From 6c3bda23ef59db869e2c1404d39c194972a799fb Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 30 Aug 2023 11:45:29 +0200
Subject: [PATCH 01/12] First attemp forecast NAO

---
 R/NAO.R          | 245 +++++++++++++++++++++++++++++++++++------------
 R/ProjectField.R |   1 -
 2 files changed, 182 insertions(+), 64 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index fb5220c..60286f4 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -11,7 +11,7 @@
 #'(obs) based on the leading EOF pattern.
 #'
 #'@param exp A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
-#'  forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+#'  hindcast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
 #'  dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
 #'  If only NAO of observational data needs to be computed, this parameter can
 #'  be left to NULL. The default value is NULL.
@@ -20,6 +20,10 @@
 #'  dimensions 'time_dim', 'ftime_dim', and 'space_dim' at least.
 #'  If only NAO of experimental data needs to be computed, this parameter can 
 #'  be left to NULL. The default value is NULL.
+#'@param exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
+#'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()] with 
+#'  If only NAO of observational data needs to be computed, this parameter can
+#'  be left to NULL. The default value is NULL.
 #'@param lat A vector of the latitudes of 'exp' and 'obs'.
 #'@param lon A vector of the longitudes of 'exp' and 'obs'.
 #'@param time_dim A character string indicating the name of the time dimension
@@ -73,6 +77,8 @@
 #'lon <- seq(-80, 40, length.out = 9) 
 #'nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 #'
+#'exp_cor <- array(rnorm(2*5*6*9), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
+#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon)
 #'# plot the NAO index
 #'  \dontrun{
 #'nao$exp <- Reorder(nao$exp, c(2, 1))
@@ -84,7 +90,7 @@
 #'@import multiApply
 #'@importFrom ClimProjDiags Subset
 #'@export
-NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
+NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sdate',
                 memb_dim = 'member', space_dim = c('lat', 'lon'),
                 ftime_dim = 'ftime', ftime_avg = 2:4, 
                 obsproj = TRUE, ncores = NULL) {
@@ -118,6 +124,18 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'obs' must have dimension names.")
     }
   }
+  if (!is.null(exp_cor)) {
+    if (!is.numeric(exp)) {
+      stop("Parameter 'exp' must be a numeric array.")
+    }
+    if (is.null(dim(exp_cor))) {
+      stop(paste0("Parameter 'exp' must have at least dimensions ",
+                  "time_dim, memb_dim, space_dim, and ftime_dim."))
+    }
+    if(any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
+      stop("Parameter 'exp' must have dimension names.")
+    }
+  }
   ## time_dim
   if (!is.character(time_dim) | length(time_dim) > 1) {
     stop("Parameter 'time_dim' must be a character string.")
@@ -153,6 +171,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       add_member_back <- FALSE
     }
   }
+  if (!is.null(exp_cor)) {
+    if (!memb_dim %in% names(dim(exp_cor))) {
+      stop("Parameter 'memb_dim' is not found in 'exp_cor' dimension.")
+    }
+  }
   ## space_dim
   if (!is.character(space_dim) | length(space_dim) != 2) {
     stop("Parameter 'space_dim' must be a character vector of 2.")
@@ -167,6 +190,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
     }
   }
+  if (!is.null(exp_cor)) {
+    if (any(!space_dim %in% names(dim(exp_cor)))) {
+      stop("Parameter 'space_dim' is not found in 'exp_cor' dimensions.")
+    }
+  }
   ## ftime_dim
   if (!is.character(ftime_dim) | length(ftime_dim) > 1) {
     stop("Parameter 'ftime_dim' must be a character string.")
@@ -181,6 +209,11 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop("Parameter 'ftime_dim' is not found in 'exp' or 'obs' dimension.")
     }
   }
+  if (!is.null(exp_cor)) {
+    if (!ftime_dim %in% names(dim(exp_cor))) {
+      stop("Parameter 'ftime_dim' is not found in 'exp_cor' dimensions.")
+    }
+  }
   ## exp and obs (2)
   if (!is.null(exp) & !is.null(obs)) {
     name_exp <- sort(names(dim(exp)))
@@ -208,11 +241,17 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       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 (!is.null(obs)) {
       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_cor)) {
+      if (max(ftime_avg) > dim(exp_cor)[ftime_dim] | min(ftime_avg) < 1) {
+        stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
+      }
+    }
   }
   ## sdate >= 2
   if (!is.null(exp)) {
@@ -234,7 +273,8 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
                   "length as the longitude dimension of 'exp' and 'obs'."))
     }
-  } else {
+  } 
+  if (!is.null(obs)) {
     if (!is.numeric(lat) | length(lat) != dim(obs)[space_dim[1]]) {
       stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
                   "length as the latitude dimension of 'exp' and 'obs'."))
@@ -244,6 +284,16 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
                   "length as the longitude dimension of 'exp' and 'obs'."))
     }
   } 
+  if (!is.null(exp_cor)) {
+    if (!is.numeric(lat) | length(lat) != dim(exp_cor)[space_dim[1]]) {
+      stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
+                  "length as the latitude dimension of 'exp_cor'."))
+    }
+    if (!is.numeric(lon) | length(lon) != dim(exp_cor)[space_dim[2]]) {
+      stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
+                  "length as the longitude dimension of 'exp_cor'."))
+    }
+  }
   stop_needed <- FALSE
   if (max(lat) > 80 | min(lat) < 20) {
     stop_needed <- TRUE
@@ -303,19 +353,34 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
       obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE)
       obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE)
     }
