diff --git a/.coverage b/.coverage
index 7bd62455cae2ec79af9cbc60f552afe43fca17d4..32823f8ffc8135d53a2545e0b6b25b62e36da705 100644
Binary files a/.coverage and b/.coverage differ
diff --git a/.coverage.DESKTOP-ATMEKSV.14516.XqtXatIx b/.coverage.DESKTOP-ATMEKSV.14516.XqtXatIx
new file mode 100644
index 0000000000000000000000000000000000000000..e1f7e6f75f5aa9a7a33884193b67b7cc22082200
Binary files /dev/null and b/.coverage.DESKTOP-ATMEKSV.14516.XqtXatIx differ
diff --git a/Outputs_Bayes_None_Calib/emcee_sampler.h5 b/Outputs_Bayes_None_Calib/emcee_sampler.h5
index 07e32bcea1754bf24449cccb9a33a415021ec432..20c29ce06ef054b40b9091828f8b2bffe8d63972 100644
Binary files a/Outputs_Bayes_None_Calib/emcee_sampler.h5 and b/Outputs_Bayes_None_Calib/emcee_sampler.h5 differ
diff --git a/examples/analytical-function/example_analytical_function.py b/examples/analytical-function/example_analytical_function.py
index 09261acba65231d3b9782d06fa10ce8c18a95085..52d8bb0042b08680255bb247c36dd9147919437f 100644
--- a/examples/analytical-function/example_analytical_function.py
+++ b/examples/analytical-function/example_analytical_function.py
@@ -25,12 +25,10 @@ matplotlib.use('agg')
 
 # import bayesvalidrox
 # Add BayesValidRox path
-sys.path.append("src/")
 sys.path.append("../../src/")
 
-import bayesvalidrox
-from bayesvalidrox import PyLinkForwardModel
-
+from bayesvalidrox.surrogate_models.inputs import Input
+from bayesvalidrox.surrogate_models.input_space import InputSpace
 from bayesvalidrox.pylink.pylink import PyLinkForwardModel
 from bayesvalidrox.surrogate_models.inputs import Input
 from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
@@ -231,8 +229,9 @@ if __name__ == "__main__":
     MetaModelOpts.ExpDesign = ExpDesign
     engine = Engine(MetaModelOpts, Model, ExpDesign)
     engine.start_engine()
-    #engine.train_sequential()
     engine.train_normal()
+    stop
+    engine.train_sequential()
 
     # Load the objects
     # with open(f"PCEModel_{Model.name}.pkl", "rb") as input:
diff --git a/src/bayesvalidrox/surrogate_models/engine.py b/src/bayesvalidrox/surrogate_models/engine.py
index 24aa0c0d7e71e9201215d4d7292377724e2869c6..022fdb63c04c6782f49148dac9a81ccdc68fcf35 100644
--- a/src/bayesvalidrox/surrogate_models/engine.py
+++ b/src/bayesvalidrox/surrogate_models/engine.py
@@ -200,7 +200,7 @@ class Engine:
         # Prepare X samples 
         # For training the surrogate use ExpDesign.X_tr, ExpDesign.X is for the model to run on 
         ExpDesign.generate_ED(ExpDesign.n_init_samples,
-                              transform=True,
+                              #transform=True,
                               max_pce_deg=np.max(MetaModel.pce_deg))
 
         # Run simulations at X 
@@ -333,7 +333,7 @@ class Engine:
             pce = True
         else:
             pce = False
-        mc_ref = True if bool(self.Model.mc_reference) else False
+        mc_ref = True if self.Model.mc_reference is not None else False
         if mc_ref:
             self.Model.read_observation('mc_ref')
 
@@ -347,14 +347,15 @@ class Engine:
         output_name = self.out_names
 
         # Setup the Sequential Design object
