[PAPER 1] refactor folder names in the paper snow load folder. add mani gelman...

[PAPER 1] refactor folder names in the paper snow load folder. add mani gelman convergence test. add the possibility to change the number of iterations for the bayesian fevd estimation

experiment/meteo_france_data/scm_models_data/visualization/study_visualization/main_study_visualizer.py

@@ -3,7 +3,7 @@ from typing import List
 
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     StudyVisualizer
-from experiment.paper_past_snow_loads.result_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+from experiment.paper_past_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
     CrocusDifferenceSnowLoad, \
     CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
     CrocusSnowDepthDifference, CrocusSnowDepthAtMaxofSwe

experiment/paper_past_snow_loads/check_mcmc_convergence_for_return_levels/gelman_convergence_test.py

+
+def compute_gelman_convergence_value(non_null_years_and_maxima, mcmc_iterations, model_class, nb_chains):
+    return 1.0
+
+# rbeta(10000, 6, 9) - 0.5
\ No newline at end of file

experiment/paper_past_snow_loads/result_mcmc_check_for_return_levels/main_bayesian_mcmc.py → experiment/paper_past_snow_loads/check_mcmc_convergence_for_return_levels/main_bayesian_mcmc.py

 import seaborn as sns
 
 import matplotlib.pyplot as plt
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSwe3Days, CrocusSnowLoadTotal
 from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
 from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
     fitted_linear_margin_estimator
@@ -31,7 +31,8 @@ def main_drawing_bayesian():
     colors = ['r', 'g', 'b', 'tab:orange']
     gev_param_name_to_color = dict(zip(GevParams.PARAM_NAMES, colors))
     gev_param_name_to_ax = dict(
-        zip(GevParams.PARAM_NAMES, [ax_trajectories, ax_trajectories, ax_trajectories_inverted]))
+        zip(GevParams.PARAM_NAMES, [ax_trajectories, ax_trajectories, ax_trajectories]))
+        # zip(GevParams.PARAM_NAMES, [ax_trajectories, ax_trajectories, ax_trajectories_inverted]))
     gev_param_name_to_label = {n: GevParams.greek_letter_from_gev_param_name(n) for n in GevParams.PARAM_NAMES}
     for j, gev_param_name in enumerate(GevParams.PARAM_NAMES[:]):
         label = gev_param_name_to_label[gev_param_name]
@@ -77,13 +78,14 @@ def main_drawing_bayesian():
 
 
 def get_return_level_bayesian_example():
-    maxima, years = CrocusSwe3Days(altitude=1800).annual_maxima_and_years('Vercors')
+    maxima, years = CrocusSnowLoadTotal(altitude=1800).annual_maxima_and_years('Vercors')
     # plt.plot(years, maxima)
     # plt.show()
     model_class = StationaryTemporalModel
     coordinates, dataset = load_temporal_coordinates_and_dataset(maxima, years)
     fitted_estimator = fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year=1958,
-                                                      fit_method=TemporalMarginFitMethod.extremes_fevd_bayesian)
+                                                      fit_method=TemporalMarginFitMethod.extremes_fevd_bayesian,
+                                                      nb_iterations_for_bayesian_fit=100000)
     return_level_bayesian = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator=fitted_estimator,
                                                                              ci_method=ConfidenceIntervalMethodFromExtremes.my_bayes,
                                                                              temporal_covariate=2017)

experiment/paper_past_snow_loads/check_mcmc_convergence_for_return_levels/main_gelman_convergence_test.py

+import pandas as pd
+from experiment.paper_past_snow_loads.paper_utils import paper_altitudes, paper_study_classes, \
+    load_altitude_to_visualizer
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+
+
+def gelman_convergence_test(mcmc_iterations, model_class, altitudes, study_class, nb_chains=3, massif_names=None):
+    """ Test if a given number of MCMC iterations is enough for the convergence of each parameter of the model_class
+     for every time series with non-null values present in the study_classes
+
+     Ideally, return a DataFrame with altitude as columns, and massif name as index that contain the R score
+     Then we could compute the max R that should ideally be below 1.2
+     """
+
+    altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names=massif_names,
+                                                         non_stationary_uncertainty=None,
+                                                         study_class=study_class, uncertainty_methods=None)
+    altitude_to_d = {}
+    for altitude, vizu in altitude_to_visualizer.items():
+        altitude_to_d[altitude] = vizu.massif_name_to_gelman_convergence_value(mcmc_iterations, model_class, nb_chains)
+    massif_names = list(altitude_to_d[altitudes[0]].keys())
+    df = pd.DataFrame(altitude_to_d, index=massif_names, columns=altitudes)
+    return df
+
+
+
+"""
+test gelman for the 4 types of models
+and the for the 3 variables considered: GSL, GSL from eurocode, GLS in 3 days
+"""
+
+if __name__ == '__main__':
+    mcmc_iterations = 1000
+    df = gelman_convergence_test(mcmc_iterations=mcmc_iterations, altitudes=paper_altitudes[:1],
+                                 study_class=paper_study_classes[0], model_class=StationaryTemporalModel,
+                                 massif_names=['Chartreuse'])
+    print(mcmc_iterations)
+    print(df.head())
+    print('Overall maxima:', df.max().max())

