From 73cc5aea61050b81fa293594b80c104edc1c08a8 Mon Sep 17 00:00:00 2001
From: kohlhaasrebecca <rebecca.kohlhaas@outlook.com>
Date: Thu, 10 Oct 2024 17:25:44 +0200
Subject: [PATCH] [fix] Sorted out_dir, changed eval_pce_model_3d to
 plot_metamodel_3d

---
 .../post_processing/post_processing.py        | 193 +++++++-----------
 1 file changed, 70 insertions(+), 123 deletions(-)

diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index a111cc150..c7b8bf98e 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -29,14 +29,21 @@ class PostProcessing:
     name : string
         Name of the PostProcessing object to be used for saving the generated files.
         The default is 'calib'.
+    out_dir : string
+        Output directory in which the images are placed. The default is ''.
         
     """
 
-    def __init__(self, engine, name='calib'):
+    def __init__(self, engine, name='calib', out_dir = ''):
         self.engine = engine
         self.name = name
         
+        self.out_dir = f'{out_dir}/Outputs_PostProcessing_{self.name}/'
         
+        # Open a pdf for the plots
+        if not os.path.exists(self.out_dir):
+            os.makedirs(self.out_dir)
+
 
     # -------------------------------------------------------------------------
     def plot_moments(self, xlabel:str='Time [s]', plot_type:str=None):
@@ -72,11 +79,6 @@ 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:
@@ -138,7 +140,7 @@ class PostProcessing:
 
             # save the current figure
             fig.savefig(
-                f'./{newpath}Mean_Std_PCE_{key}.pdf',
+                f'./{self.out_dir}Mean_Std_PCE_{key}.pdf',
                 bbox_inches='tight'
                 )
 
@@ -190,6 +192,9 @@ class PostProcessing:
             self._plot_validation()
         else:
             self._plot_validation_multi(x_values=x_values, x_axis=x_axis)
+        
+            # Zip the subdirectories
+            self.engine.Model.zip_subdirs(f'{self.engine.Model.name}valid', f'{self.engine.Model.name}valid_')
 
     # -------------------------------------------------------------------------
     def check_accuracy(self, n_samples=None, samples=None, outputs=None)-> None:
@@ -272,6 +277,7 @@ class PostProcessing:
         n_init_samples = engine.ExpDesign.n_init_samples
         n_total_samples = engine.ExpDesign.X.shape[0]
 
+        # TODO: move this one!
         newpath = f'Outputs_PostProcessing_{self.name}/seq_design_diagnostics/'
         if not os.path.exists(newpath):
             os.makedirs(newpath)
@@ -755,10 +761,6 @@ class PostProcessing:
         par_names = self.engine.ExpDesign.par_names
         x_values_orig = self.engine.ExpDesign.x_values
 
-        newpath = (f'Outputs_PostProcessing_{self.name}/')
-        if not os.path.exists(newpath):
-            os.makedirs(newpath)
-
         fig = plt.figure()
 
         for outIdx, output in enumerate(self.engine.out_names):
@@ -801,14 +803,14 @@ class PostProcessing:
                 plt.legend(loc='best', frameon=True)
 
             # Save indices
-            np.savetxt(f'./{newpath}totalsobol_' +
+            np.savetxt(f'./{self.out_dir}totalsobol_' +
                        output.replace('/', '_') + '.csv',
                        total_sobol.T, delimiter=',',
                        header=','.join(par_names), comments='')
 
             # save the current figure
             fig.savefig(
-                f'./{newpath}Sobol_indices_{output}.pdf',
+                f'./{self.out_dir}Sobol_indices_{output}.pdf',
                 bbox_inches='tight'
                 )
 
@@ -850,11 +852,6 @@ class PostProcessing:
             y_val = outputs
         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():
 
@@ -877,7 +874,7 @@ class PostProcessing:
                 plt.ylabel('Residuals')
 
                 # save the current figure
-                fig1.savefig(f'./{newpath}/Residuals_vs_Par_{i+1}.pdf',
+                fig1.savefig(f'./{self.out_dir}/Residuals_vs_Par_{i+1}.pdf',
                              bbox_inches='tight')
                 # Destroy the current plot
                 plt.close()
@@ -895,7 +892,7 @@ class PostProcessing:
             plt.ylabel('Residuals')
 
             # save the current figure
-            fig2.savefig(f'./{newpath}/Fitted_vs_Residuals.pdf',
+            fig2.savefig(f'./{self.out_dir}/Fitted_vs_Residuals.pdf',
                          bbox_inches='tight')
             # Destroy the current plot
             plt.close()
@@ -921,7 +918,7 @@ class PostProcessing:
             ax.add_artist(at)
 
             # save the current figure
-            fig3.savefig(f'./{newpath}/Hist_NormResiduals.pdf',
+            fig3.savefig(f'./{self.out_dir}/Hist_NormResiduals.pdf',
                          bbox_inches='tight')
             # Destroy the current plot
             plt.close()
@@ -937,103 +934,61 @@ class PostProcessing:
             plt.grid(True)
 
             # save the current figure
-            plt.savefig(f'./{newpath}/QQPlot_NormResiduals.pdf',
+            plt.savefig(f'./{self.out_dir}/QQPlot_NormResiduals.pdf',
                         bbox_inches='tight')
             # Destroy the current plot
             plt.close()
 
     # -------------------------------------------------------------------------
-    def eval_pce_model_3d(self):
-        # 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!')
-
-        self.n_samples = 1000
-        n_samples = self.n_samples
-
-        # 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)$')
+    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.
 
-        plt.grid()
+        Raises
+        ------
+        AttributeError
+            This function is only applicable if the MetaModel input dimension is 2.
 
-        #  Saving the figure
-        newpath = f'Outputs_PostProcessing_{self.name}/'
-        if not os.path.exists(newpath):
-            os.makedirs(newpath)
+        Returns
+        -------
+        None.
 
-        # 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()
+        """
+        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)
         
-        # 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):
@@ -1106,11 +1061,6 @@ class PostProcessing:
         y_pce_val = self.pce_out_mean
         y_val = self.model_out_dict
 
-        # Open a pdf for the plots
-        newpath = f'Outputs_PostProcessing_{self.name}/'
-        if not os.path.exists(newpath):
-            os.makedirs(newpath)
-
         fig = plt.figure()
         # Fit the data(train the model)
         for key in y_pce_val.keys():
@@ -1151,7 +1101,7 @@ class PostProcessing:
             
             # save the current figure
             plot_name = key.replace(' ', '_')
-            fig.savefig(f'./{newpath}/Model_vs_PCEModel_{plot_name}.pdf',
+            fig.savefig(f'./{self.out_dir}/Model_vs_PCEModel_{plot_name}.pdf',
                         bbox_inches='tight')
 
             # Destroy the current plot
@@ -1169,6 +1119,10 @@ class PostProcessing:
             List of x values. The default is [].
         x_axis : str, optional
             Label of the x axis. The default is "x [m]".
+            
+        Raises
+        ------
+        AttributeError: This evaluation only support PCE-type models!
 
         Returns
         -------
@@ -1182,10 +1136,6 @@ class PostProcessing:
         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']))
         marker = cycle(('x', 'd', '+', 'o', '*'))
@@ -1231,11 +1181,8 @@ class PostProcessing:
 
             # save the current figure
             plot_name = key.replace(' ', '_')
-            fig.savefig(f'./{newpath}/Model_vs_PCEModel_{plot_name}.pdf',
+            fig.savefig(f'./{self.out_dir}/Model_vs_PCEModel_{plot_name}.pdf',
                         bbox_inches='tight')
 
             # 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