-        self.SeqDes = SequentialDesign(self.MetaModel, self.Model, self.ExpDesign)
+        self.SeqDesign = SequentialDesign(self.MetaModel, self.Model, self.ExpDesign)
 
         # Handle if only one UtilityFunctions is provided
         if not isinstance(util_func, list):
             util_func = [self.ExpDesign.util_func]
 
         # Read observations or MCReference
-        # TODO: recheck the logic in this if statement
+        # TODO: recheck the logic in this if statement 
+        # ToDo: this could perhaps be moved into the Sequential Design
         if (len(self.Model.observations) != 0 or self.Model.meas_file is not None) and hasattr(self.MetaModel,
                                                                                                'Discrepancy'):
             self.observations = self.Model.read_observation()
diff --git a/src/bayesvalidrox/surrogate_models/sequential_design.py b/src/bayesvalidrox/surrogate_models/sequential_design.py
index abe3f1a13aba02d0e75efc5026e7ded13c264821..5f969c28f56cf61701da666495c332424f8b77b0 100644
--- a/src/bayesvalidrox/surrogate_models/sequential_design.py
+++ b/src/bayesvalidrox/surrogate_models/sequential_design.py
@@ -24,6 +24,9 @@ from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
 from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
 from .exploration import Exploration
 
+# TODO: These first two functions are direct copies from functions in .engine.
+#       Should they be called from there?
+
 def logpdf(x, mean, cov):
     """
     Computes the likelihood based on a multivariate normal distribution.
@@ -86,6 +89,7 @@ def subdomain(Bounds, n_new_samples):
 
     return Subdomains
 
+
 class SequentialDesign:
     """
     Contains options for choosing the next training sample iteratively.
@@ -93,12 +97,13 @@ class SequentialDesign:
     """
 
     def __init__(self, MetaMod, Model, ExpDes, parallel=False):
-        self.MetaModel = MetaMod # Surrogate should be trained
+        self.MetaModel = MetaMod  # Surrogate should be trained
         self.Model = Model
         self.ExpDesign = ExpDes
         self.parallel = parallel
 
         # Init other parameters
+        self.out_names = self.Model.Output.names
 
     def start_seqdesign(self) -> None:
         """
@@ -109,8 +114,19 @@ class SequentialDesign:
         None
 
         """
-        None
 
+        # Read observations or MCReference
+        # TODO: recheck the logic in this if statement 
+        # ToDo: this could perhaps be moved into the Sequential Design
+        if (len(self.Model.observations) != 0 or self.Model.meas_file is not None) and hasattr(self.MetaModel,
+                                                                                               'Discrepancy'):
+            self.observations = self.Model.read_observation()
+            obs_data = self.observations
+        else:
+            obs_data = []
+            # TODO: TotalSigma2 not defined if not in this else???
+            # TODO: no self.observations if in here
+            TotalSigma2 = {}
 
     # -------------------------------------------------------------------------
     def choose_next_sample(self, sigma2=None, n_candidates=5, var='DKL'):
@@ -137,8 +153,10 @@ class SequentialDesign:
         Xnew : array (n_samples, n_params)
             Selected new training point(s).
         """
+        self.start_seqdesign()
 
         # Initialization
+        # TODO: check if all of these variables are necessary, or if some should be class attributes
         Bounds = self.ExpDesign.bound_tuples
         n_new_samples = self.ExpDesign.n_new_samples
         explore_method = self.ExpDesign.explore_method
@@ -325,7 +343,6 @@ class SequentialDesign:
                 if not hasattr(self.ExpDesign, 'max_func_itr'):
                     raise AttributeError('max_func_itr not given to the experimental design')
 
-
                 # Create a sample pool for rejection sampling
                 # ToDo: remove from here, add only to BayesOptDesign option
                 MCsize = 15000
@@ -379,7 +396,7 @@ class SequentialDesign:
             # Accumulate the samples
             # ToDo: Stop assuming 2 sets of samples (should only be 1)
             finalCandidates = np.concatenate((allCandidates, candidates), axis=0)  # ToDo: Remove
