[projection snowfall] add evgam call for...

[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

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

extreme_fit/distribution/gev/evgam_fixed.R

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

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

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)
 
 
 

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
 

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

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],

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]:

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)
 

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

 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,

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 = {

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

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

extreme_fit/model/result_from_model_fit/result_from_extremes/result_from_evgam.py

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

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"""

test/test_extreme_fit/test_estimator/test_temporal_estimator/test_gev_temporal_evgam.py

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