[HYPERCUBE VISUALIZER] add nb_top_likelihood values as argument

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

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()

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):

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)