-            finalCandidates = np.unique(finalCandidates, axis=0)                   # ToDo: Remove
+            finalCandidates = np.unique(finalCandidates, axis=0)  # ToDo: Remove
 
             # Calculations take into account both exploration and exploitation
             # samples without duplicates
@@ -416,7 +433,7 @@ class SequentialDesign:
             # find an optimal point subset to add to the initial design by
             # maximization of the utility score and taking care of NaN values
             temp = totalScore.copy()
-            temp[np.isnan(totalScore)] = -np.inf                # Since we are maximizing
+            temp[np.isnan(totalScore)] = -np.inf  # Since we are maximizing
             sorted_idxtotalScore = np.argsort(temp)[::-1]
             bestIdx = sorted_idxtotalScore[:n_new_samples]
 
@@ -728,7 +745,6 @@ class SequentialDesign:
                                                       var)
         return index, -1 * U_J_d
 
-
     # -------------------------------------------------------------------------
     def util_VarBasedDesign(self, X_can, index, util_func='Entropy'):
         """
@@ -1757,4 +1773,4 @@ class SequentialDesign:
         al_unique_index = np.arange(prior_samples[:self.mc_samples + n_tp, :].shape[0])[aux2_]
         al_unique_index = al_unique_index[:self.mc_samples]
 
-        return al_unique_index
\ No newline at end of file
+        return al_unique_index
diff --git a/tests/test_BayesInference.py b/tests/test_BayesInference.py
index c31a9830e7d7e4d186936d955920885adfd39300..25c0f7ddbe71575f509fd3fec81ef7f623712bc9 100644
--- a/tests/test_BayesInference.py
+++ b/tests/test_BayesInference.py
@@ -22,6 +22,7 @@ class BayesInference:
     _plot_max_a_posteriori  Need working model to test this
     _plot_post_predictive   - x
 """
+
 import sys
 import pytest
 import numpy as np
@@ -40,8 +41,6 @@ from bayesvalidrox.bayes_inference.mcmc import MCMC
 from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
 from bayesvalidrox.bayes_inference.bayes_inference import _logpdf, _kernel_rbf
 
-import matplotlib as mpl
-mpl.rcParams.update(mpl.rcParamsDefault)
 
 
 #%% Test _logpdf
@@ -913,8 +912,7 @@ def test_plot_log_BME() -> None:
     engine = Engine(mm, mod, expdes)
 
     bi = BayesInference(engine)
-    #bi.log_BME = np.array([[0, 0.2], [0, 0.2]])
-    bi.log_BME = np.array([0, 0.2, 0, 0.2])
+    bi.log_BME = np.array([[0, 0.2], [0, 0.2]])
     bi.n_tot_measurement = 1
     bi.plot_log_BME()
 
@@ -938,7 +936,7 @@ def test_plot_log_BME_noemulator() -> None:
     engine = Engine(mm, mod, expdes)
 
     bi = BayesInference(engine)
-    bi.log_BME = np.array([0, 0.2, 0, 0.2])
+    bi.log_BME = np.array([[0, 0.2], [0, 0.2]])
     bi.n_tot_measurement = 1
     bi.emulator = False
     bi.plot_log_BME()
diff --git a/tests/test_InputSpace.py b/tests/test_InputSpace.py
index 7703f1141c70fe028444debd9b17d4ecfed12adf..40f26956376d07e59161cff6418aab1cab2c898f 100644
--- a/tests/test_InputSpace.py
+++ b/tests/test_InputSpace.py
@@ -13,11 +13,12 @@ import sys
 import pytest
 import numpy as np
 
+sys.path.append("src/")
+sys.path.append("../src/")
+
 from bayesvalidrox.surrogate_models.inputs import Input
 from bayesvalidrox.surrogate_models.input_space import InputSpace
 
-sys.path.append("src/")
-sys.path.append("../src/")
 
 
 #%% Test ExpDesign.check_valid_input
