From c8e34571a69172803f7b76d64d80357f53395bc1 Mon Sep 17 00:00:00 2001
From: faridm69 <faridmohammadi69@gmail.com>
Date: Thu, 1 Oct 2020 11:41:20 +0200
Subject: [PATCH] [surrogate] fixed bugs.

---
 .../BayesInference/BayesInference.py          |  4 +--
 .../PostProcessing/PostProcessing.py          |  2 +-
 BayesValidRox/surrogate_models/ExpDesigns.py  |  3 +-
 .../surrogate_models/RegressionFastARD.py     |  1 -
 .../surrogate_models/surrogate_models.py      | 31 ++++++++++---------
 5 files changed, 21 insertions(+), 20 deletions(-)

diff --git a/BayesValidRox/BayesInference/BayesInference.py b/BayesValidRox/BayesInference/BayesInference.py
index 0d7fc81c2..718276889 100644
--- a/BayesValidRox/BayesInference/BayesInference.py
+++ b/BayesValidRox/BayesInference/BayesInference.py
@@ -158,7 +158,7 @@ class BayesInference:
         TotalOutputs = np.empty((SampleSize,0))
         for outputType in OutputType:
             TotalOutputs = np.hstack((TotalOutputs, Outputs[outputType]))
-        
+
         # Add the std of the PCE is chosen as emulator.
         if self.emulator:
             covMatrix_PCE = np.zeros((NofMeasurements, NofMeasurements), float)
@@ -178,7 +178,7 @@ class BayesInference:
             logL = multivariate_normal.logpdf(TotalOutputs, mean=Data, cov=covMatrix)
         except:
             logL = -np.inf
-
+            
         return logL
     
     #--------------------------------------------------------------------------------------------------------
diff --git a/BayesValidRox/PostProcessing/PostProcessing.py b/BayesValidRox/PostProcessing/PostProcessing.py
index 334c9ad7d..c50be9a86 100644
--- a/BayesValidRox/PostProcessing/PostProcessing.py
+++ b/BayesValidRox/PostProcessing/PostProcessing.py
@@ -67,7 +67,7 @@ class PostProcessing:
             
             if PCEModel.DimRedMethod.lower() == 'pca':
                 PCA = PCEModel.pca[Outkey]
-                self.PCEMeans[Outkey] = PCA.inverse_transform(PCEMean)
+                self.PCEMeans[Outkey] = PCA.mean_ + np.dot(PCEMean, PCA.components_)
                 self.PCEStd[Outkey] = np.sqrt(np.dot(PCEVar, PCA.components_**2))
             else:
                 self.PCEMeans[Outkey] = PCEMean
diff --git a/BayesValidRox/surrogate_models/ExpDesigns.py b/BayesValidRox/surrogate_models/ExpDesigns.py
index 1a1ed840f..ccb9d5d7e 100644
--- a/BayesValidRox/surrogate_models/ExpDesigns.py
+++ b/BayesValidRox/surrogate_models/ExpDesigns.py
@@ -31,7 +31,8 @@ class ExpDesigns:
         self.NCandidate = 1
         self.PostSnapshot = False
         self.stepSnapshot = 2
-        self.UtilityFunction = 'None'
+        self.ExploitMethod = 'Space-filling'
+        self.UtilityFunction = 'Space-filling'
         self.ExploreMethod = 'random'
         self.NrofCandGroups = 4
         self.nReprications = 1
diff --git a/BayesValidRox/surrogate_models/RegressionFastARD.py b/BayesValidRox/surrogate_models/RegressionFastARD.py
index 880baf488..93819b75a 100755
--- a/BayesValidRox/surrogate_models/RegressionFastARD.py
+++ b/BayesValidRox/surrogate_models/RegressionFastARD.py
@@ -236,7 +236,6 @@ class RegressionFastARD(LinearModel,RegressorMixin):
             # in case of almost perfect multicollinearity between some features
             # start from feature 0
             if np.sum( XXd - X_mean**2 < np.finfo(np.float32).eps ) > 0:#>= 0:
-                print("start from feature 0.")
                 A[0]       = np.finfo(np.float16).eps
                 active[0]  = True
             
diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index 636688190..791b6dff8 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -684,8 +684,9 @@ class Metamodel:
         LCerror = AllLCerror[str(bestIdx)]
         degree = DegreeArray[np.nanargmax(scores)]
         qnorm = float(qnorm[bestqIdx])
