diff --git a/experiment/meteo_france_SCM_study/main_visualize.py b/experiment/meteo_france_SCM_study/main_visualize.py index 9110ac96520621533c500d425dd78a9d399fb30e..0a96d04258a9b8c10f765211cb0b06999376383c 100644 --- a/experiment/meteo_france_SCM_study/main_visualize.py +++ b/experiment/meteo_france_SCM_study/main_visualize.py @@ -19,7 +19,7 @@ def study_iterator(study_class, only_first_one=False, both_altitude=False, verbo if verbose: print('Loading studies....') for nb_day in nb_days: - for alti in AbstractStudy.ALTITUDES[::1]: + for alti in AbstractStudy.ALTITUDES[::-1]: if verbose: print('alti: {}, nb_day: {}'.format(alti, nb_day)) study = study_class(alti, nb_day) if is_safran_study else study_class(alti) @@ -36,11 +36,12 @@ def extended_visualization(): for study_class in SCM_EXTENDED_STUDIES[1:2]: for study in study_iterator(study_class, only_first_one=True): study_visualizer = StudyVisualizer(study) - study_visualizer.visualize_all_kde_graphs() + # study_visualizer.visualize_all_kde_graphs() + study_visualizer.visualize_experimental_law() def normal_visualization(): - for study_class in SCM_STUDIES[:1]: + for study_class in SCM_STUDIES[1:2]: for study in study_iterator(study_class, only_first_one=True): study_visualizer = StudyVisualizer(study) # study_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][0]) @@ -62,6 +63,6 @@ def complete_analysis(only_first_one=False): if __name__ == '__main__': - # normal_visualization() - extended_visualization() + normal_visualization() + # extended_visualization() # complete_analysis() diff --git a/experiment/meteo_france_SCM_study/safran/safran_visualizer.py b/experiment/meteo_france_SCM_study/safran/safran_visualizer.py index 113558758e4a6f27fcfeb0aad0e5b971d03176aa..b147faa9f99d6ae53201616aaa7689b0df4af52f 100644 --- a/experiment/meteo_france_SCM_study/safran/safran_visualizer.py +++ b/experiment/meteo_france_SCM_study/safran/safran_visualizer.py @@ -43,6 +43,10 @@ class StudyVisualizer(object): def dataset(self): return AbstractDataset(self.observations, self.coordinates) + def visualize_experimental_law(self): + plot_name = ' experimental law' + self.show_or_save_to_file(plot_name) + def visualize_all_kde_graphs(self): massif_names = self.study.safran_massif_names nb_columns = 5 diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py index 09866a6abe83064a6cd7fac38c569827cec40432..cc1ac24b4083b44da91960d0050974bdaf366223 100644 --- a/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py +++ b/extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py @@ -121,6 +121,7 @@ class AbstractMarginFunction(object): grid = [] for i, xi in enumerate(linspace): gev_param = self.get_gev_params(np.array([xi])) + assert not gev_param.has_undefined_parameters, 'This case needs to be handled during display' grid.append(gev_param.summary_dict) grid = {gev_param: [g[gev_param] for g in grid] for gev_param in GevParams.SUMMARY_NAMES} return grid, linspace diff --git a/extreme_estimator/margin_fits/extreme_params.py b/extreme_estimator/margin_fits/extreme_params.py index 01a39d3fb04a8ba143073bf980c4011546ea7da9..a20ace4f1493ebe88aa2ab28b29524b8fae33552 100644 --- a/extreme_estimator/margin_fits/extreme_params.py +++ b/extreme_estimator/margin_fits/extreme_params.py @@ -1,4 +1,5 @@ from abc import ABC +import numpy as np from extreme_estimator.margin_fits.abstract_params import AbstractParams @@ -13,4 +14,8 @@ class ExtremeParams(AbstractParams, ABC): self.location = loc self.scale = scale self.shape = shape - assert self.scale > 0 + # Scale cannot be negative + # (sometimes it happens, when we want to find a quantile for every point of a 2D map + # then it can happen that a corner point that was not used for fitting correspond to a negative scale, + # in the case we set all the parameters as equal to np.nan, and we will not display those points) + self.has_undefined_parameters = self.scale <= 0 \ No newline at end of file diff --git a/extreme_estimator/margin_fits/gev/gev_params.py b/extreme_estimator/margin_fits/gev/gev_params.py index b4d9cbc2db64732bd397bc943c64bfdc407a053d..17081517fd2aae9ebef8a28e11cce8613d8d689e 100644 --- a/extreme_estimator/margin_fits/gev/gev_params.py +++ b/extreme_estimator/margin_fits/gev/gev_params.py @@ -1,5 +1,6 @@ from extreme_estimator.extreme_models.utils import r from extreme_estimator.margin_fits.extreme_params import ExtremeParams +import numpy as np class GevParams(ExtremeParams): @@ -14,8 +15,14 @@ class GevParams(ExtremeParams): self.block_size = block_size def quantile(self, p) -> float: - return r.qgev(p, self.location, self.scale, self.shape)[0] + if self.has_undefined_parameters: + return np.nan + else: + return r.qgev(p, self.location, self.scale, self.shape)[0] @property def param_values(self): - return [self.location, self.scale, self.shape] \ No newline at end of file + if self.has_undefined_parameters: + return [np.nan for _ in range(3)] + else: + return [self.location, self.scale, self.shape] \ No newline at end of file diff --git a/extreme_estimator/margin_fits/plot/create_shifted_cmap.py b/extreme_estimator/margin_fits/plot/create_shifted_cmap.py index b0fe67fc040e6405d4c710ad71d2ca2882557677..db147f9cb422b43a8e343e7f57a63a082fe1fe85 100644 --- a/extreme_estimator/margin_fits/plot/create_shifted_cmap.py +++ b/extreme_estimator/margin_fits/plot/create_shifted_cmap.py @@ -24,6 +24,8 @@ def plot_extreme_param(ax, gev_param_name, values): midpoint = 1.0 elif vmin > 0 and vmax > 0: midpoint = 0.0 + else: + raise ValueError('Unexpected values: vmin={}, vmax={}'.format(vmin, vmax)) cmap = [plt.cm.coolwarm, plt.cm.bwr, plt.cm.seismic][1] shifted_cmap = shiftedColorMap(cmap, midpoint=midpoint, name='shifted') norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax) @@ -47,10 +49,9 @@ def get_color_rbga_shifted(ax, gev_param_name, values): def imshow_shifted(ax, gev_param_name, values, x, y): - norm, shifted_cmap = plot_extreme_param(ax, gev_param_name, values) + masked_array = np.ma.masked_where(np.isnan(values), values) + norm, shifted_cmap = plot_extreme_param(ax, gev_param_name, masked_array) shifted_cmap.set_bad(color='white') - - masked_array = values if gev_param_name != ExtremeParams.SHAPE: epsilon = 1e-2 * (np.max(values) - np.min(values)) value = np.min(values) diff --git a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_params.py b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_params.py index dab5ca5508814d9a6269addf2ec02fc24d7ed470..39d3da709d0e2ab155cdaf7f21ab601e5094bc76 100644 --- a/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_params.py +++ b/test/test_extreme_estimator/test_margin_fits/test_gev/test_gev_params.py @@ -14,6 +14,12 @@ class TestGevParams(unittest.TestCase): for quantile_name, p in gev_params.quantile_name_to_p.items(): self.assertAlmostEqual(- 1 / np.log(p), quantile_dict[quantile_name]) + def test_negative_scale(self): + gev_params = GevParams(loc=1.0, shape=1.0, scale=-1.0) + for p in [0.1, 0.5, 0.9]: + q = gev_params.quantile(p) + self.assertTrue(np.isnan(q)) + if __name__ == '__main__': unittest.main()