From a2ded98d1ebc40a99a6c27d9b61d287808894786 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Wed, 28 Nov 2018 15:00:40 +0100
Subject: [PATCH] add example of python code equivalent or r code (an example
from SpatialExtremes library)
---
example.R | 35 ++++++++++
example.py | 65 +++++++++++++++++++
.../margin_model/abstract_margin_model.py | 8 ++-
.../margin_model/smooth_margin_model.py | 10 +++
extreme_estimator/gev_params.py | 1 +
5 files changed, 117 insertions(+), 2 deletions(-)
create mode 100644 example.R
create mode 100644 example.py
diff --git a/example.R b/example.R
new file mode 100644
index 00000000..0083b6ea
--- /dev/null
+++ b/example.R
@@ -0,0 +1,35 @@
+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'])
+
diff --git a/example.py b/example.py
new file mode 100644
index 00000000..aec25196
--- /dev/null
+++ b/example.py
@@ -0,0 +1,65 @@
+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/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py b/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
index d296b19d..11349070 100644
--- a/extreme_estimator/extreme_models/margin_model/abstract_margin_model.py
+++ b/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
+
+
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 e18ff3a7..4791dcec 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/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):
diff --git a/extreme_estimator/gev_params.py b/extreme_estimator/gev_params.py
index 31cb8247..5ee8656f 100644
--- a/extreme_estimator/gev_params.py
+++ b/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])
+
--
GitLab