Newer
Older
#'@title Computing two dinamical proxies of the attractor in s2dv_cube.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function computes two dinamical proxies of the attractor:
#'The local dimension (d) and the inverse of the persistence (theta) for an
#'These two parameters will be used as a condition for the computation of dynamical
#'scores to measure predictability and to compute bias correction conditioned by
#'the dynamics with the function DynBiasCorrection
#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019).
#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability.
#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes.
#' Scientific Reports, 7-41278, 2017.
#'
#'@param data a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow
#'and the grids in ncol(dat(time,grids)
#
#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta
#'
#'@param ncores The number of cores to use in parallel computation
#'
#'@return dim and theta
#'
#'@examples
#'# Example 1: Computing the attractor using simple s2dv data
#'attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6)
#'@export
CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){
if (!inherits(data, 's2dv_cube')) {
stop("Parameter 'data' must be of the class 's2dv_cube', ",
"as output by CSTools::CST_Load.")
}
if (is.null(quanti)) {
data$data <- ProxiesAttractor(data = data$data, quanti = quanti, ncores = ncores)
return(data)
}
#'@rdname ProxiesAttractor
#'
#'@title Computing two dinamical proxies of the attractor.
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function computes two dinamical proxies of the attractor:
#'The local dimension (d) and the inverse of the persistence (theta).
#'These two parameters will be used as a condition for the computation of dynamical
#'scores to measure predictability and to compute bias correction conditioned by
#'the dynamics with the function DynBiasCorrection.
#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in:
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019).
#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability.
#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes.
#' Scientific Reports, 7-41278, 2017.
#'
#'@param data a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'.
#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta
#'@param ncores The number of cores to use in parallel computation.
#'
#'@return dim and theta
#'@import multiApply
#'
#'@examples
#'# Example 1: Computing the attractor using simple data
#'# Creating an example of matrix data(time,grids):
#'mat <- array(rnorm(36 * 40), c(time = 36, grid = 40))
#'qm <- 0.90 # imposing a threshold
#'Attractor <- ProxiesAttractor(data = mat, quanti = qm)
#'# to plot the result
#'time = c(1:length(Attractor$theta))
#'layout(matrix(c(1, 3, 2, 3), 2, 2))
#'plot(time, Attractor$dim, xlab = 'time', ylab = 'd',
#' main = 'local dimension', type = 'l')
#'plot(time, Attractor$theta, xlab = 'time', ylab = 'theta', main = 'theta')
#'plot(Attractor$dim, Attractor$theta, col = 'blue',
#' main = "Proxies of the Attractor",
#' xlab = "local dimension", ylab = "theta", lwd = 8, 'p')
ProxiesAttractor <- function(data, quanti, ncores = NULL){
if (is.null(data)) {
stop("Parameter 'data' cannot be NULL.")
}
if (is.null(quanti)) {
stop("Parameter 'quanti' is mandatory")
}
if (any(names(dim(data)) %in% 'sdate')) {
if (any(names(dim(data)) %in% 'ftime')) {
data <- MergeDims(data, merge_dims = c('ftime', 'sdate'),
rename_dim = 'time')
}
}
if (!(any(names(dim(data)) %in% 'time'))){
stop("Parameter 'data' must have a temporal dimension named 'time'.")
if (any(names(dim(data)) %in% 'lat')) {
if (any(names(dim(data)) %in% 'lon')) {
data <- MergeDims(data, merge_dims = c('lon', 'lat'),
rename_dim = 'grid')
}
}
if (any(names(dim(data)) %in% 'latitude')) {
if (any(names(dim(data)) %in% 'longitude')) {
data <- MergeDims(data, merge_dims = c('longitude', 'latitude'),
rename_dim = 'grid')
}
}
if(!(any(names(dim(data)) %in% 'grid'))){
stop("Parameter 'data' must have a spatial dimension named 'grid'.")
}
attractor <- Apply(data, target_dims = c('time', 'grid'),
fun = .proxiesattractor,
quanti = quanti , ncores = ncores)
attractor <- lapply(attractor,
FUN = function(x, dimname){
names(dim(x))[dimname] <- 'time'
return(x)},
dimname = which(names(dim(attractor[[1]])) == 'dim2'))
return(list(dim = attractor$dim, theta = attractor$theta))
}
.proxiesattractor <- function(data, quanti) {
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
# expected dimensions data: time and grid
logdista <- Apply(data, target_dims = 'grid',
fun = function(x, y){
-log(colMeans((y - as.vector(x))^2))},
y = t(data))[[1]]
#Computation of theta
Theta <- function(logdista, quanti){
#Compute the thheshold corresponding to the quantile
thresh <- quantile(logdista, quanti, na.rm = TRUE)
logdista[which(logdista == 'Inf')] <- NaN
Li <- which(as.vector(logdista) > as.numeric(thresh))
#Length of each cluster
Ti <- diff(Li)
N <- length(Ti)
q <- 1 - quanti
Si <- Ti - 1
Nc <- length(which(Si > 0))
N <- length(Ti)
theta <- (sum(q * Si) + N + Nc - sqrt(((sum(q * Si) + N + Nc)^2) -
8 * Nc * sum(q * Si))) / (2 * sum(q * Si))
#Sort the exceedances
logdista <- sort(logdista)
#Find all the Peaks over Thresholds.
findidx <- which(as.vector(logdista) > as.numeric(thresh))
if(length(findidx) < 1) {
stop("Parameter 'quanti' is too high for the length of the data provided.")
}
logextr <- logdista[findidx[[1]]:(length(logdista) - 1)]
#The inverse of the dimension is just the average of the exceedances
dim <- 1 /mean(as.numeric(logextr) - as.numeric(thresh))
return(list(dim = dim, theta = theta))
}
names(dim(logdista)) <- c('dim1', 'dim2')
proxies <- Apply(data = list(logdista = logdista),
return(list(dim = proxies$dim, theta = proxies$theta))