From 2a3c461c16cc5b64292a80f93cd8f937ab525508 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 22 Feb 2024 18:09:17 +0100
Subject: [PATCH 01/28] first version spread error ratio and significance

---
 R/ratio.R | 209 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 209 insertions(+)
 create mode 100644 R/ratio.R

diff --git a/R/ratio.R b/R/ratio.R
new file mode 100644
index 0000000..14bb36f
--- /dev/null
+++ b/R/ratio.R
@@ -0,0 +1,209 @@
+#'Compute the ratio between the ensemble spread and RMSE
+#'
+#'Compute the ratio between the standard deviation of the members around the 
+#'ensemble mean in experimental data and the RMSE between the ensemble mean of 
+#'experimental and observational data. The p-value is provided by a one-sided 
+#'Fisher's test.
+#'
+#'@param exp A named numeric array of experimental data with at least two 
+#'  dimensions 'memb_dim' and 'time_dim'.
+#'@param obs A named numeric array of observational data with at least two 
+#'  dimensions 'memb_dim' and 'time_dim'. It should have the same dimensions as
+#'  parameter 'exp' except along 'dat_dim' and 'memb_dim'.
+#'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
+#'  dimension. The default value is NULL (no dataset).
+#'@param memb_dim A character string indicating the name of the member 
+#'  dimension. It must be one dimension in 'exp' and 'obs'. The default value 
+#'  is 'member'.
+#'@param time_dim A character string indicating the name of dimension along  
+#'  which the ratio is computed. The default value is 'sdate'.
+#'@param pval A logical value indicating whether to compute or not the p-value 
+#'  of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return A list of two arrays with dimensions c(nexp, nobs, the rest of 
+#'  dimensions of 'exp' and 'obs' except memb_dim and time_dim), which nexp is
+#'  the length of dat_dim of 'exp' and nobs is the length of dat_dim of 'obs'. 
+#'  If dat_dim is NULL, nexp and nobs are omitted. \cr
+#'\item{$ratio}{
+#'  The ratio of the ensemble spread and RMSE.
+#'}
+#'\item{$p_val}{
+#'  The p-value of the one-sided Fisher's test with Ho: SD/RMSE = 1. Only present
+#'  if \code{pval = TRUE}.
+#'}
+#'
+#'@examples
+#'# Load sample data as in Load() example:
+#'example(Load)
+#'rsdrms <- RatioSDRMS(sampleData$mod, sampleData$obs, dat_dim = 'dataset')
+#'# Reorder the data in order to plot it with PlotVsLTime
+#'rsdrms_plot <- array(dim = c(dim(rsdrms$ratio)[1:2], 4, dim(rsdrms$ratio)[3]))
+#'rsdrms_plot[, , 2, ] <- rsdrms$ratio
+#'rsdrms_plot[, , 4, ] <- rsdrms$p.val
+#'\dontrun{
+#'PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", 
+#'            monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), 
+#'            listobs = c('ERSST'), biglab = FALSE, siglev = TRUE) 
+#'}
+#'
+#'@import multiApply
+#'@export
+Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member', 
+                       time_dim = 'sdate', pval = TRUE, sign = FALSE, alpha = 0.05,
+                       ncores = NULL) {
+
+  # Check inputs 
+  ## exp and obs (1)
+  if (is.null(exp) | is.null(obs)) {
+    stop("Parameter 'exp' and 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(exp) | !is.numeric(obs)) {
+    stop("Parameter 'exp' and 'obs' must be a numeric array.")
+  }
+  if (is.null(dim(exp)) | is.null(dim(obs))) {
+    stop(paste0("Parameter 'exp' and 'obs' must be array with as least two ",
+                "dimensions memb_dim and time_dim."))
+  }
+  if (any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
+      any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  ## dat_dim
+  if (!is.null(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.")
+    }
+  }
+  ## memb_dim
+  if (!is.character(memb_dim) | length(memb_dim) > 1) {
+    stop("Parameter 'memb_dim' must be a character string.")
+  }
+  if (!memb_dim %in% names(dim(exp)) & !memb_dim %in% names(dim(obs))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' nor 'obs' dimension. ",
+         "Set it as NULL if there is no member dimension.")
+  }
+  # Add [member = 1] 
+  if (memb_dim %in% names(dim(exp)) & !memb_dim %in% names(dim(obs))) {
+    dim(obs) <- c(dim(obs), 1)
+    names(dim(obs))[length(dim(obs))] <- memb_dim
+  }
+  if (!memb_dim %in% names(dim(exp)) & memb_dim %in% names(dim(obs))) {
+    dim(exp) <- c(dim(exp), 1)
+    names(dim(exp))[length(dim(exp))] <- memb_dim
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) > 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  if (!is.null(dat_dim)) {
+    name_exp <- name_exp[-which(name_exp == dat_dim)]
+    name_obs <- name_obs[-which(name_obs == dat_dim)]
+  }
+  name_exp <- name_exp[-which(name_exp == memb_dim)]
+  name_obs <- name_obs[-which(name_obs == memb_dim)]
+  if (!identical(dim(exp)[name_exp], dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all the dimensions except 'dat_dim' and 'memb_dim'."))
+  }
+  ## pval
+  if (!is.logical(pval) | length(pval) > 1) {
+    stop("Parameter 'pval' must be one logical value.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  } 
+
+
+  ###############################
+  # Calculate RatioSDRMS
+
+  # If dat_dim = NULL, insert dat dim
+  remove_dat_dim <- FALSE
+  if (is.null(dat_dim)) {
+    dat_dim <- 'dataset'
+    exp <- InsertDim(exp, posdim = 1, lendim = 1, name = 'dataset')
+    obs <- InsertDim(obs, posdim = 1, lendim = 1, name = 'dataset')
+    remove_dat_dim <- TRUE
+  }
+
+  res <- Apply(list(exp, obs), 
+               target_dims = list(c(dat_dim, memb_dim, time_dim), 
+                                  c(dat_dim, memb_dim, time_dim)),
+               pval = pval,  
+               fun = .Ratio, 
+               ncores = ncores)
+
+  if (remove_dat_dim) {
+    if (length(dim(res[[1]])) > 2) {
+      res <- lapply(res, Subset, c('nexp', 'nobs'), list(1, 1), drop = 'selected')
+    } else {
+      res <- lapply(res, as.numeric)
+    }
+  }
+
+  return(res)
+}
+
+.Ratio <- function(exp, obs, pval = TRUE, sign = FALSE, alpha = 0.05) {
+
+  # exp: [dat_exp, member, sdate]
+  # obs: [dat_obs, member, sdate]
+  nexp <- dim(exp)[1]
+  nobs <- dim(obs)[1]
+
+  # ensemble mean
+  ens_exp <- MeanDims(exp, 2, na.rm = TRUE) # [dat, sdate]
+  ens_obs <- MeanDims(obs, 2, na.rm = TRUE)
+
+  noise <- exp - InsertDim(ens_exp, 2, dim(exp)[2]) # [nexp, member, sdate]
+  enonoise <- rowSums(Eno(noise, names(dim(exp))[3]), na.rm = TRUE)
+  
+  # Create empty arrays
+  ratio <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
+  p.val <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
+
+  for (jexp in 1:nexp) {
+    for (jobs in 1:nobs) {
+      spread <- sqrt(mean(apply(exp[jexp,,], 1, var, na.rm = T), na.rm = T))
+      error <- sqrt(mean((ens_obs - ens_exp[jexp,])**2, na.rm = T))
+      ratio[jexp, jobs] <- spread/error
+
+      #effective sample size      
+      dif <- ens_exp[jexp, ] - ens_obs[jobs, ]
+      enodif <- Eno(dif)
+      
+      if (pval) {
+        F <- (enonoise[jexp] * spread[jexp]^2 / (enonoise[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
+        if (!is.na(F) & !is.na(enonoise[jexp]) & !is.na(enodif) & any(enonoise > 2) & enodif > 2) {
+          p.val[jexp, jobs] <- 1 - pf(F, enonoise[jexp] - 1, enodif - 1)
+        } else {
+          ratiosdrms[jexp, jobs] <- NA
+        }
+      }
+    }
+  }
+
+  res <- list(ratio = ratio)
+  if (pval) res <- c(res, list(p.val = p.val))
+  if (sign) {
+    sign <- p.val < alpha/2
+    res <- c(res, list(sign = sign))
+  }
+  
+  return(res)
+}
-- 
GitLab


From c3863810952fad70a0d284cad568c953c6894a47 Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Tue, 27 Feb 2024 14:41:08 +0100
Subject: [PATCH 02/28] renamed

---
 R/{ratio.R => SprErr.R} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename R/{ratio.R => SprErr.R} (100%)

diff --git a/R/ratio.R b/R/SprErr.R
similarity index 100%
rename from R/ratio.R
rename to R/SprErr.R
-- 
GitLab


From d642a820d46eb818d486c4465499dbe6b17a3d9c Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Tue, 27 Feb 2024 14:42:14 +0100
Subject: [PATCH 03/28] corrected ratio name

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

diff --git a/R/SprErr.R b/R/SprErr.R
index 14bb36f..d789eec 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -192,7 +192,7 @@ Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
         if (!is.na(F) & !is.na(enonoise[jexp]) & !is.na(enodif) & any(enonoise > 2) & enodif > 2) {
           p.val[jexp, jobs] <- 1 - pf(F, enonoise[jexp] - 1, enodif - 1)
         } else {
-          ratiosdrms[jexp, jobs] <- NA
+          ratio[jexp, jobs] <- NA
         }
       }
     }
-- 
GitLab


From 7b7772660d527b9c32a4132cd432b644af334723 Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Tue, 27 Feb 2024 15:09:25 +0100
Subject: [PATCH 04/28] included na.rm

---
 R/SprErr.R | 65 +++++++++++++++++++++++++++++-------------------------
 1 file changed, 35 insertions(+), 30 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index d789eec..a74fa64 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -1,9 +1,9 @@
 #'Compute the ratio between the ensemble spread and RMSE
 #'
-#'Compute the ratio between the standard deviation of the members around the 
+#'Compute the ratio between the spread of the members around the 
 #'ensemble mean in experimental data and the RMSE between the ensemble mean of 
-#'experimental and observational data. The p-value is provided by a one-sided 
-#'Fisher's test.
+#'experimental and observational data. The p-value and/or the statistical 
+#'significance is provided by a one-sided Fisher's test.
 #'
 #'@param exp A named numeric array of experimental data with at least two 
 #'  dimensions 'memb_dim' and 'time_dim'.
@@ -19,6 +19,8 @@
 #'  which the ratio is computed. The default value is 'sdate'.
 #'@param pval A logical value indicating whether to compute or not the p-value 
 #'  of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.
+#'@param na.rm A logical value indicating whether to remove NA values. The default 
+#'  value is TRUE.
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
@@ -50,10 +52,10 @@
 #'
 #'@import multiApply
 #'@export
-Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member', 
-                       time_dim = 'sdate', pval = TRUE, sign = FALSE, alpha = 0.05,
-                       ncores = NULL) {
-
+SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member', 
+                   time_dim = 'sdate', pval = TRUE, sign = FALSE, 
+                   alpha = 0.05, na.rm = FALSE, ncores = NULL) {
+  
   # Check inputs 
   ## exp and obs (1)
   if (is.null(exp) | is.null(obs)) {
@@ -123,15 +125,15 @@ Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
   ## ncores
   if (!is.null(ncores)) {
     if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
-      length(ncores) > 1) {
+        length(ncores) > 1) {
       stop("Parameter 'ncores' must be a positive integer.")
     }
   } 
-
-
+  
+  
   ###############################
   # Calculate RatioSDRMS
-
+  
   # If dat_dim = NULL, insert dat dim
   remove_dat_dim <- FALSE
   if (is.null(dat_dim)) {
@@ -140,14 +142,16 @@ Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
     obs <- InsertDim(obs, posdim = 1, lendim = 1, name = 'dataset')
     remove_dat_dim <- TRUE
   }
-
+  
   res <- Apply(list(exp, obs), 
                target_dims = list(c(dat_dim, memb_dim, time_dim), 
                                   c(dat_dim, memb_dim, time_dim)),
                pval = pval,  
-               fun = .Ratio, 
+               sign = sign,
+               na.rm = na.rm,
+               fun = .SprErr, 
                ncores = ncores)
-
+  
   if (remove_dat_dim) {
     if (length(dim(res[[1]])) > 2) {
       res <- lapply(res, Subset, c('nexp', 'nobs'), list(1, 1), drop = 'selected')
@@ -155,34 +159,34 @@ Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       res <- lapply(res, as.numeric)
     }
   }
-
+  
   return(res)
 }
 
-.Ratio <- function(exp, obs, pval = TRUE, sign = FALSE, alpha = 0.05) {
-
+.SprErr <- function(exp, obs, pval = TRUE, sign = FALSE, alpha = 0.05, na.rm = FALSE) {
+  
   # exp: [dat_exp, member, sdate]
   # obs: [dat_obs, member, sdate]
   nexp <- dim(exp)[1]
   nobs <- dim(obs)[1]
-
+  
   # ensemble mean
-  ens_exp <- MeanDims(exp, 2, na.rm = TRUE) # [dat, sdate]
-  ens_obs <- MeanDims(obs, 2, na.rm = TRUE)
-
+  ens_exp <- MeanDims(exp, 2, na.rm = na.rm) # [dat, sdate]
+  ens_obs <- MeanDims(obs, 2, na.rm = na.rm)
+  
   noise <- exp - InsertDim(ens_exp, 2, dim(exp)[2]) # [nexp, member, sdate]
-  enonoise <- rowSums(Eno(noise, names(dim(exp))[3]), na.rm = TRUE)
+  enonoise <- rowSums(Eno(noise, names(dim(exp))[3]), na.rm = na.rm)
   
   # Create empty arrays
   ratio <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
   p.val <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
-
+  
   for (jexp in 1:nexp) {
     for (jobs in 1:nobs) {
-      spread <- sqrt(mean(apply(exp[jexp,,], 1, var, na.rm = T), na.rm = T))
-      error <- sqrt(mean((ens_obs - ens_exp[jexp,])**2, na.rm = T))
+      spread <- sqrt(mean(apply(exp[jexp,,], 1, var, na.rm = na.rm), na.rm = na.rm))
+      error <- sqrt(mean((ens_obs - ens_exp[jexp,])**2, na.rm = na.rm))
       ratio[jexp, jobs] <- spread/error
-
+      
       #effective sample size      
       dif <- ens_exp[jexp, ] - ens_obs[jobs, ]
       enodif <- Eno(dif)
@@ -197,12 +201,13 @@ Ratio <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       }
     }
   }
-
+  
   res <- list(ratio = ratio)
-  if (pval) res <- c(res, list(p.val = p.val))
+  if (pval) {
+    res$p.val <- p.val
+  }
   if (sign) {
-    sign <- p.val < alpha/2
-    res <- c(res, list(sign = sign))
+    res$sign <- p.val <= alpha/2
   }
   
   return(res)
-- 
GitLab


From 57ae4415f15ad5b7c104deee80664c4b2090d678 Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Tue, 27 Feb 2024 15:20:04 +0100
Subject: [PATCH 05/28] documentation and checks for pval, sign and alpha

---
 R/SprErr.R | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/R/SprErr.R b/R/SprErr.R
index a74fa64..97f65f6 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -19,6 +19,11 @@
 #'  which the ratio is computed. The default value is 'sdate'.
 #'@param pval A logical value indicating whether to compute or not the p-value 
 #'  of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.
+#'@param sign A logical value indicating whether to retrieve the statistical
+#'  significance of the test Ho: ACC = 0 based on 'alpha'. The default value is
+#'  FALSE.
+#'@param alpha A numeric indicating the significance level for the statistical
+#'  significance test. The default value is 0.05. 
 #'@param na.rm A logical value indicating whether to remove NA values. The default 
 #'  value is TRUE.
 #'@param ncores An integer indicating the number of cores to use for parallel 
@@ -122,6 +127,18 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
   if (!is.logical(pval) | length(pval) > 1) {
     stop("Parameter 'pval' must be one logical value.")
   }
+  ## sign
+  if (!is.logical(sign) | length(sign) > 1) {
+    stop("Parameter 'sign' must be one logical value.")
+  }
+  # alpha
+  if (!is.numeric(alpha) | any(alpha < 0) | any(alpha > 1) | length(alpha) > 1) {
+    stop("Parameter 'alpha' must be a numeric number between 0 and 1.")
+  }
+  # na.rm
+  if (!na.rm %in% c(TRUE, FALSE)) {
+    stop("Parameter 'na.rm' must be TRUE or FALSE")
+  }
   ## ncores
   if (!is.null(ncores)) {
     if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
-- 
GitLab


From a707a06546f55819f8ae5e136c5608f85f25a116 Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Tue, 27 Feb 2024 16:14:38 +0100
Subject: [PATCH 06/28] cleaned, and identical results to easyVerification
 (significance not tested)

---
 R/SprErr.R | 29 ++++++++++++-----------------
 1 file changed, 12 insertions(+), 17 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 97f65f6..04c4728 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -191,27 +191,26 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
   ens_exp <- MeanDims(exp, 2, na.rm = na.rm) # [dat, sdate]
   ens_obs <- MeanDims(obs, 2, na.rm = na.rm)
   
-  noise <- exp - InsertDim(ens_exp, 2, dim(exp)[2]) # [nexp, member, sdate]
-  enonoise <- rowSums(Eno(noise, names(dim(exp))[3]), na.rm = na.rm)
-  
   # Create empty arrays
   ratio <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
   p.val <- array(dim = c(nexp = as.numeric(nexp), nobs = as.numeric(nobs)))  # [nexp, nobs]
   
   for (jexp in 1:nexp) {
     for (jobs in 1:nobs) {
-      spread <- sqrt(mean(apply(exp[jexp,,], 1, var, na.rm = na.rm), na.rm = na.rm))
-      error <- sqrt(mean((ens_obs - ens_exp[jexp,])**2, na.rm = na.rm))
+      
+      # spread and error
+      spread <- sqrt(mean(apply(exp[jexp,,], 2, var, na.rm = na.rm), na.rm = na.rm))
+      error <- sqrt(mean((ens_obs - ens_exp[jexp,])^2, na.rm = na.rm))
       ratio[jexp, jobs] <- spread/error
       
-      #effective sample size      
-      dif <- ens_exp[jexp, ] - ens_obs[jobs, ]
-      enodif <- Eno(dif)
+      # effective sample size
+      enospr <- sum(Eno(exp[jexp, , ] - InsertDim(ens_exp, 2, dim(exp)[2])[jexp, , ]))
+      enodif <- Eno(ens_exp[jexp, ] - ens_obs[jobs, ])
       
       if (pval) {
-        F <- (enonoise[jexp] * spread[jexp]^2 / (enonoise[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
-        if (!is.na(F) & !is.na(enonoise[jexp]) & !is.na(enodif) & any(enonoise > 2) & enodif > 2) {
-          p.val[jexp, jobs] <- 1 - pf(F, enonoise[jexp] - 1, enodif - 1)
+        F <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
+        if (!is.na(F) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
+          p.val[jexp, jobs] <- 1 - pf(F, enospr[jexp] - 1, enodif - 1)
         } else {
           ratio[jexp, jobs] <- NA
         }
@@ -220,12 +219,8 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
   }
   
   res <- list(ratio = ratio)
-  if (pval) {
-    res$p.val <- p.val
-  }
-  if (sign) {
-    res$sign <- p.val <= alpha/2
-  }
+  if (pval) {res$p.val <- p.val}
+  if (sign) {res$sign <- p.val <= alpha/2}
   
   return(res)
 }
-- 
GitLab


From 9282d408ff80fcfb4048a2c2319c9736b1764617 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Tue, 2 Apr 2024 16:47:01 +0200
Subject: [PATCH 07/28] two sided fix in SprErr

---
 R/SprErr.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 04c4728..d051d1c 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -206,11 +206,11 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       # effective sample size
       enospr <- sum(Eno(exp[jexp, , ] - InsertDim(ens_exp, 2, dim(exp)[2])[jexp, , ]))
       enodif <- Eno(ens_exp[jexp, ] - ens_obs[jobs, ])
-      
       if (pval) {
         F <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
         if (!is.na(F) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
-          p.val[jexp, jobs] <- 1 - pf(F, enospr[jexp] - 1, enodif - 1)
+          p.val[jexp, jobs] <- pf(F, enospr[jexp] - 1, enodif - 1)
+          p.val[jexp, jobs] <- 2 * min(p.val[jexp, jobs], 1 - p.val[jexp, jobs])
         } else {
           ratio[jexp, jobs] <- NA
         }
@@ -220,7 +220,7 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
   
   res <- list(ratio = ratio)
   if (pval) {res$p.val <- p.val}
-  if (sign) {res$sign <- p.val <= alpha/2}
+  if (sign) {res$sign <- p.val <= alpha}
   
   return(res)
 }
-- 
GitLab


From d7cd00ad00eb60638c9f9ebfb728ea5fa21e51da Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 3 Apr 2024 10:17:19 +0200
Subject: [PATCH 08/28] Eno dim name fixed

---
 R/SprErr.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index d051d1c..19123b0 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -204,8 +204,8 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       ratio[jexp, jobs] <- spread/error
       
       # effective sample size
-      enospr <- sum(Eno(exp[jexp, , ] - InsertDim(ens_exp, 2, dim(exp)[2])[jexp, , ]))
-      enodif <- Eno(ens_exp[jexp, ] - ens_obs[jobs, ])
+      enospr <- sum(Eno(exp[jexp, , ] - InsertDim(ens_exp, 2, dim(exp)[2])[jexp, , ], names(dim(exp))[3]))
+      enodif <- .Eno(ens_exp[jexp, ] - ens_obs[jobs, ], na.action = na.pass)
       if (pval) {
         F <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
         if (!is.na(F) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
-- 
GitLab


From 96184581c60527bbba22c0bc22f1bf04d766f087 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 22 May 2024 11:38:57 +0200
Subject: [PATCH 09/28] squared diff sprerr

---
 R/SprErr.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 19123b0..15bcba3 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -204,8 +204,8 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       ratio[jexp, jobs] <- spread/error
       
       # effective sample size
-      enospr <- sum(Eno(exp[jexp, , ] - InsertDim(ens_exp, 2, dim(exp)[2])[jexp, , ], names(dim(exp))[3]))
-      enodif <- .Eno(ens_exp[jexp, ] - ens_obs[jobs, ], na.action = na.pass)
+      enospr <- sum(Eno(apply(exp[jexp,,], 2, var, na.rm = na.rm), names(dim(exp))[3]))
+      enodif <- .Eno((ens_exp[jexp, ] - ens_obs[jobs, ])^2, na.action = na.pass)
       if (pval) {
         F <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
         if (!is.na(F) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
-- 
GitLab


From f78150bb0651258426d0c528368306ed2ae2a596 Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Mon, 3 Jun 2024 12:09:16 +0200
Subject: [PATCH 10/28] fixed sign=T when pval=F

---
 R/SprErr.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 15bcba3..d0a4640 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -206,13 +206,13 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       # effective sample size
       enospr <- sum(Eno(apply(exp[jexp,,], 2, var, na.rm = na.rm), names(dim(exp))[3]))
       enodif <- .Eno((ens_exp[jexp, ] - ens_obs[jobs, ])^2, na.action = na.pass)
-      if (pval) {
+      if (pval | sign) {
         F <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
         if (!is.na(F) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
           p.val[jexp, jobs] <- pf(F, enospr[jexp] - 1, enodif - 1)
           p.val[jexp, jobs] <- 2 * min(p.val[jexp, jobs], 1 - p.val[jexp, jobs])
         } else {
-          ratio[jexp, jobs] <- NA
+          p.val[jexp, jobs] <- NA
         }
       }
     }
-- 
GitLab


From 3900200496bd683633d2fae1daab1303b2a3b0e8 Mon Sep 17 00:00:00 2001
From: Carlos Delgado <carlos.delgado@bsc.es>
Date: Mon, 3 Jun 2024 12:18:18 +0200
Subject: [PATCH 11/28] corrected example and documentation

---
 R/SprErr.R | 21 +++++++--------------
 1 file changed, 7 insertions(+), 14 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index d0a4640..a22bdf0 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -3,7 +3,7 @@
 #'Compute the ratio between the spread of the members around the 
 #'ensemble mean in experimental data and the RMSE between the ensemble mean of 
 #'experimental and observational data. The p-value and/or the statistical 
-#'significance is provided by a one-sided Fisher's test.
+#'significance is provided by a two-sided Fisher's test.
 #'
 #'@param exp A named numeric array of experimental data with at least two 
 #'  dimensions 'memb_dim' and 'time_dim'.
@@ -37,23 +37,16 @@
 #'  The ratio of the ensemble spread and RMSE.
 #'}
 #'\item{$p_val}{
-#'  The p-value of the one-sided Fisher's test with Ho: SD/RMSE = 1. Only present
+#'  The p-value of the two-sided Fisher's test with Ho: Spread/RMSE = 1. Only present
 #'  if \code{pval = TRUE}.
 #'}
 #'
 #'@examples
-#'# Load sample data as in Load() example:
-#'example(Load)
-#'rsdrms <- RatioSDRMS(sampleData$mod, sampleData$obs, dat_dim = 'dataset')
-#'# Reorder the data in order to plot it with PlotVsLTime
-#'rsdrms_plot <- array(dim = c(dim(rsdrms$ratio)[1:2], 4, dim(rsdrms$ratio)[3]))
-#'rsdrms_plot[, , 2, ] <- rsdrms$ratio
-#'rsdrms_plot[, , 4, ] <- rsdrms$p.val
-#'\dontrun{
-#'PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", 
-#'            monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), 
-#'            listobs = c('ERSST'), biglab = FALSE, siglev = TRUE) 
-#'}
+#'exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
+#'obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
+#'sprerr1 <- SprErr(exp, obs)
+#'sprerr2 <- SprErr(exp, obs, pval=F, sign=T)
+#'sprerr3 <- SprErr(exp, obs, pval=T, sign=T)
 #'
 #'@import multiApply
 #'@export
-- 
GitLab


From e36f7bc4a1a6ddda2c95d39c1b2e28db6fe138fe Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 13 Jun 2024 14:22:45 +0200
Subject: [PATCH 12/28] Correct formatting

---
 R/SprErr.R | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index a22bdf0..2abb7eb 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -17,15 +17,15 @@
 #'  is 'member'.
 #'@param time_dim A character string indicating the name of dimension along  
 #'  which the ratio is computed. The default value is 'sdate'.
-#'@param pval A logical value indicating whether to compute or not the p-value 
+#'@param pval A logical value indicating whether to compute the p-value 
 #'  of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.
 #'@param sign A logical value indicating whether to retrieve the statistical
 #'  significance of the test Ho: ACC = 0 based on 'alpha'. The default value is
 #'  FALSE.
 #'@param alpha A numeric indicating the significance level for the statistical
 #'  significance test. The default value is 0.05. 
-#'@param na.rm A logical value indicating whether to remove NA values. The default 
-#'  value is TRUE.
+#'@param na.rm A logical value indicating whether to remove NA values. The 
+#'  default value is TRUE.
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
@@ -37,16 +37,16 @@
 #'  The ratio of the ensemble spread and RMSE.
 #'}
 #'\item{$p_val}{
-#'  The p-value of the two-sided Fisher's test with Ho: Spread/RMSE = 1. Only present
-#'  if \code{pval = TRUE}.
+#'  The p-value of the two-sided Fisher's test with Ho: Spread/RMSE = 1. Only 
+#'  present if \code{pval = TRUE}.
 #'}
 #'
 #'@examples
 #'exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
 #'obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
 #'sprerr1 <- SprErr(exp, obs)
-#'sprerr2 <- SprErr(exp, obs, pval=F, sign=T)
-#'sprerr3 <- SprErr(exp, obs, pval=T, sign=T)
+#'sprerr2 <- SprErr(exp, obs, pval = F, sign = T)
+#'sprerr3 <- SprErr(exp, obs, pval = T, sign = T)
 #'
 #'@import multiApply
 #'@export
-- 
GitLab


From 77dd3ef4d1bb21b406d0295dc8532ba2d61cdfda Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 13 Jun 2024 14:23:15 +0200
Subject: [PATCH 13/28] Generate documentation for SprErr()

---
 NAMESPACE     |  1 +
 man/SprErr.Rd | 80 +++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 81 insertions(+)
 create mode 100644 man/SprErr.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 9214a1a..587aec2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -77,6 +77,7 @@ export(Season)
 export(SignalNoiseRatio)
 export(Smoothing)
 export(Spectrum)
+export(SprErr)
 export(Spread)
 export(StatSeasAtlHurr)
 export(TPI)
diff --git a/man/SprErr.Rd b/man/SprErr.Rd
new file mode 100644
index 0000000..5fe4882
--- /dev/null
+++ b/man/SprErr.Rd
@@ -0,0 +1,80 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/SprErr.R
+\name{SprErr}
+\alias{SprErr}
+\title{Compute the ratio between the ensemble spread and RMSE}
+\usage{
+SprErr(
+  exp,
+  obs,
+  dat_dim = NULL,
+  memb_dim = "member",
+  time_dim = "sdate",
+  pval = TRUE,
+  sign = FALSE,
+  alpha = 0.05,
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numeric array of experimental data with at least two 
+dimensions 'memb_dim' and 'time_dim'.}
+
+\item{obs}{A named numeric array of observational data with at least two 
+dimensions 'memb_dim' and 'time_dim'. It should have the same dimensions as
+parameter 'exp' except along 'dat_dim' and 'memb_dim'.}
+
+\item{dat_dim}{A character string indicating the name of dataset (nobs/nexp) 
+dimension. The default value is NULL (no dataset).}
+
+\item{memb_dim}{A character string indicating the name of the member 
+dimension. It must be one dimension in 'exp' and 'obs'. The default value 
+is 'member'.}
+
+\item{time_dim}{A character string indicating the name of dimension along
+which the ratio is computed. The default value is 'sdate'.}
+
+\item{pval}{A logical value indicating whether to compute the p-value 
+of the test Ho : SD/RMSE = 1 or not. The default value is TRUE.}
+
+\item{sign}{A logical value indicating whether to retrieve the statistical
+significance of the test Ho: ACC = 0 based on 'alpha'. The default value is
+FALSE.}
+
+\item{alpha}{A numeric indicating the significance level for the statistical
+significance test. The default value is 0.05.}
+
+\item{na.rm}{A logical value indicating whether to remove NA values. The 
+default value is TRUE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A list of two arrays with dimensions c(nexp, nobs, the rest of 
+ dimensions of 'exp' and 'obs' except memb_dim and time_dim), which nexp is
+ the length of dat_dim of 'exp' and nobs is the length of dat_dim of 'obs'. 
+ If dat_dim is NULL, nexp and nobs are omitted. \cr
+\item{$ratio}{
+ The ratio of the ensemble spread and RMSE.
+}
+\item{$p_val}{
+ The p-value of the two-sided Fisher's test with Ho: Spread/RMSE = 1. Only 
+ present if \code{pval = TRUE}.
+}
+}
+\description{
+Compute the ratio between the spread of the members around the 
+ensemble mean in experimental data and the RMSE between the ensemble mean of 
+experimental and observational data. The p-value and/or the statistical 
+significance is provided by a two-sided Fisher's test.
+}
+\examples{
+exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
+obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
+sprerr1 <- SprErr(exp, obs)
+sprerr2 <- SprErr(exp, obs, pval = F, sign = T)
+sprerr3 <- SprErr(exp, obs, pval = T, sign = T)
+
+}
-- 
GitLab


From 21c1af4b5d467761802a7274f5885e7062951c25 Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 13 Jun 2024 14:42:05 +0200
Subject: [PATCH 14/28] Fix lint: change T/F to TRUE/FALSE

---
 R/SprErr.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 2abb7eb..9892d94 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -45,8 +45,8 @@
 #'exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
 #'obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
 #'sprerr1 <- SprErr(exp, obs)
-#'sprerr2 <- SprErr(exp, obs, pval = F, sign = T)
-#'sprerr3 <- SprErr(exp, obs, pval = T, sign = T)
+#'sprerr2 <- SprErr(exp, obs, pval = FALSE, sign = TRUE)
+#'sprerr3 <- SprErr(exp, obs, pval = TRUE, sign = TRUE)
 #'
 #'@import multiApply
 #'@export
-- 
GitLab


From ee9965c00b3218b98176878b94ae3fea885900bc Mon Sep 17 00:00:00 2001
From: vagudets <victoria.agudetse@bsc.es>
Date: Thu, 13 Jun 2024 15:44:29 +0200
Subject: [PATCH 15/28] Update documentation

---
 man/SprErr.Rd | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/man/SprErr.Rd b/man/SprErr.Rd
index 5fe4882..cdc647e 100644
--- a/man/SprErr.Rd
+++ b/man/SprErr.Rd
@@ -74,7 +74,7 @@ significance is provided by a two-sided Fisher's test.
 exp <- array(rnorm(30), dim = c(lat = 2, sdate = 3, member = 5))
 obs <- array(rnorm(30), dim = c(lat = 2, sdate = 3))
 sprerr1 <- SprErr(exp, obs)
-sprerr2 <- SprErr(exp, obs, pval = F, sign = T)
-sprerr3 <- SprErr(exp, obs, pval = T, sign = T)
+sprerr2 <- SprErr(exp, obs, pval = FALSE, sign = TRUE)
+sprerr3 <- SprErr(exp, obs, pval = TRUE, sign = TRUE)
 
 }
-- 
GitLab


From 15778df06c752311d952aefc7c5d88ca917df322 Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Thu, 27 Jun 2024 17:24:42 +0200
Subject: [PATCH 16/28] Fix: added "alpha = alpha" when calling .SprErr in
 SprErr.R

---
 R/SprErr.R | 1 +
 1 file changed, 1 insertion(+)

diff --git a/R/SprErr.R b/R/SprErr.R
index 9892d94..c37a460 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -158,6 +158,7 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
                                   c(dat_dim, memb_dim, time_dim)),
                pval = pval,  
                sign = sign,
+               alpha = alpha,
                na.rm = na.rm,
                fun = .SprErr, 
                ncores = ncores)
-- 
GitLab


From a7c6abc0cf9be0b1c8d17db136a8583dc3b54541 Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Mon, 1 Jul 2024 16:10:58 +0200
Subject: [PATCH 17/28] Fix: check for member dimension in SprErr.R

---
 R/SprErr.R | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index c37a460..f70ffc5 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -83,16 +83,16 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
   if (!is.character(memb_dim) | length(memb_dim) > 1) {
     stop("Parameter 'memb_dim' must be a character string.")
   }
-  if (!memb_dim %in% names(dim(exp)) & !memb_dim %in% names(dim(obs))) {
-    stop("Parameter 'memb_dim' is not found in 'exp' nor 'obs' dimension. ",
-         "Set it as NULL if there is no member dimension.")
+  if (!memb_dim %in% names(dim(exp))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' dimensions. ",
+         "'exp' must have the member dimension to compute the spread.")
   }
   # Add [member = 1] 
   if (memb_dim %in% names(dim(exp)) & !memb_dim %in% names(dim(obs))) {
     dim(obs) <- c(dim(obs), 1)
     names(dim(obs))[length(dim(obs))] <- memb_dim
   }
-  if (!memb_dim %in% names(dim(exp)) & memb_dim %in% names(dim(obs))) {
+  if (!memb_dim %in% names(dim(exp)) & memb_dim %in% names(dim(obs))) { ## check no longer needed?
     dim(exp) <- c(dim(exp), 1)
     names(dim(exp))[length(dim(exp))] <- memb_dim
   }
-- 
GitLab


From e0810d68f2416bfc18577452b90ad629920f9bdb Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Mon, 8 Jul 2024 12:33:01 +0200
Subject: [PATCH 18/28] Replace as.numeric() with as.vector() when
 remove_dat_dim = TRUE in SprErr.R

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

diff --git a/R/SprErr.R b/R/SprErr.R
index f70ffc5..945c26f 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -167,7 +167,7 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
     if (length(dim(res[[1]])) > 2) {
       res <- lapply(res, Subset, c('nexp', 'nobs'), list(1, 1), drop = 'selected')
     } else {
-      res <- lapply(res, as.numeric)
+      res <- lapply(res, as.vector)
     }
   }
   
-- 
GitLab


From 7359c2602debf0ad8dc0116f09adcce0979f6f39 Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Mon, 8 Jul 2024 13:04:54 +0200
Subject: [PATCH 19/28] Replace "F" with "f_statistic" in .SprErr

---
 R/SprErr.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 945c26f..47f3cc4 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -201,9 +201,9 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       enospr <- sum(Eno(apply(exp[jexp,,], 2, var, na.rm = na.rm), names(dim(exp))[3]))
       enodif <- .Eno((ens_exp[jexp, ] - ens_obs[jobs, ])^2, na.action = na.pass)
       if (pval | sign) {
-        F <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
-        if (!is.na(F) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
-          p.val[jexp, jobs] <- pf(F, enospr[jexp] - 1, enodif - 1)
+        f_statistic <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
+        if (!is.na(f_statistic) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
+          p.val[jexp, jobs] <- pf(f_statistic, enospr[jexp] - 1, enodif - 1)
           p.val[jexp, jobs] <- 2 * min(p.val[jexp, jobs], 1 - p.val[jexp, jobs])
         } else {
           p.val[jexp, jobs] <- NA
-- 
GitLab


From f60c39ea6bd40ccf1d49fca4389c84974eec016b Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Mon, 8 Jul 2024 16:28:07 +0200
Subject: [PATCH 20/28] Fix  'enospr' and 'spread' indexing issue in .SprErr

---
 R/SprErr.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/SprErr.R b/R/SprErr.R
index 47f3cc4..c9fbb5c 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -201,9 +201,9 @@ SprErr <- function(exp, obs, dat_dim = NULL, memb_dim = 'member',
       enospr <- sum(Eno(apply(exp[jexp,,], 2, var, na.rm = na.rm), names(dim(exp))[3]))
       enodif <- .Eno((ens_exp[jexp, ] - ens_obs[jobs, ])^2, na.action = na.pass)
       if (pval | sign) {
-        f_statistic <- (enospr[jexp] * spread[jexp]^2 / (enospr[jexp] - 1)) / (enodif * error^2 / (enodif - 1))
-        if (!is.na(f_statistic) & !is.na(enospr[jexp]) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
-          p.val[jexp, jobs] <- pf(f_statistic, enospr[jexp] - 1, enodif - 1)
+        f_statistic <- (enospr * spread^2 / (enospr - 1)) / (enodif * error^2 / (enodif - 1))
+        if (!is.na(f_statistic) & !is.na(enospr) & !is.na(enodif) & any(enospr > 2) & enodif > 2) {
+          p.val[jexp, jobs] <- pf(f_statistic, enospr - 1, enodif - 1)
           p.val[jexp, jobs] <- 2 * min(p.val[jexp, jobs], 1 - p.val[jexp, jobs])
         } else {
           p.val[jexp, jobs] <- NA
-- 
GitLab


From f8fffa24b14798ae7793bf53cb670b73ea14a1c4 Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Fri, 12 Jul 2024 15:49:42 +0200
Subject: [PATCH 21/28] Bugfix: default value of na.rm is FALSE

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

diff --git a/R/SprErr.R b/R/SprErr.R
index c9fbb5c..f89d37a 100644
--- a/R/SprErr.R
+++ b/R/SprErr.R
@@ -25,7 +25,7 @@
 #'@param alpha A numeric indicating the significance level for the statistical
 #'  significance test. The default value is 0.05. 
 #'@param na.rm A logical value indicating whether to remove NA values. The 
-#'  default value is TRUE.
+#'  default value is FALSE.
 #'@param ncores An integer indicating the number of cores to use for parallel 
 #'  computation. The default value is NULL.
 #'
-- 
GitLab


From 1a79b299538d6f6536b731bb4635d834ce79e0c1 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Mon, 15 Jul 2024 16:14:49 +0200
Subject: [PATCH 22/28] Unit test for SprErr.R

---
 tests/testthat/test-SprErr.R | 600 +++++++++++++++++++++++++++++++++++
 1 file changed, 600 insertions(+)
 create mode 100644 tests/testthat/test-SprErr.R

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
new file mode 100644
index 0000000..a116bde
--- /dev/null
+++ b/tests/testthat/test-SprErr.R
@@ -0,0 +1,600 @@
+
+library(s2dv)
+library(testthat)
+library(multiApply)
+library(ClimProjDiags) # for Subset() [commit 53ba0039 by aho]
+source("R/SprErr.R")
+source("R/Eno.R")
+
+##############################################
+# data
+##############################################
+
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(60), dim = c(member = 3, sdate = 10, lat = 2))
+set.seed(2)
+obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+
+# dat1_2
+exp1_2 <-  exp1 + 1
+obs1_2 <-  obs1 + 1
+
+# dat1_3: NAs
+exp1_3 <- exp1; exp1_3[1, 2, 1] <- NA
+obs1_3 <- obs1; obs1_3[2, 1] <- NA
+
+# dat2
+set.seed(1)
+exp2 <- array(rnorm(20), dim = c(member = 2, sdate = 10))
+set.seed(2)
+obs2 <- array(rnorm(10), dim = c(sdate = 10))
+set.seed(2)
+obs2_1 <- array(rnorm(10), dim = c(member = 1, sdate = 10))
+
+# dat3
+set.seed(1)
+exp3 <- array(rnorm(60), dim = c(member = 2, sdate = 10, dataset = 3))
+set.seed(2)
+obs3 <- array(rnorm(20), dim = c(member = 1, sdate = 10, dataset = 2))
+
+# dat3_2
+set.seed(1)
+exp3_2 <- array(rnorm(80), dim = c(member = 4, sdate = 5, dataset = 4))
+set.seed(2)
+obs3_2 <- array(rnorm(30), dim = c(member = 2, sdate = 5, dataset = 3))
+
+# dat4
+set.seed(1)
+exp4 <- array(rnorm(60), dim = c(member = 2, sdate = 10, dataset = 3, lat = 2))
+set.seed(2)
+obs4 <- array(rnorm(20), dim = c(member = 1, sdate = 10, dataset = 2, lat = 2))
+
+# dat4_2: NAs
+exp4_2 <- exp4; exp4_2[1, 2, 1, 1] <- NA
+obs4_2 <- obs4; obs4_2[1, 1:4, 1, 1] <- NA
+
+
+##############################################
+# tests
+##############################################
+
+# exp and obs (1)
+test_that("1. Input checks", {
+  expect_error(
+    SprErr(c(), obs1), ## or: SprErr(exp1, obs1, c()) ?? exp1, c() i c() obs1 - què fa quan li passes algo que no és un array? c()
+    "Parameter 'exp' and 'obs' cannot be NULL." 
+  )
+  expect_error(
+    SprErr(obs1, c()), 
+    "Parameter 'exp' and 'obs' cannot be NULL." 
+  )
+  expect_error(
+    SprErr("", obs1),
+    "Parameter 'exp' and 'obs' must be a numeric array."
+  )
+  expect_error(
+    SprErr(exp1, ""), 
+    "Parameter 'exp' and 'obs' must be a numeric array."
+  )
+  expect_error(
+    SprErr(1, obs1), 
+    "Parameter 'exp' and 'obs' must be array with as least two dimensions memb_dim and time_dim." ## si es pot es evitar el paste0 millor
+  )
+  expect_error(
+    SprErr(exp1, 1), 
+    "Parameter 'exp' and 'obs' must be array with as least two dimensions memb_dim and time_dim." 
+  )
+  expect_error(
+    SprErr(array(1), obs1), 
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+    SprErr(exp1, array(1)), 
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  # dat_dim
+  expect_error(
+    SprErr(exp1, obs1, dat_dim = 1),
+    "Parameter 'dat_dim' must be a character string."
+  ) 
+  expect_error(
+    SprErr(exp1, obs1, dat_dim = 'dat'), 
+    "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
+  )  
+  # memb_dim
+  expect_error(
+    SprErr(exp1, obs1, memb_dim = 1),
+    "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+    SprErr(exp1, obs1, memb_dim = 'check'), 
+    "Parameter 'memb_dim' is not found in 'exp' dimensions. 'exp' must have the member dimension to compute the spread."
+  )
+  # time_dim
+  expect_error(
+    SprErr(exp1, obs1, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )  
+  expect_error(
+    SprErr(exp1, obs1, time_dim = 'time'),
+    "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # exp and obs (2)
+  expect_error(
+    SprErr(exp1, array(1:9, dim = c(sdate = 9))),
+    "Parameter 'exp' and 'obs' must have same length of all the dimensions except 'dat_dim' and 'memb_dim'."
+  )
+  # pval
+  expect_error(
+    SprErr(exp1, obs1, pval = 1),
+    "Parameter 'pval' must be one logical value."
+  )  
+  # sign
+  expect_error(
+    SprErr(exp1, obs1, sign = 1),
+    "Parameter 'sign' must be one logical value."
+  )  
+  # alpha
+  expect_error(
+    SprErr(exp1, obs1, alpha = -0.05),
+    "Parameter 'alpha' must be a numeric number between 0 and 1."
+  )   
+  # na.rm
+  expect_error(
+    SprErr(exp1, obs1, na.rm = ""), 
+    "Parameter 'na.rm' must be TRUE or FALSE"
+  )  
+  # ncores
+  expect_error(
+    SprErr(exp2, obs2, ncores = 1.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+})
+
+######################################################## #
+
+
+test_that("2. Output checks: dat1", {
+  
+  ## element names
+  expect_equal(
+    names(SprErr(exp1, obs1)), ## default: pval = TRUE, sign = FALSE
+    c("ratio", "p.val")
+  )
+  expect_equal(
+    names(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)),
+    c("ratio", "sign")
+  )
+  expect_equal(
+    names(SprErr(exp1, obs1, pval = TRUE, sign = TRUE)),
+    c("ratio", "p.val", "sign")
+  )
+  
+  ## dimensions
+  expect_equal(
+    dim(SprErr(exp1, obs1)$ratio),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp1, obs1)$p.val),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$ratio), 
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$sign), 
+    c(lat = 2)
+  )
+  
+  ## values
+  expect_equal(
+    as.vector(SprErr(exp1, obs1)$ratio),
+    c(1.0646692, 0.6657522),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1)$p.val),
+    c(0.8549593, 0.2412730),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE)$sign),
+    c(FALSE, FALSE)
+  )
+  ## pval = FALSE
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$ratio),
+    c(01.0646692, 0.6657522),  
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$sign),
+    c(FALSE, FALSE)  
+  )
+  ## na.rm = TRUE
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, na.rm = TRUE)$ratio),
+    c(1.0646692, 0.6657522),  
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, na.rm = TRUE)$p.val),
+    c(0.8549593, 0.2412730),  
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE, na.rm = TRUE)$sign),
+    c(FALSE, FALSE)
+  )
+  ## alpha
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE, alpha = 0.5)$sign),
+    c(FALSE, TRUE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp1, obs1, sign = TRUE, alpha = 0.99)$sign),
+    c(TRUE, TRUE)
+  )
+  
+  # dat1_2
+  expect_equal(
+    SprErr(exp1, obs1),
+    SprErr(exp1_2, obs1_2)
+  )
+  expect_equal(
+    SprErr(exp1, obs1, sign = TRUE)$ratio,
+    SprErr(exp1_2, obs1_2)$ratio
+  )
+  expect_equal(
+    SprErr(exp1, obs1, sign = TRUE)$p.val,
+    SprErr(exp1_2, obs1_2)$p.val
+  )
+  
+  # dat1_3
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3)$ratio),
+    c(NA, 0.6657522),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3)$p.val),
+    c(NA, 0.241273),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, sign = TRUE)$sign),
+    c(NA, FALSE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, na.rm = TRUE)$ratio),
+    c(0.9656329, 0.6657522),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, na.rm = TRUE)$p.val),
+    c(0.896455, 0.241273),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp1_3, obs1_3, sign = TRUE, na.rm = TRUE)$sign),
+    c(FALSE, FALSE)
+  )
+  
+})
+
+
+##############################################
+test_that("3. Output checks: dat2", {
+  
+  expect_equal(
+    names(SprErr(exp2, obs2)),
+    c("ratio", "p.val")
+  )
+  expect_equal(
+    names(SprErr(exp2, obs2, sign = TRUE)),
+    c("ratio", "p.val","sign")
+  )
+  expect_equal(
+    dim(SprErr(exp2, obs2)$ratio),
+    NULL
+  )
+  expect_equal(
+    dim(SprErr(exp2, obs2)$p.val),
+    NULL
+  )
+  expect_equal(
+    dim(SprErr(exp2, obs2, sign = TRUE)$sign),
+    NULL
+  )
+  ## values
+  expect_equal(
+    as.vector(SprErr(exp2, obs2)$ratio),
+    c(0.6866402),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2)$p.val),
+    c(0.2779936),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp2, obs2, sign = TRUE)$sign), 
+    FALSE 
+  )
+  ## sign = TRUE
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, sign = TRUE)$ratio),
+    c(0.6866402),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, sign = TRUE)$p.val),
+    c(0.2779936),
+    tolerance = 0.0001
+  )
+  ## alpha
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, sign = TRUE, alpha = 0.99)$sign),
+    TRUE 
+  )
+  ## na.rm = TRUE
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, na.rm = TRUE)$ratio),
+    c(0.6866402),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp2, obs2, na.rm = TRUE)$p.val),
+    c(0.2779936),
+    tolerance = 0.0001
+  )
+  expect_equal( # is this expected?
+    as.vector(SprErr(exp2, obs2, na.rm = TRUE, sign = TRUE)$sign), 
+    FALSE 
+  )
+  ## other
+  expect_equal(
+    SprErr(exp2, obs2),
+    SprErr(exp2, obs2_1)
+  )
+  
+})
+
+##############################################
+test_that("4. Output checks: dat3", {
+  
+  expect_equal(
+    dim(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
+    c('nexp' = 3, 'nobs' = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp3, obs3, dat_dim = 'dataset')$p.val),
+    c('nexp' = 3, 'nobs' = 2)
+  )
+  ## values
+  expect_equal(
+    mean(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset')$p.val)[1:3],
+    c(0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE)$sign)[1:3],
+    c(FALSE, FALSE, TRUE)
+  )
+  expect_equal(
+    mean(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE)$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE)$p.val)[1:3],
+    c(0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  ## alpha
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE, alpha = 0.99)$sign)[1:3],
+    c(TRUE, TRUE, TRUE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE, alpha = 0.20)$sign)[1:3],
+    c(FALSE, TRUE, TRUE)
+  )
+  ## na.rm = TRUE
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE)$ratio),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE)$p.val)[1:3],
+    c(0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE, sign = TRUE)$sign)[1:3], 
+    c(FALSE, FALSE, TRUE) 
+  )
+  
+  # dat3_2
+  expect_equal(
+    dim(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$ratio),
+    c('nexp' = 4, 'nobs' = 3)
+  )
+  expect_equal(
+    dim(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$p.val),
+    c('nexp' = 4, 'nobs' = 3)
+  )
+  ## values
+  expect_equal(
+    mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$ratio),
+    c(1.25586),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$p.val)[1:4],
+    c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE)$sign)[1:4],
+    c(FALSE, FALSE, FALSE, FALSE)
+  )
+  expect_equal(
+    mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE)$ratio),
+    c(1.25586),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE)$p.val)[1:4],
+    c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
+    tolerance = 0.0001
+  )
+  ## alpha
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE, alpha = 0.99)$sign)[1:3],
+    c(TRUE, TRUE, TRUE)
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE, alpha = 0.70)$sign)[1:4],
+    c(TRUE, FALSE, FALSE, TRUE)
+  )
+  ## na.rm = TRUE
+  expect_equal(
+    mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE)$ratio),
+    c(1.25586),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE)$p.val)[1:4],
+    c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE, 
+                     alpha = 0.70, sign = TRUE)$sign)[1:4], 
+    c(TRUE, FALSE, FALSE, TRUE) 
+  )
+  
+})
+
+
+##############################################
+test_that("5. Output checks: dat4", {
+  
+  expect_equal(
+    dim(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio),
+    c('nexp' = 3, 'nobs' = 2, 'lat' = 2)
+  )
+  expect_equal(
+    dim(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$p.val),
+    c('nexp' = 3, 'nobs' = 2, 'lat' = 2)
+  )
+  ## values
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset')$p.val),
+    c(0.1674805),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[, , 2]),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4, obs4, dat_dim = 'dataset')$p.val[, , 2]),
+    c(0.30405979, 0.18162950, 0.01675207, 0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  # expect_equal( #############################################
+  #   SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[1],
+  #   SprErr(exp2, obs2)$ratio
+  # )
+  expect_equal(
+    SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[, , 1],
+    SprErr(exp3, obs3, dat_dim = 'dataset')$ratio
+  )
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$ratio),
+    c(0.5831841),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    mean(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$p.val),
+    c(0.1674805),
+    tolerance = 0.0001
+  )
+  
+  # dat4_2: NAs
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$ratio[, , 1]),
+    c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$p.val[, , 1]),
+    c(NA_real_, NA_real_, NA_real_, NA_real_, NA_real_, NA_real_),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', sign = TRUE)$sign[, , 1]),
+    c(NA, NA, NA, NA, NA, NA),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$ratio[, , 2]),
+    c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$p.val[, , 2]),
+    c(0.30405979, 0.18162950, 0.01675207, 0.30405979, 0.18162950, 0.01675207),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$ratio)),
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset')$p.val))
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$ratio[, , 1]),
+    c(0.4648097, 0.6571888, 0.3975804, 0.4648097, 0.6571888, 0.3975804),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$p.val[, , 1]),
+    c(0.04674814, 0.21983481, 0.01350139, 0.04351013, 0.22700305, 0.01124875),
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm =TRUE, sign = TRUE)$sign[, , 1]),
+    c(TRUE, FALSE, TRUE, TRUE, FALSE, TRUE),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$ratio)),
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$p.val))
+  )
+  expect_equal(
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$ratio)),
+    which(is.na(SprErr(exp4_2, obs4_2, dat_dim = 'dataset', na.rm = TRUE)$p.val))
+  )
+  
+})
-- 
GitLab


