diff --git a/BayesValidRox/tests/PA-A/Benchmark_PAA.py b/BayesValidRox/tests/PA-A/Benchmark_PAA.py
index f32cf5ae456dcdcf676737722a4db15675457980..9f37c18420c57338368345cf0720aefbc92d9312 100755
--- a/BayesValidRox/tests/PA-A/Benchmark_PAA.py
+++ b/BayesValidRox/tests/PA-A/Benchmark_PAA.py
@@ -36,7 +36,7 @@ plt.rc('figure', titlesize=SIZE)  # fontsize of the figure title
 try:
     import cPickle as pickle
 except ModuleNotFoundError:
-    import pickle
+    import _pickle as pickle
 
 # Add BayesValidRox path
 sys.path.insert(0,'./../../../BayesValidRox/')
@@ -334,7 +334,7 @@ if __name__ == "__main__":
 
     # Load the objects
     # with open(result_folder+'/outputs_ffpm-stokesdarcyBJ/PA_A_Bayesffpm-stokesdarcyBJ-valid.pkl', 'rb') as input:
-    #     BayesValid_BJ = pickle.load(input)
+    #     BayesValid_BJ = joblib.load(input)
 
     # Model stokesdarcy-ER
     PCEModel_ER, BayesCalib_ER, BayesValid_ER = stokesdarcy.run(paramsAvg, couplingcond='ER',
@@ -343,14 +343,14 @@ if __name__ == "__main__":
 
     # Load the objects
     # with open(result_folder+'/outputs_ffpm-stokesdarcyER/PA_A_Bayesffpm-stokesdarcyER-valid.pkl', 'rb') as input:
-    #     BayesValid_ER = pickle.load(input)
+    #     BayesValid_ER = joblib.load(input)
 
     # Model stokesdarcy-pnm with the averaged data
     PCEModel_PNM, BayesCalib_PNM, BayesValid_PNM = stokespnm.run(paramsAvg, PCEEDMethod=PCEExpDesignMethod)
 
     # Load the objects
     # with open(result_folder+'/outputs_ffpm-stokespnm/PA_A_Bayesffpm-stokespnm-valid.pkl', 'rb') as input:
-    #     BayesValid_PNM = pickle.load(input)
+    #     BayesValid_PNM = joblib.load(input)
     
     # Model stokesdarcy-pnm without the averaged data
     PCEModel_PNM_NA, BayesCalib_PNM_NA, BayesValid_PNM_NA = stokespnm.run(params, averaging=False,
diff --git a/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py b/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py
index 1da3e35ff0c5d292a3983ff8c7323d30d78f49d1..e90ef07079fa9b051d6c91c2576d16f40d963c3b 100755
--- a/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py
+++ b/BayesValidRox/tests/PA-A/ffpm_validation_stokesdarcy.py
@@ -5,7 +5,7 @@ Created on Thu Aug 13 09:53:11 2020
 
 @author: farid
 """
-import sys, os
+import sys, os, joblib
 import numpy as np
 import scipy.stats as stats
 import shutil
@@ -13,7 +13,7 @@ import pandas as pd
 try:
     import cPickle as pickle
 except ModuleNotFoundError:
-    import pickle
+    import _pickle as pickle
 
 import matplotlib
 matplotlib.use('agg')
@@ -255,11 +255,11 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
 
     ## Save PCE models
     with open('PCEModel_'+Model.Name+'.pkl', 'wb') as output:
-        pickle.dump(PCEModel, output, pickle.HIGHEST_PROTOCOL)
+        joblib.dump(PCEModel, output, 2)
 
     # Load the objects
     # with open('PCEModel_'+Model.Name+'.pkl', 'rb') as input:
-    #     PCEModel = pickle.load(input)
+    #     PCEModel = joblib.load(input)
 
     #=====================================================
     #=========  POST PROCESSING OF METAMODELS  ===========
@@ -351,11 +351,11 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
 
     #Save class objects
     with open('PA_A_Bayes'+Model.Name+'.pkl', 'wb') as output:
-        pickle.dump(BayesCalib, output, pickle.HIGHEST_PROTOCOL)
+        joblib.dump(BayesCalib, output, 2)
 
     # Load the objects
     # with open('PA_A_Bayes'+Model.Name+'.pkl', 'rb') as input:
-    #     BayesCalib = pickle.load(input)
+    #     BayesCalib = joblib.load(input)
     #=====================================================
     #=================  Visualization  ===================
     #=====================================================
@@ -428,11 +428,11 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
 
     # Save PCE models
     with open('PCEModel_'+ValidModel.Name+'.pkl', 'wb') as output:
-        pickle.dump(ValidPCEModel, output, pickle.HIGHEST_PROTOCOL)
+        joblib.dump(ValidPCEModel, output, 2)
 
     # Load the objects
     # with open('PCEModel_'+ValidModel.Name+'.pkl', 'rb') as input:
-    #     ValidPCEModel = pickle.load(input)
+    #     ValidPCEModel = joblib.load(input)
 
     #=====================================================
     #=========  POST PROCESSING OF METAMODELS  ===========
@@ -521,10 +521,8 @@ def run(params, errorPerc=0.05, couplingcond='BJ', PCEEDMethod='normal'):
     #==============  Save class objects  =================
     #=====================================================
     with open('PA_A_Bayes'+ValidModel.Name+'.pkl', 'wb') as output:
-        pickle.dump(BayesValid, output, pickle.HIGHEST_PROTOCOL)
-    # p = pickle.Pickler(open('PA_A_Bayes'+ValidModel.Name+'.pkl',"wb"))
-    # p.fast = True
-    # p.dump(BayesValid)
+        pickle.joblib(BayesValid, output, 2)
+
     #=====================================================
     #=================  Visualization  ===================
     #=====================================================
diff --git a/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py b/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py
index cd66e015c57bc52d427457f8883d7c3a6dd87d03..bd79bde60949fa7d775ba3013ab6c5e4a23c20bc 100755
--- a/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py
+++ b/BayesValidRox/tests/PA-A/ffpm_validation_stokespnm.py
@@ -5,7 +5,7 @@ Created on Thu Aug 13 09:53:11 2020
 
 @author: farid
 """
-import sys, os
+import sys, os, joblib
 import numpy as np
 import scipy.stats as stats
 import pandas as pd
@@ -13,8 +13,7 @@ import shutil
 try:
     import cPickle as pickle
 except ModuleNotFoundError:
-    import pickle
-
+    import _pickle as pickle
 
 import matplotlib
 matplotlib.use('agg')
@@ -268,11 +267,11 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
 
     # Save PCE models
     with open('PCEModel_'+Model.Name+'.pkl', 'wb') as output:
-        pickle.dump(PCEModel, output, pickle.HIGHEST_PROTOCOL)
+        joblib.dump(PCEModel, output, 2)
 
     # Load the objects
     # with open('PCEModel_'+Model.Name+'.pkl', 'rb') as input:
-    #     PCEModel = pickle.load(input)
+    #     PCEModel = joblib.load(input)
 
     #=====================================================
     #=========  POST PROCESSING OF METAMODELS  ===========
@@ -356,15 +355,15 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
     # BayesOptsCalib.MeasurementError = np.sqrt(calibMeasuError)
 
     # ----- Strat Bayesian inference -------
-    # BayesCalib = BayesOptsCalib.create_Inference()
+    BayesCalib = BayesOptsCalib.create_Inference()
 
-    # # Save class objects
-    # with open('PA_A_Bayes'+Model.Name+'.pkl', 'wb') as output:
-    #     pickle.dump(BayesCalib, output, pickle.HIGHEST_PROTOCOL)
+    # Save class objects
+    with open('PA_A_Bayes'+Model.Name+'.pkl', 'wb') as output:
+        joblib.dump(BayesCalib, output, 2)
 
     # Load the objects
-    with open('PA_A_Bayes'+Model.Name+'.pkl', 'rb') as input:
-        BayesCalib = pickle.load(input)
+    # with open('PA_A_Bayes'+Model.Name+'.pkl', 'rb') as input:
+    #     BayesCalib = joblib.load(input)
 
     #=====================================================
     #=================  Visualization  ===================
@@ -435,11 +434,11 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
 
     # Save PCE models
     with open('PCEModel_'+ValidModel.Name+'.pkl', 'wb') as output:
-        pickle.dump(ValidPCEModel, output, pickle.HIGHEST_PROTOCOL)
+        joblib.dump(ValidPCEModel, output, 2)
 
     # Load the objects
     # with open('PCEModel_'+ValidModel.Name+'.pkl', 'rb') as input:
-    #     ValidPCEModel = pickle.load(input)
+    #     ValidPCEModel = joblib.load(input)
 
     #=====================================================
     #=========  POST PROCESSING OF METAMODELS  ===========
@@ -526,7 +525,7 @@ def run(params, averaging=True,errorPerc=0.05, PCEEDMethod='normal'):
 
     # Save class objects
     with open('PA_A_Bayes'+ValidModel.Name+'.pkl', 'wb') as output:
-        pickle.dump(BayesValid, output, pickle.HIGHEST_PROTOCOL)
+        joblib.dump(BayesValid, output, 2)
 
     #=====================================================
     #=================  Visualization  ===================