diff --git a/video_prediction_tools/analysis_config/analysis_test.json b/video_prediction_tools/analysis_config/analysis_test.json
deleted file mode 100644
index 79c3e4073c71f84e3465ca2ae5e5e3ed62a161ee..0000000000000000000000000000000000000000
--- a/video_prediction_tools/analysis_config/analysis_test.json
+++ /dev/null
@@ -1,11 +0,0 @@
-
-{
-
-"results_dir": [ "/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2015to2017M01to12-64x64-3930N0000E-T2_MSL_gph500/convLSTM/20201221T181605_gong1_sunny",
-	        "/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2015to2017M01to12-64x64-3930N0000E-T2_MSL_gph500/savp/20201221T192753_gong1_savp_test1"],
-
-"metric": ["mse"],
-
-"compare_by":["model"]
-
-}
diff --git a/video_prediction_tools/main_scripts/best_model_selector.py b/video_prediction_tools/main_scripts/best_model_selector.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/video_prediction_tools/main_scripts/main_meta_postprocess.py b/video_prediction_tools/main_scripts/main_meta_postprocess.py
index 18dc83a95aaad4d8c65751f1d7054e1eb26e4201..3fda889c26ee4e91994800b7ae1ea254abf3b179 100644
--- a/video_prediction_tools/main_scripts/main_meta_postprocess.py
+++ b/video_prediction_tools/main_scripts/main_meta_postprocess.py
@@ -7,10 +7,9 @@ from __future__ import division
 from __future__ import print_function
 
 __email__ = "b.gong@fz-juelich.de"
-__author__ = "Bing Gong"
-__date__ = "2021-12-04"
-__update_date = "2022-02-02"
-
+__author__ = "Bing Gong, Yan Ji"
+__date__ = "2020-12-04"
+__update_date__ = "2022-02-02"
 import os
 from matplotlib.pylab import plt
 import json
@@ -18,47 +17,77 @@ import numpy as np
 import shutil
 import glob
 from netCDF4 import Dataset
-from  model_modules.video_prediction.metrics import *
 import xarray as xr
 
 
+def skill_score(tar_score,ref_score,best_score):
+    ss = (tar_score-ref_score) / (best_score-ref_score)
+    return ss
+
+
 class MetaPostprocess(object):
-    def __init__(self, analysis_config=None, analysis_dir=None, stochastic_ind=0, forecast_type="deterministic"):
+
+    def __init__(self, root_dir: str = "/p/project/deepacf/deeprain/video_prediction_shared_folder/",
+                 analysis_config: str = None, metric: str = "mse", exp_id: str=None, enable_skill_scores:bool=False):
         """
         This class is used for calculating the evaluation metric, analyize the models' results and make comparsion
         args:
-            analysis_config    :str, the path pointing to the aanalysis_configuration json file
+            root_dir           :str, the root directory for the shared folder
+            analysis_config    :str, the path pointing to the analysis_configuration json file
             analysis_dir       :str, the path to save the analysis results
-            forecast_type      :str, "deterministic" or "stochastic"
+            metric             :str,  based on which evalution metric for comparison, "mse","ssim", "texture"  and "acc"
+            exp_id             :str,  the given exp_id which is used as the name of postfix of the folder to store the plot
+            enable_skill_scores:bool, the
         """
-
+        self.root_dir = root_dir
         self.analysis_config = analysis_config
-        self.analysis_dir = analysis_dir
-        self.stochastic_ind=stochastic_ind
+        self.analysis_dir = os.path.join(root_dir, "meta_postprocess", exp_id)
+        self.metric = metric
+        self.exp_id = exp_id
+        self.enable_skill_scores = enable_skill_scores
+        self.models_type = []
+        self.metric_values = []  # return the shape: [num_results, persi_values, model_values]
+        self.skill_scores = []  # contain the calculated skill scores [num_results, skill_scores_values]
+
 
-    def __call__():
+    def __call__(self):
+        self.sanity_check()
         self.create_analysis_dir()
-        self.copy_analysis_config_to_analysis_dir()
+        self.copy_analysis_config()
         self.load_analysis_config()
-        self.load_results_dir_parameters()
-        self.calculate_evaluation_metrics()
-        self.save_metrics_all_dir_to_json()
-        self.make_comparsion()
+        self.get_metrics_values()
+        self.calculate_skill_scores()
+        if self.enable_skill_scores:
+            self.plot_skill_scores()
+        else:
+            self.plot_abs_scores()
 
+    def sanity_check(self):
+
+        self.available_metrics = ["mse", "ssim", "texture", "acc"]
+        if self.metric not in self.available_metrics:
+            raise ("The 'metric' must be one of the following:", available_metrics)
+        if type(self.exp_id) is not str:
+            raise ("'exp_id' must be 'str' type ")
 
     def create_analysis_dir(self):
         """
         Function to create the analysis directory if it does not exist
         """
-        if not os.path.exists(self.analysis_dir):os.makedirs(self.analysis_dir)
+        if not os.path.exists(self.analysis_dir): os.makedirs(self.analysis_dir)
+        print("1. Create analysis dir successfully: The result will be stored to the folder:", self.analysis_dir)
 
     def copy_analysis_config(self):
         """
         Copy the analysis configuration json to the analysis directory
         """
-        shutil.copy(self.analysis_config, os.path.join(self.analysis_dir,"analysis_config.json"))
-        self.analysis_config = os.path.join(self.analysis_dir,"analysis_config.json")
-
+        try:
+            shutil.copy(self.analysis_config, os.path.join(self.analysis_dir, "meta_config.json"))
+            self.analysis_config = os.path.join(self.analysis_dir, "meta_config.json")
+            print("2. Copy analysis config successs ")
+        except Exception as e:
+            print("The meta_config.json is not found in the dictory: ", self.analysis_config)
+        return None
 
     def load_analysis_config(self):
         """
@@ -66,163 +95,194 @@ class MetaPostprocess(object):
         """
         with open(self.analysis_config) as f:
             self.f = json.load(f)
-        self.metrics = self.f["metric"]
-        self.results_dirs = self.f["results_dir"]
-        self.compare_by = self.f["compare_by"]
-  
-    @staticmethod
-    def load_prediction_and_real_from_one_dir(results_dir,var="T2",stochastic_ind=0):
+
+        print("*****The following results will be compared and ploted*****")
+        [print(i) for i in self.f["results"].values()]
+        print("*******************************************************")
+        print("3. Loading analysis config success")
+
+        return None
+
+    def get_labels(self):
+        labels = list(self.f["results"].keys())
+        return labels
+
+    def get_meta_info(self):
+        """
+        get the model types meta information of the results from the options_checkpoints json file from postprocess stage
         """
-        Load the reference and prediction from nc file based on the select var
-        args:
-            var           : string, the target variable be retrieved from nc file
-            stochastic    : int   , the stochastic index 
-        return:
-            real_all      : list of list, all the reference values from all th nc files in the result directory
-            persistent_all: list of list, the persistent values from all the nc files in the result directory
-            forecast_all  : list of list, the forecasting values from all the nc files in the result directory
-            time_forecast : list of list, the forecasting timestamps from all the nc files in the results directory
-        """
-        fl_names = glob.glob(os.path.join(results_dir,"*.nc"))
-        real_all = []
-        persistent_all = []
-        forecast_all = []
-        for fl_nc in fl_names:
-            real,persistent, forecast,time_forecast = MetaPostprocess.read_values_by_var_from_nc(fl_nc,var,stochastic_ind)
-            real_all.append(real)
-            persistent_all.append(persistent)
-            forecast_all.append(forecast) 
-        return real_all,persistent_all,forecast_all, time_forecast
 
+        for i in self.f["results"].values():
+            option_checkpoints = os.path.join(i, "options_checkpoints.json")
+            with open(option_checkpoints) as f:
+                m = json.load(f)
+                self.models_type.append(m["model"])
+        return None
 
-    @staticmethod
-    def read_values_by_var_from_nc(fl_nc,var="T2",stochastic_ind=0):
+    def get_metrics_values(self):
         """
-        Function that read the referneces, persistent and forecasting values from one .nc file 
-        args:
-            fc_nc : str, the nc file full path
-            var   : str, the target variable
-            stochastic_ind: str, the stochastic index
-        Return:
-            real          : one dimension list, the reference values
-            persistent    : one dimension list, the persistent values
-            forecast      : one dimension list, the forecast values
-            time_forecast : one dimension list, the timestamps
-        """
-        #if not var in ["T2","MSL","GPH500"]: raise ValueError ("var name is not correct, should be 'T2','MSL',or 'GPH500'")
-        with Dataset(fl_nc, mode = 'r') as fl:
-           #load var prediction, real and persistent values
-           real = fl["/analysis/reference/"].variables[var][:]
-           persistent = fl["/analysis/persistent/"].variables[var][:]
-           forecast = fl["/forecasts/"+var+"/stochastic"].variables[str(stochastic_ind)][:]
-           time_forecast = fl.variables["time_forecast"][:]
-        return real, persistent, forecast, time_forecast   
+        get  the evaluation metric values of all the results, return a list  [results,persi, model]
+        """
+        self.get_meta_info()
 