From c22c2b1e2eb8a0d64d910a174037ff8c7a133e10 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Mon, 15 Jul 2024 16:37:56 +0200
Subject: [PATCH 23/28] removed source() lines in test-SprErr.R

---
 tests/testthat/test-SprErr.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
index a116bde..76bea9d 100644
--- a/tests/testthat/test-SprErr.R
+++ b/tests/testthat/test-SprErr.R
@@ -3,8 +3,7 @@ library(s2dv)
 library(testthat)
 library(multiApply)
 library(ClimProjDiags) # for Subset() [commit 53ba0039 by aho]
-source("R/SprErr.R")
-source("R/Eno.R")
+
 
 ##############################################
 # data
-- 
GitLab


From 37151473fcd5157a5707d5c6d75d45850bd98799 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Mon, 15 Jul 2024 16:55:21 +0200
Subject: [PATCH 24/28] Specified comment in test-SprErr.R

---
 tests/testthat/test-SprErr.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
index 76bea9d..dd3d0ba 100644
--- a/tests/testthat/test-SprErr.R
+++ b/tests/testthat/test-SprErr.R
@@ -2,7 +2,7 @@
 library(s2dv)
 library(testthat)
 library(multiApply)
-library(ClimProjDiags) # for Subset() [commit 53ba0039 by aho]
+library(ClimProjDiags) ## for Subset() [commit 53ba0039 by aho]
 
 
 ##############################################
