diff --git a/src/bayesvalidrox/surrogate_models/eval_rec_rule.py b/src/bayesvalidrox/surrogate_models/eval_rec_rule.py
index b583c7eb2ec58d55d19b34130812730d21a12368..730a3818d58479cf20283aa752a5ab74719c7cab 100644
--- a/src/bayesvalidrox/surrogate_models/eval_rec_rule.py
+++ b/src/bayesvalidrox/surrogate_models/eval_rec_rule.py
@@ -183,7 +183,7 @@ def eval_univ_basis(x, max_deg, poly_types, apoly_coeffs=None):
     """
     # Initilize the output array
     n_samples, n_params = x.shape
-    univ_vals = np.zeros((n_samples, n_params, max_deg+1))
+    univ_vals = np.zeros((n_samples, n_params, max_deg+1), dtype=np.float16)
 
     for i in range(n_params):
 
diff --git a/src/bayesvalidrox/surrogate_models/reg_fast_ard.py b/src/bayesvalidrox/surrogate_models/reg_fast_ard.py
index 2d9eaca8a3b0880c4cfc88a9d32316e03a208c95..1b1582860d60dd3a756e5e2ec91c213601409131 100755
--- a/src/bayesvalidrox/surrogate_models/reg_fast_ard.py
+++ b/src/bayesvalidrox/surrogate_models/reg_fast_ard.py
@@ -384,7 +384,7 @@ class RegressionFastARD(LinearModel, RegressorMixin):
             if self.normalize:
                 x = (X - self._x_mean_) / self._x_std
             var_hat = 1./self.alpha_
-            var_hat += np.sum(np.dot(x[:, self.active_], self.sigma_) *
+            var_hat += np.sum(x[:, self.active_].dot(self.sigma_) *
                               x[:, self.active_], axis=1)
             std_hat = np.sqrt(var_hat)
             return y_hat, std_hat
diff --git a/src/bayesvalidrox/surrogate_models/surrogate_models.py b/src/bayesvalidrox/surrogate_models/surrogate_models.py
index 039a5d68e783b402c53b7ebd51d11c2629f2953b..f519f39a965322cf360a6b27173ca5d77fa1089e 100644
--- a/src/bayesvalidrox/surrogate_models/surrogate_models.py
+++ b/src/bayesvalidrox/surrogate_models/surrogate_models.py
@@ -131,7 +131,7 @@ class MetaModel:
         if self.ExpDesign.method == 'sequential':
             from .sequential_design import SeqDesign
             seq_design = SeqDesign(self)
-            metamodel = seq_design.train_seq_design(Model)
+            metamodel = seq_design.train_seq_design()
 
         elif self.ExpDesign.method == 'normal':
             self.ExpDesignFlag = 'normal'
@@ -488,14 +488,14 @@ class MetaModel:
                              " match!!")
 
         # Preallocate the Psi matrix for performance
-        psi = np.ones((n_samples, n_terms))
+        psi = np.ones((n_samples, n_terms), dtype=np.float16)
         # Assemble the Psi matrix
         for m in range(basis_indices.shape[1]):
             aa = np.where(basis_indices[:, m] > 0)[0]
             try:
                 basisIdx = basis_indices[aa, m]
-                bb = np.reshape(univ_p_val[:, m, basisIdx], psi[:, aa].shape)
-                psi[:, aa] = np.multiply(psi[:, aa], bb)
+                bb = univ_p_val[:, m, basisIdx].reshape(psi[:, aa].shape)
+                psi[:, aa] *= bb
             except ValueError as err:
                 raise err
 
@@ -569,7 +569,6 @@ class MetaModel:
                                                  n_iter=1000, tol=1e-7)
 
             elif reg_method.lower() == 'lars':
-                print("lars")
                 clf_poly = lm.LassoLarsCV(fit_intercept=False)
 
             elif reg_method.lower() == 'sgdr':
@@ -1106,17 +1105,9 @@ class MetaModel:
             model_dict = self.coeffs_dict
 
         if samples is None:
-            if nsamples is None:
-                self.n_samples = 100000
-            else:
-                self.n_samples = nsamples
-
-            samples = self.ExpDesign.generate_samples(self.n_samples,
+            # Generate
+            samples = self.ExpDesign.generate_samples(nsamples,
                                                       sampling_method)
-        else:
-            self.samples = samples
-            self.n_samples = len(samples)
-
         # Transform samples
         samples = self.ExpDesign.transform(samples)
 
@@ -1151,7 +1142,6 @@ class MetaModel:
                         # with error bar
                         clf_poly = self.clf_poly[ouput][in_key]
                         y_mean, y_std = clf_poly.predict(psi, return_std=True)
-
                     except:
                         # without error bar
                         coeffs = self.coeffs_dict[ouput][in_key]