From 67e36996347512eafa252dfa997b9492e7bd306f Mon Sep 17 00:00:00 2001
From: Le Roux Erwan <erwan.le-roux@irstea.fr>
Date: Mon, 23 Mar 2020 20:53:59 +0100
Subject: [PATCH] [quantile regression project] add annual maxima osbervation
 from daily exponential. add sammple_r_function as argument to
 rmargin_from_nb_obs

---
 .../margin_model/abstract_margin_model.py     |  5 +--
 .../annual_maxima_observations.py             | 18 +++++-----
 .../daily_observations.py                     | 26 +++++++++++++++
 .../test_spatio_temporal_observations.py      | 33 +++++++++++++++++--
 4 files changed, 68 insertions(+), 14 deletions(-)
 create mode 100644 spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py

diff --git a/extreme_fit/model/margin_model/abstract_margin_model.py b/extreme_fit/model/margin_model/abstract_margin_model.py
index ef7ae1a7..b4ac6918 100644
--- a/extreme_fit/model/margin_model/abstract_margin_model.py
+++ b/extreme_fit/model/margin_model/abstract_margin_model.py
@@ -69,11 +69,12 @@ class AbstractMarginModel(AbstractModel, ABC):
         maxima_gev = self.frech2gev(maxima_frech, coordinates_values, self.margin_function_sample)
         return maxima_gev
 
-    def rmargin_from_nb_obs(self, nb_obs: int, coordinates_values: np.ndarray) -> np.ndarray:
+    def rmargin_from_nb_obs(self, nb_obs: int, coordinates_values: np.ndarray,
+                            sample_r_function='rgev') -> np.ndarray:
         maxima_gev = []
         for coordinate in coordinates_values:
             gev_params = self.margin_function_sample.get_gev_params(coordinate)
-            x_gev = r.rgev(nb_obs, **gev_params.to_dict())
+            x_gev = r(sample_r_function)(nb_obs, **gev_params.to_dict())
             maxima_gev.append(x_gev)
         return np.array(maxima_gev)
 
diff --git a/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py b/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py
index 93474165..a4dbcde9 100644
--- a/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py
+++ b/spatio_temporal_dataset/spatio_temporal_observations/annual_maxima_observations.py
@@ -9,6 +9,7 @@ from spatio_temporal_dataset.coordinates.spatio_temporal_coordinates.abstract_sp
     AbstractSpatioTemporalCoordinates
 from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations \
     import AbstractSpatioTemporalObservations
+from spatio_temporal_dataset.spatio_temporal_observations.daily_observations import DailyExp
 
 
 class AnnualMaxima(AbstractSpatioTemporalObservations):
@@ -25,7 +26,8 @@ class MarginAnnualMaxima(AnnualMaxima):
     def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates,
                       margin_model: AbstractMarginModel):
         maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs,
-                                                      coordinates_values=coordinates.coordinates_values())
+                                                      coordinates_values=coordinates.coordinates_values(),
+                                                      sample_r_function='rgev')
         df_maxima_gev = pd.DataFrame(data=maxima_gev, index=coordinates.index)
         return cls(df_maxima_gev=df_maxima_gev)
 
@@ -35,15 +37,11 @@ class DailyExpAnnualMaxima(AnnualMaxima):
     @classmethod
     def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates,
                       margin_model: AbstractMarginModel):
-        pass
-
-
-class DailyExp(AbstractSpatioTemporalObservations):
-
-    @classmethod
-    def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates,
-                      margin_model: AbstractMarginModel):
-        pass
+        # todo: to take nb_obs into accoutn i could generate nb_obs * 365 observations
+        observations = DailyExp.from_sampling(nb_obs=365, coordinates=coordinates, margin_model=margin_model)
+        df_daily_values = observations.df_maxima_gev
+        df_maxima_gev = pd.DataFrame({'0': df_daily_values.max(axis=1)}, index=df_daily_values.index)
+        return cls(df_maxima_gev=df_maxima_gev)
 
 
 class MaxStableAnnualMaxima(AnnualMaxima):
diff --git a/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py b/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py
new file mode 100644
index 00000000..9483f3b3
--- /dev/null
+++ b/spatio_temporal_dataset/spatio_temporal_observations/daily_observations.py
@@ -0,0 +1,26 @@
+import pandas as pd
+
+from extreme_fit.distribution.abstract_params import AbstractParams
+from extreme_fit.distribution.gev.gev_params import GevParams
+from extreme_fit.model.margin_model.abstract_margin_model import AbstractMarginModel
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
+from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \
+    ConsecutiveTemporalCoordinates
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+    AbstractSpatioTemporalObservations
+
+
+class DailyObservations(AbstractSpatioTemporalObservations):
+    pass
+
+
+class DailyExp(AbstractSpatioTemporalObservations):
+
+    @classmethod
+    def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates,
+                      margin_model: AbstractMarginModel):
+        exponential_values = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs,
+                                                              coordinates_values=coordinates.coordinates_values())
+        df_exponential_values = pd.DataFrame(data=exponential_values, index=coordinates.index)
+        return cls(df_maxima_gev=df_exponential_values)
diff --git a/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py b/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py
index ed6c4d40..716a5856 100644
--- a/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py
+++ b/test/test_spatio_temporal_dataset/test_spatio_temporal_observations.py
@@ -1,9 +1,17 @@
 import unittest
-import numpy as np
 
+import numpy as np
 import pandas as pd
 
-from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import AbstractSpatioTemporalObservations
+from extreme_fit.distribution.abstract_params import AbstractParams
+from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
+from extreme_fit.model.utils import set_seed_for_test
+from spatio_temporal_dataset.coordinates.temporal_coordinates.generated_temporal_coordinates import \
+    ConsecutiveTemporalCoordinates
+from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_temporal_observations import \
+    AbstractSpatioTemporalObservations
+from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import DailyExpAnnualMaxima
+from spatio_temporal_dataset.spatio_temporal_observations.daily_observations import DailyExp
 
 
 class TestSpatioTemporalObservations(unittest.TestCase):
@@ -18,5 +26,26 @@ class TestSpatioTemporalObservations(unittest.TestCase):
         self.assertTrue(np.equal(example, maxima_frech).all(), msg="{} {}".format(example, maxima_frech))
 
 
+class TestDailyObservations(unittest.TestCase):
+
+    def setUp(self) -> None:
+        super().setUp()
+        set_seed_for_test(seed=42)
+        self.coordinates = ConsecutiveTemporalCoordinates.from_nb_temporal_steps(nb_temporal_steps=10)
+        gev_param_name_to_coef_list = {
+            AbstractParams.SCALE: [1],
+        }
+        self.margin_model = StationaryTemporalModel.from_coef_list(self.coordinates, gev_param_name_to_coef_list)
+
+    def test_exponential_observations(self):
+        obs = DailyExp.from_sampling(nb_obs=1, coordinates=self.coordinates,
+                                     margin_model=self.margin_model)
+        self.assertAlmostEqual(obs.df_maxima_gev.mean()[0], 4.692385276235156)
+
+    def test_annual_maxima_observations_from_daily_observations(self):
+        obs = DailyExpAnnualMaxima.from_sampling(1, self.coordinates, self.margin_model)
+        self.assertAlmostEqual(obs.df_maxima_gev.mean()[0], 1183.8468374768636)
+
+
 if __name__ == '__main__':
     unittest.main()
-- 
GitLab