From 55041bcd3e2df555224c2b25cf5800b6434dcea1 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 7 Nov 2019 10:58:44 +0100
Subject: [PATCH] [Paper drawing] add drawing for the MCMC sampling

---
 .../scm_models_data/abstract_study.py         |   4 +
 .../__init__.py                               |   0
 .../main2_choice_to_not_use_starting_years.py |   2 +-
 ...main3_non_stationary_strength_evolution.py |   2 +-
 ..._spatial_altitude_starting_years_impact.py |   2 +-
 ...n4_common_spatial_starting_years_impact.py |   2 +-
 .../main4_individual_starting_years_impact.py |   2 +-
 .../{paper1_steps => paper1}/__init__.py      |   0
 experiment/paper1/main_drawing_bayesian.py    |  94 ++++++++++++++++++
 .../poster_EVAN2019/__init__.py               |   0
 .../poster_EVAN2019/main_poster_EVAN2019.py   |   2 +-
 .../shape_prior_check/__init__.py             |   0
 .../analyse_shape_from_some_experiment.py     |   0
 .../shape_prior_check/hist_values_shape.png   | Bin
 .../shape_from_some_experiment.txt            |   0
 .../shape_prior_check/some_experiment_EVAN.py |   2 +-
 ...snow_load_trends_percentage_significant.py |   0
 experiment/{paper1_steps => paper1}/utils.py  |   0
 .../validations/__init__.py                   |   0
 .../main0_comparison_with_observations.py     |   0
 .../main1_good_stationary_gev_fit.py          |   0
 21 files changed, 105 insertions(+), 7 deletions(-)
 rename experiment/{paper1_steps => paper1}/1 - non stationary model choice/__init__.py (100%)
 rename experiment/{paper1_steps => paper1}/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py (94%)
 rename experiment/{paper1_steps => paper1}/1 - non stationary model choice/main3_non_stationary_strength_evolution.py (95%)
 rename experiment/{paper1_steps => paper1}/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py (97%)
 rename experiment/{paper1_steps => paper1}/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py (97%)
 rename experiment/{paper1_steps => paper1}/1 - non stationary model choice/main4_individual_starting_years_impact.py (96%)
 rename experiment/{paper1_steps => paper1}/__init__.py (100%)
 create mode 100644 experiment/paper1/main_drawing_bayesian.py
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/__init__.py (100%)
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/main_poster_EVAN2019.py (99%)
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/shape_prior_check/__init__.py (100%)
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py (100%)
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/shape_prior_check/hist_values_shape.png (100%)
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt (100%)
 rename experiment/{paper1_steps => paper1}/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py (96%)
 rename experiment/{paper1_steps => paper1}/snow_load_trends_percentage_significant.py (100%)
 rename experiment/{paper1_steps => paper1}/utils.py (100%)
 rename experiment/{paper1_steps => paper1}/validations/__init__.py (100%)
 rename experiment/{paper1_steps => paper1}/validations/main0_comparison_with_observations.py (100%)
 rename experiment/{paper1_steps => paper1}/validations/main1_good_stationary_gev_fit.py (100%)

diff --git a/experiment/meteo_france_data/scm_models_data/abstract_study.py b/experiment/meteo_france_data/scm_models_data/abstract_study.py
index b16ead03..8c407129 100644
--- a/experiment/meteo_france_data/scm_models_data/abstract_study.py
+++ b/experiment/meteo_france_data/scm_models_data/abstract_study.py
@@ -102,6 +102,10 @@ class AbstractStudy(object):
     def observations_annual_maxima(self) -> AnnualMaxima:
         return AnnualMaxima(df_maxima_gev=pd.DataFrame(self.year_to_annual_maxima, index=self.study_massif_names))
 
+    def annual_maxima_and_years(self, massif_name) -> Tuple[np.ndarray, np.ndarray]:
+        df = self.observations_annual_maxima.df_maxima_gev
+        return df.loc[massif_name].values, np.array(df.columns)
+
     @cached_property
     def year_to_annual_maxima(self) -> OrderedDict:
         # Map each year to an array of size nb_massif
