[contrasting] add abstract spatio temporal polynomial model

extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R

@@ -26,7 +26,8 @@ colnames(coord) = c("X", "T")
 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, location.fun= ~sin(X) + cos(T), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
+# res = fevd_fixed(x_gev, data=coord, location.fun= ~sin(X) + cos(T), method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~T * X + X +T , method='MLE', type="GEV", verbose=FALSE, use.phi=FALSE)
 print(res)
 
 # Some display for the results

extreme_fit/estimator/margin_estimator/abstract_margin_estimator.py

@@ -28,13 +28,18 @@ class LinearMarginEstimator(AbstractMarginEstimator):
         self.margin_model = margin_model
 
     def _fit(self) -> AbstractResultFromModelFit:
-        df_coordinates_spat = self.dataset.coordinates.df_spatial_coordinates(self.train_split)
         return self.margin_model.fitmargin_from_maxima_gev(data=self.maxima_gev_train,
-                                                           df_coordinates_spat=df_coordinates_spat,
-                                                           df_coordinates_temp=self.coordinate_temp)
+                                                           df_coordinates_spat=self.df_coordinates_spat,
+                                                           df_coordinates_temp=self.df_coordinates_temp)
 
     @property
-    def coordinate_temp(self):
+    def df_coordinates_spat(self):
+        return self.dataset.coordinates.df_spatial_coordinates(split=self.train_split,
+                                                               drop_duplicates=self.margin_model.drop_duplicates)
+
+
+    @property
+    def df_coordinates_temp(self):
         return self.dataset.coordinates.df_temporal_coordinates_for_fit(split=self.train_split,
                                                                         starting_point=self.margin_model.starting_point,
                                                                         drop_duplicates=self.margin_model.drop_duplicates)

extreme_fit/model/margin_model/linear_margin_model/abstract_temporal_linear_margin_model.py

@@ -35,20 +35,23 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
 
     def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
                                   df_coordinates_temp: pd.DataFrame) -> AbstractResultFromModelFit:
-        data = data[0]
-        assert len(data) == len(df_coordinates_temp.values), 'len(data)={} != len(temp)={}'.format(len(data),
-                                                                                                   len(
-                                                                                                       df_coordinates_temp.values))
-        x = ro.FloatVector(data)
+        if len(data) == 1:
+            data = data[0]
+            assert len(data) == len(df_coordinates_temp.values), 'len(data)={} != len(temp)={}'.format(len(data),
+                                                                                                       len(
+                                                                                                           df_coordinates_temp.values))
+            x = ro.FloatVector(data)
+        else:
+            x = ro.Matrix(data)
         if self.params_class is GevParams:
             if self.fit_method == MarginFitMethod.is_mev_gev_fit:
                 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)
+                return self.extremes_fevd_fit(x, df_coordinates_temp, df_coordinates_spat)
             else:
                 raise NotImplementedError
         elif self.params_class is ExpParams:
-            return self.extreme_fevd_mle_exp_fit(x, df_coordinates_temp)
+            return self.extreme_fevd_mle_exp_fit(x, df_coordinates_temp, df_coordinates_spat)
         else:
             raise NotImplementedError
 
@@ -63,21 +66,22 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
 
     # Gev fit with extRemes package
 
-    def extremes_fevd_fit(self, x, df_coordinates_temp) -> AbstractResultFromExtremes:
+    def extremes_fevd_fit(self, x, df_coordinates_temp, df_coordinates_spat) -> AbstractResultFromExtremes:
         assert self.fit_method in FEVD_MARGIN_FIT_METHODS
         if self.fit_method == MarginFitMethod.extremes_fevd_bayesian:
             return self.extremes_fevd_bayesian_fit(x, df_coordinates_temp)
         else:
             return self.run_fevd_fixed(df_coordinates_temp=df_coordinates_temp,
+                                       df_coordinates_spat=df_coordinates_spat,
                                        method=FEVD_MARGIN_FIT_METHOD_TO_ARGUMENT_STR[self.fit_method], x=x)
 
-    def extreme_fevd_mle_exp_fit(self, x, df_coordinates_temp):
-        return self.run_fevd_fixed(df_coordinates_temp, "Exponential", x)
+    def extreme_fevd_mle_exp_fit(self, x, df_coordinates_temp, df_coordinates_spat):
+        return self.run_fevd_fixed(df_coordinates_temp, df_coordinates_spat, "Exponential", x)
 
