[SCM][NON STATIONARY TRENDS] add independent trends visualization for each massif

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

@@ -119,8 +119,8 @@ def complete_analysis(only_first_one=False):
 
 
 def trend_analysis():
-    save_to_file = True
-    only_first_one = False
+    save_to_file = False
+    only_first_one = True
     short_altitudes = [300, 1200, 2100, 3000][:1]
     full_altitude_with_at_least_2_stations = [0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200][:]
     durand_altitude = [1800]
@@ -129,14 +129,23 @@ def trend_analysis():
     study_classes = [CrocusSwe, CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][2:3]
     for study in study_iterator_global(study_classes, only_first_one=only_first_one, altitudes=altitudes):
         study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
-                                           transformation_class=normalization_class)
-        study_visualizer.visualize_temporal_trend_relevance(complete_analysis=False, multiprocessing=True)
+                                           transformation_class=normalization_class,
+                                           temporal_non_stationarity=True,
+                                           verbose=True,
+                                           multiprocessing=True,
+                                           complete_non_stationary_trend_analysis=False)
+        # study_visualizer.only_one_graph = True
+        study_visualizer.visualize_all_independent_temporal_trend()
+        # study_visualizer.visualize_temporal_trend_relevance()
 
 
 def main_run():
+    # normal_visualization(temporal_non_stationarity=True)
+    trend_analysis()
+
+
+
     # annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2100)
-    normal_visualization(temporal_non_stationarity=True)
-    # trend_analysis()
 
 
     # max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1)

experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py

@@ -15,6 +15,8 @@ from extreme_estimator.extreme_models.margin_model.linear_margin_model import \
     LinearAllParametersTwoFirstCoordinatesMarginModel, LinearAllTwoStatialCoordinatesLocationLinearMarginModel, \
     LinearStationaryMarginModel, LinearNonStationaryLocationMarginModel
 from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction
+from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import StationaryStationModel, \
+    NonStationaryStationModel
 from extreme_estimator.extreme_models.utils import OptimizationConstants
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from utils import get_display_name_from_object_type
@@ -45,7 +47,6 @@ class AbstractNonStationaryTrendTest(object):
             self._starting_point_to_estimator[starting_point] = estimator
         return self._starting_point_to_estimator[starting_point]
 
-
     def load_estimator(self, starting_point) -> Union[
         AbstractFullEstimator, AbstractMarginEstimator]:
         margin_model_class = self.stationary_margin_model_class if starting_point is None else self.non_stationary_margin_model_class
@@ -79,17 +80,25 @@ class AbstractNonStationaryTrendTest(object):
         estimator = self.get_estimator(starting_point)
         margin_function = estimator.margin_function_fitted  # type: LinearMarginFunction
         assert isinstance(margin_function, LinearMarginFunction)
-        mu_coefs = [margin_function.mu_intercept, margin_function.mu_longitude_trend, margin_function.mu_latitude_trend,
-                    margin_function.mu1_temporal_trend]
+        mu_coefs = [margin_function.mu_intercept, margin_function.mu1_temporal_trend]
+        if self.has_spatial_coordinates:
+            mu_coefs += [margin_function.mu_longitude_trend, margin_function.mu_latitude_trend]
         return dict(zip(self.mu_coef_names, mu_coefs))
 
     @property
     def mu_coef_names(self):
-        return ['mu_intercept', 'mu_longitude', 'mu_latitude', 'mu_temporal']
+        mu_coef_names = ['mu_intercept', 'mu_temporal']
+        if self.has_spatial_coordinates:
+            mu_coef_names += ['mu_longitude', 'mu_latitude']
+        return mu_coef_names
+
+    @property
+    def has_spatial_coordinates(self):
+        return self.dataset.coordinates.has_spatial_coordinates
 
     @property
     def mu_coef_colors(self):
-        return ['b', 'g', 'y', 'c']
+        return ['b', 'c', 'g', 'y', ]
 
     def visualize(self, ax, complete_analysis=True):
         years = self.years(complete_analysis)
@@ -141,17 +150,17 @@ class AbstractNonStationaryTrendTest(object):
         df_mus = pd.DataFrame(mus)
         min_mus, max_mus = df_mus.min().min(), df_mus.max().max()
         min_global, max_global = min(min_deviance_data, min_mus), max(max_deviance_data, max_mus)