diff --git a/experiment/paper1_steps/1 - non stationary model choice/__init__.py b/experiment/paper1/1 - non stationary model choice/__init__.py
similarity index 100%
rename from experiment/paper1_steps/1 - non stationary model choice/__init__.py
rename to experiment/paper1/1 - non stationary model choice/__init__.py
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py b/experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py
similarity index 94%
rename from experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py
rename to experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py
index c43a29e7..b5a2f7c5 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py	
+++ b/experiment/paper1/1 - non stationary model choice/main2_choice_to_not_use_starting_years.py	
@@ -2,7 +2,7 @@ import time
 
 from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
     Altitude_Hypercube_Year_Visualizer
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_repartition():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py b/experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py
similarity index 95%
rename from experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py
rename to experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py
index 29713ba1..6b174865 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main3_non_stationary_strength_evolution.py	
+++ b/experiment/paper1/1 - non stationary model choice/main3_non_stationary_strength_evolution.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
similarity index 97%
rename from experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
rename to experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py
index ac8a63e3..10d4fe7b 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py	
+++ b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_altitude_starting_years_impact.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py b/experiment/paper1/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_spatial_starting_years_impact.py
rename to experiment/paper1/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py
index e4e0248f..2973ad0f 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py	
+++ b/experiment/paper1/1 - non stationary model choice/main4_common_spatial_starting_years_impact.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py b/experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py
similarity index 96%
rename from experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py
rename to experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py
index c8a3fbc8..ada1e073 100644
--- a/experiment/paper1_steps/1 - non stationary model choice/main4_individual_starting_years_impact.py	
+++ b/experiment/paper1/1 - non stationary model choice/main4_individual_starting_years_impact.py	
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.gev_trend_test_two_parameters imp
 """
 Visualize the 0.99 quantile initial value and its evolution
 """
-from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
+from experiment.paper1.utils import get_full_altitude_visualizer, FULL_ALTITUDES
 
 
 def main_fast_spatial_risk_evolution():
diff --git a/experiment/paper1_steps/__init__.py b/experiment/paper1/__init__.py
similarity index 100%
rename from experiment/paper1_steps/__init__.py
rename to experiment/paper1/__init__.py
diff --git a/experiment/paper1/main_drawing_bayesian.py b/experiment/paper1/main_drawing_bayesian.py
new file mode 100644
index 00000000..713147f2
--- /dev/null
+++ b/experiment/paper1/main_drawing_bayesian.py
@@ -0,0 +1,94 @@
+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.visualization.utils import create_adjusted_axes
+from experiment.trend_analysis.univariate_test.utils import load_temporal_coordinates_and_dataset, \
+    fitted_linear_margin_estimator
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.linear_margin_model.abstract_temporal_linear_margin_model import \
+    TemporalMarginFitMethod
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
+    ExtractEurocodeReturnLevelFromMyBayesianExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+
+
+def main_drawing_bayesian():
+    return_level_bayesian = get_return_level_bayesian_example()
+    print(return_level_bayesian.result_from_fit.df_all_samples)
+    print(return_level_bayesian.result_from_fit.df_posterior_samples)
+    print(return_level_bayesian.posterior_eurocode_return_level_samples_for_temporal_covariate)
+    axes = create_adjusted_axes(1, 3)
+    # Plot the trajectories on the first axis
+    ax_trajectories = axes[0]
+    ax_trajectories_inverted = ax_trajectories.twinx()
+    df_all_samples = return_level_bayesian.result_from_fit.df_all_samples
+    iteration_step = [i + 1 for i in df_all_samples.index]
+    lns = []
+    # Last color is for the return level
+    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]))
+    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]
+        ax = gev_param_name_to_ax[gev_param_name]
+        color = gev_param_name_to_color[gev_param_name]
+        ln = ax.plot(iteration_step, df_all_samples.iloc[:, j].values, label=label, color=color)
+        lns.extend(ln)
+    ax_trajectories_inverted.set_ylim(-0.3, 10)
+    ax_trajectories.set_ylim(-0.3, 90)
+    for ax in [ax_trajectories, ax_trajectories_inverted]:
+        ax.set_xlim(min(iteration_step), max(iteration_step))
+    # labs = [l.get_label() for l in lns]
+    # ax_trajectories.legend(lns, labs, loc=0)
+    ax.axvline(x=return_level_bayesian.result_from_fit.burn_in_nb, color='k', linestyle='--')
+    ax_trajectories_inverted.legend(loc=1)
+    ax_trajectories.legend(loc=2)
+    ax_trajectories.set_xlabel("MCMC iterations")
+
+    # Plot the parameter posterior on axes 1
+    ax_parameter_posterior = axes[1]
+    ax_parameter_posterior_flip = ax_parameter_posterior.twiny()
+    gev_param_name_to_ax = dict(
+        zip(GevParams.PARAM_NAMES, [ax_parameter_posterior, ax_parameter_posterior.twiny(), ax_parameter_posterior_flip]))
+    df_posterior_samples = return_level_bayesian.result_from_fit.df_posterior_samples
+    lns = []
+    for j, gev_param_name in enumerate(GevParams.PARAM_NAMES[:]):
+        label = gev_param_name_to_label[gev_param_name]
+        color = gev_param_name_to_color[gev_param_name]
+        ax = gev_param_name_to_ax[gev_param_name]
+        ln = sns.kdeplot(df_posterior_samples.iloc[:, j], ax=ax, label=label, color=color)
+        lns.append(ln)
+    labs = [l.get_label() for l in lns]
+    ax_parameter_posterior.legend(lns, labs, loc=0)
+
+    # Plot the return level posterior on the last axes
+    ax_return_level_posterior = axes[2]
+    sns.kdeplot(return_level_bayesian.posterior_eurocode_return_level_samples_for_temporal_covariate,
+                ax=ax_return_level_posterior, color=colors[-1])
+    ax_return_level_posterior.set_xlabel("$z_p(\\theta)")
+    ax_return_level_posterior.set_ylabel("$p(z_p(\\theta)|y)")
+
+    plt.show()
+
+
+def get_return_level_bayesian_example():
+    maxima, years = CrocusSwe3Days(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)
+    return_level_bayesian = ExtractEurocodeReturnLevelFromMyBayesianExtremes(estimator=fitted_estimator,
+                                                                             ci_method=ConfidenceIntervalMethodFromExtremes.my_bayes,
+                                                                             temporal_covariate=2017)
+    return return_level_bayesian
+
+
+if __name__ == '__main__':
+    main_drawing_bayesian()
diff --git a/experiment/paper1_steps/poster_EVAN2019/__init__.py b/experiment/paper1/poster_EVAN2019/__init__.py
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/__init__.py
rename to experiment/paper1/poster_EVAN2019/__init__.py
diff --git a/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py b/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
similarity index 99%
rename from experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
rename to experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
index 063783cc..ac9ddb19 100644
--- a/experiment/paper1_steps/poster_EVAN2019/main_poster_EVAN2019.py
+++ b/experiment/paper1/poster_EVAN2019/main_poster_EVAN2019.py
@@ -7,7 +7,7 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
 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
