Commit 94aa79b3 authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[PAPER1] modify the strength criterion for the hypercube (use the 0.98...

[PAPER1] modify the strength criterion for the hypercube (use the 0.98 quantile and its evolution in 10 years)
parent e6c54f99
No related merge requests found
Showing with 219 additions and 32 deletions
+219 -32
......@@ -299,7 +299,7 @@ class AbstractStudy(object):
x, y = list(row)
massif_name = row.name
value = massif_name_to_value[massif_name]
str_value = str(round(value, 1)) if isinstance(value, str) else str(value)
str_value = str(value)
ax.text(x, y, str_value, horizontalalignment='center', verticalalignment='center', fontsize=7)
if scaled:
......
......@@ -20,11 +20,13 @@ class AbstractHypercubeVisualizer(object):
def __init__(self, tuple_to_study_visualizer: Dict[Tuple, StudyVisualizer],
trend_test_class,
nb_data_reduced_for_speed=False,
reduce_strength_array=False,
save_to_file=False,
first_starting_year=None,
last_starting_year=None,
exact_starting_year=None,
verbose=True):
self.reduce_strength_array = reduce_strength_array
self.verbose = verbose
self.save_to_file = save_to_file
self.trend_test_class = trend_test_class
......@@ -87,13 +89,17 @@ class AbstractHypercubeVisualizer(object):
)
@cached_property
def df_hypercube_trend_strength(self) -> pd.DataFrame:
def df_hypercube_trend_slope_relative_strength(self) -> pd.DataFrame:
return self._df_hypercube_trend_meta(idx=1)
@cached_property
def df_hypercube_trend_nllh(self) -> pd.DataFrame:
return self._df_hypercube_trend_meta(idx=2)
@cached_property
def df_hypercube_trend_constant_quantile(self) -> pd.DataFrame:
return self._df_hypercube_trend_meta(idx=3)
# Some properties
@property
......
......@@ -9,7 +9,9 @@ from experiment.meteo_france_data.scm_models_data.visualization.study_visualizat
SCM_STUDY_NAME_TO_COLOR, SCM_STUDY_NAME_TO_ABBREVIATION, SCM_STUDY_CLASS_TO_ABBREVIATION, SCM_STUDIES_NAMES
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
StudyVisualizer
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 extreme_estimator.margin_fits.gev.gev_params import GevParams
ALTITUDES_XLABEL = 'altitudes'
......@@ -71,11 +73,13 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
assert not isinstance(s_trend_type_percentage.index, pd.MultiIndex)
s_trend_type_percentage *= 100
series = [s_trend_type_percentage]
# # Reduce df_strength to a serie s_trend_strength
# df_strength = self.df_hypercube_trend_strength[df_bool]
# s_trend_strength = reduction_function(df_strength)
# # Group result
# series = [s_trend_type_percentage, s_trend_strength]
if self.reduce_strength_array:
# Reduce df_strength to a serie s_trend_strength
df_strength = self.df_hypercube_trend_slope_relative_strength[df_bool]
s_trend_strength = reduction_function(df_strength)
df_constant = self.df_hypercube_trend_constant_quantile[df_bool]
s_trend_constant = reduction_function(df_constant)
series.extend([s_trend_strength, s_trend_constant])
return series
def subtitle_to_reduction_function(self, reduction_function, level=None, add_detailed_plot=False, subtitle=None):
......@@ -280,6 +284,8 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
massif_to_color = {}
add_text = self.nb_rows > 1
massif_to_year = {}
massif_to_strength = {}
massif_to_constant = {}
poster_trend_types = [AbstractUnivariateTest.SIGNIFICATIVE_POSITIVE_TREND,
AbstractUnivariateTest.SIGNIFICATIVE_NEGATIVE_TREND,
AbstractUnivariateTest.NEGATIVE_TREND,
......@@ -292,15 +298,36 @@ class AltitudeHypercubeVisualizer(AbstractHypercubeVisualizer):
massif_to_color_for_trend_type = {k: color for k, v in dict(serie).items() if not np.isnan(v)}
massif_to_color.update(massif_to_color_for_trend_type)
if add_text:
massif_to_year_for_trend_type = {k: int(v) for k, v in
self.trend_type_to_series(reduction_function, isin_parameters)[
display_trend_type][1].items()
if k in massif_to_color_for_trend_type}
massif_to_year.update(massif_to_year_for_trend_type)
if self.reduce_strength_array:
massif_to_value_for_trend_type = [{k: v for k, v in
self.trend_type_to_series(reduction_function,
isin_parameters)[
display_trend_type][i].items()
if k in massif_to_color_for_trend_type} for i in [1, 2]]
massif_to_strength.update(massif_to_value_for_trend_type[0])
massif_to_constant.update(massif_to_value_for_trend_type[1])
else:
massif_to_value_for_trend_type = {k: int(v) for k, v in
self.trend_type_to_series(reduction_function,
isin_parameters)[
display_trend_type][1].items()
if k in massif_to_color_for_trend_type}
massif_to_year.update(massif_to_value_for_trend_type)
# Compute massif_to_value
if self.reduce_strength_array:
massif_name_to_value = {m: "{} {}{} / {} year(s)".format(
int(massif_to_constant[m]),
"+" if massif_to_strength[m] > 0 else "",
round(massif_to_strength[m] * massif_to_constant[m], 1),
AbstractGevChangePointTest.nb_years_for_quantile_evolution)
for m in massif_to_strength}
else:
massif_name_to_value = massif_to_year
self.study.visualize_study(None, massif_name_to_color=massif_to_color, show=False,
show_label=False, scaled=True, add_text=add_text,
massif_name_to_value=massif_to_year)
massif_name_to_value=massif_name_to_value)
print(subtitle)
title = self.set_trend_test_reparition_title(subtitle, set=False)
# row_title = self.get_title_plot(xlabel='massifs', ax_idx=i)
......
......@@ -463,11 +463,11 @@ class StudyVisualizer(VisualizationParameters):
sample_one_massif_from_each_region)
for massif_name, gev_change_point_test_results in massif_name_to_gev_change_point_test_results.items():
trend_test_res, best_idxs = gev_change_point_test_results
trend_test_res = [(a, b, c) if i in best_idxs else (np.nan, np.nan, c)
for i, (a, b, c, *_) in enumerate(trend_test_res)]
trend_test_res = [(a, b, c, d) if i in best_idxs else (np.nan, np.nan, c, np.nan)
for i, (a, b, c, d, *_) 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 == 3
assert nb_res == 4
all_massif_name_to_res = [{k: v[idx_res] for k, v in massif_name_to_trend_res.items()}
for idx_res in range(nb_res)]
return [pd.DataFrame(massif_name_to_res, index=self.starting_years_for_change_point_test).transpose()
......
......@@ -251,10 +251,10 @@ class ComparisonsVisualization(VisualizationParameters):
label += "\n {} starting in {}".format(display_trend_type, most_likely_year)
ordered_dict[TREND_TYPE_CNAME] = display_trend_type
ordered_dict['most likely year'] = most_likely_year
# Display the nllh against the starting year
# Display the deviance against the starting year
step = 1
ax2.plot(starting_years[::step], [t[3] for t in trend_test_res][::step], color=plot_color, marker='o')
ax2.plot(starting_years[::step], [t[4] for t in trend_test_res][::step], color=plot_color, marker='x')
ax2.plot(starting_years[::step], [t[4] for t in trend_test_res][::step], color=plot_color, marker='o')
ax2.plot(starting_years[::step], [t[5] for t in trend_test_res][::step], color=plot_color, marker='x')
# Plot maxima
ax.grid()
ax.plot(years, maxima, label=label, color=plot_color)
......@@ -270,7 +270,7 @@ class ComparisonsVisualization(VisualizationParameters):
res = safe_run_r_estimator(function=r('gev.fit'), xdat=ro.FloatVector(data),
use_start=True)
res = ResultFromIsmev(res, {})
gev_params = res.stationary_gev_params
gev_params = res.constant_gev_params
lim = 1.5 * max(data)
x = np.linspace(0, lim, 1000)
......
import time
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
Altitude_Hypercube_Year_Visualizer
from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
def main_fast_spatial_repartition():
for altitude in FULL_ALTITUDES[-1:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958)
vizualiser.save_to_file = False
vizualiser.visualize_massif_trend_test_one_altitude()
def main_full_spatial_repartition():
for altitude in FULL_ALTITUDES[:]:
# Compute for the most likely starting year
# vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
# vizualiser.visualize_massif_trend_test_one_altitude()
# Compute the trend for a linear trend
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958)
vizualiser.visualize_massif_trend_test_one_altitude()
def main_run():
# main_full_spatial_repartition()
main_fast_spatial_repartition()
if __name__ == '__main__':
start = time.time()
main_run()
duration = time.time() - start
print('Full run took {}s'.format(round(duration, 1)))
import time
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
Altitude_Hypercube_Year_Visualizer
"""
Visualize the 0.99 quantile initial value and its evolution
"""
from experiment.paper1_steps.utils import get_full_altitude_visualizer, FULL_ALTITUDES
def main_fast_spatial_risk_evolution():
for altitude in FULL_ALTITUDES[-1:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958, reduce_strength_array=True)
vizualiser.save_to_file = False
vizualiser.visualize_massif_trend_test_one_altitude()
def main_full_spatial_repartition():
for altitude in FULL_ALTITUDES[:]:
# Compute for the most likely starting year
# vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
# vizualiser.visualize_massif_trend_test_one_altitude()
# Compute the trend for a linear trend
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958)
vizualiser.visualize_massif_trend_test_one_altitude()
def main_run():
# main_full_spatial_repartition()
main_fast_spatial_risk_evolution()
if __name__ == '__main__':
start = time.time()
main_run()
duration = time.time() - start
print('Full run took {}s'.format(round(duration, 1)))
import time
from experiment.meteo_france_data.scm_models_data.crocus.crocus import CrocusRecentSwe
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer import \
AltitudeHypercubeVisualizer
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_hypercube_visualizer_extended import \
AltitudeHypercubeVisualizerBisExtended, QuantityHypercubeWithoutTrendExtended, \
AltitudeHypercubeVisualizerWithoutTrendExtended, QuantityHypercubeWithoutTrend
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.altitude_year_hypercube_visualizer import \
Altitude_Hypercube_Year_Visualizer
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.main_files.main_fast_hypercube_one_altitudes import \
get_fast_parameters
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.main_files.main_full_hypercube import \
get_full_parameters
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.quantity_altitude_visualizer import \
QuantityAltitudeHypercubeVisualizer
from experiment.meteo_france_data.scm_models_data.visualization.hypercube_visualization.utils_hypercube import \
load_altitude_visualizer, load_quantity_visualizer
from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.main_study_visualizer import \
ALL_ALTITUDES, SCM_STUDIES
from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevLocationChangePointTest, \
GevScaleChangePointTest
load_altitude_visualizer
from experiment.trend_analysis.univariate_test.abstract_gev_change_point_test import GevScaleChangePointTest, \
GevLocationChangePointTest
def get_fast_altitude_visualizer(altitude_hypercube_class):
altitudes, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, trend_test_class = get_fast_parameters()
study_classes = [CrocusRecentSwe]
visualizer = load_altitude_visualizer(altitude_hypercube_class, altitudes, last_starting_year,
nb_data_reduced_for_speed, only_first_one, save_to_file, study_classes,
trend_test_class)
return visualizer
def main_fast_old_spatial_repartition():
# Simply the main graph
get_fast_altitude_visualizer(Altitude_Hypercube_Year_Visualizer).visualize_massif_trend_test_one_altitude()
FULL_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
def get_full_altitude_visualizer(altitude_hypercube_class, exact_starting_year=None, altitude=900):
def get_full_altitude_visualizer(altitude_hypercube_class, exact_starting_year=None, altitude=900,
reduce_strength_array=False):
altitudes, first_starting_year, last_starting_year, nb_data_reduced_for_speed, only_first_one, save_to_file, trend_test_class = get_full_parameters(
altitude=altitude)
if exact_starting_year is not None:
first_starting_year, last_starting_year = None, None
study_classes = [CrocusRecentSwe]
trend_test_class = GevScaleChangePointTest
trend_test_class = GevLocationChangePointTest
visualizer = load_altitude_visualizer(altitude_hypercube_class, altitudes, last_starting_year,
nb_data_reduced_for_speed, only_first_one, save_to_file, study_classes,
trend_test_class, first_starting_year=first_starting_year,
exact_starting_year=exact_starting_year)
visualizer.reduce_strength_array = reduce_strength_array
return visualizer
FULL_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]
HALF_ALTITUDES = [900, 1200, 1500, 1800, 2100, 2400, 2700, 3000][::2]
def main_fast_spatial_repartition():
for altitude in FULL_ALTITUDES[-1:]:
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958)
vizualiser.save_to_file = False
vizualiser.visualize_massif_trend_test_one_altitude()
def main_full_spatial_repartition():
for altitude in FULL_ALTITUDES[:]:
# Compute for the most likely starting year
# vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude)
# vizualiser.visualize_massif_trend_test_one_altitude()
# Compute the trend for a linear trend
vizualiser = get_full_altitude_visualizer(Altitude_Hypercube_Year_Visualizer, altitude=altitude,
exact_starting_year=1958)
vizualiser.visualize_massif_trend_test_one_altitude()
def main_run():
main_full_spatial_repartition()
# main_fast_spatial_repartition()
if __name__ == '__main__':
start = time.time()
main_run()
duration = time.time() - start
print('Full run took {}s'.format(round(duration, 1)))
......@@ -21,6 +21,9 @@ from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_tempor
class AbstractGevChangePointTest(AbstractUnivariateTest):
RRunTimeError_TREND = 'R RunTimeError trend'
# I should use the quantile from the Eurocode for the buildings
quantile_for_strength = 0.98
nb_years_for_quantile_evolution = 10
def __init__(self, years, maxima, starting_year, non_stationary_model_class, gev_param_name):
super().__init__(years, maxima, starting_year)
......@@ -28,8 +31,12 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
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.coordinates = AbstractTemporalCoordinates.from_df(df,
transformation_class=CenteredScaledNormalization) # type: AbstractTemporalCoordinates
self.dataset = AbstractDataset(observations=observations, coordinates=self.coordinates)
# For the moment, we chose:
# -the 0.99 quantile (even if for building maybe we should use the 1/50 so 0.98 quantile)
# -to see the evolution between two successive years
try:
# Fit stationary model
......@@ -60,11 +67,8 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
return self.non_stationary_deviance - self.stationary_deviance
@property
def non_stationary_nllh(self):
if self.crashed:
return np.nan
else:
return self.non_stationary_estimator.result_from_fit.nllh
def non_stationary_constant_gev_params(self) -> GevParams:
return self.non_stationary_estimator.result_from_fit.constant_gev_params
@property
def stationary_deviance(self):
......@@ -80,6 +84,13 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
else:
return self.non_stationary_estimator.result_from_fit.deviance
@property
def non_stationary_nllh(self):
if self.crashed:
return np.nan
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)
......@@ -111,11 +122,23 @@ class AbstractGevChangePointTest(AbstractUnivariateTest):
return self.get_coef(self.non_stationary_estimator, AbstractCoordinates.COORDINATE_T)
@property
def test_trend_strength(self):
def test_trend_slope_strength(self):
if self.crashed:
return 0.0
else:
return self.percentage_of_change_per_year
slope = self._slope_strength()
slope *= self.nb_years_for_quantile_evolution * self.coordinates.transformed_distance_between_two_successive_years[0]
return slope
def _slope_strength(self):
raise NotImplementedError
@property
def test_trend_constant_quantile(self):
if self.crashed:
return 0.0
else:
return self.non_stationary_constant_gev_params.quantile(p=self.quantile_for_strength)
@property
def percentage_of_change_per_year(self):
......@@ -135,6 +158,10 @@ class GevLocationChangePointTest(AbstractGevChangePointTest):
super().__init__(years, maxima, starting_year,
NonStationaryLocationStationModel, GevParams.LOC)
def _slope_strength(self):
return self.non_stationary_constant_gev_params.quantile_strength_evolution_ratio(p=self.quantile_for_strength,
mu1=self.non_stationary_linear_coef)
class GevScaleChangePointTest(AbstractGevChangePointTest):
......
......@@ -105,7 +105,7 @@ class AbstractUnivariateTest(object):
return len(self.years)
@property
def test_trend_strength(self):
def test_trend_slope_strength(self):
return 0.0
@property
......
......@@ -9,7 +9,7 @@ from utils import NB_CORES
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, trend_test.non_stationary_deviance, trend_test.stationary_deviance
return trend_test.test_trend_type, trend_test.test_trend_slope_strength, trend_test.non_stationary_nllh, trend_test.test_trend_constant_quantile, trend_test.non_stationary_deviance, trend_test.stationary_deviance
def compute_gev_change_point_test_results(multiprocessing, maxima, starting_years, trend_test_class,
......
......@@ -18,7 +18,7 @@ class TemporalLinearMarginModel(LinearMarginModel):
def fitmargin_from_maxima_gev(self, data: np.ndarray, df_coordinates_spat: pd.DataFrame,
df_coordinates_temp: pd.DataFrame) -> ResultFromFit:
# Gev Fit
# Gev Fit with isMev package
assert data.shape[1] == len(df_coordinates_temp.values)
res = safe_run_r_estimator(function=r('gev.fit'), use_start=self.use_start_value,
xdat=ro.FloatVector(data[0]), y=df_coordinates_temp.values, mul=self.mul,
......
......@@ -48,6 +48,10 @@ class ResultFromFit(object):
def convergence(self) -> str:
raise NotImplementedError
@property
def constant_gev_params(self) -> GevParams:
raise NotImplementedError
class ResultFromIsmev(ResultFromFit):
......@@ -76,7 +80,7 @@ class ResultFromIsmev(ResultFromFit):
return coef_dict
@property
def stationary_gev_params(self) -> GevParams:
def constant_gev_params(self) -> GevParams:
params = {k.split('Coeff1')[0]: v for k, v in self.margin_coef_dict.items()
if 'Coeff1' in k and 'temp' not in k}
return GevParams.from_dict(params)
......
......@@ -45,6 +45,24 @@ class GevParams(ExtremeParams):
def __str__(self):
return self.to_dict().__str__()
def quantile_strength_evolution_ratio(self, p=0.99, mu1=0.0, sigma1=0.0):
"""
Compute the relative evolution of some quantile with respect to time.
(when mu1 and sigma1 can be modified with time)
:param p: level of the quantile
:param mu1: temporal slope of the location parameter
:param sigma1: temporal slope of the scale parameter
:return: A string summarizing the evolution ratio
"""
initial_quantile = self.quantile(p)
quantity_increased = mu1
if sigma1 != 0:
quantity_increased += (sigma1 / self.shape) * (1 - (- np.float_power(np.log(p), -self.shape)))
return quantity_increased / initial_quantile
# Compute some indicators (such as the mean and the variance)
def g(self, k):
# Compute the g_k parameters as defined in wikipedia
# https://fr.wikipedia.org/wiki/Loi_d%27extremum_g%C3%A9n%C3%A9ralis%C3%A9e
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment