[FIX] fix fitmaxstable, add unitary tests

example.py

-import numpy as np
-import pandas as pd
-
-from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearMarginModel, \
-    LinearAllParametersAllDimsMarginModel
-from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
-from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Schlather
-from extreme_estimator.extreme_models.utils import r, set_seed_r
-from extreme_estimator.gev_params import GevParams
-from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
-from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
-from spatio_temporal_dataset.temporal_observations.annual_maxima_observations import MaxStableAnnualMaxima
-
-print('R')
-set_seed_r()
-r("""
-n.site <- 30
-locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
-colnames(locations) <- c("lon", "lat")""")
-
-for _ in range(2):
-    set_seed_r()
-    r("""data <- rmaxstab(40, locations, cov.mod = "whitmat", nugget = 0, range = 3,
-    smooth = 0.5)""")
-    print(np.sum(r.data))
-    r("""
-    ##param.loc <- -10 + 2 * locations[,2]
-    ##TODO-IMPLEMENT SQUARE: param.scale <- 5 + 2 * locations[,1] + locations[,2]^2
-    ##param.scale <- 5 + 2 * locations[,1] + locations[,2]
-    param.shape <- rep(0.2, n.site)
-    param.loc <- rep(0.2, n.site)
-    param.scale <- rep(0.2, n.site)
-    ##Transform the unit Frechet margins to GEV
-    for (i in 1:n.site)
-    data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i], param.shape[i])""")
-    print(np.sum(r.data))
-
-print('\n\nPython')
-
-params_sample = {'range': 3, 'smooth': 0.5, 'nugget': 0.0}
-max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, params_sample=params_sample)
-df = pd.DataFrame(data=r.locations, columns=AbstractCoordinates.COORDINATE_NAMES[:2])
-coordinates = AbstractCoordinates.from_df(df)
-set_seed_r()
-maxima = MaxStableAnnualMaxima.from_sampling(nb_obs=40, max_stable_model=max_stable_model, coordinates=coordinates)
-print(np.sum(maxima.maxima_frech))
-
-# gev_param_name_to_coef_list = {
-#     GevParams.GEV_LOC: [-10, 0, 2],
-#     GevParams.GEV_SHAPE: [5, 2, 1],
-#     GevParams.GEV_SCALE: [0.2, 0, 0],
-# }
-gev_param_name_to_coef_list = {
-    GevParams.GEV_LOC: [0.2, 0, 0],
-    GevParams.GEV_SHAPE: [0.2, 0, 0],
-    GevParams.GEV_SCALE: [0.2, 0, 0],
-}
-margin_model = LinearAllParametersAllDimsMarginModel.from_coef_list(coordinates, gev_param_name_to_coef_list)
-maxima_gev = margin_model.frech2gev(maxima.maxima_frech, coordinates.coordinates_values, margin_model.margin_function_sample)
-print(np.sum(maxima_gev))
-
-# dataset = FullSimulatedDataset.from_double_sampling(nb_obs=nb_obs, margin_model=margin_model,
-#                                                     coordinates=coordinates,
-#                                                     max_stable_model=max_stable_model)
-

experiment/regression_margin/regression_margin.py

@@ -5,7 +5,7 @@ import numpy as np
 from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStepWithSmoothMargin
 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
+    LinearAllParametersAllDimsMarginModel, LinearScaleDim1MarginModel, ConstantMarginModel
 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
@@ -13,9 +13,9 @@ import matplotlib.pyplot as plt
 
 from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
 
-nb_points = 5
-nb_obs = 10
-nb_estimator = 2
+nb_points = 2
+nb_obs = 10000
+nb_estimator = 1
 show = False
 
 coordinates = LinSpaceCoordinates.from_nb_points(nb_points=nb_points)
