diff --git a/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py b/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py
index 0b3ebd762ff3d89474e8e9ce99970fae7fd6ef3e..d79fcec2e475ab3af8f594fb242ceb366a32fa3b 100644
--- a/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py
+++ b/extreme_estimator/estimator/margin_estimator/abstract_margin_estimator.py
@@ -35,9 +35,8 @@ class SmoothMarginEstimator(AbstractMarginEstimator):
def _fit(self):
maxima_gev = self.dataset.maxima_gev(split=self.train_split)
- df_coordinates = self.dataset.df_coordinates(self.train_split)
- df_coordinates_spatial = df_coordinates.loc[:, self.dataset.coordinates.coordinates_spatial_names]
- df_coordinates_temporal = df_coordinates.loc[:, self.dataset.coordinates.coordinates_temporal_names]
+ df_coordinates_spatial = self.dataset.coordinates.df_spatial_coordinates(self.train_split)
+ df_coordinates_temporal = self.dataset.coordinates.df_temporal_coordinates(self.train_split)
self._params_fitted = self.margin_model.fitmargin_from_maxima_gev(maxima_gev=maxima_gev,
df_coordinates_spatial=df_coordinates_spatial,
df_coordinates_temporal=df_coordinates_temporal)
diff --git a/extreme_estimator/extreme_models/margin_model/fitspatgev.R b/extreme_estimator/extreme_models/margin_model/fitspatgev.R
index 6f7c1f7e2d1daf9fce75a1e1df8964fadb0184ca..f8c7d264259775a5956c144f0103ad391d3ef61a 100644
--- a/extreme_estimator/extreme_models/margin_model/fitspatgev.R
+++ b/extreme_estimator/extreme_models/margin_model/fitspatgev.R
@@ -26,7 +26,7 @@ fitspatgev_2D_test <- function (n.obs){
}
# fitspatgev_test with 3D data
-fitspatgev_3D_test <- function (n.obs){
+fitspatgev_3D_test_ok <- function (n.obs){
n.site = 2
covariables <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
colnames(covariables) = c("E", "N")
@@ -47,183 +47,139 @@ fitspatgev_3D_test <- function (n.obs){
temp.cov = matrix(1:n.obs, ncol=1)
colnames(temp.cov) = c("T")
- # namedlist = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0,
- # tempCoeffLoc1=0.0, tempCoeffLoc2=1.0, tempCoeffScale1=0.1, tempCoeffShape1=1.0)
- start = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0,
- tempCoeffLoc1=0.0, tempCoeffLoc2=1.0, tempCoeffScale1=0.1, tempCoeffShape1=1.0)
-
- # res = fitspatgev(data=data, covariables=covariables, start=namedlist, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form,
- # temp.cov=temp.cov, temp.loc.form = temp_loc.form, temp.scale.form = temp_scale.form, temp.shape.form = temp_shape.form)
-
-
-
start = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0,
tempCoeffLoc1=0.0, tempCoeffScale1=0.1, tempCoeffShape1=1.0)
res = fitspatgev(data=data, covariables=covariables, start=start, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form,
temp.cov=temp.cov, temp.form.loc=temp.form.loc, temp.form.scale=temp.form.scale, temp.form.shape=temp.form.shape)
- #
- # # START PASTING CODE FROM THE FUNCTION
- # n.site <- ncol(data)
- # n.obs <- nrow(data)
- # if (n.site != nrow(covariables))
- # stop("'data' and 'covariates' doesn't match")
- # use.temp.cov <- c(!is.null(temp.form.loc), !is.null(temp.form.scale),
- # !is.null(temp.form.shape))
- # if (any(use.temp.cov) && (n.obs != nrow(temp.cov)))
- # stop("'data' and 'temp.cov' doesn't match")
- # if (any(use.temp.cov) && is.null(temp.cov))
- # stop("'temp.cov' must be supplied if at least one temporal formula is given")
- # loc.form <- update(loc.form, y ~ .)
- # scale.form <- update(scale.form, y ~ .)
- # shape.form <- update(shape.form, y ~ .)
- # if (use.temp.cov[1])
- # temp.form.loc <- update(temp.form.loc, y ~ . + 0)
- # if (use.temp.cov[2])
- # temp.form.scale <- update(temp.form.scale, y ~ . + 0)
- # if (use.temp.cov[3])
- # temp.form.shape <- update(temp.form.shape, y ~ . + 0)
- # loc.model <- modeldef(covariables, loc.form)
- # scale.model <- modeldef(covariables, scale.form)
- # shape.model <- modeldef(covariables, shape.form)
- # loc.dsgn.mat <- loc.model$dsgn.mat
- # scale.dsgn.mat <- scale.model$dsgn.mat
- # shape.dsgn.mat <- shape.model$dsgn.mat
- # loc.pen.mat <- loc.model$pen.mat
- # scale.pen.mat <- scale.model$pen.mat
- # shape.pen.mat <- shape.model$pen.mat
- # loc.penalty <- loc.model$penalty.tot
- # scale.penalty <- scale.model$penalty.tot
- # shape.penalty <- shape.model$penalty.tot
- # n.loccoeff <- ncol(loc.dsgn.mat)
- # n.scalecoeff <- ncol(scale.dsgn.mat)
- # n.shapecoeff <- ncol(shape.dsgn.mat)
- # n.pparloc <- loc.model$n.ppar
- # n.pparscale <- scale.model$n.ppar
- # n.pparshape <- shape.model$n.ppar
- # loc.names <- paste("locCoeff", 1:n.loccoeff, sep = "")
- # scale.names <- paste("scaleCoeff", 1:n.scalecoeff, sep = "")
- # shape.names <- paste("shapeCoeff", 1:n.shapecoeff, sep = "")
- # if (use.temp.cov[1]) {
- # temp.model.loc <- modeldef(temp.cov, temp.form.loc)
- # temp.dsgn.mat.loc <- temp.model.loc$dsgn.mat
- # temp.pen.mat.loc <- temp.model.loc$pen.mat
- # temp.penalty.loc <- temp.model.loc$penalty.tot
- # n.tempcoeff.loc <- ncol(temp.dsgn.mat.loc)
- # n.ppartemp.loc <- temp.model.loc$n.ppar
- # temp.names.loc <- paste("tempCoeffLoc", 1:n.tempcoeff.loc,
- # sep = "")
- # }
- # else {
- # temp.model.loc <- temp.dsgn.mat.loc <- temp.pen.mat.loc <- temp.names.loc <- NULL
- # n.tempcoeff.loc <- n.ppartemp.loc <- temp.penalty.loc <- 0
- # }
- # if (use.temp.cov[2]) {
- # temp.model.scale <- modeldef(temp.cov, temp.form.scale)
- # temp.dsgn.mat.scale <- temp.model.scale$dsgn.mat
- # temp.pen.mat.scale <- temp.model.scale$pen.mat
- # temp.penalty.scale <- temp.model.scale$penalty.tot
- # n.tempcoeff.scale <- ncol(temp.dsgn.mat.scale)
- # n.ppartemp.scale <- temp.model.scale$n.ppar
- # temp.names.scale <- paste("tempCoeffScale", 1:n.tempcoeff.scale,
- # sep = "")
- # }
- # else {
- # temp.model.scale <- temp.dsgn.mat.scale <- temp.pen.mat.scale <- temp.names.scale <- NULL
- # n.tempcoeff.scale <- n.ppartemp.scale <- temp.penalty.scale <- 0
- # }
- # if (use.temp.cov[3]) {
- # temp.model.shape <- modeldef(temp.cov, temp.form.shape)
- # temp.dsgn.mat.shape <- temp.model.shape$dsgn.mat
- # temp.pen.mat.shape <- temp.model.shape$pen.mat
- # temp.penalty.shape <- temp.model.shape$penalty.tot
- # n.tempcoeff.shape <- ncol(temp.dsgn.mat.shape)
- # n.ppartemp.shape <- temp.model.shape$n.ppar
- # temp.names.shape <- paste("tempCoeffShape", 1:n.tempcoeff.shape,
- # sep = "")
- # }
- # else {
- # temp.model.shape <- temp.dsgn.mat.shape <- temp.pen.mat.shape <- temp.names.shape <- NULL
- # n.tempcoeff.shape <- n.ppartemp.shape <- temp.penalty.shape <- 0
- # }
- # param <- c(loc.names, scale.names, shape.names, temp.names.loc,
- # temp.names.scale, temp.names.shape)
- # nllik <- function(x) x
- # body(nllik) <- parse(text = paste("-.C(C_spatgevlik, as.double(data), as.double(covariables),\n as.integer(n.site), as.integer(n.obs), as.double(loc.dsgn.mat), as.double(loc.pen.mat),\n as.integer(n.loccoeff), as.integer(n.pparloc), as.double(loc.penalty),\n as.double(scale.dsgn.mat), as.double(scale.pen.mat), as.integer(n.scalecoeff),\n as.integer(n.pparscale), as.double(scale.penalty), as.double(shape.dsgn.mat),\n as.double(shape.pen.mat), as.integer(n.shapecoeff), as.integer(n.pparshape),\n as.double(shape.penalty), as.integer(use.temp.cov), as.double(temp.dsgn.mat.loc),\n as.double(temp.pen.mat.loc), as.integer(n.tempcoeff.loc), as.integer(n.ppartemp.loc),\n as.double(temp.penalty.loc), as.double(temp.dsgn.mat.scale), as.double(temp.pen.mat.scale),\n as.integer(n.tempcoeff.scale), as.integer(n.ppartemp.scale), as.double(temp.penalty.scale),\n as.double(temp.dsgn.mat.shape), as.double(temp.pen.mat.shape), as.integer(n.tempcoeff.shape),\n as.integer(n.ppartemp.shape), as.double(temp.penalty.shape),",
- # paste("as.double(c(", paste(loc.names, collapse = ","),
- # ")), "), paste("as.double(c(", paste(scale.names,
- # collapse = ","), ")), "), paste("as.double(c(", paste(shape.names,
- # collapse = ","), ")), "), paste("as.double(c(", paste(temp.names.loc,
- # collapse = ","), ")), "), paste("as.double(c(", paste(temp.names.scale,
- # collapse = ","), ")), "), paste("as.double(c(", paste(temp.names.shape,
- # collapse = ","), ")), "), "dns = double(1), NAOK = TRUE)$dns"))
- # form.nllik <- NULL
- # for (i in 1:length(param)) form.nllik <- c(form.nllik, alist(a = ))
- # names(form.nllik) <- param
- # formals(nllik) <- form.nllik
- # if (missing(start)) {
- # loc <- scale <- shape <- rep(0, n.site)
- # for (i in 1:n.site) {
- # gev.param <- gevmle(data[, i])
- # loc[i] <- gev.param["loc"]
- # scale[i] <- gev.param["scale"]
- # shape[i] <- gev.param["shape"]
- # }
- # locCoeff <- loc.model$init.fun(loc)
- # scaleCoeff <- scale.model$init.fun(scale)
- # shapeCoeff <- shape.model$init.fun(shape)
- # locCoeff[is.na(locCoeff)] <- 0
- # scaleCoeff[is.na(scaleCoeff)] <- 0
- # shapeCoeff[is.na(shapeCoeff)] <- 0
- # scales.hat <- scale.model$dsgn.mat %*% scaleCoeff
- # if (any(scales.hat <= 0))
- # scaleCoeff[1] <- scaleCoeff[1] - 1.001 * min(scales.hat)
- # names(locCoeff) <- loc.names
- # names(scaleCoeff) <- scale.names
- # names(shapeCoeff) <- shape.names
- # if (use.temp.cov[1]) {
- # tempCoeff.loc <- rep(0, n.tempcoeff.loc)
- # names(tempCoeff.loc) <- temp.names.loc
- # }
- # else tempCoeff.loc <- NULL
- # if (use.temp.cov[2]) {
- # tempCoeff.scale <- rep(0, n.tempcoeff.scale)
- # names(tempCoeff.scale) <- temp.names.scale
- # }
- # else tempCoeff.scale <- NULL
- # if (use.temp.cov[3]) {
- # tempCoeff.shape <- rep(0, n.tempcoeff.shape)
- # names(tempCoeff.shape) <- temp.names.shape
- # }
- # else tempCoeff.shape <- NULL
- # start <- as.list(c(locCoeff, scaleCoeff, shapeCoeff,
- # tempCoeff.loc, tempCoeff.scale, tempCoeff.shape))
- # # start <- start[!(param %in% names(list(...)))]
- # }
- # if (!length(start))
- # stop("there are no parameters left to maximize over")
- # nm <- names(start)
- # l <- length(nm)
- # f <- formals(nllik)
- # names(f) <- param
- # m <- match(nm, param)
- # if (any(is.na(m)))
- # stop("'start' specifies unknown arguments")
- # formals(nllik) <- c(f[m], f[-m])
- # nllh <- function(p, ...) nllik(p, ...)
- # stop(param)
- # if (l > 1)
- # body(nllh) <- parse(text = paste("nllik(", paste("p[",
- # 1:l, "]", collapse = ", "), ", ...)"))
- # fixed.param <- list(...)[names(list(...)) %in% param]
- # if (any(!(param %in% c(nm, names(fixed.param)))))
- # stop("unspecified parameters")
- # start.arg <- c(list(p = unlist(start)), fixed.param)
print(res['fitted.values'])
}
+# fitspatgev_test with 3D data
+fitspatgev_3D_test <- function (n.obs){
+ n.site = 2
+ covariables <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2)
+ colnames(covariables) = c("E", "N")
+
+ # Generate the data
+ data <- rmaxstab(n.obs, covariables, "gauss", cov11 = 100, cov12 = 25, cov22 = 220)
+ print(data)
+
+ loc.form = loc ~ N
+ scale.form = scale ~ N
+ shape.form = shape ~ N
+
+ # Add the temporal covariates
+ temp.form.loc = loc ~ T
+ temp.form.scale = scale ~ T
+ temp.form.shape = shape ~ T
+ # temp.form.shape = NULL
+
+ temp.cov = matrix(1:n.obs, ncol=1)
+ colnames(temp.cov) = c("T")
+
+ # namedlist = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0,
+ # tempCoeffLoc1=0.0, tempCoeffLoc2=1.0, tempCoeffScale1=0.1, tempCoeffShape1=1.0)
+ start = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0,
+ tempCoeffLoc1=0.0, tempCoeffLoc2=1.0, tempCoeffScale1=0.1, tempCoeffShape1=1.0)
+
+
+
+ start = list(locCoeff1=0.0, locCoeff2=1.0, scaleCoeff1=0.1, shapeCoeff1=1.0, scaleCoeff2=0.1, shapeCoeff2=1.0,
+ tempCoeffLoc1=0.0, tempCoeffScale1=0.1, tempCoeffShape1=0.1)
+ res = fitspatgev(data=data, covariables=covariables, start=start, loc.form=loc.form, scale.form=scale.form,shape.form=shape.form,
+ temp.cov=temp.cov, temp.form.loc=temp.form.loc, temp.form.scale=temp.form.scale, temp.form.shape = temp.form.shape)
+
+ n.site <- ncol(data)
+ n.obs <- nrow(data)
+ if (n.site != nrow(covariables))
+ stop("'data' and 'covariates' doesn't match")
+ use.temp.cov <- c(!is.null(temp.form.loc), !is.null(temp.form.scale),
+ !is.null(temp.form.shape))
+ if (any(use.temp.cov) && (n.obs != nrow(temp.cov)))
+ stop("'data' and 'temp.cov' doesn't match")
+ if (any(use.temp.cov) && is.null(temp.cov))
+ stop("'temp.cov' must be supplied if at least one temporal formula is given")
+ loc.form <- update(loc.form, y ~ .)
+ scale.form <- update(scale.form, y ~ .)
+ shape.form <- update(shape.form, y ~ .)
+ print(use.temp.cov)
+ if (use.temp.cov[1])
+ temp.form.loc <- update(temp.form.loc, y ~ . + 0)
+ if (use.temp.cov[2])
+ temp.form.scale <- update(temp.form.scale, y ~ . + 0)
+ if (use.temp.cov[3])
+ temp.form.shape <- update(temp.form.shape, y ~ . + 0)
+ loc.model <- modeldef(covariables, loc.form)
+ scale.model <- modeldef(covariables, scale.form)
+ shape.model <- modeldef(covariables, shape.form)
+ loc.dsgn.mat <- loc.model$dsgn.mat
+ scale.dsgn.mat <- scale.model$dsgn.mat
+ shape.dsgn.mat <- shape.model$dsgn.mat
+ loc.pen.mat <- loc.model$pen.mat
+ scale.pen.mat <- scale.model$pen.mat
+ shape.pen.mat <- shape.model$pen.mat
+ loc.penalty <- loc.model$penalty.tot
+ scale.penalty <- scale.model$penalty.tot
+ shape.penalty <- shape.model$penalty.tot
+ n.loccoeff <- ncol(loc.dsgn.mat)
+ n.scalecoeff <- ncol(scale.dsgn.mat)
+ n.shapecoeff <- ncol(shape.dsgn.mat)
+ n.pparloc <- loc.model$n.ppar
+ n.pparscale <- scale.model$n.ppar
+ n.pparshape <- shape.model$n.ppar
+ loc.names <- paste("locCoeff", 1:n.loccoeff, sep = "")
+ scale.names <- paste("scaleCoeff", 1:n.scalecoeff, sep = "")
+ shape.names <- paste("shapeCoeff", 1:n.shapecoeff, sep = "")
+ if (use.temp.cov[1]) {
+ temp.model.loc <- modeldef(temp.cov, temp.form.loc)
+ temp.dsgn.mat.loc <- temp.model.loc$dsgn.mat
+ temp.pen.mat.loc <- temp.model.loc$pen.mat
+ temp.penalty.loc <- temp.model.loc$penalty.tot
+ n.tempcoeff.loc <- ncol(temp.dsgn.mat.loc)
+ n.ppartemp.loc <- temp.model.loc$n.ppar
+ temp.names.loc <- paste("tempCoeffLoc", 1:n.tempcoeff.loc,
+ sep = "")
+ }
+ else {
+ temp.model.loc <- temp.dsgn.mat.loc <- temp.pen.mat.loc <- temp.names.loc <- NULL
+ n.tempcoeff.loc <- n.ppartemp.loc <- temp.penalty.loc <- 0
+ }
+ if (use.temp.cov[2]) {
+ temp.model.scale <- modeldef(temp.cov, temp.form.scale)
+ temp.dsgn.mat.scale <- temp.model.scale$dsgn.mat
+ temp.pen.mat.scale <- temp.model.scale$pen.mat
+ temp.penalty.scale <- temp.model.scale$penalty.tot
+ n.tempcoeff.scale <- ncol(temp.dsgn.mat.scale)
+ n.ppartemp.scale <- temp.model.scale$n.ppar
+ temp.names.scale <- paste("tempCoeffScale", 1:n.tempcoeff.scale,
+ sep = "")
+ }
+ else {
+ temp.model.scale <- temp.dsgn.mat.scale <- temp.pen.mat.scale <- temp.names.scale <- NULL
+ n.tempcoeff.scale <- n.ppartemp.scale <- temp.penalty.scale <- 0
+ }
+ if (use.temp.cov[3]) {
+ temp.model.shape <- modeldef(temp.cov, temp.form.shape)
+ temp.dsgn.mat.shape <- temp.model.shape$dsgn.mat
+ temp.pen.mat.shape <- temp.model.shape$pen.mat
+ temp.penalty.shape <- temp.model.shape$penalty.tot
+ n.tempcoeff.shape <- ncol(temp.dsgn.mat.shape)
+ n.ppartemp.shape <- temp.model.shape$n.ppar
+ temp.names.shape <- paste("tempCoeffShape", 1:n.tempcoeff.shape,
+ sep = "")
+ }
+ else {
+ temp.model.shape <- temp.dsgn.mat.shape <- temp.pen.mat.shape <- temp.names.shape <- NULL
+ n.tempcoeff.shape <- n.ppartemp.shape <- temp.penalty.shape <- 0
+ }
+ param <- c(loc.names, scale.names, shape.names, temp.names.loc,
+ temp.names.scale, temp.names.shape)
+ print(param)
+}
+
# Boolean for python call
diff --git a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
index c9c988823a86e4644e3a05fa5184ca24df1d9f05..4196d16e51725d6a2c29e74628607772d07c69e4 100644
--- a/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
+++ b/extreme_estimator/extreme_models/margin_model/margin_function/linear_margin_function.py
@@ -95,6 +95,8 @@ class LinearMarginFunction(IndependentMarginFunction):
temporal_names = [name for name in self.coordinates.coordinates_temporal_names
if self.coefficient_name_to_dim(self.coordinates)[name] in linear_dims]
temporal_form = self.gev_param_name_to_linear_coef[gev_param_name].temporal_form_dict(temporal_names)
+ # Specifying a formula '~ 1' creates a bug in fitspatgev of SpatialExtreme R package
+ assert not any(['1' in formula for formula in temporal_form.values()])
form_dict.update(temporal_form)
return form_dict
diff --git a/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py b/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py
index 66e4067183446d086362455f4605eddcd8ffd786..8954e5e0833e3d3fc4145627bd89f7c78ea71213 100644
--- a/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py
+++ b/extreme_estimator/extreme_models/margin_model/param_function/linear_coef.py
@@ -32,7 +32,7 @@ class LinearCoef(object):
def coef_template_str(cls, gev_param_name: str, coefficient_name: str) -> str:
"""
Example of coef that can be specified
- -for the spatial covariates: locCoeff1
+ -for the spatial covariates: locCoeff
-for the temporal covariates: tempCoeffLoc
:param gev_param_name:
:param coefficient_name:
@@ -44,24 +44,41 @@ class LinearCoef(object):
else:
return 'tempCoeff' + gev_param_name.title() + '{}'
+ @staticmethod
+ def has_dependence_in_spatial_coordinates(dim_to_coefficient_name):
+ return any([coefficient_name in AbstractCoordinates.COORDINATE_SPATIAL_NAMES
+ for coefficient_name in dim_to_coefficient_name.values()])
+
+ @classmethod
+ def add_intercept_dim(cls, dims):
+ return [0] + dims
+
@classmethod
def from_coef_dict(cls, coef_dict: Dict[str, float], gev_param_name: str, linear_dims: List[int],
- dim_to_coefficient_name):
- dims = [0] + linear_dims
+ dim_to_coefficient_name: Dict[int, str]):
+ dims = cls.add_intercept_dim(linear_dims)
dim_to_coef = {}
- for j, dim in enumerate(dims, 1):
+ j = 1
+ for dim in dims:
coefficient_name = dim_to_coefficient_name[dim]
+ if coefficient_name == AbstractCoordinates.COORDINATE_T:
+ j = 1
coef = coef_dict[cls.coef_template_str(gev_param_name, coefficient_name).format(j)]
dim_to_coef[dim] = coef
+ j += 1
return cls(gev_param_name, dim_to_coef)
- def coef_dict(self, linear_dims, dim_to_coefficient_name) -> Dict[str, float]:
- dims = [0] + linear_dims
+ def coef_dict(self, linear_dims, dim_to_coefficient_name: Dict[int, str]) -> Dict[str, float]:
+ dims = self.add_intercept_dim(linear_dims)
coef_dict = {}
- for j, dim in enumerate(dims, 1):
+ j = 1
+ for dim in dims:
coefficient_name = dim_to_coefficient_name[dim]
+ if coefficient_name == AbstractCoordinates.COORDINATE_T:
+ j = 1
coef = self.dim_to_coef[dim]
coef_dict[self.coef_template_str(self.gev_param_name, coefficient_name).format(j)] = coef
+ j += 1
return coef_dict
def form_dict(self, names: List[str]) -> Dict[str, str]:
@@ -71,9 +88,10 @@ class LinearCoef(object):
def spatial_form_dict(self, coordinate_spatial_names: List[str]) -> Dict[str, str]:
"""
Example of formula that could be specified:
+ loc.form = loc ~ 1
loc.form = loc ~ coord_x
scale.form = scale ~ coord_y
- shape.form = shape ~ coord_x+coord_y
+ shape.form = shape ~ coord_x + coord_y
:return:
"""
assert all([name in AbstractCoordinates.COORDINATE_SPATIAL_NAMES for name in coordinate_spatial_names])
@@ -82,10 +100,13 @@ class LinearCoef(object):
def temporal_form_dict(self, coordinate_temporal_names: List[str]) -> Dict[str, str]:
"""
Example of formula that could be specified:
- temp.loc.form = loc ~ coord_t
- temp.scale.form = scale ~ coord_t
- temp.shape.form = shape ~ 1
+ temp.form.loc = loc ~ coord_t
+ Example of formula that could not be specified
+ temp.loc.form = shape ~ 1
:return:
"""
assert all([name in [AbstractCoordinates.COORDINATE_T] for name in coordinate_temporal_names])
- return {'temp.' + k: v for k, v in self.form_dict(coordinate_temporal_names).items()}
+ k, v = self.form_dict(coordinate_temporal_names).popitem()
+ k = 'temp.form.' + k.split('.')[0]
+ v = 'NULL' if '1' in v else v
+ return {k: v}
diff --git a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
index 57d4254ba138a9fa533758f1772927f6581c6bf5..5622ad5f3ce0c44192c0564d5e0d7d3ab5c32834 100644
--- a/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
+++ b/extreme_estimator/extreme_models/margin_model/smooth_margin_model.py
@@ -58,13 +58,25 @@ class LinearMarginModel(AbstractMarginModel):
params[(gev_param_name, dim)] = coef
return cls(coordinates, params_sample=params, params_start_fit=params)
- def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray,
- df_coordinates: pd.DataFrame) -> Dict[str, float]:
+ def fitmargin_from_maxima_gev(self, maxima_gev: np.ndarray, df_coordinates_spatial: pd.DataFrame,
+ df_coordinates_temporal: pd.DataFrame) -> Dict[str, float]:
+ # The reshaping on the line below is only valid if we have a single observation per spatio-temporal point
+ if maxima_gev.shape[1] == 1:
+ maxima_gev = maxima_gev.reshape([len(df_coordinates_temporal), len(df_coordinates_spatial)])
data = np.transpose(maxima_gev)
- covariables = get_coord(df_coordinates)
+
fit_params = get_margin_formula(self.margin_function_start_fit.form_dict)
- fit_params['start'] = r.list(**self.margin_function_start_fit.coef_dict)
- res = safe_run_r_estimator(function=r.fitspatgev, use_start=self.use_start_value, data=data, covariables=covariables, **fit_params)
+
+ # Covariables
+ covariables = get_coord(df_coordinates=df_coordinates_spatial)
+ fit_params['temp.cov'] = get_coord(df_coordinates=df_coordinates_temporal)
+
+ # Start parameters
+ coef_dict = self.margin_function_start_fit.coef_dict
+ fit_params['start'] = r.list(**coef_dict)
+
+ res = safe_run_r_estimator(function=r.fitspatgev, use_start=self.use_start_value, data=data,
+ covariables=covariables, **fit_params)
return retrieve_fitted_values(res)
diff --git a/extreme_estimator/extreme_models/utils.py b/extreme_estimator/extreme_models/utils.py
index 97c0c7eeb026426d15719e9d19808dba2e9dd63d..ec8defa76c3be4fc683cf7f5f64008af6eb3d1b6 100644
--- a/extreme_estimator/extreme_models/utils.py
+++ b/extreme_estimator/extreme_models/utils.py
@@ -16,7 +16,11 @@ r = ro.R()
numpy2ri.activate()
pandas2ri.activate()
r.library('SpatialExtremes')
-# todo: R is not reloading all the time, the SpatialExtremes, so it's quite hard to debug or print in the code...
+
+
+# Notice: R is not reloading all the time, the SpatialExtremes, so it's quite hard to debug or print in the code...
+# the best solution for debugging is to copy/paste the code module into a file that belongs to me, and then
+# I can put print & stop in the code, and I can understand where are the problems
def set_seed_r(seed=42):
r("set.seed({})".format(seed))
diff --git a/spatio_temporal_dataset/coordinates/abstract_coordinates.py b/spatio_temporal_dataset/coordinates/abstract_coordinates.py
index 0c35f5c2b1ecfc4087c163cb9970324e413621bb..5adcffcffe8bb37d249e78cfe37172d7cb89ad70 100644
--- a/spatio_temporal_dataset/coordinates/abstract_coordinates.py
+++ b/spatio_temporal_dataset/coordinates/abstract_coordinates.py
@@ -148,6 +148,8 @@ class AbstractCoordinates(object):
def nb_coordinates(self) -> int:
return len(self.coordinates_names)
+ # Spatial attributes
+
@property
def coordinates_spatial_names(self) -> List[str]:
return [name for name in self.COORDINATE_SPATIAL_NAMES if name in self.df_all_coordinates.columns]
@@ -156,6 +158,14 @@ class AbstractCoordinates(object):
def nb_coordinates_spatial(self) -> int:
return len(self.coordinates_spatial_names)
+ def df_spatial_coordinates(self, split: Split = Split.all) -> pd.DataFrame:
+ if self.nb_coordinates_spatial == 0:
+ return pd.DataFrame()
+ else:
+ return self.df_coordinates(split).loc[:, self.coordinates_spatial_names].drop_duplicates()
+
+ # Temporal attributes
+
@property
def coordinates_temporal_names(self) -> List[str]:
return [self.COORDINATE_T] if self.COORDINATE_T in self.df_all_coordinates else []
@@ -164,6 +174,12 @@ class AbstractCoordinates(object):
def nb_coordinates_temporal(self) -> int:
return len(self.coordinates_temporal_names)
+ def df_temporal_coordinates(self, split: Split = Split.all) -> pd.DataFrame:
+ if self.nb_coordinates_temporal == 0:
+ return pd.DataFrame()
+ else:
+ return self.df_coordinates(split).loc[:, self.coordinates_temporal_names].drop_duplicates()
+
# Visualization
@property
diff --git a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
index e3e1593acbb64cfcc8ef1f88681b1697b4251ab4..97c4a63ce9c89c4d6dc5dd96fc8ed04747e6a8b3 100644
--- a/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
+++ b/test/test_extreme_estimator/test_estimator/test_margin_estimators.py
@@ -2,19 +2,25 @@ import unittest
from extreme_estimator.estimator.margin_estimator.abstract_margin_estimator import SmoothMarginEstimator
from spatio_temporal_dataset.dataset.simulation_dataset import MarginDataset
-from test.test_utils import load_smooth_margin_models, load_test_1D_and_2D_spatial_coordinates
+from test.test_utils import load_smooth_margin_models, load_test_1D_and_2D_spatial_coordinates, \
+ load_test_spatiotemporal_coordinates
class TestSmoothMarginEstimator(unittest.TestCase):
DISPLAY = False
- nb_points = 2
- nb_obs = 2
- def setUp(self):
- super().setUp()
+ def test_smooth_margin_estimator_spatial(self):
+ self.nb_points = 2
+ self.nb_obs = 2
self.coordinates = load_test_1D_and_2D_spatial_coordinates(nb_points=self.nb_points)
- def test_smooth_margin_estimator(self):
+ def test_smooth_margin_estimator_spatio_temporal(self):
+ self.nb_points = 2
+ self.nb_steps = 2
+ self.nb_obs = 1
+ self.coordinates = load_test_spatiotemporal_coordinates(nb_steps=self.nb_steps, nb_points=self.nb_points)
+
+ def tearDown(self) -> None:
for coordinates in self.coordinates:
smooth_margin_models = load_smooth_margin_models(coordinates=coordinates)
for margin_model in smooth_margin_models[1:]:
@@ -33,3 +39,6 @@ class TestSmoothMarginEstimator(unittest.TestCase):
if __name__ == '__main__':
unittest.main()
+ # t = TestSmoothMarginEstimator()
+ # t.test_smooth_margin_estimator_spatio_temporal()
+ # t.tearDown()