diff --git a/src/bayesvalidrox/surrogate_models/engine.py b/src/bayesvalidrox/surrogate_models/engine.py
index a6531b978e082895fc053a7432fde38a68d2b1bc..d0b281f76544278024afa76d41afb8ca512d4816 100644
--- a/src/bayesvalidrox/surrogate_models/engine.py
+++ b/src/bayesvalidrox/surrogate_models/engine.py
@@ -362,7 +362,8 @@ class Engine:
             pce = True
         else:
             pce = False
-        mc_ref = True if bool(self.Model.mc_reference) else 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')
 
diff --git a/src/bayesvalidrox/surrogate_models/sequential_design.py b/src/bayesvalidrox/surrogate_models/sequential_design.py
index e020f689f4e919b0692578f458da6d1cce7613d2..29b2130305b803652e8a4a68c430539388a61318 100644
--- a/src/bayesvalidrox/surrogate_models/sequential_design.py
+++ b/src/bayesvalidrox/surrogate_models/sequential_design.py
@@ -23,6 +23,7 @@ from tqdm import tqdm
 from bayesvalidrox.bayes_inference.bayes_inference import BayesInference
 from bayesvalidrox.bayes_inference.discrepancy import Discrepancy
 from .exploration import Exploration
+from .surrogate_models import create_psi
 
 def logpdf(x, mean, cov):
     """
@@ -427,7 +428,7 @@ class SequentialDesign:
                     X_can = explore.closestPoints[idx]
 
                     # Calculate the maxmin score for the region of interest
-                    newSamples, maxminScore = explore.get_mc_samples(X_can)
+                    newSamples, maxminScore = explore.get_mc_samples(X_can) # TODO: Understand this line!
 
                     # select the requested number of samples
                     Xnew[i] = newSamples[np.argmax(maxminScore)]
@@ -550,6 +551,7 @@ class SequentialDesign:
             # ToDo: Check function to see what it does for scores/how it chooses points, so it gives as an output the
             #       scores. See how it works with exploration_scores.
             # Todo: Check if it is a minimization or maximization. (We think it is minimization)
+            # TODO: Move the final accumulation of scores outside of this function to match the other methods
             # ------- EXPLOITATION: ALPHABETIC -------
             Xnew = self.util_AlphOptDesign(allCandidates, var)
 
@@ -1249,12 +1251,12 @@ class SequentialDesign:
 
         # ------ Old Psi ------------
         univ_p_val = self.MetaModel.univ_basis_vals(oldExpDesignX)
-        Psi = self.MetaModel.create_psi(BasisIndices, univ_p_val)
+        Psi = create_psi(BasisIndices, univ_p_val)
 
         # ------ New candidates (Psi_c) ------------
         # Assemble Psi_c
         univ_p_val_c = self.MetaModel.univ_basis_vals(candidates)
-        Psi_c = self.MetaModel.create_psi(BasisIndices, univ_p_val_c)
+        Psi_c = create_psi(BasisIndices, univ_p_val_c)
 
         for idx in range(NCandidate):