diff --git a/BayesValidRox/surrogate_models/RegressionFastARD.py b/BayesValidRox/surrogate_models/RegressionFastARD.py
index 2f40199e11f92f499d87453796a7b430f295f17b..6aa4de770abe258a2b98318903755bff7caa1f19 100755
--- a/BayesValidRox/surrogate_models/RegressionFastARD.py
+++ b/BayesValidRox/surrogate_models/RegressionFastARD.py
@@ -281,7 +281,7 @@ class RegressionFastARD(LinearModel,RegressorMixin):
                 break
             beta    = n_samples - np.sum(active) + np.sum(Aa * Sdiag )
             beta   /= ( rss + np.finfo(np.float32).eps )
-            
+
             # update precision parameters of coefficients
             A,converged  = update_precisions(Q,S,q,s,A,active,self.tol,
                                              n_samples,False)
diff --git a/BayesValidRox/surrogate_models/aPoly_Construction.py b/BayesValidRox/surrogate_models/aPoly_Construction.py
index 1ad0179bdfa4b11eca45a8c1b680c19bde0fcf82..3fb123acf7055a12b3cc5e82a26e5fdf2446de56 100644
--- a/BayesValidRox/surrogate_models/aPoly_Construction.py
+++ b/BayesValidRox/surrogate_models/aPoly_Construction.py
@@ -1,29 +1,41 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
-"""
-% ========================================================================
-% Construction of Data-driven Orthonormal Polynomial Basis
-% Author: Dr.-Ing. habil. Sergey Oladyshkin
-% Department of Stochastic Simulation and Safety Research for Hydrosystems
-% Institute for Modelling Hydraulic and Environmental Systems
-% Universitaet Stuttgart, Pfaffenwaldring 5a, 70569 Stuttgart
-% E-mail: Sergey.Oladyshkin@iws.uni-stuttgart.de
-% http://www.iws-ls3.uni-stuttgart.de
-% The current program is using definition of arbitrary polynomial chaos expansion (aPC) which is presented in the following manuscript:
-% Oladyshkin, S. and W. Nowak. Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion. Reliability Engineering & System Safety, Elsevier, V. 106, P. 179�190, 2012. DOI: 10.1016/j.ress.2012.05.002.
-% ========================================================================
-"""
 import numpy as np
 import os
 
-def aPoly_Construction(Data, Degree,parIdx):
-    """
-      Input parameters:
-          % Data - raw data array
-          % Degree - higher degree of the orthonormal polynomial basis
-          % Poly_FileName - file name in which the orthonormal polynomial basis should be storred
+def aPoly_Construction(Data, Degree, parIdx, save=False):
     """
+    Construction of Data-driven Orthonormal Polynomial Basis
+    Author: Dr.-Ing. habil. Sergey Oladyshkin
+    Department of Stochastic Simulation and Safety Research for Hydrosystems
+    Institute for Modelling Hydraulic and Environmental Systems
+    Universitaet Stuttgart, Pfaffenwaldring 5a, 70569 Stuttgart
+    E-mail: Sergey.Oladyshkin@iws.uni-stuttgart.de
+    http://www.iws-ls3.uni-stuttgart.de
+    The current script is using definition of arbitrary polynomial chaos expansion (aPC) 
+    which is presented in the following manuscript:
+    Oladyshkin, S. and W. Nowak. Data-driven uncertainty quantification using the arbitrary 
+    polynomial chaos expansion. Reliability Engineering & System Safety, Elsevier, V. 106, P. 
+    179�190, 2012. DOI: 10.1016/j.ress.2012.05.002.
+
+    Parameters
+    ----------
+    Data : array
+        Raw data.
+    Degree : int
+        Maximum polynomial degree.
+    parIdx : TYPE
+        DESCRIPTION.
+    save : TYPE, optional
+        DESCRIPTION. The default is False.
+
+    Returns
+    -------
+    Polynomial : array
+        The coefficients of the orthonormal polynomials.
 
