[SCM][NON STATIONARITY] add non stationary trend object. to create specific plot for this topic

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

@@ -95,9 +95,18 @@ def complete_analysis(only_first_one=False):
                 study_visualizer.visualize_linear_margin_fit()
 
 
+def trend_analysis():
+    save_to_file = False
+    only_first_one = True
+    for study_class in [CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][1:2]:
+        for study in study_iterator(study_class, only_first_one=only_first_one):
+            study_visualizer = StudyVisualizer(study, save_to_file=save_to_file)
+            study_visualizer.visualize_temporal_trend_relevance(complete_analysis=False)
+
 if __name__ == '__main__':
     # annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2100)
-    normal_visualization(temporal_non_stationarity=True)
+    # normal_visualization(temporal_non_stationarity=True)
+    trend_analysis()
     # max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1)
     # extended_visualization()
     # complete_analysis()

experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py

+from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
+from scipy.stats import chi2
+from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
+    FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_estimator.extreme_models.margin_model.linear_margin_model import \
+    LinearAllParametersTwoFirstCoordinatesMarginModel, LinearAllTwoStatialCoordinatesLocationLinearMarginModel, \
+    LinearStationaryMarginModel, LinearNonStationaryLocationMarginModel
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+
+
+class AbstractNonStationaryTrendTest(object):
+    RESULT_ATTRIBUTE_METRIC = 'deviance'
+
+    def __init__(self, dataset: AbstractDataset, estimator_class,
+                 stationary_margin_model_class, non_stationary_margin_model_class):
+        self.dataset = dataset
+        self.estimator_class = estimator_class
+        self.stationary_margin_model_class = stationary_margin_model_class
+        self.non_stationary_margin_model_class = non_stationary_margin_model_class
+
+    def load_estimator(self, margin_model) -> AbstractEstimator:
+        return self.estimator_class(self.dataset, margin_model)
+
+    def get_metric(self, margin_model_class, starting_point):
+        margin_model = margin_model_class(coordinates=self.dataset.coordinates, starting_point=starting_point)
+        estimator = self.load_estimator(margin_model)   # type: AbstractEstimator
+        estimator.fit()
+        metric = estimator.result_from_fit.__getattribute__(self.RESULT_ATTRIBUTE_METRIC)
+        assert isinstance(metric, float)
+        return metric
+
+    def visualize(self, ax, complete_analysis=True):
+        # Define the year_min and year_max for the starting point
+        if complete_analysis:
+            year_min, year_max, step = 1960, 1990, 1
+        else:
+            year_min, year_max, step = 1960, 1990, 10
+        # Fit the stationary model
+        stationary_metric = self.get_metric(self.stationary_margin_model_class, starting_point=None)
+        # Fit the non stationary model
+        years = list(range(year_min, year_max + 1, step))
+        non_stationary_metrics = [self.get_metric(self.non_stationary_margin_model_class, year) for year in years]
+        difference = [m - stationary_metric for m in non_stationary_metrics]
+        # Plot some lines
+        ax.axhline(y=0, xmin=year_min, xmax=year_max)
+        # Significative line
+        significative_deviance = chi2.ppf(q=0.95, df=1)
+        ax.axhline(y=significative_deviance, xmin=year_min, xmax=year_max)
+        # todo: plot the line corresponding to the significance 0.05
+        ax.plot(years, difference, 'ro-')
+
+
+class IndependenceLocationTrendTest(AbstractNonStationaryTrendTest):
+
+    def __init__(self, dataset, coordinate_idx):
+        pass
+
+
+class ConditionalIndedendenceLocationTrendTest(AbstractNonStationaryTrendTest):
+
+    def __init__(self, dataset):
+        super().__init__(dataset=dataset,
+                         estimator_class=LinearMarginEstimator,
+                         stationary_margin_model_class=LinearStationaryMarginModel,
+                         non_stationary_margin_model_class=LinearNonStationaryLocationMarginModel)
+
+
+class MaxStableLocationTrendTest(AbstractNonStationaryTrendTest):
+
+    def __init__(self, dataset, max_stable_model):
+        super().__init__(dataset=dataset,
+                         estimator_class=FullEstimatorInASingleStepWithSmoothMargin,
+                         stationary_margin_model_class=LinearStationaryMarginModel,
+                         non_stationary_margin_model_class=LinearNonStationaryLocationMarginModel)
+        self.max_stable_model = max_stable_model
+
+    def load_estimator(self, margin_model) -> AbstractEstimator:
+        return self.estimator_class(self.dataset, margin_model, self.max_stable_model)

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

@@ -9,13 +9,15 @@ import pandas as pd
 import seaborn as sns
 
 from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
+from experiment.meteo_france_SCM_study.visualization.study_visualization.non_stationary_trends import \
+    ConditionalIndedendenceLocationTrendTest, MaxStableLocationTrendTest
 from experiment.meteo_france_SCM_study.visualization.utils import create_adjusted_axes
 from experiment.utils import average_smoothing_with_sliding_window
 from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
     FullEstimatorInASingleStepWithSmoothMargin
 from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
 from extreme_estimator.extreme_models.margin_model.linear_margin_model import LinearAllParametersAllDimsMarginModel, \
-    LinearNonStationaryMarginModel, LinearStationaryMarginModel
+    LinearNonStationaryLocationMarginModel, LinearStationaryMarginModel
 from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import \
     AbstractMarginFunction
 from extreme_estimator.extreme_models.margin_model.param_function.param_function import AbstractParamFunction
@@ -51,6 +53,8 @@ class StudyVisualizer(object):
         self._coordinates = None
         self._observations = None
 
+        self.default_covariance_function = CovarianceFunction.powexp
+
         # KDE PLOT ARGUMENTS
         self.vertical_kde_plot = vertical_kde_plot
         self.year_for_kde_plot = year_for_kde_plot
