From 2c05724bd6b196795706a9860925c8284bcdcc96 Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Thu, 1 Apr 2021 17:43:28 +0200
Subject: [PATCH] [projection snowfall] add evgam call for
abstract_temporal_linear_margin_model.py. account for log_scale in the
margin_function. rename get_coord_df to
get_r_dataframe_from_python_dataframe. create evgam_fixed.R for debugging
---
extreme_fit/distribution/gev/evgam_example.R | 18 ++++
extreme_fit/distribution/gev/evgam_fixed.R | 42 +++++++++
extreme_fit/distribution/gev/main_fevd_mle.R | 12 +--
extreme_fit/estimator/utils.py | 7 +-
.../abstract_margin_function.py | 3 +-
.../independent_margin_function.py | 6 +-
.../margin_function/linear_margin_function.py | 6 +-
.../parametric_margin_function.py | 10 +-
.../polynomial_margin_function.py | 8 +-
.../abstract_temporal_linear_margin_model.py | 28 +++++-
extreme_fit/model/margin_model/utils.py | 1 +
.../quantile_regression_model.py | 4 +-
.../abstract_result_from_model_fit.py | 5 +
.../result_from_extremes/result_from_evgam.py | 42 +++++++++
extreme_fit/model/utils.py | 17 ++--
.../test_gev_temporal_evgam.py | 93 +++++++++++++++++++
16 files changed, 269 insertions(+), 33 deletions(-)
create mode 100644 extreme_fit/distribution/gev/evgam_example.R
create mode 100644 extreme_fit/distribution/gev/evgam_fixed.R
create mode 100644 extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py
create mode 100644 test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py
diff --git a/extreme_fit/distribution/gev/evgam_example.R b/extreme_fit/distribution/gev/evgam_example.R
new file mode 100644
index 00000000..0aadfd88
--- /dev/null
+++ b/extreme_fit/distribution/gev/evgam_example.R
@@ -0,0 +1,18 @@
+# Title : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 30/03/2021
+source('evgam_fixed.R')
+library(evgam)
+data(COprcp)
+COprcp$year <- format(COprcp$date, "%Y")
+COprcp_gev <- aggregate(prcp ~ year + meta_row, COprcp, max)
+COprcp_gev <- cbind(COprcp_gev, COprcp_meta[COprcp_gev$meta_row,])
+# print(COprcp_gev)
+print('before call')
+fmla_gev2 <- list(prcp ~ s(elev, bs="cr"), ~ s(lon, fx=FALSE, k=20), ~ 1)
+m_gev2 <- evgam_fixed(fmla_gev2, data=COprcp_gev, family="gev")
+# m_gev2 <- evgam(fmla_gev2, data=COprcp_gev, family="gev")
+# summary(m_gev2)
+# print(attributes(m_gev2))
+print('good finish')
diff --git a/extreme_fit/distribution/gev/evgam_fixed.R b/extreme_fit/distribution/gev/evgam_fixed.R
new file mode 100644
index 00000000..d978021a
--- /dev/null
+++ b/extreme_fit/distribution/gev/evgam_fixed.R
@@ -0,0 +1,42 @@
+# Title : TODO
+# Objective : TODO
+# Created by: erwan
+# Created on: 01/04/2021
+library(evgam)
+attach(loadNamespace("evgam"), name = "evgam_all")
+evgam_fixed <- function(formula, data, family = "gev", correctV = TRUE, rho0 = 0, inits = NULL, outer = "bfgs", control = NULL, removeData = FALSE, trace = 0, knots = NULL, maxdata = 1e+20, maxspline = 1e+20, compact = FALSE, ald.args = list(), exi.args = list(), pp.args = list(), sandwich.args = list()) {
+ family.info <- .setup.family(family, pp.args)
+ if (is.null(family.info$lik.fns$d340)) outer <- "fd"
+ formula <- .setup.formulae(formula, family.info$npar, family.info$npar2, data, trace)
+ response.name <- attr(formula, "response.name")
+ temp.data <- .setup.data(data, response.name, formula, family, family.info$nms, removeData, exi.args, ald.args, pp.args, knots, maxdata, maxspline, compact, sandwich.args, tolower(outer), trace)
+ data <- temp.data$data
+ beta <- .setup.inner.inits(inits, temp.data$lik.data, family.info$lik.fns, family.info$npar, family)
+ lik.data <- .sandwich(temp.data$lik.data, beta)
+ if (trace > 0 & lik.data$adjust > 0) cat(paste("\n Sandwich correct lambda =", signif(lik.data$k, 3), "\n"))
+ smooths <- length(temp.data$gotsmooth) > 0
+ if (smooths) {
+ S.data <- .joinSmooth(temp.data$gams)
+ nsp <- length(attr(S.data, "Sl"))
+ if (is.null(rho0)) {
+ diagSl <- sapply(attr(S.data, "Sl"), diag)
+ rho0 <- apply(diagSl, 2, function(y) uniroot(.guess, c(-100, 100), d = attr(beta, "diagH"), s = y)$root)
+ } else {
+ if (length(rho0) == 1) rho0 <- rep(rho0, nsp)
+ }
+ lik.data$S <- .makeS(S.data, exp(rho0))
+ fit.reml <- .outer(rho0, beta, family.info$lik.fns, lik.data, S.data, control, correctV, lik.data$outer, trace)
+ sp <- exp(fit.reml$par)
+ lik.data$S <- .makeS(S.data, sp)
+ } else {
+ S.data <- NULL
+ fit.reml <- .outer.nosmooth(beta, family.info$lik.fns, lik.data, control, trace)
+ }
+ VpVc <- .VpVc(fit.reml, family.info$lik.fns, lik.data, S.data, correctV = correctV, sandwich = temp.data$sandwich, smooths = smooths, trace = trace)
+ edf <- .edf(fit.reml$beta, family.info$lik.fns, lik.data, VpVc, temp.data$sandwich)
+ names(temp.data$gams) <- family.info$nms
+ gams <- .swap(fit.reml, temp.data$gams, lik.data, VpVc, temp.data$gotsmooth, edf, smooths)
+ gams <- .finalise(gams, data, family.info$lik.fns, lik.data, S.data, fit.reml, VpVc, family, temp.data$gotsmooth, formula, response.name, removeData, edf)
+ return(gams)
+}
+
diff --git a/extreme_fit/distribution/gev/main_fevd_mle.R b/extreme_fit/distribution/gev/main_fevd_mle.R
index 972c9a64..1a97d9ad 100644
--- a/extreme_fit/distribution/gev/main_fevd_mle.R
+++ b/extreme_fit/distribution/gev/main_fevd_mle.R
@@ -17,16 +17,16 @@ x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
# N <- 50
# loc = 0; scale = 1; shape <- 0.1
# x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
-coord <- matrix(ncol=3, nrow = N)
+coord <- matrix(ncol=1, nrow = N)
coord[,1]=seq(0,N-1,1)
-coord[,2]=c(rep(0, N/2), rep(1, N/2))
-coord[,3]=c(rep(1, N/2), rep(0, N/2))
-colnames(coord) = c("T", "G1", "G2")
+# coord[,2]=c(rep(0, N/2), rep(1, N/2))
+# coord[,3]=c(rep(1, N/2), rep(0, N/2))
+colnames(coord) = c("T")
print(coord)
coord = data.frame(coord, stringsAsFactors = TRUE)
# res = fevd_fixed(x_gev, data=coord, method='MLE', verbose=TRUE, use.phi=FALSE)
-# res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, shape.fun= ~1 + (G1-G2) + I(G2), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~T, scale.fun= ~T, method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
+# res = fevd_fixed(x_gev, data=coord, shape.fun= ~1 + (G1-G2) + I(G2), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
print(res)
# Some display for the results
diff --git a/extreme_fit/estimator/utils.py b/extreme_fit/estimator/utils.py
index b42a1b64..2d8699c5 100644
--- a/extreme_fit/estimator/utils.py
+++ b/extreme_fit/estimator/utils.py
@@ -5,14 +5,17 @@ from extreme_fit.model.margin_model.linear_margin_model.linear_margin_model impo
from extreme_fit.function.margin_function.linear_margin_function import LinearMarginFunction
-def load_margin_function(estimator: AbstractEstimator, margin_model: LinearMarginModel, coef_dict=None):
+def load_margin_function(estimator: AbstractEstimator, margin_model: LinearMarginModel, coef_dict=None, log_scale=None):
if coef_dict is None:
coef_dict = estimator.result_from_model_fit.margin_coef_ordered_dict
+ if log_scale is None:
+ log_scale = estimator.result_from_model_fit.log_scale
margin_function_class = type(margin_model.margin_function)
return margin_function_class.from_coef_dict(coordinates=estimator.dataset.coordinates,
param_name_to_dims=margin_model.param_name_to_list_for_result,
coef_dict=coef_dict,
- starting_point=margin_model.starting_point)
+ starting_point=margin_model.starting_point,
+ log_scale=log_scale)
diff --git a/extreme_fit/function/margin_function/abstract_margin_function.py b/extreme_fit/function/margin_function/abstract_margin_function.py
index a271da29..8f3a7b14 100644
--- a/extreme_fit/function/margin_function/abstract_margin_function.py
+++ b/extreme_fit/function/margin_function/abstract_margin_function.py
@@ -20,8 +20,9 @@ class AbstractMarginFunction(AbstractFunction):
VISUALIZATION_RESOLUTION = 100
VISUALIZATION_TEMPORAL_STEPS = 2
- def __init__(self, coordinates: AbstractCoordinates, params_class: type = GevParams):
+ def __init__(self, coordinates: AbstractCoordinates, params_class: type = GevParams, log_scale=None):
super().__init__(coordinates)
+ self.log_scale = log_scale
self.params_class = params_class
self.mask_2D = None
diff --git a/extreme_fit/function/margin_function/independent_margin_function.py b/extreme_fit/function/margin_function/independent_margin_function.py
index 79db8d38..055c2848 100644
--- a/extreme_fit/function/margin_function/independent_margin_function.py
+++ b/extreme_fit/function/margin_function/independent_margin_function.py
@@ -14,9 +14,9 @@ class IndependentMarginFunction(AbstractMarginFunction):
IndependentMarginFunction: each parameter of the GEV are modeled independently
"""
- def __init__(self, coordinates: AbstractCoordinates, params_class: type = GevParams):
+ def __init__(self, coordinates: AbstractCoordinates, params_class: type = GevParams, log_scale=None):
"""Attribute 'param_name_to_param_function' maps each GEV parameter to its corresponding function"""
- super().__init__(coordinates, params_class)
+ super().__init__(coordinates, params_class, log_scale)
self.param_name_to_param_function = None # type: Union[None, Dict[str, AbstractParamFunction]]
def get_params(self, coordinate: np.ndarray, is_transformed: bool = True) -> GevParams:
@@ -28,6 +28,8 @@ class IndependentMarginFunction(AbstractMarginFunction):
transformed_coordinate = coordinate if is_transformed else self.transform(coordinate)
params = {param_name: param_function.get_param_value(transformed_coordinate)
for param_name, param_function in self.param_name_to_param_function.items()}
+ if self.log_scale:
+ params[GevParams.SCALE] = np.exp(params[GevParams.SCALE])
return self.params_class.from_dict(params)
def get_first_derivative_param(self, coordinate: np.ndarray, is_transformed: bool, dim: int = 0):
diff --git a/extreme_fit/function/margin_function/linear_margin_function.py b/extreme_fit/function/margin_function/linear_margin_function.py
index b989c1c6..d387341f 100644
--- a/extreme_fit/function/margin_function/linear_margin_function.py
+++ b/extreme_fit/function/margin_function/linear_margin_function.py
@@ -28,11 +28,9 @@ class LinearMarginFunction(ParametricMarginFunction):
COEF_CLASS = LinearCoef
- def __init__(self, coordinates: AbstractCoordinates, param_name_to_dims: Dict[str, List[int]],
- param_name_to_coef: Dict[str, AbstractCoef], starting_point: Union[None, int] = None,
- params_class: type = GevParams):
+ def __init__(self, *args, **kwargs):
self.param_name_to_coef = None # type: Union[None, Dict[str, LinearCoef]]
- super().__init__(coordinates, param_name_to_dims, param_name_to_coef, starting_point, params_class)
+ super().__init__(*args, **kwargs)
def load_specific_param_function(self, param_name) -> AbstractParamFunction:
return LinearParamFunction(dims=self.param_name_to_dims[param_name],
diff --git a/extreme_fit/function/margin_function/parametric_margin_function.py b/extreme_fit/function/margin_function/parametric_margin_function.py
index d08609d8..41d8f882 100644
--- a/extreme_fit/function/margin_function/parametric_margin_function.py
+++ b/extreme_fit/function/margin_function/parametric_margin_function.py
@@ -32,10 +32,11 @@ class ParametricMarginFunction(IndependentMarginFunction):
def __init__(self, coordinates: AbstractCoordinates, param_name_to_dims: Dict[str, List[int]],
param_name_to_coef: Dict[str, AbstractCoef], starting_point: Union[None, int] = None,
- params_class: type = GevParams):
+ params_class: type = GevParams,
+ log_scale=None):
# Starting point for the trend is the same for all the parameters
self.starting_point = starting_point
- super().__init__(coordinates, params_class)
+ super().__init__(coordinates, params_class, log_scale=log_scale)
self.param_name_to_dims = param_name_to_dims # type: Dict[str, List[int]]
# Check the dimension are well-defined with respect to the coordinates
@@ -89,7 +90,8 @@ class ParametricMarginFunction(IndependentMarginFunction):
@classmethod
def from_coef_dict(cls, coordinates: AbstractCoordinates, param_name_to_dims: Dict[str, List[int]],
- coef_dict: Dict[str, float], starting_point: Union[None, int] = None):
+ coef_dict: Dict[str, float], starting_point: Union[None, int] = None,
+ log_scale=None):
assert cls.COEF_CLASS is not None, 'a COEF_CLASS class attributes needs to be defined'
param_name_to_coef = {}
for param_name in GevParams.PARAM_NAMES:
@@ -97,7 +99,7 @@ class ParametricMarginFunction(IndependentMarginFunction):
coef = cls.COEF_CLASS.from_coef_dict(coef_dict=coef_dict, param_name=param_name, dims=dims,
coordinates=coordinates)
param_name_to_coef[param_name] = coef
- return cls(coordinates, param_name_to_dims, param_name_to_coef, starting_point)
+ return cls(coordinates, param_name_to_dims, param_name_to_coef, starting_point=starting_point, log_scale=log_scale)
@property
def form_dict(self) -> Dict[str, str]:
diff --git a/extreme_fit/function/margin_function/polynomial_margin_function.py b/extreme_fit/function/margin_function/polynomial_margin_function.py
index 06bbeac5..e7814ed7 100644
--- a/extreme_fit/function/margin_function/polynomial_margin_function.py
+++ b/extreme_fit/function/margin_function/polynomial_margin_function.py
@@ -14,13 +14,13 @@ class PolynomialMarginFunction(LinearMarginFunction):
def __init__(self, coordinates: AbstractCoordinates, param_name_to_dim_and_max_degree: Dict[str, List[Tuple[int, int]]],
param_name_to_coef: Dict[str, PolynomialAllCoef], starting_point: Union[None, int] = None,
- params_class: type = GevParams):
+ params_class: type = GevParams, log_scale=None):
param_name_to_dims = {}
for param_name in param_name_to_dim_and_max_degree.keys():
dims = [c[0] for c in param_name_to_dim_and_max_degree[param_name]]
param_name_to_dims[param_name] = dims
self.param_name_to_dim_and_max_degree = param_name_to_dim_and_max_degree
- super().__init__(coordinates, param_name_to_dims, param_name_to_coef, starting_point, params_class)
+ super().__init__(coordinates, param_name_to_dims, param_name_to_coef, starting_point, params_class, log_scale)
COEF_CLASS = PolynomialAllCoef
@@ -33,7 +33,7 @@ class PolynomialMarginFunction(LinearMarginFunction):
@classmethod
def from_coef_dict(cls, coordinates: AbstractCoordinates, param_name_to_dims: Dict[str, List[Tuple[int, int]]],
- coef_dict: Dict[str, float], starting_point: Union[None, int] = None):
+ coef_dict: Dict[str, float], starting_point: Union[None, int] = None, log_scale=None):
param_name_to_dim_and_max_degree = param_name_to_dims
assert cls.COEF_CLASS is not None, 'a COEF_CLASS class attributes needs to be defined'
param_name_to_coef = {}
@@ -43,5 +43,5 @@ class PolynomialMarginFunction(LinearMarginFunction):
dims=dims,
coordinates=coordinates)
param_name_to_coef[param_name] = coef
- return cls(coordinates, param_name_to_dim_and_max_degree, param_name_to_coef, starting_point)
+ return cls(coordinates, param_name_to_dim_and_max_degree, param_name_to_coef, starting_point, log_scale=log_scale)
diff --git a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
index 090cfb25..22a74298 100644
--- a/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
+++ b/extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py
@@ -1,5 +1,6 @@
import numpy as np
import pandas as pd
+from rpy2 import robjects
from extreme_fit.distribution.exp_params import ExpParams
from extreme_fit.distribution.gev.gev_params import GevParams
@@ -9,9 +10,10 @@ from extreme_fit.model.margin_model.utils import MarginFitMethod, fitmethod_to_s
from extreme_fit.model.result_from_model_fit.abstract_result_from_model_fit import AbstractResultFromModelFit
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_bayesian_extremes import \
AbstractResultFromExtremes, ResultFromBayesianExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_evgam import ResultFromEvgam
from extreme_fit.model.result_from_model_fit.result_from_extremes.result_from_mle_extremes import ResultFromMleExtremes
from extreme_fit.model.result_from_model_fit.result_from_ismev import ResultFromIsmev
-from extreme_fit.model.utils import r, ro, get_null, get_margin_formula_extremes, get_coord_df
+from extreme_fit.model.utils import r, ro, get_null, get_margin_formula_extremes, get_r_dataframe_from_python_dataframe
from extreme_fit.model.utils import safe_run_r_estimator
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -46,6 +48,8 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
return self.ismev_gev_fit(x, df_coordinates_temp)
elif self.fit_method in FEVD_MARGIN_FIT_METHODS:
return self.extremes_fevd_fit(x, df_coordinates_temp, df_coordinates_spat)
+ elif self.fit_method is MarginFitMethod.evgam:
+ return self.extremes_evgam_fit(x, df_coordinates_temp, df_coordinates_spat)
else:
raise NotImplementedError
elif self.params_class is ExpParams:
@@ -62,6 +66,26 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
sigl=self.sigl, shl=self.shl)
return ResultFromIsmev(res, self.param_name_to_list_for_result)
+ # Gev fit with evgam
+
+ def extremes_evgam_fit(self, x, df_coordinates_temp, df_coordinates_spat) -> AbstractResultFromExtremes:
+ margin_formula = get_margin_formula_extremes(self.margin_function.form_dict, transformed_as_formula=False)
+ maxima_column_name = 'Maxima'
+ formula_list = [maxima_column_name + " " + v if i == 0 else v for i, v in enumerate(margin_formula.values())]
+ formula = r.list(*[robjects.Formula(f) for f in formula_list])
+ df = pd.DataFrame({maxima_column_name: np.array(x)})
+ df = pd.concat([df, df_coordinates_spat, df_coordinates_temp], axis=1)
+ data = get_r_dataframe_from_python_dataframe(df)
+ res = safe_run_r_estimator(function=r('evgam'),
+ formula=formula,
+ data=data,
+ family=self.type_for_mle.lower(),
+ maxdata=1e10,
+ )
+ return ResultFromEvgam(res, self.param_name_to_list_for_result,
+ self.coordinates.dim_to_coordinate,
+ type_for_mle=self.type_for_mle)
+
# Gev fit with extRemes package
def extremes_fevd_fit(self, x, df_coordinates_temp, df_coordinates_spat) -> AbstractResultFromExtremes:
@@ -118,7 +142,7 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
df = df_coordinates_temp
else:
df = pd.concat([df_coordinates_spat, df_coordinates_temp], axis=1)
- y = get_coord_df(df)
+ y = get_r_dataframe_from_python_dataframe(df)
# Disable the use of log sigma parametrization
r_type_argument_kwargs = {'use.phi': False,
diff --git a/extreme_fit/model/margin_model/utils.py b/extreme_fit/model/margin_model/utils.py
index e713e18c..943ebbc0 100644
--- a/extreme_fit/model/margin_model/utils.py
+++ b/extreme_fit/model/margin_model/utils.py
@@ -8,6 +8,7 @@ class MarginFitMethod(Enum):
extremes_fevd_gmle = 3
extremes_fevd_l_moments = 4
spatial_extremes_mle = 5
+ evgam = 6
FEVD_MARGIN_FIT_METHOD_TO_ARGUMENT_STR = {
diff --git a/extreme_fit/model/quantile_model/quantile_regression_model.py b/extreme_fit/model/quantile_model/quantile_regression_model.py
index 7bf51b02..3a84f1de 100644
--- a/extreme_fit/model/quantile_model/quantile_regression_model.py
+++ b/extreme_fit/model/quantile_model/quantile_regression_model.py
@@ -3,7 +3,7 @@ from rpy2 import robjects
from extreme_fit.model.abstract_model import AbstractModel
from extreme_fit.model.result_from_model_fit.result_from_quantilreg import ResultFromQuantreg
-from extreme_fit.model.utils import r, safe_run_r_estimator, get_coord_df
+from extreme_fit.model.utils import r, safe_run_r_estimator, get_r_dataframe_from_python_dataframe
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
@@ -16,7 +16,7 @@ class AbstractQuantileRegressionModel(AbstractModel):
@property
def data(self):
- return get_coord_df(self.dataset.df_dataset)
+ return get_r_dataframe_from_python_dataframe(self.dataset.df_dataset)
@property
def first_column_of_observation(self):
diff --git a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
index 85dfbea3..9f5bd468 100644
--- a/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
+++ b/extreme_fit/model/result_from_model_fit/abstract_result_from_model_fit.py
@@ -32,6 +32,11 @@ class AbstractResultFromModelFit(object):
def all_parameters(self):
raise NotImplementedError
+ @property
+ def log_scale(self):
+ # todo: refactor, put raise NotImplementError, and fix the subfunctions for the other Result objects
+ return None
+
@property
def margin_coef_ordered_dict(self):
raise NotImplementedError
diff --git a/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py
new file mode 100644
index 00000000..8e674a52
--- /dev/null
+++ b/extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py
@@ -0,0 +1,42 @@
+import numpy as np
+from rpy2 import robjects
+
+from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_result_from_extremes import \
+ AbstractResultFromExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+ ci_method_to_method_name
+from extreme_fit.model.result_from_model_fit.utils import get_margin_coef_ordered_dict
+from extreme_fit.model.utils import r
+
+
+class ResultFromEvgam(AbstractResultFromExtremes):
+
+ def __init__(self, result_from_fit: robjects.ListVector, param_name_to_dim=None,
+ dim_to_coordinate=None,
+ type_for_mle="GEV") -> None:
+ super().__init__(result_from_fit, param_name_to_dim, dim_to_coordinate)
+ self.type_for_mle = type_for_mle
+
+ @property
+ def aic(self):
+ """Compute the aic from the list of parameters in the results,
+ find a way to comptue it directly from the result to compare"""
+ location = self.name_to_value['location']
+ a = np.array(location)
+ print(a)
+ print(len(a))
+ print('here2')
+ # 'location', 'logscale', 'shape'
+ raise NotImplementedError
+
+ @property
+ def log_scale(self):
+ return True
+
+ @property
+ def margin_coef_ordered_dict(self):
+ # print(self.name_to_value.keys())
+ # raise NotImplementedError
+ values = np.array(self.name_to_value['coefficients'])
+ return get_margin_coef_ordered_dict(self.param_name_to_dim, values, self.type_for_mle,
+ dim_to_coordinate_name=self.dim_to_coordinate)
diff --git a/extreme_fit/model/utils.py b/extreme_fit/model/utils.py
index 89a11766..30147834 100644
--- a/extreme_fit/model/utils.py
+++ b/extreme_fit/model/utils.py
@@ -4,6 +4,7 @@ import io
import os.path as op
import random
import warnings
+from collections import OrderedDict
from contextlib import redirect_stdout
from typing import Dict
@@ -25,6 +26,7 @@ pandas2ri.activate()
r.library('SpatialExtremes')
r.library('data.table')
r.library('quantreg')
+r.library('evgam')
# Desactivate temporarily warnings
default_filters = warnings.filters.copy()
warnings.filterwarnings("ignore")
@@ -151,10 +153,9 @@ def get_coord(df_coordinates: pd.DataFrame):
return coord
-def get_coord_df(df_coordinates: pd.DataFrame):
- coord = pandas2ri.py2ri_pandasdataframe(df_coordinates)
- # coord = r.transpose(coord)
- colname = df_coordinates.columns
+def get_r_dataframe_from_python_dataframe(df: pd.DataFrame):
+ coord = pandas2ri.py2ri_pandasdataframe(df)
+ colname = df.columns
coord.colnames = r.c(colname)
coord = r('data.frame')(coord, stringsAsFactors=True)
return coord
@@ -179,14 +180,18 @@ def new_coef_name_to_old_coef_names():
}
-def get_margin_formula_extremes(fit_marge_form_dict) -> Dict:
+def get_margin_formula_extremes(fit_marge_form_dict, transformed_as_formula=True) -> Dict:
v_to_str = lambda v: ' '.join(v.split()[2:]) if v != 'NULL' else ' 1'
form_dict = {
k: '~ ' + ' + '.join(
[v_to_str(fit_marge_form_dict[e]) for e in l if e in fit_marge_form_dict])
for k, l in new_coef_name_to_old_coef_names().items()
}
- return {k: robjects.Formula(v) for k, v in form_dict.items()}
+ res = OrderedDict()
+ ordered_keys = ['location.fun', 'scale.fun', 'shape.fun']
+ for k in ordered_keys:
+ res[k] = robjects.Formula(form_dict[k]) if transformed_as_formula else form_dict[k]
+ return res
# def conversion_to_FloatVector(data):
# """Convert DataFrame or numpy array into FloatVector for r"""
diff --git a/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py
new file mode 100644
index 00000000..71a69072
--- /dev/null
+++ b/test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py
@@ -0,0 +1,93 @@
+import unittest
+
+import numpy as np
+import pandas as pd
+
+from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator
+from extreme_fit.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
+from extreme_fit.model.margin_model.utils import \
+ MarginFitMethod
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel, \
+ NonStationaryLocationTemporalModel
+from extreme_fit.model.utils import r, set_seed_r
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_coordinates import \
+ AbstractTemporalCoordinates
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+ AbstractSpatioTemporalObservations
+from test.test_utils import load_non_stationary_temporal_margin_models
+
+
+class TestGevTemporalEvGam(unittest.TestCase):
+
+ def setUp(self) -> None:
+ set_seed_r()
+ r("""
+ N <- 50
+ loc = 0; scale = 2; shape <- 1
+ x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
+ start_loc = 0; start_scale = 1; start_shape = 1
+ """)
+ # Compute the stationary temporal margin with isMev
+ self.start_year = 0
+ df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: range(self.start_year, self.start_year + 50)})
+ self.coordinates = AbstractTemporalCoordinates.from_df(df)
+ df2 = pd.DataFrame(data=np.array(r['x_gev']), index=df.index)
+ observations = AbstractSpatioTemporalObservations(df_maxima_gev=df2)
+ self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
+ self.fit_method = MarginFitMethod.evgam
+
+ def test_gev_temporal_margin_fit_stationary(self):
+ # Create estimator
+ estimator = fitted_linear_margin_estimator(StationaryTemporalModel, self.coordinates, self.dataset,
+ starting_year=0,
+ fit_method=self.fit_method)
+ ref = {'loc': 0.043852626609682636, 'scale': 2.0696383929640807, 'shape': 0.8290699088428524}
+ for year in range(1, 3):
+ mle_params_estimated = estimator.function_from_fit.get_params(np.array([year])).to_dict()
+ for key in ref.keys():
+ self.assertAlmostEqual(ref[key], mle_params_estimated[key], places=3)
+
+ # def test_gev_temporal_margin_fit_nonstationary(self):
+ # # Create estimator
+ # margin_models = load_non_stationary_temporal_margin_models(self.coordinates)
+ # for margin_model in margin_models:
+ # estimator = LinearMarginEstimator(self.dataset, margin_model)
+ # estimator.fit()
+ # # Checks that parameters returned are indeed different
+ # mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
+ # mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+ # self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+ #
+ # def test_gev_temporal_margin_fit_nonstationary_with_start_point(self):
+ # # Create estimator
+ # estimator = self.fit_non_stationary_estimator(starting_point=3)
+ # self.assertNotEqual(estimator.function_from_fit.mu1_temporal_trend, 0.0)
+ # # Checks starting point parameter are well passed
+ # self.assertEqual(3, estimator.function_from_fit.starting_point)
+ # # Checks that parameters returned are indeed different
+ # mle_params_estimated_year1 = estimator.function_from_fit.get_params(np.array([1])).to_dict()
+ # mle_params_estimated_year3 = estimator.function_from_fit.get_params(np.array([3])).to_dict()
+ # self.assertEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+ # mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([5])).to_dict()
+ # self.assertNotEqual(mle_params_estimated_year5, mle_params_estimated_year3)
+ #
+ # def fit_non_stationary_estimator(self, starting_point):
+ # margin_model = NonStationaryLocationTemporalModel(self.coordinates,
+ # starting_point=starting_point + self.start_year)
+ # estimator = LinearMarginEstimator(self.dataset, margin_model)
+ # estimator.fit()
+ # return estimator
+ #
+ # def test_two_different_starting_points(self):
+ # # Create two different estimators
+ # estimator1 = self.fit_non_stationary_estimator(starting_point=3)
+ # estimator2 = self.fit_non_stationary_estimator(starting_point=28)
+ # mu1_estimator1 = estimator1.function_from_fit.mu1_temporal_trend
+ # mu1_estimator2 = estimator2.function_from_fit.mu1_temporal_trend
+ # self.assertNotEqual(mu1_estimator1, mu1_estimator2)
+
+
+if __name__ == '__main__':
+ unittest.main()
--
GitLab