[contrasting] consider also the stationary model for selection. consider a...

[contrasting] consider also the stationary model for selection. consider a level of altitude only if at least 2 altitudes are present & the reference altitude is present.

extreme_fit/model/margin_model/polynomial_margin_model/gev_altitudinal_models.py

@@ -43,7 +43,11 @@ class AbstractAltitudinalModel(AbstractSpatioTemporalPolynomialModel):
                 s += '^2'
         if self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates) in ['1', '2']:
             s += '\\sigma_t'
-            if self.dim_to_str_number(GevParams.LOC, self.coordinates.idx_temporal_coordinates) == '2':
+            if self.dim_to_str_number(GevParams.SCALE, self.coordinates.idx_temporal_coordinates) == '2':
+                s += '^2'
+        if self.dim_to_str_number(GevParams.SHAPE, self.coordinates.idx_temporal_coordinates) in ['1', '2']:
+            s += '\\zeta_t'
+            if self.dim_to_str_number(GevParams.SHAPE, self.coordinates.idx_temporal_coordinates) == '2':
                 s += '^2'
         if len(s) == 0:
             s = '0'

extreme_fit/model/margin_model/polynomial_margin_model/utils.py

@@ -105,18 +105,22 @@ ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS = [
     StationaryAltitudinal,
     AltitudinalShapeConstantTimeLocationLinear,
     AltitudinalShapeConstantTimeScaleLinear,
-    # AltitudinalShapeConstantTimeShapeLinear,
-    # AltitudinalShapeConstantTimeLocShapeLinear,
+
     AltitudinalShapeConstantTimeLocScaleLinear,
-    # AltitudinalShapeConstantTimeScaleShapeLinear,
-    # AltitudinalShapeConstantTimeLocScaleShapeLinear,
+
     # With a linear shape for the altitude (8 models)
     AltitudinalShapeLinearTimeStationary,
     AltitudinalShapeLinearTimeLocationLinear,
     AltitudinalShapeLinearTimeScaleLinear,
+
+    AltitudinalShapeLinearTimeLocScaleLinear,
+
+    # AltitudinalShapeConstantTimeShapeLinear,
+    # AltitudinalShapeConstantTimeLocShapeLinear,
+    # AltitudinalShapeConstantTimeScaleShapeLinear,
+    # AltitudinalShapeConstantTimeLocScaleShapeLinear,
     # AltitudinalShapeLinearTimeShapeLinear,
     # AltitudinalShapeLinearTimeLocShapeLinear,
-    AltitudinalShapeLinearTimeLocScaleLinear,
     # AltitudinalShapeLinearTimeScaleShapeLinear,
     # AltitudinalShapeLinearTimeLocScaleShapeLinear,
 

extreme_fit/model/result_from_model_fit/result_from_extremes/abstract_result_from_extremes.py

@@ -62,6 +62,7 @@ class AbstractResultFromExtremes(AbstractResultFromModelFit):
         if self.param_name_to_dim:
             d = {GevParams.greek_letter_from_param_name(param_name) + '1': r.c(transformed_temporal_covariate) for
                  param_name in self.param_name_to_dim.keys()}
+            print(d)
             kwargs = {
                 "vals": r.list(**d
                                )

projects/altitude_spatial_model/altitudes_fit/altitudes_studies.py

@@ -139,9 +139,11 @@ class AltitudesStudies(object):
         plot_name = 'Annual maxima of {} in {}'.format(SCM_STUDY_CLASS_TO_ABBREVIATION[self.study_class],
                                                        massif_name.replace('_', ' '))
         ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15)
+        plot_name = 'time series/' + plot_name
         ax.set_xlabel('years', fontsize=15)
         self.show_or_save_to_file(plot_name=plot_name, show=show, no_title=True)
         ax.clear()
+        plt.close()
 
     def plot_mean_maxima_against_altitude(self, massif_names=None, show=False, std=False, change=False):
         ax = plt.gca()

projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py

@@ -47,7 +47,7 @@ def main():
     # massif_names = ['Mercantour', 'Vercors', 'Ubaye']
     seasons = [Season.annual, Season.winter, Season.spring, Season.automn][:1]
 
-    main_loop(altitudes_for_groups[:1], massif_names, seasons, study_classes)
+    main_loop(altitudes_for_groups[:], massif_names, seasons, study_classes)
 
 
 def main_loop(altitudes_list, massif_names, seasons, study_classes):
@@ -68,6 +68,7 @@ def main_loop(altitudes_list, massif_names, seasons, study_classes):
             time.sleep(2)
 
 
+
 def load_visualizer_list(season, study_class, altitudes_list, massif_names):
     model_classes = ALTITUDINAL_GEV_MODELS_BASED_ON_POINTWISE_ANALYSIS
     visualizer_list = []
@@ -107,7 +108,7 @@ def plots(massif_names, season, visualizer):
     print('inner loop', season, type(studies.study))
 
     # Plot time series
-    # studies.plot_maxima_time_series(massif_names=massif_names)
+    studies.plot_maxima_time_series(massif_names=massif_names)
 
     # Plot moments
     # for std in [True, False][:]:

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitude_group.py

@@ -5,8 +5,10 @@ altitudes_for_groups = [
     [300, 600, 900],
     [1200, 1500, 1800],
     [2100, 2400, 2700],
-    [3000, 3300, 3600]
+    [3000, 3300, 3600, 3900]
 ]
+
+
 # altitudes_for_groups = [
 #     [900, 1200, 1500],
 #     [1800, 2100, 2400],
@@ -63,6 +65,7 @@ class HighAltitudeGroup(AbstractAltitudeGroup):
     def reference_altitude(self):
         return 2400
 
+
 class VeyHighAltitudeGroup(AbstractAltitudeGroup):
 
     @property

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py

@@ -21,7 +21,7 @@ from extreme_fit.model.margin_model.polynomial_margin_model.spatio_temporal_poly
 from extreme_fit.model.margin_model.utils import MarginFitMethod
 from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import \
-    get_altitude_group_from_altitudes, HighAltitudeGroup
+    get_altitude_group_from_altitudes, HighAltitudeGroup, VeyHighAltitudeGroup
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import \
     OneFoldFit
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
@@ -50,12 +50,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         # Load one fold fit
         self._massif_name_to_one_fold_fit = {}
         for massif_name in self.massif_names:
-            if any([massif_name in study.study_massif_names for study in self.studies.altitude_to_study.values()]):
+            if self.load_condition(massif_name):
                 # 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, model_classes, self.fit_method,
+                dataset_without_zero = AbstractDataset.remove_zeros(dataset.observations,
+                                                                    dataset.coordinates)
+                old_fold_fit = OneFoldFit(massif_name, dataset_without_zero, model_classes, self.fit_method,
                                           self.temporal_covariate_for_fit,
                                           type(self.altitude_group),
                                           self.display_only_model_that_pass_anderson_test)
@@ -66,7 +66,12 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
         self._max_abs_for_shape = None
 
     moment_names = ['moment', 'changes_for_moment', 'relative_changes_for_moment'][:]
-    orders = [1, 2, None][:]
+    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()
+                                if massif_name in study.study_massif_names]
+        return (self.altitude_group.reference_altitude in altitudes_for_massif) and (len(altitudes_for_massif) >= 2)
 
     @property
     def massif_name_to_one_fold_fit(self) -> Dict[str, OneFoldFit]:
