[EXTREME ESTIMATOR] add max iteration argument for any R fit. also add some...

[EXTREME ESTIMATOR] add max iteration argument for any R fit. also add some verbose for trend tests experiments.

experiment/meteo_france_SCM_study/safran/safran_variable.py

@@ -57,6 +57,9 @@ class SafranSnowfallVariable(AbstractVariable):
 
 class SafranRainfallVariable(SafranSnowfallVariable):
 
+    NAME = 'Rainfall'
+    UNIT = 'kg per m2 or mm'
+
     def __init__(self, dataset, altitude, nb_consecutive_days_of_snowfall=1, keyword='Rainf'):
         super().__init__(dataset, altitude, nb_consecutive_days_of_snowfall, keyword)
 
@@ -75,6 +78,9 @@ class SafranTotalPrecipVariable(AbstractVariable):
 
 class SafranTemperatureVariable(AbstractVariable):
 
+    NAME = 'Temperature'
+    UNIT = 'Celsius Degrees'
+
     def __init__(self, dataset, altitude, keyword='Tair'):
         super().__init__(dataset, altitude)
         # Temperature are in K, I transform them as celsius

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

+from typing import Generator, List
+
 from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
 from experiment.meteo_france_SCM_study.crocus.crocus import CrocusDepth, CrocusSwe, ExtendedCrocusDepth, \
     ExtendedCrocusSwe, CrocusDaysWithSnowOnGround
@@ -15,7 +17,7 @@ SCM_EXTENDED_STUDIES = [ExtendedSafranSnowfall, ExtendedCrocusSwe, ExtendedCrocu
 SCM_STUDY_TO_EXTENDED_STUDY = OrderedDict(zip(SCM_STUDIES, SCM_EXTENDED_STUDIES))
 
 
-def study_iterator_global(study_classes, only_first_one=False, both_altitude=False, verbose=True, altitudes=None):
+def study_iterator_global(study_classes, only_first_one=False, both_altitude=False, verbose=True, altitudes=None) -> List[AbstractStudy]:
     for study_class in study_classes:
         for study in study_iterator(study_class, only_first_one, both_altitude, verbose, altitudes):
             yield study
@@ -23,24 +25,24 @@ def study_iterator_global(study_classes, only_first_one=False, both_altitude=Fal
             break
 
 
-def study_iterator(study_class, only_first_one=False, both_altitude=False, verbose=True, altitudes=None):
+def study_iterator(study_class, only_first_one=False, both_altitude=False, verbose=True, altitudes=None) -> List[AbstractStudy]:
     all_studies = []
     is_safran_study = study_class in [SafranSnowfall, ExtendedSafranSnowfall]
     nb_days = [1] if is_safran_study else [1]
     if verbose:
-        print('Loading studies....')
+        print('\n\n\n\n\nLoading studies....', end='')
     for nb_day in nb_days:
         altis = [1800] if altitudes is None else altitudes
 
         for alti in altis:
 
             if verbose:
-                print('alti: {}, nb_day: {}'.format(alti, nb_day))
+                print('alti: {}, nb_day: {}'.format(alti, nb_day), end='')
             study = study_class(altitude=alti, nb_consecutive_days=nb_day) if is_safran_study \
                 else study_class(altitude=alti)
             massifs = study.altitude_to_massif_names[alti]
             if verbose:
-                print('{} massifs: {}'.format(len(massifs), massifs))
+                print('{} massifs: {} \n'.format(len(massifs), massifs))
             yield study
             if only_first_one and not both_altitude:
                 break
@@ -114,16 +116,17 @@ def complete_analysis(only_first_one=False):
 
 
 def trend_analysis():
-    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][-1:]
-    normalization_class = [BetweenZeroAndOneNormalization, BetweenMinusOneAndOneNormalization][1]
-    study_classes = [CrocusSwe, CrocusDepth, SafranSnowfall, SafranRainfall, SafranTemperature][:1]
+    save_to_file = True
+    only_first_one = False
+    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][:]
+    altitudes = full_altitude_with_at_least_2_stations
+    normalization_class = [None, BetweenMinusOneAndOneNormalization, BetweenZeroAndOneNormalization][-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,
                                            transformation_class=normalization_class)
-        study_visualizer.visualize_temporal_trend_relevance(complete_analysis=False)
+        study_visualizer.visualize_temporal_trend_relevance(complete_analysis=True)
 
 
 if __name__ == '__main__':

experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py

+import time
 from typing import Union
 
 from experiment.meteo_france_SCM_study.visualization.utils import align_yaxis_on_zero
@@ -11,6 +12,7 @@ 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.utils import OptimizationConstants
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from utils import get_display_name_from_object_type
 
@@ -38,7 +40,11 @@ class AbstractNonStationaryTrendTest(object):
                 estimator_name = get_display_name_from_object_type(estimator)
                 margin_model_name = get_display_name_from_object_type(margin_model)
                 print('Fitting {} with margin: {} for starting_point={}'.format(estimator_name, margin_model_name, starting_point))
+            start = time.time()
             estimator.fit()
+            duration = time.time() - start
+            if self.verbose:
+                print('Fit took {}s and was {}'.format(round(duration, 1), estimator.result_from_fit.convergence))
             self._margin_model_class_and_starting_point_to_estimator[(margin_model_class, starting_point)] = estimator
         return self._margin_model_class_and_starting_point_to_estimator[(margin_model_class, starting_point)]
 
@@ -83,7 +89,7 @@ class AbstractNonStationaryTrendTest(object):
         color_mu1 = 'c'
 
         if self.verbose:
-            print(mu1_trends)
+            print('mu1 trends:', mu1_trends, '\n')
         ax2.plot(years, mu1_trends, color_mu1 + 'o-')
         ax2.set_ylabel('mu1 parameter', color=color_mu1)
 
@@ -98,6 +104,7 @@ class AbstractNonStationaryTrendTest(object):
         # Define the year_min and year_max for the starting point
         if complete_analysis:
             year_min, year_max, step = 1960, 1990, 1
+            OptimizationConstants.USE_MAXIT = True
         else:
             year_min, year_max, step = 1960, 1990, 5
         years = list(range(year_min, year_max + 1, step))

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

@@ -47,7 +47,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,
-                 transformation_class=None):
+                 transformation_class=None,
+                 normalization_under_one_observations=True):
         self.temporal_non_stationarity = temporal_non_stationarity
         self.only_first_row = only_first_row
         self.only_one_graph = only_one_graph
@@ -55,6 +56,7 @@ class StudyVisualizer(object):
         self.study = study
         self.plot_name = None
 
+        self.normalization_under_one_observations = normalization_under_one_observations
         # Load some attributes
         self._dataset = None
         self._coordinates = None
@@ -130,6 +132,8 @@ class StudyVisualizer(object):
             self._observations = self.study.observations_annual_maxima
             if self.temporal_non_stationarity:
                 self._observations.convert_to_spatio_temporal_index(self.coordinates)
+                if self.normalization_under_one_observations:
+                    self._observations.normalize()
         return self._observations
 
     # Graph for each massif / or groups of massifs
@@ -160,13 +164,13 @@ 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):
+    def visualize_temporal_trend_relevance(self, complete_analysis, verbose=True):
         self.temporal_non_stationarity = True
-        trend_tests = [ConditionalIndedendenceLocationTrendTest(self.dataset)]
+        trend_tests = [ConditionalIndedendenceLocationTrendTest(self.dataset, verbose=verbose)]
 
         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(self.dataset, max_stable_model))
+            trend_tests.append(MaxStableLocationTrendTest(self.dataset, max_stable_model, verbose=verbose))
 
         nb_trend_tests = len(trend_tests)
         fig, axes = plt.subplots(1, nb_trend_tests, figsize=self.figsize)
@@ -176,7 +180,11 @@ class StudyVisualizer(object):
         for ax, trend_test in zip(axes, trend_tests):
             trend_test.visualize(ax, complete_analysis=complete_analysis)
 
-        self.plot_name = 'trend tests'
+        plot_name = 'trend tests'
+        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()
 
     def visualize_experimental_law(self, ax, massif_id):

extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R

@@ -70,7 +70,6 @@ rmaxstab2D <- function (n.obs){
 #     print(res['fitted.values'])
 # }
 
-fitspatgev()
 
 # rmaxstab with 1D data
 rmaxstab1D <- function (n.obs){
@@ -90,7 +89,10 @@ rmaxstab1D <- function (n.obs){
 
     # GAUSS
     namedlist = list(cov=1.0, locCoeff1=1.0, scaleCoeff1=1.0, shapeCoeff1=1.0)
-    res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist, fit.marge=TRUE, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form, iso=TRUE)
+    res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist, fit.marge=TRUE, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form, iso=TRUE, control=list(maxit=1000))
+
+    # ‘eval.max’
+    # ‘iter.max’
 
     # BROWN
     # namedlist = list(range = 3, smooth = 0.5, locCoeff1=1.0, scaleCoeff1=1.0, shapeCoeff1=1.0)
@@ -107,7 +109,7 @@ if (call_main) {
     n.obs = 500
     # rmaxstab2D(n.obs)
     # rmaxstab3D(n.obs)
-    # rmaxstab1D(n.obs)
+    rmaxstab1D(n.obs)
 
     # namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2)
     # res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist)

extreme_estimator/extreme_models/utils.py

@@ -27,6 +27,7 @@ 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
 # I can put print & stop in the code, and I can understand where are the problems
@@ -54,9 +55,18 @@ class WarningWhileRunningR(Warning):
 class WarningMaximumAbsoluteValueTooHigh(Warning):
     pass
 
+class OptimizationConstants(object):
+
+    USE_MAXIT = False
 
-def safe_run_r_estimator(function, data=None, use_start=False, threshold_max_abs_value=100,
+
+def safe_run_r_estimator(function, data=None, use_start=False, threshold_max_abs_value=100, maxit=1000000,
                          **parameters) -> robjects.ListVector:
+    if OptimizationConstants.USE_MAXIT:
+        # Add optimization parameters
+        optim_dict = {'maxit': maxit}
+        parameters['control'] = r.list(**optim_dict)
+
     # Some checks for Spatial Extremes
     if data is not None:
         # Raise warning if the maximum absolute value is above a threshold

test/test_extreme_estimator/test_extreme_models/test_safe_run_r_estimator.py

@@ -4,7 +4,7 @@ import unittest
 from extreme_estimator.extreme_models.utils import safe_run_r_estimator, WarningMaximumAbsoluteValueTooHigh
 
 
-def function(data):
+def function(data=None, control=None):
     pass