From 631f10afb9ccb2b446c4ca1e0fd6cfbf4f572771 Mon Sep 17 00:00:00 2001
From: kohlhaasrebecca <rebecca.kohlhaas@outlook.com>
Date: Wed, 9 Oct 2024 10:50:02 +0200
Subject: [PATCH] [Fix] Remove parameters MetaModel, Model, ExpDesign from
 PostProcessing

PostProcessing uses engine.* instead of self.* for MetaModel, Model and ExpDesign.
In addition it calls the engine properties instead of Model properties where possible
---
 .../post_processing/post_processing.py        | 182 ++++++++++++------
 1 file changed, 125 insertions(+), 57 deletions(-)

diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index c7b8bf98e..225c32153 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -69,9 +69,11 @@ class PostProcessing:
 
         bar_plot = True if plot_type == 'bar' else False
         meta_model_type = self.engine.MetaModel.meta_model_type
+        meta_model_type = self.engine.MetaModel.meta_model_type
 
         # Read Monte-Carlo reference
         self.mc_reference = self.engine.Model.read_observation('mc_ref')
+        self.mc_reference = self.engine.Model.read_observation('mc_ref')
 
         # Set the x values
         x_values_orig = self.engine.ExpDesign.x_values
@@ -79,8 +81,14 @@ class PostProcessing:
         # Compute the moments with the PCEModel object
         self.means, self.stds = self.engine.MetaModel.compute_moments()
 
+        # Open a pdf for the plots
+        newpath = (f'Outputs_PostProcessing_{self.name}/')
+        if not os.path.exists(newpath):
+            os.makedirs(newpath)
+
         # Plot the best fit line, set the linewidth (lw), color and
         # transparency (alpha) of the line
+        for key in self.engine.out_names:
         for key in self.engine.out_names:
             fig, ax = plt.subplots(nrows=1, ncols=2)
 
@@ -180,11 +188,14 @@ class PostProcessing:
         else:
             self.model_out_dict = self._eval_model(samples, key_str='valid')
         self.pce_out_mean, self.pce_out_std = self.engine.eval_metamodel(samples)
+        self.pce_out_mean, self.pce_out_std = self.engine.eval_metamodel(samples)
 
         try:
             key = self.engine.out_names[1]
+            key = self.engine.out_names[1]
         except IndexError:
             key = self.engine.out_names[0]
+            key = self.engine.out_names[0]
 
         n_obs = self.model_out_dict[key].shape[1]
 
@@ -217,6 +228,13 @@ class PostProcessing:
         Exception
             When neither n_samples nor samples are provided.
 
+        Returns
+        -------
+        rmse: dict
+            Root mean squared error for each output.
+        valid_error : dict
+            Validation error for each output.
+
         """
         # Set the number of samples
         if n_samples:
@@ -236,10 +254,12 @@ class PostProcessing:
 
         # Run the PCE model with the generated samples
         pce_outputs, _ = self.engine.eval_metamodel(samples)
+        pce_outputs, _ = self.engine.eval_metamodel(samples)
 
         self.rmse = {}
         self.valid_error = {}
         # Loop over the keys and compute RMSE error.
+        for key in self.engine.out_names:
         for key in self.engine.out_names:
             # Root mena square
             self.rmse[key] = mean_squared_error(outputs[key], pce_outputs[key],
@@ -259,6 +279,8 @@ class PostProcessing:
         # Save error dicts in PCEModel object
         self.engine.MetaModel.rmse = self.rmse
         self.engine.MetaModel.valid_error = self.valid_error
+        self.engine.MetaModel.rmse = self.rmse
+        self.engine.MetaModel.valid_error = self.valid_error
 
 
     # -------------------------------------------------------------------------
@@ -549,9 +571,9 @@ class PostProcessing:
         # This function currently only supports PCE/aPCE
         PCEModel = self.engine.MetaModel
         if not hasattr(PCEModel, 'meta_model_type'):
-            raise AttributeError('Sobol indices only support PCE-type models!')
+            raise AttributeError('Sobol indices currently only support PCE-type models!')
         if PCEModel.meta_model_type.lower() not in ['pce', 'apce']:
-            raise AttributeError('Sobol indices only support PCE-type models!')
+            raise AttributeError('Sobol indices currently only support PCE-type models!')
             
         # Extract the necessary variables
         basis_dict = PCEModel._basis_dict
@@ -566,6 +588,7 @@ class PostProcessing:
 
             sobol_cell_, total_sobol_ = {}, {}
 
+            for output in self.engine.out_names:
             for output in self.engine.out_names:
 
                 n_meas_points = len(coeffs_dict[f'b_{b_i+1}'][output])
@@ -754,6 +777,7 @@ class PostProcessing:
                 total_sobol_all[k][i] = v
 
         self.total_sobol = {}
+        for output in self.engine.out_names:
         for output in self.engine.out_names:
             self.total_sobol[output] = np.mean(total_sobol_all[output], axis=0)
 
@@ -763,6 +787,7 @@ class PostProcessing:
 
         fig = plt.figure()
 
+        for outIdx, output in enumerate(self.engine.out_names):
         for outIdx, output in enumerate(self.engine.out_names):
 
             # Extract total Sobol indices
@@ -846,12 +871,14 @@ class PostProcessing:
             self.n_samples = samples.shape[0]
 
         # Evaluate the original and the surrogate model
-        if outputs is None:
-            y_val = self._eval_model(samples, key_str='valid')
-        else: 
-            y_val = outputs
+        y_val = self._eval_model(samples, key_str='valid')
         y_pce_val, _ = self.engine.eval_metamodel(samples=samples)
 
+        # Open a pdf for the plots
+        newpath = f'Outputs_PostProcessing_{self.name}/'
+        if not os.path.exists(newpath):
+            os.makedirs(newpath)
+
         # Fit the data(train the model)
         for key in y_pce_val.keys():
 
@@ -940,55 +967,93 @@ class PostProcessing:
             plt.close()
 
     # -------------------------------------------------------------------------
-    def plot_metamodel_3d(self, n_samples = 10):
-        """
-        Visualize the results of a PCE MetaModel as a 3D surface over two input 
-        parameters.
-        
-        Parameters
-        ----------
-        n_samples : int
-            Number of samples that are used to generate the 3D plot.
+    def eval_pce_model_3d(self):
 
