modify r import and add set_seed_r function

extreme_estimator/extreme_models/abstract_model.py

-from extreme_estimator.extreme_models.utils import get_loaded_r
-
-
 class AbstractModel(object):
-    r = get_loaded_r()
 
     def __init__(self, params_start_fit=None, params_sample=None):
         self.default_params_start_fit = None

extreme_estimator/extreme_models/margin_model/abstract_margin_model.py

@@ -3,6 +3,7 @@ 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.utils import r
 from extreme_estimator.gev_params import GevParams
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
@@ -42,12 +43,12 @@ class AbstractMarginModel(AbstractModel):
     @classmethod
     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)
+        return cls.convert_maxima(r.gev2frech, maxima_gev, coordinates_values, margin_function)
 
     @classmethod
     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)
+        return cls.convert_maxima(r.frech2gev, maxima_frech, coordinates_values, margin_function)
 
     # Sampling methods
 
@@ -59,7 +60,7 @@ class AbstractMarginModel(AbstractModel):
         maxima_gev = []
         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())
+            x_gev = r.rgev(nb_obs, **gev_params.to_dict())
             maxima_gev.append(x_gev)
         return np.array(maxima_gev)
 

extreme_estimator/extreme_models/margin_model/margin_function/abstract_margin_function.py

@@ -12,6 +12,7 @@ class AbstractMarginFunction(object):
     def __init__(self, coordinates: AbstractCoordinates):
         self.coordinates = coordinates
         self.visualization_axes = None
+        self.dot_display = False
 
     def get_gev_params(self, coordinate: np.ndarray) -> GevParams:
         """Main method that maps each coordinate to its GEV parameters"""
@@ -19,7 +20,8 @@ class AbstractMarginFunction(object):
 
     # Visualization function
 
-    def visualize(self, axes=None, show=True):
+    def visualize(self, axes=None, show=True, dot_display=False):
+        self.dot_display = dot_display
         if axes is None:
             fig, axes = plt.subplots(3, 1, sharex='col', sharey='row')
             fig.subplots_adjust(hspace=0.4, wspace=0.4, )
@@ -46,7 +48,11 @@ class AbstractMarginFunction(object):
         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 self.dot_display:
+            ax.plot(linspace, grid[:, gev_param_idx], 'o')
+        else:
+            ax.plot(linspace, grid[:, gev_param_idx])
+
         if show:
             plt.show()
 
@@ -65,9 +71,15 @@ class AbstractMarginFunction(object):
 
     def get_grid_1D(self, x):
         # TODO: to avoid getting the value several times, I could cache the results
-        resolution = 100
+        if self.dot_display:
+            resolution = len(self.coordinates)
+            linspace = self.coordinates.coordinates_values[:, 0]
+            print('dot display')
+        else:
+            resolution = 100
+            linspace = np.linspace(x.min(), x.max(), resolution)
+
         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

extreme_estimator/extreme_models/margin_model/smooth_margin_model.py

@@ -62,6 +62,12 @@ class LinearShapeDim1MarginModel(LinearMarginModel):
         super().load_margin_functions({GevParams.GEV_SHAPE: [1]})
 
 
+class LinearScaleDim1MarginModel(LinearMarginModel):
+
+    def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):
+        super().load_margin_functions({GevParams.GEV_SCALE: [1]})
+
+
 class LinearShapeDim1and2MarginModel(LinearMarginModel):
 
     def load_margin_functions(self, margin_function_class: type = None, gev_param_name_to_linear_dims=None):

extreme_estimator/extreme_models/max_stable_model/abstract_max_stable_model.py

@@ -7,6 +7,7 @@ from rpy2.rinterface import RRuntimeError
 import rpy2.robjects as robjects
 
 from extreme_estimator.extreme_models.abstract_model import AbstractModel
+from extreme_estimator.extreme_models.utils import r
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
@@ -53,12 +54,12 @@ class AbstractMaxStableModel(AbstractModel):
             fit_params.update({k: robjects.Formula(v) for k, v in fit_marge_form_dict.items()})
         if fitmaxstab_with_one_dimensional_data:
             fit_params['iso'] = True
-        fit_params['start'] = self.r.list(**start_dict)
+        fit_params['start'] = r.list(**start_dict)
         fit_params['fit.marge'] = fit_marge
 
         # Run the fitmaxstab in R
         try:
-            res = self.r.fitmaxstab(data=data, coord=coord, **fit_params)  # type: robjects.ListVector
+            res = r.fitmaxstab(data=data, coord=coord, **fit_params)  # type: robjects.ListVector
         except RRuntimeError as error:
             raise Exception('Some R exception have been launched at RunTime: \n {}'.format(error.__repr__()))
         # todo: maybe if the convergence was not successful I could try other starting point several times
@@ -72,7 +73,7 @@ class AbstractMaxStableModel(AbstractModel):
         Return an numpy of maxima. With rows being the stations and columns being the years of maxima
         """
         maxima_frech = np.array(
-            self.r.rmaxstab(nb_obs, coordinates, *list(self.cov_mod_param.values()), **self.params_sample))
+            r.rmaxstab(nb_obs, coordinates, *list(self.cov_mod_param.values()), **self.params_sample))
         return np.transpose(maxima_frech)
 
     def remove_unused_parameters(self, start_dict, coordinate_dim):
@@ -95,5 +96,5 @@ class AbstractMaxStableModelWithCovarianceFunction(AbstractMaxStableModel):
         self.default_params_sample = {
             'range': 3,
             'smooth': 0.5,
-            'nugget': 0.5
+            'nugget': 0.0
         }

extreme_estimator/extreme_models/max_stable_model/max_stable_fit.py

@@ -3,15 +3,15 @@ import pandas as pd
 
 import numpy as np
 
-from extreme_estimator.extreme_models.utils import get_loaded_r
+from extreme_estimator.extreme_models.utils import R
 from extreme_estimator.gev.gev_mle_fit import GevMleFit
 import rpy2.robjects.numpy2ri as rpyn
 
-
 import rpy2.robjects as robjects
 
+
 def max_stable_fit():
-    r = get_loaded_r()
+    r = R().r
     r("""
     set.seed(42)
     n.obs = 50
@@ -45,8 +45,6 @@ def max_stable_fit():
 
     # res = r.fitmaxstab(data=data, coord=coord, covmod="gauss", start=namedlist, fit.marge=TRUE, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form)
 
-
-
     # m2.colnames = R.StrVector(['x', 'y'])
     # res = r.fitmaxstab()
     # mle_params = dict(r.attr(res, 'coef').items())

extreme_estimator/extreme_models/utils.py

@@ -3,19 +3,18 @@ import random
 import sys
 
 import rpy2.robjects as ro
+
 from rpy2.robjects import numpy2ri
 from rpy2.robjects import pandas2ri
 
+r = ro.R()
+numpy2ri.activate()
+pandas2ri.activate()
+r.library('SpatialExtremes')
+
 
-def get_loaded_r() -> ro.R:
-    r = ro.r
-    numpy2ri.activate()
-    pandas2ri.activate()
-    r.library('SpatialExtremes')
-    # max_int = r('.Machine$integer.max')
-    # seed = random.randrange(max_int)
-    # r("set.seed({})".format(seed))
-    return r
+def set_seed_r(seed=42):
+    r("set.seed({})".format(seed))
 
 
 def get_associated_r_file(python_filepath: str) -> str:

spatio_temporal_dataset/coordinates/spatial_coordinates/generated_spatial_coordinates.py

 import math
 import numpy as np
 import pandas as pd
+from rpy2.robjects import r
 
-from extreme_estimator.extreme_models.utils import get_loaded_r
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 import matplotlib.pyplot as plt
 
@@ -12,7 +12,6 @@ class CircleCoordinates(AbstractCoordinates):
     @classmethod
     def from_nb_points(cls, nb_points, max_radius=1.0):
         # Sample uniformly inside the circle
-        r = get_loaded_r()
         angles = np.array(r.runif(nb_points, max=2 * math.pi))
         radius = np.sqrt(np.array(r.runif(nb_points, max=max_radius)))
         df = pd.DataFrame.from_dict({cls.COORDINATE_X: radius * np.cos(angles),

spatio_temporal_dataset/coordinates/unidimensional_coordinates/coordinates_1D.py

 import pandas as pd
 
 import numpy as np
+from rpy2.robjects import r
 
-from extreme_estimator.extreme_models.utils import get_loaded_r
 from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoordinates
 
 
@@ -24,7 +24,6 @@ class UniformCoordinates(AbstractUniDimensionalCoordinates):
     @classmethod
     def from_nb_points(cls, nb_points, start=-1.0, end=1.0):
         # Sample uniformly inside the circle
-        r = get_loaded_r()
         axis_coordinates = np.array(r.runif(nb_points, min=start, max=end))
         df = pd.DataFrame.from_dict({cls.COORDINATE_X: axis_coordinates})
         return cls.from_df(df)

test/test_extreme_estimator/test_gev/test_gev_mle_fit.py

@@ -2,16 +2,15 @@ import unittest
 
 import numpy as np
 
+from extreme_estimator.extreme_models.utils import r, set_seed_r
 from extreme_estimator.gev.gev_mle_fit import GevMleFit
-from extreme_estimator.extreme_models.utils import get_loaded_r
 
 
 class TestGevMleFit(unittest.TestCase):
 
     def test_unitary_gev_mle_fit(self):
-        r = get_loaded_r()
+        set_seed_r()
         r("""
-        set.seed(42)
         N <- 50
         loc = 0; scale = 1; shape <- 1
         x_gev <- rgev(N, loc = loc, scale = scale, shape = shape)