Commit bed79371 authored by aho's avatar aho
Browse files

Merge branch 'master' into develop-ACC

parents 71f3967b 9edfc681
Pipeline #5074 passed with stage
in 3 minutes and 9 seconds
......@@ -45,4 +45,4 @@ BugReports: https://earth.bsc.es/gitlab/es/s2dv/-/issues
LazyData: true
SystemRequirements: cdo
Encoding: UTF-8
RoxygenNote: 5.0.0
RoxygenNote: 7.0.1
# s2dv 1.0.0 (Release date: 2021-)
- Add parameter 'memb_dim' and 'memb' in Corr(). They allow the existence of the member dimension
which can have different length between exp and obs, and users can choose to do the ensemble mean
first before correlation or calculate the correlation for individual member.
# s2dv 0.1.1 (Release date: 2020-11-16)
- Change the lincense to Apache License 2.0.
......
......@@ -55,6 +55,8 @@
#'@param member_dim A character string indicating the name of the member
#' dimension. The default value is 'member'. Only used if parameter 'type' is
#' 'dcpp' or 'hist'.
#'@param ncores An integer indicating the number of cores to use for parallel
#' computation. The default value is NULL.
#'
#'@return A numerical array of the AMV index with the dimensions of:
#' 1) sdate, forecast year, and member (in case of decadal predictions);
......@@ -86,7 +88,7 @@
AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lon',
mask = NULL, monini = 11, fmonth_dim = 'fmonth', sdate_dim = 'sdate',
indices_for_clim = NULL, year_dim = 'year', month_dim = 'month',
member_dim = 'member') {
member_dim = 'member', ncores = NULL) {
## Input Checks
# data
......@@ -130,6 +132,13 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo
stop(paste0("The longitude dimension of parameter 'data' must be the same",
" length of parameter 'data_lons'."))
}
# ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
length(ncores) > 1) {
stop("Parameter 'ncores' must be a positive integer.")
}
}
# mask
if (!is.null(mask)) {
if (is.array(mask) & identical(names(dim(mask)), c(lat_dim,lon_dim)) &
......@@ -141,7 +150,7 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo
return(data)
}
data <- multiApply::Apply(data = data, target_dims = c(lat_dim, lon_dim),
fun = fun_mask, mask = mask)$output1
fun = fun_mask, mask = mask, ncores = ncores)$output1
} else {
stop(paste0("Parameter 'mask' must be NULL (no mask) or a numerical array ",
"with c(lat_dim, lon_dim) dimensions and 0 in those grid ",
......@@ -209,7 +218,7 @@ AMV <- function(data, data_lats, data_lons, type, lat_dim = 'lat', lon_dim = 'lo
stop("Parameter 'member_dim' is not found in 'data' dimension.")
}
}
## Regions for AMV (Doblas-Reyes et al., 2013)
lat_min_1 <- 0; lat_max_1 <- 60
lon_min_1 <- 280; lon_max_1 <- 359.9
......
......@@ -74,7 +74,7 @@
}
## ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 |
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
length(ncores) > 1) {
stop("Parameter 'ncores' must be a positive integer.")
}
......
......@@ -136,7 +136,7 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
}
## ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 |
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
length(ncores) > 1) {
stop("Parameter 'ncores' must be a positive integer.")
}
......@@ -172,13 +172,13 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
## dat_dim: [dataset, member]
pos[i] <- which(names(dim(obs)) == dat_dim[i])
}
outrows_exp <- MeanDims(exp, pos, na.rm = FALSE) +
MeanDims(obs, pos, na.rm = FALSE)
outrows_exp <- MeanDims(exp, pos, na.rm = FALSE, ncores = ncores) +
MeanDims(obs, pos, na.rm = FALSE, ncores = ncores)
outrows_obs <- outrows_exp
for (i in 1:length(pos)) {
outrows_exp <- InsertDim(outrows_exp, pos[i], dim(exp)[pos[i]])
outrows_obs <- InsertDim(outrows_obs, pos[i], dim(obs)[pos[i]])
outrows_exp <- InsertDim(outrows_exp, pos[i], dim(exp)[pos[i]], ncores = ncores)
outrows_obs <- InsertDim(outrows_obs, pos[i], dim(obs)[pos[i]], ncores = ncores)
}
exp[which(is.na(outrows_exp))] <- NA
obs[which(is.na(outrows_obs))] <- NA
......@@ -191,7 +191,7 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
fun = .Clim,
method = method, time_dim = time_dim,
dat_dim = dat_dim, memb_dim = memb_dim,
memb = memb, na.rm = na.rm,
memb = memb, na.rm = na.rm, ncores_input = ncores,
ncores = ncores)
# Add member dimension name back
if (memb) {
......@@ -207,7 +207,7 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
fun = .Clim,
method = method, time_dim = time_dim,
dat_dim = dat_dim, ftime_dim = ftime_dim, memb_dim = memb_dim,
memb = memb, na.rm = na.rm,
memb = memb, na.rm = na.rm, ncores_input = ncores,
ncores = ncores)
} else if (method == 'NDV') {
......@@ -216,7 +216,7 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
fun = .Clim,
method = method, time_dim = time_dim,
dat_dim = dat_dim, ftime_dim = ftime_dim, memb_dim = memb_dim,
memb = memb, na.rm = na.rm,
memb = memb, na.rm = na.rm, ncores_input = ncores,
ncores = ncores)
}
......@@ -227,7 +227,7 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
.Clim <- function(exp, obs, method = 'clim',
time_dim = 'sdate', dat_dim = c('dataset', 'member'),
ftime_dim = 'ftime', memb_dim = 'member', memb = TRUE,
na.rm = TRUE) {
na.rm = TRUE, ncores_input = NULL) {
if (method == 'clim') {
# exp: [sdate, dat_dim_exp]
......@@ -269,9 +269,9 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
# exp clim
##--- NEW trend ---##
tmp_obs <- Trend(data = obs, time_dim = time_dim, interval = 1,
polydeg = 1, conf = FALSE)$trend
polydeg = 1, conf = FALSE, ncores = ncores_input)$trend
tmp_exp <- Trend(data = exp, time_dim = time_dim, interval = 1,
polydeg = 1, conf = FALSE)$trend
polydeg = 1, conf = FALSE, ncores = ncores_input)$trend
# tmp_exp: [stats, dat_dim)]
tmp_obs_mean <- apply(tmp_obs, 1, mean) #average out dat_dim (dat and member)
......@@ -331,16 +331,16 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
# Create initial data set (i.e., only first ftime)
tmp <- Subset(exp, ftime_dim, 1, drop = 'selected')
ini_exp <- InsertDim(tmp, pos_ftime, dim_ftime) #only first ftime
ini_exp <- InsertDim(tmp, pos_ftime, dim_ftime, ncores = ncores_input) #only first ftime
tmp <- Subset(obs, ftime_dim, 1, drop = 'selected')
ini_obs <- InsertDim(tmp, pos_ftime, dim_ftime) #only first ftime
ini_obs <- InsertDim(tmp, pos_ftime, dim_ftime, ncores = ncores_input) #only first ftime
#ini_: [sdate, dat_dim, ftime]
tmp_exp <- Regression(datay = exp, datax = ini_exp, reg_dim = time_dim,
na.action = na.omit,
pval = FALSE, conf = FALSE)$regression
pval = FALSE, conf = FALSE, ncores = ncores_input)$regression
tmp_obs <- Regression(datay = obs, datax = ini_obs, reg_dim = time_dim,
na.action = na.omit,
pval = FALSE, conf = FALSE)$regression
pval = FALSE, conf = FALSE, ncores = ncores_input)$regression
#tmp_: [stats = 2, dat_dim, ftime]
tmp_obs_mean <- apply(tmp_obs, c(1, length(dim(tmp_obs))), mean) #average out dat_dim (dat and member)
#tmp_obs_mean: [stats = 2, ftime]
......
......@@ -157,7 +157,7 @@ Composite <- function(data, occ, time_dim = 'time', space_dim = c('lon', 'lat'),
}
## ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 |
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
length(ncores) > 1) {
stop("Parameter 'ncores' must be a positive integer.")
}
......@@ -181,7 +181,7 @@ Composite <- function(data, occ, time_dim = 'time', space_dim = c('lon', 'lat'),
fun = .Composite,
output_dims = output_dims,
occ = occ, time_dim = time_dim, space_dim = space_dim,
K = K, lag = lag, eno = eno,
K = K, lag = lag, eno = eno, ncores_input = ncores,
ncores = ncores)
if (!is.null(fileout)) {
......@@ -192,7 +192,7 @@ Composite <- function(data, occ, time_dim = 'time', space_dim = c('lon', 'lat'),
}
.Composite <- function(data, occ, time_dim = 'time', space_dim = c('lon', 'lat'),
K = NULL, lag = 0, eno = FALSE) {
K = NULL, lag = 0, eno = FALSE, ncores_input = NULL) {
# data: [lon, lat, time]
# occ: [time]
if (is.null(K)) {
......@@ -204,12 +204,12 @@ Composite <- function(data, occ, time_dim = 'time', space_dim = c('lon', 'lat'),
pval <- array(dim = c(dim(data)[1:2], composite = K))
if (eno == TRUE) {
n_tot <- Eno(data, time_dim = time_dim)
n_tot <- Eno(data, time_dim = time_dim, ncores = ncores_input)
} else {
n_tot <- length(occ)
}
mean_tot <- MeanDims(data, dims = 3, na.rm = TRUE)
mean_tot <- MeanDims(data, dims = 3, na.rm = TRUE, ncores = ncores_input)
stdv_tot <- apply(data, c(1, 2), sd, na.rm = TRUE)
for (k in 1:K) {
......@@ -224,14 +224,14 @@ Composite <- function(data, occ, time_dim = 'time', space_dim = c('lon', 'lat'),
if (eno == TRUE) {
data_tmp <- data[, , indices]
names(dim(data_tmp)) <- names(dim(data))
n_k <- Eno(data_tmp, time_dim = time_dim)
n_k <- Eno(data_tmp, time_dim = time_dim, ncores = ncores_input)
} else {
n_k <- length(indices)
}
if (length(indices) == 1) {
composite[, , k] <- data[, , indices]
} else {
composite[, , k] <- MeanDims(data[, , indices], dims = 3, na.rm = TRUE)
composite[, , k] <- MeanDims(data[, , indices], dims = 3, na.rm = TRUE, ncores = ncores_input)
}
stdv_k <- apply(data[, , indices], c(1, 2), sd, na.rm = TRUE)
......
......@@ -19,7 +19,7 @@
#'@param exp A named numeric array of experimental data, with at least two
#' dimensions 'time_dim' and 'dat_dim'.
#'@param obs A named numeric array of observational data, same dimensions as
#' parameter 'exp' except along dat_dim.
#' parameter 'exp' except along 'dat_dim' and 'memb_dim'.
#'@param time_dim A character string indicating the name of dimension along
#' which the correlations are computed. The default value is 'sdate'.
#'@param dat_dim A character string indicating the name of dataset (nobs/nexp)
......@@ -31,6 +31,12 @@
#' be completed. The default is c(1, length(comp_dim dimension)).
#'@param method A character string indicating the type of correlation:
#' 'pearson', 'spearman', or 'kendall'. The default value is 'pearson'.
#'@param memb_dim A character string indicating the name of the member
#' dimension. It must be one dimension in 'exp' and 'obs'. If there is no
#' member dimension, set NULL. The default value is NULL.
#'@param memb A logical value indicating whether to remain 'memb_dim' dimension
#' (TRUE) or do ensemble mean over 'memb_dim' (FALSE). Only functional when
#' 'memb_dim' is not NULL. The default value is TRUE.
#'@param pval A logical value indicating whether to compute or not the p-value
#' of the test Ho: Corr = 0. The default value is TRUE.
#'@param conf A logical value indicating whether to retrieve the confidence
......@@ -42,9 +48,12 @@
#'
#'@return
#'A list containing the numeric arrays with dimension:\cr
#' c(nexp, nobs, all other dimensions of exp except time_dim).\cr
#'nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the
#'number of observation (i.e., dat_dim in obs).\cr
#' c(nexp, nobs, exp_memb, obs_memb, all other dimensions of exp except
#' time_dim and memb_dim).\cr
#'nexp is the number of experiment (i.e., 'dat_dim' in exp), and nobs is the
#'number of observation (i.e., 'dat_dim' in obs). exp_memb is the number of
#'member in experiment (i.e., 'memb_dim' in exp) and obs_memb is the number of
#'member in observation (i.e., 'memb_dim' in obs).\cr\cr
#'\item{$corr}{
#' The correlation coefficient.
#'}
......@@ -59,20 +68,37 @@
#'}
#'
#'@examples
#'# Load sample data as in Load() example:
#'# Case 1: Load sample data as in Load() example:
#'example(Load)
#'clim <- Clim(sampleData$mod, sampleData$obs)
#'corr <- Corr(clim$clim_exp, clim$clim_obs, time_dim = 'ftime', dat_dim = 'member')
#'# Renew the example when Ano and Smoothing is ready
#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
#'runmean_months <- 12
#'
#'# Smooth along lead-times
#'smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
#'smooth_ano_obs <- Smoothing(ano_obs, runmeanlen = runmean_months)
#'required_complete_row <- 3 # Discard start dates which contain any NA lead-times
#'leadtimes_per_startdate <- 60
#'corr <- Corr(MeanDims(smooth_ano_exp, 'member'),
#' MeanDims(smooth_ano_obs, 'member'),
#' comp_dim = 'ftime',
#' limits = c(ceiling((runmean_months + 1) / 2),
#' leadtimes_per_startdate - floor(runmean_months / 2)))
#'
#'# Case 2: Keep member dimension
#'corr <- Corr(smooth_ano_exp, smooth_ano_obs, memb_dim = 'member')
#'# ensemble mean
#'corr <- Corr(smooth_ano_exp, smooth_ano_obs, memb_dim = 'member', memb = FALSE)
#'
#'@rdname Corr
#'@import multiApply
#'@importFrom ClimProjDiags Subset
#'@importFrom stats cor pt qnorm
#'@export
Corr <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
comp_dim = NULL, limits = NULL,
method = 'pearson', pval = TRUE, conf = TRUE,
comp_dim = NULL, limits = NULL, method = 'pearson',
memb_dim = NULL, memb = TRUE,
pval = TRUE, conf = TRUE,
conf.lev = 0.95, ncores = NULL) {
# Check inputs
......@@ -133,6 +159,19 @@ Corr <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
if (!(method %in% c("kendall", "spearman", "pearson"))) {
stop("Parameter 'method' must be one of 'kendall', 'spearman' or 'pearson'.")
}
## memb_dim
if (!is.null(memb_dim)) {
if (!is.character(memb_dim) | length(memb_dim) > 1) {
stop("Parameter 'memb_dim' must be a character string.")
}
if (!memb_dim %in% names(dim(exp)) | !memb_dim %in% names(dim(obs))) {
stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
}
}
## memb
if (!is.logical(memb) | length(memb) > 1) {
stop("Parameter 'memb' must be one logical value.")
}
## pval
if (!is.logical(pval) | length(pval) > 1) {
stop("Parameter 'pval' must be one logical value.")
......@@ -147,7 +186,7 @@ Corr <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
}
## ncores
if (!is.null(ncores)) {
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 |
if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
length(ncores) > 1) {
stop("Parameter 'ncores' must be a positive integer.")
}
......@@ -157,9 +196,13 @@ Corr <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
name_obs <- sort(names(dim(obs)))
name_exp <- name_exp[-which(name_exp == dat_dim)]
name_obs <- name_obs[-which(name_obs == dat_dim)]
if (!is.null(memb_dim)) {
name_exp <- name_exp[-which(name_exp == memb_dim)]
name_obs <- name_obs[-which(name_obs == memb_dim)]
}
if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
"all dimension expect 'dat_dim'."))
"all dimension expect 'dat_dim' and 'memb_dim'."))
}
if (dim(exp)[time_dim] < 3) {
stop("The length of time_dim must be at least 3 to compute correlation.")
......@@ -184,48 +227,161 @@ Corr <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
}
pos <- which(names(dim(obs)) == comp_dim)
obs_sub <- Subset(obs, pos, list(limits[1]:limits[2]))
outrows <- is.na(MeanDims(obs_sub, pos, na.rm = FALSE))
outrows <- InsertDim(outrows, pos, dim(obs)[comp_dim])
outrows <- is.na(MeanDims(obs_sub, pos, na.rm = FALSE, ncores = ncores))
outrows <- InsertDim(outrows, pos, dim(obs)[comp_dim], ncores = ncores)
obs[which(outrows)] <- NA
}
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim, dat_dim),
c(time_dim, dat_dim)),
fun = .Corr,
time_dim = time_dim, method = method,
pval = pval, conf = conf, conf.lev = conf.lev,
ncores = ncores)
if (is.null(memb_dim)) {
# Define output_dims
if (conf & pval) {
output_dims <- list(corr = c('nexp', 'nobs'),
p.val = c('nexp', 'nobs'),
conf.lower = c('nexp', 'nobs'),
conf.upper = c('nexp', 'nobs'))
} else if (conf & !pval) {
output_dims <- list(corr = c('nexp', 'nobs'),
conf.lower = c('nexp', 'nobs'),
conf.upper = c('nexp', 'nobs'))
} else if (!conf & pval) {
output_dims <- list(corr = c('nexp', 'nobs'),
p.val = c('nexp', 'nobs'))
} else {
output_dims <- list(corr = c('nexp', 'nobs'))
}
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim, dat_dim),
c(time_dim, dat_dim)),
output_dims = output_dims,
fun = .Corr,
time_dim = time_dim, method = method,
pval = pval, conf = conf, conf.lev = conf.lev,
ncores = ncores)
} else {
if (!memb) { #ensemble mean
name_exp <- names(dim(exp))
margin_dims_ind <- c(1:length(name_exp))[-which(name_exp == memb_dim)]
exp <- apply(exp, margin_dims_ind, mean, na.rm = TRUE) #NOTE: remove NAs here
obs <- apply(obs, margin_dims_ind, mean, na.rm = TRUE)
# Define output_dims
if (conf & pval) {
output_dims <- list(corr = c('nexp', 'nobs'),
p.val = c('nexp', 'nobs'),
conf.lower = c('nexp', 'nobs'),
conf.upper = c('nexp', 'nobs'))
} else if (conf & !pval) {
output_dims <- list(corr = c('nexp', 'nobs'),
conf.lower = c('nexp', 'nobs'),
conf.upper = c('nexp', 'nobs'))
} else if (!conf & pval) {
output_dims <- list(corr = c('nexp', 'nobs'),
p.val = c('nexp', 'nobs'))
} else {
output_dims <- list(corr = c('nexp', 'nobs'))
}
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim, dat_dim),
c(time_dim, dat_dim)),
output_dims = output_dims,
fun = .Corr,
time_dim = time_dim, method = method,
pval = pval, conf = conf, conf.lev = conf.lev, ncores_input = ncores,
ncores = ncores)
} else { # correlation for each member
# Define output_dims
if (conf & pval) {
output_dims <- list(corr = c('nexp', 'nobs', 'exp_memb', 'obs_memb'),
p.val = c('nexp', 'nobs', 'exp_memb', 'obs_memb'),
conf.lower = c('nexp', 'nobs', 'exp_memb', 'obs_memb'),
conf.upper = c('nexp', 'nobs', 'exp_memb', 'obs_memb'))
} else if (conf & !pval) {
output_dims <- list(corr = c('nexp', 'nobs', 'exp_memb', 'obs_memb'),
conf.lower = c('nexp', 'nobs', 'exp_memb', 'obs_memb'),
conf.upper = c('nexp', 'nobs', 'exp_memb', 'obs_memb'))
} else if (!conf & pval) {
output_dims <- list(corr = c('nexp', 'nobs', 'exp_memb', 'obs_memb'),
p.val = c('nexp', 'nobs', 'exp_memb', 'obs_memb'))
} else {
output_dims <- list(corr = c('nexp', 'nobs', 'exp_memb', 'obs_memb'))
}
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim, dat_dim, memb_dim),
c(time_dim, dat_dim, memb_dim)),
output_dims = output_dims,
fun = .Corr,
time_dim = time_dim, method = method,
pval = pval, conf = conf, conf.lev = conf.lev, ncores_input = ncores,
ncores = ncores)
}
}
return(res)
}
.Corr <- function(exp, obs, time_dim = 'sdate', method = 'pearson',
conf = TRUE, pval = TRUE, conf.lev = 0.95) {
conf = TRUE, pval = TRUE, conf.lev = 0.95, ncores_input = NULL) {
if (length(dim(exp)) == 2) {
# exp: [sdate, dat_exp]
# obs: [sdate, dat_obs]
nexp <- as.numeric(dim(exp)[2])
nobs <- as.numeric(dim(obs)[2])
CORR <- array(dim = c(nexp = nexp, nobs = nobs))
eno_expand <- array(dim = c(nexp = nexp, nobs = nobs))
p.val <- array(dim = c(nexp = nexp, nobs = nobs))
# ens_mean
for (i in 1:nobs) {
CORR[, i] <- sapply(1:nexp,
function(x) {
if (any(!is.na(exp[, x])) && sum(!is.na(obs[, i])) > 2) {
cor(exp[, x], obs[, i],
use = "pairwise.complete.obs",
method = method)
} else {
CORR[, i] <- NA
}
})
# NOTE: Use sapply to replace the for loop
CORR <- sapply(1:nobs, function(i) {
sapply(1:nexp, function (x) {
if (any(!is.na(exp[, x])) && sum(!is.na(obs[, i])) > 2) { #NOTE: Is this necessary?
cor(exp[, x], obs[, i],
use = "pairwise.complete.obs",
method = method)
} else {
NA #CORR[, i] <- NA
}
})
})
if (is.null(dim(CORR))) {
CORR <- array(CORR, dim = c(1, 1))
}
} else { # member
# exp: [sdate, dat_exp, memb_exp]
# obs: [sdate, dat_obs, memb_obs]
nexp <- as.numeric(dim(exp)[2])
nobs <- as.numeric(dim(obs)[2])
exp_memb <- as.numeric(dim(exp)[3])
obs_memb <- as.numeric(dim(obs)[3])
CORR <- array(dim = c(nexp = nexp, nobs = nobs, exp_memb = exp_memb, obs_memb = obs_memb))
for (j in 1:obs_memb) {
for (y in 1:exp_memb) {
CORR[, , y, j] <- sapply(1:nobs, function(i) {
sapply(1:nexp, function (x) {
if (any(!is.na(exp[, x, y])) && sum(!is.na(obs[, i, j])) > 2) {
cor(exp[, x, y], obs[, i, j],
use = "pairwise.complete.obs",
method = method)
} else {
NA #CORR[, i] <- NA
}
})
})
}
}
}
# if (pval) {
# for (i in 1:nobs) {
# p.val[, i] <- try(sapply(1:nexp,
......@@ -240,16 +396,29 @@ cor(exp[, x], obs[, i],
if (pval | conf) {
if (method == "kendall" | method == "spearman") {
tmp <- apply(obs, 2, rank)
tmp <- apply(obs, c(1:length(dim(obs)))[-1], rank) # for memb_dim = NULL, 2; for memb_dim, c(2, 3)
names(dim(tmp))[1] <- time_dim
eno <- Eno(tmp, time_dim)
eno <- Eno(tmp, time_dim, ncores = ncores_input)
} else if (method == "pearson") {
eno <- Eno(obs, time_dim)
eno <- Eno(obs, time_dim, ncores = ncores_input)
}
for (i in 1:nexp) {
eno_expand[i, ] <- eno
if (length(dim(exp)) == 2) {
eno_expand <- array(dim = c(nexp = nexp, nobs = nobs))
for (i in 1:nexp) {
eno_expand[i, ] <- eno
}
} else { #member
eno_expand <- array(dim = c(nexp = nexp, nobs = nobs, exp_memb = exp_memb, obs_memb = obs_memb))
for (i in 1:nexp) {
for (j in 1:exp_memb) {
eno_expand[i, , j, ] <- eno