[EXTREME FIT][MAX STABLE FIT] add 6 main types of max stable model. add test

extreme_estimator/R_fit/gev_fit/gev_mle_fit.py

@@ -16,7 +16,7 @@ def mle_gev(x_gev: np.ndarray, start_loc=0, start_scale=1, start_shape=0):
     x_gev = rpyn.numpy2ri(x_gev)
     r.assign('x_gev', x_gev)
     r.assign('python_wrapping', True)
-    r.source(file="/home/erwan/Documents/projects/spatiotemporalextremes/R/gev_fit/gev_fit.R")
+    r.source(file="/home/erwan/Documents/projects/spatiotemporalextremes/extreme_estimator/gev_fit/gev_fit.extreme_estimator")
     res = r.mle_gev(loc=start_loc, scale=start_scale, shape=start_shape)
     mle_params = dict(r.attr(res, 'coef').items())
     return mle_params

extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py

+from enum import Enum
+
+from rpy2.robjects import ListVector
+from extreme_estimator.R_fit.utils import get_loaded_r
+import numpy as np
+
+
+class AbstractMaxStableModel(object):
+
+    def __init__(self, params_start_fit=None, params_sample=None):
+        self.cov_mod = None
+        self.default_params_start_fit = None
+        self.default_params_sample = None
+        self.user_params_start_fit = params_start_fit
+        self.user_params_sample = params_sample
+        self.r = get_loaded_r()
+
+    def fitmaxstab(self, maxima: np.ndarray, coord: np.ndarray, ):
+        # todo: find how to specify a startign point, understand how is it set by default
+        # res = None
+        # tries = 5
+        # nb_tries = 0
+        # while res is None and nb_tries < tries:
+        #     try:
+        #         res = self.r.fitmaxstab(maxima, coord, **self.cov_mod_param)  # type: ListVector
+        #     except rpy2.rinterface.RRuntimeError:
+        #         pass
+        #     nb_tries += 1
+        #     print(nb_tries)
+        res = self.r.fitmaxstab(np.transpose(maxima), coord, **self.cov_mod_param)  # type: ListVector
+        fitted_values = res.rx2('fitted.values')
+        fitted_values = {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
+        return fitted_values
+
+    def rmaxstab(self, nb_obs: int, coord: np.ndarray, ) -> np.ndarray:
+        """
+        Return an numpy of maxima. With rows being the stations and columns being the years of maxima
+        """
+        maxima = np.array(self.r.rmaxstab(nb_obs, coord, **self.cov_mod_param, **self.params_sample))
+        return np.transpose(maxima)
+
+    @property
+    def cov_mod_param(self):
+        return {'cov.mod': self.cov_mod}
+
+    @property
+    def params_start_fit(self):
+        return self.merge_params(default_params=self.default_params_start_fit, input_params=self.user_params_start_fit)
+
+    @property
+    def params_sample(self):
+        return self.merge_params(default_params=self.default_params_sample, input_params=self.user_params_sample)
+
+    @staticmethod
+    def merge_params(default_params, input_params):
+        merged_params = default_params.copy()
+        if input_params is not None:
+            assert isinstance(default_params, dict) and isinstance(input_params, dict)
+            assert set(input_params.keys()).issubset(set(default_params.keys()))
+            merged_params.update(input_params)
+        return merged_params
+
+
+class CovarianceFunction(Enum):
+    whitmat = 0
+    cauchy = 1
+    powexp = 2
+    bessel = 3
+
+
+class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
+
+    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
+        super().__init__(params_start_fit, params_sample)
+        assert covariance_function is not None
+        self.covariance_function = covariance_function
+        self.default_params_sample = {
+            'range': 3,
+            'smooth': 0.5,
+            'nugget': 0.5
+        }
+
+
+if __name__ == '__main__':
+     print([c.name for c in CovarianceFunction])
+     for covariance_function in CovarianceFunction:
+         print(covariance_function)
\ No newline at end of file

extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R

@@ -13,11 +13,14 @@ if (call_main) {
     print(coord)
     #  Generate the data
     # data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220)
-    data <- rmaxstab(n.obs, coord, "brown", range = 3, smooth = 0.5)
+    # data <- rmaxstab(n.obs, coord, "brown", range = 3, smooth = 0.5)
+    # data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0.0, range = 3, smooth = 0.5)
+    data <- rmaxstab(n.obs, coord, "cauchy", )
     #  Fit back the data
     print(data)
     # res = fitmaxstab(data, coord, "gauss", fit.marge=FALSE, )
-    res = fitmaxstab(data, coord, "brown")
+    # res = fitmaxstab(data, coord, "brown")
+    res = fitmaxstab(data, coord, "cauchy")
     # res = fitmaxstab(data, coord, "gauss", start=list(0,0,0))
     print(res)
     print(class(res))

extreme_estimator/R_fit/max_stable_fit/max_stable_models.py

-import rpy2
+from enum import Enum
 
-from rpy2.robjects import ListVector
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel, \
+    AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction
 
-from extreme_estimator.R_fit.utils import get_loaded_r
-import numpy as np
 
-
-def optimize_fit(max_stable_model_list, params_start_fit_list):
-    pass
-
-
-class MaxStableModel(object):
-
-    def __init__(self, params_start_fit=None, params_sample=None):
-        self.cov_mod = None
-        self.default_params_start_fit = None
-        self.default_params_sample = None
-        self.user_params_start_fit = params_start_fit
-        self.user_params_sample = params_sample
-        self.r = get_loaded_r()
-
-    def fitmaxstab(self, maxima: np.ndarray, coord: np.ndarray, ):
-        # todo: find how to specify a startign point, understand how is it set by default
-        # res = None
-        # tries = 5
-        # nb_tries = 0
-        # while res is None and nb_tries < tries:
-        #     try:
-        #         res = self.r.fitmaxstab(maxima, coord, **self.cov_mod_param)  # type: ListVector
-        #     except rpy2.rinterface.RRuntimeError:
-        #         pass
-        #     nb_tries += 1
-        #     print(nb_tries)
-        res = self.r.fitmaxstab(np.transpose(maxima), coord, **self.cov_mod_param)  # type: ListVector
-        fitted_values = res.rx2('fitted.values')
-        fitted_values = {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
-        return fitted_values
-
-    def rmaxstab(self, nb_obs: int, coord: np.ndarray, ) -> np.ndarray:
-        """
-        Return an numpy of maxima. With rows being the stations and columns being the years of maxima
-        """
-        maxima = np.array(self.r.rmaxstab(nb_obs, coord, **self.cov_mod_param, **self.params_sample))
-        return np.transpose(maxima)
-
-    @property
-    def cov_mod_param(self):
-        return {'cov.mod': self.cov_mod}
-
-    @property
-    def params_start_fit(self):
-        return self.merge_params(default_params=self.default_params_start_fit, input_params=self.user_params_start_fit)
-
-    @property
-    def params_sample(self):
-        return self.merge_params(default_params=self.default_params_sample, input_params=self.user_params_sample)
-
-    @staticmethod
-    def merge_params(default_params, input_params):
-        merged_params = default_params.copy()
-        if input_params is not None:
-            assert isinstance(default_params, dict) and isinstance(input_params, dict)
-            assert set(input_params.keys()).issubset(set(default_params.keys()))
-            merged_params.update(input_params)
-        return merged_params
-
-
-class GaussianMSP(MaxStableModel):
+class Smith(AbstractMaxStableModel):
 
     def __init__(self, params_start_fit=None, params_sample=None):
         super().__init__(params_start_fit=params_start_fit, params_sample=params_sample)
@@ -79,17 +17,45 @@ class GaussianMSP(MaxStableModel):
         }
 
 
-class BrownResick(MaxStableModel):
+class BrownResnick(AbstractMaxStableModel):
 
     def __init__(self, params_start_fit=None, params_sample=None):
         super().__init__(params_start_fit=params_start_fit, params_sample=params_sample)
         self.cov_mod = 'brown'
         self.default_params_start_fit = {}
         self.default_params_sample = {
-            "range": 3,
-            "smooth": 0.5,
+            'range': 3,
+            'smooth': 0.5,
         }
 
-class ExtremalT(MaxStableModel):
-    pass
 
+class Schlather(AbstractMaxStableModelWithCovarianceFunction):
+
+    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
+        super().__init__(params_start_fit, params_sample, covariance_function)
+        self.cov_mod = self.covariance_function.name
+        self.default_params_sample.update({})
+
+
+class Geometric(AbstractMaxStableModelWithCovarianceFunction):
+
+    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
+        super().__init__(params_start_fit, params_sample, covariance_function)
+        self.cov_mod = 'g' + self.covariance_function.name
+        self.default_params_sample .update({'sigma2': 0.5})
+
+
+class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
+
+    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
+        super().__init__(params_start_fit, params_sample, covariance_function)
+        self.cov_mod = 't' + self.covariance_function.name
+        self.default_params_sample .update({'DoF': 2})
+
+
+class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
+
+    def __init__(self, params_start_fit=None, params_sample=None, covariance_function: CovarianceFunction = None):
+        super().__init__(params_start_fit, params_sample, covariance_function)
+        self.cov_mod = 'i' + self.covariance_function.name
+        self.default_params_sample .update({'alpha': 0.5})

extreme_estimator/estimator/msp_estimator.py

 from extreme_estimator.estimator.abstract_estimator import AbstractEstimator
-from extreme_estimator.R_fit.max_stable_fit.max_stable_models import MaxStableModel
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 import numpy as np
 
 
 class MaxStableEstimator(AbstractEstimator):
 
-    def __init__(self, dataset: AbstractDataset, max_stable_model: MaxStableModel):
+    def __init__(self, dataset: AbstractDataset, max_stable_model: AbstractMaxStableModel):
         self.dataset = dataset
         self.max_stable_model = max_stable_model
         self.max_stable_params_fitted = None

extreme_estimator/robustness_plot/abstract_robustness.py

 from typing import List
 
 from extreme_estimator.estimator.msp_estimator import MaxStableEstimator
-from extreme_estimator.R_fit.max_stable_fit.max_stable_models import GaussianMSP, MaxStableModel
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import GaussianMSP, AbstractMaxStableModel
 from itertools import product
 
 from spatio_temporal_dataset.dataset.simulation_dataset import SimulatedDataset

extreme_estimator/robustness_plot/spatial_robustness.py

-from extreme_estimator.R_fit.max_stable_fit.max_stable_models import GaussianMSP, BrownResick
 from extreme_estimator.robustness_plot.abstract_robustness import DisplayItem, AbstractRobustnessPlot, \
     SpatialCoordinateClassItem, NbObservationItem, NbStationItem, MaxStableModelItem, SpatialParamsItem
 from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinates

spatio_temporal_dataset/dataset/simulation_dataset.py

-from extreme_estimator.R_fit.max_stable_fit.max_stable_models import MaxStableModel, GaussianMSP
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import AbstractMaxStableModel
 import pandas as pd
 from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 from spatio_temporal_dataset.temporal_maxima.temporal_maxima import TemporalMaxima
@@ -9,24 +9,18 @@ from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import Ci
 class SimulatedDataset(AbstractDataset):
 
     def __init__(self, temporal_maxima: TemporalMaxima, spatial_coordinates: AbstractSpatialCoordinates,
-                 max_stable_model: MaxStableModel):
+                 max_stable_model: AbstractMaxStableModel):
         super().__init__(temporal_maxima, spatial_coordinates)
         self.max_stable_model = max_stable_model
 
     @classmethod
-    def from_max_stable_sampling(cls, nb_obs: int, max_stable_model:MaxStableModel, spatial_coordinates: AbstractSpatialCoordinates):
+    def from_max_stable_sampling(cls, nb_obs: int, max_stable_model:AbstractMaxStableModel, spatial_coordinates: AbstractSpatialCoordinates):
         maxima = max_stable_model.rmaxstab(nb_obs=nb_obs, coord=spatial_coordinates.coord)
         df_maxima = pd.DataFrame(data=maxima, index=spatial_coordinates.index)
         temporal_maxima = TemporalMaxima(df_maxima=df_maxima)
         return cls(temporal_maxima=temporal_maxima, spatial_coordinates=spatial_coordinates, max_stable_model=max_stable_model)
 
 
-if __name__ == '__main__':
-    coord = CircleCoordinates.from_nb_points(nb_points=5, max_radius=1)
-    max_stable_model = GaussianMSP()
-    dataset = SimulatedDataset.from_max_stable_sampling(nb_obs=50, max_stable_model=max_stable_model, spatial_coordinates=coord)
-    print(dataset.df_dataset)
-
 
 
 

test/R/test_gev_fit.py → test/extreme_estimator/test_gev_fit.py

 import unittest
-
-from R.gev_fit.gev_marginal import frechet_unitary_transformation
-from R.gev_fit.gev_mle_fit import GevMleFit
 import rpy2.robjects as ro
 import numpy as np
 
+from extreme_estimator.R_fit.gev_fit.gev_mle_fit import GevMleFit
+
 
 class TestGevFit(unittest.TestCase):
 

test/extreme_estimator/test_max_stable_fit.py

+import unittest
+
+from extreme_estimator.R_fit.max_stable_fit.abstract_max_stable_model import \
+    AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction
+from extreme_estimator.R_fit.max_stable_fit.max_stable_models import Smith, BrownResnick, Schlather, \
+    Geometric, ExtremalT, ISchlather
+from spatio_temporal_dataset.dataset.simulation_dataset import SimulatedDataset
+from spatio_temporal_dataset.spatial_coordinates.generated_coordinates import CircleCoordinates
+
+
+class TestMaxStableFit(unittest.TestCase):
+    MAX_STABLE_CLASSES = [Smith, BrownResnick, Schlather, Geometric, ExtremalT, ISchlather]
+
+    def setUp(self):
+        super().setUp()
+        self.spatial_coord = CircleCoordinates.from_nb_points(nb_points=5, max_radius=1)
+        self.max_stable_models = []
+        for max_stable_class in self.MAX_STABLE_CLASSES:
+            print(max_stable_class)
+            if issubclass(max_stable_class, AbstractMaxStableModelWithCovarianceFunction):
+                self.max_stable_models.extend([max_stable_class(covariance_function=covariance_function)
+                                               for covariance_function in CovarianceFunction])
+            else:
+                self.max_stable_models.append(max_stable_class())
+
+    def test_sampling_fit_with_models(self, display=False):
+        for max_stable_model in self.max_stable_models:
+            dataset = SimulatedDataset.from_max_stable_sampling(nb_obs=10, max_stable_model=max_stable_model,
+                                                                spatial_coordinates=self.spatial_coord)
+            if display:
+                print(dataset.head())
+            max_stable_model.fitmaxstab(maxima=dataset.maxima, coord=dataset.coord)
+        self.assertTrue(True)
+
+
+if __name__ == '__main__':
+    unittest.main()