From e8a67417da0b8410f82e25c8d4927b9effdf038e Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 26 Oct 2022 15:53:42 +0200 Subject: [PATCH 1/4] First attempt to include window dimension --- R/CST_BiasCorrection.R | 82 +++++++++++++++++++++++++++------------ man/BiasCorrection.Rd | 4 ++ man/CST_BiasCorrection.Rd | 4 ++ 3 files changed, 66 insertions(+), 24 deletions(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 20b51082..f412ecde 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -19,6 +19,8 @@ #' dimension. By default, it is set to 'member'. #'@param sdate_dim A character string indicating the name of the start date #' dimension. By default, it is set to 'sdate'. +#'@param window_dim A character string indicating the name of the dimension to +#' enlarge the sample data. By default NULL. #'@param ncores An integer that indicates the number of cores for parallel #' computations using multiApply function. The default value is NULL. #'@return An object of class \code{s2dv_cube} containing the bias corrected @@ -46,7 +48,7 @@ #'@export CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, memb_dim = 'member', sdate_dim = 'sdate', - ncores = NULL) { + window_dim = NULL, ncores = NULL) { if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) { stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") @@ -63,6 +65,7 @@ CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, memb_dim = memb_dim, sdate_dim = sdate_dim, + window_dim = window_dim, na.rm = na.rm, ncores = ncores) pos <- match(dimnames, names(dim(BiasCorrected))) @@ -105,6 +108,8 @@ CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, #' dimension. By default, it is set to 'member'. #'@param sdate_dim A character string indicating the name of the start date #' dimension. By default, it is set to 'sdate'. +#'@param window_dim A character string indicating the name of the dimension to +#' enlarge the sample data. By default NULL. #'@param ncores An integer that indicates the number of cores for parallel #' computations using multiApply function. The default value is NULL. #' @@ -127,7 +132,7 @@ CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, #'@export BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, memb_dim = 'member', sdate_dim = 'sdate', - ncores = NULL) { + window_dim = NULL, ncores = NULL) { # Check inputs ## exp, obs if (!is.array(exp) || !is.numeric(exp)) { @@ -196,42 +201,67 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, if (memb_dim %in% names(dim(obs))) { target_dims_obs <- c(memb_dim, target_dims_obs) } - - if (is.null(exp_cor)) { - BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), - target_dims = list(target_dims_obs, + if (is.null(window_dim)) { + if (is.null(exp_cor)) { + BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), + target_dims = list(target_dims_obs, c(memb_dim, sdate_dim)), - fun = .sbc, - na.rm = na.rm, ncores = ncores)$output1 - } else { - BiasCorrected <- Apply(data = list(var_obs = obs, - var_exp = exp, - var_cor = exp_cor), + fun = .sbc, + na.rm = na.rm, ncores = ncores)$output1 + } else { + BiasCorrected <- Apply(data = list(var_obs = obs, + var_exp = exp, + var_cor = exp_cor), target_dims = list(target_dims_obs, c(memb_dim, sdate_dim), c(memb_dim, sdate_dim)), fun = .sbc, output_dims = c(memb_dim, sdate_dim), na.rm = na.rm, ncores = ncores)$output1 - } + } + } else { + if (is.null(exp_cor)) { + BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp), + target_dims = list(c(target_dims_obs, window_dim), + c(memb_dim, sdate_dim, window_dim)), + fun = .sbc, + na.rm = na.rm, ncores = ncores)$output1 + } else { + BiasCorrected <- Apply(data = list(var_obs = obs, + var_exp = exp, + var_cor = exp_cor), + target_dims = list(c(target_dims_obs, windown_dim), + c(memb_dim, sdate_dim, window_dim), + c(memb_dim, sdate_dim, window_dim)), + fun = .sbc, + output_dims = c(memb_dim, sdate_dim), + na.rm = na.rm, ncores = ncores)$output1 + } + } return(BiasCorrected) } .sbc <- function(var_obs, var_exp , var_cor = NULL, na.rm = FALSE) { - + #var_obs: [sdate, (swindow)] ---> subseasonal: syear, sday + #var_exp: [member, sdate, (swindow)] --> subseasonal: member, syear, sday + #var_cor: [member, sdate, (swindow)] --> subsesonal: member, syear, sday ntime <- dim(var_exp)[2] corrected <- NA * var_exp if (is.null(var_cor)) { for (t in 1:ntime) { # parameters - sd_obs <- sd(var_obs[-t], na.rm = na.rm) - sd_exp <- sd(var_exp[, -t], na.rm = na.rm) - clim_exp <- mean(var_exp[, -t], na.rm = na.rm) - clim_obs <- mean(var_obs[-t], na.rm = na.rm) - - # bias corrected forecast - corrected[, t] <- ((var_exp[, t] - clim_exp) * (sd_obs / sd_exp)) + clim_obs + sd_obs <- sd(Subset(var_obs, along = 1, indices = -t), na.rm = na.rm) + sd_exp <- sd(Subset(var_exp, along = 2, indices = -t), na.rm = na.rm) + clim_exp <- mean(Subset(var_exp, along = 2, indices = -t), na.rm = na.rm) + clim_obs <- mean(Subset(var_obs, along = 1, indices = -t), na.rm = na.rm) + if (length(dim(var_exp)) == 2) { + corrected[,t] <- ((Subset(var_exp, along = 2, indices = t) - clim_exp) * + (sd_obs/sd_exp)) + clim_obs + } else { + corrected[,t,] <- ((Subset(var_exp, along = 2, indices = t) - clim_exp) * + (sd_obs/sd_exp)) + clim_obs + } } } else { # parameters @@ -239,9 +269,13 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, sd_exp <- sd(var_exp, na.rm = na.rm) clim_exp <- mean(var_exp, na.rm = na.rm) clim_obs <- mean(var_obs, na.rm = na.rm) - - # bias corrected forecast - corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs + if (length(dim(var_exp)) == 2) { + # bias corrected forecast + corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs + } else { + corrected <- ((Subset(var_cor, along = 'sday', indices = 2) -clim_exp) * + (sd_obs / sd_exp)) + clim_obs + } } return(corrected) diff --git a/man/BiasCorrection.Rd b/man/BiasCorrection.Rd index d944a32f..20ff7260 100644 --- a/man/BiasCorrection.Rd +++ b/man/BiasCorrection.Rd @@ -11,6 +11,7 @@ BiasCorrection( na.rm = FALSE, memb_dim = "member", sdate_dim = "sdate", + window_dim = NULL, ncores = NULL ) } @@ -35,6 +36,9 @@ dimension. By default, it is set to 'member'.} \item{sdate_dim}{A character string indicating the name of the start date dimension. By default, it is set to 'sdate'.} +\item{window_dim}{A character string indicating the name of the dimension to +enlarge the sample data. By default NULL.} + \item{ncores}{An integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL.} } diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd index fb96babd..4e4bcee8 100644 --- a/man/CST_BiasCorrection.Rd +++ b/man/CST_BiasCorrection.Rd @@ -11,6 +11,7 @@ CST_BiasCorrection( na.rm = FALSE, memb_dim = "member", sdate_dim = "sdate", + window_dim = NULL, ncores = NULL ) } @@ -35,6 +36,9 @@ dimension. By default, it is set to 'member'.} \item{sdate_dim}{A character string indicating the name of the start date dimension. By default, it is set to 'sdate'.} +\item{window_dim}{A character string indicating the name of the dimension to +enlarge the sample data. By default NULL.} + \item{ncores}{An integer that indicates the number of cores for parallel computations using multiApply function. The default value is NULL.} } -- GitLab From 21faf01a8d1fdaed076da95f8ecb51f27ce6c377 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 26 Oct 2022 16:33:40 +0200 Subject: [PATCH 2/4] var_obs always has member dim --- R/CST_BiasCorrection.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index f412ecde..2b8da546 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -242,7 +242,7 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, } .sbc <- function(var_obs, var_exp , var_cor = NULL, na.rm = FALSE) { - #var_obs: [sdate, (swindow)] ---> subseasonal: syear, sday + #var_obs: [member, sdate, (swindow)] ---> subseasonal: member, syear, sday #var_exp: [member, sdate, (swindow)] --> subseasonal: member, syear, sday #var_cor: [member, sdate, (swindow)] --> subsesonal: member, syear, sday ntime <- dim(var_exp)[2] @@ -251,10 +251,10 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, if (is.null(var_cor)) { for (t in 1:ntime) { # parameters - sd_obs <- sd(Subset(var_obs, along = 1, indices = -t), na.rm = na.rm) + sd_obs <- sd(Subset(var_obs, along = 2, indices = -t), na.rm = na.rm) sd_exp <- sd(Subset(var_exp, along = 2, indices = -t), na.rm = na.rm) clim_exp <- mean(Subset(var_exp, along = 2, indices = -t), na.rm = na.rm) - clim_obs <- mean(Subset(var_obs, along = 1, indices = -t), na.rm = na.rm) + clim_obs <- mean(Subset(var_obs, along = 2, indices = -t), na.rm = na.rm) if (length(dim(var_exp)) == 2) { corrected[,t] <- ((Subset(var_exp, along = 2, indices = t) - clim_exp) * (sd_obs/sd_exp)) + clim_obs -- GitLab From 69c62769221dd07db31f7b0a6cd62128f904118f Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 27 Oct 2022 14:54:35 +0200 Subject: [PATCH 3/4] drop dims for the case sday --- R/CST_BiasCorrection.R | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 2b8da546..2c60f9b2 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -230,7 +230,7 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp, var_cor = exp_cor), - target_dims = list(c(target_dims_obs, windown_dim), + target_dims = list(c(target_dims_obs, window_dim), c(memb_dim, sdate_dim, window_dim), c(memb_dim, sdate_dim, window_dim)), fun = .sbc, @@ -246,9 +246,9 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, #var_exp: [member, sdate, (swindow)] --> subseasonal: member, syear, sday #var_cor: [member, sdate, (swindow)] --> subsesonal: member, syear, sday ntime <- dim(var_exp)[2] - corrected <- NA * var_exp if (is.null(var_cor)) { + corrected <- NA * var_exp for (t in 1:ntime) { # parameters sd_obs <- sd(Subset(var_obs, along = 2, indices = -t), na.rm = na.rm) @@ -269,12 +269,15 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, sd_exp <- sd(var_exp, na.rm = na.rm) clim_exp <- mean(var_exp, na.rm = na.rm) clim_obs <- mean(var_obs, na.rm = na.rm) - if (length(dim(var_exp)) == 2) { + if (length(dim(var_cor)) == 2) { # bias corrected forecast corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs } else { - corrected <- ((Subset(var_cor, along = 'sday', indices = 2) -clim_exp) * +print("HERE") +print(dim(var_cor)) + corrected <- ((Subset(var_cor, along = 3, indices = 2, drop = 'selected') - clim_exp) * (sd_obs / sd_exp)) + clim_obs +print(dim(corrected)) } } -- GitLab From b0e130d8099b1c70cc7db0829960be025a483e07 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 27 Oct 2022 17:06:32 +0200 Subject: [PATCH 4/4] fix forecast with window --- R/CST_BiasCorrection.R | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R index 2c60f9b2..6b9f5398 100644 --- a/R/CST_BiasCorrection.R +++ b/R/CST_BiasCorrection.R @@ -273,11 +273,9 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, # bias corrected forecast corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs } else { -print("HERE") -print(dim(var_cor)) + # Todo: flexible position of the window subset? corrected <- ((Subset(var_cor, along = 3, indices = 2, drop = 'selected') - clim_exp) * (sd_obs / sd_exp)) + clim_obs -print(dim(corrected)) } } -- GitLab