diff --git a/examples/OHagan-function/example_OHagan.py b/examples/OHagan-function/example_OHagan.py
index 8c7690767346db5c09e1a7b33595567c26c36a77..28010a4fcc12365779e052e425cf6ea947f4442d 100644
--- a/examples/OHagan-function/example_OHagan.py
+++ b/examples/OHagan-function/example_OHagan.py
@@ -26,15 +26,8 @@ import joblib
 import sys
 sys.path.append("../../src/bayesvalidrox/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-#from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, PostProcessing, BayesInference, Discrepancy, Engine 
+
 import matplotlib
 matplotlib.use('agg')
 
@@ -68,7 +61,7 @@ if __name__ == "__main__":
     # =====================================================
     # ======  POLYNOMIAL CHAOS EXPANSION METAMODELS  ======
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs)
+    MetaModelOpts = PCE(Inputs)
 
     # Select your metamodel method
     # 1) PCE (Polynomial Chaos Expansion) 2) aPCE (arbitrary PCE)
diff --git a/examples/analytical-function/example_analytical_function.ipynb b/examples/analytical-function/example_analytical_function.ipynb
index e63948ddc4e636a1e0eef9834c706a40f5be6890..04396bb4b57db5170296cca206bae5269b979b82 100644
--- a/examples/analytical-function/example_analytical_function.ipynb
+++ b/examples/analytical-function/example_analytical_function.ipynb
@@ -382,8 +382,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from bayesvalidrox import MetaModel\n",
-    "MetaModelOpts = MetaModel(Inputs)"
+    "from bayesvalidrox import PCE\n",
+    "MetaModelOpts = PCE(Inputs)"
    ]
   },
   {
diff --git a/examples/analytical-function/example_analytical_function_gp.py b/examples/analytical-function/example_analytical_function_gp.py
index 2d94a00a1475f45c35f15eb455dcf798090cf43f..a01ec4305d9864a8b16aa631303dbcf5189cd662 100644
--- a/examples/analytical-function/example_analytical_function_gp.py
+++ b/examples/analytical-function/example_analytical_function_gp.py
@@ -29,8 +29,7 @@ matplotlib.use('agg')
 sys.path.append("../../src/")
 print(sys.path)
 
-from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, GPESkl, PostProcessing, BayesInference, Discrepancy, \
-    Engine
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, GPESkl, PostProcessing, BayesInference, Discrepancy, Engine
 
 if __name__ == "__main__":
 
diff --git a/examples/analytical-function_sequentialdesign/example_analytical_function.ipynb b/examples/analytical-function_sequentialdesign/example_analytical_function.ipynb
index e63948ddc4e636a1e0eef9834c706a40f5be6890..04396bb4b57db5170296cca206bae5269b979b82 100644
--- a/examples/analytical-function_sequentialdesign/example_analytical_function.ipynb
+++ b/examples/analytical-function_sequentialdesign/example_analytical_function.ipynb
@@ -382,8 +382,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from bayesvalidrox import MetaModel\n",
-    "MetaModelOpts = MetaModel(Inputs)"
+    "from bayesvalidrox import PCE\n",
+    "MetaModelOpts = PCE(Inputs)"
    ]
   },
   {
diff --git a/examples/analytical-function_sequentialdesign/example_analytical_function.py b/examples/analytical-function_sequentialdesign/example_analytical_function.py
index 09261acba65231d3b9782d06fa10ce8c18a95085..8feac878e40cb728fe63bcd7ca3b1a0d7a7b98ee 100644
--- a/examples/analytical-function_sequentialdesign/example_analytical_function.py
+++ b/examples/analytical-function_sequentialdesign/example_analytical_function.py
@@ -28,19 +28,8 @@ matplotlib.use('agg')
 sys.path.append("src/")
 sys.path.append("../../src/")
 
-import bayesvalidrox
-from bayesvalidrox import PyLinkForwardModel
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, PostProcessing, BayesInference, Discrepancy, Engine
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-#from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-
-from bayesvalidrox.surrogate_models.engine import Engine
 
 if __name__ == "__main__":
 
@@ -95,7 +84,7 @@ if __name__ == "__main__":
     # =====================================================
     # ==========  DEFINITION OF THE METAMODEL  ============
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs)#, Model)
+    MetaModelOpts = PCE(Inputs)#, Model)
 
     # Select if you want to preserve the spatial/temporal depencencies
     # MetaModelOpts.dim_red_method = 'PCA'
