#'@rdname CST_Analogs
#'@title Downscaling using Analogs based on large scale fields.
#' 
#'@author M. Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Maria M. Chaves-Montero, \email{mariadm.chaves@cmcc.it }
#'@author Nuria Perez-Zanon \email{nuria.perez@bsc.es}
#'
#'@description This function perform a downscaling using Analogs. To compute 
#'the analogs, the function search for days with similar large scale conditions
#'to downscaled fields to a local scale. The large scale and the local scale 
#'regions are defined by the user. The large scale is usually given by 
#'atmospheric circulation as sea level pressure or geopotential height 
#'(Yiou et al, 2013) but the function gives the possibility to use another 
#'field. The local scale will be usually given by precipitation or temperature 
#'fields, but might be another variable.The analogs function will find the best
#'analogs based in three criterias:
#' (1) Minimum Euclidean distance in the large scale pattern (i.e. SLP)
#' (2) Minimum Euclidean distance in the large scale pattern (i.e. SLP) and 
#' minimum Euclidean distance in the local scale pattern (i.e. SLP). 
#' (3) Minimum Euclidean distance in the large scale pattern (i.e. SLP), minimum 
#' Euclidean distance in the local scale pattern (i.e. SLP) and highest 
#' rank-based (Spearman) correlation in the local variable to downscale 
#' (i.e Precipitation).
#'The search of analogs must be done in the longest dataset posible. This is 
#'important since it is necessary to have a good representation of the 
#'possible states of the field in the past, and therefore, to get better 
#'analogs. Once the search of the analogs is complete, and in order to use 
#'the three criterias the user can select a number of analogs, using parameter
#''nAnalogs' to restrict the selection of the best analogs in a short number 
#'of posibilities, the best ones.
#'This function has not constrains of specific regions, variables to downscale,
#'or data to be used (seasonal forecast data, climate projections data, 
#'reanalyses data). The regrid into a finner scale is done interpolating with 
#'CST_Load. Then, this interpolation is corrected selecting the analogs in the 
#'large and local scale in based of the observations. The function is an 
#'adapted version of the method of Yiou et al 2013.       
#'
#'@references Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard,
#' and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column 
#' from surface pressure using analogues.  Clim. Dyn., 41, 1419-1437. 
#' \email{pascal.yiou@lsce.ipsl.fr}
#'
#'@param expL an 's2dv_cube' object containing the experimental field on the 
#'large scale for which the analog is aimed. This field is used to in all the
#'criterias. If parameter 'expVar' is not provided, the function will return 
#'the expL analog. The element 'data' in the 's2dv_cube' object must have, at
#'least, latitudinal and longitudinal dimensions. The object is expect to be 
#'already subset for the desired large scale region.
#'@param obsL an 's2dv_cube' object containing the observational field on the
#'large scale. The element 'data' in the 's2dv_cube' object must have the same 
#'latitudinal and longitudinal dimensions as parameter 'expL' and a temporal 
#'dimension with the maximum number of available observations.
#'@param time_obsL a character string indicating the date of the observations 
#'in the format "dd/mm/yyyy"
#'@param time_expL a character string indicating the date of the experiments 
#'in the format "dd/mm/yyyy"
#'@param excludeTime a character string indicating the date of the observations 
#'in the format "dd/mm/yyyy" to be excluded during the search of analogs, in a 
#'forecast might be NULL, if is not a forecast can not be NULL
#'@param expVar an 's2dv_cube' object containing the experimental field on the
#'local scale, usually a different variable to the parameter 'expL'. If it is 
#'not NULL (by default, NULL), the returned field by this function will be the 
#'analog of parameter 'expVar'.
#'@param obsVar an 's2dv_cube' containing the field of the same variable as the 
#'passed in parameter 'expVar' for the same region.
#'@param region a vector of length four indicating the minimum longitude, the 
#'maximum longitude, the minimum latitude and the maximum latitude.
#'@param criteria a character string indicating the criteria to be used for the 
#'selection of analogs:
#'\itemize{
#'\item{Large_dist} minimum distance in the large scale pattern;
#'\item{Local_dist} minimum distance in the large scale pattern and minimum 
#' distance in the local scale pattern; and
#'\item{Local_cor} minimum distance in the large scale pattern, minimum 
#' distance in the local scale pattern and highest correlation in the 
#' local variable to downscale.}
#' 
#'@import multiApply
#'@import s2dverification
#'@import abind
#'@import ClimProjDiags
#'
#'@seealso code{\link{CST_Load}}, \code{\link[s2dverification]{Load}} and 
#'\code{\link[s2dverification]{CDORemap}}
#'
#'@return An 's2dv_cube' object containing the dowscaled values of the best 
#'analogs in the criteria selected. 
#'@examples
#'res <- CST_Analogs(expL = lonlat_data$exp, obsL = lonlat_data$obs)
CST_Analogs <- function(expL, obsL, time_obsL, expVar = NULL, obsVar = NULL,
                        region = NULL, criteria = "Large_dist") {
  
  if (!inherits(expL, 's2dv_cube') || !inherits(obsL, 's2dv_cube')) {
    stop("Parameter 'expL' and 'obsL' must be of the class 's2dv_cube', ",
         "as output by CSTools::CST_Load.")
  }
  
  if (!is.null(expVar) || !is.null(obsVar)) {
    if (!inherits(expVar, 's2dv_cube') || !inherits(obsVar, 's2dv_cube')) {
      stop("Parameter 'expVar' and 'obsVar' must be of the class 's2dv_cube',
           ","as output by CSTools::CST_Load.")
    }
  }
  
  timevector <- obsL$Dates$start
  
  if (!is.null(expVar)) {
    region <- c(min(expVar$lon), max(expVar$lon), 
                min(expVar$lat), max(expVar$lon))
    lonVar <- expVar$lon
    latVar <- expVar$lat
  } else {
    region <- c(min(expL$lon), max(expL$lon), 
                min(expL$lat), max(expL$lon))
    lonVar <- expL$lon
    latVar <- expL$lat
  }
  
  result <- Analogs(expL$data, obsL$data, 
                    time_obsL = timevector,
                    expVar = expVar$data, obsVar = obsVar$data,
                    criteria = criteria, 
                    lonVar = expVar$lon, latVar = expVar$lat, 
                    region = region, nAnalogs = 1, AnalogsInfo = FALSE)
  
  if (!is.null(obsVar)) {
    obsVar$data <- result$AnalogsFields
    return(obsVar)
  } else {
    obsL$data <- result$AnalogsFields
    return(obsL)
  }
}