diff --git a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
index 879e970381965110eb6699277bd6c617f47e731e..df0115f1ba7e22452ff0aa64b0de14ad77bf6439 100644
--- a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/altitude_hypercube_visualizer.py
@@ -4,7 +4,7 @@ import pandas as pd
from experiment.meteo_france_SCM_study.visualization.hypercube_visualization.abstract_hypercube_visualizer import \
AbstractHypercubeVisualizer
from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
-from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import AbstractTrendTest
+from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
from utils import get_display_name_from_object_type
@@ -27,7 +27,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
trend_type_to_series = {}
for trend_type in self.trend_types:
# Reduce df_bool df to a serie s_trend_type_percentage
- df_bool = self.df_hypercube_trend_type.isin(AbstractTrendTest.get_trend_types(trend_type))
+ df_bool = self.df_hypercube_trend_type.isin(AbstractUnivariateTest.get_trend_types(trend_type))
s_trend_type_percentage = reduction_function(df_bool)
assert isinstance(s_trend_type_percentage, pd.Series)
assert not isinstance(s_trend_type_percentage.index, pd.MultiIndex)
@@ -50,16 +50,16 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
def get_title_plot(self, xlabel, ax_idx=None):
labels = ['altitudes', 'starting years', 'massifs']
assert xlabel in labels, xlabel
- labels.remove(xlabel)
- common_txt = 'averaged on {} & {}'.format(*labels)
- if ax_idx is None:
- return common_txt
+ if ax_idx == 1:
+ return '% of change per year for the parameter value'
+ elif ax_idx == 0:
+ return '% of trend type'
else:
- assert ax_idx in [0, 1]
- if ax_idx == 0:
- return '% of trend type '
- else:
- return '% of change per year for the parameter value'
+ labels.remove(xlabel)
+ if xlabel != 'starting years':
+ labels.remove('starting years')
+ common_txt = 'averaged on {}'.format(' & '.join(labels))
+ return common_txt
def visualize_trend_test_evolution(self, reduction_function, xlabel, xlabel_values, axes=None, marker='o',
subtitle=''):
@@ -74,11 +74,11 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
ax.plot(xlabel_values, percentages_values, style + marker, label=trend_type)
if i == 0:
- # Plot the total value of significative values
- significative_values = trend_type_to_series[AbstractTrendTest.SIGNIFICATIVE_NEGATIVE_TREND][i] \
- + trend_type_to_series[AbstractTrendTest.SIGNIFICATIVE_POSITIVE_TREND][i]
- ax.plot(xlabel_values, significative_values, 'y-' + marker,
- label=AbstractTrendTest.SIGNIFICATIVE + ' trends')
+ # # Plot the total value of significative values
+ # significative_values = trend_type_to_series[AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND][i] \
+ # + trend_type_to_series[AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND][i]
+ # ax.plot(xlabel_values, significative_values, 'y-' + marker,
+ # label=AbstractUnivariateTest.SIGNIFICATIVE + ' trends')
# Global information
ax.set_ylabel(self.get_title_plot(xlabel, ax_idx=0))
@@ -144,7 +144,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
def visualize_altitude_trend_test(self, axes=None, marker='o', add_detailed_plots=False):
def altitude_reduction(df_bool, level):
# Take the mean with respect to the years
- df_bool = df_bool.mean(axis=1)
+ df_bool = df_bool.any(axis=1)
# Take the mean with respect the massifs
return df_bool.mean(level=level)
@@ -158,7 +158,7 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
def visualize_massif_trend_test(self, axes=None, add_detailed_plots=False):
def massif_reduction(df_bool, level):
# Take the mean with respect to the years
- df_bool = df_bool.mean(axis=1)
+ df_bool = df_bool.any(axis=1)
# Take the mean with respect the altitude
return df_bool.mean(level=level)
diff --git a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py
index 13f2c34c57c30f6c25c33502d9e408d81abd0937..2646e1ad01ebb246fad20a803fa71785de669bfd 100644
--- a/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py
+++ b/experiment/meteo_france_SCM_study/visualization/hypercube_visualization/main_hypercube_visualization.py
@@ -9,9 +9,8 @@ from experiment.meteo_france_SCM_study.visualization.hypercube_visualization.qua
from experiment.meteo_france_SCM_study.visualization.study_visualization.main_study_visualizer import ALL_ALTITUDES, \
SCM_STUDIES, study_iterator, study_iterator_global
from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
-from experiment.trend_analysis.univariate_trend_test.abstract_gev_trend_test import GevLocationTrendTest, \
- GevScaleTrendTest, GevShapeTrendTest
-from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import MannKendallTrendTest
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevLocationChangePointTest, \
+ GevScaleChangePointTest, GevShapeChangePointTest
from utils import get_display_name_from_object_type
@@ -33,14 +32,14 @@ from utils import get_display_name_from_object_type
# visualizer.visualize_altitude_trend_test()
-def full_trends_with_quantity_altitude_hypercube():
+def full_quantity_altitude_hypercube():
save_to_file = True
only_first_one = False
fast = False
add_detailed_plots = True
altitudes = ALL_ALTITUDES[3:-6]
study_classes = SCM_STUDIES
- for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][:]:
+ for trend_test_class in [GevLocationChangePointTest, GevScaleChangePointTest, GevShapeChangePointTest][:]:
visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False, multiprocessing=True)
for study in study_iterator_global(study_classes=study_classes, only_first_one=only_first_one,
altitudes=altitudes)]
@@ -54,13 +53,13 @@ def full_trends_with_quantity_altitude_hypercube():
visualizer.visualize_altitude_trend_test(add_detailed_plots=add_detailed_plots)
-def fast_trends_with_altitude_hypercube():
+def fast_altitude_hypercube():
save_to_file = False
only_first_one = False
fast = True
altitudes = ALL_ALTITUDES[2:4]
for study_class in SCM_STUDIES[:1]:
- for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][1:2]:
+ for trend_test_class in [GevLocationChangePointTest, GevScaleChangePointTest, GevShapeChangePointTest][1:2]:
visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False, multiprocessing=True)
for study in study_iterator(study_class=study_class, only_first_one=only_first_one,
altitudes=altitudes)]
@@ -72,13 +71,13 @@ def fast_trends_with_altitude_hypercube():
visualizer.visualize_altitude_trend_test()
-def fast_trends_with_quantity_altitude_hypercube():
+def fast_quantity_altitude_hypercube():
save_to_file = False
only_first_one = False
fast = True
altitudes = ALL_ALTITUDES[2:4]
study_classes = SCM_STUDIES[:2]
- for trend_test_class in [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][:1]:
+ for trend_test_class in [GevLocationChangePointTest, GevScaleChangePointTest, GevShapeChangePointTest][:1]:
visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False, multiprocessing=True)
for study in study_iterator_global(study_classes=study_classes, only_first_one=only_first_one,
altitudes=altitudes)]
@@ -93,9 +92,9 @@ def fast_trends_with_quantity_altitude_hypercube():
def main_run():
- # fast_trends_with_altitude_hypercube()
- # fast_trends_with_quantity_altitude_hypercube()
- full_trends_with_quantity_altitude_hypercube()
+ fast_altitude_hypercube()
+ # fast_quantity_altitude_hypercube()
+ # full_quantity_altitude_hypercube()
if __name__ == '__main__':
diff --git a/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py b/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py
index 23c2c0fa56c95bd4594266c58cae4124d380145e..dbfb52b61814050952726c5c0a3de22095dac92c 100644
--- a/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/studies_visualization/main_studies_visualizer.py
@@ -1,8 +1,8 @@
import time
from experiment.trend_analysis.abstract_score import MannKendall, WeigthedScore, MeanScore, MedianScore
-from experiment.trend_analysis.univariate_trend_test.abstract_gev_trend_test import GevLocationTrendTest, \
- GevScaleTrendTest, GevShapeTrendTest
-from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import MannKendallTrendTest
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevLocationChangePointTest, \
+ GevScaleChangePointTest, GevShapeChangePointTest
+from experiment.trend_analysis.univariate_test.abstract_univariate_test import MannKendallTrendTest
from experiment.meteo_france_SCM_study.safran.safran import ExtendedSafranTotalPrecip
from experiment.meteo_france_SCM_study.visualization.studies_visualization.studies import Studies
from experiment.meteo_france_SCM_study.visualization.studies_visualization.studies_visualizer import StudiesVisualizer, \
@@ -47,7 +47,7 @@ def altitude_trends_significant():
# altitudes = ALL_ALTITUDES[3:5]
altitudes = ALL_ALTITUDES[2:4]
for study_class in SCM_STUDIES[:1]:
- trend_test_classes = [MannKendallTrendTest, GevLocationTrendTest, GevScaleTrendTest, GevShapeTrendTest][:1]
+ trend_test_classes = [MannKendallTrendTest, GevLocationChangePointTest, GevScaleChangePointTest, GevShapeChangePointTest][:1]
visualizers = [StudyVisualizer(study, temporal_non_stationarity=True, verbose=False)
for study in study_iterator_global(study_classes=[study_class], only_first_one=only_first_one,
altitudes=altitudes)]
diff --git a/experiment/meteo_france_SCM_study/visualization/studies_visualization/studies_visualizer.py b/experiment/meteo_france_SCM_study/visualization/studies_visualization/studies_visualizer.py
index aa4b6c68742e352d99dfa33f5b9712694e0a7f1d..3b618abd6b0f53e9d0388759a618200f9c3fc7da 100644
--- a/experiment/meteo_france_SCM_study/visualization/studies_visualization/studies_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/studies_visualization/studies_visualizer.py
@@ -10,7 +10,7 @@ import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from experiment.meteo_france_SCM_study.abstract_extended_study import AbstractExtendedStudy
-from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import AbstractTrendTest
+from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
from experiment.meteo_france_SCM_study.visualization.studies_visualization.studies import \
Studies
from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
@@ -212,7 +212,7 @@ class AltitudeVisualizer(object):
# Add the marker legend
names = [get_display_name_from_object_type(c) for c in trend_test_classes]
- idx_for_positive_trend = [i for i, label in enumerate(labels) if label == AbstractTrendTest.POSITIVE_TREND]
+ idx_for_positive_trend = [i for i, label in enumerate(labels) if label == AbstractUnivariateTest.POSITIVE_TREND]
handles_ax2, labels_ax2 = [handles[i] for i in idx_for_positive_trend], names
ax2 = ax.twinx()
ax2.legend(handles_ax2, labels_ax2, loc=2)
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
index 645955e4657eede24fe0b877f4d099b6a70a153f..4f74b69f8c0bb2438ea59dcbb6c1d4b9c6c7b318 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/main_study_visualizer.py
@@ -11,6 +11,7 @@ from experiment.meteo_france_SCM_study.safran.safran import SafranSnowfall, Exte
from experiment.meteo_france_SCM_study.visualization.study_visualization.study_visualizer import StudyVisualizer
from collections import OrderedDict
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevLocationChangePointTest
from spatio_temporal_dataset.coordinates.transformed_coordinates.transformation.uniform_normalization import \
BetweenZeroAndOneNormalization, BetweenMinusOneAndOneNormalization
@@ -166,9 +167,10 @@ def trend_analysis():
verbose=True,
multiprocessing=True,
complete_non_stationary_trend_analysis=True)
- study_visualizer.visualize_all_independent_temporal_trend()
- study_visualizer.visualize_temporal_trend_relevance()
-
+ # study_visualizer.visualize_all_independent_temporal_trend()
+ # study_visualizer.visualize_temporal_trend_relevance()
+ study_visualizer.df_trend_spatio_temporal(GevLocationChangePointTest,
+ starting_years=[1958, 1980], nb_massif_for_fast_mode=2)
def maxima_analysis():
save_to_file = False
diff --git a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
index 722e28bad7cf46fa5076109d7d022bdd49648572..a5deda714f8a0c1168074bec984df70fb6143264 100644
--- a/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
+++ b/experiment/meteo_france_SCM_study/visualization/study_visualization/study_visualizer.py
@@ -11,7 +11,8 @@ import seaborn as sns
from experiment.trend_analysis.abstract_score import MeanScore, AbstractTrendScore
from experiment.meteo_france_SCM_study.abstract_study import AbstractStudy
-from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import AbstractTrendTest
+from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import AbstractGevChangePointTest
+from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
from experiment.trend_analysis.non_stationary_trends import \
ConditionalIndedendenceLocationTrendTest, MaxStableLocationTrendTest, IndependenceLocationTrendTest
from experiment.meteo_france_SCM_study.visualization.utils import create_adjusted_axes
@@ -376,8 +377,7 @@ class StudyVisualizer(object):
trend_type_and_weight = []
years, smooth_maxima = self.smooth_maxima_x_y(massif_id)
for starting_year, weight in starting_year_to_weight.items():
- test_trend_res = self.compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years)
- test_trend_type = test_trend_res[0]
+ test_trend_type = self.compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years)
trend_type_and_weight.append((test_trend_type, weight))
df = pd.DataFrame(trend_type_and_weight, columns=['trend type', 'weight'])
massif_name_to_df_trend_type[massif_name] = df
@@ -402,12 +402,17 @@ class StudyVisualizer(object):
list_args = [(smooth_maxima, starting_year, trend_test_class, years) for starting_year in
starting_years]
with Pool(NB_CORES) as p:
- trend_test_res = p.starmap(self.compute_trend_test_result, list_args)
+ trend_test_res = p.starmap(self.compute_gev_change_point_test_result, list_args)
else:
- trend_test_res = [self.compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years)
- for starting_year in starting_years]
+ trend_test_res = [
+ self.compute_gev_change_point_test_result(smooth_maxima, starting_year, trend_test_class, years)
+ 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)]
massif_name_to_trend_res[massif_name] = list(zip(*trend_test_res))
- # Build one DataFrame per trend test res
nb_res = len(list(massif_name_to_trend_res.values())[0])
assert nb_res == 2
all_massif_name_to_res = [{k: v[idx_res] for k, v in massif_name_to_trend_res.items()}
@@ -415,12 +420,16 @@ class StudyVisualizer(object):
return [pd.DataFrame(massif_name_to_res, index=starting_years).transpose()
for massif_name_to_res in all_massif_name_to_res]
+ @staticmethod
+ 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
+
@staticmethod
def compute_trend_test_result(smooth_maxima, starting_year, trend_test_class, years):
- idx = years.index(starting_year)
- # assert years[0] == starting_year, "{} {}".format(years[0], starting_year)
- trend_test = trend_test_class(years[:][idx:], smooth_maxima[:][idx:]) # type: AbstractTrendTest
- return trend_test.test_trend_type, trend_test.test_trend_strength
+ trend_test = trend_test_class(years, smooth_maxima, starting_year) # type: AbstractUnivariateTest
+ return trend_test.test_trend_type
def df_trend_test_count(self, trend_test_class, starting_year_to_weight):
"""
@@ -437,10 +446,10 @@ class StudyVisualizer(object):
df.fillna(0.0, inplace=True)
assert np.allclose(df.sum(axis=0), 100)
# Add the significant trend into the count of normal trend
- if AbstractTrendTest.SIGNIFICATIVE_POSITIVE_TREND in df.index:
- df.loc[AbstractTrendTest.POSITIVE_TREND] += df.loc[AbstractTrendTest.SIGNIFICATIVE_POSITIVE_TREND]
- if AbstractTrendTest.SIGNIFICATIVE_NEGATIVE_TREND in df.index:
- df.loc[AbstractTrendTest.NEGATIVE_TREND] += df.loc[AbstractTrendTest.SIGNIFICATIVE_NEGATIVE_TREND]
+ if AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND in df.index:
+ df.loc[AbstractUnivariateTest.POSITIVE_TREND] += df.loc[AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND]
+ if AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND in df.index:
+ df.loc[AbstractUnivariateTest.NEGATIVE_TREND] += df.loc[AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND]
return df
@cached_property
diff --git a/experiment/trend_analysis/univariate_trend_test/__init__.py b/experiment/trend_analysis/univariate_test/__init__.py
similarity index 100%
rename from experiment/trend_analysis/univariate_trend_test/__init__.py
rename to experiment/trend_analysis/univariate_test/__init__.py
diff --git a/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
similarity index 76%
rename from experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py
rename to experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
index eb09ead0f213fdfa5f8a5976a3620677149474e9..3d1b29001ff4d4c772259b3a3e2ff443fc28e360 100644
--- a/experiment/trend_analysis/univariate_trend_test/abstract_gev_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_gev_change_point_test.py
@@ -2,7 +2,7 @@ import numpy as np
import pandas as pd
from scipy.stats import chi2
-from experiment.trend_analysis.univariate_trend_test.abstract_trend_test import AbstractTrendTest
+from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import LinearMarginEstimator
from extreme_estimator.extreme_models.margin_model.param_function.linear_coef import LinearCoef
from extreme_estimator.extreme_models.margin_model.temporal_linear_margin_model import StationaryStationModel, \
@@ -19,14 +19,14 @@ from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_tempor
AbstractSpatioTemporalObservations
-class AbstractGevTrendTest(AbstractTrendTest):
+class AbstractGevChangePointTest(AbstractUnivariateTest):
RRunTimeError_TREND = 'R RunTimeError trend'
- def __init__(self, years_after_change_point, maxima_after_change_point, non_stationary_model_class, gev_param_name):
- super().__init__(years_after_change_point, maxima_after_change_point)
+ def __init__(self, years, maxima, starting_year, non_stationary_model_class, gev_param_name):
+ super().__init__(years, maxima, starting_year)
self.gev_param_name = gev_param_name
- df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: years_after_change_point})
- df_maxima_gev = pd.DataFrame(maxima_after_change_point, index=df.index)
+ df = pd.DataFrame({AbstractCoordinates.COORDINATE_T: years})
+ df_maxima_gev = pd.DataFrame(maxima, index=df.index)
observations = AbstractSpatioTemporalObservations(df_maxima_gev=df_maxima_gev)
self.coordinates = AbstractTemporalCoordinates.from_df(df, transformation_class=CenteredScaledNormalization) # type: AbstractTemporalCoordinates
self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
@@ -37,8 +37,8 @@ class AbstractGevTrendTest(AbstractTrendTest):
self.stationary_estimator.fit()
# Fit non stationary model
- self.non_stationary_estimator = LinearMarginEstimator(self.dataset,
- non_stationary_model_class(self.coordinates))
+ non_stationary_model = non_stationary_model_class(self.coordinates, starting_point=self.starting_year)
+ self.non_stationary_estimator = LinearMarginEstimator(self.dataset, non_stationary_model)
self.non_stationary_estimator.fit()
self.crashed = False
except SafeRunException:
@@ -49,6 +49,13 @@ class AbstractGevTrendTest(AbstractTrendTest):
return 2 * (self.non_stationary_estimator.result_from_fit.deviance -
self.stationary_estimator.result_from_fit.deviance)
+ @property
+ def non_stationary_nllh(self):
+ if self.crashed:
+ return -np.inf
+ else:
+ return self.non_stationary_estimator.result_from_fit.nllh
+
@property
def is_significant(self) -> bool:
return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=1)
@@ -91,29 +98,29 @@ class AbstractGevTrendTest(AbstractTrendTest):
ratio = np.abs(self.non_stationary_linear_coef) / np.abs(self.non_stationary_intercept_coef)
scaled_ratio = ratio * self.coordinates.transformed_distance_between_two_successive_years
percentage_of_change_per_year = 100 * scaled_ratio
- return percentage_of_change_per_year
+ return percentage_of_change_per_year[0]
@property
def test_sign(self) -> int:
return np.sign(self.non_stationary_linear_coef)
-class GevLocationTrendTest(AbstractGevTrendTest):
+class GevLocationChangePointTest(AbstractGevChangePointTest):
- def __init__(self, years_after_change_point, maxima_after_change_point):
- super().__init__(years_after_change_point, maxima_after_change_point,
+ def __init__(self, years, maxima, starting_year):
+ super().__init__(years, maxima, starting_year,
NonStationaryLocationStationModel, GevParams.LOC)
-class GevScaleTrendTest(AbstractGevTrendTest):
+class GevScaleChangePointTest(AbstractGevChangePointTest):
- def __init__(self, years_after_change_point, maxima_after_change_point):
- super().__init__(years_after_change_point, maxima_after_change_point,
+ def __init__(self, years, maxima, starting_year):
+ super().__init__(years, maxima, starting_year,
NonStationaryScaleStationModel, GevParams.SCALE)
-class GevShapeTrendTest(AbstractGevTrendTest):
+class GevShapeChangePointTest(AbstractGevChangePointTest):
- def __init__(self, years_after_change_point, maxima_after_change_point):
- super().__init__(years_after_change_point, maxima_after_change_point,
+ def __init__(self, years, maxima, starting_year):
+ super().__init__(years, maxima, starting_year,
NonStationaryShapeStationModel, GevParams.SHAPE)
diff --git a/experiment/trend_analysis/univariate_test/abstract_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_trend_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..195d685f8b6c43b236cda26200e70bb079d1784a
--- /dev/null
+++ b/experiment/trend_analysis/univariate_test/abstract_trend_test.py
@@ -0,0 +1,45 @@
+import random
+import warnings
+
+import matplotlib.pyplot as plt
+from collections import OrderedDict
+
+import numpy as np
+from cached_property import cached_property
+
+from experiment.trend_analysis.mann_kendall_test import mann_kendall_test
+from experiment.trend_analysis.abstract_score import MannKendall
+from experiment.trend_analysis.univariate_test.abstract_univariate_test import AbstractUnivariateTest
+
+
+class MannKendallTrendTest(AbstractUnivariateTest):
+
+ def __init__(self, years, maxima, starting_year):
+ super().__init__(years, maxima, starting_year)
+ score = MannKendall()
+ # Compute score value
+ detailed_score = score.get_detailed_score(self.years_after_starting_year, self.maxima_after_starting_year)
+ self.score_value = detailed_score[0]
+ # Compute the Mann Kendall Test
+ MK, S = mann_kendall_test(t=self.years_after_starting_year,
+ x=self.maxima_after_starting_year,
+ eps=1e-5,
+ alpha=self.SIGNIFICANCE_LEVEL,
+ Ha='upordown')
+ # Raise warning if scores are differents
+ if S != self.score_value:
+ warnings.warn('S={} is different that score_value={}'.format(S, self.score_value), WarningScoreValue)
+ self.MK = MK
+
+ @property
+ def test_sign(self) -> int:
+ return np.sign(self.score_value)
+
+ @property
+ def is_significant(self) -> bool:
+ assert 'reject' in self.MK or 'accept' in self.MK
+ return 'accept' in self.MK
+
+
+class SpearmanRhoTrendTest(AbstractUnivariateTest):
+ pass
diff --git a/experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
similarity index 60%
rename from experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py
rename to experiment/trend_analysis/univariate_test/abstract_univariate_test.py
index f08e7daa89119d6bad060d30c961040a0157ac2f..2da67e6338c107bcb11aea0319d15985dcbef0ca 100644
--- a/experiment/trend_analysis/univariate_trend_test/abstract_trend_test.py
+++ b/experiment/trend_analysis/univariate_test/abstract_univariate_test.py
@@ -5,12 +5,13 @@ import matplotlib.pyplot as plt
from collections import OrderedDict
import numpy as np
+from cached_property import cached_property
from experiment.trend_analysis.mann_kendall_test import mann_kendall_test
from experiment.trend_analysis.abstract_score import MannKendall
-class AbstractTrendTest(object):
+class AbstractUnivariateTest(object):
SIGNIFICATIVE = 'significative'
# 5 possible types of trends
NO_TREND = 'no trend'
@@ -21,10 +22,23 @@ class AbstractTrendTest(object):
SIGNIFICANCE_LEVEL = 0.05
- def __init__(self, years_after_change_point, maxima_after_change_point):
- self.years_after_change_point = years_after_change_point
- self.maxima_after_change_point = maxima_after_change_point
- assert len(self.years_after_change_point) == len(self.maxima_after_change_point)
+ def __init__(self, years, maxima, starting_year):
+ self.years = years
+ self.maxima = maxima
+ self.starting_year = starting_year
+ assert len(self.years) == len(self.maxima)
+
+ @cached_property
+ def idx_for_starting_year(self):
+ return self.years.index(self.starting_year)
+
+ @property
+ def years_after_starting_year(self):
+ return self.years[self.idx_for_starting_year:]
+
+ @property
+ def maxima_after_starting_year(self):
+ return self.maxima[self.idx_for_starting_year:]
@classmethod
def trend_type_to_style(cls):
@@ -54,11 +68,9 @@ class AbstractTrendTest(object):
else:
return plt.cm.binary
-
-
@property
def n(self):
- return len(self.years_after_change_point)
+ return len(self.years)
@property
def test_trend_strength(self):
@@ -86,7 +98,7 @@ class AbstractTrendTest(object):
raise NotImplementedError
-class ExampleRandomTrendTest(AbstractTrendTest):
+class ExampleRandomTrendTest(AbstractUnivariateTest):
@property
def test_sign(self) -> int:
@@ -99,36 +111,3 @@ class ExampleRandomTrendTest(AbstractTrendTest):
class WarningScoreValue(Warning):
pass
-
-
-class MannKendallTrendTest(AbstractTrendTest):
-
- def __init__(self, years_after_change_point, maxima_after_change_point):
- super().__init__(years_after_change_point, maxima_after_change_point)
- score = MannKendall()
- # Compute score value
- detailed_score = score.get_detailed_score(years_after_change_point, maxima_after_change_point)
- self.score_value = detailed_score[0]
- # Compute the Mann Kendall Test
- MK, S = mann_kendall_test(t=years_after_change_point,
- x=maxima_after_change_point,
- eps=1e-5,
- alpha=self.SIGNIFICANCE_LEVEL,
- Ha='upordown')
- # Raise warning if scores are differents
- if S != self.score_value:
- warnings.warn('S={} is different that score_value={}'.format(S, self.score_value), WarningScoreValue)
- self.MK = MK
-
- @property
- def test_sign(self) -> int:
- return np.sign(self.score_value)
-
- @property
- def is_significant(self) -> bool:
- assert 'reject' in self.MK or 'accept' in self.MK
- return 'accept' in self.MK
-
-
-class SpearmanRhoTrendTest(AbstractTrendTest):
- pass