diff --git a/examples/beam/example_beam.py b/examples/beam/example_beam.py
index 84b88255f53ef81ae5eb2861b4550fe0c68052ab..c9d608f1a2e018f8c3dd6563da9bb73376e2addc 100644
--- a/examples/beam/example_beam.py
+++ b/examples/beam/example_beam.py
@@ -18,16 +18,7 @@ Pfaffenwaldring 61
 import sys
 sys.path.append("../../src/bayesvalidrox/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-#from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-
-from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, PostProcessing, BayesInference, Discrepancy, Engine
 
 
 if __name__ == "__main__":
@@ -84,7 +75,7 @@ if __name__ == "__main__":
     # =====================================================
     # ==========  DEFINITION OF THE METAMODEL  ============
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs, Model)
+    MetaModelOpts = PCE(Inputs, Model)
 
     # Select if you want to preserve the spatial/temporal depencencies
     # MetaModelOpts.dim_red_method = 'PCA'
diff --git a/examples/borehole/example_borehole.py b/examples/borehole/example_borehole.py
index 869a053202cfb7343ee03a9a8f41b8b2fc4c8958..9eae7ffc536cd5b2495dc433b8143a394cfa43b1 100644
--- a/examples/borehole/example_borehole.py
+++ b/examples/borehole/example_borehole.py
@@ -27,15 +27,7 @@ import joblib
 import sys
 sys.path.append("../../src/bayesvalidrox/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-#from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, PostProcessing, BayesInference, Discrepancy, Engine 
 
 import matplotlib
 matplotlib.use('agg')
@@ -115,7 +107,7 @@ if __name__ == "__main__":
     # =====================================================
     # ======  POLYNOMIAL CHAOS EXPANSION METAMODELS  ======
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs, Model)
+    MetaModelOpts = PCE(Inputs, Model)
 
     # Select your metamodel method
     # 1) PCE (Polynomial Chaos Expansion) 2) aPCE (arbitrary PCE)
diff --git a/examples/convergence/example_trainingconvergence.py b/examples/convergence/example_trainingconvergence.py
index bfe8819bdb42dd0ba8e14bb0a839016a0d61a792..56cf407efa1daedcdce5c04f8ee04f7fed59410c 100644
--- a/examples/convergence/example_trainingconvergence.py
+++ b/examples/convergence/example_trainingconvergence.py
@@ -11,10 +11,8 @@ import numpy as np
 import matplotlib.pyplot as plt
 import random
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.engine import Engine
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, Engine
+
 import bayesvalidrox.surrogate_models.surrogate_models as sm
 
 from convergence import convergence
diff --git a/examples/ishigami/example_ishigami.py b/examples/ishigami/example_ishigami.py
index 72c9dc2d31e93601cbbd60bc40b88a31a73bd5b9..c74a6bcdbe9d8c42b0e57bce60f8e2670b987652 100644
--- a/examples/ishigami/example_ishigami.py
+++ b/examples/ishigami/example_ishigami.py
@@ -25,16 +25,8 @@ import joblib
 import sys
 sys.path.append("../../src/bayesvalidrox/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-#from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-
-from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, PostProcessing, BayesInference, Discrepancy, Engine 
+
 import matplotlib
 matplotlib.use('agg')
 
@@ -75,7 +67,7 @@ if __name__ == "__main__":
     # =====================================================
     # ======  POLYNOMIAL CHAOS EXPANSION METAMODELS  ======
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs, Model)
+    MetaModelOpts = PCE(Inputs, Model)
 
     # Select your metamodel method
     # 1) PCE (Polynomial Chaos Expansion) 2) aPCE (arbitrary PCE)
diff --git a/examples/only-model/example_analytical_function_NoSigmaMetaMod.py b/examples/only-model/example_analytical_function_NoSigmaMetaMod.py
index 9e8e68afba56ae605a484ea624c9decf3f691f57..0525ede53f9a0d12ab340520b988904881fab057 100644
--- a/examples/only-model/example_analytical_function_NoSigmaMetaMod.py
+++ b/examples/only-model/example_analytical_function_NoSigmaMetaMod.py
@@ -25,15 +25,8 @@ import joblib
 # Add BayesValidRox path
 sys.path.append("../../src/bayesvalidrox/")
 