@@ -132,7 +137,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer):
 
     @property
     def add_colorbar(self):
-        return isinstance(self.altitude_group, HighAltitudeGroup)
+        return isinstance(self.altitude_group, VeyHighAltitudeGroup)
 
     def plot_against_years(self, method_name, order):
         ax = plt.gca()

projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py

@@ -15,6 +15,10 @@ from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_m
 from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pariwise_analysis.gev_with_linear_shape_wrt_altitude import \
     AltitudinalShapeLinearTimeStationary
 from extreme_fit.model.margin_model.utils import MarginFitMethod
+from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
+    ConfidenceIntervalMethodFromExtremes
+from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
+    EurocodeConfidenceIntervalFromExtremes
 from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import AbstractAltitudeGroup, \
     DefaultAltitudeGroup
 from root_utils import classproperty
@@ -110,6 +114,17 @@ class OneFoldFit(object):
         sorted_estimators = sorted([estimator for estimator in estimators], key=lambda e: e.aic())
         return sorted_estimators
 
+    @cached_property
+    def sorted_estimators_with_stationary(self):
+        if self.only_models_that_pass_anderson_test:
+            return [e for e in self.sorted_estimators if self.goodness_of_fit_test(e)]
+        else:
+            return self._sorted_estimators_without_stationary
+
+    @property
+    def has_at_least_one_valid_model(self):
+        return len(self.sorted_estimators_with_stationary) > 0
+
     @cached_property
     def _sorted_estimators_without_stationary(self):
         return [e for e in self.sorted_estimators if not isinstance(e.margin_model, StationaryAltitudinal)]