@@ -132,6 +136,25 @@ class StudyVisualizer(object):
             'for the year {}'.format(self.year_for_kde_plot)
         self.show_or_save_to_file()
 
+    def visualize_temporal_trend_relevance(self, complete_analysis):
+        self.temporal_non_stationarity = True
+        trend_tests = [ConditionalIndedendenceLocationTrendTest(self.dataset)]
+
+        max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
+        for max_stable_model in max_stable_models[:1]:
+            trend_tests.append(MaxStableLocationTrendTest(self.dataset, max_stable_model))
+
+        nb_trend_tests = len(trend_tests)
+        fig, axes = plt.subplots(1, nb_trend_tests, figsize=self.figsize)
+        if nb_trend_tests == 1:
+            axes = [axes]
+        fig.subplots_adjust(hspace=self.subplot_space, wspace=self.subplot_space)
+        for ax, trend_test in zip(axes, trend_tests):
+            trend_test.visualize(ax, complete_analysis=complete_analysis)
+
+        self.plot_name = 'trend tests'
+        self.show_or_save_to_file()
+
     def visualize_experimental_law(self, ax, massif_id):
         # Display the experimental law for a given massif
         all_massif_data = self.get_all_massif_data(massif_id)
@@ -248,15 +271,14 @@ class StudyVisualizer(object):
         pass
 
     def visualize_linear_margin_fit(self, only_first_max_stable=False):
-        default_covariance_function = CovarianceFunction.powexp
-        margin_class = LinearNonStationaryMarginModel if self.temporal_non_stationarity else LinearStationaryMarginModel
+        margin_class = LinearNonStationaryLocationMarginModel if self.temporal_non_stationarity else LinearStationaryMarginModel
         plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
         plot_name += '\n(with {} covariance structure when a covariance is needed)'.format(
-            str(default_covariance_function).split('.')[-1])
+            str(self.default_covariance_function).split('.')[-1])
         self.plot_name = plot_name
 
         # Load max stable models
-        max_stable_models = load_test_max_stable_models(default_covariance_function=default_covariance_function)
+        max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
         if only_first_max_stable:
             # Keep only the BrownResnick model
             max_stable_models = max_stable_models[1:2]

extreme_estimator/extreme_models/margin_model/linear_margin_model.py

@@ -127,5 +127,5 @@ class LinearStationaryMarginModel(LinearAllParametersTwoFirstCoordinatesMarginMo
     pass
 
 
-class LinearNonStationaryMarginModel(LinearAllTwoStatialCoordinatesLocationLinearMarginModel):
+class LinearNonStationaryLocationMarginModel(LinearAllTwoStatialCoordinatesLocationLinearMarginModel):
     pass

extreme_estimator/extreme_models/margin_model/margin_function/parametric_margin_function.py

@@ -62,7 +62,6 @@ class ParametricMarginFunction(IndependentMarginFunction):
         raise NotImplementedError
 
     def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
-        print('here get gev', self.starting_point)
         if self.starting_point is not None:
             # Shift temporal coordinate to enable to model temporal trend with starting point
             assert self.coordinates.has_temporal_coordinates

extreme_estimator/extreme_models/result_from_fit.py

 from typing import Dict
 
+import numpy as np
 from rpy2 import robjects
 
 from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
@@ -9,6 +10,7 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 
 class ResultFromFit(object):
 
+
     def __init__(self, result_from_fit: robjects.ListVector) -> None:
         if hasattr(result_from_fit, 'names'):
             self.name_to_value = {name: result_from_fit.rx2(name) for name in result_from_fit.names}
@@ -27,6 +29,14 @@ class ResultFromFit(object):
     def margin_coef_dict(self):
         raise NotImplementedError
 
+    @property
+    def nllh(self):
+        raise NotImplementedError
+
+    @property
+    def deviance(self):
+        raise NotImplementedError
+
 
 class ResultFromIsmev(ResultFromFit):
 
@@ -47,7 +57,8 @@ class ResultFromIsmev(ResultFromFit):
             i += 1
             # Add a potential linear temporal trend
             if gev_param_name in self.gev_param_name_to_dim:
-                temporal_coef_name = LinearCoef.coef_template_str(gev_param_name, AbstractCoordinates.COORDINATE_T).format(1)
+                temporal_coef_name = LinearCoef.coef_template_str(gev_param_name,
+                                                                  AbstractCoordinates.COORDINATE_T).format(1)
                 coef_dict[temporal_coef_name] = mle_values[i]
                 i += 1
         return coef_dict
@@ -58,7 +69,7 @@ class ResultFromIsmev(ResultFromFit):
 
     @property
     def nllh(self):
-        return self.res['nllh']
+        return self.name_to_value['nllh']
 
 
 class ResultFromSpatialExtreme(ResultFromFit):
@@ -68,6 +79,10 @@ class ResultFromSpatialExtreme(ResultFromFit):
     FITTED_VALUES_NAME = 'fitted.values'
     CONVERGENCE_NAME = 'convergence'
 
+    @property
+    def deviance(self):
+        return np.array(self.name_to_value['deviance'])[0]
+
     @property
     def convergence(self) -> str:
         convergence_value = self.name_to_value[self.CONVERGENCE_NAME]

test/test_extreme_estimator/test_estimator/test_margin_estimators.py

@@ -30,6 +30,7 @@ class TestSmoothMarginEstimator(unittest.TestCase):
                 # Fit estimator
                 estimator = LinearMarginEstimator(dataset=dataset, margin_model=margin_model)
                 estimator.fit()
+                print(estimator.result_from_fit.name_to_value.keys())
                 # Plot
                 if self.DISPLAY:
                     margin_model.margin_function_sample.visualize_function(show=True)