-        ax2.set_ylim(min_global, max_global)
+        # ax2.set_ylim(min_global, max_global)
+        # if min_mus < 0.0 < max_mus:
+        #     align_yaxis_on_zero(ax2, ax)
 
+        ax.set_title(self.display_name)
         ax.set_xlabel('starting year for the linear trend of {}'.format(self.mu_coef_names[-1]))
-        if min_mus < 0.0 < max_mus:
-            align_yaxis_on_zero(ax2, ax)
-
-        title = self.display_name
-        ax.set_title(title)
         ax.grid()
-        ax.legend(loc=6)
-        ax2.legend(loc=7)
+
+        prop = {'size': 5} if not self.has_spatial_coordinates else None
+        ax.legend(loc=6, prop=prop)
+        ax2.legend(loc=7, prop=prop)
 
     def years(self, complete_analysis=True):
         # Define the year_min and year_max for the starting point
@@ -170,8 +179,16 @@ class AbstractNonStationaryTrendTest(object):
 
 class IndependenceLocationTrendTest(AbstractNonStationaryTrendTest):
 
-    def __init__(self, dataset, coordinate_idx):
-        pass
+    def __init__(self, station_name, *args, **kwargs):
+        super().__init__(*args, **kwargs,
+                         estimator_class=LinearMarginEstimator,
+                         stationary_margin_model_class=StationaryStationModel,
+                         non_stationary_margin_model_class=NonStationaryStationModel)
+        self.station_name = station_name
+
+    @property
+    def display_name(self):
+        return self.station_name
 
 
 class ConditionalIndedendenceLocationTrendTest(AbstractNonStationaryTrendTest):

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

 import os
 import os.path as op
 from collections import OrderedDict
+from copy import deepcopy
 from typing import Union
 
 import math
@@ -11,7 +12,7 @@ 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
+    ConditionalIndedendenceLocationTrendTest, MaxStableLocationTrendTest, IndependenceLocationTrendTest
 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 \
@@ -36,6 +37,9 @@ from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal
     ConsecutiveTemporalCoordinates
 from spatio_temporal_dataset.coordinates.transformed_coordinates.transformed_coordinates import TransformedCoordinates
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+    AbstractSpatioTemporalObservations
+from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
 from test.test_utils import load_test_max_stable_models
 from utils import get_display_name_from_object_type, VERSION_TIME, float_to_str_with_only_some_significant_digits
 
@@ -48,7 +52,11 @@ class StudyVisualizer(object):
                  vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
                  temporal_non_stationarity=False,
                  transformation_class=None,
+                 verbose=False,
+                 multiprocessing=False,
+                 complete_non_stationary_trend_analysis=False,
                  normalization_under_one_observations=True):
+
         self.massif_id_to_smooth_maxima = {}
         self.temporal_non_stationarity = temporal_non_stationarity
         self.only_first_row = only_first_row
@@ -58,6 +66,10 @@ class StudyVisualizer(object):
         self.plot_name = None
 
         self.normalization_under_one_observations = normalization_under_one_observations
+        self.multiprocessing = multiprocessing
+        self.verbose = verbose
+        self.complete_non_stationary_trend_analysis = complete_non_stationary_trend_analysis
+
         # Load some attributes
         self._dataset = None
         self._coordinates = None
@@ -135,10 +147,21 @@ class StudyVisualizer(object):
                 self._observations.convert_to_spatio_temporal_index(self.coordinates)
                 if self.normalization_under_one_observations:
                     self._observations.normalize()
-                self._observations.print_summary()
+                if self.verbose:
+                    self._observations.print_summary()
         return self._observations
 
-    # Graph for each massif / or groups of massifs
+    def observation_massif_id(self, massif_id):
+        annual_maxima = [data[massif_id] for data in self.study.year_to_annual_maxima.values()]
+        df_annual_maxima = pd.DataFrame(annual_maxima, index=self.temporal_coordinates.index)
+        observation_massif_id = AnnualMaxima(df_maxima_gev=df_annual_maxima)
+        if self.normalization_under_one_observations:
+            observation_massif_id.normalize()
+        if self.verbose:
+            observation_massif_id.print_summary()
+        return observation_massif_id
+
+    # PLOT FOR SEVERAL MASSIFS
 
     def visualize_massif_graphs(self, visualize_function, specified_massif_names=None):
         if self.only_one_graph:
@@ -157,29 +180,39 @@ class StudyVisualizer(object):
                 if specified_massif_names is None:
                     massif_ids = list(range(len(self.study.study_massif_names)))
                 else:
-                    massif_ids = [self.study.study_massif_names.index(massif_name) for massif_name in specified_massif_names]
+                    massif_ids = [self.study.study_massif_names.index(massif_name) for massif_name in
+                                  specified_massif_names]
                 for j, massif_id in enumerate(massif_ids):
                     row_id, column_id = j // nb_columns, j % nb_columns
                     ax = axes[row_id, column_id]
                     visualize_function(ax, massif_id)
 
-    def visualize_all_experimental_law( self):
-        self.visualize_massif_graphs(self.visualize_experimental_law)
-        self.plot_name = ' Empirical distribution \n'
-        self.plot_name += 'with data from the 23 mountain chains of the French Alps ' if self.year_for_kde_plot is None else \
-            'for the year {}'.format(self.year_for_kde_plot)
+    # TEMPORAL TREND
+
+    def visualize_all_independent_temporal_trend(self):
+        self.visualize_massif_graphs(self.visualize_independent_temporal_trend)
+        self.plot_name = ' Independent temporal trend \n'
         self.show_or_save_to_file()
 
-    def visualize_temporal_trend_relevance(self, complete_analysis, verbose=True, multiprocessing=False):
-        self.temporal_non_stationarity = True
-        trend_tests = [ConditionalIndedendenceLocationTrendTest(dataset=self.dataset, verbose=verbose,
-                                                                multiprocessing=multiprocessing)]
+    def visualize_independent_temporal_trend(self, ax, massif_id):
+        assert self.temporal_non_stationarity
+        # Create a dataset with temporal coordinates from the massif id
+        dataset_massif_id = AbstractDataset(self.observation_massif_id(massif_id), self.temporal_coordinates)
+        trend_test = IndependenceLocationTrendTest(station_name=self.study.study_massif_names[massif_id],
+                                                   dataset=dataset_massif_id, verbose=self.verbose,
+                                                   multiprocessing=self.multiprocessing)
+        trend_test.visualize(ax, self.complete_non_stationary_trend_analysis)
+
+    def visualize_temporal_trend_relevance(self):
+        assert self.temporal_non_stationarity
+        trend_tests = [ConditionalIndedendenceLocationTrendTest(dataset=self.dataset, verbose=self.verbose,
+                                                                multiprocessing=self.multiprocessing)]
 
         max_stable_models = load_test_max_stable_models(default_covariance_function=self.default_covariance_function)
         for max_stable_model in [max_stable_models[1], max_stable_models[-2]]:
             trend_tests.append(MaxStableLocationTrendTest(dataset=self.dataset,
-                                                          max_stable_model=max_stable_model, verbose=verbose,
-                                                          multiprocessing=multiprocessing))
+                                                          max_stable_model=max_stable_model, verbose=self.verbose,
+                                                          multiprocessing=self.multiprocessing))
 
         nb_trend_tests = len(trend_tests)
         fig, axes = plt.subplots(1, nb_trend_tests, figsize=self.figsize)
@@ -187,15 +220,25 @@ class StudyVisualizer(object):
             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)
+            trend_test.visualize(ax, complete_analysis=self.complete_non_stationary_trend_analysis)
 
         plot_name = 'trend tests'
-        plot_name += ' with {} applied spatially & temporally'.format(get_display_name_from_object_type(self.transformation_class))
+        plot_name += ' with {} applied spatially & temporally'.format(
+            get_display_name_from_object_type(self.transformation_class))
         if self.normalization_under_one_observations:
             plot_name += '(and maxima <= 1)'
         self.plot_name = plot_name
         self.show_or_save_to_file()
 
+    # EXPERIMENTAL LAW
+
+    def visualize_all_experimental_law(self):
+        self.visualize_massif_graphs(self.visualize_experimental_law)
+        self.plot_name = ' Empirical distribution \n'
+        self.plot_name += 'with data from the 23 mountain chains of the French Alps ' if self.year_for_kde_plot is None else \
+            'for the year {}'.format(self.year_for_kde_plot)
+        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)
@@ -281,6 +324,8 @@ class StudyVisualizer(object):
         all_massif_data = np.sort(all_massif_data)
         return all_massif_data
 