-        Raises
-        ------
-        AttributeError
-            This function is only applicable if the MetaModel input dimension is 2.
+        self.n_samples = 1000
 
-        Returns
-        -------
-        None.
+        PCEModel = self.engine.MetaModel
+        n_samples = self.n_samples
 
-        """
-        if self.engine.ExpDesign.ndim !=2:
-            raise AttributeError('This function is only applicable if the MetaModel input dimension is 2.')
-        samples = self.engine.ExpDesign.generate_samples(n_samples)
-        samples = np.sort(np.sort(samples, axis = 1), axis = 0)
-        mean, stdev = self.engine.eval_metamodel(samples = samples)
-        
-        if self.engine.emulator:
-            title = 'MetaModel'
-        else:
-            title = 'Model'
-        X,Y = np.meshgrid(samples[:,0],samples[:,1])
-        for name in self.engine.out_names:
-            for t in range(mean[name].shape[1]):
-                fig = plt.figure()
-                ax = plt.axes(projection='3d')
-                ax.plot_surface(X, Y, np.atleast_2d(mean[name][:,t]), rstride=1, cstride=1,
-                                cmap='viridis', edgecolor='none')
-                ax.set_title(title)
-                ax.set_xlabel('$x_1$')
-                ax.set_ylabel('$x_2$')
-                ax.set_zlabel('$f(x_1,x_2)$')
-        
-                plt.grid()
-        
-                # save the figure to file
-                fig.savefig(f'./{self.out_dir}/3DPlot_{title}_{name}{t}.pdf',
-                                bbox_inches='tight')
-                plt.close(fig)
+        # Create 3D-Grid
+        # TODO: Make it general
+        x = np.linspace(-5, 10, n_samples)
+        y = np.linspace(0, 15, n_samples)
+
+        X, Y = np.meshgrid(x, y)
+        PCE_Z = np.zeros((self.n_samples, self.n_samples))
+        Model_Z = np.zeros((self.n_samples, self.n_samples))
+
+        for idxMesh in range(self.n_samples):
+            sample_mesh = np.vstack((X[:, idxMesh], Y[:, idxMesh])).T
+
+            univ_p_val = PCEModel.univ_basis_vals(sample_mesh)
+
+            for Outkey, ValuesDict in PCEModel.coeffs_dict.items():
+
+                pce_out_mean = np.zeros((len(sample_mesh), len(ValuesDict)))
+                pce_out_std = np.zeros((len(sample_mesh), len(ValuesDict)))
+                model_outs = np.zeros((len(sample_mesh), len(ValuesDict)))
+
+                for Inkey, InIdxValues in ValuesDict.items():
+                    idx = int(Inkey.split('_')[1]) - 1
+                    basis_deg_ind = PCEModel.basis_dict[Outkey][Inkey]
+                    clf_poly = PCEModel.clf_poly[Outkey][Inkey]
+
+                    PSI_Val = PCEModel.create_psi(basis_deg_ind, univ_p_val)
+
+                    # Perdiction with error bar
+                    y_mean, y_std = clf_poly.predict(PSI_Val, return_std=True)
+
+                    pce_out_mean[:, idx] = y_mean
+                    pce_out_std[:, idx] = y_std
+
+                    # Model evaluation
+                    model_out_dict, _ = self.engine.Model.run_model_parallel(sample_mesh,
+                                                                 key_str='Valid3D')
+                    model_outs[:, idx] = model_out_dict[Outkey].T
+
+                PCE_Z[:, idxMesh] = y_mean
+                Model_Z[:, idxMesh] = model_outs[:, 0]
+
+        # ---------------- 3D plot for PCEModel -----------------------
+        fig_PCE = plt.figure()
+        ax = plt.axes(projection='3d')
+        ax.plot_surface(X, Y, PCE_Z, rstride=1, cstride=1,
+                        cmap='viridis', edgecolor='none')
+        ax.set_title('PCEModel')
+        ax.set_xlabel('$x_1$')
+        ax.set_ylabel('$x_2$')
+        ax.set_zlabel('$f(x_1,x_2)$')
+
+        plt.grid()
+
+        #  Saving the figure
+        newpath = f'Outputs_PostProcessing_{self.name}/'
+        if not os.path.exists(newpath):
+            os.makedirs(newpath)
+
+        # save the figure to file
+        fig_PCE.savefig(f'./{newpath}/3DPlot_PCEModel.pdf',
+                        bbox_inches='tight')
+        plt.close(fig_PCE)
+
+        # ---------------- 3D plot for Model -----------------------
+        fig_Model = plt.figure()
+        ax = plt.axes(projection='3d')
+        ax.plot_surface(X, Y, PCE_Z, rstride=1, cstride=1,
+                        cmap='viridis', edgecolor='none')
+        ax.set_title('Model')
+        ax.set_xlabel('$x_1$')
+        ax.set_ylabel('$x_2$')
+        ax.set_zlabel('$f(x_1,x_2)$')
+        plt.grid()
         
+        # Save the figure
+        fig_Model.savefig(f'./{newpath}/3DPlot_Model.pdf',
+                          bbox_inches='tight')
+        plt.close(fig_Model)
+
+        return
 
     # -------------------------------------------------------------------------
     def _get_sample(self, n_samples=None):
@@ -1003,6 +1068,7 @@ class PostProcessing:
         """
         if n_samples is None:
             n_samples = self.n_samples
