[EXTREME MODEL] add max_stable_with_spline.R & prepare test

experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py

@@ -62,6 +62,11 @@ def annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, alt
             study_visualizer.visualize_annual_mean_values(take_mean_value=take_mean_value)
 
 
+def max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1):
+    study = SafranSnowfall(altitude=altitude, nb_consecutive_days=nb_days)
+    study_visualizer = StudyVisualizer(study)
+    study_visualizer.visualize_brown_resnick_fit()
+
 def normal_visualization():
     save_to_file = False
     only_first_one = True
@@ -91,6 +96,7 @@ def complete_analysis(only_first_one=False):
 
 if __name__ == '__main__':
     # annual_mean_vizu_compare_durand_study(safran=True, take_mean_value=True, altitude=2400)
-    normal_visualization()
+    # normal_visualization()
+    max_stable_process_vizu_compare_gaume_study(altitude=1800, nb_days=1)
     # extended_visualization()
     # complete_analysis()

experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py

@@ -16,6 +16,7 @@ from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator impo
 from extreme_estimator.extreme_models.margin_model.param_function.param_function import ParamFunction
 from extreme_estimator.extreme_models.margin_model.smooth_margin_model import 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 BrownResnick
 from extreme_estimator.margin_fits.abstract_params import AbstractParams
 from extreme_estimator.margin_fits.gev.gev_params import GevParams
 from extreme_estimator.margin_fits.gev.gevmle_fit import GevMleFit
@@ -209,6 +210,10 @@ class StudyVisualizer(object):
         ax.set_xlabel('year')
         ax.set_title(self.study.safran_massif_names[massif_id])
 
+    def visualize_brown_resnick_fit(self):
+        pass
+
+
     def visualize_linear_margin_fit(self, only_first_max_stable=False):
         default_covariance_function = CovarianceFunction.cauchy
         plot_name = 'Full Likelihood with Linear marginals and max stable dependency structure'
@@ -217,9 +222,12 @@ class StudyVisualizer(object):
         self.plot_name = plot_name
         max_stable_models = load_test_max_stable_models(default_covariance_function=default_covariance_function)
         if only_first_max_stable:
-            max_stable_models = max_stable_models[:1]
+            # Keep only the BrownResnick model
+            max_stable_models = max_stable_models[1:2]
         if only_first_max_stable is None:
             max_stable_models = []
+
+
         fig, axes = plt.subplots(len(max_stable_models) + 2, len(GevParams.SUMMARY_NAMES), figsize=self.figsize)
         fig.subplots_adjust(hspace=self.subplot_space, wspace=self.subplot_space)
         margin_class = LinearAllParametersAllDimsMarginModel

extreme_estimator/extreme_models/max_stable_model/max_stable_with_spline.R

+library(SpatialExtremes)
+
+## Not run:
+##Define the coordinate of each location
+n.site <- 30
+locations <- matrix(runif(2*n.site, 0, 10), ncol = 2)
+colnames(locations) <- c("lon", "lat")
+
+##  Fitting a penalized smoothing splines for the margins with the
+##     Smith's model
+data <- rmaxstab(100, locations, cov.mod = "gauss", cov11 = 100, cov12 =
+                 25, cov22 = 220)
+
+##     And transform it to ordinary GEV margins with a non-linear
+##     function
+fun <- function(x)
+  2 * sin(pi * x / 4) + 10
+fun2 <- function(x)
+  (fun(x) - 7 ) / 15
+
+param.loc <- fun(locations[,2])
+param.scale <- fun(locations[,2])
+param.shape <- fun2(locations[,1])
+
+##Transformation from unit Frechet to common GEV margins
+for (i in 1:n.site)
+  data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i],
+param.shape[i])
+
+##Defining the knots, penalty, degree for the splines
+n.knots_x = 2
+n.knots_y = 2
+knots = quantile(locations[,1], prob=1:n.knots_x/(n.knots_x+1))
+knots2 = quantile(locations[,2], prob=1:n.knots_y/(n.knots_y+1))
+knots_tot = cbind(knots,knots2)
+print(knots)
+print(knots2)
+print(knots_tot)
+
+##Be careful the choice of the penalty (i.e. the smoothing parameter)
+##may strongly affect the result Here we use p-splines for each GEV
+##parameter - so it's really CPU demanding but one can use 1 p-spline
+##and 2 linear models.
+##A simple linear model will be clearly faster...
+loc.form <- y ~ rb(locations[,1], knots = knots, degree = 3, penalty = .5)
+scale.form <- y ~ rb(locations[,2], knots = knots2, degree = 3, penalty = .5)
+shape.form <- y ~ rb(locations, knots = knots_tot, degree = 3, penalty = .5)
+
+fitted <- fitmaxstab(data, locations, "gauss", loc.form, scale.form, shape.form,
+                     method = "BFGS")
+fitted
+
+
+## End(Not run)
\ No newline at end of file

test/test_unitary/test_fitmaxstab/test_fitmaxstab_with_margin.py

 import unittest
 
-from extreme_estimator.estimator.full_estimator.abstract_full_estimator import FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
+    FullEstimatorInASingleStepWithSmoothMargin
 from extreme_estimator.extreme_models.margin_model.smooth_margin_model import ConstantMarginModel, \
     LinearMarginModelExample
 from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model import CovarianceFunction
@@ -65,5 +66,38 @@ class TestMaxStableFitWithLinearMargin(TestUnitaryAbstract):
         self.compare()
 
 
+# class TestMaxStableFitWithSpline(TestUnitaryAbstract):
+#
+#     @property
+#     def r_output(self):
+#         TestRMaxStabWithMarginConstant.r_code()
+#         r("""
+#         # Code inspired from Johan code
+#         n.knots_x = 1
+#         n.knots_y = 2
+#         knots = quantile(locations[,1], prob=1:n.knots_x/(n.knots_x+1))
+#         knots2 = quantile(locations[,2], prob=1:n.knots_y/(n.knots_y+1))
+#         knots_tot = cbind(knots,knots2)
+#
+#         loc.form <- y ~ rb(locations[,1], knots = knots, degree = 3, penalty = .5)
+#         scale.form <- y ~ rb(locations[,2], knots = knots2, degree = 3, penalty = .5)
+#         shape.form <- y ~ rb(locations, knots = knots_tot, degree = 3, penalty = .5)
+#         """)
+#         return self.r_fitted_values_from_res_variable
+#
+#     @property
+#     def python_output(self):
+#         dataset = TestRMaxStabWithMarginConstant.python_code()
+#         max_stable_model = Schlather(covariance_function=CovarianceFunction.whitmat, use_start_value=False)
+#         margin_model = LinearMarginModelExample(dataset.coordinates)
+#         full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model,
+#                                                                     max_stable_model)
+#         full_estimator.fit()
+#         return full_estimator.fitted_values
+#
+#     def test_max_stable_fit_with_spline_margin(self):
+#         self.compare()
+
+
 if __name__ == '__main__':
     unittest.main()