diff --git a/.coverage b/.coverage
deleted file mode 100644
index 77770e96ea40081d86344c2d02b2a4fe13fed2cb..0000000000000000000000000000000000000000
Binary files a/.coverage and /dev/null differ
diff --git a/.coverage.DESKTOP-ATMEKSV.10624.XjrvRHAx b/.coverage.DESKTOP-ATMEKSV.10624.XjrvRHAx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.22780.XWOfrNvx b/.coverage.DESKTOP-ATMEKSV.22780.XWOfrNvx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.25100.XbtpWDFx b/.coverage.DESKTOP-ATMEKSV.25100.XbtpWDFx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.29056.XrBCQErx b/.coverage.DESKTOP-ATMEKSV.29056.XrBCQErx
deleted file mode 100644
index e1f7e6f75f5aa9a7a33884193b67b7cc22082200..0000000000000000000000000000000000000000
Binary files a/.coverage.DESKTOP-ATMEKSV.29056.XrBCQErx and /dev/null differ
diff --git a/.coverage.DESKTOP-ATMEKSV.34604.XhfDgthx b/.coverage.DESKTOP-ATMEKSV.34604.XhfDgthx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/.coverage.DESKTOP-ATMEKSV.37188.XIOzbGrx b/.coverage.DESKTOP-ATMEKSV.37188.XIOzbGrx
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/bayesvalidrox/bayes_inference/bayes_inference.py b/src/bayesvalidrox/bayes_inference/bayes_inference.py
index 7284348266fa88dd93a7380633437cc713e19291..c0c4cd08eb1a07f40a2a397f98cf7d46b4802368 100644
--- a/src/bayesvalidrox/bayes_inference/bayes_inference.py
+++ b/src/bayesvalidrox/bayes_inference/bayes_inference.py
@@ -961,7 +961,7 @@ class BayesInference:
         output_names = self.engine.Model.Output.names
 
         # Evaluate MetaModel on the experimental design and ValidSet
-        OutputRS, stdOutputRS = MetaModel.eval_metamodel(samples=samples)
+        OutputRS, stdOutputRS = MetaModel.eval_metamodel(samples)
 
         logLik_data = np.zeros(n_samples)
         logLik_model = np.zeros(n_samples)
diff --git a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
index 7d2cf2616ea530ceb3b95f5ca4d1d40e2f11cc11..25b1fdb5224ec42814b764ec56ad917e6a773687 100644
--- a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
+++ b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
@@ -55,8 +55,8 @@ class PCE(MetaModel):
     ----------
     input_obj : obj
         Input object with the information on the model input parameters.
-    _meta_model_type : str
-        Surrogate model types. Three surrogate model types are supported:
+    meta_model_type : str
+        PCE-surrogate model types. Two surrogate model types are supported:
         polynomial chaos expansion (`PCE`), arbitrary PCE (`aPCE`). 
         Default is PCE.
     _pce_reg_method : str
@@ -729,8 +729,6 @@ class PCE(MetaModel):
         mean_pred = {}
         std_pred = {}
         for output, values in model_dict.items():
-            print(output, values)
-
             mean = np.empty((len(samples), len(values)))
             std = np.empty((len(samples), len(values)))
             idx = 0
diff --git a/src/bayesvalidrox/surrogate_models/surrogate_models.py b/src/bayesvalidrox/surrogate_models/surrogate_models.py
index 131e59aa949da1bfb3965a640122f4d184328cb0..5a494135f7046133d2344995bffedb44cdb9197d 100644
--- a/src/bayesvalidrox/surrogate_models/surrogate_models.py
+++ b/src/bayesvalidrox/surrogate_models/surrogate_models.py
@@ -92,7 +92,6 @@ def _preprocessing_fit(fit_function):
 
         # Build the input space
         if self.InputSpace is None:
-            print(self.meta_model_type)
             self.InputSpace = InputSpace(self.input_obj, self.meta_model_type)
             n_init_samples = X.shape[0]
             self.InputSpace.n_init_samples = n_init_samples
@@ -107,14 +106,14 @@ def _preprocessing_fit(fit_function):
 
         # --- Loop through data points and fit the surrogate ---
         if self.verbose:
-            print(f"\n>>>> Training the {self._meta_model_type} metamodel "
+            print(f"\n>>>> Training the {self.meta_model_type} metamodel "
                   "started. <<<<<<\n")
 
         # Call the function for fitting
         fit_function(self, X, y, **kwargs)
 
         if self.verbose:
-            print(f"\n>>>> Training the {self._meta_model_type} metamodel"
+            print(f"\n>>>> Training the {self.meta_model_type} metamodel"
                   " sucessfully completed. <<<<<<\n")
     return decorator
 
@@ -415,6 +414,7 @@ class MetaModel:
         self.n_pca_components = None
 
         # Build general parameters
+        self._pce_deg = 1 # TODO: rename this more generally, so that it can be used for other types as well
 
         # General warnings
         if not self.is_gaussian:
diff --git a/tests/test_BayesInference.py b/tests/test_BayesInference.py
index 4b37f2cbb2776b3c881599c98a2363bc42b3f0d0..9048776c874b488446f040fae71ce03f1d973b30 100644
--- a/tests/test_BayesInference.py
+++ b/tests/test_BayesInference.py
@@ -112,7 +112,7 @@ def test_create_inference() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
 
     mod = PL()
     mod.observations = {'Z': np.array([0.45])}
diff --git a/tests/test_MCMC.py b/tests/test_MCMC.py
index e38c93b0a8cd79ad43d974e49ef6f6d90ec95559..31b24598683befbbf38c31b1cbe1e7fd22409537 100644
--- a/tests/test_MCMC.py
+++ b/tests/test_MCMC.py
@@ -54,7 +54,7 @@ def test_MCMC() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=2)
 
     mod = PL()
     mod.observations = {'Z': np.array([0.45]), 'x_values': np.array([0])}
@@ -102,7 +102,7 @@ if 0: # TODO: issue not resolved here, issue appears due to specific test setup,
     
         mm = MetaModel(inp)
         mm.fit(expdes.X, expdes.Y)
-        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=2)
     
         mod = PL()
         mod.observations = {'Z': np.array([0.45]), 'x_values': np.array([0])}
@@ -197,7 +197,7 @@ if __name__ == '__main__':
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=2)
 
     mod = PL()
     mod.observations = {'Z': np.array([0.45]), 'x_values': np.array([0])}
diff --git a/tests/test_MetaModel.py b/tests/test_MetaModel.py
index dd1eb29f26b7d7f52f9a2cf0896c61e946e86787..f2a4981ea4ee9dca1b38ad5d44ee38f043928378 100644
--- a/tests/test_MetaModel.py
+++ b/tests/test_MetaModel.py
@@ -29,7 +29,7 @@ import numpy as np
 import pytest
 import sys
 
-sys.path.append("src/")
+sys.path.append("../src/")
 
 from bayesvalidrox.surrogate_models.inputs import Input
 from bayesvalidrox.surrogate_models.input_space import InputSpace
@@ -250,7 +250,7 @@ def test_copy_meta_model_opts() -> None:
 
 #%% Test MetaModel.__select_degree
 
