CST_DynBiasCorrection.R

#'@rdname CST_DynBiasCorrection
#'@title Performing a Bias Correction conditioned by the dynamical
#'properties of the data.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function perform a bias correction conditioned by the 
#'dynamical properties of the dataset. This function internally uses the functions 
#''Predictability' to divide in terciles the two dynamical proxies 
#'computed with 'CST_ProxiesAttractor'. A bias correction
#'between the model and the observations is performed using the division into
#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'.
#'For instance, model values with lower 'dim' will be corrected with observed 
#'values with lower 'dim', and the same for theta. The function gives two options
#'of bias correction: one for 'dim' and/or one for 'theta'
#'
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., 
#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large 
#'scale atmospheric predictability.Nature Communications, 10(1), 1316. 
#'DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes. 
#' Scientific Reports, 7-41278, 2017.
#'
#'@param exp an s2v_cube object with the experiment data
#'@param obs an s2dv_cube object with the reference data 
#'@param method a character string indicating the method to apply bias 
#'correction among these ones: "PTF","RQUANT","QUANT","SSPLIN"
#'@param wetday logical indicating whether to perform wet day correction 
#'or not OR a numeric threshold below which all values are set to zero (by 
#'default is set to 'FALSE').
#'@param proxy a character string indicating the proxy for local dimension
#' 'dim' or inverse of persistence 'theta' to apply the dynamical 
#' conditioned bias correction method. 
#'@param quanti a number lower than 1 indicating the quantile to perform 
#'the computation of local dimension and theta
#'@param ncores The number of cores to use in parallel computation
#'
#'@return dynbias an s2dvcube object with a bias correction performed 
#'conditioned by local dimension 'dim' or inverse of persistence 'theta'
#'
#'@examples
#'# example 1: simple data s2dvcube style
#' set.seed(1)
#' expL <- rnorm(1:2000)
#' dim (expL) <- c(time =100,lat = 4, lon = 5)
#' obsL <- c(rnorm(1:1980),expL[1,,]*1.2)
#' dim (obsL) <- c(time = 100,lat = 4, lon = 5)
#' time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-")
#' time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-")
#' lon <- seq(-1,5,1.5)
#' lat <- seq(30,35,1.5)
#' # qm=0.98 # too high for this short dataset, it is possible that doesn't
#' # get the requirement, in that case it would be necessary select a lower qm
#' # for instance qm=0.60
#' expL <- s2dv_cube(data = expL, lat = lat, lon = lon,
#'                 Dates = list(start = time_expL, end = time_expL))
#' obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon,
#'                 Dates = list(start = time_obsL, end = time_obsL))
#' # to use DynBiasCorrection
#' dynbias1 <- DynBiasCorrection(exp = expL$data, obs = obsL$data, proxy= "dim",
#'                                  quanti = 0.6)
#' # to use CST_DynBiasCorrection
#' dynbias2 <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim",
#'                                 quanti = 0.6)
#'
#'@export
CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=FALSE,
                                 proxy = "dim", quanti,
                                 ncores = NULL) {
  if (!inherits(obs, 's2dv_cube')) {
    stop("Parameter 'obs' must be of the class 's2dv_cube', ",
         "as output by CSTools::CST_Load.")
  }
  if (!inherits(exp, 's2dv_cube')) {
    stop("Parameter 'exp' must be of the class 's2dv_cube', ",
         "as output by CSTools::CST_Load.")
  }
  exp$data <- DynBiasCorrection(exp = exp$data, obs = obs$data, method = method,
                                wetday=wetday,
                                proxy = proxy, quanti = quanti, ncores = ncores)
  return(exp)
}
#'@rdname DynBiasCorrection
#'@title Performing a Bias Correction conditioned by the dynamical
#'properties of the data.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it}
#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it}
#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}
#'
#'@description This function perform a bias correction conditioned by the 
#'dynamical properties of the dataset. This function used the functions 
#''CST_Predictability' to divide in terciles the two dynamical proxies 
#'computed with 'CST_ProxiesAttractor'. A bias correction
#'between the model and the observations is performed using the division into
#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'.
#'For instance, model values with lower 'dim' will be corrected with observed 
#'values with lower 'dim', and the same for theta. The function gives two options
#'of bias correction: one for 'dim' and/or one for 'theta'
#'
#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., 
#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large 
#'scale atmospheric predictability.Nature Communications, 10(1), 1316. 
#'DOI = https://doi.org/10.1038/s41467-019-09305-8 "
#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017).
#' Dynamical proxies of North Atlantic predictability and extremes. 
#' Scientific Reports, 7-41278, 2017.
#'
#'@param exp a multidimensional array with named dimensions with the 
#'experiment data 
#'@param obs a multidimensional array with named dimensions with the 
#'observation data
#'@param method a character string indicating the method to apply bias 
#'correction among these ones:
#'"PTF","RQUANT","QUANT","SSPLIN"
#'@param wetday logical indicating whether to perform wet day correction 
#'or not OR a numeric threshold below which all values are set to zero (by 
#'default is set to 'FALSE').
#'@param proxy a character string indicating the proxy for local dimension 
#''dim' or inverse of persistence 'theta' to apply the dynamical conditioned 
#'bias correction method. 
#'@param quanti a number lower than 1 indicating the quantile to perform the 
#'computation of local dimension and theta
#'@param ncores The number of cores to use in parallel computation
#'
#'@return a multidimensional array with named dimensions with a bias correction 
#'performed conditioned by local dimension 'dim' or inverse of persistence 'theta'
#'
#'@import multiApply
#'@importFrom s2dverification Subset
#'@import qmap
#'@examples
#'expL <- rnorm(1:2000)
#'dim (expL) <- c(time =100,lat = 4, lon = 5)
#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2)
#'dim (obsL) <- c(time = 100,lat = 4, lon = 5)
#'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT',
#'                            proxy= "dim", quanti = 0.6)
#'@export
DynBiasCorrection<- function(exp, obs, method = 'QUANT',wetday=FALSE, 
                             proxy = "dim", quanti, ncores = NULL){
  if (is.null(obs)) {
    stop("Parameter 'obs' cannot be NULL.")
  }  
  if (is.null(exp)) {
    stop("Parameter 'exp' cannot be NULL.")
  }    
  if (is.null(method)) {
    stop("Parameter 'method' cannot be NULL.")
  }  
  if (is.null(quanti)) {
    stop("Parameter 'quanti' cannot be NULL.")
  }   
  if (is.null(proxy)) {
    stop("Parameter 'proxy' cannot be NULL.")
  } 
  dims <- dim(exp)

  attractor.obs <- ProxiesAttractor(data = obs, quanti = quanti)
  predyn.obs <- Predictability(dim = attractor.obs$dim,
                               theta = attractor.obs$theta)
  attractor.exp <- ProxiesAttractor(data = exp, quanti = quanti)
  predyn.exp <- Predictability(dim = attractor.exp$dim,
                               theta = attractor.exp$theta)
 
  if (!(any(names(dim(exp)) %in% 'time'))){ 
   if (any(names(dim(exp)) %in% 'sdate')) {
    if (any(names(dim(exp)) %in% 'ftime')) {
      exp <- MergeDims(exp, merge_dims = c('ftime', 'sdate'),
                        rename_dim = 'time')
    }
   }
  }
  if (!(any(names(dim(obs)) %in% 'time'))){ 
    if (any(names(dim(obs)) %in% 'sdate')) {
      if (any(names(dim(obs)) %in% 'ftime')) {
        obs <- MergeDims(obs, merge_dims = c('ftime', 'sdate'),
                         rename_dim = 'time')
      }
    }
  }
 
  dim_exp <- dim(exp) 
  names_to_check <- names(dim_exp)[which(names(dim_exp) %in% 
                          c('time', 'lat', 'lon', 'sdate') == FALSE)]
  if (length(names_to_check) > 0) {
    dim_obs <- dim(obs)
    if (any(names(dim_obs) %in% names_to_check)) {
      if (any(dim_obs[which(names(dim_obs) %in% names_to_check)] != 
          dim_exp[which(names(dim_exp) %in% names_to_check)])) {
        for (i in names_to_check) {
          pos <- which(names(dim_obs) == i)
          names(dim(obs))[pos] <- ifelse(dim_obs[pos] != 
                                  dim_exp[which(names(dim_exp) == i)], 
                                  paste0('obs_', names(dim_obs[pos])), 
                                  names(dim(obs)[pos])) 
        }
        warning("Common dimension names with different length are renamed.")
      }
    }
  } 
  
  if (proxy == "dim") {
    adjusted <- Apply(list(exp, obs), target_dims = 'time',
                      fun = .dynbias, method, wetday,
                      predyn.exp = predyn.exp$pred.dim$pos.d,
                      predyn.obs = predyn.obs$pred.dim$pos.d,
                      ncores = ncores, output_dims = 'time')$output1
  } else if (proxy == "theta") {
    adjusted <- Apply(list(exp, obs), target_dims = 'time',
                      fun = .dynbias, method, wetday,
                      predyn.exp = predyn.exp$pred.theta$pos.t,
                      predyn.obs = predyn.obs$pred.theta$pos.t,
                      ncores = ncores, output_dims = 'time')$output1
  } else {
    stop ("Parameter 'proxy' must be set as 'dim' or 'theta'.")
  }
  
  if(any(names(dim(adjusted)) %in% 'memberObs')){
    if(dim(adjusted)['memberObs'] == 1){
    adjusted <- Subset(adjusted,along='memberObs',indices=1,drop = 'selected')
     }else{
    print('Dimension member in obs changed to memberObs')
    }
  }
  
  if(any(names(dim(adjusted)) %in% 'datasetObs')){
    if(dim(adjusted)['datasetObs'] == 1){
      adjusted <- Subset(adjusted,along='datasetObs',indices=1,drop = 'selected')
    }else{
      print('Dimension dataset in obs changed to datasetObs')
    }
  }
  return(adjusted)
}

