#'Compute forecast or observation anomalies
#'
#'This function computes anomalies from a multidimensional data array and a 
#'climatology array. 
#'
#'@param data  A numeric array with named dimensions, representing the model or
#'  observational data to be calculated the anomalies. It should involve all
#'  the dimensions in parameter 'clim', and it can have more other dimensions.
#'@param clim A numeric array with named dimensions, representing the 
#'  climatologies to be deducted from parameter 'data'. It can be generated by 
#'  Clim(). The dimensions should all be involved in parameter 'data' with the
#'  same length.
#'@param ncores An integer indicating the number of cores to use for parallel 
#'  computation. The default value is NULL.
#'
#'@return An array with same dimensions as parameter 'data'.
#'
#'@examples 
#'# Load sample data as in Load() example:
#'example(Load)
#'clim <- Clim(sampleData$mod, sampleData$obs)
#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
#'\donttest{
#'PlotAno(ano_exp, ano_obs, startDates, 
#'        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
#'}
#'@import multiApply
#'@export
  Ano <- function(data, clim, ncores = NULL) {

  # Check inputs 
  ## data
  if (is.null(data)) {
    stop("Parameter 'data' cannot be NULL.")
  }
  if (!is.numeric(data)) {
    stop("Parameter 'data' must be a numeric array.")
  }
  if (is.null(dim(data))) {  #is vector
    dim(data) <- c(length(data))
    names(dim(data)) <- 'tmp_name'
  }
  if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) {
    stop("Parameter 'data' must have dimension names.")
  }
  ## clim
  if (is.null(clim)) {
    stop("Parameter 'clim' cannot be NULL.")
  }
  if (!is.numeric(clim)) {
    stop("Parameter 'clim' must be a numeric array.")
  }
  if (is.null(dim(clim))) {  #is vector
    dim(clim) <- c(length(clim))
    names(dim(clim)) <- 'tmp_name'
  }
  if (any(is.null(names(dim(clim))))| any(nchar(names(dim(clim))) == 0)) {
    stop("Parameter 'clim' must have dimension names.")
  }
  ## data and clim
  if (!all(names(dim(clim)) %in% names(dim(data)))) {
    stop("Parameter 'data' must have all the dimensions of parameter 'clim'.")
  } else {
    pos <- names(dim(data))[match(names(dim(clim)), names(dim(data)))]
    if (any(dim(clim) != dim(data)[pos])) {
      stop("Some dimensions of parameter 'clim' have different length from parameter 'data'.")
    }
  }
  ## ncores
  if (!is.null(ncores)) {
    if (!is.numeric(ncores)) {
      stop("Parameter 'ncores' must be a positive integer.")
    } else if (ncores %% 1 != 0 | ncores <= 0 | length(ncores) > 1) {
      stop("Parameter 'ncores' must be a positive integer.")
    }
  }

  ###############################
  # Calculate Ano
  parallel_compute <- TRUE
  if (is.null(ncores)) {
    parallel_compute <- FALSE
  } else if (ncores == 1) {
    parallel_compute <- FALSE
  } 
  if (!parallel_compute) {
    target_dims_ind <- match(names(dim(clim)), names(dim(data)))
    if (any(target_dims_ind != sort(target_dims_ind))) {
      clim <- Reorder(clim, match(sort(target_dims_ind), target_dims_ind))
    }
    if (length(dim(data)) == length(dim(clim))) {
      res <- data - clim
    } else {
      target_dims_ind <- match(names(dim(clim)), names(dim(data)))
      margin_dims_ind <- c(1:length(dim(data)))[-target_dims_ind]
      res <- apply(data, margin_dims_ind, .Ano, clim)
      res <- array(res, dim = dim(data)[c(target_dims_ind, margin_dims_ind)])
    }
  } else {
    res <- Apply(list(data),
                 target_dims = names(dim(clim)),
                 output_dims = names(dim(clim)),
                 fun = .Ano,
                 clim = clim,
                 ncores = ncores)$output1
  }

  # Reorder dim back to data
  if (any(dim(res) != dim(data))) {
    res <- Reorder(res, names(dim(data)))
  }

  return(invisible(res))
  }

  .Ano <- function(data, clim) {
    data - clim
  }