-#%% Test Engine._compute_pce_moments
+#%% Test MetaModel.compute_moments
 
 def test__compute_moments() -> None:
     """
@@ -262,7 +262,7 @@ def test__compute_moments() -> None:
     inp.Marginals[0].parameters = [0, 1]
     mm = MetaModel(inp)
     mm.fit([[0.2], [0.4], [0.8]], {'Z': [[0.4], [0.2], [0.5]]})
-    mm._compute_moments()
+    mm.compute_moments()
 
 
 
diff --git a/tests/test_PolynomialChaosEmulator.py b/tests/test_PolynomialChaosEmulator.py
index fea710dd78a510b7ca19a42c08fd4d3a6314f67d..61c65889722e3a1b8206d12f9ef4b5f6f165d0ea 100644
--- a/tests/test_PolynomialChaosEmulator.py
+++ b/tests/test_PolynomialChaosEmulator.py
@@ -26,7 +26,7 @@ import numpy as np
 import pytest
 import sys
 
-sys.path.append("src/")
+sys.path.append("../src/")
 
 from bayesvalidrox.surrogate_models.inputs import Input
 from bayesvalidrox.surrogate_models.input_space import InputSpace
@@ -35,6 +35,23 @@ from bayesvalidrox.surrogate_models.polynomial_chaos import PCE
 from bayesvalidrox.surrogate_models.supplementary import gaussian_process_emulator, corr_loocv_error, create_psi
 
 
+@pytest.fixture
+def PCE_1DwithInputSpace():
+    inp = Input()
+    inp.add_marginals()
+    inp.Marginals[0].dist_type = 'normal'
+    inp.Marginals[0].parameters = [0, 1]
+    #mm = MetaModel(inp)
+    samples = np.array([[0.2], [0.8]])
+    mm = PCE(inp)
+    mm.CollocationPoints = samples
+    mm.InputSpace = InputSpace(mm.input_obj, mm.meta_model_type)
+    n_init_samples = samples.shape[0]
+    mm.InputSpace.n_init_samples = n_init_samples
+    # TODO: _pce_deg not necessarily available, generalize this!
+    mm.InputSpace.init_param_space(np.max(mm._pce_deg))
+    return mm
+
 
 #%% Test MetaMod constructor on its own
 
@@ -51,66 +68,32 @@ def test_metamod() -> None:
 
 #%% Test PCE.build_metamodel
 
-def test_build_metamodel_nosamples() -> None:
-    """
-    Build PCE without collocation samples
-    """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
-    with pytest.raises(AttributeError) as excinfo:
-        mm.build_metamodel()
-    assert str(excinfo.value) == 'Please provide samples to the metamodel before building it.'
+#def test_build_metamodel_nosamples() -> None:
+#    """
+#    Build PCE without collocation samples
+#    """
+#    inp = Input()
+#    inp.add_marginals()
+#    inp.Marginals[0].dist_type = 'normal'
+#    inp.Marginals[0].parameters = [0, 1]
+#    mm = PCE(inp)
+#    with pytest.raises(AttributeError) as excinfo:
+#        mm.build_metamodel()
+#    assert str(excinfo.value) == 'Please provide samples to the metamodel before building it.'
 
 
-def test_build_metamodel() -> None:
+def test_build_metamodel(PCE_1DwithInputSpace) -> None:
     """
     Build PCE
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     mm.CollocationPoints = np.array([[0.2], [0.8]])
     mm.build_metamodel()
 
 
-def test_build_metamodel_ninitsamples() -> None:
-    """
-    Build PCE with n_init_samples
-    """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
-    mm.CollocationPoints = np.array([[0.2], [0.8]])
-    mm.build_metamodel(n_init_samples=2)
+#%% Test PCE._generate_polynomials
 
-
-
-def test_build_metamodel_coldimerr() -> None:
-    """
-    Build PCE with wrong shape collocation samples
-    """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
-    mm.CollocationPoints = [[0.2, 0.8]]
-    with pytest.raises(AttributeError) as excinfo:
-        mm.build_metamodel()
-    assert str(
-        excinfo.value) == 'The second dimension of X should be the same size as the number of marginals in the InputObj'
-
-
-#%% Test MetaMod.generate_polynomials
-
-def test_generate_polynomials_noexp() -> None:
+def test__generate_polynomials_noexp() -> None:
     """
     Generate polynomials without ExpDeg
     """
@@ -120,11 +103,11 @@ def test_generate_polynomials_noexp() -> None:
     inp.Marginals[0].parameters = [0, 1]
     mm = PCE(inp)
     with pytest.raises(AttributeError) as excinfo:
-        mm.generate_polynomials()
+        mm._generate_polynomials()
     assert str(excinfo.value) == 'Generate or add InputSpace before generating polynomials'
 
 
-def test_generate_polynomials_nodeg() -> None:
+def test__generate_polynomials_nodeg() -> None:
     """
     Generate polynomials without max_deg
     """
@@ -137,17 +120,17 @@ def test_generate_polynomials_nodeg() -> None:
     # Setup
     mm.InputSpace = InputSpace(inp)
     mm.InputSpace.n_init_samples = 2
-    mm.InputSpace.init_param_space(np.max(mm.pce_deg))
+    mm.InputSpace.init_param_space(np.max(mm._pce_deg))
     mm.ndim = mm.InputSpace.ndim
     mm.n_params = len(mm.input_obj.Marginals)
 
     # Generate
     with pytest.raises(AttributeError) as excinfo:
-        mm.generate_polynomials()
-    assert str(excinfo.value) == 'PCE cannot generate polynomials in the given scenario!'
+        mm._generate_polynomials()
+    assert str(excinfo.value) == 'MetaModel cannot generate polynomials in the given scenario!'
 
 
-def test_generate_polynomials_deg() -> None:
+def test__generate_polynomials_deg() -> None:
     """
     Generate polynomials with max_deg
     """
@@ -160,12 +143,12 @@ def test_generate_polynomials_deg() -> None:
     # Setup
     mm.InputSpace = InputSpace(inp)
     mm.InputSpace.n_init_samples = 2
-    mm.InputSpace.init_param_space(np.max(mm.pce_deg))
+    mm.InputSpace.init_param_space(np.max(mm._pce_deg))
     mm.ndim = mm.InputSpace.ndim
     mm.n_params = len(mm.input_obj.Marginals)
 
     # Generate
-    mm.generate_polynomials(4)
+    mm._generate_polynomials(4)
 
 
 #%% Test MetaMod.add_InputSpace
@@ -236,447 +219,381 @@ def test_fit_pca() -> None:
 
 #%% Test PCE.regression
 
-def test_regression() -> None:
+def test_regression(PCE_1DwithInputSpace) -> None:
     """
     Regression without a method
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
     mm.regression(samples, outputs, psi)
 
 
