[SAFRAN] move safran directory inside experiment directory, add safran...

[SAFRAN] move safran directory inside experiment directory, add safran visualizer. add smooth marginal estimation & full estimation

safran_study/fit_safran.py → experiment/safran_study/fit_safran.py

 import pandas as pd
 
-from safran_study.safran import Safran
-from safran_study.safran_extended import ExtendedSafran
+from experiment.safran_study.safran_extended import ExtendedSafran
 from utils import VERSION
 
 

safran_study/visualize_safran.py → experiment/safran_study/main_visualize_safran.py

-from safran_study.safran import Safran
+from experiment.safran_study.safran import Safran
 from itertools import product
 
+from experiment.safran_study.safran_visualizer import SafranVisualizer
+
 
 def load_all_safran(only_first_one=False):
     all_safran = []
@@ -12,19 +14,9 @@ def load_all_safran(only_first_one=False):
     return all_safran
 
 
-def fit_mle_gev_independent(threshold=None):
-    # Dump the result in a csv
-    dfs = []
-    for safran in load_all_safran(only_first_one=True):
-        safran.visualize_margin_fits_with_cmap(threshold=threshold)
-        # path = r'/home/erwan/Documents/projects/spatiotemporalextremes/local/spatio_temporal_datasets/results/fit_mle_massif/fit_mle_gev_{}.csv'
-        # df.to_csv(path.format(VERSION))
-
-
-def fit_max_stab():
-    pass
-
-
 if __name__ == '__main__':
-    threshold = [None, 20, 40, 60][1]
-    fit_mle_gev_independent(threshold=threshold)
+    for safran in load_all_safran(only_first_one=True):
+        safran_visualizer = SafranVisualizer(safran)
+        # safran_visualizer.visualize_independent_margin_fits(threshold=[None, 20, 40, 60][1])
+        # safran_visualizer.visualize_smooth_margin_fit()
+        safran_visualizer.visualize_full_fit()

safran_study/safran.py → experiment/safran_study/safran.py

 import os
-import matplotlib as mpl
-import matplotlib.colorbar as cbar
+from typing import List
+
 import os.path as op
 from collections import OrderedDict
 
 import matplotlib.pyplot as plt
 import pandas as pd
-from matplotlib import cm
-from mpl_toolkits.axes_grid1 import make_axes_locatable
 from netCDF4 import Dataset
 
-from extreme_estimator.margin_fits.extreme_params import ExtremeParams
-from extreme_estimator.margin_fits.gev.gevmle_fit import GevMleFit
-from extreme_estimator.margin_fits.gev.gev_params import GevParams
-from extreme_estimator.margin_fits.gpd.gpd_params import GpdParams
-from extreme_estimator.margin_fits.gpd.gpdmle_fit import GpdMleFit
-from safran_study.massif import safran_massif_names_from_datasets
-from safran_study.shifted_color_map import shiftedColorMap
-from safran_study.snowfall_annual_maxima import SafranSnowfall
+from experiment.safran_study.massif import safran_massif_names_from_datasets
+from experiment.safran_study.snowfall_annual_maxima import SafranSnowfall
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 from spatio_temporal_dataset.coordinates.spatial_coordinates.abstract_spatial_coordinates import \
     AbstractSpatialCoordinates
+from spatio_temporal_dataset.spatio_temporal_observations.annual_maxima_observations import AnnualMaxima
 from utils import get_full_path, cached_property
 
 
@@ -36,103 +29,6 @@ class Safran(object):
             os.makedirs(self.result_full_path, exist_ok=True)
         df.to_csv(op.join(self.result_full_path, 'merged_array_{}_altitude.csv'.format(self.safran_altitude)))
 
