[MARGIN FUNCTION] visualize several margin functions on the same graph

experiment/regression_margin/regression_margin.py

+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()

experiment/regression_margin/visualization_margin_model.py

 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__':

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:

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)
 

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))

extreme_estimator/extreme_models/utils.py

 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
 
 

test/test_utils.py

 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][:]