experiment/paper_past_snow_loads/result_data_comparison_with_eurocode/main_comparison_with_eurocode.py → experiment/paper_past_snow_loads/discussion_data_comparison_with_eurocode/main_comparison_with_eurocode.py

@@ -6,7 +6,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
     StudyVisualizer
 import matplotlib.pyplot as plt
 
-from experiment.paper_past_snow_loads.result_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
+from experiment.paper_past_snow_loads.discussion_data_comparison_with_eurocode.crocus_study_comparison_with_eurocode import \
     CrocusDifferenceSnowLoad, \
     CrocusSnowDensityAtMaxofSwe, CrocusDifferenceSnowLoadRescaledAndEurocodeToSeeSynchronization, \
     CrocusSnowDepthAtMaxofSwe, CrocusSnowDepthDifference

experiment/paper_past_snow_loads/paper_utils.py

+from collections import OrderedDict
+
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode, \
+    CrocusSnowLoad3Days
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
+    StudyVisualizerForNonStationaryTrends
+
+paper_altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700]
+paper_study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days]
+
+
+def load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty, study_class, uncertainty_methods):
+    altitude_to_visualizer = OrderedDict()
+    for altitude in altitudes:
+        altitude_to_visualizer[altitude] = StudyVisualizerForNonStationaryTrends(
+            study=study_class(altitude=altitude), multiprocessing=True, save_to_file=True,
+            uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
+            non_stationary_contexts=non_stationary_uncertainty)
+    return altitude_to_visualizer

experiment/paper_past_snow_loads/result_trends_and_return_levels/eurocode_visualizer.py

-from typing import Dict, List, Tuple
+from typing import Dict
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 from experiment.eurocode_data.utils import EUROCODE_RETURN_LEVEL_STR, EUROCODE_ALTITUDES
-from experiment.paper_past_snow_loads.result_trends_and_return_levels.study_visualizer_for_non_stationary_trends import \
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
     StudyVisualizerForNonStationaryTrends
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
     AbstractExtractEurocodeReturnLevel
-from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
-    EurocodeConfidenceIntervalFromExtremes
 from experiment.eurocode_data.massif_name_to_departement import massif_name_to_eurocode_region
 from experiment.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes
 from root_utils import get_display_name_from_object_type

experiment/paper_past_snow_loads/result_trends_and_return_levels/main_result_trends_and_return_levels.py

-from collections import OrderedDict
-import os.path as op
-
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 
-from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_map
-from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
-    StudyVisualizer
+from experiment.paper_past_snow_loads.paper_utils import paper_study_classes, paper_altitudes
 from experiment.paper_past_snow_loads.result_trends_and_return_levels.eurocode_visualizer import \
     plot_uncertainty_massifs
-from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal, CrocusSnowLoadEurocode, \
-    CrocusSnowLoad3Days
-from experiment.paper_past_snow_loads.result_trends_and_return_levels.study_visualizer_for_non_stationary_trends import \
+from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.paper_past_snow_loads.study_visualizer_for_non_stationary_trends import \
     StudyVisualizerForNonStationaryTrends
 from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
     ConfidenceIntervalMethodFromExtremes