diff --git a/tests/test_MetaModel.py b/tests/test_MetaModel.py
index 3b67f6699e9867ef61264d2d45f4b61d24a660d8..bf6c3c39ccd8a7251cdb04e1437c583194c7c2dc 100644
--- a/tests/test_MetaModel.py
+++ b/tests/test_MetaModel.py
@@ -29,12 +29,13 @@ import numpy as np
 import pytest
 import sys
 
+sys.path.append("src/")
+
 from bayesvalidrox.surrogate_models.inputs import Input
 from bayesvalidrox.surrogate_models.input_space import InputSpace
 from bayesvalidrox.surrogate_models.surrogate_models import MetaModel, corr_loocv_error, create_psi
 from bayesvalidrox.surrogate_models.surrogate_models import gaussian_process_emulator
 
-sys.path.append("src/")
 
 
 #%% Test MetaMod constructor on its own
diff --git a/tests/test_SequentialDesign.py b/tests/test_SequentialDesign.py
new file mode 100644
index 0000000000000000000000000000000000000000..3c9c5fb731c0e53946d284669a04f7fa7736c455
--- /dev/null
+++ b/tests/test_SequentialDesign.py
@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+"""
+Test the SequentialDesign class for bayesvalidrox
+
+Tests are available for the following functions
+    logpdf               - x
+    subdomain            - x
+Engine:
+    start_seqdesign         
+    choose_next_sample
+        plotter
+    tradoff_weights      - x
+    run_util_func
+    util_VarBasedDesign
+    util_BayesianActiveDesign
+    util_BayesianDesign
+    dual_annealing
+    util_AlphOptDesign
+    _normpdf            - x    Also move outside the class?
+    _corr_factor_BME           Not used again in this class
+    _posteriorPlot      - x
+    _BME_Calculator     - x
+    _validError         - x
+    _error_Mean_Std     - x 
+    _select_indices
+
+
+"""
+
+import sys
+import pytest
+import numpy as np
+import pandas as pd
+
+sys.path.append("../src/")
+
+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.pylink.pylink import PyLinkForwardModel as PL
+from bayesvalidrox.surrogate_models.engine import Engine
+from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
+
+
+
+
+#%% Main runs
+if __name__ == '__main__':
+    inp = Input()
+    inp.add_marginals()
+    inp.Marginals[0].dist_type = 'normal'
+    inp.Marginals[0].parameters = [0, 1]
+    # prior_samples = np.swapaxes(np.array([np.random.normal(0,1,10)]),0,1)
+
+    expdes = ExpDesigns(inp)
+    expdes.init_param_space(max_deg=1)
+    expdes.n_init_samples = 2
+    expdes.n_max_samples = 4
+    expdes.X = np.array([[0], [1], [0.5]])
+    expdes.Y = {'Z': [[0.4], [0.5], [0.45]]}
+    expdes.x_values = np.array([0])  # Error in plots if this is not
+
+    mm = MetaModel(inp)
+    mm.n_params = 1
+    mm.fit(expdes.X, expdes.Y)
+    expdes.generate_ED(expdes.n_init_samples, max_pce_deg=np.max(1))
+    # y_hat, y_std = mm.eval_metamodel(prior_samples)
+
+    mod = PL()
+    mod.observations = {'Z': np.array([0.45])}
+    mod.observations = {'Z': np.array([0.45]), 'x_values': np.array([0])}  # Error if x_values not given
+    mod.Output.names = ['Z']
+    mod.n_obs = 1
+
+    engine = Engine(mm, mod, expdes)
+
+    sigma2Dict = {'Z': np.array([0.05])}
+    sigma2Dict = pd.DataFrame(sigma2Dict, columns=['Z'])
+    obsData = pd.DataFrame({'Z': np.array([0.45]), 'x_values': np.array([0])}, columns=mod.Output.names)
+    DiscrepancyOpts = Discrepancy('')
+    DiscrepancyOpts.type = 'Gaussian'
+    DiscrepancyOpts.parameters = (obsData * 0.15) ** 2
+    DiscrepancyOpts.opt_sigma = 'B'