[EXTREME FIT][MAX STABLE] fix start fit params issue for fitmaxstab. add test

extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py

+import rpy2
 from enum import Enum
 
 from rpy2.robjects import ListVector
@@ -15,19 +16,20 @@ class AbstractMaxStableModel(object):
         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
+    def fitmaxstab(self, maxima: np.ndarray, coord: np.ndarray, fit_marge=False):
+        #  Specify the fit params
+        fit_params = {
+            'fit.marge': fit_marge,
+            'start': self.r.list(**self.params_start_fit),
+        }
+        # Run the fitmaxstab in R
+        # todo: find how to specify the optim function to use
+        try:
+            res = self.r.fitmaxstab(np.transpose(maxima), coord, **self.cov_mod_param, **fit_params)  # type: ListVector
+        except rpy2.rinterface.RRuntimeError as error:
+            raise Exception('Some R exception have been launched at RunTime: {}'.format(error.__repr__()))
+        # 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')
         fitted_values = {key: fitted_values.rx2(key)[0] for key in fitted_values.names}
         return fitted_values
@@ -36,7 +38,7 @@ class AbstractMaxStableModel(object):
         """
         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))
+        maxima = np.array(self.r.rmaxstab(nb_obs, coord, *list(self.cov_mod_param.values()), **self.params_sample))
         return np.transpose(maxima)
 
     @property
@@ -53,6 +55,7 @@ class AbstractMaxStableModel(object):
 
     @staticmethod
     def merge_params(default_params, input_params):
+        assert default_params is not None, 'some default_params need to be specified'
         merged_params = default_params.copy()
         if input_params is not None:
             assert isinstance(default_params, dict) and isinstance(input_params, dict)

extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R

@@ -10,22 +10,26 @@ if (call_main) {
     n.site = 2
     coord <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
 
-    print(coord)
     #  Generate the data
-    # data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220)
+    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, "whitmat", nugget = 0.0, range = 3, smooth = 0.5)
-    data <- rmaxstab(n.obs, coord, "cauchy", )
     #  Fit back the data
-    print(data)
+    # print(data)
     # res = fitmaxstab(data, coord, "gauss", fit.marge=FALSE, )
     # 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))
-    print(names(res))
+    # res = fitmaxstab(data, coord, "whitmat", start=)
+    namedlist = list(cov11 = 1.0, cov12 = 1.2, cov22 = 2.2)
+    res = fitmaxstab(data=data, coord=coord, cov.mod="gauss", start=namedlist)
+    # res = fitmaxstab(data, coord, "whitmat", par())
+    # print(res)
+    # print(class(res))
+    # print(names(res))
+    for (name in names(res)){
+        print(name)
+        print(res[name])
+    }
     print(res['fitted.values'])
-    print(res['convergence'])
+    # print(res['convergence'])
 
 }

extreme_estimator/R_fit/max_stable_fit/max_stable_models.py

@@ -9,12 +9,12 @@ 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)
         self.cov_mod = 'gauss'
-        self.default_params_start_fit = {}
-        self.default_params_sample = {
-            'cov11': 100,
-            'cov12': 25,
-            'cov22': 220,
+        self.default_params_start_fit = {
+            'cov11': 1,
+            'cov12': 0,
+            'cov22': 1
         }
+        self.default_params_sample = self.default_params_start_fit.copy()
 
 
 class BrownResnick(AbstractMaxStableModel):
@@ -22,7 +22,10 @@ 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_start_fit = {
+            'range': 3,
+            'smooth': 0.5,
+        }
         self.default_params_sample = {
             'range': 3,
             'smooth': 0.5,
@@ -35,6 +38,7 @@ class Schlather(AbstractMaxStableModelWithCovarianceFunction):
         super().__init__(params_start_fit, params_sample, covariance_function)
         self.cov_mod = self.covariance_function.name
         self.default_params_sample.update({})
+        self.default_params_start_fit = self.default_params_sample.copy()
 
 
 class Geometric(AbstractMaxStableModelWithCovarianceFunction):
@@ -43,6 +47,7 @@ class Geometric(AbstractMaxStableModelWithCovarianceFunction):
         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})
+        self.default_params_start_fit = self.default_params_sample.copy()
 
 
 class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
@@ -51,6 +56,7 @@ class ExtremalT(AbstractMaxStableModelWithCovarianceFunction):
         super().__init__(params_start_fit, params_sample, covariance_function)
         self.cov_mod = 't' + self.covariance_function.name
         self.default_params_sample .update({'DoF': 2})
+        self.default_params_start_fit = self.default_params_sample.copy()
 
 
 class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
@@ -59,3 +65,4 @@ class ISchlather(AbstractMaxStableModelWithCovarianceFunction):
         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})
+        self.default_params_start_fit = self.default_params_sample.copy()

test/extreme_estimator/test_max_stable_fit.py

@@ -16,7 +16,6 @@ class TestMaxStableFit(unittest.TestCase):
         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])
@@ -27,9 +26,11 @@ class TestMaxStableFit(unittest.TestCase):
         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)
+            fitted_values = max_stable_model.fitmaxstab(maxima=dataset.maxima, coord=dataset.coord)
             if display:
-                print(dataset.head())
-            max_stable_model.fitmaxstab(maxima=dataset.maxima, coord=dataset.coord)
+                print(type(max_stable_model))
+                print(dataset.df_dataset.head())
+                print(fitted_values)
         self.assertTrue(True)