[contrasting] add xlabel, remove 300m from the lower altitude group. add...

[contrasting] add xlabel, remove 300m from the lower altitude group. add condition that we consider only altitude for which the massif has less than 10% of non zero annual maxima. return period is changed to 100.

extreme_data/meteo_france_data/scm_models_data/visualization/plot_utils.py

@@ -35,7 +35,8 @@ def plot_against_altitude(x_ticks, ax, massif_id, massif_name, values, altitude=
 
 
 def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_method, add_x_label=True,
-              negative_and_positive_values=True, massif_name_to_text=None, add_colorbar=True, max_abs_change=None):
+              negative_and_positive_values=True, massif_name_to_text=None, add_colorbar=True, max_abs_change=None,
+              xlabel=None):
     if max_abs_change is None:
         max_abs_change = max([abs(e) for e in massif_name_to_value.values()])
     if negative_and_positive_values:
@@ -82,5 +83,8 @@ def load_plot(cmap, graduation, label, massif_name_to_value, altitude, fit_metho
     ax.get_xaxis().set_visible(True)
     ax.set_xticks([])
     if add_x_label:
-        ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
+        if xlabel is None:
+            ax.set_xlabel('Altitude = {}m'.format(altitude), fontsize=15)
+        else:
+            ax.set_xlabel(xlabel, fontsize=10)
         # ax.set_title('Fit method is {}'.format(fit_method))

extreme_data/meteo_france_data/scm_models_data/visualization/study_visualizer.py

@@ -717,14 +717,14 @@ class StudyVisualizer(VisualizationParameters):
 
     def plot_map(self, cmap, fit_method, graduation, label, massif_name_to_value, plot_name, add_x_label=True,
                  negative_and_positive_values=True, massif_name_to_text=None, altitude=None, add_colorbar=True,
-                 max_abs_change=None):
+                 max_abs_change=None, xlabel=None):
         if altitude is None:
             altitude = self.study.altitude
         if len(massif_name_to_value) > 0:
             load_plot(cmap, graduation, label, massif_name_to_value, altitude, fitmethod_to_str(fit_method),
                       add_x_label=add_x_label, negative_and_positive_values=negative_and_positive_values,
                       massif_name_to_text=massif_name_to_text,
-                      add_colorbar=add_colorbar, max_abs_change=max_abs_change)
+                      add_colorbar=add_colorbar, max_abs_change=max_abs_change, xlabel=xlabel)
             self.plot_name = plot_name
             # self.show_or_save_to_file(add_classic_title=False, tight_layout=True, no_title=True, dpi=500)
             self.show_or_save_to_file(add_classic_title=False, no_title=True, dpi=500, tight_layout=True)

extreme_fit/distribution/gev/ci_fevd_fixed.R

@@ -424,6 +424,7 @@ ci.rl.ns.fevd.mle_fixed <- function (x, alpha = 0.05, return.period = 100, metho
             qcov.base = qcov.base)
         z.alpha <- qnorm(alpha/2, lower.tail = FALSE)
         cov.theta <- parcov.fevd(x)
+        # cov.theta <- parcov.fevd_fixed(x)
         if (is.null(cov.theta))
             stop("ci: Sorry, unable to calculate the parameter covariance matrix.  Maybe try a different method.")
         var.theta <- diag(cov.theta)
@@ -454,6 +455,60 @@ ci.rl.ns.fevd.mle_fixed <- function (x, alpha = 0.05, return.period = 100, metho
     return(out)
 }
 