+    if (!is.null(exp_cor)) {
+      exp_cor_sub <- ClimProjDiags::Subset(exp_cor, ftime_dim, ftime_avg, drop = FALSE)
+      exp_cor <- MeanDims(exp_cor_sub, ftime_dim, na.rm = TRUE)
+    }
   }
 
   # wght
   wght <- array(sqrt(cos(lat * pi/180)), dim = c(length(lat), length(lon)))
 
   if (!is.null(exp) & !is.null(obs)) {
-    res <- Apply(list(exp, obs),
-                 target_dims = list(exp = c(memb_dim, time_dim, space_dim),
-                                    obs = c(time_dim, space_dim)),
-                 fun = .NAO,
-                 lat = lat, wght = wght, 
-                 obsproj = obsproj, add_member_back = add_member_back,
-                 ncores = ncores)
+    if (is.null(exp_cor)) {
+      res <- Apply(list(exp, obs),
+                   target_dims = list(exp = c(memb_dim, time_dim, space_dim),
+                                      obs = c(time_dim, space_dim)),
+                   fun = .NAO,
+                   lat = lat, wght = wght, 
+                   obsproj = obsproj, add_member_back = add_member_back,
+                   ncores = ncores)
+    } else {
+      res <- Apply(list(exp, obs, exp_cor),
+                   target_dims = list(exp = c(memb_dim, time_dim, space_dim),
+                                      obs = c(time_dim, space_dim),
+                                      exp_cor = c(memb_dim, time_dim, space_dim)),
+                   fun = .NAO,
+                   lat = lat, wght = wght,
+                   obsproj = obsproj, add_member_back = add_member_back,
+                   ncores = ncores)
+    }
   } else if (!is.null(exp)) {
     res <- Apply(list(exp = exp),
                  target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
@@ -340,11 +405,12 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
   return(res)
 }
 
-.NAO <- function(exp = NULL, obs = NULL, lat, wght, obsproj = TRUE, add_member_back = FALSE) {
+.NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL,
+                 lat, wght, obsproj = TRUE, add_member_back = FALSE) {
   # exp: [memb_exp, sdate, lat, lon]
   # obs: [sdate, lat, lon]
+  # exp_cor: [memb, sdate=1, lat, lon]
   # wght: [lat, lon]
-
   if (!is.null(exp)) {
     ntime <- dim(exp)[2]
     nlat <- dim(exp)[3]
@@ -358,73 +424,126 @@ NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
 
   if (!is.null(obs)) NAOO.ver <- array(NA, dim = ntime)
   if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
+  if (!is.null(exp_cor)) NAOF.cor <- array(NA, dim = dim(exp_cor)[1])
 
-  for (tt in 1:ntime) {  #sdate
+  if (is.null(exp_cor)) { # cross-validation:
+    for (tt in 1:ntime) {  #sdate
 
-    if (!is.null(obs)) {
-      ## Calculate observation EOF. Excluding one forecast start year.
-      obs_sub <- obs[c(1:ntime)[-tt], , , drop = FALSE]
-      obs_EOF <- .EOF(obs_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+      if (!is.null(obs)) {
+        ## Calculate observation EOF. Excluding one forecast start year.
+        obs_sub <- obs[c(1:ntime)[-tt], , , drop = FALSE]
+        obs_EOF <- .EOF(obs_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+        ## Correct polarity of pattern.
+        # dim(obs_EOF$EOFs): [mode, lat, lon]
+        if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+          obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
+#          obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
+        }
+        ## Project observed anomalies.
+        PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+        ## Keep PCs of excluded forecast start year. Fabian.
+        NAOO.ver[tt] <- PF[tt]
+      }
+
+      if (!is.null(exp)) {
+        if (!obsproj) {
+          exp_sub <- exp[, c(1:ntime)[-tt], , , drop = FALSE]
+          # Combine 'memb' and 'sdate' to calculate EOF
+          dim(exp_sub) <- c(nmemb_exp * (ntime - 1), nlat, nlon)
+          exp_EOF <- .EOF(exp_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+
+          ## Correct polarity of pattern.
+          ##NOTE: different from s2dverification, which doesn't use mean().
+#          if (0 < exp_EOF$EOFs[1, which.min(abs(lat - 65)), ]) {
+          if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+            exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
+#            exp_EOF$PCs <- exp_EOF$PCs * sign  # not used
+          }
+
+          ### Lines below could be simplified further by computing
+          ### ProjectField() only on the year of interest... (though this is
+          ### not vital). Lauriane
+          for (imemb in 1:nmemb_exp) {
+            PF <- .ProjectField(exp[imemb, , , ], eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate, memb]
+            NAOF.ver[tt, imemb] <- PF[tt]
+          }    
+        } else {
+          ## Project forecast anomalies on obs EOF
+          for (imemb in 1:nmemb_exp) {
+            PF <- .ProjectField(exp[imemb, , , ], eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+            NAOF.ver[tt, imemb] <- PF[tt]
+          }
+        }
+      }
+
+    }  # for loop sdate
 
+    # add_member_back
+    if (add_member_back) {
+      suppressWarnings( 
+      NAOO.ver <- InsertDim(NAOO.ver, 2, 1, name = names(dim(exp))[1])
+      )
+    }
+  } else { # exp_cor provided
+    if (!is.null(obs) & obsproj) {
+      obs_EOF <- .EOF(obs, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
       ## Correct polarity of pattern.
       # dim(obs_EOF$EOFs): [mode, lat, lon]
       if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
         obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
-#        obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
+#       obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
       }
-      ## Project observed anomalies.
       PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-      ## Keep PCs of excluded forecast start year. Fabian.
-      NAOO.ver[tt] <- PF[tt]
-    }
-
-    if (!is.null(exp)) {
-      if (!obsproj) {
-        exp_sub <- exp[, c(1:ntime)[-tt], , , drop = FALSE]
-        # Combine 'memb' and 'sdate' to calculate EOF
-        dim(exp_sub) <- c(nmemb_exp * (ntime - 1), nlat, nlon)
-        exp_EOF <- .EOF(exp_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+      NAOO.ver <- PF
 
-        ## Correct polarity of pattern.
-        ##NOTE: different from s2dverification, which doesn't use mean().
-#        if (0 < exp_EOF$EOFs[1, which.min(abs(lat - 65)), ]) {
-        if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-          exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
-#          exp_EOF$PCs <- exp_EOF$PCs * sign  # not used
-        }
+      ## Project observed anomalies.
+      for (imemb in 1:dim(exp_cor)[1]) {
+        exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb,
+                                         drop = 'selected')
 
-        ### Lines below could be simplified further by computing
-        ### ProjectField() only on the year of interest... (though this is
-        ### not vital). Lauriane
-        for (imemb in 1:nmemb_exp) {
-          PF <- .ProjectField(exp[imemb, , , ], eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate, memb]
-          NAOF.ver[tt, imemb] <- PF[tt]
-        }    
-      } else {
-        ## Project forecast anomalies on obs EOF
-        for (imemb in 1:nmemb_exp) {
-          PF <- .ProjectField(exp[imemb, , , ], eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-          NAOF.ver[tt, imemb] <- PF[tt]
-        }
+        PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+        NAOF.cor[imemb] <- PF
       }
-    }
+    } else if (!is.null(exp) & !obsproj) {
+      dim(exp) <- c(nmemb_exp*ntime, nlat, nlon)
+      exp_EOF <- .EOF(exp, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+      ## Correct polarity of pattern.
+      # dim(obs_EOF$EOFs): [mode, lat, lon]
+      if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+        exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
+#       obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
+      }
+      PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+      NAOF.ver <- PF
 
-  }  # for loop sdate
+      ## Project observed anomalies.
+      for (imemb in 1:dim(exp_cor)[1]) {
+        exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb,
+                                         drop = 'selected')
 
-  # add_member_back
-  if (add_member_back) {
-    suppressWarnings( 
-    NAOO.ver <- InsertDim(NAOO.ver, 2, 1, name = names(dim(exp))[1])
-    )
+        PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+        NAOF.cor[imemb] <- PF
+      }
+    }
   }
 
   #NOTE: EOFs_obs is not returned because it's only the result of the last sdate
   #     (It is returned in s2dverification.)
-  if (!is.null(exp) & !is.null(obs)) {
-    return(list(exp = NAOF.ver, obs = NAOO.ver)) #, EOFs_obs = obs_EOF))
-  } else if (!is.null(exp)) {
-    return(list(exp = NAOF.ver))
-  } else if (!is.null(obs)) {
-    return(list(obs = NAOO.ver))
+  if (is.null(exp_cor)) {
+    if (!is.null(exp) & !is.null(obs)) {
+      return(list(exp = NAOF.ver, obs = NAOO.ver)) #, EOFs_obs = obs_EOF))
+    } else if (!is.null(exp)) {
+      return(list(exp = NAOF.ver))
+    } else if (!is.null(obs)) {
+      return(list(obs = NAOO.ver))
+    }
+  } else {
+    if (!is.null(exp) & !obsproj) {
+      return(list(exp = NAOF.ver, exp_cor = NAOF.cor))
+    } else if (!is.null(obs) & obsproj) {
+      return(list(obs = NAOO.ver, exp_cor = NAOF.cor))
+    } else {
+      stop("Unconsidered case. Please review")
+    }
   }
 }
diff --git a/R/ProjectField.R b/R/ProjectField.R
index 55e7fd2..7e73253 100644
--- a/R/ProjectField.R
+++ b/R/ProjectField.R
@@ -229,7 +229,6 @@ ProjectField <- function(ano, eof, time_dim = 'sdate', space_dim = c('lat', 'lon
   # wght: [lat, lon]
 
   ntime <- dim(ano)[1]
-
   if (length(dim(eof_mode)) == 2) {  # mode != NULL
     # Initialization of pc.ver.
     pc.ver <- array(NA, dim = ntime)  #[sdate]
-- 
GitLab


From ccac96527b80568ac9faf219ccbf55476fadc512 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 30 Aug 2023 15:53:01 +0200
Subject: [PATCH 02/12] exp_cor checks correctly placed

---
 R/NAO.R | 83 ++++++++++++++++++++++++++++++++-------------------------
 1 file changed, 46 insertions(+), 37 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index 60286f4..24f8ae3 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -125,15 +125,15 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     }
   }
   if (!is.null(exp_cor)) {
-    if (!is.numeric(exp)) {
-      stop("Parameter 'exp' must be a numeric array.")
+    if (!is.numeric(exp_cor)) {
+      stop("Parameter 'exp_cor' must be a numeric array.")
     }
     if (is.null(dim(exp_cor))) {
-      stop(paste0("Parameter 'exp' must have at least dimensions ",
+      stop(paste0("Parameter 'exp_cor' must have at least dimensions ",
                   "time_dim, memb_dim, space_dim, and ftime_dim."))
     }
     if(any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
-      stop("Parameter 'exp' must have dimension names.")
+      stop("Parameter 'exp_cor' must have dimension names.")
     }
   }
   ## time_dim
@@ -329,8 +329,8 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     if (is.null(obs)) {
       stop("Parameter 'obsproj' set to TRUE but no 'obs' provided.")
     }
-    if (is.null(exp)) {
-      .warning("parameter 'obsproj' set to TRUE but no 'exp' provided.")
+    if (is.null(exp) & is.null(exp_cor)) {
+      .warning("parameter 'obsproj' set to TRUE but no 'exp' nor 'exp_cor' provided.")
     }
   }
   ## ncores
@@ -361,46 +361,56 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
 
   # wght
   wght <- array(sqrt(cos(lat * pi/180)), dim = c(length(lat), length(lon)))
-
-  if (!is.null(exp) & !is.null(obs)) {
-    if (is.null(exp_cor)) {
-      res <- Apply(list(exp, obs),
-                   target_dims = list(exp = c(memb_dim, time_dim, space_dim),
-                                      obs = c(time_dim, space_dim)),
+  if (is.null(exp_cor)) {
+    if (!is.null(exp) & !is.null(obs)) {
+        res <- Apply(list(exp, obs),
+                     target_dims = list(exp = c(memb_dim, time_dim, space_dim),
+                                        obs = c(time_dim, space_dim)),
+                     fun = .NAO,
+                     lat = lat, wght = wght, 
+                     obsproj = obsproj, add_member_back = add_member_back,
+                     ncores = ncores)
+    } else if (!is.null(exp)) {
+      res <- Apply(list(exp = exp),
+                   target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
+                   fun = .NAO,
+                   lat = lat, wght = wght, obs = NULL, 
+                   obsproj = obsproj, add_member_back = FALSE,
+                   ncores = ncores)
+    } else if (!is.null(obs)) {
+      if (add_member_back) {
+        output_dims <- list(obs = c(time_dim, memb_dim))
+      } else {
+        output_dims <- list(obs = time_dim)
+      }
+      res <- Apply(list(obs = obs),
+                   target_dims = list(obs = c(time_dim, space_dim)),
+                   output_dims = output_dims,
                    fun = .NAO,
-                   lat = lat, wght = wght, 
+                   lat = lat, wght = wght, exp = NULL, 
                    obsproj = obsproj, add_member_back = add_member_back,
                    ncores = ncores)
-    } else {
-      res <- Apply(list(exp, obs, exp_cor),
+    }
+  } else { # exp_cor provided
+    if (!is.null(exp) & !obsproj) {
+      res <- Apply(list(exp = exp, exp_cor = exp_cor),
                    target_dims = list(exp = c(memb_dim, time_dim, space_dim),
-                                      obs = c(time_dim, space_dim),
                                       exp_cor = c(memb_dim, time_dim, space_dim)),
                    fun = .NAO,
+                   obs = NULL,
+                   lat = lat, wght = wght,
+                   obsproj = obsproj, add_member_back = add_member_back,
+                   ncores = ncores)
+    } else if (!is.null(obs) & obsproj) {
+      res <- Apply(list(obs = obs, exp_cor = exp_cor),
+                   target_dims = list(obs = c(time_dim, space_dim),
+                                      exp_cor = c(memb_dim, time_dim, space_dim)),
+                   fun = .NAO,
+                   exp = NULL,
                    lat = lat, wght = wght,
                    obsproj = obsproj, add_member_back = add_member_back,
                    ncores = ncores)
     }
-  } else if (!is.null(exp)) {
-    res <- Apply(list(exp = exp),
-                 target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
-                 fun = .NAO,
-                 lat = lat, wght = wght, obs = NULL, 
-                 obsproj = obsproj, add_member_back = FALSE,
-                 ncores = ncores)
-  } else if (!is.null(obs)) {
-    if (add_member_back) {
-      output_dims <- list(obs = c(time_dim, memb_dim))
-    } else {
-      output_dims <- list(obs = time_dim)
-    }
-    res <- Apply(list(obs = obs),
-                 target_dims = list(obs = c(time_dim, space_dim)),
-                 output_dims = output_dims,
-                 fun = .NAO,
-                 lat = lat, wght = wght, exp = NULL, 
-                 obsproj = obsproj, add_member_back = add_member_back,
-                 ncores = ncores)
   }
   return(res)
 }
@@ -515,7 +525,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       }
       PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
       NAOF.ver <- PF
-
       ## Project observed anomalies.
       for (imemb in 1:dim(exp_cor)[1]) {
         exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb,
-- 
GitLab


From e31a37a7928d870e92fb1e44e5a72c7a26434a87 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 30 Aug 2023 16:26:34 +0200
Subject: [PATCH 03/12] Fix for only exp provided

---
 R/NAO.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/NAO.R b/R/NAO.R
index 24f8ae3..16fd4e8 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -399,7 +399,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
                    fun = .NAO,
                    obs = NULL,
                    lat = lat, wght = wght,
-                   obsproj = obsproj, add_member_back = add_member_back,
+                   obsproj = obsproj, #add_member_back = add_member_back,
                    ncores = ncores)
     } else if (!is.null(obs) & obsproj) {
       res <- Apply(list(obs = obs, exp_cor = exp_cor),
-- 
GitLab


From 66e6a2bae45f8433bf639d133619791c7a460f80 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Fri, 8 Sep 2023 13:05:59 +0200
Subject: [PATCH 04/12] Trivial adjustment and fix unit test

---
 R/NAO.R                   |  6 +++---
 man/NAO.Rd                | 10 +++++++++-
 tests/testthat/test-NAO.R |  2 +-
 3 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index 16fd4e8..a34584a 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -21,7 +21,7 @@
 #'  If only NAO of experimental data needs to be computed, this parameter can 
 #'  be left to NULL. The default value is NULL.
 #'@param exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
-#'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()] with 
+#'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with 
 #'  If only NAO of observational data needs to be computed, this parameter can
 #'  be left to NULL. The default value is NULL.
 #'@param lat A vector of the latitudes of 'exp' and 'obs'.
@@ -533,6 +533,8 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
         PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
         NAOF.cor[imemb] <- PF
       }
+    } else {
+      stop("Unconsidered case. Please review.")
     }
   }
 
@@ -551,8 +553,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       return(list(exp = NAOF.ver, exp_cor = NAOF.cor))
     } else if (!is.null(obs) & obsproj) {
       return(list(obs = NAOO.ver, exp_cor = NAOF.cor))
-    } else {
-      stop("Unconsidered case. Please review")
     }
   }
 }
diff --git a/man/NAO.Rd b/man/NAO.Rd
index 999fd75..e605469 100644
--- a/man/NAO.Rd
+++ b/man/NAO.Rd
@@ -7,6 +7,7 @@
 NAO(
   exp = NULL,
   obs = NULL,
+  exp_cor = NULL,
   lat,
   lon,
   time_dim = "sdate",
@@ -20,7 +21,7 @@ NAO(
 }
 \arguments{
 \item{exp}{A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
-forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+hindcast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
 dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
 If only NAO of observational data needs to be computed, this parameter can
 be left to NULL. The default value is NULL.}
@@ -31,6 +32,11 @@ dimensions 'time_dim', 'ftime_dim', and 'space_dim' at least.
 If only NAO of experimental data needs to be computed, this parameter can 
 be left to NULL. The default value is NULL.}
 
+\item{exp_cor}{A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
+forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with 
+If only NAO of observational data needs to be computed, this parameter can
+be left to NULL. The default value is NULL.}
+
 \item{lat}{A vector of the latitudes of 'exp' and 'obs'.}
 
 \item{lon}{A vector of the longitudes of 'exp' and 'obs'.}
@@ -97,6 +103,8 @@ lat <- seq(20, 80, length.out = 6)
 lon <- seq(-80, 40, length.out = 9) 
 nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 
+exp_cor <- array(rnorm(2*5*6*9), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
+nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon)
 # plot the NAO index
  \dontrun{
 nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R
index 91b9943..bbc7832 100644
--- a/tests/testthat/test-NAO.R
+++ b/tests/testthat/test-NAO.R
@@ -153,7 +153,7 @@ test_that("1. Input checks", {
   )
   # ncores
   expect_error(
-    NAO(exp1, obs1, lat1, lon1, ncore = 3.5),
+    NAO(exp1, obs1, lat = lat1, lon = lon1, ncore = 3.5),
     "Parameter 'ncores' must be a positive integer."
   )
 
-- 
GitLab


From 59260099aa72c4ac5d6d25314281043d99c8b36f Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 20 Sep 2023 11:57:29 +0200
Subject: [PATCH 05/12] Add missing sentence in doc

---
 R/NAO.R                   | 3 ++-
 tests/testthat/test-NAO.R | 8 +++++++-
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index a34584a..0d9f9a5 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -21,7 +21,8 @@
 #'  If only NAO of experimental data needs to be computed, this parameter can 
 #'  be left to NULL. The default value is NULL.
 #'@param exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
-#'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with 
+#'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with
+#'  dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least. 
 #'  If only NAO of observational data needs to be computed, this parameter can
 #'  be left to NULL. The default value is NULL.
 #'@param lat A vector of the latitudes of 'exp' and 'obs'.
diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R
index bbc7832..f8e2b1f 100644
--- a/tests/testthat/test-NAO.R
+++ b/tests/testthat/test-NAO.R
@@ -14,7 +14,8 @@
   obs2 <- array(rnorm(108), dim = c(sdate = 3, ftime = 4, lat = 3, lon = 3))
   lat2 <- c(80, 50, 20)
   lon2 <- c(-80, 0, 40)
-
+  set.seed(3)
+  exp2_cor <- array(rnorm(72), dim = c(sdate = 1, ftime = 4, member = 2, lat = 3, lon = 3))
  
 ##############################################
 test_that("1. Input checks", {
@@ -253,6 +254,11 @@ test_that("3. dat2", {
   tolerance = 0.00001
   )
 
+  # exp_cor
+#  expect_equal(
+#  dim(NAO(exp2, obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp),
+#  c(sdate = 3, member = 2)
+#  )
 })
 
 ##############################################
-- 
GitLab


From 9db4e45b4b48eb6f0833f0fff3521b5157e957d0 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 20 Sep 2023 17:39:18 +0200
Subject: [PATCH 06/12] update docu and checks

---
 R/NAO.R | 20 +++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index a34584a..de2454e 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -22,8 +22,9 @@
 #'  be left to NULL. The default value is NULL.
 #'@param exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
 #'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with 
-#'  If only NAO of observational data needs to be computed, this parameter can
-#'  be left to NULL. The default value is NULL.
+#'  dimension 'time_dim' of length 1 (as in the case of an operational forecast).
+#'  If only NAO of reference period needs to be computed, this parameter can
+#'  be left to NULL. The default value is NULL. 
 #'@param lat A vector of the latitudes of 'exp' and 'obs'.
 #'@param lon A vector of the longitudes of 'exp' and 'obs'.
 #'@param time_dim A character string indicating the name of the time dimension
@@ -44,7 +45,8 @@
 #'  reference (TRUE) or compute the NAO by first computing the leading 
 #'  EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
 #'  year you are evaluating out), and then projecting forecast anomalies onto 
-#'  this EOF (FALSE). The default value is TRUE.
+#'  this EOF (FALSE). The default value is TRUE. If 'exp_cor' is provided, 'obs'
+#'  will be used when obsproj is TRUE and 'exp' will be used when  obsproj is FALSE.
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
@@ -135,6 +137,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     if(any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
       stop("Parameter 'exp_cor' must have dimension names.")
     }
+    if (dim(exp_cor)[time_dim] > 1) {
+      stop(paste("Parameter 'exp_cor' is expected to have length 1 on the", 
+                  time_dim, "dimension."))
+    }
   }
   ## time_dim
   if (!is.character(time_dim) | length(time_dim) > 1) {
@@ -393,6 +399,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     }
   } else { # exp_cor provided
     if (!is.null(exp) & !obsproj) {
+      if (!is.null(obs)) {
+        warning("Reference data provided in parameter 'obs' is not used when,
+                parameter 'obsproj' is set to FALSE and 'exp' is provided.")
+      }
       res <- Apply(list(exp = exp, exp_cor = exp_cor),
                    target_dims = list(exp = c(memb_dim, time_dim, space_dim),
                                       exp_cor = c(memb_dim, time_dim, space_dim)),
@@ -402,6 +412,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
                    obsproj = obsproj, #add_member_back = add_member_back,
                    ncores = ncores)
     } else if (!is.null(obs) & obsproj) {
+      if (!is.null(exp)) {
+        warning("Experimental data provided in parameter 'exp' is not used when,
+                 parameter 'obsproj' is set to TRUE and 'obs' is provided.")
+      }
       res <- Apply(list(obs = obs, exp_cor = exp_cor),
                    target_dims = list(obs = c(time_dim, space_dim),
                                       exp_cor = c(memb_dim, time_dim, space_dim)),
-- 
GitLab


From 972d06e5f1d821450b637102e8a610086f2eab15 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 21 Sep 2023 14:47:29 +0200
Subject: [PATCH 07/12] fix memb_dim

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

diff --git a/R/NAO.R b/R/NAO.R
index de2454e..9dd0854 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -21,7 +21,7 @@
 #'  If only NAO of experimental data needs to be computed, this parameter can 
 #'  be left to NULL. The default value is NULL.
 #'@param exp_cor A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
-#'  forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with 
+#'  forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
 #'  dimension 'time_dim' of length 1 (as in the case of an operational forecast).
 #'  If only NAO of reference period needs to be computed, this parameter can
 #'  be left to NULL. The default value is NULL. 
@@ -522,7 +522,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
 
       ## Project observed anomalies.
       for (imemb in 1:dim(exp_cor)[1]) {
-        exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb,
+        exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
                                          drop = 'selected')
 
         PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
@@ -541,7 +541,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       NAOF.ver <- PF
       ## Project observed anomalies.
       for (imemb in 1:dim(exp_cor)[1]) {
-        exp_sub <- ClimProjDiags::Subset(exp_cor, along = memb_dim, indices = imemb,
+        exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
                                          drop = 'selected')
 
         PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
diff --git a/man/NAO.Rd b/man/NAO.Rd
index e605469..84f445e 100644
--- a/man/NAO.Rd
+++ b/man/NAO.Rd
@@ -33,8 +33,9 @@ If only NAO of experimental data needs to be computed, this parameter can
 be left to NULL. The default value is NULL.}
 
 \item{exp_cor}{A named numeric array of the Nort Atlantic SLP (20-80N, 80W-40E)
-forecast anomalies from \cod{Ano()} or \code{Ano_CrossValid()} with 
-If only NAO of observational data needs to be computed, this parameter can
+forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
+dimension 'time_dim' of length 1 (as in the case of an operational forecast).
+If only NAO of reference period needs to be computed, this parameter can
 be left to NULL. The default value is NULL.}
 
 \item{lat}{A vector of the latitudes of 'exp' and 'obs'.}
@@ -64,7 +65,8 @@ projecting the forecast anomalies onto the leading EOF of observational
 reference (TRUE) or compute the NAO by first computing the leading 
 EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
 year you are evaluating out), and then projecting forecast anomalies onto 
-this EOF (FALSE). The default value is TRUE.}
+this EOF (FALSE). The default value is TRUE. If 'exp_cor' is provided, 'obs'
+will be used when obsproj is TRUE and 'exp' will be used when  obsproj is FALSE.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
-- 
GitLab


From 874f6d25d04f54bef0c74eaad42c8d9321ce5c93 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Fri, 22 Sep 2023 16:20:33 +0200
Subject: [PATCH 08/12] Correct exp output dim when exp and exp_cor are used

---
 R/NAO.R                   | 12 +++++++++---
 tests/testthat/test-NAO.R | 32 ++++++++++++++++++++++++++++----
 2 files changed, 37 insertions(+), 7 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index 40fc533..0f3cb7f 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -470,8 +470,11 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
 
   if (!is.null(obs)) NAOO.ver <- array(NA, dim = ntime)
   if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
-  if (!is.null(exp_cor)) NAOF.cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
-
+  if (!is.null(exp_cor)) {
+    NAOF.cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
+    #NOTE: The dimensions are flipped to fill in data correctly. Need to flip it back later.
+    if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(nmemb_exp, ntime))
+  }
   if (is.null(exp_cor)) { # cross-validation:
     for (tt in 1:ntime) {  #sdate
 
@@ -561,7 +564,10 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
 #       obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
       }
       PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-      NAOF.ver <- PF
+      NAOF.ver[, ] <- PF
+      # Flip the dimensions back
+      NAOF.ver <- aperm(NAOF.ver, 2:1)
+
       ## Project observed anomalies.
       for (imemb in 1:dim(exp_cor)[1]) {
         exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R
index 5f1eafa..7b64cac 100644
--- a/tests/testthat/test-NAO.R
+++ b/tests/testthat/test-NAO.R
@@ -263,10 +263,34 @@ test_that("3. dat2", {
   dim(NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs),
   c(sdate = 3)
   )
-#res <- NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2, ftime_avg = NULL)
-#WRONG
-#res1 <- NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)
-#res2 <- NAO(exp = exp2, obs = obs2, exp_cor = NULL, lat = lat2, lon = lon2)
+  expect_equal(
+  c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)),
+  c(0.2121340, 0.1634516),
+  tolerance = 0.0001)
+  expect_equal(
+  c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)),
+  c(0.3511294, -0.7196260, -1.5123894),
+  tolerance = 0.0001
+  )
+
+  expect_equal(
+  dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor),
+  c(sdate = 1, member = 2)
+  )
+  expect_equal(
+  dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp),
+  c(sdate = 3, member = 2)
+  )
+  expect_equal(
+  c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp),
+  c(-0.4225802, -1.3616918, 0.2808729, 0.2723077, -0.7584804, -0.7264514),
+  tolerance = 0.0001
+  )
+  expect_equal(
+  c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor),
+  c(-0.22949531, -0.06946422),
+  tolerance = 0.0001
+  )
 
 })
 
-- 
GitLab


From d00a030a9db883a9fb65cbc6c994bfe43111b4d6 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Tue, 24 Oct 2023 17:40:43 +0200
Subject: [PATCH 09/12] NAO fcst and hcst returned

---
 R/NAO.R | 88 +++++++++++++++++++++++----------------------------------
 1 file changed, 36 insertions(+), 52 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index 0f3cb7f..33adc47 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -88,7 +88,7 @@
 #'nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 #'
 #'exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
-#'nao <- NAO(exp = NULL, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon)
+#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T)
 #'# plot the NAO index
 #'  \dontrun{
 #'nao$exp <- Reorder(nao$exp, c(2, 1))
@@ -104,7 +104,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
                 memb_dim = 'member', space_dim = c('lat', 'lon'),
                 ftime_dim = 'ftime', ftime_avg = 2:4, 
                 obsproj = TRUE, ncores = NULL) {
-
   # Check inputs 
   ## exp, obs, and exp_cor (1)
   if (is.null(obs) & is.null(exp)) {
@@ -411,36 +410,17 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
                    ncores = ncores)
     }
   } else { # exp_cor provided
-    if (!is.null(exp) & !obsproj) {
-      if (!is.null(obs)) {
-        .warning("Reference data provided in parameter 'obs' is not used when ",
-                 "parameter 'obsproj' is set to FALSE and 'exp' is provided.")
-      }
-      res <- Apply(list(exp = exp, exp_cor = exp_cor),
+    if (!is.null(exp) & !is.null(obs) & !is.null(exp_cor)) {
+      res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor),
                    target_dims = list(exp = c(memb_dim, time_dim, space_dim),
+                                      obs = c(time_dim, space_dim),
                                       exp_cor = c(memb_dim, time_dim, space_dim)),
                    fun = .NAO,
-                   obs = NULL,
                    lat = lat, wght = wght,
                    obsproj = obsproj, add_member_back = add_member_back,
                    ncores = ncores)
-
-    } else if (!is.null(obs) & obsproj) {
-      if (!is.null(exp)) {
-        .warning("Experimental data provided in parameter 'exp' is not used when ",
-                 "parameter 'obsproj' is set to TRUE and 'obs' is provided.")
-      }
-      res <- Apply(list(obs = obs, exp_cor = exp_cor),
-                   target_dims = list(obs = c(time_dim, space_dim),
-                                      exp_cor = c(memb_dim, time_dim, space_dim)),
-                   fun = .NAO,
-                   exp = NULL,
-                   lat = lat, wght = wght,
-                   obsproj = obsproj, add_member_back = add_member_back,
-                   ncores = ncores)
-    } else {
-      stop("Unconsidered case. Please review inputs 'exp', 'obs', 'exp_cor', ",
-           "'obsproj' or open an issue on GitLab.")
+    } else { 
+      stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not null.")
     }
   }
 
@@ -467,13 +447,12 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     ntime_exp_cor <- dim(exp_cor)[2] # should be 1
     nmemb_exp_cor <- dim(exp_cor)[1]
   }
-
   if (!is.null(obs)) NAOO.ver <- array(NA, dim = ntime)
   if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
   if (!is.null(exp_cor)) {
     NAOF.cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
     #NOTE: The dimensions are flipped to fill in data correctly. Need to flip it back later.
-    if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(nmemb_exp, ntime))
+  #  if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(nmemb_exp, ntime))
   }
   if (is.null(exp_cor)) { # cross-validation:
     for (tt in 1:ntime) {  #sdate
@@ -535,26 +514,32 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     }
 
   } else { # exp_cor provided
-    if (!is.null(obs) & obsproj) {
-      obs_EOF <- .EOF(obs, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
-      ## Correct polarity of pattern.
-      # dim(obs_EOF$EOFs): [mode, lat, lon]
-      if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-        obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
-#       obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
-      }
-      PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-      NAOO.ver[] <- PF
-
-      ## Project observed anomalies.
+    ## Calculate observation EOF. Without cross-validation
+    obs_EOF <- .EOF(obs, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+    ## Correct polarity of pattern.
+    # dim(obs_EOF$EOFs): [mode, lat, lon]
+    if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+      obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
+      #          obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
+    }
+    ## Project observed anomalies.
+    PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs, wght = wght)  # [sdate]
+    NAOO.ver[] <- PF[1,]
+    if (obsproj) {
+      ## Project forecasts anomalies 
       for (imemb in 1:dim(exp_cor)[1]) {
         exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
                                          drop = 'selected')
-
         PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-        NAOF.cor[imemb] <- PF
+        NAOF.cor[,imemb] <- PF
       }
-    } else if (!is.null(exp) & !obsproj) {
+      for (imemb in 1:dim(exp)[1]) {
+        exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb,
+                                         drop = 'selected')
+        PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)
+        NAOF.ver[,imemb] <- PF
+      }
+    } else if (!obsproj) {
       dim(exp) <- c(nmemb_exp*ntime, nlat, nlon)
       exp_EOF <- .EOF(exp, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
       ## Correct polarity of pattern.
@@ -563,12 +548,16 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
         exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
 #       obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
       }
-      PF <- .ProjectField(exp, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-      NAOF.ver[, ] <- PF
+      dim(exp) <- c(nmemb_exp, ntime, nlat, nlon)
+      for (imemb in 1:dim(exp)[1]) {
+        exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, drop = 'selected')
+        PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+        NAOF.ver[ ,imemb] <- PF
+      }
       # Flip the dimensions back
       NAOF.ver <- aperm(NAOF.ver, 2:1)
 
-      ## Project observed anomalies.
+      ## Project hindcast anomalies on forecast
       for (imemb in 1:dim(exp_cor)[1]) {
         exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
                                          drop = 'selected')
@@ -578,7 +567,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       }
     }
   }
-
   #NOTE: EOFs_obs is not returned because it's only the result of the last sdate
   #     (It is returned in s2dverification.)
   if (is.null(exp_cor)) {
@@ -590,10 +578,6 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       return(list(obs = NAOO.ver))
     }
   } else {
-    if (!is.null(exp) & !obsproj) {
-      return(list(exp = NAOF.ver, exp_cor = NAOF.cor))
-    } else if (!is.null(obs) & obsproj) {
-      return(list(obs = NAOO.ver, exp_cor = NAOF.cor))
-    }
+    return(list(exp = NAOF.ver, obs = NAOO.ver, exp_cor = NAOF.cor))
   }
 }
-- 
GitLab


From 600907e9aae46875cf6eb37d7432203e357be2a5 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 15 Nov 2023 14:31:29 +0100
Subject: [PATCH 10/12] Tidy codes, fix unit tests

---
 R/NAO.R                   | 214 ++++++++++++++++++--------------------
 man/NAO.Rd                |  16 +--
 tests/testthat/test-NAO.R |  73 +++++++++----
 3 files changed, 161 insertions(+), 142 deletions(-)

diff --git a/R/NAO.R b/R/NAO.R
index 33adc47..2a0e118 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -9,7 +9,7 @@
 #'lead seasonal forecasts that can be plotted with PlotBoxWhisker(). It returns
 #'cross-validated PCs of the NAO index for hindcast (exp) and observations 
 #'(obs) based on the leading EOF pattern, or, if forecast (exp_cor) is provided,
-#'the NAO index for forecast and the chosen corresponding data (exp or obs).
+#'the NAO index for forecast and the corresponding data (exp and obs).
 #'
 #'@param exp A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
 #'  hindcast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
@@ -44,12 +44,12 @@
 #'  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 
-#'  EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
-#'  year you are evaluating out), and then projecting forecast anomalies onto 
-#'  this EOF (FALSE). The default value is TRUE. If 'exp_cor' is provided, 'obs'
-#'  will be used when obsproj is TRUE and 'exp' will be used when obsproj is 
-#'  FALSE.
+#'  reference (TRUE, default) or compute the NAO by first computing the leading
+#'  EOF of the forecast anomalies (in cross-validation mode, i.e. leave the 
+#'  evaluated year out), then projecting forecast anomalies onto this EOF
+#'  (FALSE). If 'exp_cor' is provided, 'obs' will be used when obsproj is TRUE
+#'  and 'exp' will be used when obsproj is FALSE, and no cross-validation is
+#'  applied.
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
@@ -141,9 +141,12 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       stop(paste0("Parameter 'exp_cor' must have at least dimensions ",
                   "time_dim, memb_dim, space_dim, and ftime_dim."))
     }
-    if(any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
+    if (any(is.null(names(dim(exp_cor)))) | any(nchar(names(dim(exp_cor))) == 0)) {
       stop("Parameter 'exp_cor' must have dimension names.")
     }
+    if (is.null(exp) || is.null(obs)) {
+      stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not provided.")
+    }
   }
   ## time_dim
   if (!is.character(time_dim) | length(time_dim) > 1) {
@@ -164,8 +167,8 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       stop("Parameter 'time_dim' is not found in 'exp_cor' dimension.")
     }
     if (dim(exp_cor)[time_dim] > 1) {
-      stop(paste("Parameter 'exp_cor' is expected to have length 1 on the ",
-                 time_dim, "dimension."))
+      stop("Parameter 'exp_cor' is expected to have length 1 in ",
+           time_dim, "dimension.")
     }
   }
 