+        self.samples = self.engine.ExpDesign.generate_samples(
         self.samples = self.engine.ExpDesign.generate_samples(
             n_samples,
             sampling_method='random')
@@ -1034,6 +1100,7 @@ class PostProcessing:
             self.n_samples = len(samples)
 
         model_outs, _ = self.engine.Model.run_model_parallel(samples, key_str=key_str)
+        model_outs, _ = self.engine.Model.run_model_parallel(samples, key_str=key_str)
 
         return model_outs
 
@@ -1129,12 +1196,9 @@ class PostProcessing:
         None.
 
         """
-        # This function currently only supports PCE/aPCE
-        PCEModel = self.engine.MetaModel
-        if not hasattr(PCEModel, 'meta_model_type'):
-            raise AttributeError('This evaluation only support PCE-type models!')
-        if PCEModel.meta_model_type.lower() not in ['pce', 'apce']:
-            raise AttributeError('This evaluation only support PCE-type models!')
+        newpath = f'Outputs_PostProcessing_{self.name}/'
+        if not os.path.exists(newpath):
+            os.makedirs(newpath)
 
         # List of markers and colors
         color = cycle((['b', 'g', 'r', 'y', 'k']))
@@ -1142,6 +1206,7 @@ class PostProcessing:
 
         fig = plt.figure()
         # Plot the model vs PCE model
+        for keyIdx, key in enumerate(self.engine.out_names):
         for keyIdx, key in enumerate(self.engine.out_names):
 
             y_pce_val = self.pce_out_mean[key]
@@ -1186,3 +1251,6 @@ class PostProcessing:
 
             # Destroy the current plot
             plt.close()
+
+        # Zip the subdirectories
+        self.engine.Model.zip_subdirs(f'{self.engine.Model.name}valid', f'{self.engine.Model.name}valid_')
-- 
GitLab