From 87596f2143f221c1cecc80dceb957d9f68b44efd Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 16 Feb 2023 16:51:17 +0100 Subject: [PATCH 1/4] Change check from integer to numeric --- R/NAO.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/NAO.R b/R/NAO.R index 6d48dba..b9b9b22 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -198,7 +198,7 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', } } ## ftime_avg - if (!is.vector(ftime_avg) | !is.integer(ftime_avg)) { + if (!is.vector(ftime_avg) | !is.numeric(ftime_avg)) { stop("Parameter 'ftime_avg' must be an integer vector.") } if (!is.null(exp)) { -- GitLab From 07b13e7590e1f3fde0490ed11dbdf7b5a1d7f088 Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 16 Mar 2023 17:52:50 +0100 Subject: [PATCH 2/4] Allow ftime_avg to be NULL --- R/NAO.R | 46 +++++++++++++++++++++------------------ man/NAO.Rd | 4 ++-- tests/testthat/test-NAO.R | 14 ++++++++++++ 3 files changed, 41 insertions(+), 23 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index b9b9b22..b535ce6 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -33,8 +33,8 @@ #'@param ftime_dim A character string indicating the name of the forecast time #' dimension of 'exp' and 'obs'. The default value is 'ftime'. #'@param ftime_avg A numeric vector of the forecast time steps to average -#' across the target period. The default value is 2:4, i.e., from 2nd to 4th -#' forecast time steps. +#' across the target period. If average is not needed, set NULL. The default +#' value is 2:4, i.e., from 2nd to 4th forecast time steps. #'@param obsproj A logical value indicating whether to compute the NAO index by #' projecting the forecast anomalies onto the leading EOF of observational #' reference (TRUE) or compute the NAO by first computing the leading @@ -198,16 +198,18 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', } } ## ftime_avg - if (!is.vector(ftime_avg) | !is.numeric(ftime_avg)) { - stop("Parameter 'ftime_avg' must be an integer vector.") - } - if (!is.null(exp)) { - if (max(ftime_avg) > dim(exp)[ftime_dim] | min(ftime_avg) < 1) { - stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.") + if (!is.null(ftime_avg)) { + if (!is.vector(ftime_avg) | !is.numeric(ftime_avg)) { + stop("Parameter 'ftime_avg' must be an integer vector.") } - } else { - if (max(ftime_avg) > dim(obs)[ftime_dim] | min(ftime_avg) < 1) { - stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.") + if (!is.null(exp)) { + if (max(ftime_avg) > dim(exp)[ftime_dim] | min(ftime_avg) < 1) { + stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.") + } + } else { + if (max(ftime_avg) > dim(obs)[ftime_dim] | min(ftime_avg) < 1) { + stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.") + } } } ## sdate >= 2 @@ -281,22 +283,24 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate', } ## 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.") } } # Average ftime - if (!is.null(exp)) { - exp_sub <- ClimProjDiags::Subset(exp, ftime_dim, ftime_avg, drop = FALSE) - exp <- MeanDims(exp_sub, ftime_dim, na.rm = TRUE) - ## Cross-validated PCs. Fabian. This should be extended to - ## nmod and nlt by simple loops. Virginie - } - if (!is.null(obs)) { - obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE) - obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE) + if (!is.null(ftime_avg)) { + if (!is.null(exp)) { + exp_sub <- ClimProjDiags::Subset(exp, ftime_dim, ftime_avg, drop = FALSE) + exp <- MeanDims(exp_sub, ftime_dim, na.rm = TRUE) + ## Cross-validated PCs. Fabian. This should be extended to + ## nmod and nlt by simple loops. Virginie + } + if (!is.null(obs)) { + obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE) + obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE) + } } # wght diff --git a/man/NAO.Rd b/man/NAO.Rd index 64b1656..71c0a55 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -50,8 +50,8 @@ of 'ano'. The default value is c('lat', 'lon').} dimension of 'exp' and 'obs'. The default value is 'ftime'.} \item{ftime_avg}{A numeric vector of the forecast time steps to average -across the target period. The default value is 2:4, i.e., from 2nd to 4th -forecast time steps.} +across the target period. If average is not needed, set NULL. The default +value is 2:4, i.e., from 2nd to 4th forecast time steps.} \item{obsproj}{A logical value indicating whether to compute the NAO index by projecting the forecast anomalies onto the leading EOF of observational diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R index 2da0d40..05fcd22 100644 --- a/tests/testthat/test-NAO.R +++ b/tests/testthat/test-NAO.R @@ -16,6 +16,7 @@ context("s2dv::NAO tests") obs2 <- array(rnorm(108), dim = c(sdate = 3, ftime = 4, lat = 3, lon = 3)) lat2 <- c(80, 50, 20) lon2 <- c(-80, 0, 40) + ############################################## test_that("1. Input checks", { @@ -208,6 +209,19 @@ test_that("2. dat1", { c(-0.1139683, 0.1056687, 0.1889449), tolerance = 0.0001 ) + expect_equal( + dim(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp), + c(sdate = 3, member = 2, dataset = 1) + ) + expect_equal( + dim(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = NULL)$exp), + c(sdate = 3, member = 2, dataset = 1, ftime = 4) + ) + expect_equal( + c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = NULL)$exp[,,1,1]), + c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp) + ) + }) ############################################## -- GitLab From d72a040226c3e2d7b176f18001df43aecb4bb37f Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 16 Mar 2023 17:56:50 +0100 Subject: [PATCH 3/4] Revise output doc --- R/NAO.R | 6 ++++-- man/NAO.Rd | 6 ++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/R/NAO.R b/R/NAO.R index b535ce6..fb5220c 100644 --- a/R/NAO.R +++ b/R/NAO.R @@ -48,11 +48,13 @@ #'A list which contains: #'\item{exp}{ #' A numeric array of forecast NAO index in verification format with the same -#' dimensions as 'exp' except space_dim and ftime_dim. +#' dimensions as 'exp' except space_dim and ftime_dim. If ftime_avg is NULL, +#' ftime_dim remains. #' } #'\item{obs}{ #' A numeric array of observed NAO index in verification format with the same -#' dimensions as 'obs' except space_dim and ftime_dim. +#' dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL, +#' ftime_dim remains. #'} #' #'@references diff --git a/man/NAO.Rd b/man/NAO.Rd index 71c0a55..999fd75 100644 --- a/man/NAO.Rd +++ b/man/NAO.Rd @@ -67,11 +67,13 @@ computation. The default value is NULL.} A list which contains: \item{exp}{ A numeric array of forecast NAO index in verification format with the same - dimensions as 'exp' except space_dim and ftime_dim. + dimensions as 'exp' except space_dim and ftime_dim. If ftime_avg is NULL, + ftime_dim remains. } \item{obs}{ A numeric array of observed NAO index in verification format with the same - dimensions as 'obs' except space_dim and ftime_dim. + dimensions as 'obs' except space_dim and ftime_dim. If ftime_avg is NULL, + ftime_dim remains. } } \description{ -- GitLab From e0fbd85e8420517b3be5c85d83c5191ca70499da Mon Sep 17 00:00:00 2001 From: aho Date: Fri, 17 Mar 2023 16:32:32 +0100 Subject: [PATCH 4/4] avoid appy() to save memory --- R/ProjectField.R | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/R/ProjectField.R b/R/ProjectField.R index 309f3ef..55e7fd2 100644 --- a/R/ProjectField.R +++ b/R/ProjectField.R @@ -237,10 +237,14 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon # Weight e.1 <- eof_mode * wght ano <- ano * InsertDim(wght, 1, ntime) + #ano <- aaply(ano, 1, '*', wght) # much heavier - na <- apply(ano, 1, mean, na.rm = TRUE) # if [lat, lon] all NA, it's NA + na <- rowMeans(ano, na.rm = TRUE) # if [lat, lon] all NA, it's NA + #na <- apply(ano, 1, mean, na.rm = TRUE) # much heavier tmp <- ano * InsertDim(e.1, 1, ntime) # [sdate, lat, lon] - pc.ver <- apply(tmp, 1, sum, na.rm = TRUE) + rm(ano) + #pc.ver <- apply(tmp, 1, sum, na.rm = TRUE) # much heavier + pc.ver <- rowSums(tmp, na.rm = TRUE) pc.ver[which(is.na(na))] <- NA } else { # mode = NULL @@ -250,7 +254,7 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon ano <- ano * InsertDim(wght, 1, ntime) dim(ano) <- c(ntime, prod(dim(ano)[2:3])) # [sdate, lat*lon] - na <- apply(ano, 1, mean, na.rm = TRUE) # if [lat, lon] all NA, it's NA + na <- rowMeans(ano, na.rm = TRUE) # if [lat, lon] all NA, it's NA na <- aperm(array(na, dim = c(ntime, dim(e.1)[1])), c(2, 1)) # Matrix multiplication e.1 [mode, lat*lon] by ano [lat*lon, sdate] -- GitLab