diff --git a/experiment/meteo_france_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py b/experiment/meteo_france_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
index 730df279670e4d0585cb33716b5bc9b0a52cf4c7..80b4737ee43ab360c4e9ddb522e8e6c0005342a3 100644
--- a/experiment/meteo_france_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
+++ b/experiment/meteo_france_data/visualization/hypercube_visualization/abstract_hypercube_visualizer.py
@@ -19,8 +19,10 @@ class AbstractHypercubeVisualizer(object):
def __init__(self, tuple_to_study_visualizer: Dict[Tuple, StudyVisualizer],
trend_test_class,
nb_data_reduced_for_speed=False,
- save_to_file=False):
+ save_to_file=False,
+ nb_top_likelihood_values=1):
self.nb_data_for_fast_mode = 7 if nb_data_reduced_for_speed else None
+ self.nb_top_likelihood_values = nb_top_likelihood_values
self.save_to_file = save_to_file
self.trend_test_class = trend_test_class
self.tuple_to_study_visualizer = tuple_to_study_visualizer # type: Dict[Tuple, StudyVisualizer]
@@ -33,7 +35,7 @@ class AbstractHypercubeVisualizer(object):
@cached_property
def starting_years(self):
- starting_years = self.study_visualizer.starting_years
+ starting_years = self.study_visualizer.starting_years[:]
if self.nb_data_for_fast_mode is not None:
starting_years = starting_years[:self.nb_data_for_fast_mode]
return starting_years
@@ -44,7 +46,8 @@ class AbstractHypercubeVisualizer(object):
@cached_property
def df_trends_spatio_temporal(self):
return [study_visualizer.df_trend_spatio_temporal(self.trend_test_class, self.starting_years,
- self.nb_data_for_fast_mode)
+ self.nb_data_for_fast_mode,
+ self.nb_top_likelihood_values)
for study_visualizer in self.tuple_to_study_visualizer.values()]
@cached_property
diff --git a/experiment/meteo_france_data/visualization/hypercube_visualization/main_hypercube_visualization.py b/experiment/meteo_france_data/visualization/hypercube_visualization/main_hypercube_visualization.py
index 57500696093b72f373a5aa589677212a54a06ef4..fcca15b5827cd205364b8c3b65cd45492f231d97 100644
--- a/experiment/meteo_france_data/visualization/hypercube_visualization/main_hypercube_visualization.py
+++ b/experiment/meteo_france_data/visualization/hypercube_visualization/main_hypercube_visualization.py
@@ -85,7 +85,8 @@ def fast_altitude_year_hypercube():
altitudes=altitudes)]
altitude_to_visualizer = OrderedDict(zip(altitudes, visualizers))
visualizer = Altitude_Hypercube_Year_Visualizer(altitude_to_visualizer, save_to_file=save_to_file,
- trend_test_class=trend_test_class, nb_data_reduced_for_speed=nb_data_reduced_for_speed)
+ trend_test_class=trend_test_class,
+ nb_data_reduced_for_speed=nb_data_reduced_for_speed)
visualizer.visualize_year_trend_test()
# visualizer.visualize_altitude_trend_test()
# visualizer.visualize_massif_trend_test()
diff --git a/experiment/meteo_france_data/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_data/visualization/study_visualization/study_visualizer.py
index a91ea95408ed81ed6b89e656b9fd6bbd4132c847..a50338a70dd3872307260aa12f9be7521d94803c 100644
--- a/experiment/meteo_france_data/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_data/visualization/study_visualization/study_visualizer.py
@@ -385,7 +385,8 @@ class StudyVisualizer(object):
massif_name_to_df_trend_type[massif_name] = df
return massif_name_to_df_trend_type
- def df_trend_spatio_temporal(self, trend_test_class, starting_years, nb_massif_for_fast_mode=None):
+ def df_trend_spatio_temporal(self, trend_test_class, starting_years, nb_massif_for_fast_mode=None,
+ nb_top_likelihood_values=1):
"""
Index are the massif
Columns are the starting year
@@ -411,9 +412,13 @@ class StudyVisualizer(object):
for starting_year in starting_years]
# Keep only the most likely starting year
# (set all the other data to np.nan so that they will not be taken into account in mean function)
- best_idx = np.argmax(trend_test_res, axis=0)[-1]
- trend_test_res = [(a, b) if i == best_idx else (np.nan, np.nan)
- for i, (a, b, _) in enumerate(trend_test_res)]
+ best_idx = list(np.argmax(trend_test_res, axis=0))[2]
+ # print(best_idx, trend_test_res)
+ best_idxs = [best_idx]
+ # todo: by doing a sorting on the deviance, I could get the nb_top_likelihood_values values
+ # best_idxs = list(np.argmax(trend_test_res, axis=0))[-nb_top_likelihood_values:]
+ trend_test_res = [(a, b) if i in best_idxs else (np.nan, np.nan)
+ for i, (a, b, *_) in enumerate(trend_test_res)]
massif_name_to_trend_res[massif_name] = list(zip(*trend_test_res))
nb_res = len(list(massif_name_to_trend_res.values())[0])
assert nb_res == 2
@@ -426,7 +431,7 @@ class StudyVisualizer(object):
def compute_gev_change_point_test_result(smooth_maxima, starting_year, trend_test_class, years):
trend_test = trend_test_class(years, smooth_maxima, starting_year) # type: AbstractGevChangePointTest
assert isinstance(trend_test, AbstractGevChangePointTest)
- return trend_test.test_trend_type, trend_test.test_trend_strength, trend_test.non_stationary_nllh
+ return trend_test.test_trend_type, trend_test.test_trend_strength, trend_test.non_stationary_nllh, trend_test.non_stationary_deviance, trend_test.stationary_deviance
@staticmethod
def compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years):
diff --git a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
index c326ac2d146e6ec6e94e810ea011c5126f7485ac..9348e9baefc1736d7bca4b00778debb032d078ae 100644
--- a/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
@@ -57,8 +57,7 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
@property
def likelihood_ratio(self):
- return 2 * (self.non_stationary_estimator.result_from_fit.deviance -
- self.stationary_estimator.result_from_fit.deviance)
+ return 2 * (self.non_stationary_deviance - self.stationary_deviance)
@property
def non_stationary_nllh(self):
@@ -67,6 +66,14 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
else:
return self.non_stationary_estimator.result_from_fit.nllh
+ @property
+ def stationary_deviance(self):
+ return self.stationary_estimator.result_from_fit.deviance
+
+ @property
+ def non_stationary_deviance(self):
+ return self.non_stationary_estimator.result_from_fit.deviance
+
@property
def is_significant(self) -> bool:
return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=1)