diff --git a/DESCRIPTION b/DESCRIPTION index 0045265cfeeefd668ffd0fc634ae9e42cac29915..6780c83d1760c794b0289f671433edacbab07636 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -50,7 +50,8 @@ Imports: ncdf4, parallel, plyr, - SpecsVerification (>= 0.5.0) + SpecsVerification (>= 0.5.0), + multiApply (>= 2.0.0) Suggests: easyVerification, testthat diff --git a/NAMESPACE b/NAMESPACE index 4e8557898444268933dfffceb4f5d7d3c8460183..532436773fa75c03b1c152975815a28a1c7ae770 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,7 +3,6 @@ export(.BrierScore) export(.Corr) export(.RMS) -export(.RMSSS) export(.RatioRMS) export(.RatioSDRMS) export(.Trend) @@ -90,6 +89,7 @@ import(graphics) import(mapproj) import(maps) import(methods) +import(multiApply) import(ncdf4) import(parallel) import(plyr) diff --git a/R/RMSSS.R b/R/RMSSS.R index 0ffe36f687f684d3f5d32f3569d79a7a0cd9f866..76d0dda94f4ec87adcd721388b0535d4bb5020b6 100644 --- a/R/RMSSS.R +++ b/R/RMSSS.R @@ -1,191 +1,210 @@ -#'Computes Root Mean Square Skill Score +#'Compute root mean square error skill score #' -#'Computes the root mean square error skill score between an array of -#'forecasts, var_exp and an array of observations, var_obs, which should -#'have the same dimensions except along posloop where the lengths can be -#'different, with the number of experiments/models for var_exp (nexp) and -#'the number of obserational datasets for var_obs (nobs).\cr -#'RMSSS computes the Root Mean Square Skill Score of each jexp in 1:nexp -#'against each jobs in 1:nobs which gives nexp x nobs RMSSS for each other -#'grid point of the matrix (each latitude/longitude/level/leadtime).\cr -#'The RMSSS are computed along the posRMS dimension which should correspond +#'Compute the root mean square error skill score (RMSSS) between an array of +#'forecast 'exp' and an array of observation 'obs'. The two arrays should +#'have the same dimensions except along memb_dim, where the length can be +#'different, with the number of experiments/models (nexp) and the number of +#'observational datasets (nobs).\cr +#'RMSSS computes the root mean square error skill score of each jexp in 1:nexp +#'against each jobs in 1:nobs which gives nexp * nobs RMSSS for each other +#'grid point of the array.\cr +#'The RMSSS are computed along the time_dim dimension which should corresponds #'to the startdate dimension.\cr -#'The p-value is optionally provided by a one-sided Fisher test.\cr\cr -#'.RMSSS provides the same functionality but taking a matrix of ensemble -#'members as input (exp). +#'The p-value is optionally provided by an one-sided Fisher test.\cr #' -#'@param var_exp Array of experimental data. -#'@param var_obs Array of observational data, same dimensions as var_exp -#' except along posloop dimension, where the length can be nobs instead of -#' nexp. -#'@param posloop Dimension nobs and nexp. -#'@param posRMS Dimension along which the RMSE are to be computed (the -#' dimension of the start dates). -#'@param pval Whether to compute or not the p-value of the test Ho : -#' RMSSS = 0. TRUE by default. -#'@param exp N by M matrix of N forecasts from M ensemble members. -#'@param obs Vector of the corresponding observations of length N. +#'@param exp A named numeric array of experimental data which contains at least +#' two dimensions for memb_dim and time_dim. +#'@param obs A named numeric array of observational data which contains at least +#' two dimensions for memb_dim and time_dim. The dimensions should be the same +#' as paramter 'exp' except the length of 'memb_dim' dimension. The order of +#' dimension can be different. +#'@param memb_dim A character string indicating the name of member (nobs/nexp) +#' dimension. The default value is 'member'. +#'@param time_dim A character string indicating the name of dimension along +#' which the RMSSS are computed. The default value is 'sdate'. +#'@param pval A logical value indicating whether to compute or not the p-value +#' of the test Ho: RMSSS = 0. If pval = TRUE, the insignificant RMSSS will +#' return NA. The default value is TRUE. +#'@param ncores An integer indicating the number of cores to use for parallel +#' computation. The default value is NULL. +#' +#'@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., memb_dim in exp), and nobs is the +#'number of observation (i.e., memb_dim in obs).\cr +#'\item{$rmsss}{ +#' The root mean square error skill score. +#'} +#'\item{$p.val}{ +#' The p-value. Only present if \code{pval = TRUE}. +#'} #' -#'@return RMSSS: Array with dimensions:\cr -#' c(length(posloop) in var_exp, length(posloop) in var_obs, 1 or 2, -#' all other dimensions of var_exp & var_obs except posRMS).\cr -#'The 3rd dimension corresponds to the RMSSS and, if \code{pval = TRUE}, -#' the p-value of the one-sided Fisher test with Ho: RMSSS = 0.\cr\cr -#'.RMSSS: -#' \itemize{ -#' \item{$rmsss}{ -#' The RMSSS. -#' } -#' \item{$p_val}{ -#' Corresponds to the p values (only present if \code{pval = TRUE}) for -#' the RMSSS. -#' } -#' } #'@keywords datagen #'@author History:\cr -#'0.1 - 2012-04 (V. Guemas, \email{vguemas@@bsc.es}) - Original code\cr -#'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@@bsc.es}) - Formatting to R CRAN\cr -#'1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter@@bsc.es}) - Adapted to veriApply() +#'0.1 - 2012-04 (V. Guemas, \email{vguemas@bsc.es}) - Original code\cr +#'1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Formatting to R CRAN\cr +#'1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter@bsc.es}) - Adapted to veriApply() +#'3.0 - 2019-12 (A. Ho, \email{an.ho@bsc.es}) - Adopt multiApply feature #'@examples -#'# Load sample data as in Load() example: -#'example(Load) -#'clim <- Clim(sampleData$mod, sampleData$obs) -#'ano_exp <- Ano(sampleData$mod, clim$clim_exp) -#'ano_obs <- Ano(sampleData$obs, clim$clim_obs) -#'rmsss <- RMSSS(Mean1Dim(ano_exp, 2), Mean1Dim(ano_obs, 2)) -#'rmsss_plot <- array(dim = c(dim(rmsss)[1:2], 4, dim(rmsss)[4])) -#'rmsss_plot[, , 2, ] <- rmsss[, , 1, ] -#'rmsss_plot[, , 4, ] <- rmsss[, , 2, ] -#' \donttest{ -#'PlotVsLTime(rmsss_plot, toptitle = "Root Mean Square Skill Score", ytitle = "", -#' monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), -#' listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1), -#' fileout = 'tos_rmsss.eps') -#' } -#'# The following example uses veriApply combined with .RMSSS instead of RMSSS -#' \dontrun{ -#'require(easyVerification) -#'RMSSS2 <- s2dverification:::.RMSSS -#'rmsss2 <- veriApply("RMSSS2", ano_exp, -#' # see ?veriApply for how to use the 'parallel' option -#' Mean1Dim(ano_obs, 2), -#' tdim = 3, ensdim = 2) -#' } +#' set.seed(1) +#' exp <- array(rnorm(15), dim = c(dat = 1, time = 3, member = 5)) +#' set.seed(2) +#' obs <- array(rnorm(6), dim = c(time = 3, member = 2, dat = 1)) +#' res <- RMSSS(exp, obs, time_dim = 'time') +#' #'@rdname RMSSS +#'@import multiApply #'@export -RMSSS <- function(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) { - # - # Enlarge var_exp & var_obs & clim to 10 dim + move posloop & posRMS - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - dimsvar <- dim(var_exp) - for (iind in 1:length(dimsvar)) { - if (iind !=posloop & dim(var_obs)[iind] != dimsvar[iind]) { - stop("var_exp & var_obs must have same dimensions except along posloop") - } +RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', + pval = TRUE, ncores = NULL) { + + # Check inputs + ## exp and obs (1) + if (is.null(exp) | is.null(obs)) { + stop("Parameter 'exp' and 'obs' cannot be NULL.") } - if (dimsvar[posRMS] < 2 ) { - stop("At least 2 values required to compute RMSE") - } - enlvarexp <- Enlarge(var_exp, 10) - enlvarobs <- Enlarge(var_obs, 10) - nexp <- dimsvar[posloop] - nobs <- dim(var_obs)[posloop] - posaperm <- numeric(10) - posaperm[1] <- posloop - posaperm[2] <- posRMS - posaperm[3:10] <- seq(1, 10)[-c(posloop, posRMS)] - permvarexp <- aperm(enlvarexp, posaperm) - permvarobs <- aperm(enlvarobs, posaperm) - dimsaperm <- dim(permvarexp) - # - # RMSSS and its pvalue computation - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - if (pval) { - nvals <- 2 - } else { - nvals <- 1 + if (!is.numeric(exp) | !is.numeric(obs)) { + stop("Parameter 'exp' and 'obs' must be a numeric array.") } - enlRMSSS <- array(dim = c(nexp, nobs, nvals, dimsaperm[3:10])) - - for (jexp in 1:nexp) { - for (jobs in 1:nobs) { - dif1 <- array(dim = dimsaperm[-1]) - dif2 <- array(dim = dimsaperm[-1]) - dif1[, , , , , , , , ] <- permvarexp[jexp, , , , , , , , - , ] - permvarobs[jobs, , , , , , , , , ] - dif2[, , , , , , , , ] <- permvarobs[jobs, , , , , , , , , ] - rms1 <- Mean1Dim(dif1 ** 2, 1, narm = TRUE) ** 0.5 - rms2 <- Mean1Dim(dif2 ** 2, 1, narm = TRUE) ** 0.5 - rms2[which(abs(rms2) <= (max(abs(rms2), na.rm = TRUE) / 1000))] <- max(abs( - rms2), na.rm = TRUE) / 1000 - enlRMSSS[jexp, jobs, 1, , , , , , , , ] <- 1 - (rms1 / rms2) - if (pval) { - eno1 <- Eno(dif1, 1) - eno2 <- Eno(dif2, 1) - F <- (eno2 * (rms2) ** 2 / (eno2 - 1)) / (eno1 * (rms1) ** 2 / (eno1 - 1)) - for (j3 in 1:dimsaperm[3]) { - for (j4 in 1:dimsaperm[4]) { - for (j5 in 1:dimsaperm[5]) { - for (j6 in 1:dimsaperm[6]) { - for (j7 in 1:dimsaperm[7]) { - for (j8 in 1:dimsaperm[8]) { - for (j9 in 1:dimsaperm[9]) { - for (j10 in 1:dimsaperm[10]) { - l1 <- eno1[j3, j4, j5, j6, j7, j8, j9, j10] - l2 <- eno2[j3, j4, j5, j6, j7, j8, j9, j10] - if (is.na(l1) == FALSE & is.na(l2) == FALSE & l1 > 2 & l2 > 2) { - enlRMSSS[jexp, jobs, 2, j3, j4, j5, j6, j7, j8, j9, - j10] <- 1 - pf(F[j3, j4, j5, j6, j7, j8, j9, - j10], l1 - 1, l2 - 1) - } else { - enlRMSSS[jexp, jobs, 1, j3, j4, j5, j6, j7, j8, j9, - j10] <- NA - } - } - } - } - } - } - } - } - } - } + if (is.null(dim(exp)) | is.null(dim(obs))) { + stop(paste0("Parameter 'exp' and 'obs' must be at least two dimensions ", + "containing time_dim and memb_dim.")) + } + if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) | + any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) { + stop("Parameter 'exp' and 'obs' must have dimension names.") + } + if(!all(names(dim(exp)) %in% names(dim(obs))) | + !all(names(dim(obs)) %in% names(dim(exp)))) { + stop("Parameter 'exp' and 'obs' must have same dimension name.") + } + ## time_dim + if (!is.character(time_dim) | length(time_dim) > 1) { + stop("Parameter 'time_dim' must be a character string.") + } + if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) { + stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.") + } + ## 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.") + } + ## pval + if (!is.logical(pval) | length(pval) > 1) { + stop("Parameter 'pval' must be one logical value.") + } + ## 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.") } + } + ## exp and obs (2) + name_exp <- sort(names(dim(exp))) + name_obs <- sort(names(dim(obs))) + 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 'memb_dim'.")) } + if (dim(exp)[time_dim] <= 2) { + stop("The length of time_dim must be more than 2 to compute RMSSS.") + } + + + ############################### + # Sort dimension + name_exp <- names(dim(exp)) + name_obs <- names(dim(obs)) + order_obs <- match(name_exp, name_obs) + obs <- s2dverification:::.aperm2(obs, order_obs) + + + ############################### + # Calculate RMSSS + + res <- Apply(list(exp, obs), + target_dims = list(c(time_dim, memb_dim), + c(time_dim, memb_dim)), + fun = .RMSSS, + time_dim = time_dim, memb_dim = memb_dim, + pval = pval, #conf = conf, conf.lev = conf.lev, + ncores = ncores) - dim(enlRMSSS) <- c(nexp, nobs, nvals, dimsvar[-c(posloop, posRMS)]) - # - # Output - # ~~~~~~~~ - # - enlRMSSS + return(res) } -#'@rdname RMSSS -#'@importFrom stats pf -#'@export -.RMSSS <- function(exp, obs, pval = TRUE) { - dif2 <- obs - dif1 <- rowMeans(exp) - obs - rms1 <- mean(dif1 ** 2, na.rm = TRUE) ** 0.5 - rms2 <- mean(dif2 ** 2, na.rm = TRUE) ** 0.5 - rms2[abs(rms2) <= abs(rms2) / 1000] <- abs(rms2) / 1000 - rmsss <- c(1 - (rms1 / rms2)) - if (pval == TRUE) { - eno1 <- Eno(dif1, 1) - eno2 <- Eno(dif2, 1) - F.stat <- (eno2 * (rms2) ** 2 / (eno2 - 1)) / (eno1 * (rms1) ** 2 / (eno1 - 1)) - if (is.na(eno1) == FALSE && is.na(eno2) == FALSE && eno1 > 2 && eno2 > 2) { - p_val <- 1 - pf(F.stat, eno1 - 1, eno2 - 1) - } + +.RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', pval = TRUE) { + # exp: [sdate, member_exp] + # obs: [sdate, member_obs] + n_exp <- as.numeric(dim(exp)[2]) + n_obs <- as.numeric(dim(obs)[2]) + n_sdate <- as.numeric(dim(exp)[1]) + + p_val <- array(dim = c(nexp = n_exp, nobs = n_obs)) + dif1 <- array(dim = c(n_sdate, n_exp, n_obs)) + names(dim(dif1)) <- c(time_dim, 'nexp', 'nobs') + +# if (conf) { +# conflow <- (1 - conf.lev) / 2 +# confhigh <- 1 - conflow +# conf_low <- array(dim = c(nexp = n_exp, nobs = n_obs)) +# conf_high <- array(dim = c(nexp = n_exp, nobs = n_obs)) +# } + + # dif1 + for (i in 1:n_obs) { + dif1[, , i] <- sapply(1:n_exp, function(x) {exp[, x] - obs[, i]}) + } + + # rms1 and eno1 + rms1 <- apply(dif1^2, c(2, 3), mean, na.rm = TRUE)^0.5 #array(dim = c(n_exp, n_obs)) + # rms2 and eno2 + rms2 <- array(colMeans(obs^2, na.rm = TRUE)^0.5, dim = c(n_obs = n_obs)) + rms2[which(abs(rms2) <= (max(abs(rms2), na.rm = TRUE) / 1000))] <- max(abs( + rms2), na.rm = TRUE) / 1000 + #rms2 above: [nobs] + rms2 <- array(rms2, dim = c(nobs = n_obs, nexp = n_exp)) + #rms2 above: [nobs, nexp] + rms2 <- s2dverification:::.aperm2(rms2, c(2, 1)) + #rms2 above: [nexp, nobs] + + # use rms1 and rms2 to calculate rmsss + rmsss <- 1 - rms1/rms2 + + ## pval and conf + if (pval) { + eno1 <- Eno(dif1, time_dim) + eno2 <- Eno(obs, time_dim) + eno2 <- array(eno2, dim = c(nobs = n_obs, nexp = n_exp)) + eno2 <- s2dverification:::.aperm2(eno2, c(2, 1)) + } + + # pval + if (pval) { + + F.stat <- (eno2 * rms2^2 / (eno2- 1)) / ((eno1 * rms1^2 / (eno1- 1))) + tmp <- !is.na(eno1) & !is.na(eno2) & eno1 > 2 & eno2 > 2 + p_val <- 1 - pf(F.stat, eno1 - 1, eno2 - 1) + p_val[which(!tmp)] <- NA + + # change not significant rmsss to NA + rmsss[which(!tmp)] <- NA + } + + # output + if (pval) { + res <- list(rmsss = rmsss, p.val = p_val) + } else { - p_val <- NA + res <- list(rmsss = rmsss) } - # - # Output - # ~~~~~~~~ - # - list(rmsss = rmsss, p_val = p_val) + return(res) } diff --git a/man/RMSSS.Rd b/man/RMSSS.Rd index f602702bf715d5adca99971942cb90a4facd47b3..240a74633342e0e8a449517443cccba652040efe 100644 --- a/man/RMSSS.Rd +++ b/man/RMSSS.Rd @@ -1,96 +1,73 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/RMSSS.R \name{RMSSS} -\alias{.RMSSS} \alias{RMSSS} -\title{Computes Root Mean Square Skill Score} +\title{Compute root mean square error skill score} \usage{ -RMSSS(var_exp, var_obs, posloop = 1, posRMS = 2, pval = TRUE) - -.RMSSS(exp, obs, pval = TRUE) +RMSSS(exp, obs, time_dim = "sdate", memb_dim = "member", pval = TRUE, + ncores = NULL) } \arguments{ -\item{var_exp}{Array of experimental data.} - -\item{var_obs}{Array of observational data, same dimensions as var_exp -except along posloop dimension, where the length can be nobs instead of -nexp.} +\item{exp}{A named numeric array of experimental data which contains at least +two dimensions for memb_dim and time_dim.} -\item{posloop}{Dimension nobs and nexp.} +\item{obs}{A named numeric array of observational data which contains at least +two dimensions for memb_dim and time_dim. The dimensions should be the same +as paramter 'exp' except the length of 'memb_dim' dimension. The order of +dimension can be different.} -\item{posRMS}{Dimension along which the RMSE are to be computed (the -dimension of the start dates).} +\item{time_dim}{A character string indicating the name of dimension along +which the RMSSS are computed. The default value is 'sdate'.} -\item{pval}{Whether to compute or not the p-value of the test Ho : -RMSSS = 0. TRUE by default.} +\item{memb_dim}{A character string indicating the name of member (nobs/nexp) +dimension. The default value is 'member'.} -\item{exp}{N by M matrix of N forecasts from M ensemble members.} +\item{pval}{A logical value indicating whether to compute or not the p-value +of the test Ho: RMSSS = 0. If pval = TRUE, the insignificant RMSSS will +return NA. The default value is TRUE.} -\item{obs}{Vector of the corresponding observations of length N.} +\item{ncores}{An integer indicating the number of cores to use for parallel +computation. The default value is NULL.} } \value{ -RMSSS: Array with dimensions:\cr - c(length(posloop) in var_exp, length(posloop) in var_obs, 1 or 2, - all other dimensions of var_exp & var_obs except posRMS).\cr -The 3rd dimension corresponds to the RMSSS and, if \code{pval = TRUE}, - the p-value of the one-sided Fisher test with Ho: RMSSS = 0.\cr\cr -.RMSSS: - \itemize{ - \item{$rmsss}{ - The RMSSS. - } - \item{$p_val}{ - Corresponds to the p values (only present if \code{pval = TRUE}) for - the RMSSS. - } - } +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., memb_dim in exp), and nobs is the +number of observation (i.e., memb_dim in obs).\cr +\item{$rmsss}{ + The root mean square error skill score. +} +\item{$p.val}{ + The p-value. Only present if \code{pval = TRUE}. +} } \description{ -Computes the root mean square error skill score between an array of -forecasts, var_exp and an array of observations, var_obs, which should -have the same dimensions except along posloop where the lengths can be -different, with the number of experiments/models for var_exp (nexp) and -the number of obserational datasets for var_obs (nobs).\cr -RMSSS computes the Root Mean Square Skill Score of each jexp in 1:nexp -against each jobs in 1:nobs which gives nexp x nobs RMSSS for each other -grid point of the matrix (each latitude/longitude/level/leadtime).\cr -The RMSSS are computed along the posRMS dimension which should correspond +Compute the root mean square error skill score (RMSSS) between an array of +forecast 'exp' and an array of observation 'obs'. The two arrays should +have the same dimensions except along memb_dim, where the length can be +different, with the number of experiments/models (nexp) and the number of +observational datasets (nobs).\cr +RMSSS computes the root mean square error skill score of each jexp in 1:nexp +against each jobs in 1:nobs which gives nexp * nobs RMSSS for each other +grid point of the array.\cr +The RMSSS are computed along the time_dim dimension which should corresponds to the startdate dimension.\cr -The p-value is optionally provided by a one-sided Fisher test.\cr\cr -.RMSSS provides the same functionality but taking a matrix of ensemble -members as input (exp). +The p-value is optionally provided by an one-sided Fisher test.\cr } \examples{ -# Load sample data as in Load() example: -example(Load) -clim <- Clim(sampleData$mod, sampleData$obs) -ano_exp <- Ano(sampleData$mod, clim$clim_exp) -ano_obs <- Ano(sampleData$obs, clim$clim_obs) -rmsss <- RMSSS(Mean1Dim(ano_exp, 2), Mean1Dim(ano_obs, 2)) -rmsss_plot <- array(dim = c(dim(rmsss)[1:2], 4, dim(rmsss)[4])) -rmsss_plot[, , 2, ] <- rmsss[, , 1, ] -rmsss_plot[, , 4, ] <- rmsss[, , 2, ] - \donttest{ -PlotVsLTime(rmsss_plot, toptitle = "Root Mean Square Skill Score", ytitle = "", - monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), - listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1), - fileout = 'tos_rmsss.eps') - } -# The following example uses veriApply combined with .RMSSS instead of RMSSS - \dontrun{ -require(easyVerification) -RMSSS2 <- s2dverification:::.RMSSS -rmsss2 <- veriApply("RMSSS2", ano_exp, - # see ?veriApply for how to use the 'parallel' option - Mean1Dim(ano_obs, 2), - tdim = 3, ensdim = 2) - } +set.seed(1) +exp <- array(rnorm(15), dim = c(dat = 1, time = 3, member = 5)) +set.seed(2) +obs <- array(rnorm(6), dim = c(time = 3, member = 2, dat = 1)) +res <- RMSSS(exp, obs, time_dim = 'time') + } \author{ History:\cr 0.1 - 2012-04 (V. Guemas, \email{vguemas@bsc.es}) - Original code\cr 1.0 - 2013-09 (N. Manubens, \email{nicolau.manubens@bsc.es}) - Formatting to R CRAN\cr 1.1 - 2017-02 (A. Hunter, \email{alasdair.hunter@bsc.es}) - Adapted to veriApply() +3.0 - 2019-12 (A. Ho, \email{an.ho@bsc.es}) - Adopt multiApply feature } \keyword{datagen} diff --git a/tests/testthat/test-RMSSS.R b/tests/testthat/test-RMSSS.R new file mode 100644 index 0000000000000000000000000000000000000000..a73f90e71ce0ac67b4e9252dadb81d94e1fba9fe --- /dev/null +++ b/tests/testthat/test-RMSSS.R @@ -0,0 +1,142 @@ +context("s2dverification::RMSSS tests") + +############################################## + # case 0 + set.seed(1) + exp0 <- array(rnorm(15), dim = c(sdate = 3, member = 5)) + set.seed(2) + obs0 <- array(rnorm(6), dim = c(sdate = 3, member = 2)) + + # case 1 + set.seed(1) + exp1 <- array(rnorm(15), dim = c(time = 3, memb = 5)) + set.seed(2) + obs1 <- array(rnorm(6), dim = c(time = 3, memb = 2)) + + # case 2 + set.seed(3) + exp2 <- array(rnorm(120), dim = c(sdate = 10, dat = 1, lon = 3, lat = 2, member = 2)) + set.seed(4) + obs2 <- array(rnorm(60), dim = c(dat = 1, sdate = 10, member = 1, lat = 2, lon = 3)) + +############################################## + +test_that("1. Input checks", { + + expect_error( + RMSSS(c(), c()), + "Parameter 'exp' and 'obs' cannot be NULL." + ) + expect_error( + RMSSS('exp', 'obs'), + "Parameter 'exp' and 'obs' must be a numeric array." + ) + expect_error( + RMSSS(c(1:10), c(2:4)), + paste0("Parameter 'exp' and 'obs' must be at least two dimensions ", + "containing time_dim and memb_dim.") + ) + expect_error( + RMSSS(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))), + "Parameter 'exp' and 'obs' must have dimension names." + ) + expect_error( + RMSSS(array(1:10, dim = c(a = 2, c = 5)), array(1:4, dim = c(a = 2, b = 2))), + "Parameter 'exp' and 'obs' must have same dimension name" + ) + expect_error( + RMSSS(exp1, obs1, time_dim = 1), + "Parameter 'time_dim' must be a character string." + ) + expect_error( + RMSSS(exp0, obs0, time_dim = 'a'), + "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension." + ) + expect_error( + RMSSS(exp0, obs0, memb_dim = NA), + "Parameter 'memb_dim' must be a character string." + ) + expect_error( + RMSSS(exp0, obs0, memb_dim = 'memb'), + "Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension." + ) + expect_error( + RMSSS(exp0, obs0, pval = c(T, T)), + "Parameter 'pval' must be one logical value." + ) + expect_error( + RMSSS(exp0, obs0, ncores = 1.4), + "Parameter 'ncores' must be a positive integer." + ) + expect_error( + RMSSS(exp = array(1:10, dim = c(sdate = 1, member = 5, a = 1)), + obs = array(1:4, dim = c(a = 1, sdate = 2, member = 2))), + "Parameter 'exp' and 'obs' must have same length of all dimension expect 'memb_dim'." + ) + expect_error( + RMSSS(exp = array(1:10, dim = c(sdate = 1, member = 5, a = 1)), + obs = array(1:4, dim = c(a = 1, sdate = 1, member = 2))), + "The length of time_dim must be more than 2 to compute RMSSS." + ) +}) + +############################################## +test_that("1. Output checks: case 1", { + + res1_1 <- RMSSS(exp1, obs1, time_dim = 'time', memb_dim = 'memb') + expect_equal( + dim(res1_1$rmsss), + c(nexp = 5, nobs = 2) + ) + expect_equal( + dim(res1_1$p.val), + c(nexp = 5, nobs = 2) + ) + expect_equal( + mean(res1_1$rmsss), + -0.5449538, + tolerance = 0.00001 + ) + + exp1_2 <- exp1 + exp1_2[2:4] <- NA + obs1_2 <- obs1 + obs1_2[1:2] <- NA + res1_2 <- RMSSS(exp1_2, obs1_2, time_dim = 'time', memb_dim = 'memb', pval = TRUE) + + expect_equal( + length(res1_2$rmsss[which(is.na(res1_2$rmsss))]), + 7 + ) + expect_equal( + range(res1_2$p.val, na.rm = T), + c(0.7159769, 0.8167073), + tolerance = 0.00001 + ) + +}) + + +############################################## +test_that("2. Output checks: case 2", { + + expect_equal( + dim(RMSSS(exp2, obs2)$rmsss), + c(nexp = 2, nobs = 1, dat = 1, lon = 3, lat = 2) + ) + expect_equal( + mean(RMSSS(exp2, obs2)$rmsss), + -0.3912208, + tolerance = 0.00001 + ) + expect_equal( + range(RMSSS(exp2, obs2)$p.val), + c(0.2627770, 0.9868412), + tolerance = 0.00001 + ) + +}) + +############################################## + +