@@ -233,6 +236,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     }
   }
   ## exp and obs (2)
+  #TODO: Add checks for exp_cor
   if (!is.null(exp) & !is.null(obs)) {
     name_exp <- sort(names(dim(exp)))
     name_obs <- sort(names(dim(obs)))
@@ -242,12 +246,12 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
       throw_error <- TRUE
     } else if (any(name_exp != name_obs)) {
       throw_error <- TRUE
-    } else if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
+    } else if (!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
       throw_error <- TRUE
     }
     if (throw_error) {
-      stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
-                  "of all the dimensions except 'memb_dim'."))
+      stop("Parameter 'exp' and 'obs' must have the same names and lengths ",
+           "of all the dimensions except 'memb_dim'.")
     }
   }
   ## ftime_avg
@@ -410,18 +414,14 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
                    ncores = ncores)
     }
   } else { # exp_cor provided
-    if (!is.null(exp) & !is.null(obs) & !is.null(exp_cor)) {
-      res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor),
-                   target_dims = list(exp = c(memb_dim, time_dim, space_dim),
-                                      obs = c(time_dim, space_dim),
-                                      exp_cor = c(memb_dim, time_dim, space_dim)),
-                   fun = .NAO,
-                   lat = lat, wght = wght,
-                   obsproj = obsproj, add_member_back = add_member_back,
-                   ncores = ncores)
-    } else { 
-      stop("Parameters 'exp' and 'obs' are required when 'exp_cor' is not null.")
-    }
+    res <- Apply(list(exp = exp, obs = obs, exp_cor = exp_cor),
+                 target_dims = list(exp = c(memb_dim, time_dim, space_dim),
+                                    obs = c(time_dim, space_dim),
+                                    exp_cor = c(memb_dim, time_dim, space_dim)),
+                 fun = .NAO,
+                 lat = lat, wght = wght,
+                 obsproj = obsproj, add_member_back = add_member_back,
+                 ncores = ncores)
   }
 
   return(res)
