[PAPER1][STARTING YEAR] write main for location & scale with common starting...

[PAPER1][STARTING YEAR] write main for location & scale with common starting year spatially & altitude

experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py

@@ -25,7 +25,10 @@ class AbstractHypercubeVisualizer(object):
                  first_starting_year=None,
                  last_starting_year=None,
                  exact_starting_year=None,
-                 verbose=True):
+                 verbose=True,
+                 sigma_for_best_year=0.0):
+        assert sigma_for_best_year >= 0.0
+        self.sigma_for_best_year = sigma_for_best_year
         self.reduce_strength_array = reduce_strength_array
         self.verbose = verbose
         self.save_to_file = save_to_file

experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py

@@ -2,6 +2,7 @@ import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
 from matplotlib.ticker import FormatStrFormatter, ScalarFormatter
+from scipy.ndimage import gaussian_filter
 
 from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.abstract_hypercube_visualizer import \
     AbstractHypercubeVisualizer
@@ -157,7 +158,8 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
             assert isinstance(serie, pd.Series)
             xlabel_values = list(serie.index)
             values = list(serie.values)
-            argmax_idx = np.argmax(values)
+            smooth_values = gaussian_filter(values, self.sigma_for_best_year)
+            argmax_idx = np.argmax(smooth_values)
             best_year = xlabel_values[argmax_idx]
             if plot_title is not None:
                 plot_title += ' (max reached in {})'.format(best_year)
@@ -168,6 +170,9 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
                 ax.plot([], [], label=ylabel, color=color)
                 linewidth = 10 if poster_plot else None
                 ax_reversed.plot(xlabel_values, values, label=ylabel, color=color, linewidth=linewidth)
+                if self.sigma_for_best_year > 0:
+                    ax_reversed.plot(xlabel_values, smooth_values, label=ylabel + ' smooth', color=color, linestyle=':',
+                                     linewidth=linewidth)
                 fontsize = 30 if poster_plot else None
                 ax_reversed.set_ylabel(ylabel, color=color, fontsize=fontsize, labelpad=-20)
                 ax_reversed.axvline(x=best_year, color=color, linestyle='--', linewidth=linewidth)
@@ -344,6 +349,8 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
         if self.reduce_strength_array:
             title += '\nEvolution of the quantile {} every {} years'.format(AbstractGevTrendTest.quantile_for_strength,
                                                                             AbstractGevTrendTest.nb_years_for_quantile_evolution)
+        else:
+            title += '\nStarting years'
         if set:
             plt.suptitle(title)
         return title

experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py

+import time
+
+from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
+    Altitude_Hypercube_Year_Visualizer, AltitudeHypercubeVisualizerWithoutTrendType
+from experiment.trend_analysis.univariate_test.gev_trend_test_one_parameter import GevScaleTrendTest, \
+    GevLocationTrendTest
+from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters import GevLocationAndScaleTrendTest
+
+"""
+Visualize the 0.99 quantile initial value and its evolution
+"""
+from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+
+
+def main_fast_spatial_risk_evolution():
+    vizualiser = get_full_altitude_visualizer(AltitudeHypercubeVisualizerWithoutTrendType, altitude=None,
+                                              reduce_strength_array=True,
+                                              trend_test_class=GevLocationAndScaleTrendTest,
+                                              offset_starting_year=28)
+    vizualiser.save_to_file = False
+    vizualiser.sigma_for_best_year = 1.0
+    res = vizualiser.visualize_year_trend_test(subtitle_specified='CrocusSwe3Days')
+    print(res)
+
+
+def main_full_spatial_risk_evolution():
+    for trend_test_class in [GevLocationAndScaleTrendTest]:
+        # Compare the risk with and without taking into account the starting year
+        vizualiser = get_full_altitude_visualizer(AltitudeHypercubeVisualizerWithoutTrendType, altitude=None,
+                                                  reduce_strength_array=True,
+                                                  trend_test_class=trend_test_class,
+                                                  offset_starting_year=20)
+        vizualiser.sigma_for_best_year = 1.0
+        res = vizualiser.visualize_year_trend_test(subtitle_specified='CrocusSwe3Days')
+        best_year = res[0][1]
+        for altitude in FULL_ALTITUDES[:]:
+            # Starting Year=1958
+            vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+                                                      exact_starting_year=1958, reduce_strength_array=True,
+                                                      trend_test_class=trend_test_class)
+            vizualiser.visualize_massif_trend_test_one_altitude()
+            # Optimal common starting year
+            vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
+                                                      exact_starting_year=best_year, reduce_strength_array=True,
+                                                      trend_test_class=trend_test_class)
+            vizualiser.visualize_massif_trend_test_one_altitude()
+
+
+def main_run():
+    main_full_spatial_risk_evolution()
+    # main_fast_spatial_risk_evolution()
+
+
+if __name__ == '__main__':
+    start = time.time()
+    main_run()
+    duration = time.time() - start
+    print('Full run took {}s'.format(round(duration, 1)))

experiment/paper1_steps/1 - non stationary model choice/main4_common_starting_years_impact.py → experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py

@@ -29,7 +29,7 @@ def main_fast_spatial_risk_evolution():
 def main_full_spatial_risk_evolution():
     # Compare the risk with and without taking into account the starting year
     for altitude in FULL_ALTITUDES[:]:
-        for trend_test_class in [GevLocationAndScaleTrendTest]:
+        for trend_test_class in [GevScaleTrendTest]:
             # Starting Year=1958
             vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
                                                       exact_starting_year=1958, reduce_strength_array=True,

experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py

@@ -16,7 +16,8 @@ def main_fast_spatial_risk_evolution():
     for altitude in [1800]:
         vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
                                                   reduce_strength_array=False,
-                                                  trend_test_class=GevLocationAndScaleTrendTest)
+                                                  trend_test_class=GevLocationAndScaleTrendTest,
+                                                  offset_starting_year=20)
         vizualiser.save_to_file = False
         vizualiser.visualize_massif_trend_test_one_altitude()
         vizualiser.reduce_strength_array = True
@@ -25,7 +26,7 @@ def main_fast_spatial_risk_evolution():
 
 def main_full_spatial_risk_evolution():
     # Compare the risk with and without taking into account the starting year
-    for altitude in FULL_ALTITUDES[:]:
+    for altitude in FULL_ALTITUDES[-2:]:
         for trend_test_class in [GevLocationAndScaleTrendTest]:
             vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
                                                       exact_starting_year=1958, reduce_strength_array=True,
@@ -33,7 +34,8 @@ def main_full_spatial_risk_evolution():
             vizualiser.visualize_massif_trend_test_one_altitude()
             vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
                                                       reduce_strength_array=True,
-                                                      trend_test_class=trend_test_class)
+                                                      trend_test_class=trend_test_class,
+                                                      offset_starting_year=20)
             vizualiser.visualize_massif_trend_test_one_altitude()
             vizualiser.reduce_strength_array = False
             vizualiser.visualize_massif_trend_test_one_altitude()