+        for i, result_dir in enumerate(self.f["results"].values()):
+            vals = MetaPostprocess.get_one_metric_values(result_dir, self.metric, self.models_type[i])
+            self.metric_values.append(vals)
+        print("4. Get metrics values success")
+        return self.metric_values
 
     @staticmethod
-    def calculate_metric_one_img(real,forecast,metric="mse"):
-         """
-         Function that calculate the evaluation metric for one image
-         args:
-              real       : one dimension list contains the reference values for one image/sample
-              forecast   : one dimension list contains the forecasting values for one image/sample
-              metric     : str, the evaluation metric
+    def get_one_metric_values(result_dir: str = None, metric: str = "mse", model: str = None):
+
+        """
+        obtain the metric values (persistence and DL model) in the "evaluation_metrics.nc" file
+        return:  list contains the evaluatioin metrics of one result. [persi,model]
         """
+        filename = 'evaluation_metrics.nc'
+        filepath = os.path.join(result_dir, filename)
+        try:
+            with xr.open_dataset(filepath) as dfiles:
+                persi = np.array(dfiles['2t_persistence_{}_bootstrapped'.format(metric)][:])
+                model = np.array(dfiles['2t_{}_{}_bootstrapped'.format(model, metric)][:])
+                print("The values for evaluation metric '{}' values are obtained from file {}".format(metric, filepath))
+                return [persi, model]
+        except Exception as e:
+            print("!! The evallution metrics retrive from the {} fails".format(filepath))
+            print(e)
 
-        if metric == "mse":
-         #compare real and persistent
-            eval_forecast = mse_imgs(real, forecast)
-        elif metric == "psnr":
-            eval_forecast = psnr_imgs(real,forecast)   
-        return  eval_forecast
-    
-    @staticmethod
-    def calculate_metric_one_dir(results_dir,var,metric="mse",forecast_type="persistent",stochastic_ind=0):
+    def calculate_skill_scores(self):
         """
-        Function that calculate the evaluation metric for one directory
-        args:
-             results_dir    : the results_dir contains the results with nc files
-             var            : the target variable want to get
-             stochastic_ind : the stochastic index would like to obtain values
-             forecast_type  : str, either "forecast" or "persistent"
-        Return: 
-              eval_forecast_all : list of list that contains the evaluation values [num_samples,forecast_timestampes]
-        """
-        real_all, persistent_all, forecast_all, self.time_forecast = MetaPostprocess.load_prediction_and_real_from_one_dir(results_dir=results_dir, var=var, stochastic_ind=stochastic_ind)
-        eval_forecast_all = []
-        #loop for real data
-        for idx in range(len(real_all)):
-            eval_forecast_per_sample_over_ts = []
-            #loop the forecast time
-            for time in range(len(real_all[idx])):
-                #loop for each sample and each timestamp
-                eval_forecast = MetaPostprocess.calculate_metric_one_img(real_all[idx][time],forecast_all[idx][time], metric=metric)
-                eval_forecast_per_sample_over_ts.append(eval_forecast)
-            eval_forecast_all.append(list(eval_forecast_per_sample_over_ts))
-        return eval_forecast_all
-    
+        calculate the skill scores
+        """
+        if self.metric_values is None:
+            raise ("metric_values should be a list but None is provided")
+
+        best_score = 0
+        if self.metric == "mse":
+            pass
+
+        elif self.metric in ["ssim", "acc", "texture"]:
+            best_score = 1
+        else:
+            raise ("The metric should be one of the following available metrics :", self.available_metrics)
+
+        if self.enable_skill_scores:
+            for i in range(len(self.metric_values)):
+                skill_val = skill_score(self.metric_values[i][1], self.metric_values[i][0], best_score)
+                self.skill_scores.append(skill_val)
+
+            return self.skill_scores
+        else:
+            return None
+
+    def get_lead_time_labels(self):
+        leadtimes = self.metric_values[0][0].shape[1]
+        leadtimelist = ["leadhour" + str(i + 1) for i in range(leadtimes)]
+        return leadtimelist
+
+    def config_plots(self):
+        self.leadtimelist = self.get_lead_time_labels()
+        self.labels = self.get_labels()
+        self.markers = self.f["markers"]
+        self.colors = self.f["colors"]
+        self.n_leadtime = len(self.leadtimelist)
+
     @staticmethod