-#import bayesvalidrox as bv
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-import matplotlib
+from bayesvalidrox import PyLinkForwardModel, Input, PCE, MetaModelEngine, PostProcessing, BayesInference, Discrepancy
+
 matplotlib.use('agg')
 
 
@@ -93,7 +86,7 @@ if __name__ == "__main__":
     # =====================================================
     # ==========  DEFINITION OF THE METAMODEL  ============
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs, Model)
+    MetaModelOpts = PCE(Inputs, Model)
 
     # Select if you want to preserve the spatial/temporal depencencies
     # MetaModelOpts.dim_red_method = 'PCA'
diff --git a/examples/only-model/example_analytical_function_noMetaMod.py b/examples/only-model/example_analytical_function_noMetaMod.py
index 7419d20a9115dfd2e5f7ba75276eb369c9164470..bae6018c2285ea2a3b59126524cce051b8e8cc10 100644
--- a/examples/only-model/example_analytical_function_noMetaMod.py
+++ b/examples/only-model/example_analytical_function_noMetaMod.py
@@ -24,12 +24,8 @@ import joblib
 # Add BayesValidRox path
 sys.path.append("../../src/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.engine import Engine
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
+from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, BayesInference, Discrepancy, Engine, BayesModelComparison
+
 from bayesvalidrox.bayes_inference.bayes_model_comparison import BayesModelComparison
 import matplotlib
 matplotlib.use('agg')
diff --git a/examples/pollution/example_pollution.py b/examples/pollution/example_pollution.py
index f2de52e2caf0196d8d05a2035950e2f28ff8fc70..02b25ec4d8813ab4154b4666813fd4083adfad27 100644
--- a/examples/pollution/example_pollution.py
+++ b/examples/pollution/example_pollution.py
@@ -26,12 +26,8 @@ import joblib
 import sys
 sys.path.append("../../src/bayesvalidrox/")
 
-from pylink.pylink import PyLinkForwardModel
-from surrogate_models.inputs import Input
-from surrogate_models.surrogate_models import MetaModel
-from post_processing.post_processing import PostProcessing
-from bayes_inference.bayes_inference import BayesInference
-from bayes_inference.discrepancy import Discrepancy
+from bayesvalidrox import PyLinkForwardModel, Input, PCE, PostProcessing, BayesInference, Discrepancy
+
 
 if __name__ == "__main__":
 
@@ -93,7 +89,7 @@ if __name__ == "__main__":
     # =====================================================
     # ======  POLYNOMIAL CHAOS EXPANSION METAMODELS  ======
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs, Model)
+    MetaModelOpts = PCE(Inputs, Model)
 
     # Select if you want to preserve the spatial/temporal depencencies
     # MetaModelOpts.dim_red_method = 'PCA'
