From 1c1ec59881d2e8a2992f29cb2b3f5942060eb74b Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 23 Feb 2021 14:53:53 +0100
Subject: [PATCH] [projections] add mean merge function for the ensemble. fix
 max abs for median/mean visualizer. Fix OneFoldFitMerge that was given empty
 list as input.

---
 .../utils_altitude_studies_visualizer.py      |  7 +++--
 .../independent_ensemble_fit.py               | 28 +++++++++++++------
 ...ld_fit_median.py => one_fold_fit_merge.py} |  9 ++++--
 ...sualizer_median.py => visualizer_merge.py} | 19 ++++++++-----
 ...ation_temporal_for_projections_ensemble.py |  8 +++---
 .../visualizer_for_projection_ensemble.py     | 23 +++++++++------
 6 files changed, 60 insertions(+), 34 deletions(-)
 rename projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/{one_fold_fit_median.py => one_fold_fit_merge.py} (58%)
 rename projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/{visualizer_median.py => visualizer_merge.py} (74%)

diff --git a/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
index 058d033c..e0a556db 100644
--- a/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
+++ b/projects/altitude_spatial_model/altitudes_fit/utils_altitude_studies_visualizer.py
@@ -27,18 +27,19 @@ def load_visualizer_list(season, study_class, altitudes_list, massif_names, mode
 
 def compute_and_assign_max_abs(visualizer_list):
     # Compute the max abs for all metrics
-    d = {}
+    method_name_and_order_to_max_abs = {}
     for method_name in AltitudesStudiesVisualizerForNonStationaryModels.moment_names:
         for order in AltitudesStudiesVisualizerForNonStationaryModels.orders:
             c = (method_name, order)
             max_abs = max([
                 max([abs(e) for e in v.method_name_and_order_to_d(method_name, order).values()
                      ]) for v in visualizer_list])
-            d[c] = max_abs
+            method_name_and_order_to_max_abs[c] = max_abs
     # Assign the max abs dictionary
     for v in visualizer_list:
-        v._method_name_and_order_to_max_abs = d
+        v._method_name_and_order_to_max_abs = method_name_and_order_to_max_abs
     # Compute the max abs for the shape parameter
     max_abs_for_shape = max([max([abs(e) for e in v.massif_name_to_shape.values()]) for v in visualizer_list])
     for v in visualizer_list:
         v._max_abs_for_shape = max_abs_for_shape
+    return method_name_and_order_to_max_abs, max_abs_for_shape
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/independent_ensemble_fit.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/independent_ensemble_fit.py
index f88b7f72..f8b2fd67 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/independent_ensemble_fit.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/independent_ensemble_fit.py
@@ -1,14 +1,18 @@
+import numpy as np
+
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
 from projects.projected_snowfall.elevation_temporal_model_for_projections.abstract_ensemble_fit import \
     AbstractEnsembleFit
-from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.visualizer_median import \
-    VisualizerMedian
+from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.visualizer_merge import \
+    VisualizerMerge
 
 
 class IndependentEnsembleFit(AbstractEnsembleFit):
     """For each gcm_rcm_couple, we create a OneFoldFit"""
+    Median_merge = 'Median'
+    Mean_merge = 'Mean'
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -28,11 +32,19 @@ class IndependentEnsembleFit(AbstractEnsembleFit):
                                                                           self.confidence_interval_based_on_delta_method,
                                                                           self.remove_physically_implausible_models)
             self.gcm_rcm_couple_to_visualizer[gcm_rcm_couple] = visualizer
-        # Load merge visualizer
+        # Load merge visualizer for various merge functions
         visualizers = list(self.gcm_rcm_couple_to_visualizer.values())
-        self.median_visualizer = VisualizerMedian(visualizers, self.models_classes, False, self.massif_names,
-                                                  self.fit_method, self.temporal_covariate_for_fit,
-                                                  self.only_models_that_pass_goodness_of_fit_test,
-                                                  self.confidence_interval_based_on_delta_method,
-                                                  self.remove_physically_implausible_models)
+        merge_function_name_to_merge_function = {
+            self.Median_merge: np.median,
+            self.Mean_merge: np.mean
+        }
+        self.merge_function_name_to_visualizer = {
+        name: VisualizerMerge(visualizers, self.models_classes, False, self.massif_names,
+                                               self.fit_method, self.temporal_covariate_for_fit,
+                                               self.only_models_that_pass_goodness_of_fit_test,
+                                               self.confidence_interval_based_on_delta_method,
+                                               self.remove_physically_implausible_models,
+                                               merge_function=merge_function)
+            for name, merge_function in merge_function_name_to_merge_function.items()
+        }
 
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_median.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_merge.py
similarity index 58%
rename from projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_median.py
rename to projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_merge.py
index ba1dcf5f..1fbcc2ba 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_median.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_merge.py
@@ -5,16 +5,19 @@ import numpy as np
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
 
 
-class OneFoldFitMedian(OneFoldFit):
+class OneFoldFitMerge(OneFoldFit):
 
-    def __init__(self, one_fold_fit_list: List[OneFoldFit], massif_name, altitude_class, temporal_covariate_for_fit):
+    def __init__(self, one_fold_fit_list: List[OneFoldFit], massif_name, altitude_class, temporal_covariate_for_fit,
+                 merge_function=np.median):
+        assert len(one_fold_fit_list) > 0
         self.one_fold_fit_list = one_fold_fit_list
         self.altitude_group = altitude_class()
         self.massif_name = massif_name
         self.temporal_covariate_for_fit = temporal_covariate_for_fit
+        self.merge_function = merge_function
 
     def get_moment(self, altitude, temporal_covariate, order=1):