-    """ Visualization methods """
-
-    def visualize(self, ax=None, massif_name_to_fill_kwargs=None, show=True):
-        if ax is None:
-            ax = plt.gca()
-        df_massif = pd.read_csv(op.join(self.map_full_path, 'massifsalpes.csv'))
-        coord_tuples = [(row_massif['idx'], row_massif[AbstractCoordinates.COORDINATE_X],
-                         row_massif[AbstractCoordinates.COORDINATE_Y])
-                        for _, row_massif in df_massif.iterrows()]
-
-        for coordinate_id in set([tuple[0] for tuple in coord_tuples]):
-            l = [coords for idx, *coords in coord_tuples if idx == coordinate_id]
-            l = list(zip(*l))
-            ax.plot(*l, color='black')
-            massif_name = self.coordinate_id_to_massif_name[coordinate_id]
-            fill_kwargs = massif_name_to_fill_kwargs[massif_name] if massif_name_to_fill_kwargs is not None else {}
-            ax.fill(*l, **fill_kwargs)
-        ax.scatter(self.massifs_coordinates.x_coordinates, self.massifs_coordinates.y_coordinates)
-        ax.axis('off')
-
-        if show:
-            plt.show()
-
-    def visualize_margin_fits_with_cmap(self, threshold=None, show=True, axes=None):
-        if threshold is None:
-            params_names = GevParams.SUMMARY_NAMES
-            df = self.df_gev_mle_each_massif
-            # todo: understand how Maurienne could be negative
-            # print(df.head())
-        else:
-            params_names = GpdParams.SUMMARY_NAMES
-            df = self.df_gpd_mle_each_massif(threshold)
-
-        if axes is None:
-            fig, axes = plt.subplots(1, len(params_names))
-            fig.subplots_adjust(hspace=1.0, wspace=1.0)
-
-        for i, gev_param_name in enumerate(params_names):
-            ax = axes[i]
-            massif_name_to_value = df.loc[gev_param_name, :].to_dict()
-            # Compute the middle point of the values for the color map
-            values = list(massif_name_to_value.values())
-            vmin, vmax = min(values), max(values)
-            try:
-                midpoint = 1 - vmax / (vmax + abs(vmin))
-            except ZeroDivisionError:
-                pass
-            # Load the shifted cmap to center on a middle point
-
-            cmap = [plt.cm.coolwarm, plt.cm.bwr, plt.cm.seismic][1]
-            if gev_param_name == ExtremeParams.SHAPE:
-                shifted_cmap = shiftedColorMap(cmap, midpoint=midpoint, name='shifted')
-                norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
-            else:
-                shifted_cmap = shiftedColorMap(cmap, midpoint=0.0, name='shifted')
-                norm = mpl.colors.Normalize(vmin=vmin-1, vmax=vmax)
-
-            m = cm.ScalarMappable(norm=norm, cmap=shifted_cmap)
-
-            massif_name_to_fill_kwargs = {massif_name: {'color': m.to_rgba(value)} for massif_name, value in
-                                          massif_name_to_value.items()}
-
-            self.visualize(ax=ax, massif_name_to_fill_kwargs=massif_name_to_fill_kwargs, show=False)
-
-            # Add colorbar
-            # plt.axis('off')
-            divider = make_axes_locatable(ax)
-            cax = divider.append_axes('right', size='5%', pad=0.05)
-            cb = cbar.ColorbarBase(cax, cmap=shifted_cmap, norm=norm)
-            cb.set_label(gev_param_name)
-
-        if show:
-            plt.show()
-
-    def visualize_cmap(self, massif_name_to_value):
-        orig_cmap = plt.cm.coolwarm
-        # shifted_cmap = shiftedColorMap(orig_cmap, midpoint=0.75, name='shifted')
-
-        massif_name_to_fill_kwargs = {massif_name: {'color': orig_cmap(value)} for massif_name, value in massif_name_to_value.items()}
-
-        self.visualize(massif_name_to_fill_kwargs=massif_name_to_fill_kwargs)
-
-    """ Statistics methods """
-
-    @property
-    def df_gev_mle_each_massif(self):
-        # Fit a margin_fits on each massif
-        massif_to_gev_mle = {massif_name: GevMleFit(self.df_annual_maxima[massif_name]).gev_params.summary_serie
-                             for massif_name in self.safran_massif_names}
-        return pd.DataFrame(massif_to_gev_mle, columns=self.safran_massif_names)
-
-    def df_gpd_mle_each_massif(self, threshold):
-        # Fit a margin fit on each massif
-        massif_to_gev_mle = {massif_name: GpdMleFit(self.df_all_snowfall_concatenated[massif_name], threshold=threshold).gpd_params.summary_serie
-                             for massif_name in self.safran_massif_names}
-        return pd.DataFrame(massif_to_gev_mle, columns=self.safran_massif_names)
-
     """ Data """
 
     @property