@@ -433,6 +433,7 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
   # obs: [sdate, lat, lon]
   # exp_cor: [memb, sdate = 1, lat, lon]
   # wght: [lat, lon]
+
   if (!is.null(exp)) {
     ntime <- dim(exp)[2]
     nlat <- dim(exp)[3]
@@ -447,30 +448,31 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
     ntime_exp_cor <- dim(exp_cor)[2] # should be 1
     nmemb_exp_cor <- dim(exp_cor)[1]
   }
-  if (!is.null(obs)) NAOO.ver <- array(NA, dim = ntime)
-  if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))
+
+  if (!is.null(obs)) nao_obs <- array(NA, dim = ntime)
+  if (!is.null(exp)) nao_exp <- array(NA, dim = c(ntime, nmemb_exp))
   if (!is.null(exp_cor)) {
-    NAOF.cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
+    nao_exp_cor <- array(NA, dim = c(ntime_exp_cor, nmemb_exp_cor))
     #NOTE: The dimensions are flipped to fill in data correctly. Need to flip it back later.
-  #  if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(nmemb_exp, ntime))
   }
-  if (is.null(exp_cor)) { # cross-validation:
-    for (tt in 1:ntime) {  #sdate
+
+  if (is.null(exp_cor)) { 
+
+    for (tt in 1:ntime) {  # cross-validation
 
       if (!is.null(obs)) {
         ## Calculate observation EOF. Excluding one forecast start year.
         obs_sub <- obs[c(1:ntime)[-tt], , , drop = FALSE]
-        obs_EOF <- .EOF(obs_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
-        ## Correct polarity of pattern.
-        # dim(obs_EOF$EOFs): [mode, lat, lon]
-        if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-          obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
-#          obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
+        EOF_obs <- .EOF(obs_sub, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
+        ## Correct polarity of pattern
+        # EOF_obs: [mode = 1, lat, lon]
+        if (0 < mean(EOF_obs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+          EOF_obs <- EOF_obs * (-1)
         }
         ## Project observed anomalies.
-        PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
+        PF <- .ProjectField(obs, eof_mode = EOF_obs[1, , ], wght = wght)  # [sdate]
         ## Keep PCs of excluded forecast start year. Fabian.
-        NAOO.ver[tt] <- PF[tt]
+        nao_obs[tt] <- PF[tt]
       }
 
       if (!is.null(exp)) {
@@ -478,106 +480,94 @@ NAO <- function(exp = NULL, obs = NULL, exp_cor = NULL, lat, lon, time_dim = 'sd
           exp_sub <- exp[, c(1:ntime)[-tt], , , drop = FALSE]
           # Combine 'memb' and 'sdate' to calculate EOF
           dim(exp_sub) <- c(nmemb_exp * (ntime - 1), nlat, nlon)
-          exp_EOF <- .EOF(exp_sub, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
+          EOF_exp <- .EOF(exp_sub, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
 
-          ## Correct polarity of pattern.
+          ## Correct polarity of pattern
           ##NOTE: different from s2dverification, which doesn't use mean().
-#          if (0 < exp_EOF$EOFs[1, which.min(abs(lat - 65)), ]) {
-          if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-            exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
-#            exp_EOF$PCs <- exp_EOF$PCs * sign  # not used
+#          if (0 < EOF_exp[1, which.min(abs(lat - 65)), ]) {
+          if (0 < mean(EOF_exp[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+            EOF_exp <- EOF_exp * (-1)
           }
 
           ### Lines below could be simplified further by computing
           ### ProjectField() only on the year of interest... (though this is
           ### not vital). Lauriane
           for (imemb in 1:nmemb_exp) {
-            PF <- .ProjectField(exp[imemb, , , ], eof_mode = exp_EOF$EOFs[1, , ], wght = wght) # [sdate, memb]
-            NAOF.ver[tt, imemb] <- PF[tt]
+            PF <- .ProjectField(exp[imemb, , , ], eof_mode = EOF_exp[1, , ], wght = wght) # [sdate, memb]
+            nao_exp[tt, imemb] <- PF[tt]
           }    
         } else {
           ## Project forecast anomalies on obs EOF
           for (imemb in 1:nmemb_exp) {
-            PF <- .ProjectField(exp[imemb, , , ], eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-            NAOF.ver[tt, imemb] <- PF[tt]
+            PF <- .ProjectField(exp[imemb, , , ], eof_mode = EOF_obs[1, , ], wght = wght)  # [sdate]
+            nao_exp[tt, imemb] <- PF[tt]
           }
         }
       }
 
     }  # for loop sdate
 
-    # add_member_back
-    if (add_member_back) {
-      suppressWarnings( 
-      NAOO.ver <- InsertDim(NAOO.ver, 2, 1, name = names(dim(exp))[1])
-      )
-    }
-
   } else { # exp_cor provided
+
     ## Calculate observation EOF. Without cross-validation
-    obs_EOF <- .EOF(obs, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
-    ## Correct polarity of pattern.
-    # dim(obs_EOF$EOFs): [mode, lat, lon]
-    if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-      obs_EOF$EOFs <- obs_EOF$EOFs * (-1)
-      #          obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
-    }
-    ## Project observed anomalies.
-    PF <- .ProjectField(obs, eof_mode = obs_EOF$EOFs, wght = wght)  # [sdate]
-    NAOO.ver[] <- PF[1,]
-    if (obsproj) {
-      ## Project forecasts anomalies 
-      for (imemb in 1:dim(exp_cor)[1]) {
-        exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
-                                         drop = 'selected')
-        PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-        NAOF.cor[,imemb] <- PF
-      }
-      for (imemb in 1:dim(exp)[1]) {
-        exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb,
-                                         drop = 'selected')
-        PF <- .ProjectField(exp_sub, eof_mode = obs_EOF$EOFs[1, , ], wght = wght)
-        NAOF.ver[,imemb] <- PF
-      }
-    } else if (!obsproj) {
-      dim(exp) <- c(nmemb_exp*ntime, nlat, nlon)
-      exp_EOF <- .EOF(exp, neofs = 1, wght = wght)  # $EOFs: [mode, lat, lon]
-      ## Correct polarity of pattern.
-      # dim(obs_EOF$EOFs): [mode, lat, lon]
-      if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
-        exp_EOF$EOFs <- exp_EOF$EOFs * (-1)
-#       obs_EOF$PCs <- obs_EOF$PCs * (-1)  # not used
-      }
-      dim(exp) <- c(nmemb_exp, ntime, nlat, nlon)
-      for (imemb in 1:dim(exp)[1]) {
-        exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, drop = 'selected')
-        PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-        NAOF.ver[ ,imemb] <- PF
+    EOF_obs <- .EOF(obs, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
+    ## Correct polarity of pattern
+    # EOF_obs: [mode, lat, lon]
+    if (0 < mean(EOF_obs[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+      EOF_obs <- EOF_obs * (-1)
+    }
+    ## Project observed anomalies
+    PF <- .ProjectField(obs, eof_mode = EOF_obs, wght = wght)  # [mode = 1, sdate]
+    nao_obs[] <- PF[1, ]
+
+    if (!obsproj) {
+      # Calculate EOF_exp
+      tmp <- array(exp, dim = c(nmemb_exp * ntime, nlat, nlon))
+      EOF_exp <- .EOF(tmp, neofs = 1, wght = wght)$EOFs  # [mode = 1, lat, lon]
+      ## Correct polarity of pattern
+      if (0 < mean(EOF_exp[1, which.min(abs(lat - 65)), ], na.rm = T)) {
+        EOF_exp <- EOF_exp * (-1)
       }
-      # Flip the dimensions back
-      NAOF.ver <- aperm(NAOF.ver, 2:1)
+      eof_mode_input <- EOF_exp[1, , ]
+    } else {
+      eof_mode_input <- EOF_obs[1, , ]
+    }
 
-      ## Project hindcast anomalies on forecast
-      for (imemb in 1:dim(exp_cor)[1]) {
-        exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
-                                         drop = 'selected')
+    # Calculate NAO_exp
+    for (imemb in 1:dim(exp)[1]) {
+      exp_sub <- ClimProjDiags::Subset(exp, along = 1, indices = imemb, 
+                                       drop = 'selected')
+      PF <- .ProjectField(exp_sub, eof_mode = eof_mode_input, wght = wght)  # [sdate]
+      nao_exp[ , imemb] <- PF
+    }
 
-        PF <- .ProjectField(exp_sub, eof_mode = exp_EOF$EOFs[1, , ], wght = wght)  # [sdate]
-        NAOF.cor[imemb] <- PF
-      }
+    # Calculate NAO_exp_cor
+    for (imemb in 1:dim(exp_cor)[1]) {
+      exp_sub <- ClimProjDiags::Subset(exp_cor, along = 1, indices = imemb,
+                                       drop = 'selected')
+      PF <- .ProjectField(exp_sub, eof_mode = eof_mode_input, wght = wght)  # [sdate]
+      nao_exp_cor[, imemb] <- PF
     }
+
   }
-  #NOTE: EOFs_obs is not returned because it's only the result of the last sdate
-  #     (It is returned in s2dverification.)
+  # add_member_back
+  if (add_member_back) {
+    memb_dim_name <- ifelse(!is.null(names(dim(exp))[1]), names(dim(exp))[1], 'member')
+    nao_obs <- InsertDim(nao_obs, 2, 1, name = memb_dim_name)
+  }
+
+  # Return results
   if (is.null(exp_cor)) {
-    if (!is.null(exp) & !is.null(obs)) {
-      return(list(exp = NAOF.ver, obs = NAOO.ver)) #, EOFs_obs = obs_EOF))
-    } else if (!is.null(exp)) {
-      return(list(exp = NAOF.ver))
-    } else if (!is.null(obs)) {
-      return(list(obs = NAOO.ver))
+    res <- NULL
+    if (!is.null(exp)) {
+      res <- c(res, list(exp = nao_exp))
     }
+    if (!is.null(obs)) {
+      res <- c(res, list(obs = nao_obs))
+    }
+    return(res)
+
   } else {
-    return(list(exp = NAOF.ver, obs = NAOO.ver, exp_cor = NAOF.cor))
+    return(list(exp = nao_exp, obs = nao_obs, exp_cor = nao_exp_cor))
   }
 }
diff --git a/man/NAO.Rd b/man/NAO.Rd
index c8ab94a..8545e6b 100644
--- a/man/NAO.Rd
+++ b/man/NAO.Rd
@@ -63,12 +63,12 @@ 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 
-reference (TRUE) or compute the NAO by first computing the leading 
-EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
-year you are evaluating out), and then projecting forecast anomalies onto 
-this EOF (FALSE). The default value is TRUE. If 'exp_cor' is provided, 'obs'
-will be used when obsproj is TRUE and 'exp' will be used when obsproj is 
-FALSE.}
+reference (TRUE, default) or compute the NAO by first computing the leading
+EOF of the forecast anomalies (in cross-validation mode, i.e. leave the 
+evaluated year out), then projecting forecast anomalies onto this EOF
+(FALSE). If 'exp_cor' is provided, 'obs' will be used when obsproj is TRUE
+and 'exp' will be used when obsproj is FALSE, and no cross-validation is
+applied.}
 
 \item{ncores}{An integer indicating the number of cores to use for parallel 
 computation. The default value is NULL.}
@@ -101,7 +101,7 @@ By default (ftime_avg = 2:4), NAO() computes the NAO index for 1-month
 lead seasonal forecasts that can be plotted with PlotBoxWhisker(). It returns
 cross-validated PCs of the NAO index for hindcast (exp) and observations 
 (obs) based on the leading EOF pattern, or, if forecast (exp_cor) is provided,
-the NAO index for forecast and the chosen corresponding data (exp or obs).
+the NAO index for forecast and the corresponding data (exp and obs).
 }
 \examples{
 # Make up synthetic data
@@ -114,7 +114,7 @@ lon <- seq(-80, 40, length.out = 9)
 nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 
 exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
-nao <- NAO(exp = NULL, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon)
+nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T)
 # plot the NAO index
  \dontrun{
 nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/tests/testthat/test-NAO.R b/tests/testthat/test-NAO.R
index 7b64cac..c3f0f75 100644
--- a/tests/testthat/test-NAO.R
+++ b/tests/testthat/test-NAO.R
@@ -6,7 +6,9 @@
   obs1 <- array(rnorm(72), dim = c(dataset = 1, member = 1, sdate = 3, ftime = 4, lat = 2, lon = 3))
   lat1 <- c(20, 80)
   lon1 <- c(40, 280, 350)
-
+  set.seed(3)
+  exp1_cor <- array(rnorm(72), dim = c(sdate = 1, ftime = 4, member = 3, lat = 2, lon = 3))
+ 
   # dat2
   set.seed(1)
   exp2 <- array(rnorm(216), dim = c(sdate = 3, ftime = 4, member = 2, lat = 3, lon = 3))
@@ -221,6 +223,38 @@ test_that("2. dat1", {
   c(NAO(exp1, obs1, lat = lat1, lon = lon1, ftime_avg = 1)$exp)
   )
 
+  # exp_cor
+  expect_equal(
+  names(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)),
+  c("exp", "obs", "exp_cor")
+  )
+  expect_equal(
+  dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp_cor),
+  c(sdate = 1, member = 3, dataset = 1)
+  )
+  expect_equal(
+  dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp),
+  c(sdate = 3, member = 2, dataset = 1)
+  )
+  expect_equal(
+  dim(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$obs),
+  c(sdate = 3, member = 1, dataset = 1)
+  )
+  expect_equal(
+  c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp_cor),
+  c(0.3896168, 0.4384543, -0.1302738),
+  tolerance = 0.0001
+  )
+  expect_equal(
+  c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$exp),
+  c(-0.1688756, -0.2658420, -0.8049575, -0.3022108, -0.3655258, 0.1237722),
+  tolerance = 0.0001
+  )
+    expect_equal(
+  c(NAO(exp = exp1, obs = obs1, exp_cor = exp1_cor, lat = lat1, lon = lon1)$obs),
+  c(-0.1762489, 0.1364694, -1.4581406),
+  tolerance = 0.0001
+  )
 
 })
 ##############################################
@@ -256,39 +290,34 @@ test_that("3. dat2", {
 
   # exp_cor
   expect_equal(
-  dim(NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor),
-  c(sdate = 1, member = 2)
+  names(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)),
+  c("exp", "obs", "exp_cor")
   )
   expect_equal(
-  dim(NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs),
-  c(sdate = 3)
+  dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor),
+  c(sdate = 1, member = 2)
   )
   expect_equal(
-  c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)),
-  c(0.2121340, 0.1634516),
-  tolerance = 0.0001)
-  expect_equal(
-  c((NAO(exp = NULL, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)),
-  c(0.3511294, -0.7196260, -1.5123894),
-  tolerance = 0.0001
+  dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp),
+  c(sdate = 3, member = 2)
   )
