From 3129af4327019a67a320aa96a4ee285036107bd4 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Tue, 23 Feb 2021 11:21:47 +0100
Subject: [PATCH] [projections] move independent_ensemble_fit.py. create median
 one fold fit & visualizer with it.

---
 ...es_visualizer_for_non_stationary_models.py | 15 +++---
 .../one_fold_analysis/one_fold_fit.py         | 27 ----------
 .../abstract_ensemble_fit.py                  |  0
 .../independent_ensemble_fit/__init__.py      |  0
 .../independent_ensemble_fit.py               | 17 ++++---
 .../one_fold_fit_median.py                    | 20 ++++++++
 .../visualizer_median.py                      | 41 +++++++++++++++
 ...ation_temporal_for_projections_ensemble.py | 51 ++++++++++---------
 .../visualizer_for_projection_ensemble.py     | 46 +++++------------
 9 files changed, 119 insertions(+), 98 deletions(-)
 rename projects/projected_snowfall/elevation_temporal_model_for_projections/{ensemble_fit => }/abstract_ensemble_fit.py (100%)
 create mode 100644 projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/__init__.py
 rename projects/projected_snowfall/elevation_temporal_model_for_projections/{ensemble_fit => independent_ensemble_fit}/independent_ensemble_fit.py (68%)
 create mode 100644 projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_median.py
 create mode 100644 projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_median.py

diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
index 5ab36a43..9c3b7f9f 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py
@@ -58,22 +58,25 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self.altitude_group = get_altitude_group_from_altitudes(self.studies.altitudes)
         self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
         self.remove_physically_implausible_models = remove_physically_implausible_models
-        # Load one fold fit
+
         self.massif_name_to_massif_altitudes = {}
+        # Load one fold fit
+        self.load_one_fold_fit()
 
+        # Cache
+        self._method_name_and_order_to_massif_name_to_value = {}
+        self._method_name_and_order_to_max_abs = {}
+        self._max_abs_for_shape = None
+
+    def load_one_fold_fit(self):
         one_fold_fit_list = [self.fit_one_fold(massif_name) for massif_name in self.massif_names]
         self._massif_name_to_one_fold_fit = {m: o for m, o in zip(self.massif_names, one_fold_fit_list) if
                                              o is not None}
-
         # Print number of massif without any validated fit
         massifs_without_any_validated_fit = [massif_name
                                              for massif_name, old_fold_fit in self._massif_name_to_one_fold_fit.items()
                                              if not old_fold_fit.has_at_least_one_valid_model]
         print('Not validated:', len(massifs_without_any_validated_fit), massifs_without_any_validated_fit)
-        # Cache
-        self._method_name_and_order_to_massif_name_to_value = {}
-        self._method_name_and_order_to_max_abs = {}
-        self._max_abs_for_shape = None
 
     def fit_one_fold(self, massif_name):
         # Load valid massif altitudes
diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
index 29fd516c..65f900e2 100644
--- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
+++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py
@@ -392,33 +392,6 @@ class OneFoldFit(object):
     def best_residuals(self):
         return self.best_estimator.sorted_empirical_standard_gumbel_quantiles(split=Split.all)
 
-    # @property
-    # def bootstrap_data(self):
-    #     start = time.time()
-    #     bootstrap = []
-    #     for _ in range(AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP):
-    #         residuals = self.best_estimator.sorted_empirical_standard_gumbel_quantiles(split=Split.all)
-    #
-    #         # yield coordinate_values_to_maxima
-    #         bootstrap.append(coordinate_values_to_maxima)
-    #     end1 = time.time()
-    #     duration = str(datetime.timedelta(seconds=end1 - start))
-    #     print('bootstrap loader duration', duration)
-    #     return bootstrap
-
-    # def bootstrap_batch_data(self, batchsize=20):
-    #     bootstrap_batch_data = []
-    #     batch = []
-    #     len_batch = 0
-    #     for bootstrap in self.bootstrap_data:
-    #         batch.append(bootstrap)
-    #         len_batch += 1
-    #         if len_batch == batchsize:
-    #             yield batch
-    #             batch = []
-    #             len_batch = 0
-    #     return bootstrap_batch_data
-
     @cached_property
     def cached_results_from_bootstrap(self):
         start = time.time()
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/abstract_ensemble_fit.py
similarity index 100%
rename from projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/abstract_ensemble_fit.py
rename to projects/projected_snowfall/elevation_temporal_model_for_projections/abstract_ensemble_fit.py
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/__init__.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/independent_ensemble_fit.py
similarity index 68%
rename from projects/projected_snowfall/elevation_temporal_model_for_projections/ensemble_fit/independent_ensemble_fit.py
rename to projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/independent_ensemble_fit.py
index d7285684..f88b7f72 100644
--- a/projects/projected_snowfall/elevation_temporal_model_for_projections/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,10 @@
-from typing import Dict, Tuple, List
-
-from extreme_fit.model.margin_model.utils import MarginFitMethod
-from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import DefaultAltitudeGroup
 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.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