-def test_regression_ols() -> None:
+def test_regression_ols(PCE_1DwithInputSpace) -> None:
     """
     Regression: ols
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='ols')
+    mm._pce_reg_method = 'ols'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_olssparse() -> None:
+def test_regression_olssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: ols and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='ols', sparsity=True)
+    mm._pce_reg_method = 'ols'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
-def test_regression_ard() -> None:
+def test_regression_ard(PCE_1DwithInputSpace) -> None:
     """
     Regression: ard
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2], [0.8]])
     outputs = np.array([0.4, 0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='ard')
+    mm._pce_reg_method = 'ard'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_ardssparse() -> None:
+def test_regression_ardssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: ard and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2], [0.8]])
     outputs = np.array([0.4, 0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='ard', sparsity=True)
+    mm._pce_reg_method = 'ard'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
-def test_regression_fastard() -> None:
+def test_regression_fastard(PCE_1DwithInputSpace) -> None:
     """
     Regression: fastard
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='fastard')
+    mm._pce_reg_method = 'fastard'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_fastardssparse() -> None:
+def test_regression_fastardssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: fastard and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='fastard', sparsity=True)
+    mm._pce_reg_method = 'fastard'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
-def test_regression_brr() -> None:
+def test_regression_brr(PCE_1DwithInputSpace) -> None:
     """
     Regression: brr
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='brr')
+    mm._pce_reg_method = 'brr'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_brrssparse() -> None:
+def test_regression_brrssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: brr and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='brr', sparsity=True)
+    mm._pce_reg_method = 'brr'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
 if 0: # Could not figure out these errors, issue most likely in chosen samples/outputs
-    def test_regression_bcs() -> None:
+    def test_regression_bcs(PCE_1DwithInputSpace) -> None:
         """
         Regression: bcs
         """
-        inp = Input()
-        inp.add_marginals()
-        inp.Marginals[0].dist_type = 'normal'
-        inp.Marginals[0].parameters = [0, 1]
-        mm = PCE(inp)
+        mm = PCE_1DwithInputSpace
         samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8], [0.9]])
         outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])
-        mm.pce_deg = 3
+        mm._pce_deg = 3
         mm.CollocationPoints = samples
-        mm.build_metamodel(n_init_samples=2)
-        BasisIndices = mm.allBasisIndices[str(mm.pce_deg)][str(1.0)]
+        mm.build_metamodel()
+        BasisIndices = mm._all_basis_indices[str(mm._pce_deg)][str(1.0)]
         univ_bas = mm.univ_basis_vals(samples)
         psi = create_psi(BasisIndices, univ_bas)
     
-        mm.regression(samples, outputs, psi, reg_method='bcs')
+        mm._pce_reg_method = 'bcs'
+        mm.regression(samples, outputs, psi)
     
     
-    def test_regression_bcsssparse() -> None:
+    def test_regression_bcsssparse(PCE_1DwithInputSpace) -> None:
         """
         Regression: bcs and sparse
         """
-        inp = Input()
-        inp.add_marginals()
-        inp.Marginals[0].dist_type = 'normal'
-        inp.Marginals[0].parameters = [0, 1]
-        mm = PCE(inp)
+        mm = PCE_1DwithInputSpace
         samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8], [0.9], [1.0]])
         outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1])
     
         mm.CollocationPoints = samples
-        mm.build_metamodel(n_init_samples=2)
-        BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+        mm.build_metamodel()
+        BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
         univ_bas = mm.univ_basis_vals(samples)
         psi = create_psi(BasisIndices, univ_bas)
     
