Commit 384abbcd authored by Le Roux Erwan's avatar Le Roux Erwan
Browse files

[MARGIN FUNCTION] add full visualization for 1D/2D margin functions

parent d130ce71
No related merge requests found
Showing with 160 additions and 55 deletions
+160 -55
import numpy as np
from extreme_estimator.estimator.full_estimator import FullEstimatorInASingleStepWithSmoothMargin
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeAxis0MarginModel, \
LinearAllParametersAllAxisMarginModel
from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith
from extreme_estimator.gev_params import GevParams
from spatio_temporal_dataset.coordinates.unidimensional_coordinates.coordinates_1D import LinSpaceCoordinates
import matplotlib.pyplot as plt
from spatio_temporal_dataset.dataset.simulation_dataset import FullSimulatedDataset
nb_points = 50
nb_obs = 100
coordinates = LinSpaceCoordinates.from_nb_points(nb_points=nb_points)
########## GENERATING THE DATA #####################
# MarginModel Linear with respect to the shape (from 0.01 to 0.02)
margin_model = LinearShapeAxis0MarginModel(coordinates=coordinates, params_sample={GevParams.GEV_SHAPE: 0.02})
max_stable_model = Smith()
dataset = FullSimulatedDataset.from_double_sampling(nb_obs=nb_obs, margin_model=margin_model,
coordinates=coordinates,
max_stable_model=max_stable_model)
# Visualize the sampling margin
dataset.margin_model.margin_function_sample.visualize_all()
# Plot a realization from the maxima gev (i.e the maxima obtained just by simulating the marginal law)
for maxima_gev in np.transpose(dataset.maxima_gev):
plt.plot(coordinates.coordinates_values, maxima_gev)
plt.show()
######### FITTING A MODEL #################
margin_model = LinearAllParametersAllAxisMarginModel(coordinates)
full_estimator = FullEstimatorInASingleStepWithSmoothMargin(dataset, margin_model, max_stable_model)
full_estimator.fit()
print(full_estimator.full_params_fitted)
# Plot the margin functions
# margin_model.margin_function_sample.visualize_2D()
import unittest
from extreme_estimator.gev_params import GevParams
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeAxis0MarginModel, \
LinearAllParametersAllAxisMarginModel
from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import CircleCoordinates
from spatio_temporal_dataset.coordinates.unidimensional_coordinates.coordinates_1D import LinSpaceCoordinates
class VisualizationMarginModel(unittest.TestCase):
DISPLAY = True
nb_points = 50
margin_model = [LinearShapeAxis0MarginModel, LinearAllParametersAllAxisMarginModel][-1]
@classmethod
def example_visualization_2D(cls):
spatial_coordinates = CircleCoordinates.from_nb_points(nb_points=cls.nb_points)
margin_model = cls.margin_model(coordinates=spatial_coordinates)
if cls.DISPLAY:
margin_model.margin_function_sample.visualize_all()
@classmethod
def example_visualization_1D(cls):
coordinates = LinSpaceCoordinates.from_nb_points(nb_points=cls.nb_points)
# MarginModel Linear with respect to the shape (from 0.01 to 0.02)
margin_model = cls.margin_model(coordinates=coordinates, params_sample={GevParams.GEV_SHAPE: 0.02})
if cls.DISPLAY:
margin_model.margin_function_sample.visualize_all()
if __name__ == '__main__':
VisualizationMarginModel.example_visualization_1D()
VisualizationMarginModel.example_visualization_2D()
......@@ -25,7 +25,7 @@ class SmoothMarginalsThenUnitaryMsp(AbstractFullEstimator):
self.margin_estimator.fit()
# Compute the maxima_frech
maxima_frech = AbstractMarginModel.gev2frech(maxima_gev=self.dataset.maxima_gev,
coordinates=self.dataset.coordinates_values,
coordinates_values=self.dataset.coordinates_values,
margin_function=self.margin_estimator.margin_function_fitted)
# Update maxima frech field through the dataset object
self.dataset.maxima_frech = maxima_frech
......@@ -47,15 +47,15 @@ class FullEstimatorInASingleStepWithSmoothMargin(AbstractFullEstimator):
self.smooth_margin_function_to_fit = margin_model.margin_function_start_fit
def _fit(self):
# todo: specify the shape of the margin
# Estimate the margin
self.max_stable_params_fitted = self.max_stable_model.fitmaxstab(
# Estimate both the margin and the max-stable structure
self.full_params_fitted = self.max_stable_model.fitmaxstab(
maxima_frech=self.dataset.maxima_frech,
df_coordinates=self.dataset.df_coordinates,
fit_marge=True,
fit_marge_form_dict=self.smooth_margin_function_to_fit.fit_marge_form_dict,
margin_start_dict=self.smooth_margin_function_to_fit.margin_start_dict
)
# Initialize
class PointwiseAndThenUnitaryMsp(AbstractFullEstimator):
......
......@@ -33,4 +33,4 @@ class SmoothMarginEstimator(AbstractMarginEstimator):
def _fit(self):
self._margin_function_fitted = self.margin_model.fitmargin_from_maxima_gev(maxima_gev=self.dataset.maxima_gev,
coordinates=self.dataset.coordinates_values)
coordinates_values=self.dataset.coordinates_values)
import numpy as np
from extreme_estimator.extreme_models.abstract_model import AbstractModel
from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function import AbstractMarginFunction
from extreme_estimator.extreme_models.margin_model.margin_function.abstract_margin_function \
import AbstractMarginFunction
from extreme_estimator.gev_params import GevParams
from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
......@@ -27,34 +28,36 @@ class AbstractMarginModel(AbstractModel):
# Conversion class methods
@classmethod
def convert_maxima(cls, convertion_r_function, maxima: np.ndarray, coordinates: np.ndarray,
margin_function: AbstractMarginFunction):
assert isinstance(coordinates, np.ndarray)
assert len(maxima) == len(coordinates)
def convert_maxima(cls, convertion_r_function, maxima: np.ndarray, coordinates_values: np.ndarray,
margin_function: AbstractMarginFunction) -> np.ndarray:
assert isinstance(coordinates_values, np.ndarray)
assert len(maxima) == len(coordinates_values)
converted_maxima = []
for x, coordinate in zip(maxima, coordinates):
for x, coordinate in zip(maxima, coordinates_values):
gev_params = margin_function.get_gev_params(coordinate)
x_gev = convertion_r_function(x, **gev_params.to_dict())
converted_maxima.append(x_gev)
return np.array(converted_maxima)
@classmethod
def gev2frech(cls, maxima_gev: np.ndarray, coordinates: np.ndarray, margin_function: AbstractMarginFunction):
return cls.convert_maxima(cls.r.gev2frech, maxima_gev, coordinates, margin_function)
def gev2frech(cls, maxima_gev: np.ndarray, coordinates_values: np.ndarray,
margin_function: AbstractMarginFunction) -> np.ndarray:
return cls.convert_maxima(cls.r.gev2frech, maxima_gev, coordinates_values, margin_function)
@classmethod
def frech2gev(cls, maxima_frech: np.ndarray, coordinates: np.ndarray, margin_function: AbstractMarginFunction):
return cls.convert_maxima(cls.r.frech2gev, maxima_frech, coordinates, margin_function)
def frech2gev(cls, maxima_frech: np.ndarray, coordinates_values: np.ndarray,
margin_function: AbstractMarginFunction) -> np.ndarray:
return cls.convert_maxima(cls.r.frech2gev, maxima_frech, coordinates_values, margin_function)
# Sampling methods
def rmargin_from_maxima_frech(self, maxima_frech: np.ndarray, coordinates: np.ndarray):
maxima_gev = self.frech2gev(maxima_frech, coordinates, self.margin_function_sample)
def rmargin_from_maxima_frech(self, maxima_frech: np.ndarray, coordinates_values: np.ndarray) -> np.ndarray:
maxima_gev = self.frech2gev(maxima_frech, coordinates_values, self.margin_function_sample)
return maxima_gev
def rmargin_from_nb_obs(self, nb_obs, coordinates):
def rmargin_from_nb_obs(self, nb_obs: int, coordinates_values: np.ndarray) -> np.ndarray:
maxima_gev = []
for coordinate in coordinates:
for coordinate in coordinates_values:
gev_params = self.margin_function_sample.get_gev_params(coordinate)
x_gev = self.r.rgev(nb_obs, **gev_params.to_dict())
maxima_gev.append(x_gev)
......@@ -62,5 +65,6 @@ class AbstractMarginModel(AbstractModel):
# Fitting methods needs to be defined in child classes
def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates_values: np.ndarray) \
-> AbstractMarginFunction:
pass
......@@ -14,12 +14,40 @@ class AbstractMarginFunction(object):
self.default_params = default_params.to_dict()
def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
"""Main function that maps each coordinate to its GEV parameters"""
"""Main method that maps each coordinate to its GEV parameters"""
pass
# Visualization function
def visualize_2D(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=False):
def visualize_all(self):
fig, axes = plt.subplots(3, 1, sharex='col', sharey='row')
fig.subplots_adjust(hspace=0.4, wspace=0.4, )
for i, gev_param_name in enumerate(GevParams.GEV_PARAM_NAMES):
ax = axes[i]
self.visualize(gev_param_name, ax, show=False)
title_str = gev_param_name
ax.set_title(title_str)
plt.show()
def visualize(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=True):
if self.coordinates.nb_columns == 1:
self.visualize_1D(gev_param_name, ax, show)
elif self.coordinates.nb_columns == 2:
self.visualize_2D(gev_param_name, ax, show)
else:
raise NotImplementedError('3D Margin visualization not yet implemented')
def visualize_1D(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=True):
x = self.coordinates.x_coordinates
grid, linspace = self.get_grid_1D(x)
gev_param_idx = GevParams.GEV_PARAM_NAMES.index(gev_param_name)
if ax is None:
ax = plt.gca()
ax.plot(linspace, grid[:, gev_param_idx])
if show:
plt.show()
def visualize_2D(self, gev_param_name=GevParams.GEV_LOC, ax=None, show=True):
x = self.coordinates.x_coordinates
y = self.coordinates.y_coordinates
grid = self.get_grid_2D(x, y)
......@@ -32,6 +60,15 @@ class AbstractMarginFunction(object):
if show:
plt.show()
def get_grid_1D(self, x):
# TODO: to avoid getting the value several times, I could cache the results
resolution = 100
grid = np.zeros([resolution, 3])
linspace = np.linspace(x.min(), x.max(), resolution)
for i, xi in enumerate(linspace):
grid[i] = self.get_gev_params(np.array([xi])).to_array()
return grid, linspace
def get_grid_2D(self, x, y):
resolution = 100
grid = np.zeros([resolution, resolution, 3])
......
......@@ -65,7 +65,8 @@ class LinearMarginFunction(IndependentMarginFunction):
start=self.default_params[gev_param_name])
# Some check on the Linear param function
if gev_param_name == GevParams.GEV_SCALE:
assert param_function.end > param_function.start, 'Impossible linear rate for Scale parameter'
assert param_function.end > 0 and param_function.start > 0, \
'Impossible start/end value for Scale parameter'
# Add the param_function to the dictionary
self.gev_param_name_to_param_function[gev_param_name] = param_function
......
......@@ -20,7 +20,7 @@ class ConstantParamFunction(ParamFunction):
class LinearOneAxisParamFunction(ParamFunction):
def __init__(self, linear_axis: int, coordinates_axis: np.ndarray, start: float, end: float = 2.0):
def __init__(self, linear_axis: int, coordinates_axis: np.ndarray, start: float, end: float = 0.01):
self.linear_axis = linear_axis
self.t_min = coordinates_axis.min()
self.t_max = coordinates_axis.max()
......@@ -36,7 +36,7 @@ class LinearOneAxisParamFunction(ParamFunction):
class LinearParamFunction(ParamFunction):
def __init__(self, linear_axes: List[int], coordinates: np.ndarray, start: float, end: float = 2.0):
def __init__(self, linear_axes: List[int], coordinates: np.ndarray, start: float, end: float = 0.01):
self.linear_one_axis_param_functions = [] # type: List[LinearOneAxisParamFunction]
self.start = start
self.end = end
......
......@@ -20,7 +20,7 @@ class LinearMarginModel(AbstractMarginModel):
default_params=GevParams.from_dict(self.params_start_fit),
gev_param_name_to_linear_axes=gev_param_name_to_linear_axis)
def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates: np.ndarray) -> AbstractMarginFunction:
def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, coordinates_values: np.ndarray) -> AbstractMarginFunction:
return self.margin_function_start_fit
......
......@@ -21,7 +21,7 @@ class AbstractMaxStableModel(AbstractModel):
return {'cov.mod': self.cov_mod}
def fitmaxstab(self, maxima_frech: np.ndarray, df_coordinates: pd.DataFrame, fit_marge=False,
fit_marge_form_dict=None, margin_start_dict=None):
fit_marge_form_dict=None, margin_start_dict=None) -> dict:
assert isinstance(maxima_frech, np.ndarray)
assert isinstance(df_coordinates, pd.DataFrame)
if fit_marge:
......
......@@ -10,6 +10,15 @@ class AbstractUniDimensionalCoordinates(AbstractCoordinates):
pass
class LinSpaceCoordinates(AbstractUniDimensionalCoordinates):
@classmethod
def from_nb_points(cls, nb_points, start=-1.0, end=1.0):
axis_coordinates = np.linspace(start, end, nb_points)
df = pd.DataFrame.from_dict({cls.COORDINATE_X: axis_coordinates})
return cls.from_df(df)
class UniformCoordinates(AbstractUniDimensionalCoordinates):
@classmethod
......
......@@ -19,8 +19,8 @@ class SimulatedDataset(AbstractDataset):
margin_model: AbstractMarginModel = None):
super().__init__(temporal_observations, coordinates)
assert margin_model is not None or max_stable_model is not None
self.margin_model = margin_model
self.max_stable_model = max_stable_model
self.margin_model = margin_model # type: AbstractMarginModel
self.max_stable_model = max_stable_model # type: AbstractMaxStableModel
class MaxStableDataset(SimulatedDataset):
......@@ -47,4 +47,5 @@ class FullSimulatedDataset(SimulatedDataset):
margin_model: AbstractMarginModel):
temporal_obs = FullAnnualMaxima.from_double_sampling(nb_obs, max_stable_model,
coordinates, margin_model)
return cls(temporal_obs, coordinates, max_stable_model)
return cls(temporal_observations=temporal_obs, coordinates=coordinates, max_stable_model=max_stable_model,
margin_model=margin_model)
......@@ -19,7 +19,7 @@ class MarginAnnualMaxima(AnnualMaxima):
@classmethod
def from_sampling(cls, nb_obs: int, coordinates: AbstractCoordinates,
margin_model: AbstractMarginModel):
maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs, coordinates=coordinates.coordinates_values)
maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=nb_obs, coordinates_values=coordinates.coordinates_values)
df_maxima_gev = pd.DataFrame(data=maxima_gev, index=coordinates.index)
return cls(df_maxima_gev=df_maxima_gev)
......@@ -41,6 +41,6 @@ class FullAnnualMaxima(MaxStableAnnualMaxima):
max_stable_annual_maxima = super().from_sampling(nb_obs, max_stable_model, coordinates)
# Compute df_maxima_gev from df_maxima_frech
maxima_gev = margin_model.rmargin_from_maxima_frech(maxima_frech=max_stable_annual_maxima.maxima_frech,
coordinates=coordinates.coordinates_values)
coordinates_values=coordinates.coordinates_values)
max_stable_annual_maxima.df_maxima_gev = pd.DataFrame(data=maxima_gev, index=coordinates.index)
return max_stable_annual_maxima
import unittest
from extreme_estimator.gev_params import GevParams
from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearShapeAxis0MarginModel
from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import CircleCoordinates
class TestLinearMarginModel(unittest.TestCase):
DISPLAY = False
def test_visualization_2D(self):
spatial_coordinates = CircleCoordinates.from_nb_points(nb_points=50)
margin_model = LinearShapeAxis0MarginModel(coordinates=spatial_coordinates)
for gev_param_name in GevParams.GEV_PARAM_NAMES:
margin_model.margin_function_sample.visualize_2D(gev_param_name=gev_param_name, show=self.DISPLAY)
# maxima_gev = margin_model.rmargin_from_nb_obs(nb_obs=10, coordinates=coordinates)
# fitted_margin_function = margin_model.fitmargin_from_maxima_gev(maxima_gev=maxima_gev, coordinates=coordinates)
if __name__ == '__main__':
unittest.main()
import unittest
from extreme_estimator.estimator.margin_estimator import SmoothMarginEstimator
from extreme_estimator.return_level_plot.spatial_2D_plot import Spatial2DPlot
from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import CircleCoordinates
from experiment.return_level_plot.spatial_2D_plot import Spatial2DPlot
from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
from test.test_utils import load_smooth_margin_models, load_test_1D_and_2D_coordinates
......
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