diff --git a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py index e938e38649a685d5d6ed8607c520e160bd055707..2209091fb44a3788b6057f6ed1d3afbc38ff2a38 100644 --- a/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py +++ b/projects/altitude_spatial_model/altitudes_fit/main_altitudes_studies.py @@ -2,35 +2,34 @@ import datetime import time from typing import List -import matplotlib +import matplotlib as mpl +mpl.rcParams['text.usetex'] = True +mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] + +import matplotlib matplotlib.use('Agg') -from extreme_data.meteo_france_data.scm_models_data.safran.safran_max_snowf import SafranSnowfall2019, \ - SafranSnowfall2020 from projects.altitude_spatial_model.altitudes_fit.plots.plot_histogram_altitude_studies import \ plot_shoe_plot_changes_against_altitude, plot_histogram_all_trends_against_altitudes, \ - plot_shoe_plot_ratio_interval_size_against_altitude, plot_histogram_all_models_against_altitudes + plot_histogram_all_models_against_altitudes from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \ AbstractExtractEurocodeReturnLevel -import matplotlib as mpl + from extreme_fit.model.utils import set_seed_for_test from projects.altitude_spatial_model.altitudes_fit.plots.plot_coherence_curves import plot_coherence_curves -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 projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitude_group import altitudes_for_groups -from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.plot_total_aic import plot_individual_aic from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day, SafranSnowfall3Days, \ - SafranSnowfall5Days, SafranSnowfall7Days, SafranDateFirstSnowfall, SafranPrecipitation1Day, \ - SafranPrecipitation3Days, SafranSnowfallCenterOnDay1dayMeanRate, SafranSnowfallNotCenterOnDay1day + SafranSnowfall5Days, SafranSnowfall7Days from extreme_data.meteo_france_data.scm_models_data.utils import Season @@ -96,13 +95,11 @@ def plot_visualizer(massif_names, visualizer): # visualizer.studies.plot_maxima_time_series(massif_names) # visualizer.studies.plot_maxima_time_series(['Vanoise']) - # Plot the results for the model that minimizes the individual aic - plot_individual_aic(visualizer) - - # Plot the results for the model that minimizes the total aic - # plot_total_aic(model_classes, visualizer) - pass - + visualizer.plot_shape_map() + visualizer.plot_moments() + visualizer.plot_qqplots() + # for std in [True, False]: + # visualizer.studies.plot_mean_maxima_against_altitude(std=std) if __name__ == '__main__': main() diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py index 9c3b7f9fb436320c51657eb45897c3cb9052b89e..ffe76956002407e7415809c06e0daa94ad1084a4 100644 --- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py +++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/altitudes_studies_visualizer_for_non_stationary_models.py @@ -191,7 +191,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): if self.temporal_covariate_for_fit is AnomalyTemperatureTemporalCovariate else ' in {}' str_for_last_year = str_for_last_year.format(self.first_one_fold_fit.covariate_after, **d_temperature) moment = moment.replace('moment', '{}{}'.format(OneFoldFit.get_moment_str(order=order), str_for_last_year)) - plot_name = '{}{} '.format(OneFoldFit.folder_for_plots, moment) + plot_name = '{} '.format(moment) if 'change' in method_name: plot_name = plot_name.replace(str_for_last_year, '') @@ -264,7 +264,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): moment = ' '.join(method_name.split('_')) moment = moment.replace('moment', '{} in {}'.format(OneFoldFit.get_moment_str(order=order), OneFoldFit.last_year)) - plot_name = '{}Model {} annual maxima of {}'.format(OneFoldFit.folder_for_plots, moment, + plot_name = 'Model {} annual maxima of {}'.format(moment, SCM_STUDY_CLASS_TO_ABBREVIATION[self.studies.study_class]) ax.set_ylabel('{} ({})'.format(plot_name, self.study.variable_unit), fontsize=15) ax.set_xlabel('altitudes', fontsize=15) @@ -321,8 +321,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): negative_and_positive_values = (max(values) > 0) and (min(values) < 0) raise NotImplementedError self.plot_map(massif_name_to_value=massif_name_to_best_coef, - label='{}Coef/{} plot for {} {}'.format(OneFoldFit.folder_for_plots, - coef_name, + label='Coef/{} plot for {} {}'.format(coef_name, SCM_STUDY_CLASS_TO_ABBREVIATION[ type(self.study)], self.study.variable_unit), @@ -387,7 +386,7 @@ class AltitudesStudiesVisualizerForNonStationaryModels(StudyVisualizer): ax.set_ylabel('{} of {} in {} ({})'.format(ylabel, SCM_STUDY_CLASS_TO_ABBREVIATION[type(self.study)], massif_name.replace('_', ' '), self.study.variable_unit)) peak_year_folder = 'Peak year ' + ylabel - plot_name = '{}{}/Peak year for {}'.format(OneFoldFit.folder_for_plots, peak_year_folder, + plot_name = '{}/Peak year for {}'.format(peak_year_folder, massif_name.replace('_', '')) self.studies.show_or_save_to_file(plot_name=plot_name, show=self.show, no_title=True, tight_layout=True) plt.close() diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py index a0451c44757d65d4c15c898937e433a5cbc436ed..525e969e6fa9df844d846d0c8f817feb3fd818a6 100644 --- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py +++ b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/one_fold_fit.py @@ -43,7 +43,6 @@ from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observat class OneFoldFit(object): SIGNIFICANCE_LEVEL = 0.05 - best_estimator_minimizes_total_aic = False return_period = 100 quantile_level = 1 - (1 / return_period) nb_years = 60 @@ -75,21 +74,12 @@ class OneFoldFit(object): # Compute sorted estimators indirectly _ = self.has_at_least_one_valid_model - # Best estimator definition - self.best_estimator_class_for_total_aic = None - # Cached object - self._folder_to_goodness = {} - def fitted_linear_margin_estimator(self, model_class, dataset): return fitted_linear_margin_estimator_short(model_class=model_class, dataset=dataset, fit_method=self.fit_method, temporal_covariate_for_fit=self.temporal_covariate_for_fit) - @classproperty - def folder_for_plots(cls): - return 'Total aic/' if cls.best_estimator_minimizes_total_aic else '' - @classmethod def get_moment_str(cls, order): if order == 1: @@ -159,7 +149,7 @@ class OneFoldFit(object): s = ' between +${}^o\mathrm{C}$ and +${}^o\mathrm{C}$' if self.temporal_covariate_for_fit is AnomalyTemperatureTemporalCovariate \ else ' between {} and {}' s = s.format(self.covariate_before, self.covariate_after, - **d_temperature) + **d_temperature) return s def relative_changes_of_moment(self, altitudes, order=1): @@ -184,7 +174,7 @@ class OneFoldFit(object): for e in estimators: coordinate_values_for_the_fit = e.coordinates_for_nllh(Split.all) coordinate_values_for_the_result = [np.array([self.altitude_group.reference_altitude, c]) - for c in self._covariate_before_and_after] + for c in self._covariate_before_and_after] coordinate_values_to_check = list(coordinate_values_for_the_fit) + coordinate_values_for_the_result has_undefined_parameters = False for coordinate in coordinate_values_to_check: @@ -218,8 +208,6 @@ class OneFoldFit(object): if self.only_models_that_pass_goodness_of_fit_test: return [e for e in self.sorted_estimators if self.goodness_of_fit_test(e) - # and self.sensitivity_of_fit_test_top_maxima(e) - # and self.sensitivity_of_fit_test_last_years(e) ] else: if not self.remove_physically_implausible_models: @@ -236,16 +224,12 @@ class OneFoldFit(object): @property def best_estimator(self): - 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] + if self.has_at_least_one_valid_model: + best_estimator = self.sorted_estimators_with_stationary[0] + return best_estimator else: - # 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 object should not have been called because ' - 'has_at_least_one_valid_model={}'.format(self.has_at_least_one_valid_model)) + raise ValueError('This object should not have been called because ' + 'has_at_least_one_valid_model={}'.format(self.has_at_least_one_valid_model)) @property def best_margin_model(self): @@ -322,52 +306,6 @@ class OneFoldFit(object): # Bootstrap based significance return self.cached_results_from_bootstrap[0] - # @property - # def goodness_of_fit_anderson_test(self): - # if self.folder_for_plots in self._folder_to_goodness: - # return self._folder_to_goodness[self.folder_for_plots] - # else: - # estimator = self.best_estimator - # goodness_of_fit_anderson_test = self.goodness_of_fit_test(estimator) - # if not goodness_of_fit_anderson_test: - # print('{} with {} does not pass the anderson test for model {}'.format(self.massif_name, - # self.folder_for_plots, - # type(self.best_margin_model))) - # self._folder_to_goodness[self.folder_for_plots] = goodness_of_fit_anderson_test - # return goodness_of_fit_anderson_test - - def sensitivity_of_fit_test_top_maxima(self, estimator: LinearMarginEstimator): - # Build the dataset without the maxima for each altitude - new_dataset = AbstractDataset.remove_top_maxima(self.dataset.observations, - self.dataset.coordinates) - # Fit the new estimator - new_estimator = fitted_linear_margin_estimator_short(model_class=type(estimator.margin_model), - dataset=new_dataset, - fit_method=self.fit_method, - temporal_covariate_for_fit=self.temporal_covariate_for_fit) - # Compare sign of change - has_not_opposite_sign = self.sign_of_change(estimator.function_from_fit) * self.sign_of_change( - new_estimator.function_from_fit) >= 0 - return has_not_opposite_sign - - def sensitivity_of_fit_test_last_years(self, estimator: LinearMarginEstimator): - # Build the dataset without the maxima for each altitude - new_dataset = AbstractDataset.remove_last_maxima(self.dataset.observations, - self.dataset.coordinates, - nb_years=10) - # Fit the new estimator - model_class = type(estimator.margin_model) - new_estimator = fitted_linear_margin_estimator_short(model_class=model_class, - dataset=new_dataset, - fit_method=self.fit_method, - temporal_covariate_for_fit=self.temporal_covariate_for_fit) - # Compare sign of change - has_not_opposite_sign = self.sign_of_change(estimator.function_from_fit) * self.sign_of_change( - new_estimator.function_from_fit) >= 0 - # if not has_not_opposite_sign: - # print('Last years', self.massif_name, model_class, self.sign_of_change(estimator), self.sign_of_change(new_estimator)) - return has_not_opposite_sign - def sign_of_change(self, function_from_fit): return_levels = [] for temporal_covariate in self._covariate_before_and_after: @@ -445,7 +383,6 @@ class OneFoldFit(object): return is_significant, altitude_and_year_to_return_level_mean_estimate, altitude_and_year_to_return_level_confidence_interval - @property def bootstrap_fitted_functions_from_fit(self): print('nb of bootstrap for confidence interval=', AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP) diff --git a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py b/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py deleted file mode 100644 index 2a50df32b9783491ae211e96f6930dfd589f572f..0000000000000000000000000000000000000000 --- a/projects/altitude_spatial_model/altitudes_fit/one_fold_analysis/plot_total_aic.py +++ /dev/null @@ -1,90 +0,0 @@ -import matplotlib as mpl -import matplotlib.pyplot as plt -import numpy as np - -from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.altitudes_studies_visualizer_for_non_stationary_models import \ - AltitudesStudiesVisualizerForNonStationaryModels - -mpl.rcParams['text.usetex'] = True -mpl.rcParams['text.latex.preamble'] = [r'\usepackage{amsmath}'] - -from extreme_data.meteo_france_data.scm_models_data.visualization.utils import create_adjusted_axes -from projects.altitude_spatial_model.altitudes_fit.one_fold_analysis.one_fold_fit import OneFoldFit -from projects.exceeding_snow_loads.utils import dpi_paper1_figure - - -def plots(visualizer: AltitudesStudiesVisualizerForNonStationaryModels): - visualizer.plot_shape_map() - visualizer.plot_moments() - visualizer.plot_qqplots() - # for std in [True, False]: - # visualizer.studies.plot_mean_maxima_against_altitude(std=std) - - -def plot_individual_aic(visualizer): - OneFoldFit.best_estimator_minimizes_total_aic = False - plots(visualizer) - - -def plot_total_aic(model_classes, visualizer): - # Compute the mean AIC for each model_class - OneFoldFit.best_estimator_minimizes_total_aic = True - model_class_to_total_aic = {model_class: 0 for model_class in model_classes} - model_class_to_name_str = {} - # model_class_to_aic_scores = {model_class: [] for model_class in model_classes} - # model_class_to_n = {model_class: [] for model_class in model_classes} - for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values(): - for model_class, estimator in one_fold_fit.model_class_to_estimator_with_finite_aic.items(): - model_class_to_total_aic[model_class] += estimator.aic() - model_class_to_name_str[model_class] = estimator.margin_model.name_str - # model_class_to_aic_scores[model_class].append(estimator.aic()) - # model_class_to_n[model_class].append(estimator.n()) - # model_class_to_mean_aic_score = {model_class: np.array(aic_scores).sum() / np.array(model_class_to_n[model_class]).sum() - # for model_class, aic_scores in model_class_to_aic_scores.items()} - # print(model_class_to_mean_aic_score) - sorted_model_class = sorted(model_classes, key=lambda m: model_class_to_total_aic[m]) - sorted_scores = [model_class_to_total_aic[model_class] for model_class in sorted_model_class] - sorted_labels = [model_class_to_name_str[model_class] for model_class in sorted_model_class] - print(sorted_model_class) - print(sorted_scores) - print(sorted_labels) - best_model_class_for_total_aic = sorted_model_class[0] - for one_fold_fit in visualizer.massif_name_to_one_fold_fit.values(): - one_fold_fit.best_estimator_class_for_total_aic = best_model_class_for_total_aic - # Plot the ranking of the model based on their total aic - plot_total_aic_repartition(visualizer, sorted_labels, sorted_scores) - # Plot the results for the model that minimizes the mean aic - plots(visualizer) - - -def plot_total_aic_repartition(visualizer, labels, scores): - """ - Plot a single trend curves - :return: - """ - scores = np.array(scores) - ax = create_adjusted_axes(1, 1) - - # parameters - width = 3 - size = 30 - legend_fontsize = 30 - labelsize = 10 - linewidth = 3 - x = [2 * width * (i + 1) for i in range(len(labels))] - ax.bar(x, scores, width=width, color='grey', edgecolor='grey', label='Total Aic', - linewidth=linewidth) - ax.legend(loc='upper right', prop={'size': size}) - ax.set_ylabel(' Total AIC score \n ' - 'i.e. sum of AIC score for all massifs ', fontsize=legend_fontsize) - ax.set_xlabel('Models', fontsize=legend_fontsize) - ax.tick_params(axis='both', which='major', labelsize=labelsize) - ax.set_xticks(x) - ax.set_ylim([scores.min(), scores.max()]) - ax.yaxis.grid() - ax.set_xticklabels(labels) - - # Save plot - visualizer.plot_name = 'Total AIC ranking' - visualizer.show_or_save_to_file(no_title=True, dpi=dpi_paper1_figure) - plt.close() diff --git a/projects/altitude_spatial_model/altitudes_fit/plot_mean_and_std.py b/projects/altitude_spatial_model/altitudes_fit/plot_mean_and_std.py deleted file mode 100644 index 20bc2fc7c6a6093833aad9e536ec0d4db2d98ea3..0000000000000000000000000000000000000000 --- a/projects/altitude_spatial_model/altitudes_fit/plot_mean_and_std.py +++ /dev/null @@ -1,8 +0,0 @@ -from extreme_data.meteo_france_data.scm_models_data.safran.safran import SafranSnowfall1Day -from projects.altitude_spatial_model.altitudes_fit.altitudes_studies import AltitudesStudies - -if __name__ == '__main__': - studies = AltitudesStudies(study_class=SafranSnowfall1Day, - altitudes=[600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600]) - for std in [True, False]: - studies.plot_mean_maxima_against_altitude(std=std) \ No newline at end of file diff --git a/projects/archive/gap_between_my_safran2019_and_safran_2019/__init__.py b/projects/archive/gap_between_my_safran2019_and_safran_2019/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/projects/altitude_spatial_model/altitudes_fit/main_gap_altitudes_studies.py b/projects/archive/gap_between_my_safran2019_and_safran_2019/main_gap_altitudes_studies.py similarity index 100% rename from projects/altitude_spatial_model/altitudes_fit/main_gap_altitudes_studies.py rename to projects/archive/gap_between_my_safran2019_and_safran_2019/main_gap_altitudes_studies.py