diff --git a/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py b/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py
index b6f05736b9d18fc9734f46fc1b9286b2d3a0d3dd..5bfe53f38905c32542fccdc55d8597a6d246d4fe 100644
--- a/extreme_estimator/R_fit/gev_fit/gev_mle_fit.py
+++ b/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
diff --git a/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
new file mode 100644
index 0000000000000000000000000000000000000000..75a4c3cccb823ee14ffb1472277d3226fd6dfd65
--- /dev/null
+++ b/extreme_estimator/R_fit/max_stable_fit/abstract_max_stable_model.py
@@ -0,0 +1,87 @@
+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
diff --git a/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R b/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
index b51112f7a53e36f9236df8eedd0e0057a7cb0513..287a186b83a731c2ff6b5aea44975952bc1b9dd9 100644
--- a/extreme_estimator/R_fit/max_stable_fit/max_stable_fit.R
+++ b/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))
diff --git a/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py b/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py
index 3aaaa7dd0e66fb2b9cf47a8aae23922ea663d0b0..2c77573f200a818e725e7c56147c0cea4d04e3d9 100644
--- a/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py
+++ b/extreme_estimator/R_fit/max_stable_fit/max_stable_models.py
@@ -1,72 +1,10 @@
-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})
diff --git a/extreme_estimator/estimator/msp_estimator.py b/extreme_estimator/estimator/msp_estimator.py
index 03c2f502a5afb33a7f4bf953e00c2a2f0ea8c84c..bb467d09985dc2d850480cc0e9990dcef47d7686 100644
--- a/extreme_estimator/estimator/msp_estimator.py
+++ b/extreme_estimator/estimator/msp_estimator.py
@@ -1,12 +1,12 @@
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
diff --git a/extreme_estimator/robustness_plot/abstract_robustness.py b/extreme_estimator/robustness_plot/abstract_robustness.py
index b40f9bd2dc28932a0e8a1933305d0d0b8d485de9..4925440748052f919a3268c57ecbc822bd3f090b 100644
--- a/extreme_estimator/robustness_plot/abstract_robustness.py
+++ b/extreme_estimator/robustness_plot/abstract_robustness.py
@@ -1,7 +1,7 @@
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
diff --git a/extreme_estimator/robustness_plot/spatial_robustness.py b/extreme_estimator/robustness_plot/spatial_robustness.py
index 3d872cd2a0ac2522038d6ac2734853c34973d665..937a66b7c51d500a495c0614ea9c09519a60b516 100644
--- a/extreme_estimator/robustness_plot/spatial_robustness.py
+++ b/extreme_estimator/robustness_plot/spatial_robustness.py
@@ -1,4 +1,3 @@
-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
diff --git a/spatio_temporal_dataset/dataset/simulation_dataset.py b/spatio_temporal_dataset/dataset/simulation_dataset.py
index bc6a2860532661a69ccda797f76e3ec36882bcc0..7312acc945598a29658bef1d76b0fb421b5cd608 100644
--- a/spatio_temporal_dataset/dataset/simulation_dataset.py
+++ b/spatio_temporal_dataset/dataset/simulation_dataset.py
@@ -1,4 +1,4 @@
-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)
-
diff --git a/test/R/test_gev_fit.py b/test/extreme_estimator/test_gev_fit.py
similarity index 90%
rename from test/R/test_gev_fit.py
rename to test/extreme_estimator/test_gev_fit.py
index c2f4b9551ee74ac39a6ecb41a9ade20be14716a2..17b963658bfe5e560ae0b5ddfd20b50599e335c9 100644
--- a/test/R/test_gev_fit.py
+++ b/test/extreme_estimator/test_gev_fit.py
@@ -1,10 +1,9 @@
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):
diff --git a/test/extreme_estimator/test_max_stable_fit.py b/test/extreme_estimator/test_max_stable_fit.py
new file mode 100644
index 0000000000000000000000000000000000000000..cca832bfa3e1454d063df35c13550882867cd5ca
--- /dev/null
+++ b/test/extreme_estimator/test_max_stable_fit.py
@@ -0,0 +1,37 @@
+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()