[SCM] add potential coordinate normalization.

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

@@ -7,6 +7,9 @@ from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, Exte
 from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
 from collections import OrderedDict
 
+from spatio_temporal_dataset.coordinates.transformed_coordinates.transformation.transformation_2D import \
+    BetweenZeroAndOne2DNormalization
+
 SCM_STUDIES = [SafranSnowfall, CrocusSwe, CrocusDepth]
 SCM_EXTENDED_STUDIES = [ExtendedSafranSnowfall, ExtendedCrocusSwe, ExtendedCrocusDepth]
 SCM_STUDY_TO_EXTENDED_STUDY = OrderedDict(zip(SCM_STUDIES, SCM_EXTENDED_STUDIES))
@@ -111,13 +114,14 @@ 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
     # [0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900, 4200, 4500, 4800] to test for others
-    altitudes = [300, 1200, 2100, 3000][:]
+    altitudes = [300, 1200, 2100, 3000][-1:]
     study_classes = [CrocusSwe, CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature]
     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)
+        study_visualizer = StudyVisualizer(study, save_to_file=save_to_file,
+                                           transformation_2D=BetweenZeroAndOne2DNormalization())
         study_visualizer.visualize_temporal_trend_relevance(complete_analysis=False)
 
 

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

-import math
 import os
 import os.path as op
 from collections import OrderedDict
+from typing import Union
 
+import math
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
@@ -16,8 +17,8 @@ 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, \
-    LinearNonStationaryLocationMarginModel, LinearStationaryMarginModel
+from extreme_estimator.extreme_models.margin_model.linear_margin_model import 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
@@ -29,6 +30,9 @@ from extreme_estimator.margin_fits.gpd.gpd_params import GpdParams
 from extreme_estimator.margin_fits.gpd.gpdmle_fit import GpdMleFit
 from spatio_temporal_dataset.coordinates.spatio_temporal_coordinates.abstract_spatio_temporal_coordinates import \
     AbstractSpatioTemporalCoordinates
+from spatio_temporal_dataset.coordinates.transformed_coordinates.transformation.transformation_2D import \
+    Transformation2D
+from spatio_temporal_dataset.coordinates.transformed_coordinates.transformed_coordinates import TransformedCoordinates
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 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
@@ -40,7 +44,8 @@ class StudyVisualizer(object):
 
     def __init__(self, study: AbstractStudy, show=True, save_to_file=False, only_one_graph=False, only_first_row=False,
                  vertical_kde_plot=False, year_for_kde_plot=None, plot_block_maxima_quantiles=False,
-                 temporal_non_stationarity=False):
+                 temporal_non_stationarity=False,
+                 transformation_2D=None):
         self.temporal_non_stationarity = temporal_non_stationarity
         self.only_first_row = only_first_row
         self.only_one_graph = only_one_graph
@@ -54,6 +59,7 @@ class StudyVisualizer(object):
         self._observations = None
 
         self.default_covariance_function = CovarianceFunction.powexp
+        self.transformation_2D = transformation_2D  # type: Union[None, Transformation2D]
 
         # KDE PLOT ARGUMENTS
         self.vertical_kde_plot = vertical_kde_plot
@@ -89,6 +95,9 @@ class StudyVisualizer(object):
     def coordinates(self):
         if self._coordinates is None:
             coordinates = self.study.massifs_coordinates
+            if self.transformation_2D is not None:
+                coordinates = TransformedCoordinates.from_coordinates(coordinates=coordinates,
+                                                                      transformation_function=self.transformation_2D)
             if self.temporal_non_stationarity:
                 # Build spatio temporal dataset from a temporal dataset
                 df_spatial = coordinates.df_spatial_coordinates()

extreme_estimator/extreme_models/utils.py

+import io
 import os.path as op
 import warnings
+from contextlib import redirect_stdout
 
 import numpy as np
 import random
@@ -20,8 +22,10 @@ r = ro.R()
 numpy2ri.activate()
 pandas2ri.activate()
 r.library('SpatialExtremes')
+default_filters = warnings.filters.copy()
+warnings.filterwarnings("ignore")
 r.library('ismev')
-
+warnings.filters = default_filters
 
 # Notice: R is not reloading all the time, the SpatialExtremes, so it's quite hard to debug or print in the code...
 # the best solution for debugging is to copy/paste the code module into a file that belongs to me, and then
@@ -43,6 +47,10 @@ def get_associated_r_file(python_filepath: str) -> str:
     return r_filepath
 
 
+class WarningWhileRunningR(Warning):
+    pass
+
+
 class WarningMaximumAbsoluteValueTooHigh(Warning):
     pass
 
