add example of python code equivalent or r code (an example from SpatialExtremes library)

example.R

+library(SpatialExtremes)
+# Title     : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 27/11/18
+# define the coordinate of each location
+set.seed(42)
+n.site <- 30
+locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
+colnames(locations) <- c("lon", "lat")
+print(locations)
+##Simulate a max-stable process - with unit Frechet margins
+data <- rmaxstab(40, locations, cov.mod = "whitmat", nugget = 0, range = 3,
+smooth = 0.5)
+##Now define the spatial model for the GEV parameters
+param.loc <- -10 + 2 * locations[,2]
+param.scale <- 5 + 2 * locations[,1] + locations[,2]^2
+param.shape <- 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])
+##Define a model for the GEV margins to be fitted
+##shape ~ 1 stands for the GEV shape parameter is constant
+##over the region
+loc.form <- loc ~ lat
+scale.form <- scale ~ lon + I(lat^2)
+shape.form <- shape ~ 1
+##Fit a max-stable process using the Schlather
+res = fitmaxstab(data, locations, "whitmat", loc.form, scale.form, shape.form)
+print(res['fitted.values'])
+## Model without any spatial structure for the GEV parameters
+## Be careful this could be *REALLY* time consuming
+# res = fitmaxstab(data, locations, "whitmat", fit.margin=TRUE)
+# print(res['fitted.values'])
+

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

extreme_estimator/extreme_models/margin_model/abstract_margin_model.py

@@ -18,8 +18,10 @@ class AbstractMarginModel(AbstractModel):
         self.margin_function_start_fit = None  # type: AbstractMarginFunction
         self.load_margin_functions()
 
-    def load_margin_functions(self, margin_function_class: type = None):
-        assert margin_function_class is not None
+    def load_margin_functions(self):
+        pass
+
+    def default_load_margin_functions(self, margin_function_class):
         self.margin_function_sample = margin_function_class(coordinates=self.coordinates,
                                                             default_params=GevParams.from_dict(self.params_sample))
         self.margin_function_start_fit = margin_function_class(coordinates=self.coordinates,
@@ -69,3 +71,5 @@ class AbstractMarginModel(AbstractModel):
     def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates_values: np.ndarray) \
             -> AbstractMarginFunction:
         pass
+
+

extreme_estimator/extreme_models/margin_model/smooth_margin_model.py

@@ -11,6 +11,8 @@ from extreme_estimator.gev_params import GevParams
 class LinearMarginModel(AbstractMarginModel):
 
     def load_margin_functions(self, gev_param_name_to_linear_dims=None):
+        assert gev_param_name_to_linear_dims is not None, 'LinearMarginModel cannot be used for sampling/fitting \n' \
+                                                          'load_margin_functions needs to be implemented in child class'
         # Load sample coef
         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)
@@ -49,6 +51,14 @@ class LinearMarginModel(AbstractMarginModel):
                                   coordinates_values: np.ndarray) -> AbstractMarginFunction:
         return self.margin_function_start_fit
 
+    @classmethod
+    def from_coef_list(cls, coordinates, gev_param_name_to_coef_list):
+        params = {}
+        for gev_param_name in GevParams.GEV_PARAM_NAMES:
+            for dim, coef in enumerate(gev_param_name_to_coef_list[gev_param_name]):
+                params[(gev_param_name, dim)] = coef
+        return cls(coordinates, params_sample=params, params_start_fit=params)
+
 
 class ConstantMarginModel(LinearMarginModel):
 

extreme_estimator/gev_params.py

@@ -27,3 +27,4 @@ class GevParams(object):
     def to_array(self) -> np.ndarray:
         gev_param_name_to_value = self.to_dict()
         return np.array([gev_param_name_to_value[gev_param_name] for gev_param_name in self.GEV_PARAM_NAMES])
+