#'Compute the Mean Bias
#'
#'The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference 
#'between the ensemble mean forecast and the observations. It is a deterministic
#'metric. Positive values indicate that the forecasts are on average too high
#'and negative values indicate that the forecasts are on average too low; however,
#'it gives no information about the typical magnitude of individual forecast errors.
#'
#'@param exp A named numerical array of the forecast with at least time 
#'  dimension.  
#'@param obs A named numerical array of the observation with at least time 
#'  dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 
#'  'dat_dim'.
#'@param time_dim A character string indicating the name of the time dimension.
#'  The default value is 'sdate'.
#'@param memb_dim A character string indicating the name of the member dimension
#'  to compute the ensemble mean; it should be set to NULL if the parameter 'exp'
#'  is already the ensemble mean. 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 numerical array of Bias with dimensions the dimensions of
#''exp' except 'time_dim' and 'memb_dim' dimensions.
#'
#'@references 
#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
#'
#'@examples
#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50))
#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50))
#'bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1)
#'
#'@import multiApply
#'@export

Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) {
  
  ## Checks
  
  
  ## Ensemble mean
  if (!is.null(memb_dim)) {
    exp <- MeanDims(exp, memb_dim, na.rm = na.rm)
  }
  
  ## Mean bias
  bias <- multiApply::Apply(data = list(exp, obs),
                            target_dims = time_dim,
                            fun = .Bias,
                            ncores = ncores)$output1
  ## Return the mean bias
  bias <- MeanDims(bias, time_dim, na.rm = na.rm)
  
  return(bias)
}

.Bias <- function(exp, obs) {
  
  bias <- exp - obs
  
  return(bias)
}