+parcov.fevd_fixed <- function (x)
+{
+    cov.theta <- NULL
+    theta.hat <- x$results$par
+    theta.names <- names(theta.hat)
+    if (is.element("log.scale", theta.names)) {
+        theta.hat[theta.names == "log.scale"] <- exp(theta.hat[theta.names ==
+            "log.scale"])
+        theta.names[theta.names == "log.scale"] <- "scale"
+        names(theta.hat) <- theta.names
+        phiU <- FALSE
+    }
+    else phiU <- x$par.models$log.scale
+    y <- datagrabber(x)
+    if (x$data.name[2] != "") {
+        xdat <- y[, 1]
+        data <- y[, -1]
+    }
+    else {
+        xdat <- y[, 1]
+        data <- NULL
+    }
+    designs <- setup.design(x)
+    if (x$method != "GMLE") {
+        hold <- try(suppressWarnings(optimHess(theta.hat, oevd,
+            gr = grlevd, o = x, des = designs, x = xdat, data = data,
+            u = x$threshold, npy = x$npy, phi = phiU, blocks = x$blocks)),
+            silent = TRUE)
+        if ((class(hold) != "try-error") && all(!is.na(hold))) {
+            cov.theta <- try(suppressWarnings(solve(hold)), silent = TRUE)
+            print(cov.theta)
+            # if (any(diag(cov.theta) <= 0))
+            #     re.do <- TRUE
+            # else re.do <- FALSE
+            re.do = TRUE
+        }
+        else re.do <- TRUE
+    }
+    else re.do <- TRUE
+    if (re.do) {
+        hold <- try(optimHess(theta.hat, oevd, o = x, des = designs,
+            x = xdat, data = data, u = x$threshold, npy = x$npy,
+            phi = phiU, blocks = x$blocks), silent = TRUE)
+        if (class(hold) != "try-error" && all(!is.na(hold))) {
+            cov.theta <- try(solve(hold), silent = TRUE)
+            if (class(cov.theta) == "try-error" || any(diag(cov.theta) <=
+                0))
+                cov.theta <- NULL
+        }
+    }
+    return(cov.theta)
+}
+
+
 return.level.ns.fevd.mle_fixed <- function (x, return.period = c(2, 20, 100), ..., alpha = 0.05,
     method = c("normal"), do.ci = FALSE, verbose = FALSE, qcov = NULL,
     qcov.base = NULL)

extreme_fit/distribution/gev/main_fevd_mle_two_covariates.R

@@ -11,8 +11,8 @@ source('ci_fevd_fixed.R')
 source('summary_fevd_fixed.R')
 # Sample from a GEV
 set.seed(42)
-N <- 50
-loc = 0; scale = 1; shape <- 1
+N <- 100
+loc = 0; scale = 1; shape <- 0.1
 x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)
 # start_loc = 0; start_scale = 1; start_shape = 1
 # N <- 50
@@ -22,14 +22,14 @@ print(N)
 coord <- matrix(ncol=2, nrow = N)
 coord[,1]=seq(0,N-1,1)
 coord[,2]=seq(0,N-1,1)
-print(coord)
+
 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= ~poly(X * T, 1, raw = TRUE),  method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
-res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 2, raw = TRUE) , method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
+res = fevd_fixed(x_gev, data=coord, location.fun= ~poly(X, 1, raw = TRUE) + poly(T, 1, raw = TRUE) , method='MLE', type="Gumbel", verbose=FALSE, use.phi=FALSE)
 # print(res)
 print(summary.fevd.mle_fixed(res))
 # print(summary(res)$AIC)
@@ -44,6 +44,18 @@ print(summary.fevd.mle_fixed(res))
 # print(m[1])
 
 
+# v = make.qcov(res, vals = list(mu0 = c(0.0), mu1 = c(0.0), mu2 = c(0.0), sigma = c(0.0)))
+v = make.qcov(res, vals = list(mu1 = c(0.0), mu2 = c(0.0)))
+
+res_ci = ci.fevd.mle_fixed(res, alpha = 0.05, type = c("return.level"),
+    return.period = 50, method = "normal", xrange = NULL, nint = 20, verbose = FALSE,
+    tscale = FALSE, return.samples = FALSE, qcov=v)
+print(res_ci)
+
+
+
+
+
 # Confidence interval staionary
 # method = "proflik"
 # res_ci = ci.fevd.mle(res, alpha = 0.05, type = c("return.level"),

projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py

@@ -43,9 +43,11 @@ class AltitudesStudies(object):
     # Dataset Loader
 
     def spatio_temporal_dataset(self, massif_name, s_split_spatial: pd.Series = None,
-                                s_split_temporal: pd.Series = None):
+                                s_split_temporal: pd.Series = None,
+                                massif_altitudes=None):
         coordinate_values_to_maxima = {}
-        massif_altitudes = self.massif_name_to_altitudes[massif_name]
+        if massif_altitudes is None:
+            massif_altitudes = self.massif_name_to_altitudes[massif_name]
         if len(massif_altitudes) == 0:
             print('{} has no data for these altitudes: {}'.format(massif_name, self.altitudes))
         for altitude in massif_altitudes:

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

