diff --git a/experiment/regression_margin/regression_margin.py b/experiment/regression_margin/regression_margin.py index bb96a5dfbae56a181a5409a10092ef6bedd1233e..a46fd2022bbbe8687164a15e7f2e19d406f98314 100644 --- a/experiment/regression_margin/regression_margin.py +++ b/experiment/regression_margin/regression_margin.py @@ -1,8 +1,11 @@ +import random + import numpy as np from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStepWithSmoothMargin -from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeAxis0MarginModel, \ - LinearAllParametersAllAxisMarginModel +from extreme_estimator.extreme_models.margin_model.margin_function.linear_margin_function import LinearMarginFunction +from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeDim1MarginModel, \ + LinearAllParametersAllDimsMarginModel from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith from extreme_estimator.gev_params import GevParams from spatio_temporal_dataset.coordinates.unidimensional_coordinates.coordinates_1D import LinSpaceCoordinates @@ -12,6 +15,7 @@ from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedData nb_points = 5 nb_obs = 10 +nb_estimator = 2 show = False coordinates = LinSpaceCoordinates.from_nb_points(nb_points=nb_points) @@ -19,25 +23,40 @@ coordinates = LinSpaceCoordinates.from_nb_points(nb_points=nb_points) ########## GENERATING THE DATA ##################### # MarginModel Linear with respect to the shape (from 0.01 to 0.02) -margin_model = LinearShapeAxis0MarginModel(coordinates=coordinates, params_sample={GevParams.GEV_SHAPE: 0.02}) +params_sample = { + (GevParams.GEV_SHAPE, 0): 0.2, + (GevParams.GEV_SHAPE, 1): 0.05, +} +margin_model = LinearShapeDim1MarginModel(coordinates=coordinates, params_sample=params_sample) max_stable_model = Smith() -dataset = FullSimulatedDataset.from_double_sampling(nb_obs=nb_obs, margin_model=margin_model, - coordinates=coordinates, - max_stable_model=max_stable_model) -if show: - # Visualize the sampling margin - dataset.margin_model.margin_function_sample.visualize_all() - # Plot a realization from the maxima gev (i.e the maxima obtained just by simulating the marginal law) - for maxima_gev in np.transpose(dataset.maxima_gev): - plt.plot(coordinates.coordinates_values, maxima_gev) - plt.show() + +# if show: +# # Plot a realization from the maxima gev (i.e the maxima obtained just by simulating the marginal law) +# for maxima_gev in np.transpose(dataset.maxima_gev): +# plt.plot(coordinates.coordinates_values, maxima_gev) +# plt.show() ######### FITTING A MODEL ################# -margin_model = LinearAllParametersAllAxisMarginModel(coordinates) -full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model, max_stable_model) -full_estimator.fit() -full_estimator.margin_function_fitted.visualize_all() + +axes = None +for _ in range(nb_estimator): + # Data part + dataset = FullSimulatedDataset.from_double_sampling(nb_obs=nb_obs, margin_model=margin_model, + coordinates=coordinates, + max_stable_model=max_stable_model) + margin_function_sample = dataset.margin_model.margin_function_sample # type: LinearMarginFunction + margin_function_sample.visualize(show=False, axes=axes) + axes = margin_function_sample.visualization_axes + + # Estimation part + margin_model_for_estimator = LinearAllParametersAllDimsMarginModel(coordinates) + full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model_for_estimator, max_stable_model) + full_estimator.fit() + full_estimator.margin_function_fitted.visualize(axes=axes, show=False) +plt.show() + +# Display all the margin on the same graph for comparison # Plot the margin functions # margin_model.margin_function_sample.visualize_2D() diff --git a/experiment/regression_margin/visualization_margin_model.py b/experiment/regression_margin/visualization_margin_model.py index 1f405c9d357c91940cb79f2f5a6c7ebe722425f6..b5f548b32dd9344d0b26cb0453de460a2752099d 100644 --- a/experiment/regression_margin/visualization_margin_model.py +++ b/experiment/regression_margin/visualization_margin_model.py @@ -1,8 +1,8 @@ import unittest from extreme_estimator.gev_params import GevParams -from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeAxis0MarginModel, \ - LinearAllParametersAllAxisMarginModel +from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeDim1MarginModel, \ + LinearAllParametersAllDimsMarginModel from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import CircleCoordinates from spatio_temporal_dataset.coordinates.unidimensional_coordinates.coordinates_1D import LinSpaceCoordinates @@ -10,14 +10,14 @@ from spatio_temporal_dataset.coordinates.unidimensional_coordinates.coordinates_ class VisualizationMarginModel(unittest.TestCase): DISPLAY = True nb_points = 50 - margin_model = [LinearShapeAxis0MarginModel, LinearAllParametersAllAxisMarginModel][-1] + margin_model = [LinearShapeDim1MarginModel, LinearAllParametersAllDimsMarginModel][-1] @classmethod def example_visualization_2D(cls): spatial_coordinates = CircleCoordinates.from_nb_points(nb_points=cls.nb_points) margin_model = cls.margin_model(coordinates=spatial_coordinates) if cls.DISPLAY: - margin_model.margin_function_sample.visualize_all() + margin_model.margin_function_sample.visualize() @classmethod def example_visualization_1D(cls): @@ -25,7 +25,7 @@ class VisualizationMarginModel(unittest.TestCase): # MarginModel Linear with respect to the shape (from 0.01 to 0.02) margin_model = cls.margin_model(coordinates=coordinates, params_sample={GevParams.GEV_SHAPE: 0.02}) if cls.DISPLAY: - margin_model.margin_function_sample.visualize_all() + margin_model.margin_function_sample.visualize() if __name__ == '__main__': diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py index 090f8f3385b8cf041029680d67e25fb5aa4ada63..8fcbc504491a971e0036fa7fff76c3db965ada5c 100644 --- a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py +++ b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py @@ -11,6 +11,7 @@ class AbstractMarginFunction(object): def __init__(self, coordinates: AbstractCoordinates): self.coordinates = coordinates + self.visualization_axes = None def get_gev_params(self, coordinate: np.ndarray) -> GevParams: """Main method that maps each coordinate to its GEV parameters""" @@ -18,17 +19,20 @@ class AbstractMarginFunction(object): # Visualization function - def visualize_all(self): - fig, axes = plt.subplots(3, 1, sharex='col', sharey='row') - fig.subplots_adjust(hspace=0.4, wspace=0.4, ) + def visualize(self, axes=None, show=True): + if axes is None: + fig, axes = plt.subplots(3, 1, sharex='col', sharey='row') + fig.subplots_adjust(hspace=0.4, wspace=0.4, ) + self.visualization_axes = axes for i, gev_param_name in enumerate(GevParams.GEV_PARAM_NAMES): ax = axes[i] - self.visualize(gev_param_name, ax, show=False) + self.visualize_single_param(gev_param_name, ax, show=False) title_str = gev_param_name ax.set_title(title_str) - plt.show() + if show: + plt.show() - def visualize(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=True): + def visualize_single_param(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=True): if self.coordinates.nb_columns == 1: self.visualize_1D(gev_param_name, ax, show) elif self.coordinates.nb_columns == 2: diff --git a/extreme_estimator/extreme_models/margin_model/param_function/param_function.py b/extreme_estimator/extreme_models/margin_model/param_function/param_function.py index df9f485a2c25c99e53fadb4e54f134407f9faeff..ad8555ef2d0bc5bc4e7bd33de8cf66eb49889ee4 100644 --- a/extreme_estimator/extreme_models/margin_model/param_function/param_function.py +++ b/extreme_estimator/extreme_models/margin_model/param_function/param_function.py @@ -26,11 +26,13 @@ class LinearOneAxisParamFunction(ParamFunction): self.t_max = coordinates[:, linear_axis].max() self.coef = coef + def get_gev_param_value_normalized(self, coordinate: np.ndarray) -> float: + return self.get_gev_param_value(coordinate) / (self.t_max - self.t_min) + def get_gev_param_value(self, coordinate: np.ndarray) -> float: t = coordinate[self.linear_axis] - t_between_zero_and_one = t / (self.t_max - self.t_min) - assert -1 <= t_between_zero_and_one <= 1, 'Out of bounds' - return t_between_zero_and_one * self.coef + assert self.t_min <= t <= self.t_max, 'Out of bounds' + return t * self.coef class LinearParamFunction(ParamFunction): @@ -40,7 +42,7 @@ class LinearParamFunction(ParamFunction): # Load each one axis linear function self.linear_one_axis_param_functions = [] # type: List[LinearOneAxisParamFunction] for linear_dim in linear_dims: - param_function = LinearOneAxisParamFunction(linear_axis=linear_dim-1, coordinates=coordinates, + param_function = LinearOneAxisParamFunction(linear_axis=linear_dim - 1, coordinates=coordinates, coef=self.linear_coef.get_coef(dim=linear_dim)) self.linear_one_axis_param_functions.append(param_function) diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py index 83e5f6d0ebd99f6be2833fc20d9039a8f0b606dd..b1115f6ca1d46b4e8f7d28f55bad4861f603a4f1 100644 --- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py +++ b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py @@ -12,31 +12,41 @@ class LinearMarginModel(AbstractMarginModel): def load_margin_functions(self, gev_param_name_to_linear_dims=None): # Load sample coef - self.default_params_sample = GevParams(1.0, 1.0, 1.0).to_dict() - linear_coef_sample = self.get_standard_linear_coef(gev_param_name_to_intercept=self.params_sample) + self.default_params_sample = self.default_param_name_and_dim_to_coef() + linear_coef_sample = self.gev_param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_sample) self.margin_function_sample = LinearMarginFunction(coordinates=self.coordinates, gev_param_name_to_linear_coef=linear_coef_sample, gev_param_name_to_linear_dims=gev_param_name_to_linear_dims) # Load start fit coef - self.default_params_start_fit = GevParams(1.0, 1.0, 1.0).to_dict() - linear_coef_start_fit = self.get_standard_linear_coef(gev_param_name_to_intercept=self.params_start_fit) + self.default_params_start_fit = self.default_param_name_and_dim_to_coef() + linear_coef_start_fit = self.gev_param_name_to_linear_coef(param_name_and_dim_to_coef=self.params_start_fit) self.margin_function_start_fit = LinearMarginFunction(coordinates=self.coordinates, gev_param_name_to_linear_coef=linear_coef_start_fit, gev_param_name_to_linear_dims=gev_param_name_to_linear_dims) @staticmethod - def get_standard_linear_coef(gev_param_name_to_intercept, slope=0.1): + def default_param_name_and_dim_to_coef() -> dict: + default_intercept = 1 + default_slope = 0.01 + gev_param_name_and_dim_to_coef = {} + for gev_param_name in GevParams.GEV_PARAM_NAMES: + gev_param_name_and_dim_to_coef[(gev_param_name, 0)] = default_intercept + for dim in [1, 2, 3]: + gev_param_name_and_dim_to_coef[(gev_param_name, dim)] = default_slope + return gev_param_name_and_dim_to_coef + + @staticmethod + def gev_param_name_to_linear_coef(param_name_and_dim_to_coef): gev_param_name_to_linear_coef = {} for gev_param_name in GevParams.GEV_PARAM_NAMES: - dim_to_coef = {dim: slope for dim in range(1, 4)} - dim_to_coef[0] = gev_param_name_to_intercept[gev_param_name] + dim_to_coef = {dim: param_name_and_dim_to_coef[(gev_param_name, dim)] for dim in [0, 1, 2, 3]} linear_coef = LinearCoef(gev_param_name=gev_param_name, dim_to_coef=dim_to_coef) gev_param_name_to_linear_coef[gev_param_name] = linear_coef return gev_param_name_to_linear_coef - - def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates_values: np.ndarray) -> AbstractMarginFunction: + def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, + coordinates_values: np.ndarray) -> AbstractMarginFunction: return self.margin_function_start_fit @@ -46,19 +56,19 @@ class ConstantMarginModel(LinearMarginModel): super().load_margin_functions({}) -class LinearShapeAxis0MarginModel(LinearMarginModel): +class LinearShapeDim1MarginModel(LinearMarginModel): def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None): super().load_margin_functions({GevParams.GEV_SHAPE: [1]}) -class LinearShapeAxis0and1MarginModel(LinearMarginModel): +class LinearShapeDim1and2MarginModel(LinearMarginModel): def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None): super().load_margin_functions({GevParams.GEV_SHAPE: [1, 2]}) -class LinearAllParametersAxis0MarginModel(LinearMarginModel): +class LinearAllParametersDim1MarginModel(LinearMarginModel): def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None): super().load_margin_functions({GevParams.GEV_SHAPE: [1], @@ -66,7 +76,7 @@ class LinearAllParametersAxis0MarginModel(LinearMarginModel): GevParams.GEV_SCALE: [1]}) -class LinearAllParametersAllAxisMarginModel(LinearMarginModel): +class LinearAllParametersAllDimsMarginModel(LinearMarginModel): def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None): all_dims = list(range(1, self.coordinates.nb_columns + 1)) diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py index 513d50c7506a84f768a2a33082d9065a9b66f26f..9b1cb54862d9e080822ab126a93f93199688d3e9 100644 --- a/extreme_estimator/extreme_models/utils.py +++ b/extreme_estimator/extreme_models/utils.py @@ -1,4 +1,6 @@ import os.path as op +import random +import sys import rpy2.robjects as ro from rpy2.robjects import numpy2ri @@ -10,6 +12,9 @@ def get_loaded_r() -> ro.R: numpy2ri.activate() pandas2ri.activate() r.library('SpatialExtremes') + # max_int = r('.Machine$integer.max') + # seed = random.randrange(max_int) + # r("set.seed({})".format(seed)) return r diff --git a/test/test_utils.py b/test/test_utils.py index 84ec34eee67f588876f86c4e507311fe744b7f44..5abfc72ae43bd69bd8c3bd5945081f7d5c8f2fb2 100644 --- a/test/test_utils.py +++ b/test/test_utils.py @@ -1,7 +1,7 @@ from extreme_estimator.estimator.full_estimator import SmoothMarginalsThenUnitaryMsp, \ FullEstimatorInASingleStepWithSmoothMargin from extreme_estimator.estimator.max_stable_estimator import MaxStableEstimator -from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearAllParametersAllAxisMarginModel, \ +from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearAllParametersAllDimsMarginModel, \ ConstantMarginModel from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import \ AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction @@ -21,7 +21,7 @@ In this case, unit test (at least on the constructor) must be ensured in the tes TEST_MAX_STABLE_MODEL = [Smith, BrownResnick, Schlather, Geometric, ExtremalT, ISchlather] TEST_1D_AND_2D_COORDINATES = [UniformCoordinates, CircleCoordinates] TEST_3D_COORDINATES = [AlpsStation3DCoordinatesWithAnisotropy] -TEST_MARGIN_TYPES = [ConstantMarginModel, LinearAllParametersAllAxisMarginModel][:] +TEST_MARGIN_TYPES = [ConstantMarginModel, LinearAllParametersAllDimsMarginModel][:] TEST_MAX_STABLE_ESTIMATOR = [MaxStableEstimator] TEST_FULL_ESTIMATORS = [SmoothMarginalsThenUnitaryMsp, FullEstimatorInASingleStepWithSmoothMargin][:]