@@ -24,33 +24,39 @@ coordinates = LinSpaceCoordinates.from_nb_points(nb_points=nb_points)
 
 # MarginModel Linear with respect to the shape (from 0.01 to 0.02)
 params_sample = {
-    (GevParams.GEV_SHAPE, 0): 0.2,
-    (GevParams.GEV_SHAPE, 1): 0.05,
+    # (GevParams.GEV_SHAPE, 0): 0.2,
+    (GevParams.GEV_LOC, 0): 10,
+    (GevParams.GEV_SHAPE, 0): 1.0,
+    (GevParams.GEV_SCALE, 0): 1.0,
 }
-margin_model = LinearShapeDim1MarginModel(coordinates=coordinates, params_sample=params_sample)
+margin_model = ConstantMarginModel(coordinates=coordinates, params_sample=params_sample)
+margin_model_for_estimator_class = [LinearAllParametersAllDimsMarginModel][-1]
 max_stable_model = Smith()
 
-# 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 #################
 
 
 axes = None
-for _ in range(nb_estimator):
+for i in range(nb_estimator):
+    print("{}/{}".format(i+1, 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)
+
+    if show and i == 0:
+        # Plot a realization from the maxima gev (i.e the maxima obtained just by simulating the marginal law)
+        for maxima in np.transpose(dataset.maxima_frech):
+            plt.plot(coordinates.coordinates_values, maxima, 'o')
+        plt.show()
+
     margin_function_sample = dataset.margin_model.margin_function_sample # type: LinearMarginFunction
-    margin_function_sample.visualize(show=False, axes=axes)
+    margin_function_sample.visualize(show=False, axes=axes, dot_display=True)
     axes = margin_function_sample.visualization_axes
 
     # Estimation part
-    margin_model_for_estimator = LinearAllParametersAllDimsMarginModel(coordinates)
+    margin_model_for_estimator = margin_model_for_estimator_class(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)

extreme_estimator/estimator/full_estimator.py

@@ -68,7 +68,7 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
     def _fit(self):
         # Estimate both the margin and the max-stable structure
         self.full_params_fitted = self.max_stable_model.fitmaxstab(
-            maxima_frech=self.dataset.maxima_frech,
+            maxima_gev=self.dataset.maxima_gev,
             df_coordinates=self.dataset.df_coordinates,
             fit_marge=True,
             fit_marge_form_dict=self.linear_margin_function_to_fit.form_dict,

extreme_estimator/extreme_models/margin_model/smooth_margin_model.py

@@ -92,6 +92,14 @@ class LinearAllParametersDim1MarginModel(LinearMarginModel):
                                        GevParams.GEV_SCALE: [1]})
 
 
+class LinearMarginModelExample(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],
+                                       GevParams.GEV_LOC: [2],
+                                       GevParams.GEV_SCALE: [1]})
+
+
 class LinearAllParametersAllDimsMarginModel(LinearMarginModel):
 
     def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):

extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py

@@ -3,11 +3,11 @@ from enum import Enum
 
 import numpy as np
 import rpy2
-from rpy2.rinterface import RRuntimeError
+from rpy2.rinterface import RRuntimeError, RRuntimeWarning
 import rpy2.robjects as robjects
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
-from extreme_estimator.extreme_models.utils import r
+from extreme_estimator.extreme_models.utils import r, set_seed_r
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
@@ -16,21 +16,23 @@ class AbstractMaxStableModel(AbstractModel):
     def __init__(self, params_start_fit=None, params_sample=None):
         super().__init__(params_start_fit, params_sample)
         self.cov_mod = None
+        self.start_value_for_fitmaxstab = True
 
     @property
     def cov_mod_param(self):
         return {'cov.mod': self.cov_mod}
 