-        mm.regression(samples, outputs, psi, reg_method='bcs', sparsity=True)
+        mm._pce_reg_method = 'bcs'
+        mm.regression(samples, outputs, psi, sparsity=True)
 
 
-def test_regression_lars() -> None:
+def test_regression_lars(PCE_1DwithInputSpace) -> None:
     """
     Regression: lars
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8], [0.9], [1.0]])
     outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='lars')
+    mm._pce_reg_method = 'lars'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_larsssparse() -> None:
+def test_regression_larsssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: lars and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8], [0.9], [1.0]])
     outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='lars', sparsity=True)
+    mm._pce_reg_method = 'lars'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
-def test_regression_sgdr() -> None:
+def test_regression_sgdr(PCE_1DwithInputSpace) -> None:
     """
     Regression: sgdr
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='sgdr')
+    mm._pce_reg_method = 'sgdr'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_sgdrssparse() -> None:
+def test_regression_sgdrssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: sgdr and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='sgdr', sparsity=True)
+    mm._pce_reg_method = 'sgdr'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
 if 0: # Could not figure out these errors, issue most likely in chosen samples/outputs
-    def test_regression_omp() -> None:
+    def test_regression_omp(PCE_1DwithInputSpace) -> None:
         """
         Regression: omp
         """
-        inp = Input()
-        inp.add_marginals()
-        inp.Marginals[0].dist_type = 'normal'
-        inp.Marginals[0].parameters = [0, 1]
-        mm = PCE(inp)
+        mm = PCE_1DwithInputSpace
         samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8], [0.9], [1.0]])
         outputs = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1])
     
         mm.CollocationPoints = samples
-        mm.build_metamodel(n_init_samples=2)
-        BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+        mm.build_metamodel()
+        BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
         univ_bas = mm.univ_basis_vals(samples)
         psi = create_psi(BasisIndices, univ_bas)
     
         mm.regression(samples, outputs, psi, reg_method='omp')
     
     
-    def test_regression_ompssparse() -> None:
+    def test_regression_ompssparse(PCE_1DwithInputSpace) -> None:
         """
         Regression: omp and sparse
         """
-        inp = Input()
-        inp.add_marginals()
-        inp.Marginals[0].dist_type = 'normal'
-        inp.Marginals[0].parameters = [0, 1]
-        mm = PCE(inp)
+        mm = PCE_1DwithInputSpace
         samples = np.array([[0.2]])
         outputs = np.array([0.5])
     
         mm.CollocationPoints = samples
-        mm.build_metamodel(n_init_samples=2)
-        BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+        mm.build_metamodel()
+        BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
         univ_bas = mm.univ_basis_vals(samples)
         psi = create_psi(BasisIndices, univ_bas)
     
         mm.regression(samples, outputs, psi, reg_method='omp', sparsity=True)
 
 
-def test_regression_vbl() -> None:
+def test_regression_vbl(PCE_1DwithInputSpace) -> None:
     """
     Regression: vbl
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='vbl')
+    mm._pce_reg_method = 'vbl'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_vblssparse() -> None:
+def test_regression_vblssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: vbl and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='vbl', sparsity=True)
+    mm._pce_reg_method = 'vbl'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
-def test_regression_ebl() -> None:
+def test_regression_ebl(PCE_1DwithInputSpace) -> None:
     """
     Regression: ebl
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='ebl')
+    mm._pce_reg_method = 'ebl'
+    mm.regression(samples, outputs, psi)
 
 
-def test_regression_eblssparse() -> None:
+def test_regression_eblssparse(PCE_1DwithInputSpace) -> None:
     """
     Regression: ebl and sparse
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.5])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    mm.regression(samples, outputs, psi, reg_method='ebl', sparsity=True)
+    mm._pce_reg_method = 'ebl'
+    mm.regression(samples, outputs, psi, sparsity=True)
 
 
 #%% Test Model.update_pce_coeffs
@@ -685,107 +602,87 @@ def test_regression_eblssparse() -> None:
 
 #%% Test PCE.univ_basis_vals
 
-def test_univ_basis_vals() -> None:
+def test_univ_basis_vals(PCE_1DwithInputSpace) -> None:
     """
     Creates univariate polynomials
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2], [0.8]])
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
+    mm.build_metamodel()
     mm.univ_basis_vals(samples)
 
 
 #%% Test PCE.adaptive_regression
 
