From 5349f19a7657f990d3d2a288737f3b010a89d9c3 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 13 May 2019 17:34:17 +0200
Subject: [PATCH] [EXTREME ESTIMATOR] add max iteration argument for any R fit.
also add some verbose for trend tests experiments.
---
.../safran/safran_variable.py | 6 +++++
.../main_study_visualizer.py | 27 ++++++++++---------
.../non_stationary_trends.py | 9 ++++++-
.../study_visualization/study_visualizer.py | 18 +++++++++----
.../max_stable_model/max_stable_fit.R | 8 +++---
extreme_estimator/extreme_models/utils.py | 12 ++++++++-
.../test_safe_run_r_estimator.py | 2 +-
7 files changed, 59 insertions(+), 23 deletions(-)
diff --git a/experiment/meteo_france_SCM_study/safran/safran_variable.py b/experiment/meteo_france_SCM_study/safran/safran_variable.py
index 1116d47c..864a6f58 100644
--- a/experiment/meteo_france_SCM_study/safran/safran_variable.py
+++ b/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
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
index 8a654900..72f2eb92 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
@@ -1,3 +1,5 @@
+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__':
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
index ea6eff9a..a77221a1 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/non_stationary_trends.py
@@ -1,3 +1,4 @@
+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))
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
index 7fb9b013..1a4b97ab 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
+++ b/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):
diff --git a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R b/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
index 744235c8..74195512 100644
--- a/extreme_estimator/extreme_models/max_stable_model/max_stable_fit.R
+++ b/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)
diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py
index fe854dda..6b4fe7ba 100644
--- a/extreme_estimator/extreme_models/utils.py
+++ b/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
diff --git a/test/test_extreme_estimator/test_extreme_models/test_safe_run_r_estimator.py b/test/test_extreme_estimator/test_extreme_models/test_safe_run_r_estimator.py
index d69cdbce..5a9e4b48 100644
--- a/test/test_extreme_estimator/test_extreme_models/test_safe_run_r_estimator.py
+++ b/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
--
GitLab