-    def fitmaxstab(self, maxima_frech: np.ndarray, df_coordinates: pd.DataFrame, fit_marge=False,
+    def fitmaxstab(self, df_coordinates: pd.DataFrame, maxima_frech: np.ndarray=None, maxima_gev: np.ndarray=None, fit_marge=False,
                    fit_marge_form_dict=None, margin_start_dict=None) -> dict:
-        assert isinstance(maxima_frech, np.ndarray)
         assert isinstance(df_coordinates, pd.DataFrame)
         if fit_marge:
             assert fit_marge_form_dict is not None
             assert margin_start_dict is not None
 
         # Prepare the data
-        data = np.transpose(maxima_frech)
+        maxima = maxima_gev if fit_marge else maxima_frech
+        assert isinstance(maxima, np.ndarray)
+        data = np.transpose(maxima)
 
         # Prepare the coord
         df_coordinates = df_coordinates.copy()
@@ -54,14 +56,19 @@ class AbstractMaxStableModel(AbstractModel):
             fit_params.update({k: robjects.Formula(v) for k, v in fit_marge_form_dict.items()})
         if fitmaxstab_with_one_dimensional_data:
             fit_params['iso'] = True
-        fit_params['start'] = r.list(**start_dict)
+        if self.start_value_for_fitmaxstab:
+            fit_params['start'] = r.list(**start_dict)
         fit_params['fit.marge'] = fit_marge
 
         # Run the fitmaxstab in R
         try:
             res = r.fitmaxstab(data=data, coord=coord, **fit_params)  # type: robjects.ListVector
-        except RRuntimeError as error:
-            raise Exception('Some R exception have been launched at RunTime: \n {}'.format(error.__repr__()))
+        except (RRuntimeError, RRuntimeWarning) as e:
+            if 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')
         # todo: maybe if the convergence was not successful I could try other starting point several times
         # Retrieve the resulting fitted values
         fitted_values = res.rx2('fitted.values')

test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py

+import unittest
+
+from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
+    LinearMarginModelExample
+from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Schlather
+from extreme_estimator.extreme_models.utils import r
+from test.test_unitary.test_rmaxstab.test_rmaxstab_with_margin import TestRMaxStabWithMarginConstant
+from test.test_unitary.test_unitary_abstract import TestUnitaryAbstract
+
+
+class TestMaxStableFitWithConstantMargin(TestUnitaryAbstract):
+
+    @property
+    def r_output(self):
+        TestRMaxStabWithMarginConstant.r_code()
+        r("""
+        shape.form = shape ~ 1
+        loc.form = loc ~ 1
+        scale.form = scale ~ 1
+        res = fitmaxstab(data, locations, "whitmat", loc.form, scale.form, shape.form)
+        """)
+        return self.r_fitted_values_from_res_variable
+
+    @property
+    def python_output(self):
+        dataset = TestRMaxStabWithMarginConstant.python_code()
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat)
+        margin_model = ConstantMarginModel(dataset.coordinates)
+        max_stable_model.start_value_for_fitmaxstab = False
+        full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
+                                                                    max_stable_model)
+        full_estimator.fit()
+        return full_estimator.full_params_fitted
+
+    def test_max_stable_fit_with_constant_margin(self):
+        self.compare()
+
+
+class TestMaxStableFitWithLinearMargin(TestUnitaryAbstract):
+
+    @property
+    def r_output(self):
+        TestRMaxStabWithMarginConstant.r_code()
+        r("""
+        loc.form <- loc ~ lat
+        ##scale.form <- scale ~ lon + I(lat^2)
+        scale.form <- scale ~ lon
+        shape.form <- shape ~ lon
+        res = fitmaxstab(data, locations, "whitmat", loc.form, scale.form, shape.form)
+        """)
+        return self.r_fitted_values_from_res_variable
+
+    @property
+    def python_output(self):
+        dataset = TestRMaxStabWithMarginConstant.python_code()
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat)
+        margin_model = LinearMarginModelExample(dataset.coordinates)
+        max_stable_model.start_value_for_fitmaxstab = False
+        full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
+                                                                    max_stable_model)
+        full_estimator.fit()
+        return full_estimator.full_params_fitted
+
+    def test_max_stable_fit_with_linear_margin(self):
+        self.compare()
+
+
+if __name__ == '__main__':
+    unittest.main()

