diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml new file mode 100644 index 0000000000000000000000000000000000000000..b1cf9976e0a017b4076e024242f9c1a9a7e61487 --- /dev/null +++ b/.gitlab-ci.yml @@ -0,0 +1,10 @@ +stages: + - build +build: + stage: build + script: + - module load R +# - module load CDO + - R CMD build --resave-data . + - R CMD check --as-cran --no-manual --run-donttest s2dverification_*.tar.gz + - R -e 'covr::package_coverage()' diff --git a/DESCRIPTION b/DESCRIPTION index f77ce597e30b6bd72d90a0934d1b2e88a221a59f..99c08ff66547282a64f34fe3ea1cfedc161aa28a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -26,7 +26,8 @@ Authors@R: c( person("Eleftheria", "Exarchou", , "eleftheria.exarchou@bsc.es", role = "ctb"), person("Ruben", "Cruz", , "ruben.cruzgarcia@bsc.es", role = "ctb"), person("Isabel", "Andreu-Burillo", , "isabel.andreu.burillo@ic3.cat", role = "ctb"), - person("Ramiro", "Saurral", , "ramiro.saurral@ic3.cat", role = "ctb")) + person("Ramiro", "Saurral", , "ramiro.saurral@ic3.cat", role = "ctb"), + person("An-Chi", "Ho", , "an.ho@bsc.es", role = "ctb")) Description: Set of tools to verify forecasts through the computation of typical prediction scores against one or more observational datasets or reanalyses (a reanalysis being a physical extrapolation of observations that relies on the equations from a model, not a pure observational dataset). Intended for seasonal to decadal climate forecasts although can be useful to verify other kinds of forecasts. The package can be helpful in climate sciences for other purposes than forecasting. Depends: maps, @@ -44,7 +45,8 @@ Imports: plyr, SpecsVerification (>= 0.5.0) Suggests: - easyVerification + easyVerification, + testthat License: LGPL-3 URL: https://earth.bsc.es/gitlab/es/s2dverification/wikis/home BugReports: https://earth.bsc.es/gitlab/es/s2dverification/issues diff --git a/R/Composite.R b/R/Composite.R index b5b54914618c441a4c4b59dd31a9a1dd42c7012e..c8294fd26c77ec2cdf7421d990b3e25f892a7cf9 100644 --- a/R/Composite.R +++ b/R/Composite.R @@ -1,39 +1,54 @@ -Composite <- function(var, occ, lag=0, eno=FALSE, fileout=NULL) { - - if ( dim(var)[3]!=length(occ) ) { stop("temporal dimension of var is not equal to length of occ") } +Composite <- function(var, occ, lag = 0, eno = FALSE, fileout = NULL) { + if ( dim(var)[3] != length(occ) ) { + stop("temporal dimension of var is not equal to length of occ.") + } K <- max(occ) composite <- array(dim = c(dim(var)[1:2], K)) tvalue <- array(dim = dim(var)[1:2]) dof <- array(dim = dim(var)[1:2]) pvalue <- array(dim = c(dim(var)[1:2], K)) - if ( eno==TRUE ) { n_tot <- Eno(var, posdim=3) } - else { n_tot <- length(occ) } - - mean_tot <- Mean1Dim(var, posdim=3, narm=TRUE) - stdv_tot <- apply(var, c(1,2), sd, na.rm=TRUE) - - for (k in 1:K) { - indices <- which(occ==k)+lag - - toberemoved=which(0>indices|indices>dim(var)[3]) - if ( length(toberemoved) > 0 ) { indices=indices[-toberemoved] } - - if ( eno==TRUE ) { n_k <- Eno(var[,,indices], posdim=3) } - else { n_k <- length(indices) } - - composite[,,k] <- Mean1Dim(var[,,indices], posdim=3, narm=TRUE) - stdv_k <- apply(var[,,indices], c(1,2), sd, na.rm=TRUE) + if (eno == TRUE) { + n_tot <- Eno(var, posdim = 3) + } else { + n_tot <- length(occ) + } + mean_tot <- Mean1Dim(var, posdim = 3, narm = TRUE) + stdv_tot <- apply(var, c(1, 2), sd, na.rm = TRUE) + + for (k in 1 : K) { + + indices <- which(occ == k) + lag + toberemoved <- which(0 > indices | indices > dim(var)[3]) + + if (length(toberemoved) > 0) { + indices <- indices[-toberemoved] + } + if (eno == TRUE) { + n_k <- Eno(var[, , indices], posdim = 3) + } else { + n_k <- length(indices) + } + if (length(indices) == 1) { + composite[, , k] <- var[, , indices] + warning(paste("Composite", k, "has length 1 and pvalue is NA.")) + } else { + composite[, , k] <- Mean1Dim(var[, , indices], posdim = 3, narm = TRUE) + } + stdv_k <- apply(var[, , indices], c(1, 2), sd, na.rm = TRUE) - tvalue[,] <- (mean_tot - composite[,,k])/sqrt(stdv_tot^2/n_tot + stdv_k^2/n_k) - dof[,] <- (stdv_tot^2/n_tot + stdv_k^2/n_k)^2/((stdv_tot^2/n_tot)^2/(n_tot - 1) + (stdv_k^2/n_k)^2/(n_k - 1)) - pvalue[,,k] <- 2*pt(-abs(tvalue[,]), df=dof[,]) + tvalue <- (mean_tot - composite[, , k]) / + sqrt(stdv_tot ^ 2 / n_tot + stdv_k ^ 2 / n_k) + dof <- (stdv_tot ^ 2 / n_tot + stdv_k ^ 2 / n_k) ^ 2 / + ((stdv_tot ^ 2 / n_tot) ^ 2 / (n_tot - 1) + + (stdv_k ^ 2 / n_k) ^ 2 / (n_k - 1)) + pvalue[, , k] <- 2 * pt(-abs(tvalue), df = dof) } - if ( is.null(fileout)==FALSE ) { - output <- list(composite=composite, pvalue=pvalue) - save(output,file=paste(fileout,'.sav',sep='')) + if (is.null(fileout) == FALSE) { + output <- list(composite = composite, pvalue = pvalue) + save(output, file = paste(fileout, '.sav', sep = '')) } invisible(list(composite = composite, pvalue = pvalue)) diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000000000000000000000000000000000000..19e87e1d0acf31320fbbbe91b14573bd1ba0d99d --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(s2dverification) + +test_check("s2dverification") diff --git a/tests/testthat/test-Composite.R b/tests/testthat/test-Composite.R new file mode 100644 index 0000000000000000000000000000000000000000..42c003a00d87f2323a14d72ddad5a48b7fb4c521 --- /dev/null +++ b/tests/testthat/test-Composite.R @@ -0,0 +1,34 @@ +context("Generic tests") +test_that("Sanity checks", { + + expect_error( + Composite(var = array(1:20, dim = c(2, 5, 2)), c(1, 1, 0)), + "temporal dimension of var is not equal to length of occ.") + + expect_warning( + Composite(var = array(1:40, dim = c(2, 5, 4)), c(1, 2, 2, 2)), + "Composite 1 has length 1 and pvalue is NA.") + + var <- array(rep(c(1, 3, 2, 1, 2), 8), dim = c(x = 2, y = 4, time = 5)) + occ <- c(1, 2, 2, 2, 1) + output <- c(x = 2, y = 4, 2) #dim(asd$composite) + expect_equal( + dim(Composite(var, occ)$composite), + output + ) + output <- c(1.5, 2.0, 2.5, 2.0) + expect_equal( + Composite(var, occ)$composite[1, , 1], + output + ) + + var <- array(rep(c(1, 3, 2, 1, 2), 8), dim = c(x = 2, y = 4, time = 5)) + occ <- c(1, 2, 2, 3, 3) + output <- array(as.numeric(rep(NA, 8)), dim = c(2, 4)) + expect_equal( + Composite(var, occ)$pvalue[, , 1], + output + ) + +}) +