#'Compute the North Atlantic Oscillation (NAO) Index
#'
#'Compute the North Atlantic Oscillation (NAO) index based on the leading EOF 
#'of the sea level pressure (SLP) anomalies over the north Atlantic region 
#'(20N-80N, 80W-40E). The PCs are obtained by projecting the forecast and 
#'observed anomalies onto the observed EOF pattern (Pobs) or the forecast 
#'anomalies onto the EOF pattern of the other years of the forecast (Pmod). 
#'By default (ftime_avg = 2:4) NAO() computes the NAO index for 1-month 
#'lead seasonal forecasts that can be plotted with PlotBoxWhisker(). It returns
#'cross-validated PCs of the NAO index for forecast (exp) and observations 
#'(obs) based on the leading EOF pattern.
#'
#'@param exp A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
#'  forecast anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
#'  dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
#'  If only NAO of observational data needs to be computed, this parameter can
#'  be left to NULL. The default value is NULL.
#'@param obs A named numeric array of North Atlantic SLP (20N-80N, 80W-40E) 
#'  observed anomalies from \code{Ano()} or \code{Ano_CrossValid()} with 
#'  dimensions 'time_dim', 'memb_dim', 'ftime_dim', and 'space_dim' at least.
#'  If only NAO of experimental data needs to be computed, this parameter can 
#'  be left to NULL. The default value is NULL.
#'@param lat A vector of the latitudes of 'exp' and 'obs'.
#'@param lon A vector of the longitudes of 'exp' and 'obs'.
#'@param time_dim A character string indicating the name of the time dimension
#' of 'exp' and 'obs'. The default value is 'sdate'. 
#'@param memb_dim A character string indicating the name of the member 
#'  dimension of 'exp' and 'obs'. The default value is 'member'.
#'@param space_dim A vector of two character strings. The first is the dimension
#'  name of latitude of 'ano' and the second is the dimension name of longitude
#'  of 'ano'. The default value is c('lat', 'lon').
#'@param ftime_dim A character string indicating the name of the forecast time 
#'  dimension of 'exp' and 'obs'. The default value is 'ftime'.
#'@param ftime_avg A numeric vector of the forecast time steps to average
#'  across the target period. The default value is 2:4, i.e., from 2nd to 4th 
#'  forecast time steps.
#'@param obsproj A logical value indicating whether to compute the NAO index by
#'  projecting the forecast anomalies onto the leading EOF of observational 
#'  reference (TRUE) or compute the NAO by first computing the leading 
#'  EOF of the forecast anomalies (in cross-validation mode, i.e. leaving the 
#'  year you are evaluating out), and then projecting forecast anomalies onto 
#'  this EOF (FALSE). 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 
#'A list which contains:
#'\item{exp}{
#'  A numeric array of forecast NAO index in verification format with the same 
#'  dimensions as 'exp' except space_dim and ftime_dim.
#'  }
#'\item{obs}{
#'  A numeric array of observed NAO index in verification format with the same
#'  dimensions as 'obs' except space_dim and ftime_dim.
#'}
#'
#'@references
#'Doblas-Reyes, F.J., Pavan, V. and Stephenson, D. (2003). The skill of 
#'  multi-model seasonal forecasts of the wintertime North Atlantic 
#'  Oscillation. Climate Dynamics, 21, 501-514. 
#'  DOI: 10.1007/s00382-003-0350-4
#'
#'@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)
#'# No example data is available over NAO region, so in this example we will 
#'# tweak the available data. In a real use case, one can Load() the data over 
#'# the NAO region directly.
#'sampleData$lon[] <- c(40, 280, 340)
#'sampleData$lat[] <- c(20, 80)
#'  }
#'
#'# Now ready to compute the EOFs and project on, for example, the first 
#'# variability mode.
#'ano <- Ano_CrossValid(sampleData$mod, sampleData$obs)
#'# Note that computing the NAO over the region for which there is available 
#'# example data is not the full NAO area: NAO() will raise a warning.
#'nao <- NAO(ano$exp, ano$obs, sampleData$lat, sampleData$lon)
#'# Finally plot the NAO index
#'  \dontrun{
#'nao$exp <- Reorder(nao$exp, c(2, 1))
#'nao$obs <- Reorder(nao$obs, c(2, 1))
#'PlotBoxWhisker(nao$exp, nao$obs, "NAO index, DJF", "NAO index (PC1) TOS",
#'        monini = 12, yearini = 1985, freq = 1, "Exp. A", "Obs. X")
#'  }
#'
#'@import multiApply
#'@importFrom ClimProjDiags Subset
#'@export
NAO <- function(exp = NULL, obs = NULL, lat, lon, time_dim = 'sdate',
                memb_dim = 'member', space_dim = c('lat', 'lon'),
                ftime_dim = 'ftime', ftime_avg = 2:4, 
                obsproj = TRUE, ncores = NULL) {

  # Check inputs 
  ## exp and obs (1)
  if (is.null(obs) & is.null(exp)) {
    stop("Parameter 'exp' and 'obs' cannot both be NULL.")
  }
  if (!is.null(exp)) {
    if (!is.numeric(exp)) {
      stop("Parameter 'exp' must be a numeric array.")
    }
    if (is.null(dim(exp))) {
      stop(paste0("Parameter 'exp' and must have at least dimensions ",
                  "time_dim, memb_dim, space_dim, and ftime_dim."))
    }
    if(any(is.null(names(dim(exp)))) | any(nchar(names(dim(exp))) == 0)) {
      stop("Parameter 'exp' must have dimension names.")
    }
  }
  if (!is.null(obs)) {
    if (!is.numeric(obs)) {
      stop("Parameter 'obs' must be a numeric array.")
    }
    if (is.null(dim(obs))) {
      stop(paste0("Parameter 'obs' and must have at least dimensions ",
                  "time_dim, memb_dim, space_dim, and ftime_dim."))
    }
    if(any(is.null(names(dim(obs)))) | any(nchar(names(dim(obs))) == 0)) {
      stop("Parameter 'obs' must have dimension names.")
    }
  }
  ## time_dim
  if (!is.character(time_dim) | length(time_dim) > 1) {
    stop("Parameter 'time_dim' must be a character string.")
  }
  if (!is.null(exp)) {
    if (!time_dim %in% names(dim(exp))) {
      stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  if (!is.null(obs)) {
    if (!time_dim %in% names(dim(obs))) {
      stop("Parameter 'time_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 (!is.null(exp)) {
    if (!memb_dim %in% names(dim(exp))) {
      stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  if (!is.null(obs)) {
    if (!memb_dim %in% names(dim(obs))) {
      stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  ## space_dim
  if (!is.character(space_dim) | length(space_dim) != 2) {
    stop("Parameter 'space_dim' must be a character vector of 2.")
  }
  if (!is.null(exp)) {
    if (any(!space_dim %in% names(dim(exp)))) {
      stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  if (!is.null(obs)) {
    if (any(!space_dim %in% names(dim(obs)))) {
      stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  ## ftime_dim
  if (!is.character(ftime_dim) | length(ftime_dim) > 1) {
    stop("Parameter 'ftime_dim' must be a character string.")
  }
  if (!is.null(exp)) {
    if (!ftime_dim %in% names(dim(exp))) {
      stop("Parameter 'ftime_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  if (!is.null(obs)) {
    if (!ftime_dim %in% names(dim(obs))) {
      stop("Parameter 'ftime_dim' is not found in 'exp' or 'obs' dimension.")
    }
  }
  ## exp and obs (2)
  if (!is.null(exp) & !is.null(obs)) {
    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)]
    if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
      stop(paste0("Parameter 'exp' and 'obs' must have the same length of ",
                  "all dimensions except 'memb_dim'."))
    }
  }
  ## ftime_avg
  if (!is.vector(ftime_avg) | !is.integer(ftime_avg)) {
    stop("Parameter 'ftime_avg' must be an integer vector.")
  }
  if (!is.null(exp)) {
    if (max(ftime_avg) > dim(exp)[ftime_dim] | min(ftime_avg) < 1) {
      stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
    }
  } else {
    if (max(ftime_avg) > dim(obs)[ftime_dim] | min(ftime_avg) < 1) {
      stop("Parameter 'ftime_avg' must be within the range of ftime_dim length.")
    }
  }
  ## sdate >= 2
  if (!is.null(exp)) {
    if (dim(exp)[time_dim] < 2) {
      stop("The length of time_dim must be at least 2.")
    }
  } else {
    if (dim(obs)[time_dim] < 2) {
      stop("The length of time_dim must be at least 2.")
    }
  }
  ## lat and lon
  if (!is.null(exp)) {
    if (!is.numeric(lat) | length(lat) != dim(exp)[space_dim[1]]) {
      stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
                  "length as the latitude dimension of 'exp' and 'obs'."))
    }
    if (!is.numeric(lon) | length(lon) != dim(exp)[space_dim[2]]) {
      stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
                  "length as the longitude dimension of 'exp' and 'obs'."))
    }
  } else {
    if (!is.numeric(lat) | length(lat) != dim(obs)[space_dim[1]]) {
      stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
                  "length as the latitude dimension of 'exp' and 'obs'."))
    }
    if (!is.numeric(lon) | length(lon) != dim(obs)[space_dim[2]]) {
      stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
                  "length as the longitude dimension of 'exp' and 'obs'."))
    }
  } 
  stop_needed <- FALSE
  if (tail(lat, 1) < 70 | tail(lat, 1) > 90 | 
      head(lat, 1) > 30 | head(lat, 1) < 10) {
    stop_needed <- TRUE
  }
  #NOTE: different from s2dverification
  # lon is not used in the calculation actually. EOF only uses lat to do the
  # weight. So we just need to ensure the data is in this region, regardless
  # the order. 
  if (any(lon < 0)) {  #[-180, 180] 
    if (!(min(lon) > -90 & min(lon) < -70 & max(lon) < 50 & max(lon) > 30)) {
      stop_needed <- TRUE
    }
  } else {  #[0, 360]
    if (any(lon >= 50 & lon <= 270)) {
      stop_needed <- TRUE
    } else {
      lon_E <- lon[which(lon < 50)]
      lon_W <- lon[-which(lon < 50)]
      if (max(lon_E) < 30 | min(lon_W) > 290) {
        stop_needed <- TRUE
      }
    }
  }
  if (stop_needed) {
    stop(paste0("The typical domain used to compute the NAO is 20N-80N, ",
                "80W-40E. 'lat' or 'lon' is out of range."))
  }
  ## obsproj
  if (!is.logical(obsproj)  | length(obsproj) > 1) {
    stop("Parameter 'obsproj' must be either TRUE or FALSE.")
  }
  if (obsproj) {
    if (is.null(obs)) {
      stop("Parameter 'obsproj' set to TRUE but no 'obs' provided.")
    }
    if (is.null(exp)) {
      .warning("parameter 'obsproj' set to TRUE but no 'exp' provided.")
    }
  }
  ## 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.")
    }
  }

  #-------- Average ftime -----------
  if (!is.null(exp)) {
    exp_sub <- ClimProjDiags::Subset(exp, ftime_dim, ftime_avg, drop = FALSE)
    exp <- MeanDims(exp_sub, ftime_dim, na.rm = TRUE)
    ## Cross-validated PCs. Fabian. This should be extended to
    ## nmod and nlt by simple loops. Virginie
  }

  if (!is.null(obs)) {
    obs_sub <- ClimProjDiags::Subset(obs, ftime_dim, ftime_avg, drop = FALSE)
    obs <- MeanDims(obs_sub, ftime_dim, na.rm = TRUE)
  }

  if (!is.null(exp) & !is.null(obs)) {
    res <- Apply(list(exp, obs),
                 target_dims = list(exp = c(memb_dim, time_dim, space_dim),
                                    obs = c(memb_dim, time_dim, space_dim)),
                 fun = .NAO,
                 obsproj = obsproj, lat = lat, lon = lon,
                 ncores = ncores)
  } else if (!is.null(exp)) {
    res <- Apply(list(exp = exp),
                 target_dims = list(exp = c(memb_dim, time_dim, space_dim)),
                 fun = .NAO,
                 obsproj = obsproj, lat = lat, lon = lon, obs = NULL,
                 ncores = ncores)
  } else if (!is.null(obs)) {
    res <- Apply(list(obs = obs),
                 target_dims = list(obs = c(memb_dim, time_dim, space_dim)),
                 fun = .NAO,
                 obsproj = obsproj, lat = lat, lon = lon, exp = NULL,
                 ncores = ncores)
  }
  return(res)
}

.NAO <- function(exp = NULL, obs = NULL, lat, lon,
                 obsproj = TRUE, ncores = NULL) {
  # exp: [memb_exp, sdate, lat, lon]
  # obs: [memb_obs, sdate, lat, lon]
  if (!is.null(exp)) {
    ntime <- dim(exp)[2]
    nlat <- dim(exp)[3]
    nlon <- dim(exp)[4]
    nmemb_exp <- dim(exp)[1]
    nmemb_obs <- dim(obs)[1]
  } else {
    ntime <- dim(obs)[2]
    nlat <- dim(obs)[3]
    nlon <- dim(obs)[4]
    nmemb_obs <- dim(obs)[1]
  }

  if (!is.null(obs)) NAOO.ver <- array(NA, dim = c(ntime, nmemb_obs))
  if (!is.null(exp)) NAOF.ver <- array(NA, dim = c(ntime, nmemb_exp))

  for (tt in 1:ntime) {  #sdate

    if (!is.null(obs)) {
      ## Observed EOF excluding one forecast start year.
      obs_sub <- ClimProjDiags::Subset(obs, 2, c(1:ntime)[-tt], drop = FALSE)
      obs_EOF <- EOF(obs_sub, lat = lat, lon = lon, time_dim = names(ntime),
                     space_dim = c(names(nlat), names(nlon)), neofs = 1)

      ## Correct polarity of pattern.
      #NOTE: different from s2dverification
      # dim(obs_EOF$EOFs): [mode, lat, lon, member]
      for (imemb in 1:nmemb_obs) {
        if (0 < mean(obs_EOF$EOFs[1, which.min(abs(lat - 65)), , ], na.rm = T)) {
          obs_EOF$EOFs[1, , , imemb] <- obs_EOF$EOFs[1, , , imemb] * (-1)
        }
      }
#      obs_EOF$PCs <- obs_EOF$PCs * sign  # not used

      ## Project observed anomalies.
      PF <- ProjectField(obs, eof = obs_EOF, time_dim = names(ntime),
                         space_dim = c(names(nlat), names(nlon)), mode = 1)
      NAOO.ver[tt, ] <- PF[tt, ]
      ## Keep PCs of excluded forecast start year. Fabian.
    }

    if (!is.null(exp)) {
      if (!obsproj) {
        exp_sub <- ClimProjDiags::Subset(exp, 2, c(1:ntime)[-tt], drop = FALSE)
        #NOTE: different from s2dverification. Here, 'member' is considered.
        exp_EOF <- EOF(exp_sub, lat = lat, lon = lon, time_dim = names(ntime),
                       space_dim = c(names(nlat), names(nlon)), neofs = 1)

        ## Correct polarity of pattern.
        #NOTE: different from s2dverification
        for (imemb in 1:nmemb_exp) {
          if (0 < mean(exp_EOF$EOFs[1, which.min(abs(lat - 65)), , imemb], na.rm = T)) {
          exp_EOF$EOFs[1, , , imemb] <- exp_EOF$EOFs[1, , , imemb] * (-1)
          }
        }
#        exp_EOF$PCs <- exp_EOF$PCs * sign  # not used

        ### Lines below could be simplified further by computing
        ### ProjectField() only on the year of interest... (though this is
        ### not vital). Lauriane
        PF <- ProjectField(exp, eof = exp_EOF, time_dim = names(ntime),
                           space_dim = c(names(nlat), names(nlon)), mode = 1)
        NAOF.ver[tt, ] <- PF[tt, ]
        
      } else {
      ## Project forecast anomalies on obs EOF
      #NOTE: Because obs and exp have different nmemb, do ensemble mean to 
      #      obs_EOF$EOFs first, then expand the memb dim to be the same as exp.
      obs_EOF$EOFs <- apply(obs_EOF$EOFs, c(1, 2, 3), mean, na.rm = T)
      obs_EOF$EOFs <- array(obs_EOF$EOFs, dim = c(dim(obs_EOF$EOFs), as.numeric(nmemb_exp)))
      names(dim(obs_EOF$EOFs))[4] <- names(nmemb_obs)
      PF <- ProjectField(exp, obs_EOF, mode = 1)
      NAOF.ver[tt, ] <- PF[tt, ]
      }
    }
  }
  #NOTE: EOFs_obs is not returned because it's only the result of the last sdate
  #     (It is returned in s2dverification.)
  if (!is.null(exp) & !is.null(obs)) {
    return(list(exp = NAOF.ver, obs = NAOO.ver)) #, EOFs_obs = obs_EOF))
  } else if (!is.null(exp)) {
    return(list(exp = NAOF.ver))
  } else if (!is.null(obs)) {
    return(list(obs = NAOO.ver))
  }
}