-def test_adaptive_regression_fewsamples() -> None:
+def test_adaptive_regression_fewsamples(PCE_1DwithInputSpace) -> None:
     """
     Adaptive regression, no specific method, too few samples given
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.8])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
+    mm.build_metamodel()
 
     # Evaluate the univariate polynomials on InputSpace
     mm.univ_p_val = mm.univ_basis_vals(mm.CollocationPoints)
 
     with pytest.raises(AttributeError) as excinfo:
-        mm.adaptive_regression(outputs, 0)
+        mm.adaptive_regression(samples, outputs, 0)
     assert str(excinfo.value) == ('There are too few samples for the corrected loo-cv error. Fit surrogate on at least as '
                            'many samples as parameters to use this')
 
 
-def test_adaptive_regression() -> None:
+def test_adaptive_regression(PCE_1DwithInputSpace) -> None:
     """
     Adaptive regression, no specific method
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.0], [0.1]])
     outputs = np.array([0.0, 0.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
+    mm.build_metamodel()
 
     # Evaluate the univariate polynomials on InputSpace
     mm.univ_p_val = mm.univ_basis_vals(mm.CollocationPoints)
-    mm.adaptive_regression(outputs, 0)
+    mm.adaptive_regression(samples, outputs, 0)
 
 
-def test_adaptive_regression_verbose() -> None:
+def test_adaptive_regression_verbose(PCE_1DwithInputSpace) -> None:
     """
     Adaptive regression, no specific method, verbose output
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.0], [0.1]])
     outputs = np.array([0.0, 0.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
+    mm.build_metamodel()
 
     # Evaluate the univariate polynomials on InputSpace
     mm.univ_p_val = mm.univ_basis_vals(mm.CollocationPoints)
-    mm.adaptive_regression(outputs, 0, True)
+    mm.adaptive_regression(samples, outputs, 0, True)
 
 
-def test_adaptive_regression_ols() -> None:
+def test_adaptive_regression_ols(PCE_1DwithInputSpace) -> None:
     """
     Adaptive regression, ols
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = PCE(inp)
+    mm = PCE_1DwithInputSpace
     samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8],
                         [0.9], [1.0]])
     outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
+    mm.build_metamodel()
 
     # Evaluate the univariate polynomials on InputSpace
     mm.univ_p_val = mm.univ_basis_vals(mm.CollocationPoints)
-    mm.pce_reg_method = 'ols'
-    mm.adaptive_regression(outputs, 0)
+    mm._pce_reg_method = 'ols'
+    mm.adaptive_regression(samples, outputs, 0)
 
 
 #%% Test PCE.pca_transformation
@@ -915,9 +812,9 @@ def test_copy_meta_model_opts() -> None:
 
 #%% Test PCE.__select_degree
 
-#%% Test Engine._compute_moments
+#%% Test PCE.compute_moments
 
-def test__compute_moments() -> None:
+def test_compute_moments() -> None:
     """
     Compute moments of a pce-surrogate
     """
@@ -927,10 +824,10 @@ def test__compute_moments() -> None:
     inp.Marginals[0].parameters = [0, 1]
     mm = PCE(inp)
     mm.fit([[0.2], [0.4], [0.8]], {'Z': [[0.4], [0.2], [0.5]]})
-    mm._compute_moments()
+    mm.compute_moments()
 
 
-def test__compute_moments_pca() -> None:
+def test_compute_moments_pca() -> None:
     """
     Compute moments of a pce-surrogate with pca
     """
@@ -941,7 +838,7 @@ def test__compute_moments_pca() -> None:
     mm = PCE(inp)
     mm.dim_red_method = 'pca'
     mm.fit([[0.2], [0.8]], {'Z': [[0.4, 0.4], [0.5, 0.6]]})
-    mm._compute_moments()
+    mm.compute_moments()
 
 
 #%% Test PCE.update_metamodel
