one_fold_fit.py

import datetime
import math
import time
from itertools import chain
from multiprocessing import Pool

import numpy as np
from cached_property import cached_property
from scipy.stats import chi2
from sklearn.utils import resample

from extreme_fit.distribution.gev.gev_params import GevParams
from extreme_fit.distribution.gumbel.gumbel_gof import goodness_of_fit_anderson
from extreme_fit.estimator.margin_estimator.utils import fitted_linear_margin_estimator_short
from extreme_fit.function.param_function.polynomial_coef import PolynomialAllCoef, PolynomialCoef
from extreme_fit.model.margin_model.linear_margin_model.temporal_linear_margin_models import StationaryTemporalModel
from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models import StationaryAltitudinal
from extreme_fit.model.margin_model.polynomial_margin_model.gev_altitudinal_models_only_altitude_and_scale import \
    AltitudinalOnlyScale, StationaryAltitudinalOnlyScale
from extreme_fit.model.margin_model.polynomial_margin_model.gumbel_altitudinal_models import \
    StationaryGumbelAltitudinal, AbstractGumbelAltitudinalModel
from extreme_fit.model.margin_model.polynomial_margin_model.models_based_on_pariwise_analysis.gev_with_linear_shape_wrt_altitude import \
    AltitudinalShapeLinearTimeStationary
from extreme_fit.model.margin_model.utils import MarginFitMethod
from extreme_fit.model.result_from_model_fit.result_from_extremes.abstract_extract_eurocode_return_level import \
    AbstractExtractEurocodeReturnLevel
from extreme_fit.model.result_from_model_fit.result_from_extremes.confidence_interval_method import \
    ConfidenceIntervalMethodFromExtremes
from extreme_fit.model.result_from_model_fit.result_from_extremes.eurocode_return_level_uncertainties import \
    EurocodeConfidenceIntervalFromExtremes
from extreme_trend.one_fold_fit.altitude_group import DefaultAltitudeGroup, altitudes_for_groups
from root_utils import NB_CORES, batch
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
from spatio_temporal_dataset.coordinates.temporal_coordinates.abstract_temporal_covariate_for_fit import \
    TimeTemporalCovariate
from spatio_temporal_dataset.coordinates.temporal_coordinates.temperature_covariate import \
    AnomalyTemperatureWithSplineTemporalCovariate
from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima


class OneFoldFit(object):
    SIGNIFICANCE_LEVEL = 0.05
    return_period = 100
    quantile_level = 1 - (1 / return_period)

    def __init__(self, massif_name: str, dataset: AbstractDataset, models_classes,
                 first_year, last_year,
                 fit_method=MarginFitMethod.extremes_fevd_mle,
                 temporal_covariate_for_fit=None,
                 altitude_group=None,
                 only_models_that_pass_goodness_of_fit_test=True,
                 confidence_interval_based_on_delta_method=False,
                 remove_physically_implausible_models=False,
                 climate_coordinates_with_effects=None
                 ):
        self.first_year = first_year
        self.last_year = last_year
        self.years_of_difference = last_year - first_year
        self.remove_physically_implausible_models = remove_physically_implausible_models
        self.confidence_interval_based_on_delta_method = confidence_interval_based_on_delta_method
        self.only_models_that_pass_goodness_of_fit_test = only_models_that_pass_goodness_of_fit_test
        self.altitude_group = altitude_group
        self.massif_name = massif_name
        self.dataset = dataset
        self.models_classes = models_classes
        self.fit_method = fit_method
        self.temporal_covariate_for_fit = temporal_covariate_for_fit
        self.climate_coordinates_with_effects = climate_coordinates_with_effects

        # Fit Estimators
        self.model_class_to_estimator = {}
        for model_class in models_classes:
            self.model_class_to_estimator[model_class] = self.fitted_linear_margin_estimator(model_class, self.dataset)
        # Compute sorted estimators indirectly
        _ = self.has_at_least_one_valid_model

    def fitted_linear_margin_estimator(self, model_class, dataset):
        return fitted_linear_margin_estimator_short(model_class=model_class,
                                                    dataset=dataset,
                                                    fit_method=self.fit_method,
                                                    temporal_covariate_for_fit=self.temporal_covariate_for_fit,
                                                    drop_duplicates=False,
                                                    climate_coordinates_with_effects=self.climate_coordinates_with_effects)

    @classmethod
    def get_moment_str(cls, order):
        if order == 1:
            return 'mean'
        elif order == 2:
            return 'std'
        elif order is None:
            return '{}-year return levels'.format(cls.return_period)

    def get_moment_for_plots(self, altitudes, order=1, covariate_before=None, covariate_after=None):
        return [self.get_moment(altitudes[0], covariate_after, order)]

    def get_moment(self, altitude, temporal_covariate, order=1):
        gev_params = self.get_gev_params(altitude, temporal_covariate)
        if order == 1:
            return gev_params.mean
        elif order == 2:
            return gev_params.std
        elif order is None:
            return gev_params.return_level(return_period=self.return_period)
        elif order in GevParams.PARAM_NAMES:
            return gev_params.to_dict()[order]
        else:
            raise NotImplementedError

    def get_gev_params(self, altitude, year):
        coordinate = self.get_coordinate(altitude, year)
        gev_params = self.best_margin_function_from_fit.get_params(coordinate, is_transformed=False)
        return gev_params

    def moment(self, altitudes, order=1):
        return [self.get_moment(altitude, self.covariate_after, order) for altitude in altitudes]

    @property
    def change_in_return_level_for_reference_altitude(self) -> float:
        return self.changes_of_moment(altitudes=[self.altitude_plot], order=None)[0]

    @property
    def relative_change_in_return_level_for_reference_altitude(self) -> float:
        return self.relative_changes_of_moment(altitudes=[self.altitude_plot], order=None)[0]

    def changes_of_moment(self, altitudes, order=1, covariate_before=None, covariate_after=None):
        covariate_after, covariate_before = self.set_covariate_before_and_after(covariate_after, covariate_before)
        changes = []
        for altitude in altitudes:
            mean_after = self.get_moment(altitude, covariate_after, order)
            mean_before = self.get_moment(altitude, covariate_before, order)
            change = mean_after - mean_before
            changes.append(change)
        return changes

    def set_covariate_before_and_after(self, covariate_after, covariate_before):
        if covariate_before is None:
            covariate_before = self.covariate_before
        if covariate_after is None:
            covariate_after = self.covariate_after
        return covariate_after, covariate_before

    @property
    def covariate_before(self):
        return self._covariate_before_and_after[0]

    @property
    def covariate_after(self):
        return self._covariate_before_and_after[1]

    @property
    def _covariate_before_and_after(self):
        if self.temporal_covariate_for_fit in [None, TimeTemporalCovariate]:
            return self.first_year, self.last_year
        elif self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate:
            # In 2020, we are roughly at 1 degree. Thus it natural to see the augmentation from 1 to 2 degree.
            return 1, 2
        else:
            raise NotImplementedError

    @property
    def between_covariate_str(self):
        d_temperature = {'C': '{C}'}
        s = ' between +${}^o\mathrm{C}$ and +${}^o\mathrm{C}$' if self.temporal_covariate_for_fit is AnomalyTemperatureWithSplineTemporalCovariate \
            else ' between {} and {}'
        s = s.format(self.covariate_before, self.covariate_after,
                     **d_temperature)
        return s

    def relative_changes_of_moment(self, altitudes, order=1, covariate_before=None, covariate_after=None):
        covariate_after, covariate_before = self.set_covariate_before_and_after(covariate_after, covariate_before)
        relative_changes = []
        for altitude in altitudes:
            mean_after = self.get_moment(altitude, covariate_after, order)
            mean_before = self.get_moment(altitude, covariate_before, order)
            relative_change = 100 * (mean_after - mean_before) / np.abs(mean_before)
            relative_changes.append(relative_change)
        return relative_changes

    # Minimizing the AIC and some properties

    @cached_property
    def sorted_estimators_with_aic(self):
        return self._sorted_estimators_with_method_name(method_name='aic')

    def method_name_to_best_estimator(self, method_names):
        return {self._sorted_estimators_with_method_name(method_name) for method_name in method_names}

    def _sorted_estimators_with_method_name(self, method_name):
        estimators = self.estimators_quality_checked
        try:
            sorted_estimators = sorted([estimator for estimator in estimators],
                                       key=lambda e: e.__getattribute__(method_name))
        except AssertionError as e:
            print('Error for:\n', self.massif_name, self.altitude_group)
            raise e
        return sorted_estimators

    @cached_property
    def estimators_quality_checked(self):
        estimators = list(self.model_class_to_estimator.values())
        if self.remove_physically_implausible_models:
            # Remove wrong shape
            estimators = [e for e in estimators if -0.5 < self._compute_shape_for_reference_altitude(e) < 0.5]
            # Remove models with undefined parameters for the coordinate of interest
            well_defined_estimators = []
            for e in estimators:
                coordinate_values_for_the_fit = e.coordinates_for_nllh

                if isinstance(self.altitude_group, DefaultAltitudeGroup):
                    coordinate_values_for_the_result = []
                else:
                    coordinate_values_for_the_result = [np.array([self.altitude_group.reference_altitude, c])
                                                        for c in self._covariate_before_and_after]
                coordinate_values_to_check = list(coordinate_values_for_the_fit) + coordinate_values_for_the_result
                has_undefined_parameters = False
                for coordinate in coordinate_values_to_check:
                    gev_params = e.margin_function_from_fit.get_params(coordinate)
                    if gev_params.has_undefined_parameters:
                        has_undefined_parameters = True
                        break
                if not has_undefined_parameters:
                    well_defined_estimators.append(e)
                else:
                    print('undefined', e)
            estimators = well_defined_estimators

            if len(estimators) == 0:
                print(self.massif_name, " has only implausible models")
        #  Apply the goodness of fit
        if self.only_models_that_pass_goodness_of_fit_test:
            estimators = [e for e in estimators if self.goodness_of_fit_test(e)]
        if not (self.remove_physically_implausible_models or self.only_models_that_pass_goodness_of_fit_test):
            assert len(estimators) == len(self.models_classes)
        return estimators

    def get_coordinate(self, altitude, year):
        if isinstance(self.altitude_group, DefaultAltitudeGroup):
            if not isinstance(altitude, list):
                coordinate = np.array([year])
            else:
                coordinate = [year] + altitude
        else:
            coordinate = np.array([altitude, year])
        return coordinate

    def _compute_shape_for_reference_altitude(self, estimator):
        coordinate = self.get_coordinate(self.altitude_plot, self.covariate_after)
        gev_params = estimator.margin_function_from_fit.get_params(coordinate, is_transformed=False)
        shape = gev_params.shape
        return shape

    @property
    def has_at_least_one_valid_model(self):
        return len(self.sorted_estimators_with_aic) > 0

    @property
    def model_class_to_estimator_with_finite_aic(self):
        return {type(estimator.margin_model): estimator for estimator in self.sorted_estimators_with_aic}

    @property
    def best_estimator(self):
        if self.has_at_least_one_valid_model:
            best_estimator = self.sorted_estimators_with_aic[0]
            return best_estimator
        else:
            raise ValueError('This object should not have been called because '
                             'has_at_least_one_valid_model={}'.format(self.has_at_least_one_valid_model))

    @property
    def best_margin_model(self):
        return self.best_estimator.margin_model

    @property
    def best_margin_function_from_fit(self):
        return self.best_estimator.margin_function_from_fit

    @property
    def best_shape(self):
        return self.get_gev_params(altitude=self.altitude_plot, year=self.last_year).shape

    @property
    def altitude_plot(self):
        return self.altitude_group.reference_altitude

    def best_coef(self, param_name, dim, degree):
        try:
            coef = self.best_margin_function_from_fit.param_name_to_coef[param_name]  # type: PolynomialAllCoef
            if coef.dim_to_polynomial_coef is None:
                return coef.intercept
            else:
                coef = coef.dim_to_polynomial_coef[dim]  # type: PolynomialCoef
                coef = coef.idx_to_coef[degree]
            return coef
        except (TypeError, KeyError):
            return None

    @property
    def model_names(self):
        return [e.margin_model.name_str for e in self.sorted_estimators_with_aic]

    @property
    def best_name(self):
        name = self.best_estimator.margin_model.name_str
        latex_command = 'textbf' if self.is_significant else 'textrm'
        best_name = '$\\' + latex_command + '{' + name + '}$'
        if self.is_significant:
            best_name = '\\underline{' + best_name + '}'
        return best_name

    # Significant

    @property
    def stationary_estimator(self):
        try:
            if isinstance(self.best_estimator.margin_model, AbstractGumbelAltitudinalModel):
                return self.model_class_to_estimator_with_finite_aic[StationaryGumbelAltitudinal]
            elif isinstance(self.best_estimator.margin_model, AltitudinalOnlyScale):
                return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinalOnlyScale]
            elif isinstance(self.best_estimator.margin_model, AltitudinalShapeLinearTimeStationary):
                return self.model_class_to_estimator_with_finite_aic[AltitudinalShapeLinearTimeStationary]
            elif isinstance(self.best_estimator.margin_model, AltitudinalShapeLinearTimeStationary):
                return self.model_class_to_estimator_with_finite_aic[AltitudinalShapeLinearTimeStationary]
            else:
                if isinstance(self.altitude_group, DefaultAltitudeGroup):
                    return self.model_class_to_estimator_with_finite_aic[StationaryTemporalModel]
                else:
                    return self.model_class_to_estimator_with_finite_aic[StationaryAltitudinal]
        except KeyError:
            raise KeyError('A stationary model must be added either in the list of models, or here in the method')

    @property
    def likelihood_ratio(self):
        return self.stationary_estimator.deviance - self.best_estimator.deviance

    @property
    def degree_freedom_chi2(self):
        return self.best_estimator.nb_params - self.stationary_estimator.nb_params

    @cached_property
    def is_significant(self) -> bool:
        if isinstance(self.altitude_group, DefaultAltitudeGroup):
            # Likelihood ratio based significance
            print("likelihood ratio based significance")
            stationary_model_classes = [StationaryAltitudinal, StationaryGumbelAltitudinal,
                                        AltitudinalShapeLinearTimeStationary]
            if any([isinstance(self.best_estimator.margin_model, c)
                    for c in stationary_model_classes]):
                return False
            else:
                return self.likelihood_ratio > chi2.ppf(q=1 - self.SIGNIFICANCE_LEVEL, df=self.degree_freedom_chi2)
        else:
            # Bootstrap based significance
            return self.cached_results_from_bootstrap[0]

    def sign_of_change(self, function_from_fit):
        return_levels = []
        for temporal_covariate in self._covariate_before_and_after:
            coordinate = np.array([self.altitude_plot, temporal_covariate])
            return_level = function_from_fit.get_params(
                coordinate=coordinate,
                is_transformed=False).return_level(return_period=self.return_period)
            return_levels.append(return_level)
        return 100 * (return_levels[1] - return_levels[0]) / return_levels[0]

    def goodness_of_fit_test(self, estimator):
        quantiles = estimator.sorted_empirical_standard_gumbel_quantiles()
        goodness_of_fit_anderson_test = goodness_of_fit_anderson(quantiles, self.SIGNIFICANCE_LEVEL)
        return goodness_of_fit_anderson_test

    def standard_gumbel_quantiles(self, n=None):
        standard_gumbel_distribution = GevParams(loc=0, scale=1, shape=0)
        if n is None:
            n = len(self.dataset.coordinates)
        standard_gumbel_quantiles = [standard_gumbel_distribution.quantile(i / (n + 1)) for i in range(1, n + 1)]
        return standard_gumbel_quantiles

    def best_confidence_interval(self, altitude, year) -> EurocodeConfidenceIntervalFromExtremes:
        coordinate = self.get_coordinate(altitude, year)
        if self.confidence_interval_based_on_delta_method:
            EurocodeConfidenceIntervalFromExtremes.quantile_level = self.quantile_level
            return EurocodeConfidenceIntervalFromExtremes.from_estimator_extremes(
                estimator_extremes=self.best_estimator,
                ci_method=ConfidenceIntervalMethodFromExtremes.ci_mle,
                coordinate=coordinate)
        else:
            key = (altitude, year)
            mean_estimate = self.cached_results_from_bootstrap[1][key]
            confidence_interval = self.cached_results_from_bootstrap[2][key]
            return EurocodeConfidenceIntervalFromExtremes(mean_estimate, confidence_interval)

    def get_return_level(self, function_from_fit, coordinate):
        return function_from_fit.get_params(coordinate).return_level(self.return_period)

    @cached_property
    def best_residuals(self):
        return self.best_estimator.sorted_empirical_standard_gumbel_quantiles()

    @cached_property
    def cached_results_from_bootstrap(self):
        start = time.time()
        bootstrap_fitted_functions = self.bootstrap_fitted_functions_from_fit
        end1 = time.time()
        duration = str(datetime.timedelta(seconds=end1 - start))
        print('Fit duration', duration)

        # First result - Compute the significance
        sign_of_changes = [self.sign_of_change(f) for f in bootstrap_fitted_functions]
        if self.sign_of_change(self.best_margin_function_from_fit) > 0:
            is_significant = np.quantile(sign_of_changes, self.SIGNIFICANCE_LEVEL) > 0
        else:
            is_significant = np.quantile(sign_of_changes, 1 - self.SIGNIFICANCE_LEVEL) < 0

        # Second result - Compute some dictionary for the return level
        altitude_and_year_to_return_level_mean_estimate = {}
        altitude_and_year_to_return_level_confidence_interval = {}
        altitudes = altitudes_for_groups[self.altitude_group.group_id - 1]
        years = [1959, 2019]
        for year in years:
            for altitude in altitudes:
                key = (altitude, year)
                coordinate = self.get_coordinate(altitude, year)
                mean_estimate = self.get_return_level(self.best_margin_function_from_fit, coordinate)
                bootstrap_return_levels = [self.get_return_level(f, coordinate) for f in
                                           bootstrap_fitted_functions]
                confidence_interval = tuple([np.quantile(bootstrap_return_levels, q)
                                             for q in AbstractExtractEurocodeReturnLevel.bottom_and_upper_quantile])
                altitude_and_year_to_return_level_mean_estimate[key] = mean_estimate
                altitude_and_year_to_return_level_confidence_interval[key] = confidence_interval

        return is_significant, altitude_and_year_to_return_level_mean_estimate, altitude_and_year_to_return_level_confidence_interval

    @property
    def return_level_last_temporal_coordinate(self):
        df_temporal_covariate = self.dataset.coordinates.df_temporal_coordinates_for_fit(
            temporal_covariate_for_fit=self.temporal_covariate_for_fit,
            drop_duplicates=False)
        last_temporal_coordinate = df_temporal_covariate.loc[:, AbstractCoordinates.COORDINATE_T].max()
        # todo: améliorer the last temporal coordinate. on recupère la liste des rights_limits, puis on prend la valeur juste au dessus ou égale."""
        print('last temporal coordinate', last_temporal_coordinate)
        altitude = self.altitude_group.reference_altitude
        coordinate = self.get_coordinate(altitude, last_temporal_coordinate)
        return self.get_return_level(self.best_margin_function_from_fit, coordinate)

    @property
    def bootstrap_fitted_functions_from_fit(self):
        print('nb of bootstrap for confidence interval=', AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP)
        multiprocess = None
        idxs = list(range(AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP))

        if multiprocess is None:
            start = time.time()
            with Pool(NB_CORES) as p:
                batchsize = math.ceil(AbstractExtractEurocodeReturnLevel.NB_BOOTSTRAP / NB_CORES)
                list_functions_from_fit = p.map(self.fit_batch_bootstrap_estimator, batch(idxs, batchsize=batchsize))
                functions_from_fit = list(chain.from_iterable(list_functions_from_fit))

        elif multiprocess:
            print('multiprocessing')
            start = time.time()
            with Pool(NB_CORES) as p:
                functions_from_fit = p.map(self.fit_one_bootstrap_estimator, idxs)

        else:
            start = time.time()

            functions_from_fit = []
            for idx in idxs:
                estimator = self.fit_one_bootstrap_estimator(idx)
                functions_from_fit.append(estimator)

        end1 = time.time()
        duration = str(datetime.timedelta(seconds=end1 - start))
        print('Multiprocessing duration', duration)
        return functions_from_fit

    def fit_batch_bootstrap_estimator(self, idxs):
        list_function_from_fit = []
        for idx in idxs:
            list_function_from_fit.append(self.fit_one_bootstrap_estimator(idx))
        return list_function_from_fit

    def fit_one_bootstrap_estimator(self, idx):
        resample_residuals = resample(self.best_residuals)
        coordinate_values_to_maxima = self.best_estimator. \
            coordinate_values_to_maxima_from_standard_gumbel_quantiles(standard_gumbel_quantiles=resample_residuals)

        coordinates = self.dataset.coordinates
        observations = AnnualMaxima.from_coordinates(coordinates, coordinate_values_to_maxima)
        dataset = AbstractDataset(observations=observations, coordinates=coordinates)
        model_class = type(self.best_margin_model)

        function_from_fit = self.fitted_linear_margin_estimator(model_class, dataset).margin_function_from_fit

        return function_from_fit