From b3e27c1b6473ef6f64a54fb6ba4213bf0434baf5 Mon Sep 17 00:00:00 2001 From: Le Roux Erwan <erwan.le-roux@irstea.fr> Date: Fri, 10 May 2019 13:40:38 +0200 Subject: [PATCH] [SCM] add potential coordinate normalization. --- .../main_study_visualizer.py | 12 ++++-- .../study_visualization/study_visualizer.py | 17 ++++++-- extreme_estimator/extreme_models/utils.py | 43 +++++++++++-------- .../coordinates/abstract_coordinates.py | 4 +- .../transformation/transformation_2D.py | 6 +++ .../test_dataset.py | 23 +++++++++- 6 files changed, 76 insertions(+), 29 deletions(-) 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 53144ef8..818c4be0 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 @@ -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) 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 799ad66a..457bf951 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 @@ -1,8 +1,9 @@ -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() diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py index 8ca53f60..fe854dda 100644 --- a/extreme_estimator/extreme_models/utils.py +++ b/extreme_estimator/extreme_models/utils.py @@ -1,5 +1,7 @@ +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 diff --git a/spatio_temporal_dataset/coordinates/abstract_coordinates.py b/spatio_temporal_dataset/coordinates/abstract_coordinates.py index 7b6ae41b..91c85d25 100644 --- a/spatio_temporal_dataset/coordinates/abstract_coordinates.py +++ b/spatio_temporal_dataset/coordinates/abstract_coordinates.py @@ -1,5 +1,5 @@ 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: diff --git a/spatio_temporal_dataset/coordinates/transformed_coordinates/transformation/transformation_2D.py b/spatio_temporal_dataset/coordinates/transformed_coordinates/transformation/transformation_2D.py index 4344c669..d92ec4ae 100644 --- a/spatio_temporal_dataset/coordinates/transformed_coordinates/transformation/transformation_2D.py +++ b/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 diff --git a/test/test_spatio_temporal_dataset/test_dataset.py b/test/test_spatio_temporal_dataset/test_dataset.py index 9aea31b2..481314b9 100644 --- a/test/test_spatio_temporal_dataset/test_dataset.py +++ b/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) -- GitLab