diff --git a/tests/test_SequentialDesign.py b/tests/test_SequentialDesign.py
index 52e54be03a22ec15f49bbbe4b9c4a6dd950ace00..f922a77d7b78feee8ca01de8976ed4bd036607a7 100644
--- a/tests/test_SequentialDesign.py
+++ b/tests/test_SequentialDesign.py
@@ -170,7 +170,7 @@ def test_choose_next_sample() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -200,7 +200,7 @@ def test_choose_next_sample_da_spaceparallel() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -231,7 +231,7 @@ def test_choose_next_sample_da_spacenoparallel() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -262,7 +262,7 @@ def test_choose_next_sample_loo_space() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -292,7 +292,7 @@ def test_choose_next_sample_vor_space() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -327,7 +327,7 @@ def test_choose_next_sample_latin_space() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -357,7 +357,7 @@ def test_choose_next_sample_latin_alphD() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -387,7 +387,7 @@ def test_choose_next_sample_latin_alphK() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -417,7 +417,7 @@ def test_choose_next_sample_latin_alphA() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -448,7 +448,7 @@ def test_choose_next_sample_latin_VarALM() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -479,7 +479,7 @@ def test_choose_next_sample_latin_VarEIGF() -> None:
 
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -512,7 +512,7 @@ if 0:
     
         mm = MetaModel(inp)
         mm.fit(expdes.X, expdes.Y)
-        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
         mod = PL()
         
         engine = Engine(mm, mod, expdes)
@@ -543,7 +543,7 @@ def test_choose_next_sample_latin_BODMI() -> None:
     expdes.util_func = 'MI'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -576,7 +576,7 @@ def test_choose_next_sample_vor_BODMI() -> None:
     expdes.util_func = 'MI'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -614,7 +614,7 @@ def test_choose_next_sample_latin_BODALC() -> None:
     expdes.util_func = 'ALC'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -647,7 +647,7 @@ def test_choose_next_sample_latin_BODDKL() -> None:
     expdes.util_func = 'DKL'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -680,7 +680,7 @@ def test_choose_next_sample_latin_BODDPP() -> None:
     expdes.util_func = 'DPP'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -713,7 +713,7 @@ def test_choose_next_sample_latin_BODAPP() -> None:
     expdes.util_func = 'APP'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -746,7 +746,7 @@ def test_choose_next_sample_latin_BODMI_() -> None:
     expdes.util_func = 'MI'
     mm = MetaModel(inp)
     mm.fit(expdes.X, expdes.Y)
-    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
     mod = PL()
     
     engine = Engine(mm, mod, expdes)
@@ -779,7 +779,7 @@ if 0:
         expdes.util_func = 'BME'
         mm = MetaModel(inp)
         mm.fit(expdes.X, expdes.Y)
-        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
         mod = PL()
         
         engine = Engine(mm, mod, expdes)
@@ -813,7 +813,7 @@ if 0:
         expdes.util_func = 'DKL'
         mm = MetaModel(inp)
         mm.fit(expdes.X, expdes.Y)
-        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
         mod = PL()
         
         engine = Engine(mm, mod, expdes)
@@ -847,7 +847,7 @@ if 0:
         expdes.util_func = 'infEntropy'
         mm = MetaModel(inp)
         mm.fit(expdes.X, expdes.Y)
-        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm.pce_deg))
+        expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(mm._pce_deg))
         mod = PL()
         
         engine = Engine(mm, mod, expdes)
diff --git a/tests/test_supplementary.py b/tests/test_supplementary.py
index 7537b0ca38d2644748ef626650bc02cf7036324f..1eea5b70ba398cc27219c6d9184d253a30301094 100644
--- a/tests/test_supplementary.py
+++ b/tests/test_supplementary.py
@@ -14,7 +14,7 @@ import numpy as np
 import pytest
 import sys
 
-sys.path.append("src/")
+sys.path.append("../src/")
 
 from bayesvalidrox.surrogate_models.inputs import Input
 from bayesvalidrox.surrogate_models.input_space import InputSpace
@@ -22,95 +22,101 @@ from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
 from bayesvalidrox.surrogate_models.polynomial_chaos import PCE
 from bayesvalidrox.surrogate_models.supplementary import gaussian_process_emulator, corr_loocv_error, create_psi
 
+@pytest.fixture
+def PCE_withInputSpace():
+    inp = Input()
+    inp.add_marginals()
+    inp.Marginals[0].dist_type = 'normal'
+    inp.Marginals[0].parameters = [0, 1]
+    #mm = MetaModel(inp)
+    samples = np.array([[0.2], [0.8]])
+    mm = PCE(inp)
+    mm.CollocationPoints = samples
+    mm.InputSpace = InputSpace(mm.input_obj, mm.meta_model_type)
+    n_init_samples = samples.shape[0]
+    mm.InputSpace.n_init_samples = n_init_samples
+    mm.InputSpace.init_param_space(np.max(mm._pce_deg))
+    return mm
+    
+
 
 #%% Test create_psi
 