-from root_utils import VERSION_TIME
 
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
@@ -42,12 +35,8 @@ def intermediate_result(altitudes, massif_names=None,
     :return:
     """
     # Load altitude to visualizer
-    altitude_to_visualizer = OrderedDict()
-    for altitude in altitudes:
-        altitude_to_visualizer[altitude] = StudyVisualizerForNonStationaryTrends(
-            study=study_class(altitude=altitude), multiprocessing=True, save_to_file=True,
-            uncertainty_massif_names=massif_names, uncertainty_methods=uncertainty_methods,
-            non_stationary_contexts=non_stationary_uncertainty)
+    altitude_to_visualizer = load_altitude_to_visualizer(altitudes, massif_names, non_stationary_uncertainty,
+                                                         study_class, uncertainty_methods)
     # Plot trends
     max_abs_tdrl = max([visualizer.max_abs_tdrl for visualizer in altitude_to_visualizer.values()])
     for visualizer in altitude_to_visualizer.values():
@@ -57,15 +46,14 @@ def intermediate_result(altitudes, massif_names=None,
     return altitude_to_visualizer
 
 
+
 def major_result():
-    altitudes = [300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700]
     uncertainty_methods = [ConfidenceIntervalMethodFromExtremes.my_bayes,
                            ConfidenceIntervalMethodFromExtremes.ci_mle][:]
     massif_names = None
     non_stationary_uncertainty = [False, True][:]
-    study_classes = [CrocusSnowLoadTotal, CrocusSnowLoadEurocode, CrocusSnowLoad3Days]
-    for study_class in study_classes[2:]:
-        intermediate_result(altitudes, massif_names, non_stationary_uncertainty, uncertainty_methods, study_class)
+    for study_class in paper_study_classes[2:]:
+        intermediate_result(paper_altitudes, massif_names, non_stationary_uncertainty, uncertainty_methods, study_class)
 
 
 if __name__ == '__main__':

experiment/paper_past_snow_loads/result_trends_and_return_levels/study_visualizer_for_non_stationary_trends.py → experiment/paper_past_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -11,6 +11,8 @@ from experiment.meteo_france_data.plot.create_shifted_cmap import get_shifted_ma
 from experiment.meteo_france_data.scm_models_data.abstract_study import AbstractStudy
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     StudyVisualizer
+from experiment.paper_past_snow_loads.check_mcmc_convergence_for_return_levels.gelman_convergence_test import \
+    compute_gelman_convergence_value
 from experiment.trend_analysis.abstract_score import MeanScore
 from experiment.trend_analysis.univariate_test.abstract_gev_trend_test import AbstractGevTrendTest
 from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
@@ -147,7 +149,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def massif_name_to_tdrl_value(self):
         return {m: t.time_derivative_of_return_level for m, t in
                 self.massif_name_to_minimized_aic_non_stationary_trend_test.items()}
-    
+
     @cached_property
     def cmap(self):
         return get_shifted_map(-self._max_abs_tdrl, self._max_abs_tdrl)
@@ -169,8 +171,8 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     # Part 2 - Uncertainty return level plot
 
-    def massif_name_to_model_class(self, massif_name, non_stationary_model):
-        if not non_stationary_model:
+    def massif_name_and_non_stationary_context_to_model_class(self, massif_name, non_stationary_context):
+        if not non_stationary_context:
             return StationaryTemporalModel
         else:
             return self.massif_name_to_minimized_aic_non_stationary_trend_test[massif_name].unconstrained_model_class
@@ -179,16 +181,12 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def nb_contexts(self):
         return len(self.non_stationary_contexts)
 
-    @property
-    def nb_uncertainty_method(self):
-        return len(self.uncertainty_methods)
-
-    def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_model) \
+    def all_massif_name_to_eurocode_uncertainty_for_minimized_aic_model_class(self, ci_method, non_stationary_context) \
             -> Dict[str, Tuple[int, EurocodeConfidenceIntervalFromExtremes]]:
         # Compute for the uncertainty massif names
         arguments = [
             [self.massif_name_to_non_null_years_and_maxima[m],
-             self.massif_name_to_model_class(m, non_stationary_model),
+             self.massif_name_and_non_stationary_context_to_model_class(m, non_stationary_context),
              ci_method, self.effective_temporal_covariate,
              self.massif_name_to_eurocode_quantile_level_in_practice[m]
              ] for m in self.uncertainty_massif_names]
@@ -219,3 +217,15 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
 
     def model_name_to_uncertainty_method_to_ratio_above_eurocode(self):
         assert self.uncertainty_massif_names == self.study.study_massif_names
+
+    # Some checks with Gelman convergence diagnosis
+
+    def massif_name_to_gelman_convergence_value(self, mcmc_iterations, model_class, nb_chains):
+        arguments = [(self.massif_name_to_non_null_years_and_maxima[m], mcmc_iterations, model_class, nb_chains)
+                     for m in self.uncertainty_massif_names]
+        if self.multiprocessing:
+            with Pool(NB_CORES) as p:
+                res = p.starmap(compute_gelman_convergence_value, arguments)
+        else:
+            res = [compute_gelman_convergence_value(*argument) for argument in arguments]
+        return dict(zip(self.uncertainty_massif_names, res))

experiment/trend_analysis/univariate_test/utils.py

@@ -20,8 +20,8 @@ def load_temporal_coordinates_and_dataset(maxima, years):
     return coordinates, dataset
 
 
-def fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method):
-    model = model_class(coordinates, starting_point=starting_year, fit_method=fit_method)
+def fitted_linear_margin_estimator(model_class, coordinates, dataset, starting_year, fit_method, **model_kwargs):
+    model = model_class(coordinates, starting_point=starting_year, fit_method=fit_method, **model_kwargs)
     estimator = LinearMarginEstimator(dataset, model)
     estimator.fit()
     return estimator

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

@@ -24,8 +24,10 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
     """Linearity only with respect to the temporal coordinates"""
 
     def __init__(self, coordinates: AbstractCoordinates, use_start_value=False, params_start_fit=None,
-                 params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit):
+                 params_sample=None, starting_point=None, fit_method=TemporalMarginFitMethod.is_mev_gev_fit,
+                 nb_iterations_for_bayesian_fit=5000):
         super().__init__(coordinates, use_start_value, params_start_fit, params_sample, starting_point)
+        self.nb_iterations_for_bayesian_fit = nb_iterations_for_bayesian_fit
         assert isinstance(fit_method, TemporalMarginFitMethod)
         self.fit_method = fit_method
 
@@ -73,7 +75,7 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
                                    method='Bayesian',
                                    priorFun="fevdPriorCustom",
                                    priorParams=r.list(q=r.c(6), p=r.c(9)),
-                                   iter=5000,
+                                   iter=self.nb_iterations_for_bayesian_fit,
                                    **r_type_argument_kwargs
                                    )
         return ResultFromBayesianExtremes(res, self.margin_function_start_fit.gev_param_name_to_dims)