@@ -134,8 +149,13 @@ class OneFoldFit(object):
         if self.best_estimator_minimizes_total_aic and self.best_estimator_class_for_total_aic is not None:
             return self.model_class_to_estimator[self.best_estimator_class_for_total_aic]
         else:
-            if self.has_at_least_one_valid_non_stationary_model:
-                best_estimator = self.sorted_estimators_without_stationary[0]
+            # Without stationary
+            # if self.has_at_least_one_valid_non_stationary_model:
+            #     best_estimator = self.sorted_estimators_without_stationary[0]
+            #     return best_estimator
+            # With stationary
+            if self.has_at_least_one_valid_model:
+                best_estimator = self.sorted_estimators_with_stationary[0]
                 return best_estimator
             else:
                 raise ValueError('This should not happen')
@@ -169,7 +189,10 @@ class OneFoldFit(object):
     def best_name(self):
         name = self.best_estimator.margin_model.name_str
         latex_command = 'textbf' if self.is_significant else 'textrm'
-        return '$\\' + latex_command + '{' + name + '}$'
+        best_name = '$\\' + latex_command + '{' + name + '}$'
+        if self.is_significant:
+            best_name = '\\underline{' + best_name + '}'
+        return best_name
 
     # Significant
 
@@ -181,6 +204,8 @@ class OneFoldFit(object):
             return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinalOnlyScale]
         elif isinstance(self.best_estimator.margin_model, AltitudinalShapeLinearTimeStationary):
             return self.model_class_to_estimator_with_finite_aic[AltitudinalShapeLinearTimeStationary]
+        elif isinstance(self.best_estimator.margin_model, AltitudinalShapeLinearTimeStationary):
+            return self.model_class_to_estimator_with_finite_aic[AltitudinalShapeLinearTimeStationary]
         else:
             return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinal]
 
@@ -194,7 +219,7 @@ class OneFoldFit(object):
 
     @property
     def is_significant(self) -> bool:
-        stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal]
+        stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal, AltitudinalShapeLinearTimeStationary]
         if any([isinstance(self.best_estimator.margin_model, c)
                 for c in stationary_model_classes]):
             return False
@@ -233,3 +258,8 @@ class OneFoldFit(object):
             empirical_quantiles.append(maximum_standardized)
         empirical_quantiles = sorted(empirical_quantiles)
         return empirical_quantiles
+
+    # def best_confidence_interval(self):
+    #     EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(self.best_estimator,
+    #                                                                    ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
+    #                                                                    temporal_covariate=np.array([2019, self.altitude_plot]),)
\ No newline at end of file