test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py

+import unittest
+
+from extreme_estimator.estimator.max_stable_estimator import MaxStableEstimator
+from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Schlather
+from extreme_estimator.extreme_models.utils import r
+from spatio_temporal_dataset.dataset.simulation_dataset import MaxStableDataset
+from test.test_unitary.test_rmaxstab.test_rmaxstab_without_margin import TestRMaxStab
+from test.test_unitary.test_unitary_abstract import TestUnitaryAbstract
+
+
+class TestMaxStableFitWithoutMargin(TestUnitaryAbstract):
+
+    @property
+    def r_output(self):
+        TestRMaxStab.r_code()
+        r("""res = fitmaxstab(data, locations, "whitmat")""")
+        return self.r_fitted_values_from_res_variable
+
+    @property
+    def python_output(self):
+        coordinates, max_stable_model = TestRMaxStab.python_code()
+        dataset = MaxStableDataset.from_sampling(nb_obs=40, max_stable_model=max_stable_model, coordinates=coordinates)
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat)
+        max_stable_model.start_value_for_fitmaxstab = False
+        max_stable_estimator = MaxStableEstimator(dataset, max_stable_model)
+        max_stable_estimator.fit()
+        return max_stable_estimator.max_stable_params_fitted
+
+    def test_max_stable_fit_without_margin(self):
+        self.compare()
+
+
+if __name__ == '__main__':
+    unittest.main()

test/test_unitary/test_rmaxstab/test_rmaxstab_with_margin.py

+import unittest
+
+import numpy as np
+
+from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
+    LinearAllParametersAllDimsMarginModel
+from extreme_estimator.extreme_models.utils import r
+from extreme_estimator.gev_params import GevParams
+from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
+from test.test_unitary.test_rmaxstab.test_rmaxstab_without_margin import TestRMaxStab
+from test.test_unitary.test_unitary_abstract import TestUnitaryAbstract
+
+
+class TestRMaxStabWithMarginConstant(TestUnitaryAbstract):
+
+    @classmethod
+    def r_code(cls):
+        TestRMaxStab.r_code()
+        r("""
+        param.shape <- rep(0.2, n.site)
+        param.loc <- rep(0.2, n.site)
+        param.scale <- rep(0.2, n.site)""")
+        r("""
+        for (i in 1:n.site)
+        data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i], param.shape[i])
+        """)
+
+
+    @classmethod
+    def python_code(cls):
+        coordinates, max_stable_model = TestRMaxStab.python_code()
+        # Load margin model
+        gev_param_name_to_coef_list = {
+            GevParams.GEV_LOC: [0.2],
+            GevParams.GEV_SHAPE: [0.2],
+            GevParams.GEV_SCALE: [0.2],
+        }
+        margin_model = ConstantMarginModel.from_coef_list(coordinates, gev_param_name_to_coef_list)
+        # Load dataset
+        dataset = FullSimulatedDataset.from_double_sampling(nb_obs=40, margin_model=margin_model,
+                                                            coordinates=coordinates,
+                                                            max_stable_model=max_stable_model)
+
+        return dataset
+
+    @property
+    def r_output(self):
+        self.r_code()
+        return np.sum(r.data)
+
+    @property
+    def python_output(self):
+        dataset = self.python_code()
+        return np.sum(dataset.maxima_gev)
+
+    def test_rmaxstab_with_constant_margin(self):
+        self.compare()
+
+
+class TestRMaxStabWithLinearMargin(TestUnitaryAbstract):
+
+    @classmethod
+    def r_code(cls):
+        TestRMaxStab.r_code()
+        r("""
+        param.loc <- -10 + 2 * locations[,2]
+        ##TODO-IMPLEMENT SQUARE: param.scale <- 5 + 2 * locations[,1] + locations[,2]^2
+        param.shape <- 2 -3 * locations[,1]
+        param.scale <- 5 + 2 * locations[,1] + locations[,2]
+        ##Transform the unit Frechet margins to GEV
+        for (i in 1:n.site)
+        data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i], param.shape[i])
+        """)
+
+    @classmethod
+    def python_code(cls):
+        coordinates, max_stable_model = TestRMaxStab.python_code()
+        # Load margin model
+        gev_param_name_to_coef_list = {
+            GevParams.GEV_LOC: [-10, 0, 2],
+            GevParams.GEV_SHAPE: [2, -3, 0],
+            GevParams.GEV_SCALE: [5, 2, 1],
+        }
+        margin_model = LinearAllParametersAllDimsMarginModel.from_coef_list(coordinates, gev_param_name_to_coef_list)
+        # Load dataset
+        dataset = FullSimulatedDataset.from_double_sampling(nb_obs=40, margin_model=margin_model,
+                                                            coordinates=coordinates,
+                                                            max_stable_model=max_stable_model)
+        return dataset
+
+    @property
+    def r_output(self):
+        self.r_code()
+        return np.sum(r.data)
+
+    @property
+    def python_output(self):
+        dataset = self.python_code()
+        return np.sum(dataset.maxima_gev)
+
+    def test_rmaxstab_with_linear_margin(self):
+        self.compare()
+
+
+if __name__ == '__main__':
+    unittest.main()