@@ -9,7 +9,6 @@ from projects.altitude_spatial_model.altitudes_fit.plot_histogram_altitude_studi
 mpl.rcParams['text.usetex'] = True
 mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}']
 
-
 from projects.altitude_spatial_model.altitudes_fit.utils_altitude_studies_visualizer import load_visualizer_list
 
 from extreme_fit.model.margin_model.polynomial_margin_model.utils import \
@@ -29,13 +28,20 @@ from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_s
 def main():
     study_classes = [SafranSnowfall1Day, SafranSnowfall3Days, SafranSnowfall5Days, SafranSnowfall7Days][:1]
     # study_classes = [SafranPrecipitation1Day][:1]
-
-    massif_names = None
-    # massif_names = ['Mercantour', 'Vercors', 'Ubaye']
-
     seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
 
-    main_loop(altitudes_for_groups[:], massif_names, seasons, study_classes)
+    fast = None
+    if fast is None:
+        massif_names = None
+        altitudes_list = altitudes_for_groups[:2]
+    elif fast:
+        massif_names = ['Mercantour', 'Vercors', 'Ubaye']
+        altitudes_list = altitudes_for_groups[:2]
+    else:
+        massif_names = None
+        altitudes_list = altitudes_for_groups
+
+    main_loop(altitudes_list, massif_names, seasons, study_classes)
 
 
 def main_loop(altitudes_list, massif_names, seasons, study_classes):
@@ -59,7 +65,7 @@ def plot_visualizers(massif_names, visualizer_list):
 
 def plot_visualizer(massif_names, visualizer):
     # Plot time series
-    # visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
+    visualizer.studies.plot_maxima_time_series(massif_names=massif_names)
     # Plot moments against altitude
     # for std in [True, False][:]:
     #     for change in [True, False, None]:

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py

@@ -2,26 +2,13 @@ from enum import Enum
 
 # The order is important
 altitudes_for_groups = [
-    [300, 600, 900],
+    [300, 600, 900][1:],
     [1200, 1500, 1800],
     [2100, 2400, 2700],
     [3000, 3300, 3600, 3900]
 ]
 
 
-# altitudes_for_groups = [
-#     [900, 1200, 1500],
-#     [1800, 2100, 2400],
-#     [2700, 3000, 3300]
-# ]
-
-# altitudes_for_groups = [
-#     [600, 900, 1200, 1500, 1800],
-#     [1500, 1800, 2100, 2400, 2700],
-#     [2400, 2700, 3000, 3300, 3600]
-# ]
-
-
 class AbstractAltitudeGroup(object):
 
     @property
@@ -32,6 +19,18 @@ class AbstractAltitudeGroup(object):
     def reference_altitude(self):
         raise NotImplementedError
 
+    @property
+    def xlabel(self):
+        # warning: the label could not correspond to all massifs, some might have been fitted with less data
+        # idx = get_index_group_from_reference_altitude(reference_altitude)
+        # min_altitude, *_, max_altitude = altitudes_for_groups[idx]
+        i = self.reference_altitude // 1000
+        min_altitude, max_altitude = 1000 * i, 1000 * (i + 1)
+        return 'Altitude = {} m\n' \
+               'Estimated with maxima between {} m and {} m'.format(self.reference_altitude,
+                                                                        min_altitude,
+                                                                        max_altitude)
+
 
 class LowAltitudeGroup(AbstractAltitudeGroup):
 
@@ -87,7 +86,6 @@ class DefaultAltitudeGroup(AbstractAltitudeGroup):
     def reference_altitude(self):
         return 500
 
-
 def get_altitude_group_from_altitudes(altitudes):
     s = set(altitudes)
     if s == set(altitudes_for_groups[0]):

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -51,12 +51,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         # Load one fold fit
         self._massif_name_to_one_fold_fit = {}
         for massif_name in self.massif_names:
-            if self.load_condition(massif_name):
+            # Load valid massif altitudes
+            massif_altitudes = self.get_massif_altitudes(massif_name)
+            if self.load_condition(massif_altitudes):
                 # Load dataset