@@ -64,22 +72,23 @@ def safe_run_r_estimator(function, data=None, use_start=False, threshold_max_abs
     # Then if it crashes, use start value
     run_successful = False
     res = None
-    while not run_successful:
-        current_parameter = parameters.copy()
-        if not use_start and 'start' in current_parameter:
-            current_parameter.pop('start')
-        try:
-            res = function(**current_parameter)  # type:
-            run_successful = True
-        except (RRuntimeError, RRuntimeWarning) as e:
-            if not use_start:
-                use_start = True
-                continue
-            elif isinstance(e, RRuntimeError):
-                raise Exception('Some R exception have been launched at RunTime: \n {}'.format(e.__repr__()))
-            if isinstance(e, RRuntimeWarning):
-                print(e.__repr__())
-                print('WARNING')
+    f = io.StringIO()
+    with redirect_stdout(f):
+        while not run_successful:
+            current_parameter = parameters.copy()
+            if not use_start and 'start' in current_parameter:
+                current_parameter.pop('start')
+            try:
+                res = function(**current_parameter)  # type:
+                run_successful = True
+            except (RRuntimeError, RRuntimeWarning) as e:
+                if not use_start:
+                    use_start = True
+                    continue
+                elif isinstance(e, RRuntimeError):
+                    raise Exception('Some R exception have been launched at RunTime: \n {}'.format(e.__repr__()))
+                if isinstance(e, RRuntimeWarning):
+                    warnings.warn(e.__repr__(), WarningWhileRunningR)
     return res
 
 

spatio_temporal_dataset/coordinates/abstract_coordinates.py

 import os.path as op
-from typing import List, Tuple
+from typing import List, Tuple, Union
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -50,7 +50,7 @@ class AbstractCoordinates(object):
         # Slicing attributes
         self.s_split_spatial = s_split_spatial  # type: pd.Series
         self.s_split_temporal = s_split_temporal  # type: pd.Series
-        self.slicer = None  # type: AbstractSlicer
+        self.slicer = None  # type: Union[None, AbstractSlicer]
 
         # Load the slicer
         if slicer_class is TemporalSlicer:

spatio_temporal_dataset/coordinates/transformed_coordinates/transformation/transformation_2D.py

@@ -40,3 +40,9 @@ class BetweenZeroAndOne2DNormalization(Uniform2DNormalization):
         s_coord_shifted = s_coord - self.min_coord
         s_coord_scaled = s_coord_shifted / (self.max_coord - self.min_coord)
         return s_coord_scaled
+
+
+class BetweenMinusOneAndOne2DNormalization(BetweenZeroAndOne2DNormalization):
+
+    def uniform_normalization(self, s_coord: pd.Series) -> pd.Series:
+        pass

test/test_spatio_temporal_dataset/test_dataset.py

@@ -52,7 +52,7 @@ class TestSpatioTemporalDataset(unittest.TestCase):
                                                    margin_model=smooth_margin_model,
                                                    coordinates=self.coordinates)
 
-    def test_spatio_temporal_array(self):
+    def test_spatio_temporal_array_wrt_time(self):
         # The test could have been on a given station. But we decided to do it for a given time step.
         self.load_dataset(nb_obs=1)
 
@@ -71,7 +71,26 @@ class TestSpatioTemporalDataset(unittest.TestCase):
         self.assertTrue(equality, msg='v1={} is different from v2={}'.format(observation_at_time_0_v1,
                                                                              observation_at_time_0_v2))
 
-    def test_spatio_temporal_case_to_resolve(self):
+    def test_spatio_temporal_array_wrt_space(self):
+        # The test could have been on a given station. But we decided to do it for a given time step.
+        self.load_dataset(nb_obs=1)
+
+        # Load observation for time 0
+        ind_station_0 = self.dataset.coordinates.ind_of_df_all_coordinates(
+            coordinate_name=AbstractCoordinates.COORDINATE_X,
+            value=-1)
+        observation_at_station_0_v1 = self.dataset.observations.df_maxima_gev.loc[ind_station_0].values.flatten()
+
+        # Load observation correspond to time 0
+        maxima_gev = self.dataset.maxima_gev_for_spatial_extremes_package()
+        maxima_gev = np.transpose(maxima_gev)
+        self.assertEqual(maxima_gev.shape, (3, 2))
+        observation_at_time_0_v2 = maxima_gev[0, :]
+        equality = np.equal(observation_at_station_0_v1, observation_at_time_0_v2).all()
+        self.assertTrue(equality, msg='v1={} is different from v2={}'.format(observation_at_station_0_v1,
+                                                                             observation_at_time_0_v2))
+
+    def test_spatio_temporal_array_with_multiple_observations(self):
         # In this case, we must check that the observations are the same
         self.load_dataset(nb_obs=2)