From 2d4ac21fcb25d06dfe2ff82f28c14acef4c71ef0 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Mon, 13 Dec 2021 18:16:19 +0100
Subject: [PATCH 1/8] Add area weighting (undone)

---
 R/ACC.R | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/R/ACC.R b/R/ACC.R
index 0f1040c..9e92da4 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -293,6 +293,9 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
                                     c(space_dim, dat_dim)),
                  fun = .ACC,
                  dat_dim = dat_dim, avg_dim = avg_dim,
+#------------NEW--------------
+                 lat = lat,
+#------------NEW_END--------------
                  conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
                  ncores_input = ncores,
                  ncores = ncores)
@@ -321,6 +324,9 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
                                     c(space_dim, avg_dim, dat_dim)),
                  fun = .ACC,
                  dat_dim = dat_dim, avg_dim = avg_dim,
+#------------NEW--------------
+                 lat = lat,
+#------------NEW_END--------------
                  conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
                  ncores_input = ncores,
                  ncores = ncores)
@@ -409,6 +415,15 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
 
       if (is.null(avg_dim)) {
         # ACC
+#-----------NEW--------------
+        # area weighted
+        wt <- cos(lat * pi/180)
+        wt <- rep(wt, as.numeric(dim(exp_sub)[2])) ## dim(exp): [lat, lon]
+        exp_sub <- as.vector(t(exp_sub)) ## first longitudes and then latitudes
+        obs_sub <- as.vector(t(obs_sub)) ## first longitudes and then latitudes
+        exp_sub <- sqrt(wt) * (exp_sub - mean(exp_sub))
+        obs_sub <- sqrt(wt) * (obs_sub - mean(obs_sub))
+#----------NEW_END-----------
         top <- sum(exp_sub*obs_sub, na.rm = TRUE)  #a number
         bottom <- sqrt(sum(exp_sub^2, na.rm = TRUE) * sum(obs_sub^2, na.rm = TRUE))
         acc[iexp, iobs] <- top/bottom #a number
-- 
GitLab


From 56462824cc0ba8bea4f4807e9cd7a91e4d5b4788 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 16 Dec 2021 14:35:22 +0100
Subject: [PATCH 2/8] Add area weighting and allow dat_dim = NULL

---
 R/ACC.R    | 281 +++++++++++++++++++++++++++++++++--------------------
 man/ACC.Rd |  47 +++++----
 2 files changed, 202 insertions(+), 126 deletions(-)

diff --git a/R/ACC.R b/R/ACC.R
index 9e92da4..340dd77 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -1,37 +1,40 @@
-#'Compute the anomaly correlation coefficient between the forecast and corresponding observation
+#'Compute the spatial anomaly correlation coefficient between the forecast and corresponding observation
 #'
-#'Calculate the anomaly correlation coefficient for the ensemble mean of each 
-#'model and the corresponding references over a spatial domain. It can return a
-#'forecast time series if the data contain forest time dimension, and also the 
-#'start date mean if the data contain start date dimension. 
-#'The domain of interest can be specified by providing the list 
-#'of longitudes/latitudes (lon/lat) of the data together with the corners 
-#'of the domain: lonlatbox = c(lonmin, lonmax, latmin, latmax).
+#'Calculate the spatial anomaly correlation coefficient (ACC) for the ensemble
+#'mean of each model and the corresponding references over a spatial domain. It
+#'can return a forecast time series if the data contain forest time dimension, 
+#'and also the ACC mean over one dimension, e.g., start date dimension.
+#'The domain of interest can be specified by providing the list of longitudes/
+#'latitudes of the data together with the corners of the domain: lonlatbox = 
+#'c(lonmin, lonmax, latmin, latmax). The data will be adjusted to have a spatial
+#'mean of zero, then area weighting is applied. The formula is referenced from 
+#'Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).
 #'
 #'@param exp A numeric array of experimental anomalies with named dimensions.
-#'  It must have at least 'dat_dim' and 'space_dim'.
+#'  The dimension must have at least 'space_dim'.
 #'@param obs A numeric array of observational anomalies with named dimensions.
-#'  It must have the same dimensions as 'exp' except the length of 'dat_dim' 
+#'  The dimension should be the same as 'exp' except the length of 'dat_dim' 
 #'  and 'memb_dim'.
 #'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
-#'  dimension. The default value is 'dataset'.
+#'  dimension. The default value is 'dataset'. If there is no dataset 
+#'  dimension, set NULL.
 #'@param space_dim A character string vector of 2 indicating the name of the
 #'  latitude and longitude dimensions (in order) along which ACC is computed. 
 #'  The default value is c('lat', 'lon').
 #'@param avg_dim A character string indicating the name of the dimension to be
-#'  averaged. It must be one of 'time_dim'. The mean ACC is calculated along 
-#'  averaged. If no need to calculate mean ACC, set as NULL. The default value 
-#'  is 'sdate'.
+#'  averaged, which is usually the time dimension. If no need to calculate mean
+#'  ACC, set as NULL. The default value is 'sdate'.
 #'@param memb_dim A character string indicating the name of the member 
 #'  dimension. If the data are not ensemble ones, set as NULL. The default 
 #'  value is 'member'.
-#'@param lat A vector of the latitudes of the exp/obs grids. Only required when
-#'  the domain of interested is specified. The default value is NULL.
+#'@param lat A vector of the latitudes of the exp/obs grids. It is used for
+#'  area weighting and when the domain of interested 'lonlatbox' is specified.
 #'@param lon A vector of the longitudes of the exp/obs grids. Only required when
-#'  the domain of interested is specified. The default value is NULL.
+#'  the domain of interested 'lonlatbox' is specified. The default value is 
+#'  NULL.
 #'@param lonlatbox A numeric vector of 4 indicating the corners of the domain of
-#'  interested: c(lonmin, lonmax, latmin, latmax). Only required when the domain
-#'  of interested is specified. The default value is NULL.
+#'  interested: c(lonmin, lonmax, latmin, latmax). The default value is NULL 
+#'  and the whole data will be used.
 #'@param conf A logical value indicating whether to retrieve the confidence 
 #'  intervals or not. The default value is TRUE.
 #'@param conftype A charater string of "parametric" or "bootstrap". 
@@ -55,7 +58,8 @@
 #'\item{acc}{
 #'  The ACC with the dimensions c(nexp, nobs, the rest of the dimension except 
 #'  space_dim and memb_dim). 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), and nobs is the number of observation (i.e., dat_dim in obs). If 
+#'  dat_dim is NULL, nexp and nobs are omitted.
 #'}
 #'\item{conf.lower (if conftype = "parametric") or acc_conf.lower (if 
 #'      conftype = "bootstrap")}{
@@ -113,7 +117,8 @@
 #'PlotACC(res_bootstrap, startDates)
 #'  }
 #'@references Joliffe and Stephenson (2012). Forecast Verification: A 
-#'  Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.
+#'  Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.; 
+#'  Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).
 #'@import multiApply
 #'@importFrom abind abind
 #'@importFrom stats qt qnorm quantile
@@ -146,11 +151,14 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
     stop("Parameter 'exp' and 'obs' must have same dimension names.")
   }
   ## dat_dim
-  if (!is.character(dat_dim) | length(dat_dim) > 1) {
-    stop("Parameter 'dat_dim' must be a character string.")
-  }
-  if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
-    stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
+  if (!is.null(dat_dim)) {
+    if (!is.character(dat_dim) | length(dat_dim) > 1) {
+      stop("Parameter 'dat_dim' must be a character string or NULL.")
+    }
+    if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+      stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.",
+           "Set it as NULL if there is no dataset dimension.")
+    }
   }
   ## space_dim
   if (!is.character(space_dim) | length(space_dim) != 2) {
@@ -178,11 +186,12 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
     }
   }
   ## lat
-  if (!is.null(lat)) {
-    if (!is.numeric(lat) | length(lat) != dim(exp)[space_dim[1]]) {
-      stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
-                  "length as the latitude dimension of 'exp' and 'obs'."))
-    }
+  if (is.null(lat)) {
+    stop("Parameter 'lat' cannot be NULL. It is required for area weighting.")
+  }
+  if (!is.numeric(lat) | length(lat) != dim(exp)[space_dim[1]]) {
+    stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
+                "length as the latitude dimension of 'exp' and 'obs'."))
   }
   ## lon
   if (!is.null(lon)) {
@@ -196,16 +205,22 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
     if (!is.numeric(lonlatbox) | length(lonlatbox) != 4) {
       stop("Parameter 'lonlatbox' must be a numeric vector of 4.")
     }
-  }
-  ## lat, lon, and lonlatbox
-  if (!is.null(lon) & !is.null(lat) & !is.null(lonlatbox)) {
+    if (is.null(lon)) {
+      stop("Parameter 'lat' and 'lon' are required if 'lonlatbox' is specified.")
+    }
     select_lonlat <- TRUE
-  } else if (is.null(lon) & is.null(lat) & is.null(lonlatbox)) {
-    select_lonlat <- FALSE
   } else {
-    stop(paste0("Parameters 'lon', 'lat', and 'lonlatbox' must be used or be ",
-                "NULL at the same time."))
+    select_lonlat <- FALSE
   }
+#  ## lat, lon, and lonlatbox
+#  if (!is.null(lon) & !is.null(lat) & !is.null(lonlatbox)) {
+#    select_lonlat <- TRUE
+#  } else if (is.null(lon) & is.null(lat) & is.null(lonlatbox)) {
+#    select_lonlat <- FALSE
+#  } else {
+#    stop(paste0("Parameters 'lon', 'lat', and 'lonlatbox' must be used or be ",
+#                "NULL at the same time."))
+#  }
   ## conf
   if (!is.logical(conf) | length(conf) > 1) {
     stop("Parameter 'conf' must be one logical value.")
@@ -219,7 +234,8 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
       stop("Parameter 'memb_dim' cannot be NULL when parameter 'conftype' is 'bootstrap'.")
     }
     ## conf.lev
-    if (!is.numeric(conf.lev) | conf.lev < 0 | conf.lev > 1 | length(conf.lev) > 1) {
+    if (!is.numeric(conf.lev) | any(conf.lev < 0) | any(conf.lev > 1) | 
+        length(conf.lev) > 1) {
       stop("Parameter 'conf.lev' must be a numeric number between 0 and 1.")
     }
   }
@@ -229,8 +245,8 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   }
   ## ncores
   if (!is.null(ncores)) {
-    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores < 0 |
-      length(ncores) > 1) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
       stop("Parameter 'ncores' must be a positive integer.")
     }
   }
@@ -293,9 +309,7 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
                                     c(space_dim, dat_dim)),
                  fun = .ACC,
                  dat_dim = dat_dim, avg_dim = avg_dim,
-#------------NEW--------------
                  lat = lat,
-#------------NEW_END--------------
                  conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
                  ncores_input = ncores,
                  ncores = ncores)
@@ -324,9 +338,7 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
                                     c(space_dim, avg_dim, dat_dim)),
                  fun = .ACC,
                  dat_dim = dat_dim, avg_dim = avg_dim,
-#------------NEW--------------
                  lat = lat,
-#------------NEW_END--------------
                  conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
                  ncores_input = ncores,
                  ncores = ncores)
@@ -354,23 +366,27 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
  return(res)
 }
 
-.ACC <- function(exp, obs,  dat_dim = 'dataset', #space_dim = c('lat', 'lon'),
-                 avg_dim = 'sdate', #memb_dim = NULL,
-                 lon = NULL, lat = NULL, lonlatbox = NULL,
+.ACC <- function(exp, obs, dat_dim = 'dataset', avg_dim = 'sdate', lat,
                  conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE,
                  ncores_input = NULL) {
-
-# if (is.null(avg_dim)) 
-  # exp: [space_dim, dat_exp]
-  # obs: [space_dim, dat_obs]
-# if (!is.null(avg_dim)) 
-  # exp: [space_dim, avg_dim, dat_exp]
-  # obs: [space_dim, avg_dim, dat_obs]
-
   # .ACC() should use all the spatial points to calculate ACC. It returns [nexp, nobs].
-
-  nexp <- as.numeric(dim(exp)[length(dim(exp))])
-  nobs <- as.numeric(dim(obs)[length(dim(obs))])
+  # If dat_dim = NULL, it returns a number.
+
+  # if (is.null(avg_dim)) 
+  ## exp: [lat, lon, (dat_exp)]
+  ## obs: [lat, lon, (dat_obs)]
+  # if (!is.null(avg_dim)) 
+  ## exp: [lat, lon, avg_dim, (dat_exp)]
+  ## obs: [lat, lon, avg_dim, (dat_obs)]
+
+  # Add dat_dim temporarily if dat_dim = NULL
+  if (is.null(dat_dim)) {
+    nexp <- 1
+    nobs <- 1
+  } else {
+    nexp <- as.numeric(dim(exp)[length(dim(exp))])
+    nobs <- as.numeric(dim(obs)[length(dim(obs))])
+  }
 
   if (is.null(avg_dim)) {
     acc <- array(dim = c(nexp = nexp, nobs = nobs))
@@ -378,34 +394,67 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
     if (conf) {
       conf.upper <- array(dim = c(nexp = nexp, nobs = nobs))
       conf.lower <- array(dim = c(nexp = nexp, nobs = nobs))
-      if (conftype == 'bootstrap') {
-         ndraw <- 100
-         acc_draw <- array(dim = c(nexp = nexp, nobs = nobs, ndraw)) 
-      }
     }
     
   } else {
-        acc <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
-        names(dim(acc))[3] <- avg_dim
-        macc <- array(dim = c(nexp = nexp, nobs = nobs))
-        if (pval) p.val <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
+    acc <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[avg_dim]))
+    names(dim(acc))[3] <- avg_dim
+    macc <- array(dim = c(nexp = nexp, nobs = nobs))
+    if (pval) p.val <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[avg_dim]))
     if (conf) {
-        conf.upper <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
-        conf.lower <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
-      if (conftype == 'bootstrap') {
-         ndraw <- 100
-         acc_draw <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1], ndraw))
-         macc_draw <- array(dim = c(nexp = nexp, nobs = nobs, ndraw))
-      }
+        conf.upper <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[avg_dim]))
+        conf.lower <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[avg_dim]))
     }
   }  
   
+  # centralize & area weighted
+  ## spatial centralization for each [avg_dim, dat]
+  dim_exp <- dim(exp)
+  dim_obs <- dim(obs)
+  wt <- cos(lat * pi/180)
+  wt <- rep(wt, times = prod(dim_exp[2:length(dim_exp)]))
+
+  if (is.null(avg_dim)) { # [lat, lon, (dat)]
+    if (is.null(dat_dim)) { #[lat, lon]
+      # turn exp and obs into vector, first latitudes and then longitudes
+      exp <- as.vector(exp)
+      obs <- as.vector(obs)
+      exp <- array(sqrt(wt) * (exp - mean(exp, na.rm = TRUE)), dim = dim_exp)
+      obs <- array(sqrt(wt) * (obs - mean(obs, na.rm = TRUE)), dim = dim_obs)
+    } else { #[lat, lon, dat]
+      # exp
+      exp <- array(exp, dim = c(prod(dim_exp[1:2]), dim_exp[3]))
+      mean_exp <- rep(apply(exp, 2, mean, na.rm = TRUE), each = prod(dim_exp[1:2]))
+      exp <- array(sqrt(wt) * (as.vector(exp) - mean_exp), dim = dim_exp)
+      # obs
+      obs <- array(obs, dim = c(prod(dim_obs[1:2]), dim_obs[3]))
+      mean_obs <- rep(apply(obs, 2, mean, na.rm = TRUE), each = prod(dim_obs[1:2]))
+      obs <- array(sqrt(wt) * (as.vector(obs) - mean_obs), dim = dim_obs)
+    }
+  } else {  #[lat, lon, avg_dim, (dat)]
+    # exp
+    exp <- array(exp, dim = c(prod(dim_exp[1:2]), dim_exp[3:length(dim_exp)]))
+    mean_exp <- apply(exp, 2:length(dim(exp)), mean, na.rm = TRUE)  # [avg_dim, (dat)]
+    mean_exp <- rep(as.vector(mean_exp), each = prod(dim_exp[1:2]))
+    exp <- array(sqrt(wt) * (as.vector(exp) - mean_exp), dim = dim_exp)
+    # obs
+    obs <- array(obs, dim = c(prod(dim_obs[1:2]), dim_obs[3:length(dim_obs)]))
+    mean_obs <- apply(obs, 2:length(dim(obs)), mean, na.rm = TRUE)  # [avg_dim, (dat)]
+    mean_obs <- rep(as.vector(mean_obs), each = prod(dim_obs[1:2]))
+    obs <- array(sqrt(wt) * (as.vector(obs) - mean_obs), dim = dim_obs)
+  }
+
   # Per-paired exp and obs. NAs should be in the same position in both exp and obs
   for (iobs in 1:nobs) {
     for (iexp in 1:nexp) {
-      exp_sub <- ClimProjDiags::Subset(exp, dat_dim, iexp, drop = 'selected')
-      obs_sub <- ClimProjDiags::Subset(obs, dat_dim, iobs, drop = 'selected')
-      # dim: [space_dim]
+      if (!is.null(dat_dim)) {
+        exp_sub <- ClimProjDiags::Subset(exp, dat_dim, iexp, drop = 'selected')
+        obs_sub <- ClimProjDiags::Subset(obs, dat_dim, iobs, drop = 'selected')
+      } else {
+        exp_sub <- exp
+        obs_sub <- obs
+      }
+      # dim: [lat, lon, (avg_dim)]
 
       # Variance(iexp) should not take into account any point 
       # that is not available in iobs and therefore not accounted for 
@@ -413,29 +462,26 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
       exp_sub[is.na(obs_sub)] <- NA
       obs_sub[is.na(exp_sub)] <- NA
 
+      # area weights
+      wt <- cos(lat * pi/180)
+
       if (is.null(avg_dim)) {
         # ACC
-#-----------NEW--------------
-        # area weighted
-        wt <- cos(lat * pi/180)
-        wt <- rep(wt, as.numeric(dim(exp_sub)[2])) ## dim(exp): [lat, lon]
-        exp_sub <- as.vector(t(exp_sub)) ## first longitudes and then latitudes
-        obs_sub <- as.vector(t(obs_sub)) ## first longitudes and then latitudes
-        exp_sub <- sqrt(wt) * (exp_sub - mean(exp_sub))
-        obs_sub <- sqrt(wt) * (obs_sub - mean(obs_sub))
-#----------NEW_END-----------
-        top <- sum(exp_sub*obs_sub, na.rm = TRUE)  #a number
+        top <- sum(exp_sub * obs_sub, na.rm = TRUE)  #a number
         bottom <- sqrt(sum(exp_sub^2, na.rm = TRUE) * sum(obs_sub^2, na.rm = TRUE))
         acc[iexp, iobs] <- top/bottom #a number
         # handle bottom = 0
         if (is.infinite(acc[iexp, iobs])) acc[iexp, iobs] <- NA
+
         # pval and conf
         if (pval | conf) {
           if (conftype == "parametric") {
-            # calculate effective sample size along space_dim
-            # combine space_dim into one dim first
-            obs_tmp <- array(obs_sub, dim = c(space = length(obs_sub)))
-            eno <- Eno(obs_tmp, 'space', ncores = ncores_input)  # a number
+            # calculate effective sample size
+#            # combine space_dim into one dim first
+#            obs_tmp <- array(obs_sub, dim = c(space = length(obs_sub)))
+#            eno <- Eno(obs_tmp, 'space', ncores = ncores_input)  # a number
+            eno <- .Eno(as.vector(obs_sub), na.action = na.pass)
+
             if (pval) {
               t <- qt(conf.lev, eno - 2)  # a number
               p.val[iexp, iobs] <- sqrt(t^2 / (t^2 + eno - 2))
@@ -448,18 +494,22 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
         }
 
       } else {  #avg_dim is not NULL
+        # exp_sub and obs_sub: [lat, lon, avg_dim]
+
         # MACC
-        top <- sum(exp_sub*obs_sub, na.rm = TRUE)  #a number
+#        ## area weighted
+        top <- sum(exp_sub * obs_sub, na.rm = TRUE)  #a number
         bottom <- sqrt(sum(exp_sub^2, na.rm = TRUE) * sum(obs_sub^2, na.rm = TRUE))
         macc[iexp, iobs] <- top/bottom #a number
+
         # handle bottom = 0
         if (is.infinite(macc[iexp, iobs])) macc[iexp, iobs] <- NA
+
         # ACC
-        for (i in 1:dim(acc)[3]) {   #NOTE: use sapply!!!
-          exp_sub_i <- ClimProjDiags::Subset(exp_sub, avg_dim, i, drop = 'selected')
-          obs_sub_i <- ClimProjDiags::Subset(obs_sub, avg_dim, i, drop = 'selected')
-          #dim: [space_dim]
-          top <- sum(exp_sub_i*obs_sub_i, na.rm = TRUE)  #a number
+        for (i in 1:dim(acc)[3]) {
+          exp_sub_i <- exp_sub[, , i]
+          obs_sub_i <- obs_sub[, , i]
+          top <- sum(exp_sub_i * obs_sub_i, na.rm = TRUE)  #a number
           bottom <- sqrt(sum(exp_sub_i^2, na.rm = TRUE) * sum(obs_sub_i^2, na.rm = TRUE))
           acc[iexp, iobs, i] <- top/bottom #a number
           # handle bottom = 0
@@ -471,9 +521,10 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
           if (conftype == "parametric") {
             # calculate effective sample size along space_dim 
             # combine space_dim into one dim first
-            obs_tmp <- array(obs_sub, dim = c(space = prod(dim(obs_sub)[-length(dim(obs_sub))]), 
-                                              dim(obs_sub)[length(dim(obs_sub))]))
-            eno <- Eno(obs_tmp, 'space', ncores = ncores_input)  # a vector of avg_dim
+            obs_tmp <- array(obs_sub, 
+                             dim = c(space = prod(dim(obs_sub)[1:2]), 
+                                     dim(obs_sub)[3]))
+            eno <- apply(obs_tmp, 2, .Eno, na.action = na.pass)  # a vector of avg_dim
             if (pval) {
               t <- qt(conf.lev, eno - 2)  # a vector of avg_dim
               p.val[iexp, iobs, ] <- sqrt(t^2 / (t^2 + eno - 2))
@@ -491,8 +542,29 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   }
 
 #------------------------------------------------
-  
-
+  # Remove nexp and nobs if dat_dim = NULL
+  if (is.null(dat_dim)) {
+    if (is.null(avg_dim)) {
+      acc <- as.vector(acc)
+      if (conf) {
+        conf.lower <- as.vector(conf.lower)
+        conf.upper <- as.vector(conf.upper)
+      }   
+      if (pval) {
+        p.val <- as.vector(p.val)
+      }
+    } else {
+      dim(acc) <- dim(acc)[3:length(dim(acc))]
+      macc <- as.vector(macc)
+      if (conf) {
+        dim(conf.lower) <- dim(conf.lower)[3:length(dim(conf.lower))]
+        dim(conf.upper) <- dim(conf.upper)[3:length(dim(conf.upper))]
+      }
+      if (pval) {
+          dim(p.val) <- dim(p.val)[3:length(dim(p.val))]
+      }
+    }
+  }
 
   # Return output
   if (is.null(avg_dim)) {
@@ -524,9 +596,8 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
 }
 
 
-.ACC_bootstrap <- function(exp, obs,  dat_dim = 'dataset', #space_dim = c('lat', 'lon'),
-                           avg_dim = 'sdate', memb_dim = NULL,
-                           lon = NULL, lat = NULL, lonlatbox = NULL,
+.ACC_bootstrap <- function(exp, obs,  dat_dim = 'dataset', 
+                           avg_dim = 'sdate', memb_dim = NULL, lat = NULL,
                            conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE, 
                            ncores_input = NULL) {
 # if (is.null(avg_dim)) 
diff --git a/man/ACC.Rd b/man/ACC.Rd
index d1389fd..a9039d6 100644
--- a/man/ACC.Rd
+++ b/man/ACC.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/ACC.R
 \name{ACC}
 \alias{ACC}
-\title{Compute the anomaly correlation coefficient between the forecast and corresponding observation}
+\title{Compute the spatial anomaly correlation coefficient between the forecast and corresponding observation}
 \usage{
 ACC(
   exp,
@@ -23,37 +23,38 @@ ACC(
 }
 \arguments{
 \item{exp}{A numeric array of experimental anomalies with named dimensions.
-It must have at least 'dat_dim' and 'space_dim'.}
+The dimension must have at least 'space_dim'.}
 
 \item{obs}{A numeric array of observational anomalies with named dimensions.
-It must have the same dimensions as 'exp' except the length of 'dat_dim' 
+The dimension should be the same as 'exp' except the length of 'dat_dim' 
 and 'memb_dim'.}
 
 \item{dat_dim}{A character string indicating the name of dataset (nobs/nexp) 
-dimension. The default value is 'dataset'.}
+dimension. The default value is 'dataset'. If there is no dataset 
+dimension, set NULL.}
 
 \item{space_dim}{A character string vector of 2 indicating the name of the
 latitude and longitude dimensions (in order) along which ACC is computed. 
 The default value is c('lat', 'lon').}
 
 \item{avg_dim}{A character string indicating the name of the dimension to be
-averaged. It must be one of 'time_dim'. The mean ACC is calculated along 
-averaged. If no need to calculate mean ACC, set as NULL. The default value 
-is 'sdate'.}
+averaged, which is usually the time dimension. If no need to calculate mean
+ACC, set as NULL. The default value is 'sdate'.}
 
 \item{memb_dim}{A character string indicating the name of the member 
 dimension. If the data are not ensemble ones, set as NULL. The default 
 value is 'member'.}
 
-\item{lat}{A vector of the latitudes of the exp/obs grids. Only required when
-the domain of interested is specified. The default value is NULL.}
+\item{lat}{A vector of the latitudes of the exp/obs grids. It is used for
+area weighting and when the domain of interested 'lonlatbox' is specified.}
 
 \item{lon}{A vector of the longitudes of the exp/obs grids. Only required when
-the domain of interested is specified. The default value is NULL.}
+the domain of interested 'lonlatbox' is specified. The default value is 
+NULL.}
 
 \item{lonlatbox}{A numeric vector of 4 indicating the corners of the domain of
-interested: c(lonmin, lonmax, latmin, latmax). Only required when the domain
-of interested is specified. The default value is NULL.}
+interested: c(lonmin, lonmax, latmin, latmax). The default value is NULL 
+and the whole data will be used.}
 
 \item{conf}{A logical value indicating whether to retrieve the confidence 
 intervals or not. The default value is TRUE.}
@@ -82,7 +83,8 @@ A list containing the numeric arrays:\cr
 \item{acc}{
  The ACC with the dimensions c(nexp, nobs, the rest of the dimension except 
  space_dim and memb_dim). 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), and nobs is the number of observation (i.e., dat_dim in obs). If 
+ dat_dim is NULL, nexp and nobs are omitted.
 }
 \item{conf.lower (if conftype = "parametric") or acc_conf.lower (if 
      conftype = "bootstrap")}{
@@ -113,13 +115,15 @@ A list containing the numeric arrays:\cr
 }
 }
 \description{
-Calculate the anomaly correlation coefficient for the ensemble mean of each 
-model and the corresponding references over a spatial domain. It can return a
-forecast time series if the data contain forest time dimension, and also the 
-start date mean if the data contain start date dimension. 
-The domain of interest can be specified by providing the list 
-of longitudes/latitudes (lon/lat) of the data together with the corners 
-of the domain: lonlatbox = c(lonmin, lonmax, latmin, latmax).
+Calculate the spatial anomaly correlation coefficient (ACC) for the ensemble
+mean of each model and the corresponding references over a spatial domain. It
+can return a forecast time series if the data contain forest time dimension, 
+and also the ACC mean over one dimension, e.g., start date dimension.
+The domain of interest can be specified by providing the list of longitudes/
+latitudes of the data together with the corners of the domain: lonlatbox = 
+c(lonmin, lonmax, latmin, latmax). The data will be adjusted to have a spatial
+mean of zero, then area weighting is applied. The formula is referenced from 
+Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).
 }
 \examples{
  \dontshow{
@@ -151,5 +155,6 @@ PlotACC(res_bootstrap, startDates)
 }
 \references{
 Joliffe and Stephenson (2012). Forecast Verification: A 
- Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.
+ Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.; 
+ Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).
 }
-- 
GitLab


From 7f0cd2a474c9992774bafe161824b23f9eecd29a Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 16 Dec 2021 15:04:02 +0100
Subject: [PATCH 3/8] Delete unneeded lines and simplified code

---
 R/ACC.R    | 142 ++++++++++++++++++-----------------------------------
 man/ACC.Rd |   4 +-
 2 files changed, 51 insertions(+), 95 deletions(-)

diff --git a/R/ACC.R b/R/ACC.R
index 340dd77..55220d7 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -105,8 +105,8 @@
 #'clim <- Clim(sampleData$mod, sampleData$obs)
 #'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 #'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-#'acc <- ACC(ano_exp, ano_obs)
-#'acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+#'acc <- ACC(ano_exp, ano_obs, lat = sampleData$lat)
+#'acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap', lat = sampleData$lat)
 #'# Combine acc results for PlotACC
 #'res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
 #'             dim = c(dim(acc$acc), 4))
@@ -212,15 +212,6 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   } else {
     select_lonlat <- FALSE
   }
-#  ## lat, lon, and lonlatbox
-#  if (!is.null(lon) & !is.null(lat) & !is.null(lonlatbox)) {
-#    select_lonlat <- TRUE
-#  } else if (is.null(lon) & is.null(lat) & is.null(lonlatbox)) {
-#    select_lonlat <- FALSE
-#  } else {
-#    stop(paste0("Parameters 'lon', 'lat', and 'lonlatbox' must be used or be ",
-#                "NULL at the same time."))
-#  }
   ## conf
   if (!is.logical(conf) | length(conf) > 1) {
     stop("Parameter 'conf' must be one logical value.")
@@ -303,72 +294,50 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
     obs <- MeanDims(obs, memb_dim, na.rm = TRUE)
   }
 
-  if (is.null(avg_dim)) {
-    res <- Apply(list(exp, obs),
-                 target_dims = list(c(space_dim, dat_dim),
-                                    c(space_dim, dat_dim)),
-                 fun = .ACC,
-                 dat_dim = dat_dim, avg_dim = avg_dim,
-                 lat = lat,
-                 conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
-                 ncores_input = ncores,
-                 ncores = ncores)
-
-    if (conftype == 'bootstrap') {
-      res_conf <- Apply(list(exp_ori, obs_ori),
-                 target_dims = list(c(memb_dim, dat_dim, space_dim),
-                                    c(memb_dim, dat_dim, space_dim)),
-                 fun = .ACC_bootstrap,
-                 dat_dim = dat_dim, memb_dim = memb_dim, avg_dim = avg_dim,
-                 conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
-                 ncores_input = ncores,
-                 ncores = ncores)
-      #NOTE: pval?
-      res <- list(acc = res$acc, 
-                  acc_conf.lower = res_conf$acc_conf.lower, 
-                  acc_conf.upper = res_conf$acc_conf.upper,
-                  macc = res$macc,
-                  macc_conf.lower = res_conf$macc_conf.lower, 
-                  macc_conf.upper = res_conf$macc_conf.upper)
-    }
 
+  if (is.null(avg_dim)) {
+    target_dims <- list(c(space_dim, dat_dim), c(space_dim, dat_dim))
   } else {
-    res <- Apply(list(exp, obs),
-                 target_dims = list(c(space_dim, avg_dim, dat_dim),
-                                    c(space_dim, avg_dim, dat_dim)),
+    target_dims <- list(c(space_dim, avg_dim, dat_dim), c(space_dim, avg_dim, dat_dim))
+  }
+  res <- Apply(list(exp, obs),
+                 target_dims = target_dims,
                  fun = .ACC,
                  dat_dim = dat_dim, avg_dim = avg_dim,
                  lat = lat,
                  conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
-                 ncores_input = ncores,
                  ncores = ncores)
 
-    if (conftype == 'bootstrap') {
-      res_conf <- Apply(list(exp_ori, obs_ori),
-                        target_dims = list(c(memb_dim, dat_dim, avg_dim, space_dim),
-                                           c(memb_dim, dat_dim, avg_dim, space_dim)),
-                        fun = .ACC_bootstrap,
-                        dat_dim = dat_dim, memb_dim = memb_dim, avg_dim = avg_dim,
-                        conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
-                        ncores_input = ncores,
-                        ncores = ncores)
-      res <- list(acc = res$acc, 
-                  acc_conf.lower = res_conf$acc_conf.lower, 
-                  acc_conf.upper = res_conf$acc_conf.upper,
-                  macc = res$macc,
-                  macc_conf.lower = res_conf$macc_conf.lower,
-                  macc_conf.upper = res_conf$macc_conf.upper)
-
+  # If bootstrap, calculate confidence level
+  if (conftype == 'bootstrap') {
+    if (is.null(avg_dim)) {
+      target_dims_bs <- list(c(memb_dim, dat_dim, space_dim),
+                             c(memb_dim, dat_dim, space_dim))
+    } else {
+      target_dims_bs <- list(c(memb_dim, dat_dim, avg_dim, space_dim),
+                             c(memb_dim, dat_dim, avg_dim, space_dim))
     }
-
-  }
-
- return(res)
+    res_conf <- Apply(list(exp_ori, obs_ori),
+                      target_dims = target_dims_bs, 
+                      fun = .ACC_bootstrap, 
+                      dat_dim = dat_dim, memb_dim = memb_dim, avg_dim = avg_dim,
+                      lat = lat,
+                      conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
+                      ncores = ncores)
+    #NOTE: pval?
+    res <- list(acc = res$acc,
+                acc_conf.lower = res_conf$acc_conf.lower,
+                acc_conf.upper = res_conf$acc_conf.upper,
+                macc = res$macc,
+                macc_conf.lower = res_conf$macc_conf.lower,
+                macc_conf.upper = res_conf$macc_conf.upper)
+  }
+
+  return(res)
 }
 
 .ACC <- function(exp, obs, dat_dim = 'dataset', avg_dim = 'sdate', lat,
-                 conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE,
-                 ncores_input = NULL) {
+                 conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE) {
   # .ACC() should use all the spatial points to calculate ACC. It returns [nexp, nobs].
   # If dat_dim = NULL, it returns a number.
 
@@ -414,24 +383,13 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   wt <- cos(lat * pi/180)
   wt <- rep(wt, times = prod(dim_exp[2:length(dim_exp)]))
 
-  if (is.null(avg_dim)) { # [lat, lon, (dat)]
-    if (is.null(dat_dim)) { #[lat, lon]
-      # turn exp and obs into vector, first latitudes and then longitudes
-      exp <- as.vector(exp)
-      obs <- as.vector(obs)
-      exp <- array(sqrt(wt) * (exp - mean(exp, na.rm = TRUE)), dim = dim_exp)
-      obs <- array(sqrt(wt) * (obs - mean(obs, na.rm = TRUE)), dim = dim_obs)
-    } else { #[lat, lon, dat]
-      # exp
-      exp <- array(exp, dim = c(prod(dim_exp[1:2]), dim_exp[3]))
-      mean_exp <- rep(apply(exp, 2, mean, na.rm = TRUE), each = prod(dim_exp[1:2]))
-      exp <- array(sqrt(wt) * (as.vector(exp) - mean_exp), dim = dim_exp)
-      # obs
-      obs <- array(obs, dim = c(prod(dim_obs[1:2]), dim_obs[3]))
-      mean_obs <- rep(apply(obs, 2, mean, na.rm = TRUE), each = prod(dim_obs[1:2]))
-      obs <- array(sqrt(wt) * (as.vector(obs) - mean_obs), dim = dim_obs)
-    }
-  } else {  #[lat, lon, avg_dim, (dat)]
+  if (is.null(avg_dim) & is.null(dat_dim)) {  #[lat, lon]
+    # turn exp and obs into vector, first latitudes and then longitudes
+    exp <- as.vector(exp)
+    obs <- as.vector(obs)
+    exp <- array(sqrt(wt) * (exp - mean(exp, na.rm = TRUE)), dim = dim_exp)
+    obs <- array(sqrt(wt) * (obs - mean(obs, na.rm = TRUE)), dim = dim_obs)
+  } else { # [lat, lon, dat], [lat, lon, avg_dim], or [lat, lon, avg_dim, dat]
     # exp
     exp <- array(exp, dim = c(prod(dim_exp[1:2]), dim_exp[3:length(dim_exp)]))
     mean_exp <- apply(exp, 2:length(dim(exp)), mean, na.rm = TRUE)  # [avg_dim, (dat)]
@@ -477,9 +435,6 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
         if (pval | conf) {
           if (conftype == "parametric") {
             # calculate effective sample size
-#            # combine space_dim into one dim first
-#            obs_tmp <- array(obs_sub, dim = c(space = length(obs_sub)))
-#            eno <- Eno(obs_tmp, 'space', ncores = ncores_input)  # a number
             eno <- .Eno(as.vector(obs_sub), na.action = na.pass)
 
             if (pval) {
@@ -497,7 +452,6 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
         # exp_sub and obs_sub: [lat, lon, avg_dim]
 
         # MACC
-#        ## area weighted
         top <- sum(exp_sub * obs_sub, na.rm = TRUE)  #a number
         bottom <- sqrt(sum(exp_sub^2, na.rm = TRUE) * sum(obs_sub^2, na.rm = TRUE))
         macc[iexp, iobs] <- top/bottom #a number
@@ -530,8 +484,10 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
               p.val[iexp, iobs, ] <- sqrt(t^2 / (t^2 + eno - 2))
             }
             if (conf) {
-              conf.upper[iexp, iobs, ] <- tanh(atanh(acc[iexp, iobs, ]) + qnorm(1 - (1 - conf.lev) / 2) / sqrt(eno - 3))
-              conf.lower[iexp, iobs, ] <- tanh(atanh(acc[iexp, iobs, ]) + qnorm((1 - conf.lev) / 2) / sqrt(eno - 3))
+              conf.upper[iexp, iobs, ] <- tanh(atanh(acc[iexp, iobs, ]) + 
+                                          qnorm(1 - (1 - conf.lev) / 2) / sqrt(eno - 3))
+              conf.lower[iexp, iobs, ] <- tanh(atanh(acc[iexp, iobs, ]) + 
+                                          qnorm((1 - conf.lev) / 2) / sqrt(eno - 3))
             }
           }
         }
@@ -597,9 +553,9 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
 
 
 .ACC_bootstrap <- function(exp, obs,  dat_dim = 'dataset', 
-                           avg_dim = 'sdate', memb_dim = NULL, lat = NULL,
-                           conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE, 
-                           ncores_input = NULL) {
+                           avg_dim = 'sdate', memb_dim = NULL, lat,
+                           conf = TRUE, conftype = "parametric", conf.lev = 0.95, 
+                           pval = TRUE) {
 # if (is.null(avg_dim)) 
   # exp: [memb_exp, dat_exp, space_dim]
   # obs: [memb_obs, dat_obs, space_dim]
@@ -657,7 +613,7 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
    
     #calculate the ACC of the randomized field
     tmpACC <- .ACC(drawexp, drawobs, conf = FALSE, pval = FALSE, avg_dim = avg_dim,
-                   ncores_input = ncores_input)
+                   lat = lat)
     if (is.null(avg_dim)) {
       acc_draw[, , jdraw] <- tmpACC$acc
     } else {
diff --git a/man/ACC.Rd b/man/ACC.Rd
index a9039d6..979afcb 100644
--- a/man/ACC.Rd
+++ b/man/ACC.Rd
@@ -141,8 +141,8 @@ sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2)
 clim <- Clim(sampleData$mod, sampleData$obs)
 ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-acc <- ACC(ano_exp, ano_obs)
-acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+acc <- ACC(ano_exp, ano_obs, lat = sampleData$lat)
+acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap', lat = sampleData$lat)
 # Combine acc results for PlotACC
 res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
             dim = c(dim(acc$acc), 4))
-- 
GitLab


From 560ee5c8a3b548c62eab771a630c659805e32c8d Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 16 Dec 2021 15:28:24 +0100
Subject: [PATCH 4/8] Small docfix

---
 R/ACC.R    | 6 ++++--
 man/ACC.Rd | 3 ++-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/R/ACC.R b/R/ACC.R
index 55220d7..ab8185c 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -78,7 +78,8 @@
 #'\item{macc}{
 #'  The mean anomaly correlation coefficient with dimensions
 #'  c(nexp, nobs, the rest of the dimension except space_dim, memb_dim, and 
-#'  avg_dim). Only present if 'avg_dim' is not NULL.
+#'  avg_dim). Only present if 'avg_dim' is not NULL. If dat_dim is NULL, nexp 
+#'  and nobs are omitted.
 #'}
 #'\item{macc_conf.lower}{
 #'  The lower confidence interval of MACC with the same dimensions as MACC. 
@@ -182,7 +183,8 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
       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.")
+      stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.",
+           "Set it as NULL if there is no member dimension.")
     }
   }
   ## lat
diff --git a/man/ACC.Rd b/man/ACC.Rd
index 979afcb..5bd17f6 100644
--- a/man/ACC.Rd
+++ b/man/ACC.Rd
@@ -103,7 +103,8 @@ A list containing the numeric arrays:\cr
 \item{macc}{
  The mean anomaly correlation coefficient with dimensions
  c(nexp, nobs, the rest of the dimension except space_dim, memb_dim, and 
- avg_dim). Only present if 'avg_dim' is not NULL.
+ avg_dim). Only present if 'avg_dim' is not NULL. If dat_dim is NULL, nexp 
+ and nobs are omitted.
 }
 \item{macc_conf.lower}{
  The lower confidence interval of MACC with the same dimensions as MACC. 
-- 
GitLab


From 75cf9835164545d02a2424afcdb665529db3bcb4 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 16 Dec 2021 15:32:32 +0100
Subject: [PATCH 5/8] Remove unneeded line

---
 R/ACC.R | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/R/ACC.R b/R/ACC.R
index ab8185c..113d7ac 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -422,9 +422,6 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
       exp_sub[is.na(obs_sub)] <- NA
       obs_sub[is.na(exp_sub)] <- NA
 
-      # area weights
-      wt <- cos(lat * pi/180)
-
       if (is.null(avg_dim)) {
         # ACC
         top <- sum(exp_sub * obs_sub, na.rm = TRUE)  #a number
-- 
GitLab


From 257c6ec62705749cb1d579fa5cabdf06ce2cfb65 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Mon, 3 Jan 2022 18:23:50 +0100
Subject: [PATCH 6/8] Change space_dim to lat_dim and lon_dim

---
 R/ACC.R                   |  84 ++++++++++++++++-----------
 man/ACC.Rd                |  18 ++++--
 tests/testthat/test-ACC.R | 117 ++++++++++++++++++++++++++------------
 3 files changed, 144 insertions(+), 75 deletions(-)

diff --git a/R/ACC.R b/R/ACC.R
index 113d7ac..ec71d03 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -11,16 +11,19 @@
 #'Wilks (2011; section 7.6.4; https://doi.org/10.1016/B978-0-12-385022-5.00008-7).
 #'
 #'@param exp A numeric array of experimental anomalies with named dimensions.
-#'  The dimension must have at least 'space_dim'.
+#'  The dimension must have at least 'lat_dim' and 'lon_dim'.
 #'@param obs A numeric array of observational anomalies with named dimensions.
 #'  The dimension should be the same as 'exp' except the length of 'dat_dim' 
 #'  and 'memb_dim'.
 #'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
 #'  dimension. The default value is 'dataset'. If there is no dataset 
 #'  dimension, set NULL.
-#'@param space_dim A character string vector of 2 indicating the name of the
-#'  latitude and longitude dimensions (in order) along which ACC is computed. 
-#'  The default value is c('lat', 'lon').
+#'@param lat_dim A character string indicating the name of the latitude
+#'  dimension of 'exp' and 'obs' along which ACC is computed. The default value
+#'  is 'lat'. 
+#'@param lon_dim A character string indicating the name of the longitude
+#'  dimension of 'exp' and 'obs' along which ACC is computed. The default value
+#'  is 'lon'. 
 #'@param avg_dim A character string indicating the name of the dimension to be
 #'  averaged, which is usually the time dimension. If no need to calculate mean
 #'  ACC, set as NULL. The default value is 'sdate'.
@@ -57,7 +60,7 @@
 #'A list containing the numeric arrays:\cr 
 #'\item{acc}{
 #'  The ACC with the dimensions c(nexp, nobs, the rest of the dimension except 
-#'  space_dim and memb_dim). nexp is the number of experiment (i.e., dat_dim in
+#'  lat_dim, lon_dim and memb_dim). 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). If 
 #'  dat_dim is NULL, nexp and nobs are omitted.
 #'}
@@ -77,7 +80,7 @@
 #'}
 #'\item{macc}{
 #'  The mean anomaly correlation coefficient with dimensions
-#'  c(nexp, nobs, the rest of the dimension except space_dim, memb_dim, and 
+#'  c(nexp, nobs, the rest of the dimension except lat_dim, lon_dim, memb_dim, and 
 #'  avg_dim). Only present if 'avg_dim' is not NULL. If dat_dim is NULL, nexp 
 #'  and nobs are omitted.
 #'}
@@ -125,8 +128,8 @@
 #'@importFrom stats qt qnorm quantile
 #'@importFrom ClimProjDiags Subset
 #'@export
-ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
-                avg_dim = 'sdate', memb_dim = 'member', 
+ACC <- function(exp, obs, dat_dim = 'dataset', lat_dim = 'lat', lon_dim = 'lon',
+                space_dim = c('lat', 'lon'), avg_dim = 'sdate', memb_dim = 'member', 
                 lat = NULL, lon = NULL, lonlatbox = NULL, 
                 conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE,
                 ncores = NULL) {
@@ -141,7 +144,7 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   }
   if (is.null(dim(exp)) | is.null(dim(obs))) {
     stop(paste0("Parameter 'exp' and 'obs' must have at least dimensions ",
-                "dat_dim and space_dim."))
+                "lat_dim and lon_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)) {
@@ -161,12 +164,25 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
            "Set it as NULL if there is no dataset dimension.")
     }
   }
-  ## space_dim
-  if (!is.character(space_dim) | length(space_dim) != 2) {
-    stop("Parameter 'space_dim' must be a character vector of 2.")
+  ## space_dim (deprecated)
+  if (!missing("space_dim")) {
+    warning("Parameter 'space_dim' is deprecated. Use 'lat_dim' and 'lon_dim' instead.")
+    lat_dim <- space_dim[1]
+    lon_dim <- space_dim[2]
   }
-  if (any(!space_dim %in% names(dim(exp))) | any(!space_dim %in% names(dim(obs)))) {
-    stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
+  ## lat_dim
+  if (!is.character(lat_dim) | length(lat_dim) != 1) {
+    stop("Parameter 'lat_dim' must be a character string.")
+  }
+  if (!lat_dim %in% names(dim(exp)) | !lat_dim %in% names(dim(obs))) {
+    stop("Parameter 'lat_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## lon_dim
+  if (!is.character(lon_dim) | length(lon_dim) != 1) {
+    stop("Parameter 'lon_dim' must be a character string.")
+  }
+  if (!lon_dim %in% names(dim(exp)) | !lon_dim %in% names(dim(obs))) {
+    stop("Parameter 'lon_dim' is not found in 'exp' or 'obs' dimension.")
   }
   ## avg_dim
   if (!is.null(avg_dim)) {
@@ -191,13 +207,13 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   if (is.null(lat)) {
     stop("Parameter 'lat' cannot be NULL. It is required for area weighting.")
   }
-  if (!is.numeric(lat) | length(lat) != dim(exp)[space_dim[1]]) {
+  if (!is.numeric(lat) | length(lat) != dim(exp)[lat_dim]) {
     stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
                 "length as the latitude dimension of 'exp' and 'obs'."))
   }
   ## lon
   if (!is.null(lon)) {
-    if (!is.numeric(lon) | length(lon) != dim(exp)[space_dim[2]]) {
+    if (!is.numeric(lon) | length(lon) != dim(exp)[lon_dim]) {
       stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
                   "length as the longitude dimension of 'exp' and 'obs'."))
     }
@@ -246,13 +262,15 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   ## exp and obs (2)
   name_exp <- sort(names(dim(exp)))
   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(dat_dim)) {
+    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])) {
+  if (!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
     stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
                 "all the dimensions expect 'dat_dim' and 'memb_dim'."))
   }
@@ -282,8 +300,8 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
                     (lonlatbox[1] > lonlatbox[2] & (lon > lonlatbox[1] | lon < lonlatbox[2])))
     indlat <- which(lat >= lonlatbox[3] & lat <= lonlatbox[4])
 
-    exp <- ClimProjDiags::Subset(exp, space_dim, list(indlat, indlon), drop = FALSE)
-    obs <- ClimProjDiags::Subset(obs, space_dim, list(indlat, indlon), drop = FALSE)
+    exp <- ClimProjDiags::Subset(exp, c(lat_dim, lon_dim), list(indlat, indlon), drop = FALSE)
+    obs <- ClimProjDiags::Subset(obs, c(lat_dim, lon_dim), list(indlat, indlon), drop = FALSE)
   }
 
    # Ensemble mean
@@ -298,9 +316,9 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
 
 
   if (is.null(avg_dim)) {
-    target_dims <- list(c(space_dim, dat_dim), c(space_dim, dat_dim))
+    target_dims <- list(c(lat_dim, lon_dim, dat_dim), c(lat_dim, lon_dim, dat_dim))
   } else {
-    target_dims <- list(c(space_dim, avg_dim, dat_dim), c(space_dim, avg_dim, dat_dim))
+    target_dims <- list(c(lat_dim, lon_dim, avg_dim, dat_dim), c(lat_dim, lon_dim, avg_dim, dat_dim))
   }
   res <- Apply(list(exp, obs),
                  target_dims = target_dims,
@@ -313,11 +331,11 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
   # If bootstrap, calculate confidence level
   if (conftype == 'bootstrap') {
     if (is.null(avg_dim)) {
-      target_dims_bs <- list(c(memb_dim, dat_dim, space_dim),
-                             c(memb_dim, dat_dim, space_dim))
+      target_dims_bs <- list(c(memb_dim, dat_dim, lat_dim, lon_dim),
+                             c(memb_dim, dat_dim, lat_dim, lon_dim))
     } else {
-      target_dims_bs <- list(c(memb_dim, dat_dim, avg_dim, space_dim),
-                             c(memb_dim, dat_dim, avg_dim, space_dim))
+      target_dims_bs <- list(c(memb_dim, dat_dim, avg_dim, lat_dim, lon_dim),
+                             c(memb_dim, dat_dim, avg_dim, lat_dim, lon_dim))
     }
     res_conf <- Apply(list(exp_ori, obs_ori),
                       target_dims = target_dims_bs, 
@@ -472,8 +490,8 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
         # pval and conf
         if (pval | conf) {
           if (conftype == "parametric") {
-            # calculate effective sample size along space_dim 
-            # combine space_dim into one dim first
+            # calculate effective sample size along lat_dim and lon_dim 
+            # combine lat_dim and lon_dim into one dim first
             obs_tmp <- array(obs_sub, 
                              dim = c(space = prod(dim(obs_sub)[1:2]), 
                                      dim(obs_sub)[3]))
@@ -556,11 +574,11 @@ ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
                            conf = TRUE, conftype = "parametric", conf.lev = 0.95, 
                            pval = TRUE) {
 # if (is.null(avg_dim)) 
-  # exp: [memb_exp, dat_exp, space_dim]
-  # obs: [memb_obs, dat_obs, space_dim]
+  # exp: [memb_exp, dat_exp, lat, lon]
+  # obs: [memb_obs, dat_obs, lat, lon]
 # if (!is.null(avg_dim)) 
-  # exp: [memb_exp, dat_exp, avg_dim, space_dim]
-  # obs: [memb_obs, dat_obs, avg_dim, space_dim]
+  # exp: [memb_exp, dat_exp, avg_dim, lat, lon]
+  # obs: [memb_obs, dat_obs, avg_dim, lat, lon]
 
   nexp <- as.numeric(dim(exp)[2])
   nobs <- as.numeric(dim(obs)[2])
diff --git a/man/ACC.Rd b/man/ACC.Rd
index 5bd17f6..1402808 100644
--- a/man/ACC.Rd
+++ b/man/ACC.Rd
@@ -8,6 +8,8 @@ ACC(
   exp,
   obs,
   dat_dim = "dataset",
+  lat_dim = "lat",
+  lon_dim = "lon",
   space_dim = c("lat", "lon"),
   avg_dim = "sdate",
   memb_dim = "member",
@@ -23,7 +25,7 @@ ACC(
 }
 \arguments{
 \item{exp}{A numeric array of experimental anomalies with named dimensions.
-The dimension must have at least 'space_dim'.}
+The dimension must have at least 'lat_dim' and 'lon_dim'.}
 
 \item{obs}{A numeric array of observational anomalies with named dimensions.
 The dimension should be the same as 'exp' except the length of 'dat_dim' 
@@ -33,9 +35,13 @@ and 'memb_dim'.}
 dimension. The default value is 'dataset'. If there is no dataset 
 dimension, set NULL.}
 
-\item{space_dim}{A character string vector of 2 indicating the name of the
-latitude and longitude dimensions (in order) along which ACC is computed. 
-The default value is c('lat', 'lon').}
+\item{lat_dim}{A character string indicating the name of the latitude
+dimension of 'exp' and 'obs' along which ACC is computed. The default value
+is 'lat'.}
+
+\item{lon_dim}{A character string indicating the name of the longitude
+dimension of 'exp' and 'obs' along which ACC is computed. The default value
+is 'lon'.}
 
 \item{avg_dim}{A character string indicating the name of the dimension to be
 averaged, which is usually the time dimension. If no need to calculate mean
@@ -82,7 +88,7 @@ computation. The default value is NULL.}
 A list containing the numeric arrays:\cr 
 \item{acc}{
  The ACC with the dimensions c(nexp, nobs, the rest of the dimension except 
- space_dim and memb_dim). nexp is the number of experiment (i.e., dat_dim in
+ lat_dim, lon_dim and memb_dim). 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). If 
  dat_dim is NULL, nexp and nobs are omitted.
 }
@@ -102,7 +108,7 @@ A list containing the numeric arrays:\cr
 }
 \item{macc}{
  The mean anomaly correlation coefficient with dimensions
- c(nexp, nobs, the rest of the dimension except space_dim, memb_dim, and 
+ c(nexp, nobs, the rest of the dimension except lat_dim, lon_dim, memb_dim, and 
  avg_dim). Only present if 'avg_dim' is not NULL. If dat_dim is NULL, nexp 
  and nobs are omitted.
 }
diff --git a/tests/testthat/test-ACC.R b/tests/testthat/test-ACC.R
index 7631401..5c36725 100644
--- a/tests/testthat/test-ACC.R
+++ b/tests/testthat/test-ACC.R
@@ -11,6 +11,8 @@ context("s2dv::ACC tests")
   set.seed(2)
   obs1 <- array(rnorm(30), dim = c(dataset = 1, member = 1, sdate = 5,
                                    ftime = 1, lat = 2, lon = 3))
+  lat1 <- c(30, 35)
+  lon1 <- c(0, 5, 10)
   # dat2
   set.seed(1)
   exp2 <- array(rnorm(60), dim = c(dataset = 2, sdate = 5,
@@ -21,10 +23,12 @@ context("s2dv::ACC tests")
   set.seed(2)
   na <- floor(runif(2, min = 1, max = 30))
   obs2[na] <- NA
+  lat2 <- c(30, 35)
+  lon2 <- c(0, 5, 10)
 
 ##############################################
 test_that("1. Input checks", {
-
+  # exp and obs (1)
   expect_error(
   ACC(c(), c()),
   "Parameter 'exp' and 'obs' cannot be NULL."
@@ -36,7 +40,7 @@ test_that("1. Input checks", {
   expect_error(
   ACC(c(1:10), c(2:4)),
   paste0("Parameter 'exp' and 'obs' must have at least dimensions ",
-         "dat_dim and space_dim.")
+         "lat_dim and lon_dim.")
   )
   expect_error(
   ACC(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
@@ -46,6 +50,7 @@ test_that("1. Input checks", {
   ACC(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"
   )
+  # dat_dim
   expect_error(
   ACC(exp1, obs1, dat_dim = 1),
   "Parameter 'dat_dim' must be a character string."
@@ -54,14 +59,30 @@ test_that("1. Input checks", {
   ACC(exp1, obs1, dat_dim = 'a'),
   "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
   )
+  # space_dim (deprecated)
+  expect_warning(
+  ACC(exp1, obs1, space_dim = c('lat', 'lon'), lat = c(1, 2)),
+  "Parameter 'space_dim' is deprecated. Use 'lat_dim' and 'lon_dim' instead."
+  )
+  # lat_dim
+  expect_error(
+  ACC(exp1, obs1, lat_dim = 1),
+  "Parameter 'lat_dim' must be a character string."
+  )
+  expect_error(
+  ACC(exp1, obs1, lat_dim = 'latitude'),
+  "Parameter 'lat_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # lon_dim
   expect_error(
-  ACC(exp1, obs1, space_dim = c('lon')),
-  "Parameter 'space_dim' must be a character vector of 2."
+  ACC(exp1, obs1, lon_dim = 2),
+  "Parameter 'lon_dim' must be a character string."
   )
   expect_error(
-  ACC(exp1, obs1, space_dim = c('lon', 'lev')),
-  "Parameter 'space_dim' is not found in 'exp' or 'obs' dimension."
+  ACC(exp1, obs1, lon_dim = 'lonn'),
+  "Parameter 'lon_dim' is not found in 'exp' or 'obs' dimension."
   )
+  # avg_dim
   expect_error(
   ACC(exp1, obs1, avg_dim = 1),
   "Parameter 'avg_dim' must be a character string."
@@ -70,6 +91,7 @@ test_that("1. Input checks", {
   ACC(exp1, obs1, avg_dim = c('lev')),
   "Parameter 'avg_dim' is not found in 'exp' or 'obs' dimension."
   )
+  # memb_dim
   expect_error(
   ACC(exp1, obs1, memb_dim = TRUE),
   "Parameter 'memb_dim' must be a character string."
@@ -78,52 +100,57 @@ test_that("1. Input checks", {
   ACC(exp1, obs1, memb_dim = 'memb'),
   "Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension."
   )
+  # lat
   expect_error(
   ACC(exp1, obs1, lat = c(1, 2, 3)),
   paste0("Parameter \'lat\' must be a numeric vector with the same ",
          "length as the latitude dimension of \'exp\' and \'obs\'.")
   )
+  # lon
   expect_error(
-  ACC(exp1, obs1, lon = c(1, 3)),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 3)),
   paste0("Parameter \'lon\' must be a numeric vector with the same ",
          "length as the longitude dimension of \'exp\' and \'obs\'.")
   )
+  # lonlatbox
   expect_error(
-  ACC(exp1, obs1, lonlatbox = c(-90, 90)),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), lonlatbox = c(-90, 90)),
   "Parameter 'lonlatbox' must be a numeric vector of 4."
   )
+  # conf
   expect_error(
-  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3)),
-  paste0("Parameters 'lon', 'lat', and 'lonlatbox' must be used or be ",
-         "NULL at the same time.")
-  )
-  expect_error(
-  ACC(exp1, obs1, conf = 1),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), conf = 1),
   "Parameter 'conf' must be one logical value."
   )
+  # conftype
   expect_error(
-  ACC(exp1, obs1, conftype = 'a'),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), conftype = 'a'),
   "Parameter \'conftype\' must be either \'parametric\' or \'bootstrap\'."
   )
   expect_error(
-  ACC(exp1, obs1, memb_dim = NULL, conftype = 'bootstrap'),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), memb_dim = NULL, conftype = 'bootstrap'),
   "Parameter 'memb_dim' cannot be NULL when parameter 'conftype' is 'bootstrap'."
   )
+  # conf.lev
   expect_error(
-  ACC(exp1, obs1, conf.lev = -1),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), conf.lev = -1),
   "Parameter 'conf.lev' must be a numeric number between 0 and 1."
   )
+  # pval
   expect_error(
-  ACC(exp1, obs1, pval = 'TRUE'),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), pval = 'TRUE'),
   "Parameter 'pval' must be one logical value."
   )
+  # ncores
   expect_error(
-  ACC(exp1, obs1, ncores = 1.5),
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3), ncores = 1.5),
   "Parameter 'ncores' must be a positive integer."
   )
+  # exp and obs (2)
   expect_error(
   ACC(exp = array(1:10, dim = c(dataset = 5, member = 1, lat = 2, lon = 1)),
       obs = array(1:4, dim = c(dataset = 2, member = 2, lat = 1, lon = 1)),
+      lat = c(1, 2), lon = c(1),
       avg_dim = NULL),
   "Parameter 'exp' and 'obs' must have same length of all the dimensions expect 'dat_dim' and 'memb_dim'."
   )
@@ -134,66 +161,84 @@ test_that("1. Input checks", {
 test_that("2. Output checks: dat1", {
 
   expect_equal(
-    dim(ACC(exp1, obs1)$acc),
+    dim(ACC(exp1, obs1, lat = lat1, lon = lon1)$acc),
     c(nexp = 1, nobs = 1, sdate = 5, ftime = 1)
   )
   expect_equal(
-    names(ACC(exp1, obs1)),
+    names(ACC(exp1, obs1, lat = lat1)),
     c("acc", "conf.lower", "conf.upper", "p.val", "macc")
   )
   expect_equal(
-    mean(ACC(exp1, obs1)$acc),
-    -0.2001352,
+    as.vector(ACC(exp1, obs1, lat = lat1)$acc),
+    c(-0.006384684, -0.664580451, -0.637292180, 0.320974249, -0.453246995),
     tolerance = 0.00001
   )
   expect_equal(
-    as.vector(ACC(exp1, obs1)$p.val),   
+    as.vector(ACC(exp1, obs1, lat = lat1)$p.val),   
     c(0.7292993, 0.7292993, 0.7292993, 0.7292993, 0.7292993),
     tolerance = 0.00001    
   )
   expect_equal(
-    as.vector(ACC(exp1, obs1)$conf.lower),
-    c(-0.8595534, -0.9644555, -0.9408508, -0.6887500, -0.7619374),
+    as.vector(ACC(exp1, obs1, lat = lat1)$conf.lower),
+    c(-0.8137296, -0.9589397, -0.9549511, -0.6633950, -0.9246771),
     tolerance = 0.00001
   )
   expect_equal(
-    as.vector(ACC(exp1, obs1)$conf.upper),
-    c(0.7493799, 0.2515608, 0.4759707, 0.8890967, 0.8517117),
+    as.vector(ACC(exp1, obs1, lat = lat1)$conf.upper),
+    c(0.8093704, 0.3190711, 0.3609563, 0.8984881, 0.5668074),
     tolerance = 0.00001
   )
   expect_equal(
-    names(ACC(exp1, obs1, avg_dim = NULL)),
+    names(ACC(exp1, obs1, lat = lat1, avg_dim = NULL)),
     c("acc", "conf.lower", "conf.upper", "p.val")
   )
   expect_equal(
-    dim(ACC(exp1, obs1, dat_dim = 'member', memb_dim = NULL)$acc),
+    dim(ACC(exp1, obs1, lat = lat1, dat_dim = 'member', memb_dim = NULL)$acc),
     c(nexp = 2, nobs = 1, sdate = 5, dataset = 1, ftime = 1)
   )
   expect_equal(
-    names(ACC(exp1, obs1, conf = FALSE)),
+    names(ACC(exp1, obs1, lat = lat1, conf = FALSE)),
     c("acc", "p.val", "macc")
   )
   expect_equal(
-    names(ACC(exp1, obs1, pval = FALSE)),
+    names(ACC(exp1, obs1, lat = lat1, pval = FALSE)),
     c("acc", "conf.lower", "conf.upper", "macc")
   )
   expect_equal(
-    names(ACC(exp1, obs1, conf = FALSE, pval = FALSE)),
+    names(ACC(exp1, obs1, lat = lat1, conf = FALSE, pval = FALSE)),
     c("acc", "macc")
   )
   expect_equal(
-    as.vector(ACC(exp1, obs1, conf = FALSE, avg_dim = NULL, conf.lev = 0.9)$p.val),
+    as.vector(ACC(exp1, obs1, lat = lat1, conf = FALSE, avg_dim = NULL, conf.lev = 0.9)$p.val),
     c(0.6083998, 0.6083998, 0.6083998, 0.6083998, 0.6083998),
     tolerance = 0.00001
   )
   expect_equal(
-    mean(ACC(exp1, obs1, lat = c(10, 20), lon = c(20, 30, 40), lonlatbox = c(20, 30, 10, 20))$acc),
-    -0.1681097,
+    as.vector(ACC(exp1, obs1, lat = c(10, 20), lon = c(20, 30, 40), lonlatbox = c(20, 30, 10, 20))$acc),
+    c(0.3555372, -0.5205672, -0.3403301, 0.2690648, -0.5146910),
     tolerance = 0.00001
   )
+
 })
 
 
+##############################################
 
+test_that("3. Output checks: dat2", {
 
+expect_equal(
+  dim(ACC(exp2, obs2, lat = lat2, lon = lon2, memb_dim = NULL)$acc),
+  c(nexp = 2, nobs = 1, sdate = 5, ftime = 1)
+)
+expect_equal(
+  as.vector(ACC(exp2, obs2, lat = lat2, memb_dim = NULL)$acc)[3:7],
+  c(-0.3601880, -0.5624773, -0.4603762, -0.6997169, -0.1336961),
+  tolerance = 0.00001
+)
+expect_equal(
+  mean(ACC(exp2, obs2, lat = lat2, memb_dim = NULL)$acc),
+  -0.1484762,
+  tolerance = 0.00001
+)
 
+})
-- 
GitLab


From 3af3504e5d11c63b7da1291c86fc53f9bc5225f3 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Mon, 3 Jan 2022 18:32:18 +0100
Subject: [PATCH 7/8] Revise example due to new development of ACC()

---
 R/PlotACC.R    | 4 ++--
 man/PlotACC.Rd | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/PlotACC.R b/R/PlotACC.R
index 3bffa68..a674ff6 100644
--- a/R/PlotACC.R
+++ b/R/PlotACC.R
@@ -65,8 +65,8 @@
 #'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 #'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
 
-#'acc <- ACC(ano_exp, ano_obs)
-#'acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+#'acc <- ACC(ano_exp, ano_obs, lat = sampleData$lat)
+#'acc_bootstrap <- ACC(ano_exp, ano_obs, lat = sampleData$lat, conftype = 'bootstrap')
 #'# Combine acc results for PlotACC
 #'res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
 #'             dim = c(dim(acc$acc), 4))
diff --git a/man/PlotACC.Rd b/man/PlotACC.Rd
index cd2b357..2764de3 100644
--- a/man/PlotACC.Rd
+++ b/man/PlotACC.Rd
@@ -110,8 +110,8 @@ sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2)
 clim <- Clim(sampleData$mod, sampleData$obs)
 ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-acc <- ACC(ano_exp, ano_obs)
-acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+acc <- ACC(ano_exp, ano_obs, lat = sampleData$lat)
+acc_bootstrap <- ACC(ano_exp, ano_obs, lat = sampleData$lat, conftype = 'bootstrap')
 # Combine acc results for PlotACC
 res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
             dim = c(dim(acc$acc), 4))
-- 
GitLab


From 3507a8d7e9310e5759e9d7d00fafdaa38e520c4f Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Thu, 20 Jan 2022 09:44:58 +0100
Subject: [PATCH 8/8] Add 'space_dim' back to documentation

---
 R/ACC.R    | 4 ++++
 man/ACC.Rd | 5 +++++
 2 files changed, 9 insertions(+)

diff --git a/R/ACC.R b/R/ACC.R
index ec71d03..44c724c 100644
--- a/R/ACC.R
+++ b/R/ACC.R
@@ -24,6 +24,10 @@
 #'@param lon_dim A character string indicating the name of the longitude
 #'  dimension of 'exp' and 'obs' along which ACC is computed. The default value
 #'  is 'lon'. 
+#'@param space_dim A character string vector of 2 indicating the name of the
+#'  latitude and longitude dimensions (in order) along which ACC is computed. 
+#'  The default value is c('lat', 'lon'). This argument has been deprecated.
+#'  Use 'lat_dim' and 'lon_dim' instead.
 #'@param avg_dim A character string indicating the name of the dimension to be
 #'  averaged, which is usually the time dimension. If no need to calculate mean
 #'  ACC, set as NULL. The default value is 'sdate'.
diff --git a/man/ACC.Rd b/man/ACC.Rd
index 1402808..7df6abb 100644
--- a/man/ACC.Rd
+++ b/man/ACC.Rd
@@ -43,6 +43,11 @@ is 'lat'.}
 dimension of 'exp' and 'obs' along which ACC is computed. The default value
 is 'lon'.}
 
+\item{space_dim}{A character string vector of 2 indicating the name of the
+latitude and longitude dimensions (in order) along which ACC is computed. 
+The default value is c('lat', 'lon'). This argument has been deprecated.
+Use 'lat_dim' and 'lon_dim' instead.}
+
 \item{avg_dim}{A character string indicating the name of the dimension to be
 averaged, which is usually the time dimension. If no need to calculate mean
 ACC, set as NULL. The default value is 'sdate'.}
-- 
GitLab