-                dataset = studies.spatio_temporal_dataset(massif_name=massif_name)
-                dataset_without_zero = AbstractDataset.remove_zeros(dataset.observations,
-                                                                    dataset.coordinates)
-                old_fold_fit = OneFoldFit(massif_name, dataset_without_zero, model_classes, self.fit_method,
+                dataset = studies.spatio_temporal_dataset(massif_name=massif_name, massif_altitudes=massif_altitudes)
+                old_fold_fit = OneFoldFit(massif_name, dataset, model_classes, self.fit_method,
                                           self.temporal_covariate_for_fit,
                                           type(self.altitude_group),
                                           self.display_only_model_that_pass_anderson_test)
@@ -69,10 +69,25 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
     moment_names = ['moment', 'changes_for_moment', 'relative_changes_for_moment'][:]
     orders = [1, 2, None][2:]
 
-    def load_condition(self, massif_name):
-        altitudes_for_massif = [altitude for altitude, study in self.studies.altitude_to_study.items()
+    def get_massif_altitudes(self, massif_name):
+        valid_altitudes = [altitude for altitude, study in self.studies.altitude_to_study.items()
                                 if massif_name in study.study_massif_names]
-        return (self.altitude_group.reference_altitude in altitudes_for_massif) and (len(altitudes_for_massif) >= 2)
+        massif_altitudes = []
+        for altitude in valid_altitudes:
+            study = self.studies.altitude_to_study[altitude]
+            annual_maxima = study.massif_name_to_annual_maxima[massif_name]
+            percentage_of_non_zeros = 100 * np.count_nonzero(annual_maxima) / len(annual_maxima)
+            if percentage_of_non_zeros > 90:
+                massif_altitudes.append(altitude)
+            else:
+                print(massif_name, altitude, percentage_of_non_zeros)
+        return massif_altitudes
+
+    def load_condition(self, massif_altitudes):
+        # At least two altitudes for the estimated, including the reference altitude
+        reference_altitude_is_in_altitudes = (self.altitude_group.reference_altitude in massif_altitudes)
+        at_least_two_altitudes = (len(massif_altitudes) >= 2)
+        return reference_altitude_is_in_altitudes and at_least_two_altitudes
 
     @property
     def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
@@ -150,6 +165,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                       add_colorbar=add_colorbar,
                       max_abs_change=max_abs_change,
                       massif_name_to_text=massif_name_to_text,
+                      xlabel=self.altitude_group.xlabel,
                       )
 
     @property
@@ -250,6 +266,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
                       add_colorbar=self.add_colorbar,
                       max_abs_change=self._max_abs_for_shape,
                       fit_method=self.fit_method,
+                      xlabel=self.altitude_group.xlabel,
                       )
 
     def plot_altitude_for_the_peak(self):
@@ -374,24 +391,25 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
 
         nb_valid_massif_names = len(valid_massif_names)
         nbs = np.zeros(4)
-        relative_changes = [[], []]
+        relative_changes = []
         for one_fold in [one_fold for m, one_fold in self.massif_name_to_one_fold_fit.items()
                                    if m in valid_massif_names]:
-            if one_fold.change_for_reference_altitude == 0:
-                continue
             # Compute relative changes
-            relative_changes[0].append(one_fold.relative_change_for_reference_altitude)
-            if one_fold.is_significant:
-                relative_changes[1].append(one_fold.relative_change_for_reference_altitude)
+            relative_changes.append(one_fold.relative_change_in_return_level_for_reference_altitude)
+            # Compute nb of non stationary models
+            if one_fold.change_in_return_level_for_reference_altitude == 0:
+                continue
             # Compute nbs
-            idx = 0 if one_fold.change_for_reference_altitude < 0 else 2
+            idx = 0 if one_fold.change_in_return_level_for_reference_altitude < 0 else 2
             nbs[idx] += 1
             if one_fold.is_significant:
                 nbs[idx + 1] += 1
 
         percents = 100 * nbs / nb_valid_massif_names
-        mean_relative_changes = [np.mean(l) for l in relative_changes]
-        return [nb_valid_massif_names] + mean_relative_changes + list(percents)
+        mean_relative_change = np.mean(relative_changes)
+        median_relative_change = np.median(relative_changes)
+        print('mean', mean_relative_change, relative_changes)
+        return [nb_valid_massif_names, mean_relative_change, median_relative_change] + list(percents)
 
     def get_valid_names(self, massif_names):
         valid_massif_names = set(self.massif_name_to_one_fold_fit.keys())

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -28,7 +28,7 @@ from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
 class OneFoldFit(object):
     SIGNIFICANCE_LEVEL = 0.05
     best_estimator_minimizes_total_aic = False