diff --git a/examples/pollution/example_valid_pollution.py b/examples/pollution/example_valid_pollution.py
index 74b87787b741a0c990c6d1ccbf06596b43256a57..ea3d45ecd4e240fa453a827cf3fef0d8d7349beb 100644
--- a/examples/pollution/example_valid_pollution.py
+++ b/examples/pollution/example_valid_pollution.py
@@ -23,15 +23,8 @@ import joblib
 import sys
 sys.path.append("../../src/bayesvalidrox/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-#from bayesvalidrox.surrogate_models.meta_model_engine import MetaModelEngine
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox import PyLinkForwardModel, Input, PCE, PostProcessing, BayesInference, Discrepancy, Engine, ExpDesigns
+
 
 if __name__ == "__main__":
 
@@ -93,7 +86,7 @@ if __name__ == "__main__":
     # =====================================================
     # ======  POLYNOMIAL CHAOS EXPANSION METAMODELS  ======
     # =====================================================
-    MetaModelOpts = MetaModel(Inputs)
+    MetaModelOpts = PCE(Inputs)
 
     # Select if you want to preserve the spatial/temporal depencencies
     MetaModelOpts.dim_red_method = 'PCA'
diff --git a/examples/principal_component_analysis/example_principalcomponentanalysis.py b/examples/principal_component_analysis/example_principalcomponentanalysis.py
index 69cb4aa8fa2846e4d20a2cdb16e749716f4b8097..bc36876a1e39444511713639669c95546f84087e 100644
--- a/examples/principal_component_analysis/example_principalcomponentanalysis.py
+++ b/examples/principal_component_analysis/example_principalcomponentanalysis.py
@@ -24,15 +24,8 @@ import joblib
 # Add BayesValidRox path
 sys.path.append("../../src/")
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.engine import Engine
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
-from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.bayes_inference.bayes_model_comparison import BayesModelComparison
+from bayesvalidrox import PyLinkForwardModel, Input, PCE, PostProcessing, BayesInference, Discrepancy, Engine, ExpDesigns, BayesModelComparison
+
 import matplotlib
 matplotlib.use('agg')
 
@@ -104,7 +97,7 @@ if __name__ == "__main__":
     # ------------------------------------------------
     # ------------- PCE Specification ----------------
     # ------------------------------------------------
-    MetaModelOpts = MetaModel(Inputs)
+    MetaModelOpts = PCE(Inputs)
 
     # Select if you want to preserve the spatial/temporal depencencies
     MetaModelOpts.dim_red_method = 'PCA'
diff --git a/examples/sinusoidal-function/test.ipynb b/examples/sinusoidal-function/test.ipynb
deleted file mode 100644
index 0ff49e1e71323ec44763301fed8e84353a84a173..0000000000000000000000000000000000000000
--- a/examples/sinusoidal-function/test.ipynb
+++ /dev/null
@@ -1,207 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "There are no estimations of surrogate uncertainty available for the chosen regression options. This might lead to issues in later steps.\n",
-      "\n",
-      " Now the forward model needs to be run!\n",
-      "\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "Computing orth. polynomial coeffs: 100%|##########| 1/1 [00:01<00:00,  1.61s/it]\n",
-      "Running forward model valid: 100%|██████████| 3/3 [00:00<00:00,  4.28it/s]\n",
-      "Running forward model valid: 100%|██████████| 1000/1000 [00:01<00:00, 775.69it/s]\n",
-      "Running forward model validSet: 100%|██████████| 3/3 [00:00<00:00,  4.29it/s]\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n",
-      ">>>>> Errors of Z <<<<<\n",
-      "\n",
-      "Index  |  RMSE   |  Validation Error\n",
-      "-----------------------------------\n",
-      "1  |  1.319e-07  |  4.873e-13\n"
-     ]
-    }
-   ],
-   "source": [
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "import sys\n",
-    "import matplotlib\n",
-    "import matplotlib.pyplot as plt\n",
-    "matplotlib.use('agg')\n",
-    "\n",
-    "sys.path.append(\"../../src/\")\n",
-    "\n",
-    "from bayesvalidrox import PyLinkForwardModel, Input, ExpDesigns, PCE, PostProcessing, Discrepancy, Engine\n",
-    "\n",
-    "if __name__ == \"__main__\":\n",
-    "    \n",
-    "    # =====================================================\n",
-    "    # =============   COMPUTATIONAL MODEL  ================\n",
-    "    # =====================================================\n",
-    "    Model=PyLinkForwardModel()\n",
-    "\n",
-    "    Model.link_type = 'Function'\n",
-    "    Model.py_file = 'sinusoidal_function'\n",
-    "    Model.name = 'SinusoidalFunc'\n",
-    "\n",
-    "    Model.Output.names = ['Z']\n",
-    "\n",
-    "    # =====================================================\n",
-    "    # =========   PROBABILISTIC INPUT MODEL  ==============\n",
-    "    # =====================================================\n",
-    "    # Define the uncertain parameters with their mean and\n",
-    "    # standard deviation\n",
-    "    Inputs = Input()\n",
-    "\n",
-    "    Inputs.add_marginals()\n",
-    "    Inputs.Marginals[0].name = 'X'\n",
-    "    Inputs.Marginals[0].dist_type = 'uniform'\n",
-    "    Inputs.Marginals[0].parameters = [0,10]\n",
-    "\n",
-    "    # =====================================================\n",
-    "    # ==========  DEFINITION OF THE METAMODEL  ============\n",
-    "    # =====================================================\n",
-    "    MetaModelOpts = PCE(Inputs)\n",
-    "\n",
-    "    MetaModelOpts.meta_model_type = 'aPCE'\n",
-    "\n",
-    "    # ------------------------------------------------\n",
-    "    # ------------- PCE Specification ----------------\n",
-    "    # ------------------------------------------------\n",
-    "    MetaModelOpts._pce_reg_method = 'OLS'\n",
-    "\n",
-    "    MetaModelOpts.bootstrap_method = 'fast'\n",
-    "    MetaModelOpts.n_bootstrap_itrs = 1\n",
-    "\n",
-    "    MetaModelOpts._pce_deg = 15\n",
-    "\n",
-    "    MetaModelOpts._pce_q_norm = 1\n",
-    "\n",
-    "    # ------------------------------------------------\n",
-    "    # ------ Experimental Design Configuration -------\n",
-    "    # ------------------------------------------------\n",
-    "    ExpDesign = ExpDesigns(Inputs)\n",
-    "\n",
-    "    #ExpDesign.method = 'sequential'\n",
-    "    ExpDesign.n_init_samples = 100\n",
-    "\n",
-    "    ExpDesign.sampling_method = 'latin_hypercube' #is this correct? \n",
-    "\n",
-    "    ExpDesign.n_new_samples = 1\n",
-    "    ExpDesign.n_max_samples = 3\n",
-    "    ExpDesign.mod_LOO_threshold = 1e-16\n",
-    "\n",
-    "    ExpDesign.tradeoff_scheme = None\n",
-    "    ExpDesign.explore_method = 'latin_hypercube' #is this correct? How do I choose this?\n",
-    "\n",
-    "    ExpDesign.n_canddidate = 1000\n",
-    "    ExpDesign.n_cand_groups = 4\n",
-    "\n",
-    "    # -------- Exploitation ------\n",
-    "    ExpDesign.exploit_method = 'Space-filling'\n",
-    "\n",
-    "    #BayesOptDesign -> when data is available \n",
-    "\n",
-    "    # VarBasedOptDesign -> when data is not available\n",
-    "    ExpDesign.util_func = 'ALM'\n",
-    "\n",
-    "    ExpDesign.valid_samples = np.load(\"./data/valid_samples.npy\")\n",
-    "    ExpDesign.valid_model_runs = {'Z': np.load(\"./data/valid_samples.npy\")}\n",
-    "\n",
-    "    # >>>>>>>>>>>>>>>>>>>>>> Build Surrogate <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n",
-    "    MetaModelOpts.ExpDesign = ExpDesign\n",
-    "    egine = Engine(MetaModelOpts, Model, ExpDesign)\n",
-    "    egine.start_engine()\n",
-    "    egine.train_normal()\n",
-    "\n",
-    "    # =====================================================\n",
-    "    # =========  POST PROCESSING OF METAMODELS  ===========\n",
-    "    # =====================================================\n",
-    "    PostPCE = PostProcessing(egine)\n",
-    "\n",
-    "    #Plot to check validation visually\n",
-    "    PostPCE.valid_metamodel(n_samples=3)\n",
-    "\n",
-    "    #Check the quality of your regression model \n",
-    "    PostPCE.check_reg_quality()\n",
-    "\n",
-    "    #Compute and print RMSE error\n",
-    "    PostPCE.check_accuracy(n_samples=3)\n",
-    "\n",
-    "    PostPCE.plot_moments()\n",
-    "\n",
-    "    # if MetaModelOpts.ExpDesign.method == 'sequential':\n",
-    "    #     PostPCE.plot_seq_design_diagnostics()\n",
-    "\n",
-    "    total_sobol = PostPCE.sobol_indices()\n",
-    "\n",
-    "    X = list(ExpDesign.X)\n",
-    "    Y = list(ExpDesign.Y.values())\n",
-    "\n",
-    "    X = [float(x[0]) for x in X]\n",
-    "    Y = [float(y) for sublist in Y for y in sublist]\n",
-    "\n",
-    "    #print(f\"x-values {len(X)} \\n y-values {len(Y)}\")\n",
-    "\n",
-    "    plt.plot(X, Y, 'g', label='PCE predictor - LSTSQ')\n",
-    "    plt.scatter(X, Y, label='training data')\n",
-    "    #plt.plot(x_, f, 'm', label='function')\n",
-    "    plt.title('PCE surrogate - prediction accuracy')\n",
-    "    plt.legend()\n",
-    "    plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.plot(X, Y, 'g', label='PCE predictor - LSTSQ')\n",
-    "plt.scatter(X, Y, label='training data')\n",
-    "#plt.plot(x_, f, 'm', label='function')\n",
-    "plt.title('PCE surrogate - prediction accuracy')\n",
-    "plt.legend()\n",
-    "plt.show()"
-   ]
-  }
- ],
- "metadata": {
-  "kernelspec": {
-   "display_name": "meta_models",
-   "language": "python",
-   "name": "python3"
-  },
-  "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.10.14"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/examples/umbridge_tsunamitutorial/example_umbridge_testmodel.py b/examples/umbridge_tsunamitutorial/example_umbridge_testmodel.py
index 22a0f750e96c4b7bf14913eb97f0040fb360a257..5f813b3224e76a6c8da799c750991904566dba85 100644
--- a/examples/umbridge_tsunamitutorial/example_umbridge_testmodel.py
+++ b/examples/umbridge_tsunamitutorial/example_umbridge_testmodel.py
@@ -12,21 +12,15 @@ An example for using UM-Bridge with an explicitly defined model can be found in
 import joblib
 import numpy as np 
 import umbridge
-import bayesvalidrox as bv
 
-from bayesvalidrox.pylink.pylink import PyLinkForwardModel
-from bayesvalidrox.surrogate_models.inputs import Input
-from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
-from bayesvalidrox.surrogate_models.surrogate_models import MetaModel
-from bayesvalidrox.post_processing.post_processing import PostProcessing
-from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox import PyLinkForwardModel, Input, PCE, PostProcessing, Engine, ExpDesigns
 
 
 if __name__ == '__main__':
     
     # This model has 2 inputs and four outputs
     n_prior_sample = 1000
-    priors = bv.Input()
+    priors = Input()
     priors.add_marginals()
     priors.Marginals[0].name = 'x'
     priors.Marginals[0].input_data = np.random.uniform(50,150,n_prior_sample)  
@@ -41,7 +35,7 @@ if __name__ == '__main__':
     model.x_values = np.array([0])
     
     # Create the surrogate
-    surrogate_opts = MetaModel(priors)
+    surrogate_opts = PCE(priors)
     
     # Select the surrogate type and properties
     surrogate_opts.meta_model_type = 'aPCE'
diff --git a/examples/umbridge_tsunamitutorial/example_umbridge_tsunamitutorial.py b/examples/umbridge_tsunamitutorial/example_umbridge_tsunamitutorial.py
index c456787d0c5d9d9714083b406895b8096a4d16ce..0d9912b94ece22370f2deed22e2045adaeeeabb9 100644
--- a/examples/umbridge_tsunamitutorial/example_umbridge_tsunamitutorial.py
+++ b/examples/umbridge_tsunamitutorial/example_umbridge_tsunamitutorial.py
@@ -19,12 +19,13 @@ import pandas as pd
 import umbridge
 
 import bayesvalidrox as bv
+from bayesvalidrox import PCE, PostProcessing, Engine, ExpDesigns, Input, PyLinkForwardModel
 from tsunami_model import tsunami_model
 
 if __name__ == '__main__':
     # This model has 2 inputs and four outputs
     n_prior_sample = 1000
-    priors = bv.Input()
+    priors = Input()
     priors.add_marginals()
     priors.Marginals[0].name = 'x'
     priors.Marginals[0].input_data = np.random.uniform(50,150,n_prior_sample)  
@@ -33,7 +34,7 @@ if __name__ == '__main__':
     priors.Marginals[1].input_data = np.random.uniform(50,150,n_prior_sample)
     
     # Define the model - level 0
-    model0 = bv.PyLinkForwardModel()
+    model0 = PyLinkForwardModel()
     model0.link_type = 'function'
     model0.py_file = 'tsunami_model'
     model0.name = 'tsunami_model'
@@ -41,7 +42,7 @@ if __name__ == '__main__':
     #model.observations = data_dict1_times
     
     # Define the model - level 1
-    model1 = bv.PyLinkForwardModel()
+    model1 = PyLinkForwardModel()
     model1.link_type = 'function'
     model1.py_file = 'tsunami_model1'
     model1.name = 'tsunami_model1'
@@ -49,7 +50,7 @@ if __name__ == '__main__':
     #model.observations = data_dict1_times
     
     # Create the surrogate
-    surrogate_opts0 = bv.MetaModel(priors, model0)
+    surrogate_opts0 = PCE(priors, model0)
     
     # Select your metamodel method
     surrogate_opts0.meta_model_type = 'aPCE'
@@ -58,21 +59,21 @@ if __name__ == '__main__':
     surrogate_opts0.pce_q_norm = 0.4#1.0
 
     # Define ExpDesign - this is the same for both models       
-    ExpDesign0 = bv.ExpDesigns(priors)
+    ExpDesign0 = ExpDesigns(priors)
     ExpDesign0.method = 'normal'
     ExpDesign0.n_init_samples = 50
     ExpDesign0.sampling_method = 'latin-hypercube'
     ExpDesign1 = copy.deepcopy(ExpDesign0)
     
     # Start and run the engine
-    engine0 = bv.Engine(surrogate_opts0, model0, ExpDesign0)
+    engine0 = Engine(surrogate_opts0, model0, ExpDesign0)
     engine0.start_engine()
     engine0.train_normal(parallel = False)
     print('Surrogate 0 completed')
     print('')
     
     # Create the surrogate
-    surrogate_opts1 = bv.MetaModel(priors, model1)
+    surrogate_opts1 = PCE(priors, model1)
     
     # Select your metamodel method
     surrogate_opts1.meta_model_type = 'aPCE'
@@ -81,7 +82,7 @@ if __name__ == '__main__':
     surrogate_opts1.pce_q_norm = 0.4#1.0
 
     # Start and run the engine
-    engine1 = bv.Engine(surrogate_opts1, model1, ExpDesign1)
+    engine1 = Engine(surrogate_opts1, model1, ExpDesign1)
     engine1.start_engine()
     engine1.train_normal(parallel = False)
     print('Surrogate 1 completed')
@@ -96,7 +97,7 @@ if __name__ == '__main__':
         joblib.dump(engine1.MetaModel, output, 2)
     
     # Post processing on model 1
-    L2_PostPCE = bv.PostProcessing(engine1)
+    L2_PostPCE = PostProcessing(engine1)
     L2_PostPCE.plot_moments(plot_type='line')
     # Plot to check validation visually.
     L2_PostPCE.valid_metamodel(n_samples=1)
@@ -118,12 +119,12 @@ if __name__ == '__main__':
     # Bayesian Inference
     # Set uncertainty
     obsData = pd.DataFrame(true_data_nox, columns=model0.Output.names)
-    DiscrepancyOpts = bv.Discrepancy('')
+    DiscrepancyOpts = Discrepancy('')
     DiscrepancyOpts.type = 'Gaussian'
     DiscrepancyOpts.parameters = (obsData*0.15)**2
     
     # Parameter estimation / single model validation via TOM for model 1
-    BayesOpts = bv.BayesInference(engine1)
+    BayesOpts = BayesInference(engine1)
     BayesOpts.emulator= True
     BayesOpts.plot_post_pred = True
     #BayesOpts.inference_method = 'rejection'
@@ -167,7 +168,7 @@ if __name__ == '__main__':
         }
 
     # Run model comparison
-    BayesOpts = bv.BayesModelComparison(
+    BayesOpts = BayesModelComparison(
         justifiability=True,
         n_bootstarp=10,
         just_n_meas=2
diff --git a/examples/umbridge_tsunamitutorial/testmodel.pk1 b/examples/umbridge_tsunamitutorial/testmodel.pk1
new file mode 100644
index 0000000000000000000000000000000000000000..6021864ade548d6adb3202eaa24924932eef03d7
Binary files /dev/null and b/examples/umbridge_tsunamitutorial/testmodel.pk1 differ
diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index 0d7fbad16b8d25482f1d350466cb8e713373d157..1b244516d39a9bf478dc165911f3b3717b85d9e7 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -1146,7 +1146,7 @@ class PostProcessing:
             # the total number of explanatory variables in the model
             # (not including the constant term)
             length_list = []
-            for key, value in PCEModel.coeffs_dict['b_1'][key].items():
+            for key, value in PCEModel._coeffs_dict['b_1'][key].items():
                 length_list.append(len(value))
             n_predictors = min(length_list)
             n_samples = x_val.shape[0]