-        return np.median([o.get_moment(altitude, temporal_covariate, order) for o in self.one_fold_fit_list])
+        return self.merge_function([o.get_moment(altitude, temporal_covariate, order) for o in self.one_fold_fit_list])
 
 
 
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_median.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_merge.py
similarity index 74%
rename from projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_median.py
rename to projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_merge.py
index 30a346fc..907cccf2 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_median.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_merge.py
@@ -1,15 +1,17 @@
 import copy
 from typing import Dict, List
 
+import numpy as np
+
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
-from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.one_fold_fit_median import \
-    OneFoldFitMedian
+from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.one_fold_fit_merge import \
+    OneFoldFitMerge
 
 
-class VisualizerMedian(AltitudesStudiesVisualizerForNonStationaryModels):
+class VisualizerMerge(AltitudesStudiesVisualizerForNonStationaryModels):
 
     def __init__(self, visualizers: List[AltitudesStudiesVisualizerForNonStationaryModels],
                  model_classes,
@@ -19,7 +21,9 @@ class VisualizerMedian(AltitudesStudiesVisualizerForNonStationaryModels):
                  temporal_covariate_for_fit=None,
                  display_only_model_that_pass_anderson_test=True,
                  confidence_interval_based_on_delta_method=False,
-                 remove_physically_implausible_models=False):
+                 remove_physically_implausible_models=False,
+                 merge_function=np.median):
+        self.merge_function = merge_function
         self.visualizers = visualizers
         super().__init__(studies=visualizers[0].studies, model_classes=model_classes, show=show, massif_names=massif_names,
                          fit_method=fit_method, temporal_covariate_for_fit=temporal_covariate_for_fit,
@@ -32,9 +36,10 @@ class VisualizerMedian(AltitudesStudiesVisualizerForNonStationaryModels):
         for massif_name in self.massif_names:
             one_fold_fit_list = [v.massif_name_to_one_fold_fit[massif_name] for v in self.visualizers
                                  if massif_name in v.massif_name_to_one_fold_fit]
-            one_fold_fit_merge = OneFoldFitMedian(one_fold_fit_list, massif_name,
-                                                  type(self.altitude_group), self.temporal_covariate_for_fit)
-            self._massif_name_to_one_fold_fit[massif_name] = one_fold_fit_merge
+            if len(one_fold_fit_list) > 0:
+                one_fold_fit_merge = OneFoldFitMerge(one_fold_fit_list, massif_name,
+                                                     type(self.altitude_group), self.temporal_covariate_for_fit)
+                self._massif_name_to_one_fold_fit[massif_name] = one_fold_fit_merge
 
     @property
     def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
index 546011ae..268b38fc 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/main_elevation_temporal_for_projections_ensemble.py
@@ -40,16 +40,16 @@ def main():
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = False
-    scenarios = rcp_scenarios[:1] if fast is None else [AdamontScenario.rcp45]
+    fast = None
+    scenarios = rcp_scenarios[:1] if fast is False else [AdamontScenario.rcp26]
 
     for scenario in scenarios:
         gcm_rcm_couples = get_gcm_rcm_couples(scenario)
         if fast is None:
-            massif_names = ['Vanoise', 'Vercors']
+            massif_names = None
             gcm_rcm_couples = gcm_rcm_couples[:]
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-            altitudes_list = altitudes_for_groups[1:2]
+            altitudes_list = altitudes_for_groups[:1]
         elif fast:
             AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
             massif_names = None
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
index b13c2af4..9d51602b 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/visualizer_for_projection_ensemble.py
@@ -60,24 +60,29 @@ class MetaVisualizerForProjectionEnsemble(object):
             for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
                 visualizer_list.extend(list(ensemble_fit.gcm_rcm_couple_to_visualizer.values()))
             # Potentially I could add more visualizer here...
-            compute_and_assign_max_abs(visualizer_list)
+            method_name_and_order_to_max_abs, max_abs_for_shape = compute_and_assign_max_abs(visualizer_list)
+            # Assign the same max abs for the 
+            for ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
+                for v in ensemble_fit.merge_function_name_to_visualizer.values():
+                    v._max_abs_for_shape = max_abs_for_shape
+                    v._method_name_and_order_to_max_abs = method_name_and_order_to_max_abs
             # Plot
             self.plot_independent()
 
     def plot_independent(self):
         with_significance = False
         # Aggregated at gcm_rcm_level plots
-        gcm_rcm_couples = self.gcm_rcm_couples + [None]
-        if None in gcm_rcm_couples:
-            assert gcm_rcm_couples[-1] is None
-        for gcm_rcm_couple in gcm_rcm_couples:
-            visualizer_list = [independent_ensemble_fit.gcm_rcm_couple_to_visualizer[gcm_rcm_couple]
-                               if gcm_rcm_couple is not None else independent_ensemble_fit.median_visualizer
+        merge_keys = [IndependentEnsembleFit.Median_merge, IndependentEnsembleFit.Mean_merge]
+        keys = self.gcm_rcm_couples + merge_keys
+        for key in keys:
+            visualizer_list = [independent_ensemble_fit.gcm_rcm_couple_to_visualizer[key]
+                               if key in self.gcm_rcm_couples
+                               else independent_ensemble_fit.merge_function_name_to_visualizer[key]
                                for independent_ensemble_fit in self.ensemble_fits(IndependentEnsembleFit)
                                ]
-            if gcm_rcm_couple is None:
+            if key in merge_keys:
                 for v in visualizer_list:
-                    v.studies.study.gcm_rcm_couple = ("Median", "merge")
+                    v.studies.study.gcm_rcm_couple = (key, "merge")
             for v in visualizer_list:
                 v.plot_moments()
             plot_histogram_all_trends_against_altitudes(self.massif_names, visualizer_list, with_significance=with_significance)
-- 
GitLab