-from experiment.paper1_steps.utils import get_full_altitude_visualizer
+from experiment.paper1.utils import get_full_altitude_visualizer
 
 POSTER_ALTITUDES = [900, 1800, 2700]
 import matplotlib as mpl
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/__init__.py b/experiment/paper1/poster_EVAN2019/shape_prior_check/__init__.py
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/__init__.py
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/__init__.py
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py b/experiment/paper1/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/analyse_shape_from_some_experiment.py
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/hist_values_shape.png b/experiment/paper1/poster_EVAN2019/shape_prior_check/hist_values_shape.png
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/hist_values_shape.png
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/hist_values_shape.png
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt b/experiment/paper1/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt
similarity index 100%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/shape_from_some_experiment.txt
diff --git a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py b/experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
similarity index 96%
rename from experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
rename to experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
index d93577b4..9d0175bf 100644
--- a/experiment/paper1_steps/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
+++ b/experiment/paper1/poster_EVAN2019/shape_prior_check/some_experiment_EVAN.py
@@ -9,7 +9,7 @@ from experiment.trend_analysis.univariate_test.abstract_comparison_non_stationar
 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
-from experiment.paper1_steps.utils import get_full_altitude_visualizer
+from experiment.paper1.utils import get_full_altitude_visualizer
 
 POSTER_ALTITUDES = [900, 1800, 2700]
 import matplotlib as mpl
diff --git a/experiment/paper1_steps/snow_load_trends_percentage_significant.py b/experiment/paper1/snow_load_trends_percentage_significant.py
similarity index 100%
rename from experiment/paper1_steps/snow_load_trends_percentage_significant.py
rename to experiment/paper1/snow_load_trends_percentage_significant.py
diff --git a/experiment/paper1_steps/utils.py b/experiment/paper1/utils.py
similarity index 100%
rename from experiment/paper1_steps/utils.py
rename to experiment/paper1/utils.py
diff --git a/experiment/paper1_steps/validations/__init__.py b/experiment/paper1/validations/__init__.py
similarity index 100%
rename from experiment/paper1_steps/validations/__init__.py
rename to experiment/paper1/validations/__init__.py
diff --git a/experiment/paper1_steps/validations/main0_comparison_with_observations.py b/experiment/paper1/validations/main0_comparison_with_observations.py
similarity index 100%
rename from experiment/paper1_steps/validations/main0_comparison_with_observations.py
rename to experiment/paper1/validations/main0_comparison_with_observations.py
diff --git a/experiment/paper1_steps/validations/main1_good_stationary_gev_fit.py b/experiment/paper1/validations/main1_good_stationary_gev_fit.py
similarity index 100%
rename from experiment/paper1_steps/validations/main1_good_stationary_gev_fit.py
rename to experiment/paper1/validations/main1_good_stationary_gev_fit.py
-- 
GitLab