.dynbias <- function(exp, obs, method, wetday,predyn.exp, predyn.obs) {
   result <- array(rep(NA, length(exp)))
   res <- lapply(1:3, function(x) {
     exp_sub <- exp[predyn.exp[[x]]]
     obs_sub <- obs[predyn.obs[[x]]]
     adjust <- .qbiascorrection(exp_sub, obs_sub, method,wetday)
     result[predyn.exp[[x]]] <<- adjust
     return(NULL)
   })
   return(result)
}   
.qbiascorrection <- function(expX, obsX, method,wetday) {
  ## functions fitQmap and doQmap
  if (method == "PTF") {
    qm.fit <- fitQmap(obsX, expX, method = "PTF", transfun = "expasympt",
                      cost = "RSS", wet.day = wetday)
    qmap <- doQmap(expX, qm.fit)
  } else if (method == "QUANT") {
    qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01, wet.day = wetday)
    qmap <- doQmap(expX, qm.fit, type = "tricub")
  } else if (method == "RQUANT") {
    qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01,wet.day = wetday)
    qmap <- doQmap(expX, qm.fit, type = "linear")
  } else if (method == "SSPLIN") {
    qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN",wet.day = wetday)
    qmap <- doQmap(expX, qm.fit)
  } else {
    stop ("Parameter 'method' doesn't match any of the available methods.")
  }
  return(qmap)
}