-    return_period = 50
+    return_period = 100
 
     def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
                  fit_method=MarginFitMethod.extremes_fevd_mle, temporal_covariate_for_fit=None,
@@ -89,12 +89,12 @@ class OneFoldFit(object):
         return [self.get_moment(altitude, year, order) for altitude in altitudes]
 
     @property
-    def change_for_reference_altitude(self) -> float:
-        return self.changes_for_moment(altitudes=[self.altitude_plot])[0]
+    def change_in_return_level_for_reference_altitude(self) -> float:
+        return self.changes_for_moment(altitudes=[self.altitude_plot], order=None)[0]
 
     @property
-    def relative_change_for_reference_altitude(self) -> float:
-        return self.relative_changes_for_moment(altitudes=[self.altitude_plot])[0]
+    def relative_change_in_return_level_for_reference_altitude(self) -> float:
+        return self.relative_changes_for_moment(altitudes=[self.altitude_plot], order=None)[0]
 
     def changes_for_moment(self, altitudes, year=2019, nb_years=50, order=1):
         changes = []

projects/altitude_spatial_model/altitudes_fit/plot_histogram_altitude_studies.py

@@ -6,6 +6,7 @@ import matplotlib.pyplot as plt
 
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \
     AltitudesStudiesVisualizerForNonStationaryModels
+from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit
 
 
 def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: List[
@@ -32,8 +33,8 @@ def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: L
         x_shifted = [tick_middle + width * shift / 2 for shift in range(-3, 5, 2)]
         percentages = model_name_to_percentages[model_name]
         colors = ['white', 'yellow', 'orange', 'red']
-        labels = ['{} m (\% out of {} massifs)'.format(v.altitude_group.reference_altitude,
-                                                      len(v.get_valid_names(massif_names))) for v in visualizer_list]
+        labels = ['{} m - {} m (\% out of {} massifs)'.format(1000 * i, 1000 * (i+1),
+                                                      len(v.get_valid_names(massif_names))) for i, v in enumerate(visualizer_list)]
         for x, color, percentage, label in zip(x_shifted, colors, percentages, labels):
             ax.bar([x], [percentage], width=width, label=label,
                    linewidth=linewidth, edgecolor='black', color=color)
@@ -53,9 +54,12 @@ def plot_histogram_all_models_against_altitudes(massif_names, visualizer_list: L
     plt.close()
 
 
-def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list):
+def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list: List[
+    AltitudesStudiesVisualizerForNonStationaryModels]):
+    visualizer = visualizer_list[0]
+
     all_trends = [v.all_trends(massif_names) for v in visualizer_list]
-    nb_massifs, mean_relative_change, mean_significant_relative_change, *all_l = zip(*all_trends)
+    nb_massifs, mean_relative_changes, median_relative_changes, *all_l = zip(*all_trends)
 
     plt.close()
     ax = plt.gca()
@@ -87,17 +91,14 @@ def plot_histogram_all_trends_against_altitudes(massif_names, visualizer_list):
 
     # PLot mean relative change against altitude
     ax_twinx = ax.twinx()
-    colors = ['grey', 'black']
-    relative_changes = [mean_relative_change, mean_significant_relative_change]
-    labels = ['All massifs with a trend', 'Only massifs with a significant trend']
-    for color, y, label in zip(colors, relative_changes, labels):
-        print(x, y)
-        ax_twinx.plot(x, y, label=label, color=color)
+    ax_twinx.plot(x, mean_relative_changes, label='Mean relative change', color='black')
+    ax_twinx.plot(x, median_relative_changes, label='Median relative change', color='grey')
     ax_twinx.legend(loc='upper right', prop={'size': size})
+    ylabel = 'Relative change of {}-year return levels ({})'.format(OneFoldFit.return_period,
+                                                            visualizer.study.variable_unit)
+    ax_twinx.set_ylabel(ylabel, fontsize=legend_fontsize)
 
-    ax_twinx.set_ylabel('Mean relative change average on massif with trends', fontsize=legend_fontsize)
 
-    visualizer = visualizer_list[0]
     visualizer.plot_name = 'All trends'
     visualizer.show_or_save_to_file(add_classic_title=False, no_title=True)