[snow projection] add rank analysis & relative bias analysis

projects/projected_snowfall/comparison_with_scm/comparison_historical_visualizer.py

@@ -49,6 +49,26 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
                 pass
         return np.array(values), gcm_rcm_couples
 
+    def massif_name_to_relative_bias(self, plot_maxima=True):
+        study_method = self.get_study_method(plot_maxima)
+        massif_name_to_relative_bias = {}
+        for massif_name in self.massif_names:
+            values, gcm_rcm_couples = self.get_values(study_method, massif_name)
+            if len(gcm_rcm_couples) == len(self.adamont_studies.gcm_rcm_couples):
+                mean_values = np.mean(values, axis=1)
+                bias_in_the_mean = (mean_values - mean_values[0])[1:]
+                relative_bias = 100 * np.mean(bias_in_the_mean) / mean_values[0]
+                massif_name_to_relative_bias[massif_name] = relative_bias
+        return massif_name_to_relative_bias
+
+    def compute_real_bias_list_in_the_mean(self, massif_name, relative, study_method):
+        values, gcm_rcm_couples = self.get_values(study_method, massif_name)
+        mean_values = np.mean(values, axis=1)
+        bias_in_the_mean = (mean_values - mean_values[0])[1:]
+        if relative:
+            bias_in_the_mean *= 100 / mean_values[0]
+        return bias_in_the_mean, gcm_rcm_couples
+
     def compute_bias_list_in_the_mean(self, massif_name, relative, study_method):
         values, gcm_rcm_couples = self.get_values(study_method, massif_name)
         mean_values = np.mean(values, axis=1)
@@ -71,17 +91,15 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
     def massif_name_to_mean_bias_in_the_mean(self, plot_maxima=True, relative=False):
         return {m: np.mean(l) for m, l in self.massif_name_to_bias_list_in_the_mean(plot_maxima, relative).items()}
 
-    @property
-    def massif_name_to_rank(self):
+    def massif_name_to_rank(self, plot_maxima=True):
         massif_name_to_rank = {}
-        for massif_name, bias_list in self.massif_name_to_bias_list_in_the_mean(plot_maxima=True).items():
+        for massif_name, bias_list in self.massif_name_to_bias_list_in_the_mean(plot_maxima).items():
             # Count the number of bias negative
             nb_of_negative = sum([b < 0 for b in bias_list])
             # Rank starts to 0
             massif_name_to_rank[massif_name] = float(nb_of_negative)
         return massif_name_to_rank
 
-
     # Map gcm_rcm_couple to bias list (ordered by the massif_names)
 
     @cached_property
@@ -170,12 +188,13 @@ class ComparisonHistoricalVisualizer(StudyVisualizer):
         plt.close()
 
     def plot_map_with_the_rank(self):
-        massif_name_to_value = self.massif_name_to_rank
+        massif_name_to_value = self.massif_name_to_rank(plot_maxima=True)
         max_abs_change = self.adamont_studies.nb_ensemble_members
         ylabel = 'Rank of the mean maxima\n' \
                  'for the period {}-{} and the scenario {}\n' \
                  'which is between 0 (lowest) and {} (largest)'.format(self.study.year_min, self.study.year_max,
-                                                               scenario_to_str(self.study.scenario), max_abs_change)
+                                                                       scenario_to_str(self.study.scenario),
+                                                                       max_abs_change)
         plot_name = op.join('rank', ylabel)
         self.plot_map(cmap=plt.cm.coolwarm, graduation=1,
                       label=ylabel,

projects/projected_snowfall/comparison_with_scm/comparison_plot.py

+from collections import Counter
+
 import matplotlib.pyplot as plt
 from typing import Dict
 
@@ -13,38 +15,66 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
 def individual_plot(visualizer: ComparisonHistoricalVisualizer):
     # visualizer.adamont_studies.plot_maxima_time_series_adamont(visualizer.massif_names, visualizer.scm_study)
     # visualizer.shoe_plot_bias_maxima_comparison()
-    for relative in [True, False]:
-        visualizer.plot_map_with_the_mean_bias_in_the_mean(relative)
-
-    visualizer.plot_map_with_the_rank()
+    # for relative in [True, False]:
+    #     visualizer.plot_map_with_the_mean_bias_in_the_mean(relative)
+    # visualizer.plot_map_with_the_rank()
     # for relative
     # for plot_maxima in [True, False][:1]:
     #     visualizer.plot_comparison(plot_maxima)
+    pass
 
 
 def collective_plot(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]):
     visualizer = list(altitude_to_visualizer.values())[0]