-    def reshape_eval_to_one_dim(values):
-        return np.array(values).flatten()
-
-    def calculate_metric_all_dirs(self, forecast_type="persistent", metric="mse"):
-        """
-        Return the evaluation metrics for persistent and forecasing model over forecasting timestampls
-        
-        return:
-               eval_forecast : list, the evaluation metric values for persistent  with respect to the dimenisons [results_dir,samples,timestampe]
-               forecast_type : str, either "forecast" or "persistent"
-
-        """
-        if forecast_type not in ["persistent","forecast"]: raise ValueError("forecast type should be either 'persistent' or 'forecast'")
-        eval_forecast_all_dirs = []
-        for results_dir in self.results_dirs:
-            eval_forecast_all = MetaPostprocess.calculate_metric_one_dir(results_dir,var,metric="mse",forecast_type="persistent",stochastic_ind=0)
-            eval_forecast_all_dirs.append(eval_forecast_all)
-
-        times = list(range(len(self.time_forecast)))
-        samples = list(range(len(real_all)))
-        print("shape of list",np.array(eval_forecast_all_dirs).shape)
-        evals_forecast = xr.DataArray(eval_forecast_all_dirs, coords=[self.results_dirs, samples , times], dims=["results_dirs", "samples","time_forecast"])
-        return evals_forecast
-
-
-
-
-    def load_results_dir_parameters(self,compare_by="model"):
-        self.compare_by_values = []
-        for results_dir in self.results_dirs:
-            with open(os.path.join(results_dir, "options_checkpoints.json")) as f:
-                self.options = json.loads(f.read())
-                print("self.options:",self.options)
-                #if self.compare_by == "model":
-                self.compare_by_values.append(self.options[compare_by])
-
-    
-    def plot_results(self,one_persistent=True):
-        """
-        Plot the mean and vars for the user-defined metrics
-        """
-        self.load_results_dir_parameters(compare_by="model")
-        evals_forecast = self.calculate_metric_all_dirs(is_persistent=False,metric="mse")
-        t = evals_forecast["time_forecast"]
-        mean_forecast = evals_forecast.groupby("time_forecast").mean(dim="samples").values
-        var_forecast = evals_forecast.groupby("time_forecast").var(dim="samples").values
-        print("mean_foreast",mean_forecast)
-        x = np.array(self.compare_by_values)
-        y = np.array(mean_forecast)
-        e = np.array(var_forecast)
-       
-       # plt.errorbar(t,y[0],e[0],label="convlstm")
-        plt.errorbar(t,y[0],e[0],label="savp")
+    def map_ylabels(metric):
+        if metric == "mse":
+            ylabel = "MSE"
+        elif metric == "acc":
+            ylabel = "ACC"
+        elif metric == "ssim":
+            ylabel = "SSIM"
+        elif metric == "texture":
+            ylabel = "Ratio of gradient ($r_G$)"
+        else:
+            raise ("The metric is not correct!")
+        return ylabel
+
+    def plot_abs_scores(self):
+        self.config_plots()
+
+        fig = plt.figure(figsize = (8, 6))
+        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
+        for i in range(len(self.metric_values)):
+            score_plot = np.nanquantile(self.metric_values[i][1], 0.5, axis = 0)
+            plt.plot(np.arange(1, 1 + self.n_leadtime), score_plot, label = self.labels[i], color = self.colors[i],
+                     marker = self.markers[i], markeredgecolor = 'k', linewidth = 1.2)
+            plt.fill_between(np.arange(1, 1 + self.n_leadtime),
+                             np.nanquantile(self.metric_values[i][1], 0.95, axis = 0),
+                             np.nanquantile(self.metric_values[i][1], 0.05, axis = 0), color = self.colors[i],
+                             alpha = 0.2)
+
+            if self.models_type[i] == "convLSTM":
+                score_plot = np.nanquantile(self.metric_values[i][0], 0.5, axis = 0)
+                plt.plot(np.arange(1, 1 + self.n_leadtime), score_plot, label = "Persi_cv" + str(i),
+                         color = self.colors[i], marker = "D", markeredgecolor = 'k', linewidth = 1.2)
+                plt.fill_between(np.arange(1, 1 + self.n_leadtime),
+                                 np.nanquantile(self.metric_values[i][0], 0.95, axis = 0),
+                                 np.nanquantile(self.metric_values[i][0], 0.05, axis = 0), color = "b", alpha = 0.2)
+
+        plt.yticks(fontsize = 16)
+        plt.xticks(np.arange(1, 13), np.arange(1, 13, 1), fontsize = 16)
+        legend = ax.legend(loc = 'upper right', bbox_to_anchor = (1.46, 0.95),
+                           fontsize = 14)  # 'upper right', bbox_to_anchor=(1.38, 0.8),
+        ylabel = MetaPostprocess.map_ylabels(self.metric)
+        ax.set_xlabel("Lead time (hours)", fontsize = 21)
+        ax.set_ylabel(ylabel, fontsize = 21)
+        fig_path = os.path.join(self.analysis_dir, self.metric + "_abs_values.png")
+        # fig_path = os.path.join(prefix,fig_name)
+        plt.savefig(fig_path, bbox_inches = "tight")
         plt.show()
