#'@rdname CST_ProxiesAttractor
#'@title Computing two dinamical proxies of the attractor in s2dvcube.
#'
#'@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 computes two dinamical proxies of the attractor: 
#'The local dimension (d) and the inverse of the persistence (theta) for an
#''s2dvcube' object.
#'These two parameters will be used as a condition for the computation of dynamical 
#'scores to measure predictability and to compute bias correction conditioned  by
#'the dynamics with the function DynBiasCorrection  
#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in 
#'@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 data data to create the attractor. Must be a matrix with the timesteps in nrow 
#'and the grids in ncol(dat(time,grids)
# 
#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to
#' be provided as a list of length 2 with longitudes and latitudes 
#' 
#'@param ncores The number of cores to use in parallel computation
#' 
#'@return dim and theta
#'@import s2dverification
#'@import fields
#'@import multiApply
#'
#'@examples
#'# Example 1: Computing the attractor using simple s2dv data
#'expL <- rnorm(1:200)
#'dim(expL) <- c(member=10,lat = 4, lon = 5)
#'obsL <- c(rnorm(1:180),expL[1,,]*1.2)
#'dim(obsL) <- c(time = 10,lat = 4, lon = 5)
#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-")
#'time_expL <- time_obsL[1]
#'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
#'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))
#'
#'attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm)
#'
CST_ProxiesAttractor <- function(data,quanti,ncores=NULL){
  if (is.null(data)) {
    stop("Parameter 'data' is mandatory")
  }  
  if (!inherits(data, 's2dv_cube')) {
    stop("Parameter 'data' must be of the class 's2dv_cube', ",
         "as output by CSTools::CST_Load.")
  }
  if (is.null(quanti)) {
    stop("Parameter 'quanti' is mandatory")
  }        
  if (any(names(dim(data$data)) %in% 'sdate')) {
    if (any(names(dim(data$data)) %in% 'ftime')) {
      data$data <- MergeDims(data$data, merge_dims = c('ftime', 'sdate'),
                             rename_dim = 'time')
    }
  }
  if (!(any(names(dim(data$data)) %in% 'time'))){
    stop('dimension time is mandatory')
  }
  if (any(names(dim(data$data)) %in% 'lat')) {
    if (any(names(dim(data$data)) %in% 'lon')) {
      data$data <- MergeDims(data$data, merge_dims = c('lon', 'lat'),
                             rename_dim = 'grid')
    }
  }
  if(!(any(names(dim(data$data)) %in% 'grid'))){
    stop('dimension grid is mandatory')
  }  
  
  attractor <-Apply(data$data,target_dims=c('time','grid'), fun=ProxiesAttractor,
                    quanti=quanti,iplot=FALSE, ncores=ncores)
  # rename dimensions
  attractor <- lapply(attractor, FUN=function(x,dimname){names(dim(x))[dimname] <- 'time'
  return(x)}, dimname = which(names(dim(attractor[[1]]))=='dim2'))
  
  return(list(dim=attractor$dim,theta=attractor$theta))
}