-    def run_fevd_fixed(self, df_coordinates_temp, method, x):
+    def run_fevd_fixed(self, df_coordinates_temp, df_coordinates_spat, method, x):
         if self.fit_method == MarginFitMethod.extremes_fevd_l_moments:
             assert self.margin_function.is_a_stationary_model
-        r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp)
+        r_type_argument_kwargs, y = self.extreme_arguments(df_coordinates_temp, df_coordinates_spat)
         res = safe_run_r_estimator(function=r('fevd_fixed'),
                                    x=x,
                                    data=y,
@@ -107,13 +111,22 @@ class AbstractTemporalLinearMarginModel(LinearMarginModel):
                                    )
         return ResultFromBayesianExtremes(res, self.param_name_to_list_for_result)
 
-    def extreme_arguments(self, df_coordinates_temp):
+    def extreme_arguments(self, df_coordinates_temp, df_coordinates_spat=None):
+
+        # Load parameters
+
+        if df_coordinates_spat is None or df_coordinates_spat.empty:
+            df = df_coordinates_temp
+        else:
+            df = pd.concat([df_coordinates_spat, df_coordinates_temp], axis=1)
+            print(df.shape)
+            print(df.head())
+        y = get_coord_df(df)
+
         # Disable the use of log sigma parametrization
         r_type_argument_kwargs = {'use.phi': False,
                                   'verbose': False}
-        # Load parameters
         r_type_argument_kwargs.update(get_margin_formula_extremes(self.margin_function.form_dict))
-        y = get_coord_df(df_coordinates_temp)
         return r_type_argument_kwargs, y
 
     # Default arguments for all methods

extreme_fit/model/margin_model/polynomial_margin_model/spatio_temporal_polynomial_model.py

+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import PolynomialMarginModel
+from extreme_fit.model.margin_model.utils import MarginFitMethod
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+
+
+class AbstractSpatioTemporalPolynomialModel(PolynomialMarginModel):
+
+    def __init__(self, coordinates: AbstractCoordinates, params_user=None, starting_point=None,
+                 fit_method=MarginFitMethod.extremes_fevd_mle, nb_iterations_for_bayesian_fit=5000,
+                 params_initial_fit_bayesian=None, type_for_MLE="GEV", params_class=GevParams, max_degree=2):
+        super().__init__(coordinates, params_user, starting_point, fit_method, nb_iterations_for_bayesian_fit,
+                         params_initial_fit_bayesian, type_for_MLE, params_class, max_degree)
+        self.drop_duplicates = False
+
+
+class NonStationaryLocationSpatioTemporalLinearityModel(AbstractSpatioTemporalPolynomialModel):
+
+    def load_margin_function(self, param_name_to_dims=None):
+        return super().load_margin_function({GevParams.LOC: [
+            (self.coordinates.idx_temporal_coordinates, 1),
+            (self.coordinates.idx_x_coordinates, 1),
+        ]})

projects/altitude_spatial_model/altitudes_fit/two_fold_detail_fit.py

@@ -69,7 +69,8 @@ class TwoFoldModelFit(object):
         self.model_class = model_class
         self.fit_method = fit_method
         self.estimators = [fitted_linear_margin_estimator_short(model_class=self.model_class, dataset=dataset,
-                                                 fit_method=self.fit_method) for dataset in two_fold_datasets] # type: List[LinearMarginEstimator]
+                                                                fit_method=self.fit_method)
+                           for dataset in two_fold_datasets]  # type: List[LinearMarginEstimator]
         self.estimator_fold_1 = self.estimators[0]
         self.estimator_fold_2 = self.estimators[1]
 

spatio_temporal_dataset/coordinates/abstract_coordinates.py

@@ -210,11 +210,12 @@ class AbstractCoordinates(object):
     def has_spatial_coordinates(self) -> bool:
         return self.nb_spatial_coordinates > 0
 
-    def df_spatial_coordinates(self, split: Split = Split.all, transformed=True) -> pd.DataFrame:
+    def df_spatial_coordinates(self, split: Split = Split.all, transformed=True, drop_duplicates=True) -> pd.DataFrame:
         if self.nb_spatial_coordinates == 0:
             return pd.DataFrame()
         else:
-            return self.df_coordinates(split, transformed).loc[:, self.spatial_coordinates_names].drop_duplicates()
+            df = self.df_coordinates(split, transformed).loc[:, self.spatial_coordinates_names]
+            return df.drop_duplicates() if drop_duplicates else df
 
     def nb_stations(self, split: Split = Split.all) -> int:
         return len(self.df_spatial_coordinates(split))