-        plt.savefig(os.path.join(self.analysis_dir,self.metrics[0]+".png"))
         plt.close()
-       
-          
+        print("The plot saved to {}".format(fig_path))
+
+    def plot_skill_scores(self):
+        """
+        Plot the skill scores once the enable_skill is True
+        """
+        self.config_plots()
+        fig = plt.figure(figsize = (8, 6))
+        ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
+        for i in range(len(self.skill_scores)):
+            if self.models_type[i] == "convLSTM":
+                c = "r"
+            elif self.models_type[i] == "savp":
+                c = "b"
+            else:
+                raise ("current only support convLSTM and SAVP for plotinig the skil scores")
+
+            plt.boxplot(self.skill_scores[i], positions = np.arange(1, self.n_leadtime + 1), medianprops = {'color': c},
+                        capprops = {'color': c}, boxprops = {'color': c}, showfliers = False)
+            score_plot = np.nanquantile(self.skill_scores[i], 0.5, axis = 0)
+            plt.plot(np.arange(1, 1 + self.n_leadtime), score_plot, color = c, linewidth = 1.2, label = self.labels[i])
+
+        legend = ax.legend(loc = 'upper right', bbox_to_anchor = (1.46, 0.95), fontsize = 14)
+        plt.yticks(fontsize = 16)
+        plt.xticks(np.arange(1, 13), np.arange(1, 13, 1), fontsize = 16)
+        ax.set_xlabel("Lead time (hours)", fontsize = 21)
+        ax.set_ylabel("Skill scores of {}".format(self.metric), fontsize = 21)
+        fig_path = os.path.join(self.analysis_dir, self.metric + "_skill_scores.png")
+        plt.savefig(fig_path, bbox_inches = "tight")
+        plt.show()
+        plt.close()
+        print("The plot saved to {}".format(fig_path))
+
 
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--analysis_config", type=str, required=True, help="The path points to the  meta_postprocess configuration file.",
+                        default="../meta_postprocess_config/meta_config.json")
+    parser.add_argument("--metric", help="Based on which the models are compared, the value should be in one of [mse,ssim,acc,texture]",default="mse")
+    parser.add_argument("--exp_id", help="The experiment id which will be used as postfix of the output directory",default="exp1")
+    parser.add_argument("--enable_skill_scores", help="compared by skill scores or the absolute evaluation values",default=True)
+    args = parser.parse_args()
 
+    meta = MetaPostprocess(analysis_config=args.analysis_config, metric=args.metric, exp_id=args.metric,
+                           enable_skill_scores=args.enable_skill_scores)
+    meta()
 
 
+if __name__ == '__main__':
+    main()
 
-   
diff --git a/video_prediction_tools/meta_postprocess_config/meta_config.json b/video_prediction_tools/meta_postprocess_config/meta_config.json
new file mode 100644
index 0000000000000000000000000000000000000000..96ed1f2d9b1725d39d4841971f4731186bac7c2a
--- /dev/null
+++ b/video_prediction_tools/meta_postprocess_config/meta_config.json
@@ -0,0 +1,12 @@
+{
+"results":
+{"ConvLSTM_cv1":"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM/20210712T104840_gong1_convLSTM_cv11_epoch20",
+"ConvLSTM_cv2":"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM/20210709T225929_gong1_convLSTM_epoch20",
+"ConvLSTM_cv3":"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/convLSTM/20210712T110410_gong1_convLSTM_cv13_epoch20",
+"SAVP_cv1":"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210709T224456_gong1_savp_cv11_best_model",
+"SAVP_cv2":"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210901T090059_gong1_savp_cv12/",
+"SAVP_cv3":"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210709T225131_gong1_savp_cv13_best_model"},
+"markers":["X","X","X","P","P","P"],
+"colors":["r","g","b","r","g","b"]
+
+}