+    altitudes = list(altitude_to_visualizer.keys())
+    print('Summary for rank for altitudes:', altitudes)
+    for plot_maxima in [True, False]:
+        study_setting = '{} {} {} {}-{}'.format('maxima' if plot_maxima else 'daily', ' values for ', scenario_to_str(visualizer.study.scenario),
+              visualizer.study.year_min, visualizer.study.year_max)
+        print(study_setting)
+        summary_rank(altitude_to_visualizer, visualizer, plot_maxima)
+        visualizer.plot_name = 'Histogram of the rank analysis for {}'.format(study_setting)
+        visualizer.show_or_save_to_file()
+        plt.close()
+
+        summary_bias(altitude_to_visualizer, visualizer, plot_maxima)
+
+
+def summary_bias(altitude_to_visualizer, visualizer, plot_maxima):
+    relative_bias_list = []
+    for v in altitude_to_visualizer.values():
+        relative_bias = np.array(list(v.massif_name_to_relative_bias(plot_maxima).values()))
+        relative_bias_list.extend(relative_bias)
+    mean_relative_bias = np.round(np.mean(relative_bias_list), 1)
+    print('Average relative bias w.r.t to the mean:', mean_relative_bias, '\%\n')
+
+
+def summary_rank(altitude_to_visualizer, visualizer, plot_maxima):
     count_number_of_total_massifs = 0
     count_number_of_time_the_reanalysis_is_the_smallest = 0
     count_number_of_time_the_reanalysis_is_the_biggest = 0
-    altitudes = list(altitude_to_visualizer.keys())
     all_ranks = []
     for v in altitude_to_visualizer.values():
-        ranks = np.array(list(v.massif_name_to_rank.values()))
+        ranks = np.array(list(v.massif_name_to_rank(plot_maxima).values()))
         count_number_of_total_massifs += len(ranks)
         count_number_of_time_the_reanalysis_is_the_smallest += sum(ranks == 0.0)
         all_ranks.extend(ranks)
-    print(scenario_to_str(visualizer.study.scenario), visualizer.study.year_min, visualizer.study.year_max)
-    print('Summary for rank for altitudes:', altitudes)
-    print('Mean ranks:', np.mean(all_ranks))
-    print('percentages of time reanalysis is the biggest:', 
-          100 * count_number_of_time_the_reanalysis_is_the_biggest / count_number_of_total_massifs)
-    print('number of time reanalysis is the biggest:', count_number_of_time_the_reanalysis_is_the_biggest,
-          ' out of ', count_number_of_total_massifs, ' time series')
-    print('percentages of time reanalysis is the smallest:', 
-          100 * count_number_of_time_the_reanalysis_is_the_smallest / count_number_of_total_massifs)
-    print('number of time reanalysis is the smallest:', count_number_of_time_the_reanalysis_is_the_smallest,
-          ' out of ', count_number_of_total_massifs, ' time series')
+    print('Average ranks w.r.t to the mean:', np.round(np.mean(all_ranks), 1))
+    print('percentages of time reanalysis is the biggest:',
+          np.round(100 * count_number_of_time_the_reanalysis_is_the_biggest / count_number_of_total_massifs, 1), '\%')
+    print('percentages of time reanalysis is the smallest:',
+          np.round(100 * count_number_of_time_the_reanalysis_is_the_smallest / count_number_of_total_massifs, 1), '\%')
+
+    # Plot histogram.
+    ax = plt.gca()
+    c = Counter(all_ranks)
+    x = list(range(15))
+    y = 100 * np.array([c[i] for i in x]) / len(all_ranks)
+    ax.bar(x, y, width=0.5, edgecolor='black')
+    ax.set_ylabel('Percentage of time series (\%)')
+    ax.set_xlabel('Rank w.r.t to the mean')
+
+
 
 
 def bias_of_the_mean_with_the_altitude(altitude_to_visualizer: Dict[int, ComparisonHistoricalVisualizer]):

projects/projected_snowfall/comparison_with_scm/main_comparison_on_rcp_85_extended.py

@@ -18,8 +18,7 @@ from projects.projected_snowfall.comparison_with_scm.comparison_historical_visua
     ComparisonHistoricalVisualizer
 
 
-def main():
-    fast = 7
+def main(fast):
     year_max = 2019
     # Set the year_min and year_max for the comparison
     if fast is 1:
@@ -40,16 +39,6 @@ def main():
         year_min = 2006
         altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
     elif fast is 5:
-        year_max = 2019
-        massif_names = None
-        year_min = 2013
-        altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
-    elif fast is 6:
-        year_max = 2012
-        massif_names = None
-        year_min = 2006
-        altitudes = [600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600][:]
-    elif fast is 7:
         year_max = 2019
         massif_names = None
         year_min = 1982
@@ -97,4 +86,6 @@ def load_visualizer(altitude, massif_names, year_min, year_max) -> ComparisonHis
 
 
 if __name__ == '__main__':
-    main()
+    fast_list = [2, 4, 6][1:]
+    for fast in fast_list:
+        main(fast)