From eac70e8eed51315bfe10e131e8672886b6af4dca Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 31 Aug 2022 17:24:30 +0200 Subject: [PATCH 1/6] Parameter out_len to .predictability to fix Apply warnings --- R/CST_ProxiesAttractor.R | 49 +++++++++++++++--------------- R/Predictability.R | 64 +++++++++++++++++++++++----------------- 2 files changed, 61 insertions(+), 52 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 518968fa..04776e74 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -24,6 +24,8 @@ #'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 ncores The number of cores to use in parallel computation #' @@ -34,7 +36,8 @@ #'attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6) #' #'@export -CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ +CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), + time_dims = c('sdate', 'ftime'), ncores = NULL){ if (!inherits(data, 's2dv_cube')) { stop("Parameter 'data' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") @@ -70,6 +73,8 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ #' #'@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 ncores The number of cores to use in parallel computation. #' #'@return dim and theta @@ -94,38 +99,31 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ #' #'@export -ProxiesAttractor <- function(data, quanti, ncores = NULL){ +ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), + time_dims = c('sdate', 'ftime'), 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 (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% 'grid'))){ - stop("Parameter 'data' must have a spatial dimension named 'grid'.") + if (length(time_dims) == 2) { + data <- MergeDims(data, merge_dims = time_dims, rename = 'time') + time_dims <- 'time' + } else if (length(time_dims) > 2) { + stop("Parameter 'time_dims' maximum length 2.") } - attractor <- Apply(data, target_dims = c('time', 'grid'), + attractor <- Apply(data, target_dims = c(time_dims, spatial_dims), fun = .proxiesattractor, quanti = quanti , ncores = ncores) # rename dimensions @@ -134,6 +132,7 @@ ProxiesAttractor <- function(data, quanti, ncores = NULL){ names(dim(x))[dimname] <- 'time' return(x)}, dimname = which(names(dim(attractor[[1]])) == 'dim2')) + return(list(dim = attractor$dim, theta = attractor$theta)) } diff --git a/R/Predictability.R b/R/Predictability.R index 12cd6b41..b2274b48 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -61,14 +61,6 @@ #'@export #' Predictability<- function(dim, theta, ncores = NULL) { - # if (!inherits(dim, 's2dv_cube')) { - # stop("Parameter 'dim' must be of the class 's2dv_cube', ", - # "as output by CSTools::CST_Load.") - # } - # if (!inherits(theta, 's2dv_cube')) { - # stop("Parameter 'theta' must be of the class 's2dv_cube', ", - # "as output by CSTools::CST_Load.") - # } if (any(names(dim(dim)) %in% 'sdate')) { if (any(names(dim(dim)) %in% 'ftime')) { dim <- MergeDims(dim, merge_dims = c('ftime', 'sdate'), @@ -88,29 +80,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) @@ -118,7 +123,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) @@ -128,6 +133,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)] @@ -150,10 +156,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)) } -- GitLab From a80160aaa44fe031177cbdabc32e72e48191424d Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 1 Sep 2022 17:04:51 +0200 Subject: [PATCH 2/6] order dimensions --- R/CST_ProxiesAttractor.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 04776e74..193071eb 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -99,8 +99,8 @@ CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), #' #'@export -ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), - time_dims = c('sdate', 'ftime'), ncores = NULL){ +ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), + time_dims = c('ftime', 'sdate'), ncores = NULL){ if (is.null(data)) { stop("Parameter 'data' cannot be NULL.") } -- GitLab From efd398e454f4d9dd1ba05775e7cbc85babec8d63 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 1 Sep 2022 17:41:43 +0200 Subject: [PATCH 3/6] na.rm for logdista --- R/CST_ProxiesAttractor.R | 23 ++++++++++++++++------- man/CST_ProxiesAttractor.Rd | 15 ++++++++++++++- man/ProxiesAttractor.Rd | 15 ++++++++++++++- 3 files changed, 44 insertions(+), 9 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 193071eb..64d83a30 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -26,7 +26,7 @@ #'@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 @@ -37,7 +37,8 @@ #' #'@export CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), - time_dims = c('sdate', 'ftime'), ncores = NULL){ + 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.") @@ -46,7 +47,10 @@ CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), stop("Parameter 'quanti' cannot be NULL.") } - data$data <- ProxiesAttractor(data = data$data, quanti = quanti, ncores = ncores) + data$data <- ProxiesAttractor(data = data$data, quanti = quanti, + spatial_dims = spatial_dims, + time_dims = time_dims, + na.rm = na.rm, ncores = ncores) return(data) } @@ -75,6 +79,7 @@ CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), #'@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 @@ -100,7 +105,8 @@ CST_ProxiesAttractor <- function(data, quanti, spatial_dims = c('lat', 'lon'), #'@export ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), - time_dims = c('ftime', 'sdate'), ncores = NULL){ + time_dims = c('ftime', 'sdate'), + na.rm = TRUE, ncores = NULL){ if (is.null(data)) { stop("Parameter 'data' cannot be NULL.") } @@ -123,9 +129,12 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), } else if (length(time_dims) > 2) { stop("Parameter 'time_dims' maximum length 2.") } + if (any(is.na(data))) { + warning("Parameter 'data' contains NAs and they will be removed.") + } 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){ @@ -136,11 +145,11 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), 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', 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/man/CST_ProxiesAttractor.Rd b/man/CST_ProxiesAttractor.Rd index ecf9b9da..66784547 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}{a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow @@ -12,6 +19,12 @@ 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} +\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} } \value{ diff --git a/man/ProxiesAttractor.Rd b/man/ProxiesAttractor.Rd index 768ba736..e15d4ab0 100644 --- a/man/ProxiesAttractor.Rd +++ b/man/ProxiesAttractor.Rd @@ -4,13 +4,26 @@ \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 attractor. It requires a temporal dimension named 'time' and spatial 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{ -- GitLab From 73f46b67c9333e498d86b2920ff940d5b3da6130 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 1 Sep 2022 18:09:00 +0200 Subject: [PATCH 4/6] proxies attractor return array with input dimension names --- R/CST_ProxiesAttractor.R | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 64d83a30..a3f7782a 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -124,6 +124,7 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), stop("Paramenter 'spatial_dims' maximum length 2.") } if (length(time_dims) == 2) { + split_time <- dim(data)[time_dims] data <- MergeDims(data, merge_dims = time_dims, rename = 'time') time_dims <- 'time' } else if (length(time_dims) > 2) { @@ -138,10 +139,18 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), # 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[2], + new_dim_name = names(split_time[2])) + attractor <- lapply(attractor, function(x, name) { + names(dim(x))[which(names(dim(x)) == 'time')] <- + names(split_time)[1] + return(x)}) + } return(list(dim = attractor$dim, theta = attractor$theta)) } -- GitLab From b051ec74dc8032671e46ac5b1bc5ef6f8cfe2f08 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 1 Sep 2022 18:31:28 +0200 Subject: [PATCH 5/6] split dims proxiesattractor --- R/CST_ProxiesAttractor.R | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index a3f7782a..41b855aa 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -125,8 +125,8 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), } if (length(time_dims) == 2) { split_time <- dim(data)[time_dims] - data <- MergeDims(data, merge_dims = time_dims, rename = 'time') - time_dims <- 'time' + 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.") } @@ -145,10 +145,10 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), if (length(split_time) == 2) { attractor <- lapply(attractor, SplitDim, split_dim = time_dims, indices = NULL, freq = split_time[2], - new_dim_name = names(split_time[2])) + new_dim_name = names(split_time[1])) attractor <- lapply(attractor, function(x, name) { - names(dim(x))[which(names(dim(x)) == 'time')] <- - names(split_time)[1] + names(dim(x))[which(names(dim(x)) == time_dims)] <- + names(split_time)[2] return(x)}) } return(list(dim = attractor$dim, theta = attractor$theta)) @@ -156,7 +156,8 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), .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, na.rm = na.rm))}, y = t(data))[[1]] -- GitLab From 055929389dbdaa5cae1814cccb3b3c698d99812e Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 1 Sep 2022 19:09:04 +0200 Subject: [PATCH 6/6] improve split --- R/CST_ProxiesAttractor.R | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 41b855aa..454013a9 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -124,7 +124,16 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), stop("Paramenter 'spatial_dims' maximum length 2.") } if (length(time_dims) == 2) { - split_time <- dim(data)[time_dims] + 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) { @@ -144,11 +153,11 @@ ProxiesAttractor <- function(data, quanti, spatial_dims = c('lon', 'lat'), 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[2], - new_dim_name = names(split_time[1])) + 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)[2] + names(split_time)[1] return(x)}) } return(list(dim = attractor$dim, theta = attractor$theta)) -- GitLab