+    """
+    
     # Initialization
     d=Degree #Degree of polinomial expansion
     dd=d+1 #Degree of polinomial for roots defenition
@@ -115,7 +127,8 @@ def aPoly_Construction(Data, Degree,parIdx):
     #     Polynomial[:,k] = Polynomial[:,k]/(MeanOfData**(k))
     
     # Save the Polynomial matrix
-    os.makedirs('polyCoeffs',exist_ok=True)
-    np.save('polyCoeffs/par_{}.npy'.format(parIdx+1), Polynomial)
+    if save:
+        os.makedirs('polyCoeffs',exist_ok=True)
+        np.save('polyCoeffs/par_{}.npy'.format(parIdx+1), Polynomial)
     
     return Polynomial
\ No newline at end of file
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index dcf927cff59de44d4d733205bf167733310794ef..c265d99a0069690002cdc9cacffafa1c904927ed 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -794,8 +794,6 @@ class Metamodel:
         self.LCerror = self.AutoVivification()
             
         if self.ExpDesignFlag != 'sequential':
-            self.allBasisIndices = self.AutoVivification()
-            
             # Read observations or MCReference
             if Model.MeasurementFile is not None or len(Model.Observations) != 0:
                 self.Observations = Model.read_Observation()
@@ -809,15 +807,15 @@ class Metamodel:
         degNew = range(1,maxDeg+1)[np.argmin(abs(M_uptoMax(maxDeg)-ndim*nSamples*d))]
         degNew = degNew if self.ExpDesignFlag == 'sequential' else self.MaxPceDegree
         self.q = np.array(self.q) if not np.isscalar(self.q) else np.array([self.q])
-        if degNew > self.MinPceDegree:
+        if degNew > self.MinPceDegree and self.RegMethod.lower() != 'fastard':
             self.DegreeArray = np.arange(self.MinPceDegree,degNew+1)
         else:
             self.DegreeArray = np.array([degNew])
 
         # Generate all basis indices
+        self.allBasisIndices = self.AutoVivification()
         for deg in self.DegreeArray:
             if deg not in np.fromiter(self.allBasisIndices.keys(), dtype=float):
-                # self.allBasisIndices = self.AutoVivification()
                 # Generate the polynomial basis indices
                 for qidx, q in enumerate(self.q):
                     self.allBasisIndices[str(deg)][str(q)] = self.PolyBasisIndices(degree=deg, q=q)
diff --git a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
index d00381c5181b4b39cc2d69ebf6c05b7d27c30c15..cb3b4221d44ede2844a47cebfabd031686f60bb7 100755
--- a/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
+++ b/BayesValidRox/tests/AnalyticalFunction/Test_AnalyticalFunction.py
@@ -88,7 +88,7 @@ if __name__ == "__main__":
     MetaModelOpts = Metamodel(Inputs)
     
     # Select if you want to preserve the spatial/temporal depencencies
-    MetaModelOpts.DimRedMethod = 'PCA'
+    # MetaModelOpts.DimRedMethod = 'PCA'
     # MetaModelOpts.varPCAThreshold = 100 #99.999
     
     # Select your metamodel method
@@ -118,7 +118,7 @@ if __name__ == "__main__":
     MetaModelOpts.q = np.linspace(0.3,0.8,3) if ndim<5 else 0.75
 
     # Print summary of the regression results
-    MetaModelOpts.DisplayFlag = True
+    # MetaModelOpts.DisplayFlag = True
     
     # ------------------------------------------------
     # ------ Experimental Design Configuration -------
@@ -128,8 +128,8 @@ if __name__ == "__main__":
     MetaModelOpts.addExpDesign()
     
     # One-shot (normal) or Sequential Adaptive (sequential) Design
-    MetaModelOpts.ExpDesign.Method = 'normal'
-    NrofInitSamples = 350 #5*ndim #50
+    MetaModelOpts.ExpDesign.Method = 'sequential'
+    NrofInitSamples = 5*ndim #50
     MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
     
     # Sampling methods
@@ -229,7 +229,7 @@ if __name__ == "__main__":
     # Compute and print RMSE error
     PostPCE.accuracyCheckMetaModel(nSamples=200)
     
-    # sys.exit('STOP!!')
+    sys.exit('STOP!!')
     #=====================================================
     #========  Bayesian inference with Emulator ==========
     #=====================================================
@@ -248,8 +248,7 @@ if __name__ == "__main__":
     BayesOpts.MCMCnwalkers = 300
     # BayesOpts.MCMCinitSamples = np.zeros((1,10)) + 1e-3 * np.random.randn(200, 10)
     # BayesOpts.MCMCnSteps = 1000
-    import emcee
-    BayesOpts.MCMCmoves = emcee.moves.KDEMove()
+    # import emcee
     # BayesOpts.MCMCmoves = [(emcee.moves.DEMove(), 0.2),
     #                         (emcee.moves.WalkMove(), 0.3),
     #                         (emcee.moves.KDEMove(), 0.3),
diff --git a/BayesValidRox/tests/PollutionTest/Pollution_Test.py b/BayesValidRox/tests/PollutionTest/Pollution_Test.py
index b5c7a458b50465e8a49c0c1242a4a05902e999d6..888d43121ccf6b626991a54b40d835e9c2ff76eb 100644
--- a/BayesValidRox/tests/PollutionTest/Pollution_Test.py
+++ b/BayesValidRox/tests/PollutionTest/Pollution_Test.py
@@ -106,7 +106,7 @@ if __name__ == "__main__":
     MetaModelOpts = Metamodel(Inputs)
     
     # Select if you want to preserve the spatial/temporal depencencies
-    MetaModelOpts.DimRedMethod = 'PCA'
+    # MetaModelOpts.DimRedMethod = 'PCA'
     # MetaModelOpts.varPCAThreshold = 99.5 #99.999
     
      # Select your metamodel method
@@ -131,13 +131,13 @@ if __name__ == "__main__":
     # error (or the highest score=1-LOO)estimator is chosen as the final 
     # metamodel.
     MetaModelOpts.MinPceDegree = 1
-    MetaModelOpts.MaxPceDegree = 8 #12
+    MetaModelOpts.MaxPceDegree = 10 #8
     
     # q-quasi-norm 0<q<1 (default=1)
-    # MetaModelOpts.q = 0.75 #np.linspace(0.3,0.8, 3)
+    MetaModelOpts.q = 0.95 #np.linspace(0.3,0.8, 3)
     
     # Print summary of the regression results
-    MetaModelOpts.DisplayFlag = True
+    # MetaModelOpts.DisplayFlag = True
     
     # ------ Experimental Design --------
     # Generate an experimental design of size NrExpDesign based on a latin 
@@ -145,8 +145,8 @@ if __name__ == "__main__":
     MetaModelOpts.addExpDesign()
     
     # One-shot (normal) or Sequential Adaptive (sequential) Design
-    MetaModelOpts.ExpDesign.Method = 'normal'
-    NrofInitSamples = 100 #200
+    MetaModelOpts.ExpDesign.Method = 'sequential'
+    NrofInitSamples = 20 #200
     MetaModelOpts.ExpDesign.NrSamples = NrofInitSamples
     # Sampling methods
     # 1) random 2) latin_hypercube 3) sobol 4) halton 5) hammersley 6) chebyshev(FT) 
@@ -264,13 +264,13 @@ if __name__ == "__main__":
     
     # Select the inference method
     BayesOpts.SamplingMethod = 'MCMC'
-    BayesOpts.MCMCnwalkers = 100
+    BayesOpts.MCMCnwalkers = 300
     # BayesOpts.MCMCnSteps = 500
-    import emcee
-    BayesOpts.MCMCmoves = [(emcee.moves.DEMove(), 0.4),
-                           (emcee.moves.WalkMove(), 0.1),
-                           (emcee.moves.KDEMove(), 0.1),
-                           (emcee.moves.DESnookerMove(), 0.4),]
+    # import emcee
+    # BayesOpts.MCMCmoves = [(emcee.moves.KDEMove(), 1.0)]
+    #                        (emcee.moves.WalkMove(), 0.1),
+    #                        (emcee.moves.DEMove(), 0.1),
+    #                        (emcee.moves.DESnookerMove(), 0.4),]
     
     BayesOpts.PlotPostDist = True
     BayesOpts.PlotPostPred = True