Commit 2029e59f authored by aho's avatar aho
Browse files

Merge branch 'develop-BrierScore' into 'master'

Develop BrierScore and UltimateBrier

See merge request !60
parents 4afa1dc5 41a66d62
Pipeline #5733 passed with stage
in 10 minutes and 28 seconds
......@@ -37,6 +37,7 @@ Imports:
ncdf4,
NbClust,
multiApply (>= 2.1.1),
SpecsVerification (>= 0.5.0),
easyNCDF
Suggests:
easyVerification,
......
......@@ -5,6 +5,7 @@ export(AMV)
export(AnimateMap)
export(Ano)
export(Ano_CrossValid)
export(BrierScore)
export(CDORemap)
export(Clim)
export(Cluster)
......@@ -60,10 +61,12 @@ export(Spectrum)
export(TPI)
export(ToyModel)
export(Trend)
export(UltimateBrier)
export(clim.colors)
export(clim.palette)
import(GEOmap)
import(NbClust)
import(SpecsVerification)
import(bigmemory)
import(geomapdata)
import(graphics)
......
#'Compute Brier score, its decomposition, and Brier skill score
#'
#'Compute the Brier score (BS) and the components of its standard decompostion
#'with the two within-bin components described in Stephenson et al., (2008). It
#'also returns the bias-corrected decomposition of the BS (Ferro and Fricker,
#'2012). BSS has the climatology as the reference forecast.
#'
#'@param exp A vector or a numeric array with named dimensions. It should be
#' the predicted probabilities which are within the range [0, 1] if memb_dim
#' doesn't exist. If it has memb_dim, the value should be 0 or 1, and the
#' predicted probabilities will be computed by ensemble mean. The dimensions
#' must at least have 'time_dim'.
#' range [0, 1].
#'@param obs A numeric array with named dimensions of the binary observations
#' (0 or 1). The dimension must be the same as 'exp' except memb_dim, which is
#' optional. If it has 'memb_dim', then the length must be 1. The length of
#' 'dat_dim' can be different from 'exp' if it has.
#'@param thresholds A numeric vector used to bin the forecasts. The default
#' value is \code{seq(0.1, 0.9, 0.1)}, which means that the bins are
#' \code{[0, 0.1), [0.1, 0.2), ... [0.9, 1]}.
#'@param time_dim A character string indicating the name of dimension along
#' which Brier score is computed. The default value is 'sdate'.
#'@param dat_dim A character string indicating the name of dataset dimension in
#' 'exp' and 'obs'. The length of this dimension can be different between
#' 'exp' and 'obs'. The default value is NULL.
#'@param memb_dim A character string of the name of the member dimension in
#' 'exp' (and 'obs', optional). The function will do the ensemble mean
#' over this dimension. If there is no member dimension, set NULL. The default
#' value is NULL.
#'@param ncores An integer indicating the number of cores to use for parallel
#' computation. The default value is NULL.
#'
#'@return
#'A list that contains:
#'\item{$rel}{standard reliability}
#'\item{$res}{standard resolution}
#'\item{$unc}{standard uncertainty}
#'\item{$bs}{Brier score}
#'\item{$bs_check_res}{rel - res + unc}
#'\item{$bss_res}{res - rel / unc}
#'\item{$gres}{generalized resolution}
#'\item{$bs_check_gres}{rel - gres + unc}
#'\item{$bss_gres}{gres - rel / unc}
#'\item{$rel_bias_corrected}{bias - corrected rel}
#'\item{$gres_bias_corrected}{bias - corrected gres}
#'\item{$unc_bias_corrected}{bias - corrected unc}
#'\item{$bss_bias_corrected}{gres_bias_corrected - rel_bias_corrected / unc_bias_corrected}
#'\item{$nk}{number of forecast in each bin}
#'\item{$fkbar}{average probability of each bin}
#'\item{$okbar}{relative frequency that the observed event occurred}
#'The data type and dimensions of the items depend on if the input 'exp' and
#''obs' are:\cr
#'(a) Vectors\cr
#'(b) Arrays with 'dat_dim' specified\cr
#'(c) Arrays with no 'dat_dim' specified\cr
#'Items 'rel', 'res', 'unc', 'bs', 'bs_check_res', 'bss_res', 'gres',
#''bs_check_gres', 'bss_gres', 'rel_bias_corrected', 'gres_bias_corrected',
#''unc_bias_corrected', and 'bss_bias_corrected' are (a) a number (b) an array
#'with dimensions c(nexp, nobs, all the rest dimensions in 'exp' and 'obs'
#'expect 'time_dim' and 'memb_dim') (c) an array with dimensions of
#''exp' and 'obs' except 'time_dim' and 'memb_dim'\cr
#'Items 'nk', 'fkbar', and 'okbar' are (a) a vector of length of bin number
#'determined by 'threshold' (b) an array with dimensions c(nexp, nobs,
#'no. of bins, all the rest dimensions in 'exp' and 'obs' expect 'time_dim' and
#''memb_dim') (c) an array with dimensions c(no. of bin, all the rest dimensions
#'in 'exp' and 'obs' expect 'time_dim' and 'memb_dim')
#'
#'@references
#'Wilks (2006) Statistical Methods in the Atmospheric Sciences.\cr
#'Stephenson et al. (2008). Two extra components in the Brier score decomposition.
#' Weather and Forecasting, 23: 752-757.\cr
#'Ferro and Fricker (2012). A bias-corrected decomposition of the BS.
#' Quarterly Journal of the Royal Meteorological Society, DOI: 10.1002/qj.1924.
#'
#'@examples
#'# Inputs are vectors
#'exp <- runif(10)
#'obs <- round(exp)
#'x <- BrierScore(exp, obs)
#'
#'# Inputs are arrays
#'example(Load)
#'clim <- Clim(sampleData$mod, sampleData$obs)
#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
#'bins_ano_exp <- ProbBins(ano_exp, thr = c(1/3, 2/3))
#'bins_ano_obs <- ProbBins(ano_obs, thr = c(1/3, 2/3))
#'res <- BrierScore(bins_ano_exp, MeanDims(bins_ano_obs, 'member'), memb_dim = 'member')
#'
#'@import multiApply
#'@export
BrierScore <- function(exp, obs, thresholds = seq(0.1, 0.9, 0.1), time_dim = 'sdate',
dat_dim = NULL, memb_dim = NULL, 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 vector or a numeric array.")
}
if (is.null(dim(exp))) { #is vector
dim(exp) <- c(length(exp))
names(dim(exp)) <- time_dim
}
if (is.null(dim(obs))) { #is vector
dim(obs) <- c(length(obs))
names(dim(obs)) <- 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.")
}
if (any(!obs %in% c(0, 1))) {
stop("Parameter 'obs' must be binary events (0 or 1).")
}
## thresholds
if (!is.numeric(thresholds) | !is.vector(thresholds)) {
stop("Parameter 'thresholds' must be a numeric vector.")
}
if (any(thresholds <= 0 | thresholds >= 1)) {
stop("Parameter 'thresholds' must be between 0 and 1 as the bin-breaks.")
}
## 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' and 'obs' dimension.")
}
## 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' and 'obs' dimension.")
}
}
## memb_dim
if (!is.null(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))) {
stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
}
if (memb_dim %in% names(dim(obs))) {
if (dim(obs)[memb_dim] != 1) {
stop("The length of parameter 'memb_dim' in 'obs' must be 1.")
}
}
}
## exp and obs (2)
if (is.null(memb_dim)) {
if (max(exp) > 1 | min(exp) < 0) {
stop("Parameter 'exp' must be within [0, 1] range.")
}
} else {
if (any(!exp %in% c(0, 1))) {
stop("Parameter 'exp' must be 0 or 1 if it has memb_dim.")
}
}
name_exp <- sort(names(dim(exp)))
name_obs <- sort(names(dim(obs)))
if (!is.null(memb_dim)) {
name_exp <- name_exp[-which(name_exp == memb_dim)]
if (memb_dim %in% name_obs) {
name_obs <- name_obs[-which(name_obs == memb_dim)]
}
}
if (!is.null(dat_dim)) {
name_exp <- name_exp[-which(name_exp == dat_dim)]
name_obs <- name_obs[-which(name_obs == dat_dim)]
}
if (any(name_exp != name_obs)) {
stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
"of all the dimensions expect 'dat_dim' and 'memb_dim'."))
}
if (!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
"of all the dimensions expect 'dat_dim' and 'memb_dim'."))
}
## 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 Brier score
## ensemble mean
if (!is.null(memb_dim)) {
exp <- MeanDims(exp, memb_dim)
if (memb_dim %in% names(dim(obs))) {
obs <- MeanDims(obs, memb_dim)
}
}
if (is.null(dat_dim)) {
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim),
c(time_dim)),
fun = .BrierScore,
thresholds = thresholds,
ncores = ncores)
} else {
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim, dat_dim),
c(time_dim, dat_dim)),
fun = .BrierScore,
thresholds = thresholds,
ncores = ncores)
}
return(res)
}
.BrierScore <- function(exp, obs, thresholds = seq(0.1, 0.9, 0.1)) {
# exp: [sdate] or [sdate, nexp]
# obs: [sdate] or [sdate, nobs]
if (length(dim(exp)) == 2) {
nexp <- as.numeric(dim(exp)[2])
nobs <- as.numeric(dim(obs)[2])
exp_ori <- exp
obs_ori <- obs
# Create empty arrays
arr_rel <- arr_res <- arr_unc <- arr_bs <- arr_bs_check_res <- arr_bss_res <-
arr_gres <- arr_bs_check_gres <- arr_bss_gres <- arr_rel_bias_corrected <-
arr_gres_bias_corrected <- arr_unc_bias_corrected <- arr_bss_bias_corrected <-
array(dim = c(nexp = nexp, nobs = nobs))
arr_nk <- arr_fkbar <- arr_okbar <-
array(dim = c(nexp = nexp, nobs = nobs, bin = length(thresholds) + 1))
} else {
nexp <- 1
nobs <- 1
}
for (n_exp in 1:nexp) {
for (n_obs in 1:nobs) {
if (exists('exp_ori')) {
exp <- exp_ori[, n_exp]
obs <- obs_ori[, n_obs]
}
n <- length(exp)
nbins <- length(thresholds) + 1 # Number of bins
bins <- vector('list', nbins) #as.list(paste("bin", 1:nbins, sep = ""))
for (i in 1:nbins) {
if (i == 1) {
bins[[i]] <- list(which(exp >= 0 & exp < thresholds[i]))
} else if (i == nbins) {
bins[[i]] <- list(which(exp >= thresholds[i - 1] & exp <= 1))
} else {
bins[[i]] <- list(which(exp >= thresholds[i - 1] & exp < thresholds[i]))
}
}
fkbar <- okbar <- nk <- array(0, dim = nbins)
for (i in 1:nbins) {
nk[i] <- length(bins[[i]][[1]])
fkbar[i] <- sum(exp[bins[[i]][[1]]]) / nk[i]
okbar[i] <- sum(obs[bins[[i]][[1]]]) / nk[i]
}
#-----in old .BrierScore()---------
# fkbar[fkbar == Inf] <- 0
# okbar[is.nan(okbar)] <- 0
#----------------------------------
obar <- sum(obs) / length(obs)
relsum <- ressum <- term1 <- term2 <- 0
for (i in 1:nbins) {
if (nk[i] > 0) {
relsum <- relsum + nk[i] * (fkbar[i] - okbar[i])^2
ressum <- ressum + nk[i] * (okbar[i] - obar)^2
for (j in 1:nk[i]) {
term1 <- term1 + (exp[bins[[i]][[1]][j]] - fkbar[i])^2
term2 <- term2 + (exp[bins[[i]][[1]][j]] - fkbar[i]) * (obs[bins[[i]][[1]][j]] - okbar[i])
}
}
}
rel <- relsum / n
res <- ressum / n
unc <- obar * (1 - obar)
bs <- sum((exp - obs)^2) / n
bs_check_res <- rel - res + unc
bss_res <- (res - rel) / unc
gres <- res - term1 * (1 / n) + term2 * (2 / n) # Generalized resolution
bs_check_gres <- rel - gres + unc # BS using GRES
bss_gres <- (gres - rel) / unc # BSS using GRES
# Estimating the bias-corrected components of the BS
term3 <- array(0, nbins)
for (i in 1:nbins) {
term3[i] <- (nk[i] / (nk[i] - 1)) * okbar[i] * (1 - okbar[i])
}
term_a <- sum(term3, na.rm = T) / n
term_b <- (obar * (1 - obar)) / (n - 1)
rel_bias_corrected <- rel - term_a
gres_bias_corrected <- gres - term_a + term_b
if (rel_bias_corrected < 0 || gres_bias_corrected < 0) {
rel_bias_corrected2 <- max(rel_bias_corrected, rel_bias_corrected - gres_bias_corrected, 0)
gres_bias_corrected2 <- max(gres_bias_corrected, gres_bias_corrected - rel_bias_corrected, 0)
rel_bias_corrected <- rel_bias_corrected2
gres_bias_corrected <- gres_bias_corrected2
}
unc_bias_corrected <- unc + term_b
bss_bias_corrected <- (gres_bias_corrected - rel_bias_corrected) / unc_bias_corrected
#if (round(bs, 8) == round(bs_check_gres, 8) & round(bs_check_gres, 8) == round((rel_bias_corrected - gres_bias_corrected + unc_bias_corrected), 8)) {
# cat("No error found \ n")
# cat("BS = REL - GRES + UNC = REL_lessbias - GRES_lessbias + UNC_lessbias \ n")
#}
# Add name for nk, fkbar, okbar
names(dim(nk)) <- 'bin'
names(dim(fkbar)) <- 'bin'
names(dim(okbar)) <- 'bin'
if (exists('exp_ori')) {
arr_rel[n_exp, n_obs] <- rel
arr_res[n_exp, n_obs] <- res
arr_unc[n_exp, n_obs] <- unc
arr_bs[n_exp, n_obs] <- bs
arr_bs_check_res[n_exp, n_obs] <- bs_check_res
arr_bss_res[n_exp, n_obs] <- bss_res
arr_gres[n_exp, n_obs] <- gres
arr_bs_check_gres[n_exp, n_obs] <- bs_check_gres
arr_bss_gres[n_exp, n_obs] <- bss_gres
arr_rel_bias_corrected[n_exp, n_obs] <- rel_bias_corrected
arr_gres_bias_corrected[n_exp, n_obs] <- gres_bias_corrected
arr_unc_bias_corrected[n_exp, n_obs] <- unc_bias_corrected
arr_bss_bias_corrected[n_exp, n_obs] <- bss_bias_corrected
arr_nk[n_exp, n_obs, ] <- nk
arr_fkbar[n_exp, n_obs, ] <- fkbar
arr_okbar[n_exp, n_obs, ] <- okbar
}
}
}
if (exists('exp_ori')) {
res_list <- list(rel = arr_rel, res = arr_res, unc = arr_unc, bs = arr_bs,
bs_check_res = arr_bs_check_res, bss_res = arr_bss_res,
gres = arr_gres, bs_check_gres = arr_bs_check_gres,
bss_gres = arr_bss_gres, rel_bias_corrected = arr_rel_bias_corrected,
gres_bias_corrected = arr_gres_bias_corrected,
unc_bias_corrected = arr_unc_bias_corrected,
bss_bias_corrected = arr_bss_bias_corrected, nk = arr_nk,
fkbar = arr_fkbar, okbar = arr_okbar) #bins = list(bins),
} else {
res_list <- list(rel = rel, res = res, unc = unc, bs = bs, bs_check_res = bs_check_res,
bss_res = bss_res, gres = gres, bs_check_gres = bs_check_gres,
bss_gres = bss_gres, rel_bias_corrected = rel_bias_corrected,
gres_bias_corrected = gres_bias_corrected,
unc_bias_corrected = unc_bias_corrected,
bss_bias_corrected = bss_bias_corrected, nk = nk, fkbar = fkbar,
okbar = okbar) #bins = list(bins),
}
return(invisible(res_list))
}
#'Compute Brier scores
#'
#'Interface to compute probabilistic scores (Brier Score, Brier Skill Score)
#'from the forecast and observational data anomalies. It provides six types
#'to choose.
#'
#'@param exp A numeric array of forecast anomalies with named dimensions that
#' at least include 'dat_dim', 'memb_dim', and 'time_dim'. It can be provided
#' by \code{Ano()}.
#'@param obs A numeric array of observational reference anomalies with named
#' dimensions that at least include 'dat_dim' and 'time_dim'. If it has
#' 'memb_dim', the length must be 1. The dimensions should be consistent with
#' 'exp' except 'dat_dim' and 'memb_dim'. It can be provided by \code{Ano()}.
#'@param dat_dim A character string indicating the name of the dataset
#' dimension in 'exp' and 'obs'. The default value is 'dataset'.
#'@param memb_dim A character string indicating the name of the member
#' dimension in 'exp' (and 'obs') for ensemble mean calculation. The default
#' value is 'member'.
#'@param time_dim A character string indicating the dimension along which to
#' compute the probabilistic scores. The default value is 'sdate'.
#'@param quantile A logical value to decide whether a quantile (TRUE) or a
#' threshold (FALSE) is used to estimate the forecast and observed
#' probabilities. If 'type' is 'FairEnsembleBS' or 'FairEnsembleBSS', it must
#' be TRUE. The default value is TRUE.
#'@param thr A numeric vector to be used in probability calculation (for 'BS',
#' 'FairStartDatesBS', 'BSS', and 'FairStartDatesBSS') and binary event
#' judgement (for 'FairEnsembleBS' and 'FairEnsembleBSS'). It is as
#' quantiles if 'quantile' is TRUE or as thresholds if 'quantile' is FALSE.
#' The default value is \code{c(0.05, 0.95)} for 'quantile = TRUE'.
#'@param type A character string of the desired score type. It can be the
#' following values:
#'\itemize{
#' \item{'BS': Simple Brier Score. Use SpecsVerification::BrierDecomp inside.}
#' \item{'FairEnsembleBS': Corrected Brier Score computed across ensemble
#' members. Use SpecsVerification::FairBrier inside.}
#' \item{'FairStartDatesBS': Corrected Brier Score computed across starting
#' dates. Use s2dv:::.BrierScore inside.}
#' \item{'BSS': Simple Brier Skill Score. Use s2dv:::.BrierScore inside.}
#' \item{'FairEnsembleBSS': Corrected Brier Skill Score computed across
#' ensemble members. Use SpecsVerification::FairBrierSs inside.}
#' \item{'FairStartDatesBSS': Corrected Brier Skill Score computed across
#' starting dates. Use s2dv:::.BrierScore inside.}
#'}
#' The default value is 'BS'.
#'@param decomposition A logical value to determine whether the decomposition
#' of the Brier Score should be provided (TRUE) or not (FALSE). It is only
#' used when 'type' is 'BS' or 'FairStartDatesBS'. 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
#'If 'type' is 'BS' or 'FairStartDatesBS' and 'decomposition' is TRUE, the
#'output is a list of 4 arrays (see details below.) In other cases, the output
#'is an array of Brier scores or Brier skill scores. All the arrays have the
#'same dimensions:
#'c(nexp, nobs, no. of bins, the rest dimensions of 'exp' except 'time_dim' and
#''memb_dim'). 'nexp' and 'nobs' is the length of dataset dimension in 'exp'
#'and 'obs' respectively.\cr
#'The list of 4 includes:
#' \itemize{
#' \item{$bs: Brier Score}
#' \item{$rel: Reliability component}
#' \item{$res: Resolution component}
#' \item{$unc: Uncertainty component}
#' }
#'
#'@examples
#' \dontshow{
#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
#' c('observation'), startDates,
#' leadtimemin = 1,
#' leadtimemax = 4,
#' output = 'lonlat',
#' latmin = 27, latmax = 48,
#' lonmin = -12, lonmax = 40)
#' }
#'sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
#'sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2)
#'clim <- Clim(sampleData$mod, sampleData$obs)
#'exp <- Ano(sampleData$mod, clim$clim_exp)
#'obs <- Ano(sampleData$obs, clim$clim_obs)
#'bs <- UltimateBrier(exp, obs)
#'bss <- UltimateBrier(exp, obs, type = 'BSS')
#'
#'@import SpecsVerification plyr multiApply
#'@export
UltimateBrier <- function(exp, obs, dat_dim = 'dataset', memb_dim = 'member', time_dim = 'sdate',
quantile = TRUE, thr = c(5/100, 95/100), type = 'BS',
decomposition = TRUE, 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 vector or a numeric array.")
}
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.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))) {
stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
}
if (!memb_dim %in% names(dim(obs))) {
# Insert memb_dim into obs for the ease of later calculation
obs <- InsertDim(obs, posdim = 2, lendim = 1, name = memb_dim)
} else if (dim(obs)[memb_dim] != 1) {
stop("The length of parameter 'memb_dim' in 'obs' must be 1.")
}
## 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)))
name_exp <- name_exp[-which(name_exp == memb_dim)]
name_obs <- name_obs[-which(name_obs == memb_dim)]
name_exp <- name_exp[-which(name_exp == dat_dim)]
name_obs <- name_obs[-which(name_obs == dat_dim)]
if (any(name_exp != name_obs)) {
stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
"of all the dimensions expect 'dat_dim' and 'memb_dim'."))
}
if (!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
stop(paste0("Parameter 'exp' and 'obs' must have the same names and lengths ",
"of all the dimensions expect 'dat_dim' and 'memb_dim'."))
}
## quantile
if (!is.logical(quantile) | length(quantile) > 1) {
stop("Parameter 'quantile' must be one logical value.")
}
## thr
if (!is.numeric(thr) | !is.vector(thr)) {
stop("Parameter 'thr' must be a numeric vector.")
}
if (quantile) {
if (!all(thr < 1 & thr > 0)) {
stop("Parameter 'thr' must be between 0 and 1 when quantile is TRUE.")
}
}
if (!quantile & (type %in% c('FairEnsembleBSS', 'FairEnsembleBS'))) {
stop("Parameter 'quantile' must be TRUE if 'type' is 'FairEnsembleBSS' or 'FairEnsembleBS'.")
}
## type
if (!(type %in% c("BS", "BSS", "FairEnsembleBS", "FairEnsembleBSS", "FairStartDatesBS", "FairStartDatesBSS"))) {
stop("Parameter 'type' must be one of 'BS', 'BSS', 'FairEnsembleBS', 'FairEnsembleBSS', 'FairStartDatesBS' or 'FairStartDatesBSS'.")
}
## decomposition
if (!is.logical(decomposition) | length(decomposition) > 1) {
stop("Parameter 'decomposition' 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 UltimateBrier
if (type %in% c('FairEnsembleBSS', 'FairEnsembleBS')) {
res <- Apply(list(exp, obs),
target_dims = list(c(time_dim, dat_dim, memb_dim),