From 8b581b73ccc70d3918b57271ffe1b56b98402a85 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 29 Jul 2019 18:37:04 +0200
Subject: [PATCH] [PAPER1][STARTING YEAR] write main for location & scale with
 common starting year spatially & altitude

---
 .../abstract_hypercube_visualizer.py          |  5 +-
 .../altitude_hypercube_visualizer.py          |  9 ++-
 ..._spatial_altitude_starting_years_impact.py | 58 +++++++++++++++++++
 ...4_common_spatial_starting_years_impact.py} |  2 +-
 .../main4_individual_starting_years_impact.py |  8 ++-
 5 files changed, 76 insertions(+), 6 deletions(-)
 create mode 100644 experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
 rename experiment/paper1_steps/1 - non stationary model choice/{main4_common_starting_years_impact.py => main4_common_spatial_starting_years_impact.py} (97%)

diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
index 8bd9527e..728090c1 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
+++ b/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
diff --git a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py b/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
index 5db8641c..eba9f36a 100644
--- a/experiment/meteo_france_data/scm_models_data/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
+++ b/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
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py b/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
new file mode 100644
index 00000000..ac8a63e3
--- /dev/null
+++ b/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py	
@@ -0,0 +1,58 @@
+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)))
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_common_starting_years_impact.py b/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py
similarity index 97%
rename from experiment/paper1_steps/1 - non stationary model choice/main4_common_starting_years_impact.py
rename to experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py
index 3f08aa4a..e07bdd64 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_common_starting_years_impact.py	
+++ b/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,
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py b/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py
index 3c051130..22368c51 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py	
+++ b/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()
-- 
GitLab