diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index e99677176490c53be8b881e190f1168ab532f39e..554121de2682e8f76457a453b6ef0858520adb1a 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -21,11 +21,15 @@ #'Dynamical proxies of North Atlantic predictability and extremes. #'Scientific Reports, 7-41278, 2017. #' -#'@param data An 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. +#'@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 spatial_dims could be 1 or 2 they would be merged +#'@param time_dims could be 1 or more they would be merged +#'@param na.rm logical indicating whether to remove the NAs when calculating the logaritmic distance (default TRUE). +#'@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 @@ -39,8 +43,9 @@ #'attractor <- CST_ProxiesAttractor(data = data, quanti = 0.6) #'@import multiApply #'@export -CST_ProxiesAttractor <- function(data, quanti, ncores = NULL) { - # Check 's2dv_cube' +CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), + time_dims = c('ftime', 'sdate'), + na.rm = TRUE, ncores = NULL){ if (!inherits(data, 's2dv_cube')) { stop("Parameter 'data' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") @@ -51,7 +56,9 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL) { } data$data <- ProxiesAttractor(data = data$data, quanti = quanti, - ncores = ncores) + spatial_dims = spatial_dims, + time_dims = time_dims, + na.rm = na.rm, ncores = ncores) return(data) } @@ -77,11 +84,11 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL) { #' 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 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 spatial_dims could be 1 or 2. If more than 1, they will be merged. +#'@param time_dims could be 1 or more. If more than 1, they will be merged +#'@param na.rm logical indicating whether to remove the NAs when calculating the logaritmic distance. #'@param ncores The number of cores to use in parallel computation. #' #'@return dim and theta @@ -98,54 +105,72 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL) { #' main = 'local dimension', type = 'l') #'@import multiApply #'@export -ProxiesAttractor <- function(data, quanti, ncores = NULL){ + +ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), + time_dims = c('ftime', 'sdate'), + na.rm = TRUE, 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 (is.null(spatial_dims)) { + data <- InsertDim(data, posdim = 1, lendim = 1, name = 'grid') + spatial_dims <- 'grid' + } + if (length(spatial_dims) == 2) { + data <- MergeDims(data, merge_dims = spatial_dims, rename = 'grid') + spatial_dims <- 'grid' + } else if (length(spatial_dims) > 2) { + stop("Paramenter 'spatial_dims' maximum length 2.") } - 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 (length(time_dims) == 2) { + pos1 <- which(names(dim(data)) == time_dims[1]) + pos2 <- which(names(dim(data)) == time_dims[2]) + #split_time <- dim(data)[time_dims] + if (pos1 < pos2) { + time_dims <- time_dims[c(1, 2)] + split_time <- dim(data)[c(pos1, pos2)] + } else { + time_dims <- time_dims[c(2, 1)] + split_time <- dim(data)[c(pos2, pos1)] + } + data <- MergeDims(data, merge_dims = time_dims, rename = 'temporals') + time_dims <- 'temporals' + } else if (length(time_dims) > 2) { + stop("Parameter 'time_dims' maximum length 2.") } - if(!(any(names(dim(data)) %in% 'grid'))){ - stop("Parameter 'data' must have a spatial dimension named 'grid'.") + if (any(is.na(data))) { + warning("Parameter 'data' contains NAs and they will be removed.") } - attractor <- Apply(data, target_dims = c('time', 'grid'), + attractor <- Apply(data, target_dims = c(time_dims, spatial_dims), fun = .proxiesattractor, - quanti = quanti , ncores = ncores) + quanti = quanti , na.rm = na.rm, ncores = ncores) # rename dimensions attractor <- lapply(attractor, FUN = function(x, dimname){ - names(dim(x))[dimname] <- 'time' + names(dim(x))[dimname] <- time_dims return(x)}, dimname = which(names(dim(attractor[[1]])) == 'dim2')) + if (length(split_time) == 2) { + attractor <- lapply(attractor, SplitDim, split_dim = time_dims, + indices = NULL, freq = split_time[1], + new_dim_name = names(split_time[2])) + attractor <- lapply(attractor, function(x, name) { + names(dim(x))[which(names(dim(x)) == time_dims)] <- + names(split_time)[1] + return(x)}) + } return(list(dim = attractor$dim, theta = attractor$theta)) } -.proxiesattractor <- function(data, quanti) { +.proxiesattractor <- function(data, quanti, na.rm) { # expected dimensions data: time and grid - logdista <- Apply(data, target_dims = 'grid', + grid_dim <- names(dim(data))[2] + logdista <- Apply(data, target_dims = grid_dim, fun = function(x, y){ - -log(colMeans((y - as.vector(x))^2))}, + -log(colMeans((y - as.vector(x))^2, na.rm = na.rm))}, y = t(data))[[1]] #Computation of theta diff --git a/R/Predictability.R b/R/Predictability.R index 680666df92b5c779e79c7488f5d6906672f88eb3..d93c244220be12fe50dd1bea2f2be9742b1325c1 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -58,8 +58,8 @@ #' attractor <- ProxiesAttractor(dat = m, quanti = 0.60) #' predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) #'@export -Predictability <- function(dim, theta, ncores = NULL) { - +#' +Predictability<- function(dim, theta, ncores = NULL) { if (any(names(dim(dim)) %in% 'sdate')) { if (any(names(dim(dim)) %in% 'ftime')) { dim <- MergeDims(dim, merge_dims = c('ftime', 'sdate'), @@ -79,29 +79,42 @@ Predictability <- function(dim, theta, ncores = NULL) { if (!(any(names(dim(theta)) %in% 'time'))){ stop("Parameter 'data' must have a temporal dimension named 'time'.") } - + # The number of elements in each tercile can vary and รง + # Apply requires to merge them being of the same length + # there will be NAs in the outputs + length_theta <- Apply(list(theta), target_dims = 'time', + function(x) { + qt1 <- quantile(x, 0.33, na.rm = TRUE) + qt3 <- quantile(x, 0.66, na.rm = TRUE) + th3 <- which(x >= qt3) + th1 <- which(x <= qt1) + th2 <- which(x > qt1 & x < qt3) + max(length(th1), length(th2), length(th3))}, + ncores = ncores)$output1 + length_theta <- max(length_theta) pred <- Apply(list(dim, theta), target_dims = 'time', - fun = .predictability, + fun = .predictability, out_len = length_theta, ncores = ncores) dim(pred$dyn_scores) <- dim(theta) - return(list(pred.dim = list(qdim = list(pred$qdim.d1,pred$qdim.d2,pred$qdim.d3), - pos.d = list(pred$pos.d1,pred$pos.d2,pred$pos.d3)), - pred.theta = list(qtheta = list(pred$qtheta.th1,pred$qtheta.th2,pred$qtheta.th3), - pos.t = list(pred$pos.th1,pred$pos.th2,pred$pos.th3)), - dyn_scores = pred$dyn_scores)) + return(list(pred.dim = list(qdim = list(d1 = pred$qdim.d1, + d2 = pred$qdim.d2, + d3 = pred$qdim.d3), + pos.d = list(d1 = pred$pos.d1, + d2 = pred$pos.d2, + d3 = pred$pos.d3)), + pred.theta = list(qtheta = list(th1 = pred$qtheta.th1, + th2 = pred$qtheta.th2, + th3 = pred$qtheta.th3), + pos.t = list(th1 = pred$pos.th1, + th2 = pred$pos.th2, + th3 = pred$pos.th3)), + dyn_scores = pred$dyn_scores)) } -.predictability <- function(dim, theta) { - if (is.null(dim)) { - stop("Parameter 'dim' cannot be NULL.") - } - if (is.null(theta)) { - stop("Parameter 'theta' cannot be NULL.") - } +.predictability <- function(dim, theta, out_len) { if (length(dim) != length(theta)) { stop("Parameters 'dim' and 'theta' must have the same length") } pos <- c(1:length(dim)) - # dim qd1 <- quantile(dim, 0.33, na.rm = TRUE) qd3 <- quantile(dim, 0.66, na.rm = TRUE) @@ -109,7 +122,7 @@ Predictability <- function(dim, theta, ncores = NULL) { d1 <- which(dim <= qd1) d2 <- which(dim > qd1 & dim < qd3) qdim <- list(d1 = dim[d1], d2 = dim[d2], d3 = dim[d3]) - pos.d <- list(d1=d1,d2=d2,d3=d3) + pos.d <- list(d1 = d1, d2 = d2, d3 = d3) #theta qt1 <- quantile(theta, 0.33, na.rm = TRUE) @@ -119,6 +132,7 @@ Predictability <- function(dim, theta, ncores = NULL) { th2 <- which(theta > qt1 & theta < qt3) qtheta <- list(th1 = theta[th1], th2 = theta[th2], th3 = theta[th3]) pos.t <- list(th1 = th1, th2 = th2, th3 = th3) + #scores position d1th1 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th1)] @@ -141,10 +155,14 @@ Predictability <- function(dim, theta, ncores = NULL) { dyn_scores[which(pos %in% d2th3)]<- 3/9 dyn_scores[which(pos %in% d3th2)]<- 2/9 dyn_scores[which(pos %in% d3th3)]<- 1/9 - + + th1 <- c(th1, rep(NA, out_len - length(th1))) + th2 <- c(th2, rep(NA, out_len - length(th2))) + th3 <- c(th3, rep(NA, out_len - length(th3))) + return(list(qdim.d1 = dim[d1], qdim.d2 = dim[d2], qdim.d3 = dim[d3], - pos.d1 = d1, pos.d2 = d2, pos.d3 =d3, - qtheta.th1 = theta[th1], qtheta.th2 = theta[th2], qtheta.th3 = theta[th3], + pos.d1 = d1, pos.d2 = d2, pos.d3 = d3, + qtheta.th1 = theta[th1], qtheta.th2 = theta[th2], qtheta.th3 = theta[th3], pos.th1 = th1, pos.th2 = th2, pos.th3 = th3, dyn_scores = dyn_scores)) } diff --git a/man/CST_ProxiesAttractor.Rd b/man/CST_ProxiesAttractor.Rd index 58f949af388ee6c572421514b885d1e487e43f46..bb5df1214cc57126b31722337d44dcaa52c17c45 100644 --- a/man/CST_ProxiesAttractor.Rd +++ b/man/CST_ProxiesAttractor.Rd @@ -4,7 +4,14 @@ \alias{CST_ProxiesAttractor} \title{Computing two dinamical proxies of the attractor in s2dv_cube.} \usage{ -CST_ProxiesAttractor(data, quanti, ncores = NULL) +CST_ProxiesAttractor( + data, + quanti, + spatial_dims = c("lat", "lon"), + time_dims = c("ftime", "sdate"), + na.rm = TRUE, + ncores = NULL +) } \arguments{ \item{data}{An s2dv_cube object with the data to create the attractor. Must be @@ -13,7 +20,17 @@ a matrix with the timesteps in nrow and the grids in ncol(dat(time,grids)} \item{quanti}{A number lower than 1 indicating the quantile to perform the computation of local dimension and theta.} +<<<<<<< HEAD +\item{spatial_dims}{could be 1 or 2 they would be merged} + +\item{time_dims}{could be 1 or more they would be merged} + +\item{na.rm}{logical indicating whether to remove the NAs when calculating the logaritmic distance (default TRUE).} + +\item{ncores}{The number of cores to use in parallel computation} +======= \item{ncores}{The number of cores to use in parallel computation.} +>>>>>>> bd3fdc56cf575ad2f16f8d7b6024b618511de2ba } \value{ dim and theta diff --git a/man/ProxiesAttractor.Rd b/man/ProxiesAttractor.Rd index ffa1b36b18654fd83119cbfec3079abd6f0dde7d..5d728ec3926326acf3c0217d1d5b5150a8fb63d1 100644 --- a/man/ProxiesAttractor.Rd +++ b/man/ProxiesAttractor.Rd @@ -4,7 +4,14 @@ \alias{ProxiesAttractor} \title{Computing two dinamical proxies of the attractor.} \usage{ -ProxiesAttractor(data, quanti, ncores = NULL) +ProxiesAttractor( + data, + quanti, + spatial_dims = c("lon", "lat"), + time_dims = c("ftime", "sdate"), + na.rm = TRUE, + ncores = NULL +) } \arguments{ \item{data}{A multidimensional array with named dimensions to create the @@ -14,6 +21,12 @@ dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'.} \item{quanti}{A number lower than 1 indicating the quantile to perform the computation of local dimension and theta} +\item{spatial_dims}{could be 1 or 2. If more than 1, they will be merged.} + +\item{time_dims}{could be 1 or more. If more than 1, they will be merged} + +\item{na.rm}{logical indicating whether to remove the NAs when calculating the logaritmic distance.} + \item{ncores}{The number of cores to use in parallel computation.} } \value{