-
   expect_equal(
-  dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor),
-  c(sdate = 1, member = 2)
+  dim(NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs),
+  c(sdate = 3)
   )
   expect_equal(
-  dim(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp),
-  c(sdate = 3, member = 2)
+  c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp_cor)),
+  c(0.2121340, 0.1634516),
+  tolerance = 0.0001
   )
   expect_equal(
-  c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp),
-  c(-0.4225802, -1.3616918, 0.2808729, 0.2723077, -0.7584804, -0.7264514),
+  c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$exp)),
+  c(0.01457391, 0.06668166, 0.20193275, -0.20154315, -0.49487925, -0.04181974),
   tolerance = 0.0001
-  )
+  ) 
   expect_equal(
-  c(NAO(exp = exp2, obs = NULL, exp_cor = exp2_cor, lat = lat2, lon = lon2, obsproj = F)$exp_cor),
-  c(-0.22949531, -0.06946422),
+  c((NAO(exp = exp2, obs = obs2, exp_cor = exp2_cor, lat = lat2, lon = lon2)$obs)),
+  c(0.3511294, -0.7196260, -1.5123894),
   tolerance = 0.0001
   )
 
-- 
GitLab


From 5e75233a47c0182f498f6b3cc43c8e0b5fe24910 Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 15 Nov 2023 14:57:25 +0100
Subject: [PATCH 11/12] Change T to TRUE

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