test/test_unitary/test_rmaxstab/test_rmaxstab_without_margin.py

+import unittest
+
+import numpy as np
+import pandas as pd
+
+from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Schlather
+from extreme_estimator.extreme_models.utils import r
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.temporal_observations.annual_maxima_observations import MaxStableAnnualMaxima
+from test.test_unitary.test_unitary_abstract import TestUnitaryAbstract
+
+
+class TestRMaxStab(TestUnitaryAbstract):
+
+    @classmethod
+    def r_code(cls):
+        r("""data <- rmaxstab(40, locations, cov.mod = "whitmat", nugget = 0, range = 3, smooth = 0.5)""")
+
+    @classmethod
+    def python_code(cls):
+        # Load coordinate object
+        df = pd.DataFrame(data=r.locations, columns=AbstractCoordinates.COORDINATE_NAMES[:2])
+        coordinates = AbstractCoordinates.from_df(df)
+        # Load max stable model
+        params_sample = {'range': 3, 'smooth': 0.5, 'nugget': 0}
+        max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, params_sample=params_sample)
+        return coordinates, max_stable_model
+
+    @property
+    def r_output(self):
+        self.r_code()
+        return np.sum(r.data)
+
+    @property
+    def python_output(self):
+        coordinates, max_stable_model = self.python_code()
+        m = MaxStableAnnualMaxima.from_sampling(nb_obs=40, max_stable_model=max_stable_model, coordinates=coordinates)
+        # TODO: understand why the array are not in the same order
+        return np.sum(m.maxima_frech)
+
+    def test_rmaxstab(self):
+        self.compare()
+
+
+if __name__ == '__main__':
+    unittest.main()

test/test_unitary/test_unitary_abstract.py

+import unittest
+
+from extreme_estimator.extreme_models.utils import set_seed_r, r
+
+
+class TestUnitaryAbstract(unittest.TestCase):
+
+    def setUp(self):
+        r("""
+        n.site <- 30
+        locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
+        colnames(locations) <- c("lon", "lat")
+        """)
+
+    @property
+    def r_fitted_values_from_res_variable(self):
+        res = r.res
+        fitted_values = res.rx2('fitted.values')
+        fitted_values = {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
+        return fitted_values
+
+    @property
+    def r_output(self):
+        pass
+
+    @property
+    def python_output(self):
+        pass
+
+    def compare(self):
+        set_seed_r()
+        r_output = self.r_output
+        set_seed_r()
+        python_output = self.python_output
+        self.assertEqual(r_output, python_output, msg="python: {}, r: {}".format(python_output, r_output))
+
+