+    #
+
     def visualize_all_mean_and_max_graphs(self):
         # Compute the order of massif names
         massif_name_to_score = {}
@@ -289,7 +334,8 @@ class StudyVisualizer(object):
             massif_name_to_score[massif_name] = score
         ordered_massif_names = sorted(self.study.study_massif_names[:], key=lambda s: massif_name_to_score[s])
         self.visualize_massif_graphs(self.visualize_mean_and_max_graph, specified_massif_names=ordered_massif_names)
-        self.plot_name = ' smoothing values temporally with sliding window of size {}'.format(self.window_size_for_smoothing)
+        self.plot_name = ' smoothing values temporally with sliding window of size {}'.format(
+            self.window_size_for_smoothing)
         self.show_or_save_to_file()
 
     def visualize_mean_and_max_graph(self, ax, massif_id):
@@ -313,7 +359,6 @@ class StudyVisualizer(object):
         ax.set_title(self.study.study_massif_names[massif_id])
         ax.xaxis.set_ticks(x[2::10])
 
-
     def smooth_maxima_x_y(self, massif_id):
         if massif_id not in self.massif_id_to_smooth_maxima:
             tuples_x_y = [(year, annual_maxima[massif_id]) for year, annual_maxima in

extreme_estimator/extreme_models/result_from_fit.py

@@ -8,8 +8,10 @@ from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
-class ResultFromFit(object):
+def convertFloatVector_to_float(f):
+    return np.array(f)[0]
 
+class ResultFromFit(object):
 
     def __init__(self, result_from_fit: robjects.ListVector) -> None:
         if hasattr(result_from_fit, 'names'):
@@ -73,7 +75,16 @@ class ResultFromIsmev(ResultFromFit):
 
     @property
     def nllh(self):
-        return self.name_to_value['nllh']
+        return convertFloatVector_to_float(self.name_to_value['nllh'])
+
+    @property
+    def deviance(self):
+        return - 2 * self.nllh
+
+    @property
+    def convergence(self) -> str:
+        return convertFloatVector_to_float(self.name_to_value['conv']) == 0
+
 
 
 class ResultFromSpatialExtreme(ResultFromFit):

spatio_temporal_dataset/spatio_temporal_observations/abstract_spatio_temporal_observations.py

@@ -3,6 +3,8 @@ import pandas as pd
 import numpy as np
 
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
+    AbstractTemporalCoordinates
 from spatio_temporal_dataset.slicer.abstract_slicer import df_sliced, AbstractSlicer
 from spatio_temporal_dataset.slicer.split import Split
 
@@ -91,6 +93,24 @@ class AbstractSpatioTemporalObservations(object):
         if df is not None:
             return pd.DataFrame(np.concatenate([df[c].values for c in df.columns]), index=index)
 
+    def convert_to_temporal_index(self, temporal_coordinates: AbstractTemporalCoordinates, spatial_idx: int):
+        assert len(self.index) > len(temporal_coordinates) and len(self.index) % len(temporal_coordinates) == 0
+        spatial_len = len(self.index) // len(temporal_coordinates)
+        assert 0 <= spatial_idx < spatial_len
+        # Build ind to select the observations of interest
+        ind = np.zeros(spatial_len, dtype=bool)
+        ind[spatial_idx] = True
+        ind = np.concatenate([ind for _ in range(len(temporal_coordinates))])
+        self.df_maxima_frech = self.loc_df(self.df_maxima_frech, ind, temporal_coordinates.index)
+        self.df_maxima_gev = self.loc_df(self.df_maxima_gev, ind, temporal_coordinates.index)
+
+    @staticmethod
+    def loc_df(df, ind, new_index):
+        if df is not None:
+            df = df.loc[ind]
+            df.index = new_index
+            return df
+
     def maxima_gev(self, split: Split = Split.all, slicer: AbstractSlicer = None) -> np.ndarray:
         return df_sliced(self.df_maxima_gev, split, slicer).values
 
@@ -118,5 +138,3 @@ class AbstractSpatioTemporalObservations(object):
     @_df_maxima.setter
     def _df_maxima(self, value):
         self.__df_maxima = value
-
-