-from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.abstract_ensemble_fit import \
+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
 
 
 class IndependentEnsembleFit(AbstractEnsembleFit):
@@ -32,4 +28,11 @@ 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
+        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)
 
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_median.py
new file mode 100644
index 00000000..ba1dcf5f
--- /dev/null
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/one_fold_fit_median.py
@@ -0,0 +1,20 @@
+from typing import List
+
+import numpy as np
+
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
+
+
+class OneFoldFitMedian(OneFoldFit):
+
+    def __init__(self, one_fold_fit_list: List[OneFoldFit], massif_name, altitude_class, temporal_covariate_for_fit):
+        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
+
+    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])
+
+
+
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_median.py
new file mode 100644
index 00000000..30a346fc
--- /dev/null
+++ b/projects/projected_snowfall/elevation_temporal_model_for_projections/independent_ensemble_fit/visualizer_median.py
@@ -0,0 +1,41 @@
+import copy
+from typing import Dict, List
+
+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
+
+
+class VisualizerMedian(AltitudesStudiesVisualizerForNonStationaryModels):
+
+    def __init__(self, visualizers: List[AltitudesStudiesVisualizerForNonStationaryModels],
+                 model_classes,
+                 show=False,
+                 massif_names=None,
+                 fit_method=MarginFitMethod.extremes_fevd_mle,
+                 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):
+        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,
+                         display_only_model_that_pass_anderson_test=display_only_model_that_pass_anderson_test,
+                         confidence_interval_based_on_delta_method=confidence_interval_based_on_delta_method,
+                         remove_physically_implausible_models=remove_physically_implausible_models)
+
+    def load_one_fold_fit(self):
+        self._massif_name_to_one_fold_fit = {}
+        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
+
+    @property
+    def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
+        return self._massif_name_to_one_fold_fit
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 74d80cc3..546011ae 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
@@ -17,16 +17,13 @@ from extreme_fit.model.utils import set_seed_for_test
 from extreme_data.meteo_france_data.adamont_data.adamont.adamont_snowfall import AdamontSnowfall
 from extreme_data.meteo_france_data.adamont_data.adamont_scenario import AdamontScenario, get_gcm_rcm_couples, \
     rcp_scenarios
-from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.independent_ensemble_fit import \
+from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.independent_ensemble_fit import \
     IndependentEnsembleFit
 from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
     AnomalyTemperatureTemporalCovariate, TimeTemporalCovariate
 
 matplotlib.use('Agg')
 
-from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
-    plot_shoe_plot_changes_against_altitude, plot_histogram_all_trends_against_altitudes
-
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
     AbstractExtractEurocodeReturnLevel
 
@@ -36,33 +33,37 @@ from extreme_data.meteo_france_data.scm_models_data.utils import Season
 
 
 def main():
+    start = time.time()
     study_classes = [AdamontSnowfall][:1]
-    scenario = AdamontScenario.rcp45
-    gcm_rcm_couples = get_gcm_rcm_couples(scenario)
     ensemble_fit_class = [IndependentEnsembleFit]
     temporal_covariate_for_fit = [TimeTemporalCovariate, AnomalyTemperatureTemporalCovariate][1]
     set_seed_for_test()
     AbstractExtractEurocodeReturnLevel.ALPHA_CONFIDENCE_INTERVAL_UNCERTAINTY = 0.2
 
-    fast = True
-    if fast is None:
-        massif_names = None
-        gcm_rcm_couples = gcm_rcm_couples[:2]
-        AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-        altitudes_list = altitudes_for_groups[:2]
-    elif fast:
-        AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
-        massif_names = None
-        gcm_rcm_couples = [('EC-EARTH', 'RACMO22E')]
-        altitudes_list = altitudes_for_groups[1:2]
-    else:
-        massif_names = None
-        altitudes_list = altitudes_for_groups[:]
-
-    assert isinstance(gcm_rcm_couples, list)
-    start = time.time()
-    main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_class, scenario,
-              temporal_covariate_for_fit)
+    fast = False
+    scenarios = rcp_scenarios[:1] if fast is None else [AdamontScenario.rcp45]
+
+    for scenario in scenarios:
+        gcm_rcm_couples = get_gcm_rcm_couples(scenario)
+        if fast is None:
+            massif_names = ['Vanoise', 'Vercors']
+            gcm_rcm_couples = gcm_rcm_couples[:]
+            AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
+            altitudes_list = altitudes_for_groups[1:2]
+        elif fast:
+            AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP = 10
+            massif_names = None
+            gcm_rcm_couples = [('EC-EARTH', 'RACMO22E')]
+            altitudes_list = altitudes_for_groups[1:2]
+        else:
+            massif_names = None
+            altitudes_list = altitudes_for_groups[:]
+
+        assert isinstance(gcm_rcm_couples, list)
+
+        main_loop(gcm_rcm_couples, altitudes_list, massif_names, study_classes, ensemble_fit_class, scenario,
+                  temporal_covariate_for_fit)
+
     end = time.time()
     duration = str(datetime.timedelta(seconds=end - start))
     print('Total duration', duration)
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 79c75e4d..b13c2af4 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
@@ -1,39 +1,16 @@
-from collections import Counter
-from math import ceil, floor
-from typing import List, Dict
+from typing import List
 
-import matplotlib
-import matplotlib.pyplot as plt
-
-import numpy as np
-from cached_property import cached_property
-
-from extreme_data.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
-from extreme_data.meteo_france_data.scm_models_data.visualization.create_shifted_cmap import get_shifted_map, \
-    get_colors, ticks_values_and_labels_for_percentages, get_half_colormap, ticks_values_and_labels_for_positive_value, \
-    get_inverse_colormap, get_cmap_with_inverted_blue_and_green_channels, remove_the_extreme_colors
-from extreme_data.meteo_france_data.scm_models_data.visualization.main_study_visualizer import \
-    SCM_STUDY_CLASS_TO_ABBREVIATION, ALL_ALTITUDES_WITHOUT_NAN
-from extreme_data.meteo_france_data.scm_models_data.visualization.plot_utils import plot_against_altitude
-from extreme_data.meteo_france_data.scm_models_data.visualization.study_visualizer import StudyVisualizer
-from extreme_fit.distribution.gev.gev_params import GevParams
-from extreme_fit.function.margin_function.abstract_margin_function import AbstractMarginFunction
-from extreme_fit.function.param_function.linear_coef import LinearCoef
 from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
     AbstractSpatioTemporalPolynomialModel
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import \
-    get_altitude_group_from_altitudes, HighAltitudeGroup, VeyHighAltitudeGroup, MidAltitudeGroup
-from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
-    OneFoldFit
+    get_altitude_group_from_altitudes
 from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \
     plot_histogram_all_trends_against_altitudes, plot_shoe_plot_changes_against_altitude
 from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import compute_and_assign_max_abs
-from projects.projected_snowfall.elevation_temporal_model_for_projections.ensemble_fit.independent_ensemble_fit import \
+from projects.projected_snowfall.elevation_temporal_model_for_projections.independent_ensemble_fit.independent_ensemble_fit import \
     IndependentEnsembleFit
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
-from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 
 
 class MetaVisualizerForProjectionEnsemble(object):
@@ -89,17 +66,20 @@ class MetaVisualizerForProjectionEnsemble(object):
 
     def plot_independent(self):
         with_significance = False
-        # Individual plots
-        for independent_ensemble_fit in self.ensemble_fits(IndependentEnsembleFit):
-            print(independent_ensemble_fit)
-            for c, v in independent_ensemble_fit.gcm_rcm_couple_to_visualizer.items():
-                print(c, v.altitude_group)
-                v.plot_moments()
         # Aggregated at gcm_rcm_level plots
-        for gcm_rcm_couple in self.gcm_rcm_couples:
+        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
                                for independent_ensemble_fit in self.ensemble_fits(IndependentEnsembleFit)
                                ]
+            if gcm_rcm_couple is None:
+                for v in visualizer_list:
+                    v.studies.study.gcm_rcm_couple = ("Median", "merge")
+            for v in visualizer_list:
+                v.plot_moments()
             plot_histogram_all_trends_against_altitudes(self.massif_names, visualizer_list, with_significance=with_significance)
             for relative in [True, False]:
                 plot_shoe_plot_changes_against_altitude(self.massif_names, visualizer_list, relative=relative, with_significance=with_significance)
-- 
GitLab