From 219ab0c431a00e38e86e467773deb403930a1e9c Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 28 Nov 2018 18:09:09 +0100
Subject: [PATCH] [FIX] fix fitmaxstable, add unitary tests
---
example.py | 65 -----------
.../regression_margin/regression_margin.py | 36 +++---
extreme_estimator/estimator/full_estimator.py | 2 +-
.../margin_model/smooth_margin_model.py | 8 ++
.../abstract_max_stable_model.py | 23 ++--
test/test_unitary/__init__.py | 0
test/test_unitary/test_fitmaxstab/__init__.py | 0
.../test_fitmaxstab_with_margin.py | 71 ++++++++++++
.../test_fitmaxstab_without_margin.py | 35 ++++++
test/test_unitary/test_rmaxstab/__init__.py | 0
.../test_rmaxstab_with_margin.py | 106 ++++++++++++++++++
.../test_rmaxstab_without_margin.py | 47 ++++++++
test/test_unitary/test_unitary_abstract.py | 37 ++++++
13 files changed, 341 insertions(+), 89 deletions(-)
delete mode 100644 example.py
create mode 100644 test/test_unitary/__init__.py
create mode 100644 test/test_unitary/test_fitmaxstab/__init__.py
create mode 100644 test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
create mode 100644 test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py
create mode 100644 test/test_unitary/test_rmaxstab/__init__.py
create mode 100644 test/test_unitary/test_rmaxstab/test_rmaxstab_with_margin.py
create mode 100644 test/test_unitary/test_rmaxstab/test_rmaxstab_without_margin.py
create mode 100644 test/test_unitary/test_unitary_abstract.py
diff --git a/example.py b/example.py
deleted file mode 100644
index aec25196..00000000
--- a/example.py
+++ /dev/null
@@ -1,65 +0,0 @@
-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)
-
diff --git a/experiment/regression_margin/regression_margin.py b/experiment/regression_margin/regression_margin.py
index a46fd202..ee49fa93 100644
--- a/experiment/regression_margin/regression_margin.py
+++ b/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)
diff --git a/extreme_estimator/estimator/full_estimator.py b/extreme_estimator/estimator/full_estimator.py
index 35d9ca96..20eda73f 100644
--- a/extreme_estimator/estimator/full_estimator.py
+++ b/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,
diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
index 4791dcec..4fb13aa5 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/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):
diff --git a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py b/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
index b051b687..f9e1235c 100644
--- a/extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py
+++ b/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')
diff --git a/test/test_unitary/__init__.py b/test/test_unitary/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/test/test_unitary/test_fitmaxstab/__init__.py b/test/test_unitary/test_fitmaxstab/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
new file mode 100644
index 00000000..c0ca47c8
--- /dev/null
+++ b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py
@@ -0,0 +1,71 @@
+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()
diff --git a/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py
new file mode 100644
index 00000000..9256c0ae
--- /dev/null
+++ b/test/test_unitary/test_fitmaxstab/test_fitmaxstab_without_margin.py
@@ -0,0 +1,35 @@
+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()
diff --git a/test/test_unitary/test_rmaxstab/__init__.py b/test/test_unitary/test_rmaxstab/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/test/test_unitary/test_rmaxstab/test_rmaxstab_with_margin.py b/test/test_unitary/test_rmaxstab/test_rmaxstab_with_margin.py
new file mode 100644
index 00000000..fac7b393
--- /dev/null
+++ b/test/test_unitary/test_rmaxstab/test_rmaxstab_with_margin.py
@@ -0,0 +1,106 @@
+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()
diff --git a/test/test_unitary/test_rmaxstab/test_rmaxstab_without_margin.py b/test/test_unitary/test_rmaxstab/test_rmaxstab_without_margin.py
new file mode 100644
index 00000000..c97a0726
--- /dev/null
+++ b/test/test_unitary/test_rmaxstab/test_rmaxstab_without_margin.py
@@ -0,0 +1,47 @@
+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()
diff --git a/test/test_unitary/test_unitary_abstract.py b/test/test_unitary/test_unitary_abstract.py
new file mode 100644
index 00000000..acf305e2
--- /dev/null
+++ b/test/test_unitary/test_unitary_abstract.py
@@ -0,0 +1,37 @@
+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))
+
+
--
GitLab