From c060a7455d662918349e7616ab096e4a1a9df296 Mon Sep 17 00:00:00 2001 From: Nicolau Manubens Date: Tue, 7 Feb 2017 18:17:43 +0100 Subject: [PATCH 1/6] Attempt to fix the bug. Vero to test it. --- R/ProbBins.R | 143 ++++++++++++++++++++++++++++----------------------- 1 file changed, 80 insertions(+), 63 deletions(-) diff --git a/R/ProbBins.R b/R/ProbBins.R index c6e4fec1..0f529faa 100644 --- a/R/ProbBins.R +++ b/R/ProbBins.R @@ -1,9 +1,8 @@ -ProbBins <- function(ano, fcyr, thr, quantile=TRUE, posdates = 3, - posdim = 2, compPeriod= "Full period") { +ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, + posdim = 2, compPeriod = "Full period") { # define dimensions #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ nbdim <- length(dim(ano)) - nfcyr<-length(fcyr) if (nbdim < 7){ ano <- Enlarge(ano, 7) @@ -13,78 +12,96 @@ ProbBins <- function(ano, fcyr, thr, quantile=TRUE, posdates = 3, posdim <- setdiff(posdim, posdates) nbpos <- length(posdim) #permute dimensions in ano - ano <- aperm(ano, c(posdates, posdim, setdiff(seq(1,7,1), c(posdates, posdim)))) + if (posdates != 1 || posdim != 2) { + ano <- aperm(ano, c(posdates, posdim, setdiff(seq(1,7,1), c(posdates, posdim)))) + posdates <- 1 + posdim <- 2 + } dimano <- dim(ano) - nsdates=dimano[1] + nsdates <- dimano[1] #calculate the number of elements on posdim dimension in ano - nmemb=1 + nmemb <- 1 if (nbpos > 0){ for (idim in 2:(nbpos+1)){ - nmemb=nmemb*dimano[idim] + nmemb <- nmemb*dimano[idim] } } - - # separate forecast and hindcast - #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - fore <- array(ano[fcyr, , , , , , ], dim = c(nfcyr, - dimano[2:7])) - # the members and startdates are combined in one dimension - sample_fore <- array(fore, dim=c(nfcyr*nmemb, dimano[(nbpos+2):7])) - - if(compPeriod=="Full period"){ - hind <- ano - sample <- array(hind, dim=c(nsdates*nmemb, dimano[(nbpos+2):7])) + + if (fcyr == 'all') { + fcyr <- 1:nsdates } - - if (compPeriod=="Without fcyr"){ - hind <- array(ano[-fcyr, , , , , , ], dim = c(nsdates-nfcyr, + nfcyr <- length(fcyr) + + if (compPeriod == "Cross-validation") { + result <- array(dim = c(nfcyr, length(thr) + 1, 1, nmemb, dimano[(nbpos + 2):nbdim])) + store_indices <- as.list(rep(TRUE, length(dim(result)))) + for (iyr in fcyr) { + store_indices[[1]] <- iyr + result <- do.call("[<-", c(list(x = result), + store_indices, + list(value = ProbBins(ano, iyr, thr, quantile, + posdates, posdim, + "Without fcyr")) + ) + ) + } + return(result) + } else if (compPeriod %in% c("Full period", "Without fcyr")) { + # separate forecast and hindcast + #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + fore <- array(ano[fcyr, , , , , , ], dim = c(nfcyr, dimano[2:7])) - sample <- array(hind, dim=c((nsdates-nfcyr)*nmemb, dimano[(nbpos+2):7])) - } - - #quantiles for each grid point and experiment - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # the members and startdates are combined in one dimension + sample_fore <- array(fore, dim=c(nfcyr*nmemb, dimano[(nbpos+2):7])) + + if(compPeriod == "Full period") { + hind <- ano + sample <- array(hind, dim=c(nsdates*nmemb, dimano[(nbpos+2):7])) + } else if (compPeriod == "Without fcyr") { + hind <- array(ano[-fcyr, , , , , , ], dim = c(nsdates-nfcyr, + dimano[2:7])) + sample <- array(hind, dim=c((nsdates-nfcyr)*nmemb, dimano[(nbpos+2):7])) + } - if (quantile==TRUE){ - qum <- apply(sample, seq(2,7-nbpos,1), FUN=quantile,probs=thr,na.rm=TRUE,names=FALSE,type=8) - }else{ - qum<-array(thr,dim=c(length(thr), dimano[(nbpos+2):7])) - } - - # This function assign the values to a category which is limited by the thresholds - # It provides binary information + #quantiles for each grid point and experiment + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + if (quantile==TRUE){ + qum <- apply(sample, seq(2,7-nbpos,1), FUN=quantile,probs=thr,na.rm=TRUE,names=FALSE,type=8) + }else{ + qum<-array(thr,dim=c(length(thr), dimano[(nbpos+2):7])) + } - counts <- function (dat, nbthr){ - thr <- dat[1:nbthr] - data <- dat[nbthr+1:(length(dat)-nbthr)] - prob <- array(NA, dim=c(nbthr+1, length(dat)-nbthr)) - prob[1,]=1*(data <= thr[1]) - if(nbthr!=1){ - for (ithr in 2:(nbthr)){ - prob[ithr,]=1*((data > thr[ithr - 1]) & (data <= thr[ithr])) + # This function assign the values to a category which is limited by the thresholds + # It provides binary information + + counts <- function (dat, nbthr){ + thr <- dat[1:nbthr] + data <- dat[nbthr+1:(length(dat)-nbthr)] + prob <- array(NA, dim=c(nbthr+1, length(dat)-nbthr)) + prob[1,]=1*(data <= thr[1]) + if(nbthr!=1){ + for (ithr in 2:(nbthr)){ + prob[ithr,]=1*((data > thr[ithr - 1]) & (data <= thr[ithr])) + } } + prob[nbthr+1,]=1*(data > thr[nbthr]) + return(prob) } - prob[nbthr+1,]=1*(data > thr[nbthr]) - return(prob) - } - - # The thresholds and anomalies are combined to use apply - data <- abind(qum, sample_fore, along = 1) - - # PBF:Probabilistic bins of a forecast. - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # This array contains zeros and ones that indicate the category where your forecast is. - - PBF <- array(apply(data, seq(2,7-nbpos,1), FUN=counts, nbthr=length(thr)), - dim=c(length(thr)+1, nfcyr, nmemb, dimano[(nbpos+2):nbdim])) - - return(PBF) - + + # The thresholds and anomalies are combined to use apply + data <- abind(qum, sample_fore, along = 1) + + # PBF:Probabilistic bins of a forecast. + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # This array contains zeros and ones that indicate the category where your forecast is. - if (compPeriod == "cross-validation"){ - for (iyr in 1:fcyr){ - ProbBins(ano,fcyr=iyr,thr=thr,posdates=posdates,posdim=posdim) - } + PBF <- array(apply(data, seq(2,7-nbpos,1), FUN=counts, nbthr=length(thr)), + dim=c(length(thr)+1, nfcyr, nmemb, dimano[(nbpos+2):nbdim])) + + return(PBF) + } else { + stop("Parameter 'compPeriod' must be one of 'Full period', 'Without fcyr' or 'Cross-validation'.") } } -- GitLab From 0aa2f79760e8e127e9c6bcb45e9a0a71dd1cbeb8 Mon Sep 17 00:00:00 2001 From: Nicolau Manubens Date: Tue, 7 Feb 2017 18:43:19 +0100 Subject: [PATCH 2/6] Enhancement. --- R/ProbBins.R | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/R/ProbBins.R b/R/ProbBins.R index 0f529faa..c4b11d60 100644 --- a/R/ProbBins.R +++ b/R/ProbBins.R @@ -33,17 +33,11 @@ ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, nfcyr <- length(fcyr) if (compPeriod == "Cross-validation") { - result <- array(dim = c(nfcyr, length(thr) + 1, 1, nmemb, dimano[(nbpos + 2):nbdim])) - store_indices <- as.list(rep(TRUE, length(dim(result)))) + result <- NULL for (iyr in fcyr) { store_indices[[1]] <- iyr - result <- do.call("[<-", c(list(x = result), - store_indices, - list(value = ProbBins(ano, iyr, thr, quantile, - posdates, posdim, - "Without fcyr")) - ) - ) + result <- abind(result, ProbBins(ano, iyr, thr, quantile, + posdates, posdim, "Without fcyr"), 2) } return(result) } else if (compPeriod %in% c("Full period", "Without fcyr")) { -- GitLab From 7e22d4e40904cd74507abe987b34e2b88d5633bb Mon Sep 17 00:00:00 2001 From: Nicolau Manubens Date: Tue, 7 Feb 2017 18:47:40 +0100 Subject: [PATCH 3/6] Doc enhancements. --- man/ProbBins.Rd | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/man/ProbBins.Rd b/man/ProbBins.Rd index 96ec0576..a1c5ac50 100644 --- a/man/ProbBins.Rd +++ b/man/ProbBins.Rd @@ -7,7 +7,7 @@ Computes probabilistic information of a forecast relative to a threshold or a qu Compute probabilistic bins of a set of forecast years ('fcyr') relative to the forecast climatology over the whole period of anomalies, optionally excluding the selected forecast years ('fcyr') or the forecast year for which the probabilistic bins are being computed (see 'compPeriod'). } \usage{ -ProbBins(ano, fcyr, thr, quantile = TRUE, posdates = 3, posdim = 2, +ProbBins(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, posdim = 2, compPeriod = "Full period") } \arguments{ @@ -16,8 +16,8 @@ Array of anomalies from Ano().\cr Must be of dimension (nexp/nobs, nmemb, nsdates, nleadtime, nlat, nlon) } \item{fcyr}{ -Indices of the forecast years of the anomalies of which to compute the probabilistic bins.\cr -Ex: c(1:5), c(1, 4) or 4. +Indices of the forecast years of the anomalies which to compute the probabilistic bins for, or 'all' to compute the bins for all the years.\cr +Ex: c(1:5), c(1, 4), 4 or 'all'. } \item{thr}{ Values used as thresholds to bin the anomalies. @@ -33,12 +33,12 @@ Position of the dimension in \code{ano} that corresponds to the start dates (def Position of the dimension in \code{ano} which will be combined with 'posdates' to compute the quantiles (default = 2, ensemble members). } \item{compPeriod}{ -Three options: "Full period"/"Without fcyr"/"cross-validation" (The probabilities are computed with the terciles based on ano/ano with all 'fcyr's removed/cross-validation). The default is "Full period". +Three options: "Full period"/"Without fcyr"/"Cross-validation" (The probabilities are computed with the terciles based on ano/ano with all 'fcyr's removed/cross-validation). The default is "Full period". } } \value{ Matrix with probabilistic information and dimensions:\cr -c(length('thr'+1), nfcyr, nmemb/nparam, nmod/nexp/nobs, nltime, nlat, nlon)\cr +c(length('thr') + 1, length(fcyr), nmemb/nparam, nmod/nexp/nobs, nltime, nlat, nlon)\cr The values along the first dimension take values 0 or 1 depending on which of the 'thr'+1 cathegories the forecast/observation at the corresponding grid point, time step, member and starting date belongs to. } \examples{ @@ -74,6 +74,7 @@ PB <- ProbBins(ano_exp, fcyr = 3, thr = c(1/3, 2/3), quantile = TRUE, posdates = \author{ History:\cr 1.0 - 2013 (F.Lienert) - Original code\cr -2.0 - 2014-03 (N. Gonzalez and V.Torralba, \email{veronica.torralba@ic3.cat}) - Debugging +2.0 - 2014-03 (N. Gonzalez and V. Torralba, \email{veronica.torralba at bsc.es}) - Debugging +2.1 - 2017-02 (V. Torralba and N. Manubens, \email{veronica.torralba at bsc.es}) - Fix bug with cross-validation } \keyword{datagen} -- GitLab From fbedcd25dd1dd75e96b890e64b839e26147b3938 Mon Sep 17 00:00:00 2001 From: Nicolau Manubens Date: Tue, 7 Feb 2017 18:50:26 +0100 Subject: [PATCH 4/6] Small fix. --- R/ProbBins.R | 1 - 1 file changed, 1 deletion(-) diff --git a/R/ProbBins.R b/R/ProbBins.R index c4b11d60..1651a01f 100644 --- a/R/ProbBins.R +++ b/R/ProbBins.R @@ -35,7 +35,6 @@ ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, if (compPeriod == "Cross-validation") { result <- NULL for (iyr in fcyr) { - store_indices[[1]] <- iyr result <- abind(result, ProbBins(ano, iyr, thr, quantile, posdates, posdim, "Without fcyr"), 2) } -- GitLab From 82a6087a67de1017b6048341da4438f4f19116c2 Mon Sep 17 00:00:00 2001 From: Nicolau Manubens Date: Tue, 7 Feb 2017 19:08:35 +0100 Subject: [PATCH 5/6] Small fix. --- R/ProbBins.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/ProbBins.R b/R/ProbBins.R index 1651a01f..2e2495fb 100644 --- a/R/ProbBins.R +++ b/R/ProbBins.R @@ -36,7 +36,8 @@ ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, result <- NULL for (iyr in fcyr) { result <- abind(result, ProbBins(ano, iyr, thr, quantile, - posdates, posdim, "Without fcyr"), 2) + posdates, posdim, "Without fcyr"), + along = 2) } return(result) } else if (compPeriod %in% c("Full period", "Without fcyr")) { -- GitLab From 4af7600172005d16ca7e2f383daf0b0097365548 Mon Sep 17 00:00:00 2001 From: Nicolau Manubens Date: Wed, 8 Feb 2017 11:46:52 +0100 Subject: [PATCH 6/6] Now returning dimension names. --- R/ProbBins.R | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/R/ProbBins.R b/R/ProbBins.R index 2e2495fb..4aee0817 100644 --- a/R/ProbBins.R +++ b/R/ProbBins.R @@ -13,7 +13,10 @@ ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, nbpos <- length(posdim) #permute dimensions in ano if (posdates != 1 || posdim != 2) { - ano <- aperm(ano, c(posdates, posdim, setdiff(seq(1,7,1), c(posdates, posdim)))) + dimnames_backup <- names(dim(ano)) + perm <- c(posdates, posdim, (1:7)[-c(posdates, posdim)]) + ano <- aperm(ano, perm) + names(dim(ano)) <- dimnames_backup[perm] posdates <- 1 posdim <- 2 } @@ -35,9 +38,16 @@ ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, if (compPeriod == "Cross-validation") { result <- NULL for (iyr in fcyr) { - result <- abind(result, ProbBins(ano, iyr, thr, quantile, - posdates, posdim, "Without fcyr"), - along = 2) + if (is.null(result)) { + result <- ProbBins(ano, iyr, thr, quantile, + posdates, posdim, "Without fcyr") + } else { + dimnames_backup <- names(dim(result)) + result <- abind(result, ProbBins(ano, iyr, thr, quantile, + posdates, posdim, "Without fcyr"), + along = 2) + names(dim(result)) <- dimnames_backup + } } return(result) } else if (compPeriod %in% c("Full period", "Without fcyr")) { @@ -94,6 +104,7 @@ ProbBins <- function(ano, fcyr = 'all', thr, quantile = TRUE, posdates = 3, PBF <- array(apply(data, seq(2,7-nbpos,1), FUN=counts, nbthr=length(thr)), dim=c(length(thr)+1, nfcyr, nmemb, dimano[(nbpos+2):nbdim])) + names(dim(PBF)) <- c('bin', 'sdate', 'member', names(dim(ano))[(nbpos+2):nbdim]) return(PBF) } else { stop("Parameter 'compPeriod' must be one of 'Full period', 'Without fcyr' or 'Cross-validation'.") -- GitLab