diff --git a/R/NAO.R b/R/NAO.R
index 2a0e118..0810f2b 100644
--- a/R/NAO.R
+++ b/R/NAO.R
@@ -88,7 +88,7 @@
 #'nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 #'
 #'exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
-#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T)
+#'nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = TRUE)
 #'# plot the NAO index
 #'  \dontrun{
 #'nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/man/NAO.Rd b/man/NAO.Rd
index 8545e6b..8115fa2 100644
--- a/man/NAO.Rd
+++ b/man/NAO.Rd
@@ -114,7 +114,7 @@ lon <- seq(-80, 40, length.out = 9)
 nao <- NAO(exp = exp, obs = obs, lat = lat, lon = lon)
 
 exp_cor <- array(rnorm(540), dim = c(member = 2, sdate = 1, ftime = 5, lat = 6, lon = 9))
-nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = T)
+nao <- NAO(exp = exp, obs = obs, exp_cor = exp_cor, lat = lat, lon = lon, obsproj = TRUE)
 # plot the NAO index
  \dontrun{
 nao$exp <- Reorder(nao$exp, c(2, 1))
-- 
GitLab


From 8fb2a2e8576e4efc87c926e5070e3384eefc545d Mon Sep 17 00:00:00 2001
From: aho <an.ho@bsc.es>
Date: Wed, 15 Nov 2023 15:27:35 +0100
Subject: [PATCH 12/12] Fix example for new NAO dev

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

diff --git a/R/PlotBoxWhisker.R b/R/PlotBoxWhisker.R
index 2ddcec0..9b4c88e 100644
--- a/R/PlotBoxWhisker.R
+++ b/R/PlotBoxWhisker.R
@@ -90,7 +90,7 @@
 #'ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
 #'ano_exp <- array(ano$exp, dim = dim(ano$exp)[-2])
 #'ano_obs <- array(ano$obs, dim = dim(ano$obs)[-2])
-#'nao <- NAO(ano_exp, ano_obs, sampleData$lat, sampleData$lon)
+#'nao <- NAO(exp = ano_exp, obs = ano_obs, lat = sampleData$lat, lon = sampleData$lon)
 #'# Finally plot the nao index
 #'  \dontrun{
 #'nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/R/PlotStereoMap.R b/R/PlotStereoMap.R
index 4b4fbd2..a2fdb85 100644
--- a/R/PlotStereoMap.R
+++ b/R/PlotStereoMap.R
@@ -170,8 +170,10 @@
 #'data <- matrix(rnorm(100 * 50), 100, 50)
 #'x <- seq(from = 0, to = 360, length.out = 100)
 #'y <- seq(from = -90, to = 90, length.out = 50)
+#'  \dontrun{
 #'PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50,
 #'              toptitle = "This is the title")
+#' }
 #'@import mapproj
 #'@importFrom grDevices dev.cur dev.new dev.off gray
 #'@importFrom stats median
diff --git a/man/PlotBoxWhisker.Rd b/man/PlotBoxWhisker.Rd
index 9c5a3f4..e868160 100644
--- a/man/PlotBoxWhisker.Rd
+++ b/man/PlotBoxWhisker.Rd
@@ -123,7 +123,7 @@ sampleData$lat[] <- c(20, 80)
 ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
 ano_exp <- array(ano$exp, dim = dim(ano$exp)[-2])
 ano_obs <- array(ano$obs, dim = dim(ano$obs)[-2])
-nao <- NAO(ano_exp, ano_obs, sampleData$lat, sampleData$lon)
+nao <- NAO(exp = ano_exp, obs = ano_obs, lat = sampleData$lat, lon = sampleData$lon)
 # Finally plot the nao index
  \dontrun{
 nao$exp <- Reorder(nao$exp, c(2, 1))
diff --git a/man/PlotStereoMap.Rd b/man/PlotStereoMap.Rd
index 1b7f166..10e7503 100644
--- a/man/PlotStereoMap.Rd
+++ b/man/PlotStereoMap.Rd
@@ -283,6 +283,8 @@ compatible with figure layouts if colour bar is disabled.
 data <- matrix(rnorm(100 * 50), 100, 50)
 x <- seq(from = 0, to = 360, length.out = 100)
 y <- seq(from = -90, to = 90, length.out = 50)
+ \dontrun{
 PlotStereoMap(data, x, y, latlims = c(60, 90), brks = 50,
              toptitle = "This is the title")
 }
+}
-- 
GitLab