-def test_create_psi() -> None:
+def test_create_psi(PCE_withInputSpace) -> None:
     """
     Create psi-matrix
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = MetaModel(inp)
+    mm = PCE_withInputSpace#(inp, samples)
     samples = np.array([[0.2], [0.8]])
-    mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    #mm.InputSpace = InputSpace(mm.input_obj, mm.meta_model_type)
+    #n_init_samples = samples.shape[0]
+    #mm.InputSpace.n_init_samples = n_init_samples
+    # TODO: _pce_deg not necessarily available, generalize this!
+    #mm.InputSpace.init_param_space(np.max(mm._pce_deg))
+    mm.build_metamodel()
+    
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     create_psi(BasisIndices, univ_bas)
 
 #%% Test corr_loocv_error
 
-def test_corr_loocv_error_nosparse() -> None:
+def test_corr_loocv_error_nosparse(PCE_withInputSpace) -> None:
     """
     Corrected loocv error
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = MetaModel(inp)
+    mm = PCE_withInputSpace
     samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7],
                         [0.8], [0.9], [1.0]])
     outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    outs = mm.regression(samples, outputs, psi, reg_method='ebl')
+    outs = mm.regression(samples, outputs, psi)
     corr_loocv_error(outs['clf_poly'], outs['sparePsi'], outs['coeffs'],
                      outputs)
 
 
-def test_corr_loocv_error_singley() -> None:
+def test_corr_loocv_error_singley(PCE_withInputSpace) -> None:
     """
     Corrected loocv error
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = MetaModel(inp)
+    mm = PCE_withInputSpace
     samples = np.array([[0.2]])
     outputs = np.array([0.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
 
-    outs = mm.regression(samples, outputs, psi, reg_method='ols')
+    outs = mm.regression(samples, outputs, psi)
     corr_loocv_error(outs['clf_poly'], outs['sparePsi'], outs['coeffs'],
                      outputs)
 
 
-def test_corr_loocv_error_sparse() -> None:
+def test_corr_loocv_error_sparse(PCE_withInputSpace) -> None:
     """
     Corrected loocv error from sparse results
     """
-    inp = Input()
-    inp.add_marginals()
-    inp.Marginals[0].dist_type = 'normal'
-    inp.Marginals[0].parameters = [0, 1]
-    mm = MetaModel(inp)
+    mm = PCE_withInputSpace
     samples = np.array([[0.0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7],
                         [0.8], [0.9], [1.0]])
     outputs = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1])
 
     mm.CollocationPoints = samples
-    mm.build_metamodel(n_init_samples=2)
-    BasisIndices = mm.allBasisIndices[str(1)][str(1.0)]
+    mm.build_metamodel()
+    BasisIndices = mm._all_basis_indices[str(1)][str(1.0)]
     univ_bas = mm.univ_basis_vals(samples)
     psi = create_psi(BasisIndices, univ_bas)
+    mm._pce_reg_method='ebl'
 
-    outs = mm.regression(samples, outputs, psi, reg_method='ebl',
-                         sparsity=True)
+    outs = mm.regression(samples, outputs, psi, sparsity=True)
     corr_loocv_error(outs['clf_poly'], outs['sparePsi'], outs['coeffs'],
                      outputs)
 
@@ -147,7 +153,7 @@ def test_gaussian_process_emulator_autosel() -> None:
     inp.add_marginals()
     inp.Marginals[0].dist_type = 'normal'
     inp.Marginals[0].parameters = [0, 1]
-    gaussian_process_emulator([[0.2], [0.8]], [0.4, 0.5], autoSelect=True)
+    gaussian_process_emulator([[0.2], [0.8]], [0.4, 0.5], auto_select=True)
 
 
 def test_gaussian_process_emulator_varidx() -> None:
@@ -158,5 +164,5 @@ def test_gaussian_process_emulator_varidx() -> None:
     inp.add_marginals()
     inp.Marginals[0].dist_type = 'normal'
     inp.Marginals[0].parameters = [0, 1]
-    gaussian_process_emulator([[0.2], [0.8]], [0.4, 0.5], varIdx=1)
+    gaussian_process_emulator([[0.2], [0.8]], [0.4, 0.5], var_idx=1)