@@ -145,6 +41,10 @@ class Safran(object):
     def df_annual_maxima(self):
         return pd.DataFrame(self.year_to_annual_maxima, index=self.safran_massif_names).T
 
+    @property
+    def observations_annual_maxima(self):
+        return AnnualMaxima(df_maxima_gev=self.df_annual_maxima.T)
+
     """ Load some attributes only once """
 
     @cached_property
@@ -166,7 +66,7 @@ class Safran(object):
         return year_to_snowfall
 
     @property
-    def safran_massif_names(self):
+    def safran_massif_names(self) -> List[str]:
         # Load the names of the massif as defined by SAFRAN
         return safran_massif_names_from_datasets(self.year_to_dataset_ordered_dict.values())
 
@@ -196,12 +96,38 @@ class Safran(object):
         df_centroid = pd.read_csv(op.join(self.map_full_path, 'coordonnees_massifs_alpes.csv'))
         # Assert that the massif names are the same between SAFRAN and the coordinate file
         assert not set(self.safran_massif_names).symmetric_difference(set(df_centroid['NOM']))
+        df_centroid.set_index('NOM', inplace=True)
+        df_centroid = df_centroid.loc[self.safran_massif_names]
         return df_centroid
 
     @property
-    def coordinate_id_to_massif_name(self):
+    def coordinate_id_to_massif_name(self) -> dict:
         df_centroid = self.load_df_centroid()
-        return dict(zip(df_centroid['id'], df_centroid['NOM']))
+        return dict(zip(df_centroid['id'], df_centroid.index))
+
+
+    """ Visualization methods """
+
+    def visualize(self, ax=None, massif_name_to_fill_kwargs=None, show=True):
+        if ax is None:
+            ax = plt.gca()
+        df_massif = pd.read_csv(op.join(self.map_full_path, 'massifsalpes.csv'))
+        coord_tuples = [(row_massif['idx'], row_massif[AbstractCoordinates.COORDINATE_X],
+                         row_massif[AbstractCoordinates.COORDINATE_Y])
+                        for _, row_massif in df_massif.iterrows()]
+
+        for coordinate_id in set([tuple[0] for tuple in coord_tuples]):
+            l = [coords for idx, *coords in coord_tuples if idx == coordinate_id]
+            l = list(zip(*l))
+            ax.plot(*l, color='black')
+            massif_name = self.coordinate_id_to_massif_name[coordinate_id]
+            fill_kwargs = massif_name_to_fill_kwargs[massif_name] if massif_name_to_fill_kwargs is not None else {}
+            ax.fill(*l, **fill_kwargs)
+        ax.scatter(self.massifs_coordinates.x_coordinates, self.massifs_coordinates.y_coordinates)
+        ax.axis('off')
+
+        if show:
+            plt.show()
 
     """ Some properties """
 

safran_study/safran_extended.py → experiment/safran_study/safran_extended.py

@@ -2,7 +2,7 @@ from collections import OrderedDict
 
 import pandas as pd
 
-from safran_study.safran import Safran
+from experiment.safran_study.safran import Safran
 from utils import cached_property
 
 

experiment/safran_study/safran_visualizer.py

+import matplotlib as mpl
+import matplotlib.colorbar as cbar
+import matplotlib.pyplot as plt
+import pandas as pd
+from matplotlib import cm
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+from extreme_estimator.estimator.full_estimator.abstract_full_estimator import \
+    FullEstimatorInASingleStepWithSmoothMargin
+from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import SmoothMarginEstimator
+from extreme_estimator.estimator.max_stable_estimator.abstract_max_stable_estimator import MaxStableEstimator
+from extreme_estimator.extreme_models.margin_model.smooth_margin_model import LinearAllParametersAllDimsMarginModel
+from extreme_estimator.extreme_models.max_stable_model.max_stable_models import ExtremalT, Smith
+from extreme_estimator.margin_fits.extreme_params import ExtremeParams
+from extreme_estimator.margin_fits.gev.gev_params import GevParams
+from extreme_estimator.margin_fits.gev.gevmle_fit import GevMleFit
+from extreme_estimator.margin_fits.gpd.gpd_params import GpdParams
+from extreme_estimator.margin_fits.gpd.gpdmle_fit import GpdMleFit
+from experiment.safran_study.safran import Safran
+from experiment.safran_study.shifted_color_map import shiftedColorMap
+from spatio_temporal_dataset.dataset.abstract_dataset import AbstractDataset
+
+
+class SafranVisualizer(object):
+
+    def __init__(self, safran: Safran, show=True):
+        self.safran = safran
+        self.show = show
+
+    @property
+    def observations(self):
+        return self.safran.observations_annual_maxima
+
+    @property
+    def coordinates(self):
+        return self.safran.massifs_coordinates
+
+    @property
+    def dataset(self):
+        return AbstractDataset(self.observations, self.coordinates)
+
+    def visualize_smooth_margin_fit(self):
+        margin_model = LinearAllParametersAllDimsMarginModel(coordinates=self.coordinates)
+        estimator = SmoothMarginEstimator(dataset=self.dataset, margin_model=margin_model)
+        estimator.fit()
+        estimator.margin_function_fitted.visualize(show=self.show)
+
+    def visualize_full_fit(self):
+        max_stable_model = Smith()
+        margin_model = LinearAllParametersAllDimsMarginModel(coordinates=self.coordinates)
+        estimator = FullEstimatorInASingleStepWithSmoothMargin(self.dataset, margin_model, max_stable_model)
+        estimator.fit()
+        estimator.margin_function_fitted.visualize(show=self.show)
+
+    def visualize_independent_margin_fits(self, threshold=None, axes=None):
+        if threshold is None:
+            params_names = GevParams.SUMMARY_NAMES
+            df = self.df_gev_mle_each_massif
+            # todo: understand how Maurienne could be negative
+            # print(df.head())
+        else:
+            params_names = GpdParams.SUMMARY_NAMES
+            df = self.df_gpd_mle_each_massif(threshold)
+
+        if axes is None:
+            fig, axes = plt.subplots(1, len(params_names))
+            fig.subplots_adjust(hspace=1.0, wspace=1.0)
+
+        for i, gev_param_name in enumerate(params_names):
+            ax = axes[i]
+            massif_name_to_value = df.loc[gev_param_name, :].to_dict()
+            # Compute the middle point of the values for the color map
+            values = list(massif_name_to_value.values())
+            vmin, vmax = min(values), max(values)
+            try:
+                midpoint = 1 - vmax / (vmax + abs(vmin))
+            except ZeroDivisionError:
+                pass
+            # Load the shifted cmap to center on a middle point
+
+            cmap = [plt.cm.coolwarm, plt.cm.bwr, plt.cm.seismic][1]
+            if gev_param_name == ExtremeParams.SHAPE:
+                shifted_cmap = shiftedColorMap(cmap, midpoint=midpoint, name='shifted')
+                norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
+            else:
+                shifted_cmap = shiftedColorMap(cmap, midpoint=0.0, name='shifted')
+                norm = mpl.colors.Normalize(vmin=vmin - 1, vmax=vmax)
+
+            m = cm.ScalarMappable(norm=norm, cmap=shifted_cmap)
+
+            massif_name_to_fill_kwargs = {massif_name: {'color': m.to_rgba(value)} for massif_name, value in
+                                          massif_name_to_value.items()}
+
+            self.safran.visualize(ax=ax, massif_name_to_fill_kwargs=massif_name_to_fill_kwargs, show=False)
+
+            # Add colorbar
+            # plt.axis('off')
+            divider = make_axes_locatable(ax)
+            cax = divider.append_axes('right', size='5%', pad=0.05)
+            cb = cbar.ColorbarBase(cax, cmap=shifted_cmap, norm=norm)
+            cb.set_label(gev_param_name)
+
+        if self.show:
+            plt.show()
+
+    def visualize_cmap(self, massif_name_to_value):
+        orig_cmap = plt.cm.coolwarm
+        # shifted_cmap = shiftedColorMap(orig_cmap, midpoint=0.75, name='shifted')
+
+        massif_name_to_fill_kwargs = {massif_name: {'color': orig_cmap(value)} for massif_name, value in
+                                      massif_name_to_value.items()}
+
+        self.safran.visualize(massif_name_to_fill_kwargs=massif_name_to_fill_kwargs)
+
+    """ Statistics methods """
+
+    @property
+    def df_gev_mle_each_massif(self):
+        # Fit a margin_fits on each massif
+        massif_to_gev_mle = {massif_name: GevMleFit(self.safran.df_annual_maxima[massif_name]).gev_params.summary_serie
+                             for massif_name in self.safran.safran_massif_names}
+        return pd.DataFrame(massif_to_gev_mle, columns=self.safran.safran_massif_names)
+
+    def df_gpd_mle_each_massif(self, threshold):
+        # Fit a margin fit on each massif
+        massif_to_gev_mle = {massif_name: GpdMleFit(self.safran.df_all_snowfall_concatenated[massif_name],
+                                                    threshold=threshold).gpd_params.summary_serie
+                             for massif_name in self.safran.safran_massif_names}
+        return pd.DataFrame(massif_to_gev_mle, columns=self.safran.safran_massif_names)

spatio_temporal_dataset/dataset/abstract_dataset.py

@@ -16,7 +16,7 @@ from spatio_temporal_dataset.spatio_temporal_observations.abstract_spatio_tempor
 class AbstractDataset(object):
 
     def __init__(self, observations: AbstractSpatioTemporalObservations, coordinates: AbstractCoordinates):
-        assert pd.Index.equals(observations.index, coordinates.index)
+        assert pd.Index.equals(observations.index, coordinates.index), '\n{}\n{}'.format(observations.index, coordinates.index)
         self.observations = observations  # type: AbstractSpatioTemporalObservations
         self.coordinates = coordinates  # type: AbstractCoordinates
         self.subset_id_to_column_idxs = None  # type: Dict[int, List[int]]

test/test_experiment/test_safran/test_safran_visualizer.py

+# import unittest
+#
+# from experiment.safran_study.safran import Safran
+# from experiment.safran_study.safran_visualizer import SafranVisualizer
+#
+#
+# class TestSafranVisualizer(unittest.TestCase):
+#     DISPLAY = False
+#
+#     def setUp(self):
+#         super().setUp()
+#         self.safran = Safran(1800, 1)
+#         self.safran_visualizer = SafranVisualizer(self.safran, show=self.DISPLAY)
+#
+#     def tearDown(self) -> None:
+#         self.assertTrue(True)
+#
+#     def test_safran_smooth_margin_estimator(self):
+#         self.safran_visualizer.visualize_smooth_margin_fit()
+#
+#     def test_safran_independent_margin_fits(self):
+#         self.safran_visualizer.visualize_independent_margin_fits()
+#
+#     def test_safran_full_estimator(self):
+#         self.safran_visualizer.visualize_full_fit()
+#
+#
+# if __name__ == '__main__':
+#     unittest.main()

test/test_utils.py

@@ -7,7 +7,7 @@ from extreme_estimator.extreme_models.max_stable_model.abstract_max_stable_model
     AbstractMaxStableModelWithCovarianceFunction, CovarianceFunction
 from extreme_estimator.extreme_models.max_stable_model.max_stable_models import Smith, BrownResnick, Schlather, \
     Geometric, ExtremalT, ISchlather
-from safran_study.safran import Safran
+from experiment.safran_study.safran import Safran
 from spatio_temporal_dataset.coordinates.spatial_coordinates.alps_station_3D_coordinates import \
     AlpsStation3DCoordinatesWithAnisotropy
 from spatio_temporal_dataset.coordinates.spatial_coordinates.generated_spatial_coordinates import \