@@ -302,7 +303,11 @@ class AbstractCoordinates(object):
     def idx_temporal_coordinates(self):
         return self.coordinates_names.index(self.COORDINATE_T)
 
-    # Spatio temporal attributes
+    @property
+    def idx_x_coordinates(self):
+        return self.coordinates_names.index(self.COORDINATE_X)
+
+# Spatio temporal attributes
 
     @property
     def has_spatio_temporal_coordinates(self) -> bool:

test/test_extreme_fit/test_estimator/test_gev_spatio_temporal_extremes_mle.py

+import unittest
+
+import numpy as np
+import pandas as pd
+
+from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.polynomial_margin_model.polynomial_margin_model import \
+    NonStationaryQuadraticLocationModel, \
+    NonStationaryQuadraticScaleModel, NonStationaryQuadraticLocationGumbelModel, NonStationaryQuadraticScaleGumbelModel
+from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_polynomial_model import \
+    NonStationaryLocationSpatioTemporalLinearityModel
+from extreme_trend.abstract_gev_trend_test import fitted_linear_margin_estimator
+from extreme_fit.model.margin_model.utils import \
+    MarginFitMethod
+from extreme_fit.model.utils import r, set_seed_r
+from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
+from projects.altitude_spatial_model.altitudes_fit.two_fold_datasets_generator import TwoFoldDatasetsGenerator
+from projects.altitude_spatial_model.altitudes_fit.two_fold_fit import TwoFoldFit
+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_projects.test_contrasting.test_two_fold_fit import load_two_fold_fit
+
+
+class TestGevTemporalQuadraticExtremesMle(unittest.TestCase):
+
+    def get_estimator_fitted(self, model_class):
+        altitudes = [900, 1200]
+        study_class = SafranSnowfall1Day
+        studies = AltitudesStudies(study_class, altitudes, year_max=2019)
+        two_fold_datasets_generator = TwoFoldDatasetsGenerator(studies, nb_samples=1, massif_names=['Vercors'])
+        model_family_name_to_model_class = {'Non stationary': [model_class]}
+        two_fold_fit = TwoFoldFit(two_fold_datasets_generator=two_fold_datasets_generator,
+                                  model_family_name_to_model_classes=model_family_name_to_model_class,
+                                  fit_method=MarginFitMethod.extremes_fevd_mle)
+        massif_fit = two_fold_fit.massif_name_to_massif_fit['Vercors']
+        sample_fit = massif_fit.sample_id_to_sample_fit[0]
+        model_fit = sample_fit.model_class_to_model_fit[model_class]  # type: TwoFoldModelFit
+        estimator = model_fit.estimator_fold_1
+        print(estimator.dataset)
+        return estimator
+
+    # def test_location_spatio_temporal_linearity(self):
+    #     # Create estimator
+    #     # estimator = fitted_linear_margin_estimator(model_class,
+    #     #                                            self.coordinates, self.dataset,
+    #     #                                            starting_year=0,
+    #     #                                            fit_method=self.fit_method)
+    #     # # 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([21])).to_dict()
+    #     # mle_params_estimated_year5 = estimator.function_from_fit.get_params(np.array([41])).to_dict()
+    #     # self.assertNotEqual(mle_params_estimated_year1, mle_params_estimated_year3)
+    #     # self.assertNotEqual(mle_params_estimated_year3, mle_params_estimated_year5)
+    #     # # Assert the relationship for the location is indeed quadratic
+    #     # diff1 = mle_params_estimated_year1[quadratic_param] - mle_params_estimated_year3[quadratic_param]
+    #     # diff2 = mle_params_estimated_year3[quadratic_param] - mle_params_estimated_year5[quadratic_param]
+    #     # self.assertNotAlmostEqual(diff1, diff2)
+    #     estimator = self.get_estimator_fitted(NonStationaryLocationSpatioTemporalLinearityModel)
+    #     # Assert that indicators are correctly computed
+    #     self.assertAlmostEqual(estimator.result_from_model_fit.nllh, estimator.nllh())
+    #     self.assertAlmostEqual(estimator.result_from_model_fit.aic, estimator.aic())
+    #     self.assertAlmostEqual(estimator.result_from_model_fit.bic, estimator.bic())
+
+
+if __name__ == '__main__':
+    unittest.main()