-
-        if verbose:
+        
+        # ------------------ Summary of results ------------------
+        if self.DisplayFlag:
             # Create PSI_Sparse by removing redundent terms
             nnz_idx = np.nonzero(coeffs)[0]
             BasisIndices_Sparse = PolynomialDegrees[nnz_idx]
@@ -697,8 +698,8 @@ class Metamodel:
                                                 round(len(BasisIndices_Sparse)/P, 3)))
             print('Final ModLOO error estimate: %.3e'%(1-max(scores)))
             print('\n'+'-'*50)
-        # ------------------ Summary of results ------------------
-        if self.DisplayFlag:
+        
+        if verbose:
             print('='*50)
             print(' '* 10 +' Summary of results ')
             print('='*50)
@@ -867,9 +868,10 @@ class Metamodel:
         
         # Transform via Principal Component Analysis
         from sklearn.decomposition import IncrementalPCA as sklearnPCA
-        # from sklearn.decomposition import TruncatedSVD as TruncatedSVD
+        # from sklearn.decomposition import TruncatedSVD as sklearnPCA
+
         n_samples, n_features = Output.shape
-        covar_matrix = sklearnPCA(n_components=None)#, svd_solver='full')
+        covar_matrix = sklearnPCA(n_components=None) #, svd_solver='full')
         covar_matrix.fit(Output)
         var = np.cumsum(np.round(covar_matrix.explained_variance_ratio_, decimals=5)*100)
         try:
@@ -878,14 +880,14 @@ class Metamodel:
             selected_n_components = min(n_samples, n_features)
         
         nComponents = min(n_samples, n_features, selected_n_components)
-
-        pca = sklearnPCA(n_components=nComponents)#, svd_solver='full')
+        pca = sklearnPCA(n_components=nComponents) #, svd_solver='full')
         OutputMatrix = pca.fit_transform(Output)
 
+        # pca.mean_ = np.zeros(Output.shape[1])
         return pca, OutputMatrix
 
     #--------------------------------------------------------------------------------------------------------
-    def GaussianProcessEmulator(self, X, y, autoSelect=True, varIdx=None):
+    def GaussianProcessEmulator(self, X, y, autoSelect=False, varIdx=None):
         
         from sklearn.gaussian_process import GaussianProcessRegressor
         from sklearn.gaussian_process.kernels import (RBF, Matern, RationalQuadratic,
@@ -914,7 +916,7 @@ class Metamodel:
             gp = gp[np.argmax(BME)]    
         
         else:
-            gp = GaussianProcessRegressor(kernel=kernels[-1], n_restarts_optimizer=3,
+            gp = GaussianProcessRegressor(kernel=kernels[1], n_restarts_optimizer=3,
                                           normalize_y=True)
             gp.fit(X, y)
             
@@ -1074,7 +1076,7 @@ class Metamodel:
         #                 prevBasisIndices = self.BasisDict[key]["y_"+str(idx)] if idx != 0 else None
                             
         #                 outDict = self.adaptiveSparseaPCE(CollocationPoints,OutputMatrix[:,idx],
-        #                                           idx,prevBasisIndices,verbose=True)
+        #                                           idx,prevBasisIndices)
                         
                         
                         
@@ -2183,7 +2185,7 @@ class Metamodel:
                 PCEOutputs_mean[:, idx] = y_mean
                 PCEOutputs_std[:, idx] = y_std
                 idx += 1
-            
+
             # TODO: Plot PCE vs Components
             # PCA = self.pca[Outkey]
             # self.plot_components(PCA, Y[Outkey], PCEOutputs_mean, PCEOutputs_std, plotname='PCA_PCE_'+Outkey)
@@ -2251,13 +2253,12 @@ class Metamodel:
         np.fill_diagonal(covMatrix_PCE, varPCE)
         
         # Aggregate the cov matrices
-        covMatrix = covMatrix + covMatrix_PCE
+        covMatrix += covMatrix_PCE
         
         # Compute likelihood
         self.Likelihoods = stats.multivariate_normal.pdf(TotalOutputs, 
                                                          mean=TotalData, 
-                                                         cov=covMatrix)#,
-                                                         # allow_singular=True)
+                                                         cov=covMatrix)
         
         return self.Likelihoods       
     
-- 
GitLab