#'@rdname CST_ProxiesAttractor
#'@title Computing two dinamical proxies of the attractor in s2dv_cube.
#'
#'@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
#''s2dv_cube' 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 a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow 
#'and the grids in ncol(dat(time,grids)
# 
#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta
#' 
#'@param ncores The number of cores to use in parallel computation
#' 
#'@return dim and theta
#'
#'@examples
#'# Example 1: Computing the attractor using simple s2dv data
#'attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6)
#'
#'@export
CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){
  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' cannot be NULL.")
  }        
  
  data$data <- ProxiesAttractor(data = data$data, quanti = quanti, ncores = ncores)

  return(data)
}
#'@rdname ProxiesAttractor
#'
#'@title Computing two dinamical proxies of the attractor.
#'@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). 
#'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 a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'. 
#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta
#'@param ncores The number of cores to use in parallel computation.
#'
#'@return dim and theta
#'
#'@import multiApply
#'
#'@examples
#'# Example 1: Computing the attractor using simple data
#'# Creating an example of matrix data(time,grids):
#'mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) 
#'qm <- 0.90 # imposing a threshold
#'Attractor <- ProxiesAttractor(data = mat, quanti = qm)
#'# to plot the result
#'time = c(1:length(Attractor$theta))
#'layout(matrix(c(1, 3, 2, 3), 2, 2))
#'plot(time, Attractor$dim, xlab = 'time', ylab = 'd',
#'     main = 'local dimension', type = 'l')
#'plot(time, Attractor$theta, xlab = 'time', ylab = 'theta', main = 'theta') 
#'plot(Attractor$dim, Attractor$theta, col = 'blue', 
#'     main = "Proxies of the Attractor",
#'     xlab = "local dimension", ylab = "theta", lwd = 8, 'p')
#'
#'@export

ProxiesAttractor <- function(data, quanti, ncores = NULL){
  if (is.null(data)) {
    stop("Parameter 'data' cannot be NULL.")
  }  
  if (is.null(quanti)) {
    stop("Parameter 'quanti' is mandatory")
  }    
   if (any(names(dim(data)) %in% 'sdate')) {
    if (any(names(dim(data)) %in% 'ftime')) {
      data <- MergeDims(data, merge_dims = c('ftime', 'sdate'),
                             rename_dim = 'time')
    }
  }
  if (!(any(names(dim(data)) %in% 'time'))){
    stop("Parameter 'data' must have a temporal dimension named 'time'.")
  }
  if (any(names(dim(data)) %in% 'lat')) {
    if (any(names(dim(data)) %in% 'lon')) {
      data <- MergeDims(data, merge_dims = c('lon', 'lat'),
                             rename_dim = 'grid')
    }
  }
  if (any(names(dim(data)) %in% 'latitude')) {
    if (any(names(dim(data)) %in% 'longitude')) {
      data <- MergeDims(data, merge_dims = c('longitude', 'latitude'),
                             rename_dim = 'grid')
    }
  }
  if(!(any(names(dim(data)) %in% 'grid'))){
    stop("Parameter 'data' must have a spatial dimension named 'grid'.")
  }
  attractor <- Apply(data, target_dims = c('time', 'grid'),
                     fun = .proxiesattractor,
                     quanti = quanti , 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))
}

.proxiesattractor <- function(data, quanti) {
  # expected dimensions data: time and grid
  logdista <- Apply(data, target_dims =  'grid',
                    fun = function(x, y){
                      -log(colMeans((y - as.vector(x))^2))},
                    y = t(data))[[1]]
  
  #Computation of theta
  Theta <- function(logdista, quanti){
    #Compute the thheshold corresponding to the quantile
    thresh  <- quantile(logdista, quanti, na.rm = TRUE)
    logdista[which(logdista == 'Inf')] <-  NaN
    Li <- which(as.vector(logdista) > as.numeric(thresh))
    #Length of each cluster
    Ti <- diff(Li)
    N <-  length(Ti)
    q <-  1 - quanti
    Si <- Ti - 1
    Nc <- length(which(Si > 0))
    N <-  length(Ti)
    theta <- (sum(q * Si) + N + Nc - sqrt(((sum(q * Si) + N + Nc)^2) - 
             8 * Nc * sum(q * Si))) / (2 * sum(q * Si))
    #Sort the exceedances  
    logdista <- sort(logdista)
    #Find all the Peaks over Thresholds.
    findidx <- which(as.vector(logdista) > as.numeric(thresh))
    if(length(findidx) < 1) {
      stop("Parameter 'quanti' is too high for the length of the data provided.")
    }
    logextr <- logdista[findidx[[1]]:(length(logdista) - 1)]
    #The inverse of the dimension is just the average of the exceedances
    dim <- 1 /mean(as.numeric(logextr) - as.numeric(thresh))
    return(list(dim = dim, theta = theta))
  }
  names(dim(logdista)) <- c('dim1', 'dim2')
  proxies <- Apply(data = list(logdista = logdista),
                  target_dims = list('dim1'), fun = Theta, quanti = quanti)
  
  return(list(dim = proxies$dim, theta = proxies$theta))
}