-- 
GitLab


From f6cad5a720f9c7ce98a69cc7c3a66f5b09bdba76 Mon Sep 17 00:00:00 2001
From: ARIADNA BATALLA FERRES <ariadna.batalla@bsc.es>
Date: Mon, 15 Jul 2024 16:59:59 +0200
Subject: [PATCH 25/28] Edited comments in test-SprErr.R

---
 tests/testthat/test-SprErr.R | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
index dd3d0ba..7cdbe09 100644
--- a/tests/testthat/test-SprErr.R
+++ b/tests/testthat/test-SprErr.R
@@ -2,7 +2,7 @@
 library(s2dv)
 library(testthat)
 library(multiApply)
-library(ClimProjDiags) ## for Subset() [commit 53ba0039 by aho]
+library(ClimProjDiags) ## for Subset() [view commit 53ba0039 by aho]
 
 
 ##############################################
@@ -61,7 +61,7 @@ obs4_2 <- obs4; obs4_2[1, 1:4, 1, 1] <- NA
 # exp and obs (1)
 test_that("1. Input checks", {
   expect_error(
-    SprErr(c(), obs1), ## or: SprErr(exp1, obs1, c()) ?? exp1, c() i c() obs1 - què fa quan li passes algo que no és un array? c()
+    SprErr(c(), obs1), 
     "Parameter 'exp' and 'obs' cannot be NULL." 
   )
   expect_error(
@@ -78,7 +78,7 @@ test_that("1. Input checks", {
   )
   expect_error(
     SprErr(1, obs1), 
-    "Parameter 'exp' and 'obs' must be array with as least two dimensions memb_dim and time_dim." ## si es pot es evitar el paste0 millor
+    "Parameter 'exp' and 'obs' must be array with as least two dimensions memb_dim and time_dim."
   )
   expect_error(
     SprErr(exp1, 1), 
@@ -158,7 +158,7 @@ test_that("2. Output checks: dat1", {
   
   ## element names
   expect_equal(
-    names(SprErr(exp1, obs1)), ## default: pval = TRUE, sign = FALSE
+    names(SprErr(exp1, obs1)),
     c("ratio", "p.val")
   )
   expect_equal(
@@ -523,7 +523,8 @@ test_that("5. Output checks: dat4", {
     c(0.30405979, 0.18162950, 0.01675207, 0.30405979, 0.18162950, 0.01675207),
     tolerance = 0.0001
   )
-  # expect_equal( #############################################
+  ## TODO: remove?
+  # expect_equal( 
   #   SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[1],
   #   SprErr(exp2, obs2)$ratio
   # )
-- 
GitLab


From 5496826f93479c2a23795bd5af68d97dbe6c965e Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Thu, 25 Jul 2024 13:23:29 +0200
Subject: [PATCH 26/28] Clean comments and format of test-SprErr.R

---
 tests/testthat/test-SprErr.R | 53 ++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 26 deletions(-)

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
index 7cdbe09..ca9d7c5 100644
--- a/tests/testthat/test-SprErr.R
+++ b/tests/testthat/test-SprErr.R
@@ -1,8 +1,7 @@
-
 library(s2dv)
 library(testthat)
 library(multiApply)
-library(ClimProjDiags) ## for Subset() [view commit 53ba0039 by aho]
+library(ClimProjDiags)
 
 
 ##############################################
@@ -151,12 +150,12 @@ test_that("1. Input checks", {
   )
 })
 
-######################################################## #
+##############################################
 
 
 test_that("2. Output checks: dat1", {
   
-  ## element names
+  # element names
   expect_equal(
     names(SprErr(exp1, obs1)),
     c("ratio", "p.val")
@@ -169,8 +168,7 @@ test_that("2. Output checks: dat1", {
     names(SprErr(exp1, obs1, pval = TRUE, sign = TRUE)),
     c("ratio", "p.val", "sign")
   )
-  
-  ## dimensions
+  # dimensions
   expect_equal(
     dim(SprErr(exp1, obs1)$ratio),
     c(lat = 2)
@@ -187,8 +185,7 @@ test_that("2. Output checks: dat1", {
     dim(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$sign), 
     c(lat = 2)
   )
-  
-  ## values
+  # values
   expect_equal(
     as.vector(SprErr(exp1, obs1)$ratio),
     c(1.0646692, 0.6657522),
@@ -203,7 +200,7 @@ test_that("2. Output checks: dat1", {
     as.vector(SprErr(exp1, obs1, sign = TRUE)$sign),
     c(FALSE, FALSE)
   )
-  ## pval = FALSE
+  # pval = FALSE
   expect_equal(
     as.vector(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$ratio),
     c(01.0646692, 0.6657522),  
@@ -213,7 +210,7 @@ test_that("2. Output checks: dat1", {
     as.vector(SprErr(exp1, obs1, pval = FALSE, sign = TRUE)$sign),
     c(FALSE, FALSE)  
   )
-  ## na.rm = TRUE
+  # na.rm = TRUE
   expect_equal(
     as.vector(SprErr(exp1, obs1, na.rm = TRUE)$ratio),
     c(1.0646692, 0.6657522),  
@@ -228,7 +225,7 @@ test_that("2. Output checks: dat1", {
     as.vector(SprErr(exp1, obs1, sign = TRUE, na.rm = TRUE)$sign),
     c(FALSE, FALSE)
   )
-  ## alpha
+  # alpha
   expect_equal(
     as.vector(SprErr(exp1, obs1, sign = TRUE, alpha = 0.5)$sign),
     c(FALSE, TRUE)
@@ -284,8 +281,9 @@ test_that("2. Output checks: dat1", {
   
 })
 
-
 ##############################################
+
+
 test_that("3. Output checks: dat2", {
   
   expect_equal(
@@ -308,7 +306,7 @@ test_that("3. Output checks: dat2", {
     dim(SprErr(exp2, obs2, sign = TRUE)$sign),
     NULL
   )
-  ## values
+  # values
   expect_equal(
     as.vector(SprErr(exp2, obs2)$ratio),
     c(0.6866402),
@@ -323,7 +321,7 @@ test_that("3. Output checks: dat2", {
     as.vector(SprErr(exp2, obs2, sign = TRUE)$sign), 
     FALSE 
   )
-  ## sign = TRUE
+  # sign = TRUE
   expect_equal(
     as.vector(SprErr(exp2, obs2, sign = TRUE)$ratio),
     c(0.6866402),
@@ -334,12 +332,12 @@ test_that("3. Output checks: dat2", {
     c(0.2779936),
     tolerance = 0.0001
   )
-  ## alpha
+  # alpha
   expect_equal(
     as.vector(SprErr(exp2, obs2, sign = TRUE, alpha = 0.99)$sign),
     TRUE 
   )
-  ## na.rm = TRUE
+  # na.rm = TRUE
   expect_equal(
     as.vector(SprErr(exp2, obs2, na.rm = TRUE)$ratio),
     c(0.6866402),
@@ -350,11 +348,11 @@ test_that("3. Output checks: dat2", {
     c(0.2779936),
     tolerance = 0.0001
   )
-  expect_equal( # is this expected?
+  expect_equal(
     as.vector(SprErr(exp2, obs2, na.rm = TRUE, sign = TRUE)$sign), 
     FALSE 
   )
-  ## other
+  # other
   expect_equal(
     SprErr(exp2, obs2),
     SprErr(exp2, obs2_1)
@@ -363,6 +361,8 @@ test_that("3. Output checks: dat2", {
 })
 
 ##############################################
+
+
 test_that("4. Output checks: dat3", {
   
   expect_equal(
@@ -373,7 +373,7 @@ test_that("4. Output checks: dat3", {
     dim(SprErr(exp3, obs3, dat_dim = 'dataset')$p.val),
     c('nexp' = 3, 'nobs' = 2)
   )
-  ## values
+  # values
   expect_equal(
     mean(SprErr(exp3, obs3, dat_dim = 'dataset')$ratio),
     c(0.5831841),
@@ -403,7 +403,7 @@ test_that("4. Output checks: dat3", {
     c(0.30405979, 0.18162950, 0.01675207),
     tolerance = 0.0001
   )
-  ## alpha
+  # alpha
   expect_equal(
     as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE, alpha = 0.99)$sign)[1:3],
     c(TRUE, TRUE, TRUE)
@@ -412,7 +412,7 @@ test_that("4. Output checks: dat3", {
     as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', sign = TRUE, alpha = 0.20)$sign)[1:3],
     c(FALSE, TRUE, TRUE)
   )
-  ## na.rm = TRUE
+  # na.rm = TRUE
   expect_equal(
     as.vector(SprErr(exp3, obs3, dat_dim = 'dataset', na.rm = TRUE)$ratio),
     c(0.7006396, 0.6277856, 0.4211269, 0.7006396, 0.6277856, 0.4211269),
@@ -437,7 +437,7 @@ test_that("4. Output checks: dat3", {
     dim(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$p.val),
     c('nexp' = 4, 'nobs' = 3)
   )
-  ## values
+  # values
   expect_equal(
     mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset')$ratio),
     c(1.25586),
@@ -462,7 +462,7 @@ test_that("4. Output checks: dat3", {
     c(0.6927309, 0.7390035, 0.8834023, 0.4421531),
     tolerance = 0.0001
   )
-  ## alpha
+  # alpha
   expect_equal(
     as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE, alpha = 0.99)$sign)[1:3],
     c(TRUE, TRUE, TRUE)
@@ -471,7 +471,7 @@ test_that("4. Output checks: dat3", {
     as.vector(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', sign = TRUE, alpha = 0.70)$sign)[1:4],
     c(TRUE, FALSE, FALSE, TRUE)
   )
-  ## na.rm = TRUE
+  # na.rm = TRUE
   expect_equal(
     mean(SprErr(exp3_2, obs3_2, dat_dim = 'dataset', na.rm = TRUE)$ratio),
     c(1.25586),
@@ -490,8 +490,9 @@ test_that("4. Output checks: dat3", {
   
 })
 
+############################################## 
+
 
-##############################################
 test_that("5. Output checks: dat4", {
   
   expect_equal(
@@ -502,7 +503,7 @@ test_that("5. Output checks: dat4", {
     dim(SprErr(exp4, obs4, dat_dim = 'dataset', sign = TRUE)$p.val),
     c('nexp' = 3, 'nobs' = 2, 'lat' = 2)
   )
-  ## values
+  # values
   expect_equal(
     mean(SprErr(exp4, obs4, dat_dim = 'dataset')$ratio),
     c(0.5831841),
-- 
GitLab


From 1fa20d4e69c95445f367bb6fc64f98f806a43b84 Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Thu, 25 Jul 2024 13:24:52 +0200
Subject: [PATCH 27/28] Remove expect_equal(SprErr(exp4, obs4, dat_dim =
 'dataset')$ratio[1], SprErr(exp2, obs2)$ratio) in test-SprErr.R

---
 tests/testthat/test-SprErr.R | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
index ca9d7c5..bb03574 100644
--- a/tests/testthat/test-SprErr.R
+++ b/tests/testthat/test-SprErr.R
@@ -524,11 +524,6 @@ test_that("5. Output checks: dat4", {
     c(0.30405979, 0.18162950, 0.01675207, 0.30405979, 0.18162950, 0.01675207),
     tolerance = 0.0001
   )
-  ## TODO: remove?
-  # expect_equal( 
-  #   SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[1],
-  #   SprErr(exp2, obs2)$ratio
-  # )
   expect_equal(
     SprErr(exp4, obs4, dat_dim = 'dataset')$ratio[, , 1],
     SprErr(exp3, obs3, dat_dim = 'dataset')$ratio
-- 
GitLab


From 621ed770988dd61307bc3372c9cebc37777ecc5e Mon Sep 17 00:00:00 2001
From: abatalla <ariadna.batalla@bsc.es>
Date: Thu, 25 Jul 2024 15:39:12 +0200
Subject: [PATCH 28/28] Format fix test-SprErr.R

---
 tests/testthat/test-SprErr.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/testthat/test-SprErr.R b/tests/testthat/test-SprErr.R
index bb03574..65ff728 100644
--- a/tests/testthat/test-SprErr.R
+++ b/tests/testthat/test-SprErr.R
@@ -57,8 +57,9 @@ obs4_2 <- obs4; obs4_2[1, 1:4, 1, 1] <- NA
 # tests
 ##############################################
 
-# exp and obs (1)
 test_that("1. Input checks", {
+
+  # exp and obs (1)
   expect_error(
     SprErr(c(), obs1), 
     "Parameter 'exp' and 'obs' cannot be NULL." 
-- 
GitLab