[paper 1] add some options on study visualizer to handle years without snow

papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/main_qqplot_for_big_shapes.py

@@ -17,6 +17,8 @@ def get_tuple_ordered_by_shape(fast=False):
                                                                               select_only_acceptable_shape_parameter=False,
                                                                               multiprocessing=True,
                                                                               save_to_file=True,
+                                                                              fit_gev_only_on_non_null_maxima=False,
+                                                                              fit_only_time_series_with_ninety_percent_of_non_null_values=False,
                                                                               show=False)
                               for altitude in altitudes}
     # Extract all the values

papers/exceeding_snow_loads/check_mle_convergence_for_trends/qqplot/plot_qqplot.py

@@ -7,6 +7,8 @@ import pandas as pd
 from matplotlib.ticker import PercentFormatter
 
 from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusSnowLoadTotal
+from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
+    ALL_ALTITUDES_WITHOUT_NAN
 from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
     TemporalMarginFitMethod
 from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
@@ -183,20 +185,21 @@ if __name__ == '__main__':
     300 Haut_Var-Haut_Verdon 0.8446498197950775
 
     """
-    altitudes = [300, 600, 900, 1200, 1500, 1800][:]
+    altitudes = [300, 600, 900, 1200, 1500, 1800][:3]
     # altitudes = ALL_ALTITUDES_WITHOUT_NAN
     # altitudes = [900, 1800, 2700]
     altitude_to_visualizer = {altitude: StudyVisualizerForNonStationaryTrends(CrocusSnowLoadTotal(altitude=altitude),
                                                                               select_only_acceptable_shape_parameter=True,
                                                                               fit_method=TemporalMarginFitMethod.extremes_fevd_mle,
                                                                               multiprocessing=True,
+                                                                              fit_gev_only_on_non_null_maxima=True,
                                                                               save_to_file=True,
                                                                               show=False)
                               for altitude in altitudes}
 
     # plot_qqplot_wrt_standard_gumbel(altitude_to_visualizer)
-    # plot_hist_psnow(altitude_to_visualizer)
-    plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(altitude_to_visualizer, nb_worst_examples=None)
+    plot_hist_psnow(altitude_to_visualizer)
+    # plot_intensity_against_gumbel_quantile_for_time_series_with_missing_zeros(altitude_to_visualizer, nb_worst_examples=None)
     # plot_exceedance_psnow(altitude_to_visualizer)
     # non_stationarity_psnow(altitude_to_visualizer)
 

papers/exceeding_snow_loads/check_mle_convergence_for_trends/without_maximum/study_visualizer_for_fit_witout_maximum.py

@@ -18,7 +18,7 @@ class StudyVisualizerForFitWithoutMaximum(StudyVisualizerForNonStationaryTrends)
 
     @cached_property
     def massif_name_to_maximum_index_for_non_null_values(self) -> Tuple[Dict, Dict]:
-        d = super().massif_name_to_non_null_years_and_maxima
+        d = super().massif_name_to_years_and_maxima_for_model_fitting
         d_without_maximum = {}
         d_maximum = {}
         for m, (years, maxima) in d.items():
@@ -35,7 +35,7 @@ class StudyVisualizerForFitWithoutMaximum(StudyVisualizerForNonStationaryTrends)
         return self.massif_name_to_maximum_index_for_non_null_values[1]
 
     @cached_property
-    def massif_name_to_non_null_years_and_maxima(self):
+    def massif_name_to_years_and_maxima_for_model_fitting(self):
         return self.massif_name_to_maximum_index_for_non_null_values[0]
 
     def maximum_value_test(self):

papers/exceeding_snow_loads/study_visualizer_for_non_stationary_trends.py

@@ -52,12 +52,17 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
                  relative_change_trend_plot=True,
                  non_stationary_trend_test_to_marker=None,
                  fit_method=TemporalMarginFitMethod.extremes_fevd_mle,
-                 select_only_acceptable_shape_parameter=True):
+                 select_only_acceptable_shape_parameter=True,
+                 fit_gev_only_on_non_null_maxima=True,
+                 fit_only_time_series_with_ninety_percent_of_non_null_values=True,
+                 ):
         super().__init__(study, show, save_to_file, only_one_graph, only_first_row, vertical_kde_plot,
                          year_for_kde_plot, plot_block_maxima_quantiles, temporal_non_stationarity,
                          transformation_class, verbose, multiprocessing, complete_non_stationary_trend_analysis,
                          normalization_under_one_observations, score_class)
         # Add some attributes
+        self.fit_only_time_series_with_ninety_percent_of_non_null_values = fit_only_time_series_with_ninety_percent_of_non_null_values
+        self.fit_gev_only_on_non_null_maxima = fit_gev_only_on_non_null_maxima
         self.select_only_acceptable_shape_parameter = select_only_acceptable_shape_parameter
         self.fit_method = fit_method
         self.non_stationary_trend_test_to_marker = non_stationary_trend_test_to_marker
@@ -107,15 +112,22 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     def massif_name_to_eurocode_quantile_level_in_practice(self):
         """Due to missing data, the the eurocode quantile which 0.98 if we have all the data
         correspond in practice to the quantile psnow x 0.98 of the data where there is snow"""
-        return {m: 1 - ((1 - EUROCODE_QUANTILE) / p_snow) for m, p_snow in self.massif_name_to_psnow.items()}
+        if self.fit_gev_only_on_non_null_maxima:
+            return {m: 1 - ((1 - EUROCODE_QUANTILE) / p_snow) for m, p_snow in self.massif_name_to_psnow.items()}
+        else:
+            return {m: EUROCODE_QUANTILE for m in self.massif_name_to_psnow.keys()}
 
     @cached_property
-    def massif_name_to_non_null_years_and_maxima(self):
-        # return self.massif_name_to_years_and_maxima
-        d = {}
-        for m, (years, maxima) in self.massif_name_to_years_and_maxima.items():
-            mask = np.nonzero(maxima)
-            d[m] = (years[mask], maxima[mask])
+    def massif_name_to_years_and_maxima_for_model_fitting(self):
+        if self.fit_gev_only_on_non_null_maxima:
+            d = {}
+            for m, (years, maxima) in self.massif_name_to_years_and_maxima.items():
+                mask = np.nonzero(maxima)
+                d[m] = (years[mask], maxima[mask])
+        else:
+            d = self.massif_name_to_years_and_maxima
+            if self.fit_only_time_series_with_ninety_percent_of_non_null_values:
+                d = {m: v for m, v in d.items() if self.massif_name_to_psnow[m] >= 0.90}
         return d
 
     @property
@@ -137,7 +149,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
         massif_name_to_trend_test_that_minimized_aic = {}
         massif_name_to_stationary_trend_test_that_minimized_aic = {}
         massif_name_to_gumbel_trend_test_that_minimized_aic = {}
-        for massif_name, (x, y) in self.massif_name_to_non_null_years_and_maxima.items():
+        for massif_name, (x, y) in self.massif_name_to_years_and_maxima_for_model_fitting.items():
             quantile_level = self.massif_name_to_eurocode_quantile_level_in_practice[massif_name]
             all_trend_test = [
                 t(years=x, maxima=y, starting_year=starting_year, quantile_level=quantile_level,
@@ -317,7 +329,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
             -> Dict[str, EurocodeConfidenceIntervalFromExtremes]:
         # Compute for the uncertainty massif names
         arguments = [
-            [self.massif_name_to_non_null_years_and_maxima[m],
+            [self.massif_name_to_years_and_maxima_for_model_fitting[m],
              self.massif_name_and_model_subset_to_model_class(m, model_subset_for_uncertainty),
              ci_method, self.effective_temporal_covariate,
              self.massif_name_to_eurocode_quantile_level_in_practice[m]
@@ -353,7 +365,7 @@ class StudyVisualizerForNonStationaryTrends(StudyVisualizer):
     # 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)
+        arguments = [(self.massif_name_to_years_and_maxima_for_model_fitting[m], mcmc_iterations, model_class, nb_chains)
                      for m in self.uncertainty_massif_names]
         if self.multiprocessing:
             with Pool(NB_CORES) as p: