rp.R

# package casnet ----
#
# Recurrence Plots ----
#
# rp functions

#' Create a Distance Matrix
#'
#' @param y1 A numeric vector or time series
#' @param y2 A numeric vector or time series for cross recurrence
#' @param emDim The embedding dimensions
#' @param emLag The embedding lag
#' @param emRad The threshold (emRad) to apply to the distance matrix to create a binary or weighted matrix. If `NULL`, an unthresholded matrix will be created (default = `NULL`)
#' @param theiler Use a `theiler` window around the main diagonal (Line of Identity/Synchronisation) to remove auto-correlations at short time-lags:
#' * `0` will include the main diagonal in all RQA measure calculations.
#' * `1` will remove the main diagonal from all RQA measure calculations.
#' * `NA` (default), will check if the matrix is symmetrical , if so, it will remove the diagonal by setting `theiler = 1` (Line of Identity, Auto-RQA), if it is not symmetrical (Line of Synchronisation, Cross-RQA) it will set `theiler = 0`.
#' * A value greater than `1` will remove that many diagonals around and including the diagonal from all RQA measure calculations. So `theiler = 2` means exclude `2` diagonals around the main diagonal, including the main diagonal itself: `[-1,0,1]`.
#' If `theiler` is a numeric vector of `length(theiler) == 2` it is possible to exclude an asymmetrical window. The values are interpreted as end points in a sequence of diagonal ID's, e.g. `theiler = c(-1,5)` will exclude `[-1,0,1,2,3,4,5]`. If `length(theiler) > 2`, the values will be considered individual diagonal ID's, e.g. `theiler = c(-3,-1,0,2,5)`, will exclude only those specific ID's. Also see the note.
#' @param to.ts Should `y1` and `y2` be converted to time series objects?
#' @param order.by If `to.ts = TRUE`, pass a vector of the same length as `y1` and `y2`. It will be used as the time index, if `NA` the vector indices will be used to represent time.
#' @param to.sparse Should sparse matrices be used?
#' @param weighted If `FALSE` a binary matrix will be returned. If `TRUE` every value larger than `emRad` will be `0`, but values smaller than `emRad` will be retained (default = `FALSE`)
#' @param weightedBy After setting values smaller than `emRad` to `0`, what should the recurrent values represent? The default is to use the state space similarity (distance/proximity) values as weights (`"si"`). Other option are `"rt"` for *recurrence time* and `"rf"` for *recurrence time frequency* (default = `"si"`)
#' @param method Distance measure to use. Any option that is valid for argument `method` of [proxy::dist()]. Type `proxy::pr_DB$get_entries()` to see a list of all the options. Common methods are: `"Euclidean", "Manhattan", "Minkowski", "Chebysev"` (or the same but shorter: `"L2","L1","Lp", "max"` distance). To use the shape based distance for phase-based recurrence use `"SBD"` (default = `"Euclidean"`)
#' @param rescaleDist Should the distance matrix be rescaled? Options are "none", "maxDist" to create a unit scale, "meanScale" to creat z-scores based on the mean distance. (default = `"none"`)
#' @param targetValue A value passed to `est_radius(...,type="fixed", targetMeasure="RR")` if `is.na(emRad)==TRUE`.
#' @param chromatic Perform a chromatic RQA. This assumes the recurring values represent the labels of an unordered categorical variable (default = `FALSE`)
#' @param returnMeasures Should the output of [rp_measures()] be returned as an attribute `"measures"` to the matrix? If `silent = FALSE` results will also be output to the console. (default = `FALSE`)
#' @param doPlot Plot the matrix by calling [rp_plot()] with default settings
#' @param doEmbed If `FALSE`, a distance matrix will be returned that is not embedded by `emDim` and `emLag` (Multidimensional RQA). If `y1` and/or `y2` are data frames, the columns will be used as the state space dimensions (default = `TRUE`)
#' @param silent Silent-ish mode
#' @param ... Any parameters to pass to [rp_plot()] if `doPlot = TRUE`
#'
#' @note The calculation of the (C)RQA measures in [casnet] can be different from other packages. For example, depending on the value of `theiler` the main diagonal can be included or excluded from the calculations, whereas some software will always include the diagonal.
#'
#' @return A (Coss-) Recurrence matrix with attributes:
#'
#' * `emdims1` and `emdims2` - A matrix of surrogate dimensions
#' * `emdims1.name` and `emdims2.name` - Names of surrogate dimensions
#' * `method` and `call` - The distance `method` used by [proxy::dist()]
#' * `weighted` - Whether a weighted matrix is returned
#' * `emDim`, `emLag` and `emRad` - The embedding parameters
#' * `AUTO` - Whether the matrix represents AUTO recurrence
#'
#'
#' @export
#'
#' @family Distance matrix operations (recurrence plot)
#'
#'
rp <- function(y1, y2 = NULL,
               emDim = 1,
               emLag = 1,
               emRad = NULL,
               theiler = NA,
               to.ts = NULL,
               order.by = NULL,
               to.sparse = TRUE,
               weighted = FALSE,
               weightedBy = "si",
               method = c("Euclidean","SBD")[1],
               rescaleDist = c("none","maxDist","meanDist")[1],
               targetValue  = .05,
               chromatic = FALSE,
               returnMeasures = FALSE,
               doPlot = FALSE,
               doEmbed = TRUE,
               silent = TRUE,
               ...){


  if(is.null(y2)){
    y2 <- y1
    attributes(y2) <- attributes(y1)
    AUTO <- TRUE
  }

  if(method == "SBD"){
    checkPkg("dtwclust")
  }

  atlist <- attributes(y1)
  if(any(names(atlist) %in% c("emDims1","emDims2","emDims1.name","emDims2.name"))){
    stop("Input is already a recurrence matrix created by 'rp()'. To manually create a binary matrix from a distance matrix use 'mat_di2bi()'. To create a weighted matrix use 'mat_di2we()'")
  }


  if(chromatic){
    if(any(is.data.frame(y1),is.data.frame(y2))){
      stop("Multidimensional Chromatic RQA has not been invented (yet)!")
    }
    uniqueID <- as.numeric_discrete(x = c(y1,y2), sortUnique = TRUE, keepNA = TRUE)
    y1  <- uniqueID[1:length(y1)]
    y2  <- uniqueID[(length(y2)+1):length(uniqueID)]

    chromaDims <- list(y1 = uniqueID[1:length(y1)],
                       y2 = uniqueID[(length(y2)+1):length(uniqueID)])

    chromaNames <- unique(uniqueID)
    nList <- list()
    for(n in seq_along(chromaNames)){
      nList[[n]] <- names(uniqueID)[which(uniqueID%in%chromaNames[n])]
    }
    names(chromaNames) <- unique(unlist(nList))

  } else {
    chromaNames <- NA
    chromaDims  <- NA
  }

  if(!is.data.frame(y1)){
    id <- deparse(substitute(y1))
    y1 <- as.data.frame(y1)
    if(length(id)==NCOL(y1)){
      colnames(y1) <- id
    }
  }

  if(!is.data.frame(y2)){
    id <- deparse(substitute(y2))
    y2 <- as.data.frame(y2)
    if(length(id)==NCOL(y2)){
      colnames(y2) <- id
    }
  }

  if(any(!doEmbed,chromatic)){
    emDim <- 1
    emLag <- 0
    if(chromatic){
      emRad <- 0
    }
  }

  et1 <- ts_embed(y1, emDim=emDim, emLag=emLag, silent = silent)
  et2 <- ts_embed(y2, emDim=emDim, emLag=emLag, silent = silent)

  dist_method <- return_error(proxy::pr_DB$get_entry(method))
  if("error"%in%class(dist_method$value)){
    stop("Unknown distance metric!\nUse proxy::pr_DB$get_entries() to see a list of valid options.")
  } else {

    dmat <- proxy::dist(x = et1,
                        y = et2,
                        method = method,
                        diag = TRUE)

    if(method == "SBD"){
      message("Make sure all phase space dimensions are on the same scale when using Shape Based Distance.")
      dmat <- Matrix::as.matrix(dmat)
    }
  }

  # Remove proxy class
  dmat <- unclass(dmat)

  # Check a time index was requested
  if(!is.null(to.ts)){
    if(is.null(order.by)){
      order.by <- lubridate::as_datetime(1:NROW(dmat), origin = lubridate::ymd_hms(Sys.time()))
    }
    dmat <-  switch(to.ts,
                    "xts" =  xts::xts(dmat, order.by = lubridate::as_datetime(order.by)),
                    "zoo" =  zoo::zoo(dmat, order.by = lubridate::as_datetime(order.by))
    )
  }
  # Add names if ordered by time vector
  if(!is.null(order.by)){
    colnames(dmat) <- paste(order.by)
    rownames(dmat) <- paste(order.by)
  }

  if(!chromatic){
    if(rescaleDist=="maxDist"){dmat <- dmat/max(dmat,na.rm = TRUE)}
    if(rescaleDist=="meanDist"){dmat <- (dmat-mean(dmat, na.rm=TRUE))/sd(dmat, na.rm = TRUE)}
  }

  if(!is.null(emRad)){
    if(is.na(emRad)){
      suppressMessages(emRad <- est_radius(RM = dmat, emDim = emDim, emLag = emLag, targetValue = targetValue, silent = TRUE))
      if(emRad$Converged){
        emRad <- emRad$Radius
      } else {
        emRad <- est_radius(RM = dmat, emDim = emDim, emLag = emLag, targetValue = targetValue, tol = 0.1, silent = TRUE, radiusOnFail = "percentile")$Radius
      }
    }
    if(chromatic){
      dmat <- mat_di2ch(dmat, y = et1, emRad = emRad, convMat = to.sparse)
    } else {
      if(weighted){
        dmat <- mat_di2we(dmat, emRad = emRad, convMat = to.sparse)
      } else {
        dmat <- mat_di2bi(dmat, emRad = emRad, convMat = to.sparse)
      }
    }
  }

  dmat <- rp_checkfix(dmat, checkAUTO = TRUE, fixAUTO = TRUE, checkSPARSE = TRUE, fixSPARSE = TRUE, checkS4 = TRUE, fixS4 = TRUE)
  suppressMessages(dmat <- setTheiler(RM = dmat, theiler = theiler, chromatic = chromatic))
  theiler <- attributes(dmat)$theiler

  if(returnMeasures){
    rpOut <-  rp_measures(RM = dmat,
                          emRad = emRad%00%NA,
                          silent = silent,
                          theiler = theiler,
                          chromatic = chromatic)
  } else {
    rpOut <- NA
  }

  if(to.sparse){
    attributes(dmat)$emDims1  <- et1
    attributes(dmat)$emDims2  <- et2
    attributes(dmat)$emDims1.name <- colnames(y1)
    attributes(dmat)$emDims2.name <- colnames(y2)
    attributes(dmat)$embedded <- doEmbed
    attributes(dmat)$emLag <- emLag
    attributes(dmat)$emDim <- emDim
    attributes(dmat)$emRad <- emRad%00%NA
    attributes(dmat)$measures <- rpOut
    attributes(dmat)$weighted <- weighted
    attributes(dmat)$weightedBy <- weightedBy
    attributes(dmat)$chromatic <- chromatic
    attributes(dmat)$chromaNames <- chromaNames
    attributes(dmat)$chromaDims <- chromaDims
    attributes(dmat)$theiler <- theiler
  } else {
    attr(dmat,"emDims1") <- et1
    attr(dmat,"emDims2") <- et2
    attr(dmat,"emDims1.name") <- colnames(y1)
    attr(dmat,"emDims2.name") <- colnames(y2)
    attr(dmat,"weighted") <- weighted
    attr(dmat,"embedded") <- doEmbed
    attr(dmat,"emLag") <- emLag
    attr(dmat,"emDim") <- emDim
    attr(dmat,"emRad") <- emRad%00%NA
    attr(dmat,"measures") <- rpOut
    attr(dmat,"weighted") <- weighted
    attr(dmat,"weightedBy") <- weightedBy
    attr(dmat,"chromatic") <- chromatic
    attr(dmat,"chromaNames") <- chromaNames
    attr(dmat,"chromaDims") <- chromaDims
    attr(dmat,"theiler") <- theiler
  }

  if(is.null(attributes(dmat))){stop("lost attributes")}

  if(doPlot){

    if(rp_size(RM = dmat)$rp_size_total>1024){
      warning("Plotting will take a long time. Consider running rp_plot(..., courseGrain = TRUE)", immediate. = TRUE)
    }
    rp_plot(RM = dmat, courseGrain = FALSE)

    # if(...length()>0){
    #   dotArgs <- list(...)
    #   nameOK  <- names(dotArgs)%in%methods::formalArgs(rp_plot)
    #   dotArgs <- dotArgs[nameOk]
    # } else {
    #   # Plot with defaults
    #   dotArgs  <- formals(rp_plot)
    #   #nameOk   <- rep(TRUE,length(dotArgs))
    # }
    #
    # if(!"RM"%in%names(dotArgs)|length(dotArgs$RM)==1){
    #   dotArgs$RM <- dmat
    # }
    # if(!is.na(emRad%00%NA)&(!"radiusValue"%in%names(dotArgs))){
    #   dotArgs$radiusValue <- emRad
    # }
    # do.call(rp_plot, as.list(dotArgs))
  }

  return(dmat)
}


#' Get (C)RQA measures based on a binary matrix
#'
#' A zoo of measures based on singular recurrent points, diagonal, vertical and horizontal line structures (anisotropic) will be caluclated.
#'
#' @aliases crqa_rp
#' @inheritParams rp
#' @param RM A distance matrix (set `emRad = NA` to estimate a radius), or a matrix of zeroes and ones
#' @param emRad Threshold for distance value that counts as a recurrence (ignored is `RM` is a binary matrix)
#' @param DLmin Minimal diagonal line length (default = `2`)
#' @param VLmin Minimal vertical line length (default = `2`)
#' @param HLmin Minimal horizontal line length (default = `2`)
#' @param DLmax Maximal diagonal line length (default = length of diagonal -1)
#' @param VLmax Maximal vertical line length (default = length of diagonal -1)
#' @param HLmax Maximal horizontal line length (default = length of diagonal -1)
#' @param AUTO Auto-recurrence? (default = `FALSE`)
#' @param chromatic Force chromatic RQA? If `NA` the value of the `RM` attribute `"chromatic"` will be used, if present (default = `NA`)
#' @param anisotropyHV Return anisotropy ratio measures based on Horizontal and Vertical lines. The ratios are calculated as `(horizontal - vertical) / (horizontal + vertical)`. So a value of 0 means no anisotropy, negative ratios indicate the measures based on vertical lines had  higher values, positive ratios indicate the measures based on horizontal lines had higher values  (default = `FALSE`)
#' @param asymmetryUL Return asymmetry ratio measures based on Upper and Lower triangles. The ratios are calculated as `(upper - lower) / (upper + lower)`. So a value of 0 means no asymmetry, negative ratios indicate the measures based on the lower triangle had the higher values, positive ratios indicate measures based on the upper triangle had higher values (default = `FALSE`)
#' @param returnUL Return the (C)RQA values for the upper and lower triangle on which asymmetry ratio calculations are based? (default = `FALSE`)
#' @param recurrenceTimes Return measures based on 'white lines', the recurrence times (default = `FALSE`)
#' @param matrices Return matrices? (default = `FALSE`)
#' @param Nboot How many bootstrap replications? (default = `NULL`)
#' @param CL Confidence limit for bootstrap results (default = `.95`)
#' @param targetValue A value passed to `est_radius(...,type="fixed", targetMeasure="RR", tol = .2)` if `is.na(emRad)==TRUE`, it will estimate a radius (default = `.05`).
#' @param use_rqa_par Use function [rqa_par()] if matrix dims >  (default = `FALSE`)
#' @param silent Do not display any messages (default = `TRUE`)
#'
#' @return A list object containing (C)RQA measures (and matrices if requested)
#'
#' @export
#'
#' @family Recurrence Quantification Analysis
#'
#'
rp_measures <- function(RM,
                        emRad = NA,
                        DLmin = 2,
                        VLmin = 2,
                        HLmin = 2,
                        DLmax = length(Matrix::diag(RM)),
                        VLmax = length(Matrix::diag(RM)),
                        HLmax = length(Matrix::diag(RM)),
                        AUTO      = NULL,
                        theiler   = NA,
                        chromatic = NA,
                        anisotropyHV = FALSE,
                        asymmetryUL = FALSE,
                        returnUL = FALSE,
                        recurrenceTimes = FALSE,
                        matrices  = FALSE,
                        Nboot     = NULL,
                        CL        = .95,
                        targetValue = .05,
                        use_rqa_par  = FALSE,
                        silent     = TRUE){


  # Input should be a distance matrix, or a matrix of zeroes and ones with emRad = NULL, output is a list
  # Fred Hasselman - August 2013

  # check auto-recurrence and make sure Matrix has sparse triplet representation
  RM <- rp_checkfix(RM, checkAUTO = TRUE, fixAUTO = TRUE, checkTSPARSE = TRUE, fixTSPARSE = TRUE)

  if(is.null(AUTO)){
    AUTO <- attr(RM,"AUTO")
  }

  #require(parallel)
  if(is.na(theiler)){
    if(AUTO){
      theiler <- 1
      message("Auto-RQA, not including diagonal, theiler set to 1...")
    } else {
      theiler <- 0
      message("Cross-RQA, including diagonal, theiler set to 0...")
    }
  }




  if(is.na(chromatic)){
    if(!is.null(attr(RM,"chromatic"))){
      chromatic <- attr(RM,"chromatic")
    } else {
      chromatic <- FALSE
    }
  }

  if(chromatic){
    if(abs(diff(range(RM)))==1){
      message("There is only 1 ctategory, chromatic (C)RQA is not sensible. Setting chromatic to FALSE...")
      chromatic <- FALSE
    }
  }

  if(is.na(emRad)){
    if(!is.null(attributes(RM)$emRad)){
      emRad <- attributes(RM)$emRad
    } else {
      # Check for attributes
      if(is.na(targetValue)){
        targetValue <- .05
      }
      if(!is.null(attributes(RM)$emDim)){
        emDim <- attributes(RM)$emDim
      } else {
        emdim <- 1
      }
      if(!is.null(attributes(RM)$emLag)){
        emLag <- attributes(RM)$emLag
      } else {
        emLag <- 1
      }

      emRad <- est_radius(RM = RM, emDim = emDim, emLag = emLag, targetValue = targetValue, tol = .2, radiusOnFail = "percentile", silent = silent)
      if(emRad$Converged){
        emRad <- emRad$Radius
      } else {
        emRad <- stats::sd(RM,na.rm = TRUE)
      }
    } # not emRad attr
  }

  if(use_rqa_par){
    if(max(dim(RM))>2046){

      y1 <- attributes(RM)$emDims1
      y2 <- attributes(RM)$emDims2

      out <- rqa_par(y1 = y1,
                     y2 = y2,
                     emRad = emRad,
                     DLmin = DLmin,
                     VLmin = VLmin,
                     HLmin = HLmin,
                     DLmax = DLmax,
                     VLmax = VLmax,
                     HLmax = HLmax,
                     theiler = theiler,
                     AUTO  = AUTO,
                     chromatic = chromatic,
                     anisotropyHV = anisotropyHV,
                     asymmetryUL = asymmetryUL,
                     recurrenceTimes = recurrenceTimes,
                     weighted = FALSE,
                     weightedBy = "si",
                     method = method,
                     rescaleDist = c("none","maxDist","meanDist")[1],
                     targetValue  = targetValue,
                     doEmbed = FALSE,
                     silent = silent,
                     ...)

      return(out)
    }
  } else {


  #uval <- unique(as.vector(RM))
  if(!all(as.vector(RM)==0|as.vector(RM)==1)){

    if(chromatic){

      # prepare data
      if(any(grepl("Matrix",class(RM)))){
        RM     <- rp_checkfix(RM, checkTSPARSE = TRUE, fixTSPARSE = TRUE)
        #meltRP <- data.frame(Var1 = (RM@i+1), Var2 = (RM@j+1), value = as.numeric(RM@x))
        meltRP <-  mat_mat2ind(Matrix::as.matrix(RM))
      } else {
        meltRP <-  mat_mat2ind(as.matrix(RM))
      }

      if(!is.null(attr(RM,"chromaNames"))){
        chromaNumbers <- attr(RM,"chromaNames")
        Cind <- sort(chromaNumbers, index.return = TRUE)
        chromaNames <- names(chromaNumbers)[Cind$ix]
        chromaNumbers <- chromaNumbers[Cind$ix]
      } else {
        chromas <- meltRP %>% dplyr::filter(.data$value!=0)
        chromaNumbers <- as.numeric_discrete(sort(unique(chromas$value)))
      }

      chromaNumbers <- stats::na.exclude(chromaNumbers)

      rm(meltRP)

      if(NROW(chromaNumbers)>=(NROW(RM)/2)){
        warning(paste("Chromatic RQA will have a large number of categories:",NROW(chromaNumbers)))
      }

      chromaList <- list()
      matrixList <- list()

      for(i in seq_along(chromaNumbers)){

        RMtmp <- RM
        RMtmp[RM!=i] <- 0
        RMtmp[RM==i] <- 1

        suppressMessages(tmpOut <- rp_calc(RMtmp,
                                           emRad = emRad,
                                           DLmin = DLmin,
                                           VLmin = VLmin,
                                           HLmin = HLmin,
                                           DLmax = DLmax,
                                           VLmax = VLmax,
                                           HLmax = HLmax,
                                           theiler = theiler,
                                           AUTO  = AUTO,
                                           chromatic = chromatic,
                                           anisotropyHV = anisotropyHV,
                                           asymmetryUL = asymmetryUL,
                                           recurrenceTimes = recurrenceTimes,
                                           matrices  = matrices)
        )

        if(matrices){
          chromaList[[i]] <- tmpOut$crqaMeasures
          matrixList[[i]] <- tmpOut$crqaMatrices
        } else {
          chromaList[[i]] <- tmpOut
        }
        rm(RMtmp, tmpOut)
      }

      names(chromaList) <- names(chromaNumbers)

      if(is.na(length(chromaList)%00%NA)){
        stop("Chromalist is empty")
      }

      if(matrices){
        names(matrixList) <- names(chromaNumbers)
        out <- list(crqaMeasures = plyr::ldply(chromaList, .id = "chroma"),
                    crqaMatrices =  matrixList)
      } else {
        out <- plyr::ldply(chromaList, .id = "chroma")
      }

    } else { #IF CHROMATIC
      if(!is.null(emRad)){
        RM <- mat_di2bi(RM,emRad)
      } else {
        if(!chromatic){
          stop("Expecting a binary (0,1) matrix.\nUse 'est_radius()', or set 'chromatic = TRUE'")
        }
      }
    }
  }
  #rm(RM)

  if(!chromatic){

      suppressMessages(out <- rp_calc(RM,
                                      emRad = emRad,
                                      DLmin = DLmin,
                                      VLmin = VLmin,
                                      HLmin = HLmin,
                                      DLmax = DLmax,
                                      VLmax = VLmax,
                                      HLmax = HLmax,
                                      theiler = theiler,
                                      AUTO  = AUTO,
                                      chromatic = chromatic,
                                      anisotropyHV = anisotropyHV,
                                      asymmetryUL = asymmetryUL,
                                      recurrenceTimes = recurrenceTimes))
  }

    # } else {
    #
    #   suppressMessages(out <- rp_calc(RM,
    #                                   emRad = emRad,
    #                                   DLmin = DLmin,
    #                                   VLmin = VLmin,
    #                                   HLmin = HLmin,
    #                                   DLmax = DLmax,
    #                                   VLmax = VLmax,
    #                                   HLmax = HLmax,
    #                                   theiler = theiler,
    #                                   AUTO  = AUTO,
    #                                   chromatic = chromatic,
    #                                   anisotropyHV = anisotropyHV,
    #                                   asymmetryUL = asymmetryUL,
    #                                   recurrenceTimes = recurrenceTimes,
    #                                   matrices  = matrices))
    # }
  #}

  } # if !use_rqa_par


  if(matrices){
      tab <- out$crqaMeasures
    } else {
      tab <- out
    }

  if(chromatic){
    Nrows <- length(chromaNames)
  } else {
    Nrows <- 1
  }

    outTable <- list(`Global Measures` = data.frame(Global = "Matrix",
                                                    `Max points` = tab$RP_max,
                                                    `N points` = tab$RP_N,
                                                    RR = tab$RR,
                                                    Singular = tab$SING_N%00%NA,
                                                    Divergence = tab$DIV_dl,
                                                    Repetitiveness = tab$REP_av),
                     `Line-based Measures` = data.frame(`Lines`    = rep(c("Diagonal", "Vertical", "Horizontal", "V+H Total"), each = Nrows),
                                                        `N lines`  = c(tab$N_dl,tab$N_vl, tab$N_hl, tab$N_hv),
                                                        `N points` = c(tab$N_dlp%00%NA,tab$N_vlp%00%NA,tab$N_hlp%00%NA, tab$N_hvp%00%NA),
                                                        `Measure`  = rep(c("DET","V LAM","H LAM", "V+H LAM"), each = Nrows),
                                                        `Rate`     = c(tab$DET, tab$LAM_vl, tab$LAM_hl, tab$LAM_hv),
                                                        `Mean`     = c(tab$MEAN_dl, tab$TT_vl, tab$TT_vl, tab$TT_hv),
                                                        `Max.`     = c(tab$MAX_dl, tab$MAX_vl, tab$MAX_hl, tab$MAX_hv),
                                                        `ENT`      = c(tab$ENT_dl, tab$ENT_vl, tab$ENT_hl, tab$ENT_hv),
                                                        `ENT_rel` = c(tab$ENTrel_dl, tab$ENTrel_vl, tab$ENTrel_hl, tab$ENTrel_hv),
                                                        `CoV`   = c(tab$CoV_dl, tab$CoV_vl, tab$CoV_hl, tab$CoV_hv)))

    if(chromatic){
      rownames(outTable$`Global Measures`) <- chromaNames
      rownames(outTable$`Line-based Measures`) <- paste(1:NROW(outTable$`Line-based Measures`), rep(chromaNames, 4))
    }

    if(anisotropyHV){

      outTable$`Horizontal/Vertical line anisotropy` <-  data.frame(`Ratio` = "H/V line measures",
                                                                    `N lines`  = tab$ratios.hv.Nlines_ani,
                                                                    `N points` = tab$ratios.hv.N_hlp%00%NA/tab$ratios.hv.N_vlp%00%NA,
                                                                    `Measure` = "LAM",
                                                                    `Rate`    = tab$ratios.hv.LAM_ani,
                                                                    `Mean`    = tab$ratios.hv.MEAN_ani,
                                                                    `Max`     = tab$ratios.hv.MAX_ani,
                                                                    `ENT`     = tab$ratios.hv.ENT_ani)
      if(chromatic){
      rownames(outTable$`Horizontal/Vertical line anisotropy`) <- chromaNames
      }
    }

    if(asymmetryUL){

      outTable$`Upper/Lower triangle asymmetry` <- list(
        `Global Measures` =  data.frame(`Global Ratio` = "U/L of points",
                                        `N points` = tab$ratios.ul.Npoints_asym,
                                        RR = tab$ratios.ul.RR_asym,
                                        Singular = tab$ratios.ul.SING_N_asym,
                                        Divergence = tab$ratios.ul.DIV_asym,
                                        Repetetiveness = tab$ratios.ul.REP_asym),
        `Line-based Measures` = data.frame(
          `Line ratio` = rep(c("D lines", "V lines", "H lines"), each = Nrows),
          `N lines`  = c(tab$ratios.ul.NDlines_asym,
                         tab$ratios.ul.NVlines_asym,
                         tab$ratios.ul.NHlines_asym),
          `N points` = c(tab$ratios.ul.NDpoints_asym,
                         tab$ratios.ul.NVpoints_asym,
                         tab$ratios.ul.NHpoints_asym),
          `Measure` = rep(c("DET","V LAM", "H LAM"), each = Nrows),
          `Rate`    = c(tab$ratios.ul.DET_asym, tab$ratios.ul.LAM_vl_asym, tab$ratios.ul.LAM_hl_asym),
          `Mean`    = c(tab$ratios.ul.MEAN_dl_asym, tab$ratios.ul.MEAN_vl_asym, tab$ratios.ul.MEAN_hl_asym),
          `Max`     = c(tab$ratios.ul.MAX_dl_asym, tab$ratios.ul.MAX_vl_asym, tab$ratios.ul.MAX_hl_asym),
          `ENT`     = c(tab$ratios.ul.ENT_dl_asym, tab$ratios.ul.ENT_vl_asym, tab$ratios.ul.ENT_hl_asym))
      )

      if(chromatic){
        id <- which(names(outTable)%in%"Upper/Lower triangle asymmetry")
        rownames(outTable[[id]]$`Global Measures`) <- chromaNames
        rownames(outTable[[id]]$`Line-based Measures`) <- paste(1:NROW(outTable[[id]]$`Line-based Measures`), rep(chromaNames, 3))
      }
    }

    if(!silent){
      cat("\n~~~o~~o~~casnet~~o~~o~~~\n")
      if(chromatic){
        cat(paste0("Chromatic RQA with categories: ", paste(chromaNames, collapse = ", "),"\n"))
      }
      cat(" Global Measures\n")
      print(format(outTable$`Global Measures`, digits = 3))
      cat(paste("\n\n Line-based Measures\n"))
      print(format(outTable$`Line-based Measures`, digits = 3))

      if(anisotropyHV){
        id <- which(names(outTable)%in%"Horizontal/Vertical line anisotropy")
        cat(paste0("\n\n",names(outTable)[id],"\n\n"))
        print(format(outTable$`Horizontal/Vertical line anisotropy`, digits = 3))
      }

      if(asymmetryUL){
        id <- which(names(outTable)%in%"Upper/Lower triangle asymmetry")
        cat(paste0("\n\n",names(outTable)[id],"\n\n"))
        cat(" Global Measures\n")
        print(format(outTable[[id]]$`Global Measures`, digits = 3))
        cat(paste("\n\n Line-based Measures\n"))
        print(format(outTable[[id]]$`Line-based Measures`, digits = 3))
      }
      cat("\n~~~o~~o~~casnet~~o~~o~~~\n")
    }

    attr(out,"measuresTable") <- outTable

    return(invisible(out))
}

  # tb <- c("| Max rec. points (post theiler) | N rec. points | RR | Singular points | Divergence | Repetiteveness | Anisotropy |",
  #         "|--------------------------------|---------------|----|-----------------|------------|----------------|------------|".
  #         )
  #
  #cat(c("Global Recurrence Measures\n",paste0(tb,"\n")))


  #   if(is.null(Nboot)){Nboot = 1}
  #
  #   out <- rp_measures_measures(RM,
  #                            emRad= emRad,
  #                            DLmin = DLmin,
  #                            VLmin = VLmin,
  #                            HLmin = HLmin,
  #                            DLmax = DLmax,
  #                            VLmax = VLmax,
  #                            HLmax = HLmax,
  #                            AUTO  = AUTO,
  #                            chromatic = chromatic,
  #                            matrices  = matrices,
  #                            doHalf = doHalf,
  #                            Nboot  = Nboot,
  #                            CL     = CL)
  #
  #   dfori <- dplyr::gather(out, key = measure, value = value)
  #
  #   col.ind <- dplyr::tbl_df(index(RM))
  #   row.ind <- dplyr::tbl_df(sample(index(RM),size=nrow(RM)))
  #
  #   if(Nboot>1){cat(paste0("Bootstrapping Recurrence Matrix... ",Nboot," iterations.\n"))
  #     bootout <- col.ind  %>%
  #       bootstrap(Nboot) %>%
  #       do(rp_measures_measures(RM[row.ind,unlist(.)],
  #                            emRad= emRad,
  #                            DLmin = DLmin,
  #                            VLmin = VLmin,
  #                            HLmin = HLmin,
  #                            DLmax = DLmax,
  #                            VLmax = VLmax,
  #                            HLmax = HLmax,
  #                            AUTO  = AUTO,
  #                            chromatic = chromatic,
  #                            matrices  = matrices,
  #                            doHalf = doHalf))
  #
  #     dfrepl <- gather(bootout, key = measure, value = value, -replicate)
  #
  #     if(length(CL)==1){
  #       ci.lo <- (1-CL)/2
  #       ci.hi <- CL + ci.lo
  #     } else {
  #       ci.lo <- CL[1]
  #       ci.hi <- CL[2]
  #     }
  #
  #     rqout <-  dfrepl %>% dplyr::group_by(measure) %>%
  #       summarize(
  #         val     = NA,
  #         ci.low  = stats::quantile(value, ci.lo, na.rm = TRUE),
  #         ci.high = stats::quantile(value, ci.hi, na.rm = TRUE),
  #         mean    = mean(value, na.rm = TRUE),
  #         sd      = stats::sd(value, na.rm = TRUE),
  #         var     = stats::var(value, na.rm = TRUE),
  #         N       = n(),
  #         se      = stats::sd(value, na.rm = TRUE)/sqrt(n()),
  #         median  = mean(value, na.rm = TRUE),
  #         mad     = stats::mad(value, na.rm = TRUE)
  #       )
  #
  #     for(m in 1:nrow(rqout)){
  #       rqout$val[rqout$measure%in%dfori$measure[m]] <- dfori$value[m]
  #     }
  #   } else {
  #     rqout <- dfori
  #   }


#' Get (C)RQA measures from a Recurrence Matrix
#'
#' Use `rp_measures`
#'
#' @inheritParams rp_measures
#'
#' @family Recurrence Quantification Analysis
#'
#' @export
#'
#' @keywords internal
#'
rp_measures_main <- function(RM,
                             emRad = NULL,
                             DLmin = 2,
                             VLmin = 2,
                             HLmin = 2,
                             DLmax = length(Matrix::diag(RM)),
                             VLmax = length(Matrix::diag(RM)),
                             HLmax = length(Matrix::diag(RM)),
                             AUTO      = NULL,
                             theiler = NULL,
                             chromatic = FALSE,
                             matrices  = FALSE,
                             doHalf    = FALSE,
                             Nboot     = NULL,
                             CL        = .95,
                             doParallel = FALSE){


  if(is.null(Nboot)){
    Nboot <- 1
  }

  if(Nboot>1|doParallel){
    checkPkg("parallel")
    mc.cores <- parallel::detectCores()-1
    if(Nboot<mc.cores) mc.cores <- Nboot
  }

  NRows <- NROW(RM)
  NCols <- NCOL(RM)

  if(doParallel){

    tstart <- Sys.time()
    cl  <- parallel::makeCluster(mc.cores)
    out <- parallel::mclapply(1, function(i){
      #rp_prep(matrix(RM[ceiling(NCols*NRows*stats::runif(NCols*NRows))], ncol=NCols, nrow = NRows),
      rp_prep(RP = RM,
              emRad= emRad,
              DLmin = DLmin,
              VLmin = VLmin,
              HLmin = HLmin,
              DLmax = DLmax,
              VLmax = VLmax,
              HLmax = HLmax,
              AUTO  = AUTO,
              chromatic = chromatic,
              matrices  = matrices,
              doHalf    = doHalf)
    },
    mc.cores = mc.cores
    )
    parallel::stopCluster(cl)
    tend  <- Sys.time()
  } else {
    out <- rp_prep(RP = RM,
                   emRad= emRad,
                   DLmin = DLmin,
                   VLmin = VLmin,
                   HLmin = HLmin,
                   DLmax = DLmax,
                   VLmax = VLmax,
                   HLmax = HLmax,
                   AUTO  = AUTO,
                   chromatic = chromatic,
                   matrices  = matrices,
                   doHalf    = doHalf)
  }


  dfori <- tidyr::gather(as.data.frame(out), key = "measure", value = "value")

  if(Nboot>1){

    cat(paste0("Bootstrapping Recurrence Matrix... ",Nboot," iterations...\n"))
    cat(paste0("Estimated duration: ", round((difftime(tend,tstart, units = "mins")*Nboot)/max((round(mc.cores/2)-1),1), digits=1)," min.\n"))

    tstart <-Sys.time()
    cl      <- parallel::makeCluster(mc.cores)
    bootOut <-  parallel::mclapply(1:Nboot, function(i){
      replicate <- as.data.frame(rp_prep(matrix(RM[ceiling(NCols*NRows*stats::runif(NCols*NRows))],
                                                ncol=NCols, nrow = NRows),
                                         emRad= emRad,
                                         DLmin = DLmin,
                                         VLmin = VLmin,
                                         HLmin = HLmin,
                                         DLmax = DLmax,
                                         VLmax = VLmax,
                                         HLmax = HLmax,
                                         AUTO  = AUTO,
                                         chromatic = chromatic,
                                         matrices  = matrices,
                                         doHalf    = doHalf))
      replicate$replicate = i
      return(replicate)
    },
    mc.cores = mc.cores
    )
    parallel::stopCluster(cl)
    tend <- Sys.time()
    cat(paste0("Actual duration: ", round(difftime(tend,tstart, units = "mins"), digits=1)," min.\n"))

    dfrepl <- plyr::ldply(bootOut)
    dfrepl <- tidyr::gather(dfrepl, key = "measure", value = "value", -replicate)

    if(length(CL)==1){
      ci.lo <- (1-CL)/2
      ci.hi <- CL + ci.lo
    } else {
      ci.lo <- CL[1]
      ci.hi <- CL[2]
    }

    rqout <-  dfrepl %>%
      dplyr::group_by(dfrepl$measure) %>%
      dplyr::summarise(
        val          = NA,
        ci.lower     = stats::quantile(dfrepl$value%00%NA, ci.lo, na.rm = TRUE),
        ci.upper     = stats::quantile(dfrepl$value%00%NA, ci.hi, na.rm = TRUE),
        BOOTmean     = mean(dfrepl$value%00%NA, na.rm = TRUE),
        BOOTsd       = stats::sd(dfrepl$value%00%NA, na.rm = TRUE),
        BOOTse       = stats::sd(dfrepl$value%00%NA, na.rm = TRUE)/sqrt(Nboot),
        BOOTvar      = stats::var(dfrepl$value%00%NA, na.rm = TRUE),
        BOOTmedian   = stats::median(dfrepl$value%00%NA, na.rm = TRUE),
        BOOTmad      = stats::mad(dfrepl$value%00%NA, na.rm = TRUE),
        BOOTn        = Nboot
      )

    for(m in 1:nrow(rqout)){
      rqout$val[rqout$measure%in%dfori$measure[m]] <- dfori$value[m]

    }
  } else {
    rqout <- dfori %>% tidyr::spread(dfori$measure,dfori$value)
  }
  return(rqout)
}




#' Diagonal Recurrence Profile
#'
#' @aliases crqa_diagprofile
#' @inheritParams rp_measures
#' @param RM A binary recurrence matrix
#' @param diagWin Window around the line of synchrony
#' @param xname Label for x-axis
#' @param yname Label for y-axis
#' @param doShuffle Should a shuffled baseline be calculated (default = `FALSE`)
#' @param y1 The original `y1` time series
#' @param y2 The original `y2` time series
#' @param shuffleWhich Which of the time series should be shuffled: 'y1' or 'y2'? (default = 'y2')
#' @param Nshuffle How many shuffled versions to make up the baseline? The default is `19`, which is the minimum for a one-sided surrogate test.
#' @param AUTO Auto-recurrence? (default = `FALSE`)
#' @param doEmbed If `doShuffle = TRUE`, should the data in y1 and y2 be considered embedded time series? The temporal order of all columns in `y2` will be randomly shuffled in the same way, keeping coordinates together (default = `FALSE`)
#' @param chromatic Force chromatic RQA? (default = `FALSE`)
#' @param matrices Return matrices? (default = `FALSE`)
#' @param doPlot Plot (default = `TRUE`)
#' @param returnOnlyPlot Don't plot to graphics device, but do return the plot (default = `FALSE`)
#'
#' @return A plot and/or the data for the plot
#'
#' @export
#'
rp_diagProfile <- function(RM,
                           diagWin = NULL,
                           xname = "X-axis",
                           yname = "Y-axis",
                           theiler = 0,
                           DLmin = 2,
                           VLmin = 2,
                           HLmin = 2,
                           DLmax = length(Matrix::diag(RM)),
                           VLmax = length(Matrix::diag(RM)),
                           HLmax = length(Matrix::diag(RM)),
                           doShuffle = FALSE,
                           y1        = NULL,
                           y2        = NULL,
                           shuffleWhich = "y1",
                           Nshuffle  = 19,
                           doEmbed   = TRUE,
                           AUTO      = NULL,
                           chromatic = FALSE,
                           matrices  = FALSE,
                           doPlot    = TRUE,
                           returnOnlyPlot = FALSE){


  # crqa_results_xy <– crqa(ts1 = lorData$x, ts2 = lorData$y, delay = 9, embed = 4, rescale = 2, radius = 20, normalize = 2, mindiagline = 2, minvertline = 2, tw = 0, whiteline = FALSE, recpt = FALSE, side = “both”) # compute cross-recurrence plot
  #

  if(!all(as.vector(RM)==0|as.vector(RM)==1)){
    stop("Expecting a binary (0,1) matrix.")
  }

  if(is.null(AUTO)){
    checkAUTO <- TRUE
    fixAUTO <- TRUE
  } else {
    checkAUTO <- AUTO
    fixAUTO <- AUTO
  }
  # check auto-recurrence and make sure Matrix has sparse triplet representation
  RM <- rp_checkfix(RM, checkAUTO = checkAUTO, fixAUTO = checkAUTO, checkTSPARSE = TRUE, fixTSPARSE = TRUE)

  if(is.null(AUTO)){
    AUTO <- attr(RM,"AUTO")
  }


  RM <- setTheiler(RM = RM, theiler = theiler, silent = TRUE)

  # if(is.na(theiler%00%NA)){
  #   if(!is.null(attributes(RM)$theiler)){
  #     if(attributes(RM)$theiler>0){
  #       message(paste0("Value found in attribute 'theiler'... assuming a theiler window of size:",attributes(RM)$theiler,"was already removed."))
  #     }
  #   }
  #   theiler <- 0
  # } else {
  #   if(theiler < 0){
  #     theiler <- 0
  #   }
  # }

  if(is.numeric(diagWin)){
    if(diagWin>0){
      diagWin <- seq(-diagWin,diagWin)
    } else {
      diagWin <- seq(-1*NCOL(RM),NCOL(RM))
    }
  } else {
    stop("Diagonal window must be 1 positive integer!!")
  }

  if(doShuffle){
    if(any(is.null(y1),is.null(y2))){
      stop("Need time series (y1 and y2) in order to do the shuffle!!")
    } else {
      cat("Calculating diagonal recurrence profiles... \n")

      if(!is.data.frame(y1)){
        y1 <- data.frame(y1)
      }
      if(!is.data.frame(y2)){
        y2 <- data.frame(y2)
      }

        if(shuffleWhich %in% "y1"){
         TSrnd <- plyr::llply(0:Nshuffle, function(r){
           if(r==0){
             return(y1)
           } else {
           return(data.frame(y1[sample.int(n = NROW(y1), size = NROW(y1)),]))
             }
         })
        }
        if(shuffleWhich %in% "y2"){
          TSrnd <- plyr::llply(0:Nshuffle, function(r){
            if(r==0){
              return(y2)
            } else {
            return(data.frame(y2[sample.int(n = NROW(y2), size = NROW(y2)),]))
            }
        })
       }

      # else {
      #   if(shuffleWhich %in% "y1"){
      #   TSrnd <- tseries::surrogate(x=y2, ns=Nshuffle, fft=FALSE, amplitude = FALSE)
      #   if(Nshuffle==1){
      #     TSrnd <- matrix(TSrnd,ncol=1)
      # }
    }
    emDim <- attr(RM,"emDim")
    emLag <- attr(RM,"emLag")
    emRad <- attr(RM,"emRad")
  } else {
    Nshuffle <- 0
    TSrnd<- list(RM)
  }

  #out <- vector(mode = "list", length = Nshuffle+1)

  out <- plyr::llply(TSrnd, function(r){

    if(doShuffle){
      if(shuffleWhich %in% "y1"){
        RMtmp <- rp(y1 = r, y2 = y2, emDim = emDim, emLag = emLag, emRad = emRad, to.sparse = TRUE)
      } else {
        RMtmp <- rp(y1 = y1, y2 = r, emDim = emDim, emLag = emLag, emRad = emRad, to.sparse = TRUE)
      }
    } else {
      RMtmp <- r
    }

    B <- rp_nzdiags(RMtmp, removeNZ = FALSE, d = diagWin)
    rm(RMtmp)

    diagID <- 1:NCOL(B)
    names(diagID) <- colnames(B)
    if(length(diagWin)<NCOL(B)){
      cID <- which(colnames(B)%in%diagWin)
      B <- B[,cID]
      diagID <- seq_along(cID)
      names(diagID) <- colnames(B)
    }

    #winRR <- sum(B==1, na.rm = TRUE)
    df <- plyr::ldply(diagID, function(i){
      data.frame(index = i, RR = sum(B[,i]==1, na.rm = TRUE)/ (NROW(B)-abs(as.numeric(colnames(B)[i]))))
    }, .id = "Diagonal")

    df$group  <- 1
    df$labels <- paste(df$Diagonal)
    # df$labels[df$Diagonal==0] <- ifelse(AUTO,"LOI","LOS")
    df$labels <- factor(df$labels,levels = df$labels, ordered = TRUE)
    return(df)
  })


  if(doShuffle){
    names(out) <- c("obs",paste0("Shuffled", 1:Nshuffle))
  } else {
    names(out) <- "obs"
  }

  # for(r in seq_along(out)){
  #
  #   if(r==1){
  #     RMd <- RM
  #   } else {
  #     RMd <- rp(y1 = y1, y2 = TSrnd[,(r-1)],emDim = emDim, emLag = emLag, emRad = emRad, to.sparse = TRUE)
  #   }
    #rp_lineDist(RM,d = diagWin, matrices = TRUE)
    # B <- rp_nzdiags(RMd, removeNZ = FALSE, d = diagWin)
    # rm(RMd)
    #
    # diagID <- 1:NCOL(B)
    # names(diagID) <- colnames(B)
    # if(length(diagWin)<NCOL(B)){
    #   cID <- which(colnames(B)%in%diagWin)
    #   B <- B[,cID]
    #   diagID <- seq_along(cID)
    #   names(diagID) <- colnames(B)
    # }
    #
    # #winRR <- sum(B==1, na.rm = TRUE)
    # df <- plyr::ldply(diagID, function(i){
    #   data.frame(index = i, RR = sum(B[,i]==1, na.rm = TRUE)/ (NROW(B)-abs(as.numeric(colnames(B)[i]))))
    # }, .id = "Diagonal")
    #
    # df$group  <- 1
    # df$labels <- paste(df$Diagonal)
    # # df$labels[df$Diagonal==0] <- ifelse(AUTO,"LOI","LOS")
    # df$labels <- factor(df$labels,levels = df$labels, ordered = TRUE)
    #
    # out[[r]] <- df
    # #rm(df,B,cID,diagID)
  #}

  dy      <- plyr::ldply(out)
  if(doShuffle){
    df_shuf <- dplyr::filter(dy, .data$.id!="obs")
  } else {
    df_shuf <- dy
  }

  dy_m <- df_shuf %>%
    dplyr::group_by(.data$Diagonal,.data$labels) %>%
    dplyr::summarise(`Mean Shuffled` = mean(.data$RR), sdRRrnd = stats::sd(.data$RR)) %>%
    dplyr::ungroup()

  if(Nshuffle==1){
    dy_m$sdRRrnd <- 0
  }

  dy_m$ciHI <- dy_m$`Mean Shuffled` + 1.96*(dy_m$sdRRrnd/sqrt(Nshuffle))
  dy_m$ciLO <- dy_m$`Mean Shuffled` - 1.96*(dy_m$sdRRrnd/sqrt(Nshuffle))

  obs <- dy[dy$.id=="obs",]
  if(doShuffle){
    if(!all(obs$Diagonal%in%dy_m$Diagonal)){
      notInd <- which(!obs$Diagonal%in%dy_m$Diagonal)
      tmp <- matrix(nrow = NROW(dy_m), ncol = NCOL(dy_m), dimnames = list(NULL,colnames(dy_m)))
      for(r in 1:NROW(dy_m)){
        if(!r%in%notInd){
          tmp[r,] <- obs[r,]
        }
      }
      obs <- tmp
      rm(tmp)
    }
  }
  dy_m$Observed <- obs$RR

  df <- tidyr::gather(dy_m, key = "variable", value = "RR", -c(.data$Diagonal,.data$sdRRrnd, .data$labels, .data$ciLO, .data$ciHI))
  df$Diagonal <- as.numeric(df$Diagonal)

  if(doPlot){

    # Diags <- as.numeric(levels(df$labels))
    # Diags[is.na(Diags)] <- 0

    # if(length(diagWin)>21){
    #   ext <- max(min(abs(Diags),na.rm = TRUE),abs(max(Diags,na.rm = TRUE)))-1
    #   breaks <- which(Diags%in%(c(seq(-ext,-1,length.out = 10),0,seq(1,ext,length.out = 10))))
    #   labels <- sort(unique(c(Diags[breaks],0)))
    #   breaks <- sort(unique(c(breaks,stats::median(breaks))))
    # } else {
    #   breaks <- seq_along(Diags)
    #   labels <- sort(unique(c(Diags[breaks],0)))
    # }


    # if(which.min(c(length(breaks),length(labels)))==1){
    #   labels <- labels[1:length(breaks)]
    # } else {
    #   breaks <- breaks[1:length(labels)]
    # }

    breaks <- seq_along(as.numeric(levels(df$labels)))
    labels <- sort(unique(c(as.numeric(levels(df$labels))[breaks],0)))

    if(length(breaks)>21){
      ID <- round(seq(1,length(breaks),length.out = 21))
      if(any(labels[ID]==0)){
        labels <- labels[ID]
        breaks <- ID
      } else {
        breaks <- sort(c(which(labels==0),ID))
        labels <- labels[breaks]
      }
    }

    x1 <- round((which.min(as.numeric(paste(df$Diagonal))))+(length(diagWin)*.1))
    x2 <- round((which.max(as.numeric(paste(df$Diagonal))))-(length(diagWin)*.1))
    yL <- max(as.numeric(paste(df$RR)),na.rm = TRUE)+0.1
    # col <- c("ciHI" = "grey70", "ciLO" = "grey70", "meanRRrnd" = "grey40","y_obs" = "black")
    if(Nshuffle>0){
      col <- c("Mean Shuffled" = "grey40","Observed" = "black")
      leg <- paste0("Diagonal Profile\n(N surrogates = ",Nshuffle,")")
    } else {
      col <- c("Observed" = "black")
      leg <- paste0("Diagonal Profile")
    }

    #siz <- c("ciHI" = .5, "ciLO" = .5, "meanRRrnd" = .5,"y_obs" = 1)
    #g<- ggplot(df, aes_(x = ~labels, y = ~RR, colour = ~variable))
    #   geom_line() + theme_bw()
    maxData <- max(df$RR, na.rm = TRUE)
    g <- ggplot2::ggplot(df, ggplot2::aes_(x=~Diagonal)) +
      ggplot2::geom_vline(xintercept = df$Diagonal[df$labels=="0"][1], size=1, colour = "grey80")
      if(doShuffle){
        maxData <- max(c(maxData,df$ciHI), na.rm = TRUE)
        g <- g + ggplot2::geom_ribbon(ggplot2::aes_(ymin=~ciLO, ymax=~ciHI), alpha=0.3)
      }
    g <- g + ggplot2::geom_line(ggplot2::aes_(y=~RR, colour = ~variable), size = .5) +
      #ggplot2::geom_line(ggplot2::aes_(y=~y_obs), colour = "black", size = 1) +
      ggplot2::scale_y_continuous("Recurrence Rate",limits = c(0,max(c(.5,maxData)))) +
      ggplot2::scale_x_continuous(name = paste0(xname, " << recurrences due to >> ", yname),
                                  breaks = breaks, labels = labels, expand = c(0,0)) +
      # ggplot2::geom_label(x=x1,y=.8,label=paste0("Recurrences due to\n ",xname),hjust="left", inherit.aes = FALSE, size = 3) +
      # ggplot2::geom_label(x=x2,y=.8,label=paste0("Recurrences due to\n ",yname),hjust="right", inherit.aes = FALSE, size = 3) +
      ggplot2::scale_colour_manual(leg, values = col) +
      ggplot2::theme_bw() +
      theme(panel.grid.minor = element_blank(),
            axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5, size = rel(.6)),
            axis.title.x = element_text(size = rel(.8)))

    if(!returnOnlyPlot){
      print(g)
    }
    if(doShuffle){
      df <- tidyr::spread(df,key = .data$variable, value = .data$RR)
    }

    return(invisible(list(plot = g, data = df)))

  } else {

    if(doShuffle){df <- tidyr::spread(df,key = .data$variable, value = .data$RR)}
    return(df)

  }
}


#' Fast (C)RQA (command line crp)
#'
#' This function will run the [commandline Recurrence Plots](http://tocsy.pik-potsdam.de/commandline-rp.php) executable provided by Norbert Marwan.
#'
#' @aliases crqa_cl
#' @param y1 Time series 1
#' @param y2 Time series 2 for Cross Recurrence Analysis (default = `NULL`)
#' @param emDim Embedding dimensions (default = `1`)
#' @param emLag Embedding lag (default = `1`)
#' @param emRad Radius on distance matrix (default = `1`)
#' @param DLmin Minimum length of diagonal structure to be considered a line (default = `2`)
#' @param VLmin Minimum length of vertical structure to be considered a line (default = `2`)
#' @param theiler Theiler window (default = `0`)
#' @param win Window to calculate the (C)RQA (default = minimum of length of `y1` or `y2`)
#' @param step Stepsize for sliding windows (default = size of `win`, so no sliding window)
#' @param JRP Wether to calculate a Joint Recurrence Plot (default = `FALSE`)
#' @param distNorm One of "EUCLIDEAN" (default), `"MAX", "MIN"`, or `"OP"` for an Order Pattern recurrence matrix
#' @param standardise Standardise data: `"none"` (default), `"mean.sd"`, or `"median.mad"`
#' @param returnMeasures Return the (C)RQA measures? (default = `TRUE`)
#' @param returnRPvector Return the recurrent points in a dataframe? (default = `FALSE`)
#' @param returnLineDist Return the distribution of diagonal and horizontal line length distances (default = `FALSE`)
#' @param doPlot Produce a plot of the recurrence matrix by calling [rp_plot()], values can be `"rp"` (the thresholded recurrence matrix),`"distmat"` (the unthresholded recurrence matrix) or `"noplot"` (default = `"noplot"`)
#' @param path_to_rp Path to the command line executable (default = path set during installation, use `getOption("casnet.path_to_rp")` to see)
#' @param saveOut Save the output to files? If `TRUE` and `path_out = NA`, the current working directory will be used (default = `FALSE`)
#' @param path_out Path to save output if `saveOut = TRUE` (default = `NULL`)
#' @param file_ID A file ID which will be a prefix to to the filename if `saveOut = TRUE` (default = `NULL`, an integer will be added tot the file name to ensure unique files)
#' @param silent Do not display any messages (default = `TRUE`)
#' @param surrogateTest Perform surrogate tests. If `TRUE`, will run surrogate tests using default settings for a two-sided test of \eqn{H_0: The data generating process is a rescaled linear Gaussian process} at \eqn{\alpha = .05} (arguments `ns = 39, fft = TRUE, amplitude = TRUE`)
#' @param targetValue A value passed to `est_radius(...,type="fixed", targetMeasure="RR")` if `is.na(emRad)==TRUE`. This is useful for windowed analysis, it will estimate a new radius for each window.
#' @param useParallel Speed up calculations by using the parallel processing options provided by `parallel` to assign a seperate process/core for each window in windowed (C)RQA analysis and [parallel::detectCores()] with`logical = TRUE` to decide on the available cores (default = `FALSE`)
#' @param ... Additional parameters (currently not used)
#'
#' @details The `rp` executable is installed when the function is called for the first time and is renamed to `rp`, from a platform specific filename downloaded from http://tocsy.pik-potsdam.de/commandline-rp.php or extracted from an archive located in the directory:
#' `...\\casnet\\commandline_rp\\`
#' The file is copied to the directory: `...\\casnet\\exec\\`
#' The latter location is stored as an option and can be read by calling `getOption("casnet.path_to_rp")`.
#'
#' @section Troubleshooting:
#' Some notes on resolving errors with `rp`.The script will first try to download the correct executable, if that fails, the copy will have... failed. It should be relatively easy to get `rp_cl()` working though, by using some custom settings:
#'
#'
#' + *Copy failed* - Every time the function `rp_cl()` is called it will check whether a log file `rp_instal_log.txt` is present in the `...\\casnet\\exec\\` directory. If you delete the `rp_instal_log.txt` file, and call the function, another attempt will be made to download a copy of the executable.
#'
#' + *Copy still fails and/or no permission to copy* - If you cannot acces the directory `...\\casnet\\commandline_rp\\`, download the appropriate executable from the [Commandline Recurrence Plots](http://tocsy.pik-potsdam.de/commandline-rp.php) page and copy to a directory you do have the rights to: *execute* commands, *write* and *read* files. Make sure you rename the file to `rp` (`rp.exe` on Windows OS). Then, either pass the path to `rp` as the argument `path_to_rp` in the `rp_cl(.., path_to_rp = "YOUR_PATH")` function call, or, as a more permanent solution, set the `path_to_rp` option by calling `options(casnet.path_to_rp="YOUR_PATH")`.
#'
#' + *Error in execution of `rp`* - This can have a variety of causes, the `rp` executable is called using [system2()] and makes use of the [normalizePath()] function with argument `mustWork = FALSE`. Problems caused by specific OS, machine, or, locale problems (e.g. the `winslash` can be reported as an \href{https://github.com/FredHasselman/casnet/issues}{issue on Github}). One execution error that occurs when the OS is not recognised properly can be resolved by chekcing `getOption("casnet.rp_prefix")`. On Windows OS this should return an empty character vector, on Linux or macOS it should return `"./"`. You can manually set the correct prefix by calling `options(casnet.rp_prefix="CORRECT OS PREFIX")` and fill in the prefix that is correct for your OS
#'
#'
#' @return A list object containing 1-3 elements, depending on arguments requesting output.
#'
#' * `rqa_measures` - A list of the (C)RQA measures returned if `returnMeasures = TRUE`:
#'      1. RR = Recurrence rate
#'      2. DET = Determinism
#'      3. DET_RR = Ratio DET/RR
#'      4. LAM = Laminarity
#'      5. LAM_DET = Ratio LAM/DET
#'      6. L_max = maximal diagonal line length
#'      7. L_mean = mean diagonal line length
#'      8. L_entr = Entropy of diagonal line length distribution
#'      9. DIV =  Divergence (1/L_max)
#'      10. V_max = maximal vertical line length
#'      11. TT = Trapping time
#'      12. V_entr = Entropy of vertical line length distribution
#'      13. T1 = Recurrence times 1st type
#'      14. T2 = Recurrence times 2nd type
#'      15. W_max = Max interval length
#'      16. W_mean = Mean of interval lengths
#'      17. W_entr = Entropy of interval length distribution
#'      18. W_prob = Probability of interval
#'      19. F_min = F min
#' * `rqa_rpvector` - The radius thresholded distance matrix (recurrence matrix), which can be visualised as a recurrence plot by calling [rp_plot()]. If a sliding window analysis is conducted this will be a list of matrices and could potentially grow too large to handle. It is recommended you save the output to disk by setting `saveOut = TRUE`.
#' * `rqa_diagdist` - The distribution of diagonal line lengths
#'
#'
#' @note The platform specific `rp` command line executables were created by Norbert Marwan and obtained under a Creative Commons License from the website of the Potsdam Institute for Climate Impact Research at [http://tocsy.pik-potsdam.de/](http://tocsy.pik-potsdam.de/).
#'
#' The full copyright statement on the website is as follows:
#'
#' (C) 2004-2017 SOME RIGHTS RESERVED
#'
#' University of Potsdam, Interdisciplinary Center for Dynamics of Complex Systems, Germany
#'
#' Potsdam Institute for Climate Impact Research, Transdisciplinary Concepts and Methods, Germany
#'
#' This work is licensed under a [Creative Commons Attribution-NonCommercial-NoDerivs 2.0 Germany License](https://creativecommons.org/licenses/by-nc-nd/2.0/de/).
#'
#' More information about recurrence analysis can be found on the [Recurrence Plot](http://www.recurrence-plot.tk) website.
#'
#' @family Recurrence Quantification Analysis

#' @export
#'
rp_cl            <- function(y1,
                             y2      = NULL,
                             emDim   = 1,
                             emLag   = 1,
                             emRad   = NA,
                             DLmin   = 2,
                             VLmin   = 2,
                             theiler = 0,
                             win            = min(length(y1),ifelse(is.null(y2),(length(y1)+1), length(y2)), na.rm = TRUE),
                             step           = win,
                             JRP            = FALSE,
                             distNorm       = c("EUCLIDEAN", "MAX", "MIN", "OP")[[1]],
                             standardise    = c("none","mean.sd","median.mad")[1],
                             returnMeasures = TRUE,
                             returnRPvector = FALSE,
                             returnLineDist = FALSE,
                             doPlot         = c("noplot","rp","distmat")[[1]],
                             path_to_rp     = getOption("casnet.path_to_rp"),
                             saveOut        = FALSE,
                             path_out       = NULL,
                             file_ID        = NULL,
                             silent         = TRUE,
                             surrogateTest  = FALSE,
                             targetValue    = .05,
                             useParallel    = FALSE,
                             ...){

  if(!file.exists(normalizePath(file.path(getOption("casnet.path_to_rp"),"rp_install_log.txt"), mustWork = FALSE))){
    set_command_line_rp()
  }

  os <- get_os()
  sysinf <- Sys.info()
  if(os%in%"osx"&as.numeric(strsplit(sysinf[which(names(sysinf)%in%"release")],"[.]")[[1]][1])>=19){
    stop("As of macOS Catalina, 32-bit applications are no longer supported, please use rp_measures() to conduct Recurrence Quantification Analysis.")
  }

  cat("\n~~~o~~o~~casnet~~o~~o~~~\n")

  dotArgs <- list(...)
  if(is.na(dotArgs%00%NA)){dotArgs<-list(nothing=NA)}

  if(!is.null(y2)){
    y2 <- zoo::as.zoo(y2)
    N1 <- NROW(y1)
    N2 <- NROW(y2)
    if(N1>N2){
      y2[N2:(N1+(N1-N2))] <- 0
    }
    if(N1<N2){
      y1[N1:(N2+(N2-N1))] <- 0
    }
    df <- cbind.data.frame(y1=unclass(y1),y2=unclass(y2))
    df <- df[stats::complete.cases(df),]
  } else {
    if(any(is.na(y1))){
      y1 <- y1[!is.na(y1)]
      if(!silent){message("Removed NAs from timeseries y1 before (C)RQA")}
    }
    df <- cbind.data.frame(y1=unclass(y1),y2=NA)
  }

  # Begin input checks
  if(is.null(path_out)){
    if(saveOut){
      path_out <- getwd()
    } else {
      path_out <- tempdir(check = TRUE)
    }
  } else {
    path_out <- normalizePath(path_out, mustWork = TRUE)
  }

  if(is.null(y2)){
    cat("\nPerforming auto-RQA\n")
    if(theiler<1){theiler=1}
  } else {
    cat("\nPerforming cross-RQA\n")
  }

  if(surrogateTest){
    surrogateSeries <- tseries::surrogate(y1,ns=39, fft = TRUE, amplitude = TRUE)
  }

  windowedAnalysis <- FALSE
  if(win==NROW(df)|step==(NROW(df))){
    win <- step <- NROW(df)
    wIndex <- seq(1,NROW(df))
    if(is.na(emRad)){cat(paste("Radius will be calulated to yield targetValue =",targetValue,"RR...\n"))}
  } else {
    windowedAnalysis = TRUE
    wIndex <- seq(1,NROW(df)-win, by = step)
    cat(paste("\nCalculating recurrence measures in a window of size",win,"taking steps of",step,"data points. \n"))
    if(is.na(emRad)){cat(paste("Radius will be re-calculated to yield targetValue =",targetValue,"RR for each window...\n\n"))}
  }


  if(windowedAnalysis){

    # Check window length vs. embedding
    if(win < (emLag*emDim)){

      stop(paste0("The size of win = ", win, " must be larger than the product of emLag = ",emLag, " and emDim = ", emDim))

    } else {

      # Adjust time series lengths
      if((dplyr::last(wIndex)+win-NROW(df))>0){

        yy1 <- ts_trimfill(x=seq(1,dplyr::last(wIndex)+win),y=df[,1])
        yy2 <- rep(NA,length(yy1))
        if(!all(is.na(df[,2]))){
          yy2 <- ts_trimfill(x=seq(1,dplyr::last(wIndex)+win),y=df[,2])
        }
        df <- cbind.data.frame(y1=yy1,y2=yy2)
        if(!silent){message("\nAdded 0s to vector to make window and stepsize combination fit to time series length\n\n ")}
      } else {
        if(!silent){message("\nWindow and stepsize combination do not fit to the full length of the time series!\n\n ")}
      }

      if(silent){cat("\nsst! [it's a silent-ish mode]\n\n")}

      wIndices <- plyr::llply(wIndex, function(w){seq(w,w+(win))})
      names(wIndices) <- paste0("window: ",seq_along(wIndices)," | start: ",wIndex," | stop: ",wIndex+win)

    }
  } else {

    wIndices <- list(wIndex)
    names(wIndices) <- paste0("window: 1 | start: ",wIndex[1]," | stop: ",wIndex[NROW(df)])

  } # If windowed

  #cl <- parallel::makeCluster(mc.cores)
  #parallel::clusterExport(cl = cl, c("rp_cl_main","wIndices","df","y1","y2","emDim","emLag","emRad","DLmin","VLmin","theiler", "win","step","JRP","distNorm","returnMeasures","returnRPvector","returnLineDist","doPlot","path_to_rp","saveOut","path_out","file_ID","silent","..."))

  # Create the data, for all windows,need this for parallel, but also speeds up processing in general.
  dfList <- plyr::llply(wIndices, function(ind){cbind(y1 = df[ind,1],y2 = df[ind,2])})

  if(useParallel){

    if(names(dotArgs)%in%"logical"){logical <- logical} else {logical <- TRUE}

    cat("\n...using parallel processing...\n")

    # Only return measures
    returnRPvector <- FALSE
    returnLineDist <- FALSE

    # How many cores is optimal?
    cores_available <- parallel::detectCores(logical = logical)-1
    if(cores_available>1){
      if(length(wIndices)%%2==0){odd_windows <- FALSE} else {odd_windows <- TRUE}
      if(cores_available%%2==0){odd_cores <- FALSE} else {odd_cores <- TRUE}
      if(odd_windows&!odd_cores){cores_available<-cores_available-1}
      if(!odd_windows&odd_cores){cores_available<-cores_available-1}
    } else {
      cores_available <- 1
    }

    cl       <- parallel::makeCluster(cores_available)

    # parallel::clusterEvalQ(cl, library(callr))
    parallel::clusterEvalQ(cl,library(utils))
    parallel::clusterEvalQ(cl,library(plyr))
    parallel::clusterEvalQ(cl,library(tidyverse))
    parallel::clusterEvalQ(cl,library(pROC))
    parallel::clusterEvalQ(cl,library(Matrix))
    parallel::clusterEvalQ(cl,library(casnet))

    # parallel::clusterExport(cl, varlist = c("data","emDim","emLag","emRad","DLmin","VLmin","theiler","win","step","JRP","distNorm","returnMeasures","returnRPvector","returnLineDist","doPlot","path_to_rp", "saveOut","path_out","file_ID","silent","targetValue", "useParallel"))

    # cluster_library(c("devtools","utils","plyr","dplyr","tidyr","Matrix","pROC")) %>%
    # cluster_assign_value("rp_cl_main", rp_cl_main) %>%
    # cluster_assign_value("est_radius", est_radius) %>%
    # cluster_assign_value("emDim", emDim)  %>%
    # cluster_assign_value("emLag", emLag)  %>%
    # cluster_assign_value("emRad", emRad)  %>%
    # cluster_assign_value("DLmin", DLmin)  %>%
    # cluster_assign_value("VLmin", VLmin)  %>%
    # cluster_assign_value("theiler", theiler)  %>%
    # cluster_assign_value("JRP", JRP)  %>%
    # cluster_assign_value("distNorm", distNorm)  %>%
    # cluster_assign_value("returnMeasures", returnMeasures)  %>%
    # cluster_assign_value("returnRPvector", returnRPvector)  %>%
    # cluster_assign_value("returnLineDist", returnLineDist)  %>%
    # cluster_assign_value("doPlot", doPlot)  %>%
    # cluster_assign_value("path_to_rp", path_to_rp)  %>%
    # cluster_assign_value("saveOut", saveOut)  %>%
    # cluster_assign_value("path_out", path_out)  %>%
    # cluster_assign_value("file_ID", file_ID)  %>%
    # cluster_assign_value("silent", silent)  %>%
    # cluster_assign_value("targetValue", targetValue) %>%
    # cluster_assign_value("useParallel", useParallel)

    #  parallel::clusterExport(cl = cl, c("rp_cl_main","wIndices","df","y1","y2","emDim","emLag","emRad","DLmin","VLmin","theiler", "win","step","JRP","distNorm","returnMeasures","returnRPvector","returnLineDist","doPlot","path_to_rp","saveOut","path_out","file_ID","silent","..."))


    start <- proc.time()
    wList <- parallel::parLapply(cl,dfList,function(df){rp_cl_main(data = df,
                                                                   emDim          = emDim,
                                                                   emLag          = emLag,
                                                                   emRad          = emRad,
                                                                   DLmin          = DLmin,
                                                                   VLmin          = VLmin,
                                                                   theiler        = theiler,
                                                                   JRP            = JRP,
                                                                   distNorm       = distNorm,
                                                                   returnMeasures = returnMeasures,
                                                                   returnRPvector = returnRPvector,
                                                                   returnLineDist = returnLineDist,
                                                                   doPlot    = doPlot,
                                                                   path_to_rp     = path_to_rp,
                                                                   saveOut        = saveOut,
                                                                   path_out       = path_out,
                                                                   file_ID        = file_ID,
                                                                   silent         = silent,
                                                                   targetValue    = targetValue,
                                                                   useParallel    = useParallel)})

    parallel::stopCluster(cl)
    time_elapsed_parallel <- proc.time() - start # End clock

    cat("\nCompleted in:\n")
    print(time_elapsed_parallel)

  } else {

    cat("\n...using sequential processing...\n")

    start <- proc.time()

    wList <- plyr::llply(dfList, function(data){
      rp_cl_main(
        data           = data,
        emDim          = emDim,
        emLag          = emLag,
        emRad          = emRad,
        DLmin          = DLmin,
        VLmin          = VLmin,
        theiler        = theiler,
        JRP            = JRP,
        distNorm       = distNorm,
        returnMeasures = returnMeasures,
        returnRPvector = returnRPvector,
        returnLineDist = returnLineDist,
        doPlot    = doPlot,
        path_to_rp     = path_to_rp,
        saveOut        = saveOut,
        path_out       = path_out,
        file_ID        = file_ID,
        silent         = silent,
        targetValue    = targetValue)},
      .progress = plyr::progress_text(char = "o~o"))

    time_elapsed_parallel <- proc.time() - start # End clock

    cat(paste("\nCompleted in:\n"))
    print(time_elapsed_parallel)


  } # useParallel
  rm(dfList)

  rqa_measures <- rqa_rpvector <- rqa_diagdist <- rqa_horidist <- NA
  if(useParallel){
    if(is.list(wList)){
      rqa_measures <-  plyr::ldply(wList)
    } else {
      rqa_measures <-  tidyr::unnest(wList)
    }
  } else {
    rqa_measures <-plyr::ldply(wList, function(l) l$measures) # %>% dplyr::mutate(win = win, step = step, index = attr(wlist, "index"))
    rqa_rpvector <-plyr::ldply(wList, function(l) l$rpMAT) # %>% dplyr::mutate(win = win, step = step, index = attr(wlist, "index"))
    rqa_diagdist <-plyr::ldply(wList, function(l) l$diag_disthist) # %>% dplyr::mutate(win = win, step = step, index = attr(wlist, "index"))
    rqa_horidist <-plyr::ldply(wList, function(l) l$hori_disthist) # %>% dplyr::mutate(win = win, step = step, index = attr(wlist, "index"))
  }

  wPlot <- which(doPlot%in%c("noplot","rp","distmat"))

  if(wPlot>1){
    if(wPlot==2){
      if(windowedAnalysis){
        plotList <- plyr::llply(wIndices, function(ind) rp(y1 = df[ind,1], y2 = df[ind,2], emDim = emDim, emLag = emLag, emRad = emRad, doPlot = TRUE))
        cowplot::plot_grid(plotList)
      } else {
        rp(y1 = df[,1], y2 = df[,2], emDim = emDim, emLag = emLag, emRad = emRad, doPlot = TRUE)
      }
    }
    if(wPlot==3){
      if(windowedAnalysis){
        plotList <- plyr::llply(wIndices, function(ind) rp(y1 = df[ind,1],y2 = df[ind,2], emDim = emDim, emLag = emLag, doPlot = TRUE))
        cowplot::plot_grid(plotList)
      } else {
        rp(y1 = df[,1],y2 = df[,2], emDim = emDim, emLag = emLag, doPlot = TRUE)
      }
    }

  }

  justData <- FALSE
  if(win!=NROW(df)){
    if(!silent){message("Windowed (c)rqa: Returning only (c)rqa measures, use saveOut = TRUE, to get the distribution files")}
    justData <- TRUE
  }

  if(returnMeasures&!returnRPvector&!returnLineDist){
    justData<- TRUE
  }

  out <- list(rqa_measures = rqa_measures,
              rqa_rpvector = rqa_rpvector,
              rqa_diagdist = rqa_diagdist,
              rqa_horidist = rqa_horidist)
  if(saveOut){saveRDS(out,paste0(path_out,"CRQA_out",file_ID,".rds"))}

  cat("\n~~~o~~o~~casnet~~o~~o~~~\n")

  if(justData){
    return(rqa_measures)
  } else {
    return(out[c(returnMeasures,returnRPvector,returnLineDist,returnLineDist)])
  }
}


#' rp_cl_main
#'
#' @inheritParams rp_cl
#'
#' @keywords internal
#'
#' @export
#'
rp_cl_main <- function(data,
                       emDim  = 1,
                       emLag  = 1,
                       emRad  = NA,
                       DLmin = 2,
                       VLmin = 2,
                       theiler = 0,
                       win     = min(length(y1),ifelse(is.null(y2),(length(y1)+1), length(y2)), na.rm = TRUE),
                       step    = step,
                       JRP     = FALSE,
                       distNorm       = c("EUCLIDEAN", "MAX", "MIN", "OP")[[1]],
                       returnMeasures = TRUE,
                       returnRPvector = FALSE,
                       returnLineDist = FALSE,
                       doPlot = c("noplot","rp","distmat")[[1]],
                       path_to_rp = getOption("casnet.path_to_rp"),
                       saveOut    = FALSE,
                       path_out   = NULL,
                       file_ID    = NULL,
                       silent     = TRUE,
                       targetValue  = .05,
                       useParallel = FALSE,
                       ...){


  if(useParallel){
    #data <- unnest(data)
    returnRPvector = FALSE
    returnLineDist = FALSE
  }
  y1     <- data[,1]
  y2     <- data[,2]

  fixedRR <- FALSE
  if(is.na(emRad)){fixedRR=TRUE}

  file_ID <- file_ID%00%0

  RQAmeasures <- list(
    RR      = 'Recurrence rate',
    DET     = 'Determinism',
    DET_RR  = 'Ratio DET/RR',
    LAM     = 'Laminarity',
    LAM_DET = 'Ratio LAM/DET',
    L_max   = 'maximal diagonal line length',
    L_mean  = 'mean diagonal line length',
    L_entr  = 'Entropy of diagonal line length distribution',
    DIV     = 'Divergence (1/L_max)',
    V_max   = 'maximal vertical line length',
    TT      = 'Trapping time',
    V_entr  = 'Entropy of vertical line length distribution',
    T1      = 'Recurrence times 1st type',
    T2      = 'Recurrence times 2nd type',
    W_max   = 'Max interval',
    W_mean  = 'Mean of interval',
    W_entr  = 'Entropy interval distribution',
    W_prob  = 'Prob.of intervals',
    F_min   = 'F min'
  )

  if(doPlot%in%"distmat"){-1 * emRad}

  tmpd  <- tempdir()
  tmpf1 <- tempfile(tmpdir = tmpd, fileext = ".dat")
  utils::write.table(as.data.frame(y1), tmpf1, col.names = FALSE, row.names = FALSE, sep = "\t")

  #  fileSep <- ifelse(any(path_out%in%"/"),"/","\\")
  file_ID=1
  plotOUT     <- file.path(normalizePath(path_out,mustWork = TRUE),paste0("RQAplot_",     file_ID, ".txt"))
  measOUT     <- normalizePath(file.path(path_out,paste0("RQAmeasures_", file_ID, ".txt")), mustWork = FALSE)
  histOUTdiag <- normalizePath(file.path(path_out,paste0("RQAhist_diag_",file_ID, ".txt")), mustWork = FALSE)
  histOUThori <- normalizePath(file.path(path_out,paste0("RQAhist_hori_",file_ID, ".txt")), mustWork = FALSE)

  if(any(is.null(y2))|any(is.na(y2%00%NA))){

    if(is.na(emRad)){
      if(!is.na(targetValue)){
        emRad <- est_radius(y1 = y1,
                            emDim = emDim,
                            emLag = emLag,
                            targetValue = targetValue,
                            radiusOnFail = "percentile", tol = .2, silent = silent)
        if(emRad$Converged){
          emRad <- emRad$Radius
        } else {
          emRad <-  emRad <- stats::sd(c(y1,y2),na.rm = T)
        }
      }
    }
    opts <- paste("-i", shQuote(tmpf1),
                  "-r", shQuote(plotOUT),
                  "-o", shQuote(measOUT),
                  "-p", shQuote(histOUTdiag),
                  "-q", shQuote(histOUThori),
                  "-m", emDim,
                  "-t", emLag,
                  "-e", emRad,
                  "-l", DLmin,
                  "-v", VLmin,
                  "-w", theiler,
                  "-n", shQuote(distNorm),
                  ifelse(silent,"-s",""))
  } else {
    if(is.na(emRad)){
      if(!is.na(targetValue)){
        emRad <- est_radius(y1 = y1, y2 = y2, emDim = emDim, emLag = emLag, targetValue = targetValue, tol = .2, radiusOnFail = "percentile", silent = silent)
        if(emRad$Converged){
          emRad <- emRad$Radius
        } else {
          emRad <- stats::sd(c(y1,y2),na.rm = T)
        }
      }
    }
    tmpf2 <- tempfile(tmpdir = tmpd, fileext = ".dat")
    utils::write.table(as.data.frame(y2), tmpf2, col.names = FALSE, row.names = FALSE, sep = "\t")
    opts <- paste("-i", shQuote(tmpf1),
                  "-j", shQuote(tmpf2),
                  "-r", shQuote(plotOUT),
                  "-o", shQuote(measOUT),
                  "-p", shQuote(histOUTdiag),
                  "-q", shQuote(histOUThori),
                  "-m", emDim,
                  "-t", emLag,
                  "-e", emRad,
                  "-l", DLmin,
                  "-v", VLmin,
                  "-w", theiler,
                  "-n", shQuote(distNorm),
                  ifelse(silent,"-s",""))
  }

  #closeAllConnections()

  ## RCMD
  #callr::rcmd_copycat(cmd = paste0(getOption("casnet.rp_command")), cmdargs = opts, wd = file.path(normalizePath(path_to_rp, mustWork = FALSE)), show = silent, echo = TRUE)

  system2(command = file.path(path_to_rp,getOption("casnet.rp_command")), args = opts)

  measures     <- return_error(utils::read.delim(normalizePath(gsub("[']+","",measOUT)),header=TRUE))
  rpMAT        <- return_error(utils::read.delim(normalizePath(gsub("[']+","",plotOUT)),header=TRUE))
  disthistDiag <- return_error(utils::read.delim(normalizePath(gsub("[']+","",histOUTdiag)), header=FALSE, sep = " "))
  disthistHori <- return_error(utils::read.delim(normalizePath(gsub("[']+","",histOUThori)), header=FALSE, sep = " "))

  if(all(is.null(measures$warning),is.data.frame(measures$value))){
    measures <- measures$value
  } else {
    measures <- rbind.data.frame(rep(NA,length(RQAmeasures)))
  }
  colnames(measures) <- gsub("(#RR)|(X.RR)","RR",colnames(measures))
  measures <- cbind.data.frame(measures, emDim=emDim, emLag=emLag, emRad=emRad, DLmin=DLmin, VLmin=VLmin, distNorm =distNorm)

  if(all(is.null(rpMAT$warning),is.data.frame(rpMAT$value))){
    rpMAT <- rpMAT$value
  } else {
    rpMAT <- data.frame(y1=NA,y2=NA,dist=NA)
  }

  if(any(is.null(disthistDiag$warning),is.data.frame(grepl("Error",paste(disthistDiag$value))))){
    disthistDiag <- disthistDiag$value
  } else {
    disthistDiag <- data.frame(line.length=NA,freq=NA)
  }

  if(any(is.null(disthistHori$warning),is.data.frame(grepl("Error",paste(disthistHori$value))))){
    disthistHori <- disthistHori$value
  } else {
    disthistHori <- data.frame(line.length=NA,freq=NA)
  }

  if(!silent){cat(paste0("[ID ",file_ID,"] Analysis completed... "))}
  if(!saveOut){
    file.remove(measOUT)
    file.remove(plotOUT)
    file.remove(histOUTdiag)
    file.remove(histOUThori)
    if(!silent){cat("temporary files removed ...\n")}
  } else {
    if(!silent){
      cat("files saved ...\n")
      cat(measOUT,"\n",plotOUT,"\n",histOUTdiag,"\n",histOUThori,"\n")
    }
  }

  if(useParallel){
    return(tibble::as_tibble(measures))
  } else{
    return(list(measures = measures,
                rpMAT    = rpMAT,
                diag_disthist = disthistDiag,
                hori_disthist = disthistHori))
  }
}



#' rp_nzdiags
#'
#' @description Get all nonzero diagonals of a binary matrix, or, diagonals specified as a vector by argument `d`.
#'
#' @param RM A binary (0,1) matrix.
#' @param d An optional vector of diagonals to extract.
#' @param returnVectorList Return list
#' @param returnNZtriplets Return a dataframe with coordinates of only nonzero elements in diagonals (default = `FALSE`)
#' @param removeNZ Remove nonzero diagonals if `TRUE`. If `FALSE` returns the full diagonals matrix. Use e.g. to plot diagonal recurrence profiles (default = `TRUE`)
#' @param silent Silent-ish mode
#'
#' @author Fred Hasselman
#'
#' @return A matrix object with nonzero diagonals as columns and/or a dataframe with coordinates of nonzero diagonal elements
#'
#' @export
#'
#' @family Distance matrix operations (recurrence plot)
#'
rp_nzdiags <- function(RM=NULL, d=NULL, returnVectorList=TRUE, returnNZtriplets=FALSE, removeNZ = TRUE, silent = TRUE){
  # Loosely based on MATLAB function spdiags() by Rob Schreiber - Copyright 1984-2014 The MathWorks, Inc.

  # if(!is.na(win)){
  #   if(length(win)==1){
  #     win <- c(-win,win)
  #   }
  #   if(length(win)==2){
  #       win <- sort(win)
  #   }
  #   if(!length(win)%[]%c(1,2)){
  #     stop("Windowsize must be NA, or a 1, or 2 element vector!")
  #     }
  # }

  if(grepl("matrix",class(RM)[1],ignore.case = TRUE)){

    if(all(RM>0)){warning("All matrix elements are nonzero.")}

    s  <- Sys.time()

    nzdiagsM <- methods::as(RM, "TsparseMatrix") #methods::as(RM, "dgTMatrix")
    nzdiags  <- data.frame(row   = nzdiagsM@i,
                           col   = nzdiagsM@j,
                           value = nzdiagsM@x,
                           ndiag = (nzdiagsM@i)-(nzdiagsM@j)) #(nzdiagsM@j)-(nzdiagsM@i))
    nzdiags <- dplyr::arrange(nzdiags,nzdiags$ndiag)

    if(is.null(d)){
      d <-  c(-1*rev(1:(NCOL(nzdiagsM)-1)),0,1:(NCOL(nzdiagsM)-1))
    }
    rm(nzdiagsM)

    if(removeNZ){
      nd <- unique(nzdiags$ndiag)
      # Get diagonals which have nonzero elements
      d <- nd[nd%in%sort(as.vector(d))]
    }

    indMat  <- col(RM)-row(RM)
    indMatV <- as.vector(indMat)

    #cat("\nStart extraction of nonzero diagonals\n")
    m  <- NROW(RM)
    n  <- NCOL(RM)
    p  <- length(d)
    if(is.logical(RM)){
      B <- matrix(FALSE,nrow = max(c(m,n), na.rm = TRUE), ncol = p)
    } else {
      B <- matrix(0, nrow = max(c(m,n), na.rm = TRUE), ncol = p)
    }
    colnames(B) <- paste(d)

    for(i in seq_along(d)){
      B[(nzdiags$row[nzdiags$ndiag==d[i]]+1), i] <- nzdiags$value[nzdiags$ndiag==d[i]]
    }

    zID <- which(Matrix::colSums(Matrix::as.matrix(B))==0)
    if(length(zID)>0){
      if(removeNZ){
        B  <- B[,-zID]
        nzdiags <- nzdiags[nzdiags$ndiag%in%zID,]
      }
    }
    #dl <- plyr::llply(seq_along(toList),function(dd){RM[indMat==toList[[dd]]]}, .progress = plyr::progress_text(char = "o~o"))
    e  <- Sys.time()

    if(!silent){cat(paste0("\nTime: ", round(difftime(e,s, units = "mins"), digits=1)," min.\n"))}

    if(returnNZtriplets){
      return(list(B = B, nzdiags = nzdiags))
    } else {
      return(B)
    }
  } else {
    stop("Input to rp_nzdiags is not a matrix-like object.")
  }
}



rp_nzdiags_matlab <- function(RP,d=NULL){
  # Loosely based on MATLAB function spdiags() by Rob Schreiber - Copyright 1984-2014 The MathWorks, Inc.
  #require(Matrix)

  if(grepl("matrix",class(RP),ignore.case = TRUE)){

    A <- RP

    if(all(A>0)){warning("All matrix elements are nonzero.")}

    # create an indicator for all diagonals in the matrix
    ind   <- col(A)-row(A)

    # Split the matrix!
    spd <- split(A, ind)

    if(is.null(d)){

      # Get diagonals which have nonzero elements
      keepID <-plyr::ldply(spd, function(di) any(di>0))
      nzdiag <- spd[keepID$V1]
      # Indices of nonzero diagonals
      d      <- as.numeric(keepID$.id[keepID$V1])

    } else {

      # Diagonals are specified
      d <- sort(as.vector(d))
      keepID <- names(spd)%in%d
      nzdiag <- spd[keepID]
    }

    # Deal with uneven rows and cols
    m  <- NROW(A)
    n  <- NCOL(A)
    p  <- length(d)
    if(is.logical(A)){
      B <- matrix(FALSE,nrow = min(c(m,n), na.rm = TRUE), ncol = p)
    } else {
      B <- matrix(0, nrow = min(c(m,n), na.rm = TRUE), ncol = p)
    }

    for(i in seq_along(d)){
      B[zoo::index(nzdiag[[i]]),i] <- nzdiag[[i]]
    }
    colnames(B) <- paste(d)

    return(B)
  }
}


rp_nzdiags_chroma <- function(RP, d=NULL){
  # Loosely based on MATLAB function spdiags() by Rob Schreiber - Copyright 1984-2014 The MathWorks, Inc.
  #require(Matrix)

  if(grepl("matrix",class(RP),ignore.case = TRUE)){

    A <- RP

    if(all(A>0)){warning("All matrix elements are nonzero.")}

    # create an indicator for all diagonals in the matrix
    ind   <- col(A)-row(A)

    # Split the matrix!
    spd <- split(A, ind)

    if(is.null(d)){

      # Get diagonals which have nonzero elements
      keepID <-plyr::ldply(spd, function(di) any(di>0))
      nzdiag <- spd[keepID$V1]
      # Indices of nonzero diagonals
      d      <- as.numeric(keepID$.id[keepID$V1])

    } else {

      # Diagonals are specified
      d <- sort(as.vector(d))
      keepID <-  names(spd)%in%d
      nzdiag <- spd[keepID]
    }

    # Deal with uneven rows and cols
    m  <- NROW(A)
    n  <- NCOL(A)
    p  <- length(d)
    if(is.logical(A)){
      B <- matrix(FALSE,nrow = min(c(m,n), na.rm = TRUE), ncol = p)
    } else {
      B <- matrix(0, nrow = min(c(m,n), na.rm = TRUE), ncol = p)
    }

    for(i in seq_along(d)){
      B[zoo::index(nzdiag[[i]]),i] <- nzdiag[[i]]
    }
    colnames(B) <- paste(d)

    return(B)
  }
}


#' Line length distributions
#'
#' @inheritParams rp
#' @inheritParams rp_measures
#' @param d Vector of diagonals to be extracted from matrix `RP` before line length distributions are calculated. A one element vector will be interpreted as a windowsize, e.g., `d = 50` will extract the diagonal band `-50:50`. A two element vector will be interpreted as a band, e.g. `d = c(-50,100)` will extract diagonals `-50:100`. If `length(d) > 2`, the numbers will be interpreted to refer to individual diagonals, `d = c(-50,50,100)` will extract diagonals `-50,50,100`. If `length(d)` is `NULL`, 1 or 2, the theiler window is applied before diagonals are extracted. The theiler window is ignored if `length(d)>2`, or if it is larger than the matrix or band indicated by parameter `d`. A warning will be given is a theiler window was already applied to the matrix.
#' @param recurrenceTimes Relevant for Recurrence Time analysis: Return the distribution of 0 valued segments in nonzero diagonals/verticals/horizontals. This indicates the time between subsequent line structures.
#'
#' @description Extract lengths of diagonal, vertical and horizontal line segments from a recurrence matrix.
#'
#' @details Based on the Matlab function `linedists` by Stefan Schinkel, Copyright (C) 2009 Stefan Schinkel, University of Potsdam, http://www.agnld.uni-potsdam.de
#'
#' References:
#' S. Schinkel, N. Marwan, O. Dimigen & J. Kurths (2009):
#' "Confidence Bounds of recurrence-based complexity measures
#' Physics Letters A,  373(26), pp. 2245-2250
#'
#' Copyright (C) 2009 Stefan Schinkel, University of Potsdam
#' <http://www.agnld.uni-potsdam.de>
#'
#' @author Fred Hasselman
#' @return A list object with distributions of line lengths. If `matrices = TRUE` datafr are returned whose columns represent the nonzero diagonals, verticals, or, horizontals.
#'
#' @export
#'
#' @family Distance matrix operations (recurrence plot)
#'
rp_lineDist <- function(RM,
                        DLmin = 2,
                        VLmin = 2,
                        HLmin = 2,
                        DLmax = length(Matrix::diag(RM))-1,
                        VLmax = length(Matrix::diag(RM))-1,
                        HLmax = length(Matrix::diag(RM))-1,
                        d         = NULL,
                        theiler   = NA,
                        recurrenceTimes    = FALSE,
                        AUTO      = NULL,
                        chromatic = FALSE,
                        matrices  = FALSE){

  # For boot()
  # RP <- RP[indices,]

  if(!all(as.vector(RM)==0|as.vector(RM)==1)){stop("Matrix should be a binary (0,1) matrix!!")}

  if(recurrenceTimes){
    RM <- Matrix::Matrix(1-RM,sparse = TRUE)
  }

  if(length(d)<=2){
    suppressMessages(RM <- setTheiler(RM, theiler))
  } else {
    if(!is.null(attributes(RM)$theiler)){
      message(paste0("Value found in attribute 'theiler'... assuming a theiler window was already applied to the matrix."))
    }
  }

  if(!is.null(d)){
    if(length(d)==1){d <- seq(-d,d)}
    if(length(d)==2){d <- seq(min(d),max(d))}
  }

  B <- rp_nzdiags(RM)
  V <- Matrix::as.matrix(RM)[,colSums(Matrix::as.matrix(RM))>0]

  RPt <- Matrix::t(RM)

  H <- Matrix::as.matrix(RPt)[,colSums(Matrix::as.matrix(RPt))>0]
  rm(RPt)

  if(NCOL(B)==0|is.null(dim(B))){

    B <- matrix(0)
  }
  if(NCOL(V)==0|is.null(dim(V))){
    V <- matrix(0)
  }
  if(NCOL(H)==0|is.null(dim(H))){
    H <- matrix(0)
  }

  # TEST
  # if(invert){
  #   IDremove <- which(B[1,]==1)
  #   for(i in IDremove){
  #     l1 <- which(diff(c(B[1,IDremove[i]],B[,IDremove[i]]))!=0)
  #     if(l1 > DLmin){
  #       B[1:l1,IDremove[i]] <- 0
  #     }
  #   }
  # }

  # Get diagonal lines & pad with zeros
  diagonals   <- rbind.data.frame(rep(0,dim(B)[2]),
                                  B,
                                  rep(0,dim(B)[2])
  )

  # get nonzero vertical Lines & pad with zeros
  verticals <- rbind.data.frame(rep(x = 0, times = dim(V)[2]),
                                V,
                                rep(x= 0, times = dim(V)[2])
  )

  # get nonzero horizontal Lines & pad with zeros
  horizontals <- rbind.data.frame(rep(0,dim(H)[2]),
                                  H,
                                  rep(0,dim(H)[2])
  )

  # Get indices of line lengths
  diagonals.ind   <- tidyr::gather(diagonals,   key = "diagonal",   value = "segment")
  verticals.ind   <- tidyr::gather(verticals,   key = "vertical",   value = "segment")
  horizontals.ind <- tidyr::gather(horizontals, key = "horizontal", value = "segment")

  D <- diagonals.ind$segment
  names(D) <- paste0(diagonals.ind$diagonal,ifelse(recurrenceTimes,"DT","D"))
  V <- verticals.ind$segment
  names(V) <- paste0(verticals.ind$vertical,ifelse(recurrenceTimes,"VT","V"))
  H <- horizontals.ind$segment
  names(H) <- paste0(horizontals.ind$horizontal,ifelse(recurrenceTimes,"HT","H"))

  # Get consecutive nonzero segments from indices, their difference is the segment length
  # We added a row of 0s so we'll get sequences of -1, 1, -1
  diagonals.dist   <- sort(which(diff(D)==-1)-which(diff(D)==1))
  verticals.dist   <- sort(which(diff(V)==-1)-which(diff(V)==1))
  horizontals.dist <- sort(which(diff(H)==-1)-which(diff(H)==1))

  diagonals.dist   <- diagonals.dist[diagonals.dist%[]%c(DLmin,DLmax)]
  verticals.dist   <- verticals.dist[verticals.dist%[]%c(VLmin,VLmax)]
  horizontals.dist <- horizontals.dist[horizontals.dist%[]%c(HLmin,HLmax)]

  if(length(diagonals.dist)==0){diagonals.dist <- NA}
  if(length(verticals.dist)==0){verticals.dist <- NA}
  if(length(horizontals.dist)==0){horizontals.dist <- NA}

  return(list(diagonals.dist   = diagonals.dist,
              verticals.dist   = verticals.dist,
              horizontals.dist = horizontals.dist,
              diagonals.mat = diagonals[-c(1,NROW(diagonals)),],
              verticals.mat = verticals[-c(1,NROW(verticals)),],
              horizontals.mat = horizontals[-c(1,NROW(horizontals)),])[c(TRUE,TRUE,TRUE,matrices,matrices,matrices)]
  )
}


#' Copy Matrix Attributes
#'
#' Simple attribute copy used in `casnet` to convert between `matrix` and `Matrix` classes and back.
#'
#' @param source Source matrix
#' @param target Target matrix
#' @param source_remove Remove these attribute fields from the source before copying.
#'
#' @return The target matrix with attributes copied from the source matrix.
#' @export
#'
rp_copy_attributes <- function(source, target, source_remove = c("names", "row.names", "class","dim", "dimnames","x")){
  attrList <- attributes(source)
  attrList[source_remove] <- NULL
  attributes(target) <- c(attributes(target), attrList)
  return(target)
}

#' Check a Recurrence Matrix
#'
#' @param RM RM
#' @param checkS4 checkS4
#' @param checkAUTO checkAUTO
#' @param checkSPARSE checkSPARSE
#' @param fixS4 fixS4
#' @param fixAUTO fixAUTO
#' @param fixSPARSE fixSPARSE
#'
#' @return A checked and/or fixed matrix
#' @export
#' @keywords internal
#'
rp_checkfix <- function(RM, checkS4 = TRUE, checkAUTO = TRUE, checkSPARSE = FALSE, checkTSPARSE = FALSE, fixS4 = FALSE, fixAUTO = TRUE, fixSPARSE = FALSE, fixTSPARSE = FALSE){

  dummy <- Matrix::Matrix(matrix(c(0,0)))
  dummy <- rp_copy_attributes(source = RM, target =dummy)

  # Always check S4
  if(!checkS4){checkS4 <- TRUE}

  yesS4      <- FALSE
  yesAUTO    <- FALSE
  yesSPARSE  <- FALSE
  yesTSPARSE <- FALSE

  if(checkS4){
    yesS4 <- isS4(RM)
  }
  if(!yesS4&fixS4){
    RM <- Matrix::Matrix(RM, sparse = TRUE)
    yesS4 <- isS4(RM)
  }

  # check auto-recurrence
  if(checkAUTO){
    if(yesS4){
      yesAUTO <- Matrix::isSymmetric(RM)
    } else {
      yesAUTO <- identical(as.vector(RM[lower.tri(RM)]),as.vector(t(RM)[lower.tri(RM)]))
    }
  }
  if(fixAUTO){
    attributes(RM)$AUTO    <- yesAUTO
    attributes(dummy)$AUTO <- yesAUTO
  }

  if(checkTSPARSE){
    if(class(RM)[1]%in%names(methods::getClass("TsparseMatrix")@subclasses)){
      yesTSPARSE <- TRUE
    }
  }
  if(fixTSPARSE&!yesTSPARSE){
    if(!yesS4){
      RM <- Matrix::Matrix(RM, sparse = TRUE)
    }
    Mtype <- gsub("CMatrix","TMatrix",class(RM)[1])
    eval(parse(text=paste0("RM <- methods::as(RM,'",Mtype,"')")))
  }


  RM <- rp_copy_attributes(source = dummy, target = RM, source_remove = c("Dimnames", "i", "class","Dim", "p","x","factors"))

  return(RM)
}


#' Plot (thresholded) distance matrix as a recurrence plot
#'
#' @param RM A distance matrix or recurrence matrix, preferably generated by function [rp] or [rn].
#' @param plotDimensions Should the state vectors be plotted if they are available as attributes of RM (default = `TRUE`)
#' @param plotDimensionLegend If `plotDimensions = TRUE` plot a simple legend (default = `FALSE`)
#' @param plotMeasures Print common (C)RQA measures in the plot if the matrix is binary (default = `FALSE`)
#' @param plotRadiusRRbar The `Radius-RR-bar` is a colour-bar guide plotted with an unthresholded distance matrix indicating a number of `RR` values one would get if a certain distance threshold were chosen (default = `TRUE`)
#' @param drawGrid Draw a grid on the recurrence plot (default = `FALSE`)
#' @param drawDiagonal One usually omits the main diagonal, the Line Of Incidence (LOI) from the calculation of Auto-RQA measures., it is however common to plot it. Set to `FALSE` to omit the LOI from an Auto-Recurrence Plot (default = `TRUE`)
#' @param diagColour Colour of the main diagonal (default = `"grey50"`)
#' @param drawNA Draw `NA` values in the recurrence plot? If `FALSE` then `NA` values will be considered non-recurring (default = `FALSE`)
#' @param markEpochsLOI Pass a factor whose levels indicate different epochs or phases in the time series and use the line of identity to represent the levels by different colours (default = `NULL`)
#' @param radiusValue If `plotMeasures = TRUE` and RM is an unthresholded matrix, this value will be used to calculate recurrence measures. If `plotMeasures = TRUE` and RM is already a binary recurrence matrix, pass the radius that was used as a threshold to create the matrix for display purposes. If `plotMeasures = TRUE` and `radiusValue = NA`, function `est_radius()` will be called with default settings (find a radius that yields `.05` recurrence rate). If `plotMeasures = FALSE` this setting will be ignored.
#' @param courseGrain Reduce the size of the matrix before plotting (greatly improves speed). If `plotMeasures = TRUE` coursegraining takes place after calculation of the CRQA measures. See [mat_coursegrain()] for details (default = `TRUE`)
#' @param maxSize The maximum size of the matrix (cols x rows) above which the coursegraining function [mat_coursegrain] will be implemented if `courseGrain = TRUE`. This will speed up the plotting process
#' @param title A title for the plot
#' @param xlabel An x-axis label
#' @param ylabel An y-axis label
#' @param plotSurrogate Should a 2-panel comparison plot based on surrogate time series be added? If `RM` has attributes `y1` and `y2` containing the time series data (i.e. it was created by a call to [rp]), the following options are available: "RS" (random shuffle), "RP" (randomised phases), "AAFT" (amplitude adjusted fourier transform). If no timeseries data is included, the columns will be shuffled.  NOTE: This is not a surrogate test, just 1 surrogate is created from `y1`. (default = `FALSE`)
#' @param returnOnlyObject Return the ggplot object only, do not draw the plot (default = `TRUE`)
#'
#' @return A nice plot of the recurrence matrix.
#' @export
#'
#' @family Distance matrix operations (recurrence plot)
#'
rp_plot <- function(RM,
                    plotDimensions = FALSE,
                    plotDimensionLegend = FALSE,
                    plotMeasures = FALSE,
                    plotRadiusRRbar = TRUE,
                    drawGrid = FALSE,
                    drawDiagonal = TRUE,
                    diagColour = "grey50",
                    drawNA = FALSE,
                    markEpochsLOI = NULL,
                    radiusValue = NA,
                    courseGrain = TRUE,
                    maxSize = 1000^2,
                    title = "",
                    xlabel = "",
                    ylabel="",
                    plotSurrogate = NA,
                    returnOnlyObject = FALSE){

  useGtable <- TRUE

  if(is.null(attr(RM,"chromatic"))){
    chromatic <- FALSE
    attr(RM,"chromatic") <- chromatic
  } else {
    chromatic <- attr(RM,"chromatic")
  }

  if(!all(stats::na.exclude(as.vector(RM))%in%c(0,1))&!chromatic){
    unthresholded <- TRUE
  } else {
    unthresholded <- FALSE
  }

  #  else {
  #   # This is just to set the attribute in case
  #   if(unthresholded|plotMeasures){
  #     radiusValue <- est_radius(RM,silent = TRUE)$Radius
  #     attr(RM,"emRad") <- radiusValue
  #   }
  # }

  # Check if we need to display CRQA measures for the original matrix

  # Get Radius attribute (if present), otherwise set it (if not present)
  if(is.na(radiusValue%00%NA)){
    if(!is.na(attr(RM,"emRad")%00%NA)){
      radiusValue <- attr(RM,"emRad")
    }
  } else {
    if(is.numeric(radiusValue)&radiusValue>=0){
      if(is.na(attr(RM,"emRad")%00%NA)){
        attr(RM,"emRad") <- radiusValue
      }
    } else {
      radiusValue <- NA
    }
  } # is.na(radiusvalue)

  # check auto-recurrence and make sure Matrix has sparse triplet representation
  AUTO <- attr(RM,"AUTO")
  if(!is.na(AUTO%00%NA)){
    RM   <- rp_checkfix(RM, checkAUTO = AUTO, fixAUTO = AUTO, checkSPARSE = TRUE)
  } else {
    RM   <- rp_checkfix(RM, checkAUTO = TRUE, fixAUTO = TRUE, checkSPARSE = TRUE)
    AUTO <- attr(RM,"AUTO")
  }

  # Get CRQA measures if requested
  if(plotMeasures){
      if(!is.na(attr(RM,"emRad")%00%NA)){
        radiusValue <- attr(RM,"emRad")
      } else {
        radiusValue <- est_radius(RM, silent = TRUE)$Radius
        attr(RM,"emRad") <- radiusValue
      }

    if(unthresholded){
      rpOUT   <- rp_measures(RM, emRad = radiusValue, AUTO = AUTO)
    } else {
      rpOUT   <- rp_measures(RM, AUTO = AUTO)
    }
  } # plotMeasures


  # Size check
  reduced <- FALSE
  if((rp_size(RM=RM)$rp_size_total>=maxSize)&courseGrain){
    message("NOTE: To speed up the plotting process, the RP will represent a coursegrained matrix. Set argument 'courseGrain = FALSE' to see the full matrix.")
    if((rp_size(RM=RM)$rp_size_total/2)>=maxSize){
      target_height <- round(sqrt(maxSize))
      target_width  <- round(sqrt(maxSize))
    } else {
      target_height <- NROW(RM)/2
      target_width  <- NCOL(RM)/2
    }


    RM <- mat_coursegrain(RM, target_height = target_height, target_width = target_width)
    reduced <- TRUE
  }

  colvec <- c("#FFFFFF00","#000000")
  names(colvec) <- c("0","1")


  # Make sure we can draw a diagonal if requested
  maxDist <- max(RM, na.rm = TRUE)
  if(drawDiagonal&AUTO){
    if(all(stats::na.exclude(as.vector(RM))%in%c(0,1))){
      RM <- bandReplace(RM,0,0,1)
    } else {
      if(!chromatic){
        RM <- bandReplace(RM,0,0,(maxDist+1))
      }
    }
  }


  # prepare data
  #if(attr(RM,"package")%00%""%in%"Matrix"){
  if(any(grepl("Matrix",class(RM)))){
    RM     <- rp_checkfix(RM, checkTSPARSE = TRUE, fixTSPARSE = TRUE)
    #meltRP <- data.frame(Var1 = (RM@i+1), Var2 = (RM@j+1), value = as.numeric(RM@x))
    meltRP <-  mat_mat2ind(Matrix::as.matrix(RM))
  } else {
    meltRP <-  mat_mat2ind(as.matrix(RM))
  }

  hasNA <- FALSE

  if(chromatic){
    if(is.null(attr(RM,"chromaNames"))){
      chromaNumbers <- sort(unique(meltRP$value))[sort(unique(meltRP$value))>0]
      names(chromaNumbers) <- c("No recurrence", paste("Cat.", chromaNames))
      attr(RM,"chromaNames") <- chromaNumbers
      chromaNames <- names(chromaNumbers)
    } else {
      chromaNumbers <- attr(RM,"chromaNames")
      Cind <- sort(chromaNumbers, index.return = TRUE)
      chromaNames <- names(chromaNumbers)[Cind$ix]
      chromaNumbers <- chromaNumbers[Cind$ix]
    }
    if(NA%in%chromaNames){
      hasNA <- TRUE
    }
  }


  if(drawNA){
    if(!is.null(attr(RM,"NAij"))){
      hasNA <- TRUE
      NAid <- data.frame(attr(RM,"NAij"))
      NAid$value <- NA
      colnames(NAid)[1:2] <- c("Var1","Var2")
      #meltRP$value[meltRP$Var1==NAid[,2]&meltRP$Var2==NAid[,1]] <- NA
      meltRP <- rbind(meltRP,NAid)
    } else {
      message("Cannot find any NA to draw!")
      drawNA <- FALSE
      hasNA <- FALSE
    }
  } else {
    if(any(is.na(meltRP$value))){
      hasNA <- FALSE
      meltRP$value[is.na(meltRP$value)] <- (maxDist + 1)
    }
  }


  # Unthresholded START ----
  showL <- FALSE
  if(unthresholded){

    if(!is.null(attr(RM,"weighted"))){
      if((attr(RM,"weighted")%00%FALSE)&plotRadiusRRbar){
        warning("Can't show the radius vs. RR bar on a weighted Recurrence Plot!")
        plotRadiusRRbar <- FALSE
      }
    }

    if(!is.null(markEpochsLOI)){
      warning("Can't show epochs on an unthresholded Recurrence Plot!")
    }

    # if(chromatic){
    #   plotRadiusRRbar <- FALSE
    # }

  } else { # unthresholded

    # THRESHOLDED = TRUE
    unthresholded <- FALSE
    plotRadiusRRbar <- FALSE

    if(!is.null(markEpochsLOI)){
      if(is.factor(markEpochsLOI)&length(markEpochsLOI)==max(c(NROW(RM),NCOL(RM)))){
        start <- max(meltRP$value, na.rm = TRUE) + 1
        #cpal <- paletteer::paletteer_d(package = "rcartocolor",palette = "Safe", n = nlevels(markEpochsLOI))
        cpal <- getColours(Ncols = nlevels(markEpochsLOI))
        #cpal <- paletteer::paletteer_d(package = "ggthemes",palette = "tableau_colorblind10",n = nlevels(markEpochsLOI),direction = 1)

        if(hasNA){
          colvec <- c("#FFFFFF00","#000000","#FF0000", cpal)
          names(colvec) <- c("0","1","NA",levels(markEpochsLOI))
          vals <- 1:3
        } else {
          colvec <- c("#FFFFFF00","#000000", cpal) #viridis::viridis_pal()(nlevels(markEpochsLOI)))
          names(colvec) <- c("0","1",levels(markEpochsLOI))
          vals <- 1:2
        }
        N <- max(c(NROW(RM),NCOL(RM)))
        LOIvals <- as.numeric_discrete(markEpochsLOI, sortUnique = TRUE)
        for(i in 1:N){
          j <- i
          meltRP$value[meltRP$Var1==i&meltRP$Var2==j] <- start + LOIvals[i]
        }

        # uni <- sort(unique(meltRP$value))
        # LOIvals <- as.numeric_discrete()
        # if(length(uni)==length(names(colvec))){
        #   which()
        #
        # }
        meltRP$value <- factor(meltRP$value, labels = names(colvec)) #as.character(c(as.numeric(colvec)[vals],(LOIvals+start)))
        showL <- TRUE
      } else {
        warning("Variable passed to 'markEpochsLOI' is not a factor or doesn't have correct length.")
      }
    } else {

      if(chromatic){

        if("No recurrence"%in%chromaNames){
          chromaNames[which("No recurrence"%in%chromaNames)] <- "0"
        } else {
          if(!"0"%in%chromaNames){
          chromaNames <- c("0",chromaNames)
          }
        }

        if(!0%in%chromaNumbers){
          chromaNumbers <- c(0,chromaNumbers)
          names(chromaNumbers)[1] <- "0"
        }

        colvec <- getColours(length(chromaNames))
        colvec[which(chromaNames%in%"0")] <- "#FFFFFF00"

        if(hasNA){
          if(NA%in%chromaNames){
            colvec[which(chromaNames%in%NA)] <-"#FF0000"
            names(colvec) <- chromaNames
          } else {
            colvec <- c(colvec,"#FF0000")
            names(colvec) <- c(chromaNames,NA)
            if(!NA%in%chromaNumbers){
              chromaNumbers <- c(chromaNumbers, NA)
            }
          }
        } else {
          names(colvec) <- chromaNames
        }

        if(!drawNA){
          chromaNames   <- stats::na.exclude(chromaNames)
          chromaNumbers <- stats::na.exclude(chromaNumbers)
        }

        # CHECK VECTOR LENGTHS
        allEqual <- FALSE
        if(length(unique(meltRP$value))==length(chromaNames)){
          if(length(unique(meltRP$value))==length(chromaNumbers)){
            allEqual <- TRUE
          }
        }

        if(allEqual){
          meltRP$value <- factor(meltRP$value, levels = as.numeric(chromaNumbers), labels = chromaNames, exclude = NULL)
        } else {
          #message("Check values and chromaNames...")

          excludeVals <- unique(meltRP$value)[!unique(meltRP$value)%in%chromaNumbers]
          meltRP$value <- factor(meltRP$value,
                                 levels = c(as.numeric(chromaNumbers), excludeVals),
                                 labels = c(chromaNames, as.character(excludeVals)),
                                 exclude = c(as.character(excludeVals),""))
        }

        showL <- FALSE
        tmp <- data.frame(x = 1:length(chromaNames), y = 1:length(chromaNames), cols = chromaNames)

        grL <- ggplot(tmp, aes_(x = ~x, y = ~y)) +
          geom_point(aes_(colour = ~cols, fill = ~cols), pch=22, size =10) +
          scale_fill_manual(name  = "Phases", breaks = chromaNames,
                            values = colvec,
                            na.translate = TRUE,
                            na.value = "#FF0000") +
          scale_colour_manual(name  = "Phases", breaks = chromaNames,
                              values = colvec,
                              na.translate = TRUE,
                              na.value = "#FF0000")
        grLegend <- cowplot::get_legend(grL)
        rm(grL)

      } else {

      colvec <- c("#FFFFFF00","#000000")
      names(colvec) <- c("0","1")
      meltRP$value <- factor(meltRP$value, levels = c(0,1), labels = c("0","1"))
      }
    }
  } # unthresholded = FALSE

  ## Unthresholded END ##

  # Main plot ----
  gRP <-  ggplot2::ggplot(ggplot2::aes_(x=~Var1, y=~Var2, fill = ~value), data= meltRP) +
    ggplot2::geom_raster(hjust = 0, vjust=0, show.legend = showL)

  if(drawDiagonal){
    gRP <- gRP + ggplot2::geom_abline(slope = 1,colour = diagColour, size = 1)
  }

  ## Unthresholded START ----
  if(unthresholded){

    # #barValue <- 0.05
    # if(is.na(radiusValue)){
    #   barValue <- mean(meltRP$value, na.rm = TRUE)
    #   } else {
    #   barValue <- est_radius(RM, targetValue = radiusValue, silent = TRUE, maxIter = 100, radiusOnFail = "percentile")$Radius
    # }

    if(!is.na(radiusValue)){
      barValue <- radiusValue
    } else {
      barValue <- mean(meltRP$value, na.rm = TRUE)
    }

    if(attr(RM,"weighted")){

      gRP <- gRP + ggplot2::scale_fill_gradient(low      = "white",
                                                high     = "red3",
                                                na.value = "#FF0000",
                                                space    = "Lab",
                                                name     = "")

    } else {




      gRP <- gRP + ggplot2::scale_fill_gradient2(low      = "red3",
                                                 high     = "steelblue",
                                                 mid      = "white",
                                                 na.value = "#FF0000",
                                                 midpoint = barValue*1.1, #mean(meltRP$value, na.rm = TRUE),
                                                 limit    = c(min(meltRP$value, na.rm = TRUE),max(meltRP$value, na.rm = TRUE)),
                                                 space    = "Lab",
                                                 name     = "")

    } # Weighted


    ## RadiusRRbar START ----
    if(plotRadiusRRbar){
      # Create a custom legend
      distrange  <- round(seq(0,max(RM,na.rm = TRUE),length.out=7),2)
      resol      <- sort(unique(round(as.vector(RM),2)))
      if(length(resol)<7){
        resol <- distrange
      }
      if(length(resol)>100){
        resol <- round(seq(0,max(RM,na.rm = TRUE),length.out=100),2)
      }
      resol <- resol %>% tibble::as_tibble() %>% dplyr::mutate(y= seq(exp(0),exp(1),length.out=NROW(resol)), x=0.5)

      distrange <- plyr::ldply(c(0.001, 0.005, 0.01, 0.05, 0.1, 0.5), function(t){
        suppressWarnings(suppressMessages(est_radius(RM,targetValue = t, silent = TRUE, maxIter = 100, radiusOnFail = "percentile")))
      })

      if(AUTO){
      maxD <- suppressMessages(max(RM[lower.tri(RM)], na.rm = TRUE))
      } else {
      maxD <- suppressMessages(max(RM, na.rm = TRUE))
      }
      RecScale <- data.frame(RR=distrange$Measure,epsilon=distrange$Radius)
      RecScale <- RecScale %>%
        dplyr::add_row(epsilon=mean(c(0,distrange$Radius[1])),RR=mean(c(0,distrange$Measure[1])),.before = 1) %>%
        dplyr::add_row(epsilon=maxD,RR=1)

      resol$y <- elascer(x = resol$y,lo = min(log(RecScale$RR),na.rm = TRUE), hi = max(log(RecScale$RR),na.rm = TRUE))
      resol <- resol[-1,]

      if(!is.na(radiusValue)){
        barValue <- round(RecScale$RR[which(round(RecScale$epsilon,4)>=radiusValue)[1]],4)
        barValue <- resol$value[which(resol$y>=log(barValue))[1]]
      } else {
        barValue <- mean(meltRP$value, na.rm = TRUE)
      }

      gDist <-  ggplot2::ggplot(resol,ggplot2::aes_(x=~x,y=~y,fill=~value)) +
        ggplot2::geom_tile(show.legend = FALSE) +
        ggplot2::scale_y_continuous(name = "Recurrence Rate", breaks = log(round(RecScale$RR,3)), labels = paste(round(RecScale$RR,3)), sec.axis = dup_axis(name=expression(paste("recurrence threshold",~ epsilon)), labels = paste(round(RecScale$epsilon,3)))) +
        ggplot2::scale_fill_gradient2(low      = "red3",
                                      high     = "steelblue",
                                      mid      = "white",
                                      na.value = "grey",
                                      midpoint =  barValue * 1.1,#mean(meltRP$value, na.rm = TRUE),
                                      #limit    = c(min(RecScale$epsilon),max(RecScale$epsilon, na.rm = TRUE)),
                                      space    = "Lab",
                                      name     = "") +
        ggplot2::coord_equal(1, expand = FALSE) +
        ggplot2::theme_bw() +
        ggplot2::theme(panel.background = element_blank(),
                       panel.grid.major  = element_blank(),
                       panel.grid.minor  = element_blank(),
                       legend.background = element_blank(),
                       legend.key = element_blank(),
                       panel.border = element_blank(),
                       axis.text.y  =  element_text(size = 8),
                       # axis.text.y.left  =  element_text(size = 8),
                       # axis.text.y.right =  element_text(size = 8),
                       axis.text.x  = element_blank(),
                       axis.ticks.x = element_blank(),
                       axis.title.x = element_blank(),
                       axis.title.y = element_text(hjust = 0, size = 10),
                       # axis.title.y.left = element_text(hjust = 0, size = 10),
                       # axis.title.y.right = element_text(hjust = 0, size = 10),
                       plot.margin = margin(0,5,0,5, unit = "pt"))
    }

  }
  # else { # unthresholded
  #
  # warning("This is a binary matrix. Cannot plot the Radius vs. RR bar!")
  #
  # }
  ##RadiusRRbar END##


  ## Create Theme START ----
  rptheme <- ggplot2::theme_bw() + ggplot2::theme(panel.background = element_blank(),
                                                  panel.grid.minor  = element_blank(),
                                                  panel.border = element_rect("grey50",fill=NA),
                                                  legend.key = element_rect(colour = "grey90"),
                                                  axis.ticks = element_blank(),
                                                  axis.text = element_blank(),
                                                  plot.margin = margin(0,0,0,0))

  if(drawGrid){
    rptheme <- rptheme +  theme(panel.grid.minor  = element_blank(),
                                panel.grid.major  = element_line("grey80",size = .1),
                                panel.ontop = unthresholded)
  } else {
    rptheme <- rptheme +  theme(panel.grid.major  = element_blank(),
                                panel.grid.minor  = element_blank(),
                                panel.ontop = FALSE)
  } #drawGrid


  if(showL){
    rptheme <- rptheme +  theme(legend.position = c(1.2,1),
                                legend.direction = "vertical",
                                legend.background = element_rect(fill="grey90"),
                                legend.justification = "top")
  }

  if(plotDimensions){
    rptheme <- rptheme + theme(axis.title.y =element_blank(),
                               axis.title.x =element_blank())
  }

  if(!is.null(markEpochsLOI)){
    if(!unthresholded){
      gRP <- gRP +  ggplot2::scale_fill_manual(name  = "Key:",
                                               values = colvec,
                                               na.translate = TRUE ,
                                               na.value = scales::muted("slategray4"),
                                               guide = "legend",
                                               limits = levels(meltRP$value))
    }

  } else {
    if(!unthresholded){

      if(chromatic){

      gRP <- gRP +  ggplot2::scale_fill_manual(name  = "", breaks = sort(unique(meltRP$value)),
                                               values = colvec,
                                               na.translate = TRUE ,
                                               na.value = scales::muted("slategray4"))

      } else {

      gRP <- gRP +  ggplot2::scale_fill_manual(name  = "", breaks = c(0,1),
                                               values = colvec,
                                               na.translate = TRUE ,
                                               na.value = scales::muted("slategray4"),
                                               guide = "none")
      }
    }
  } # markEpochsLOI


  if(plyr::is.discrete(meltRP$Var1)){
    gRP <- gRP + ggplot2::scale_x_discrete(breaks=meltRP$Var1,expand = c(0,0))
  } else {
    gRP <- gRP + ggplot2::scale_x_continuous(breaks=meltRP$Var1,expand = c(0,0))
  }
  if(plyr::is.discrete(meltRP$Var2)){
    gRP <- gRP + ggplot2::scale_y_discrete(breaks=meltRP$Var2,expand = c(0,0))
  } else {
    gRP <- gRP + ggplot2::scale_y_continuous(breaks=meltRP$Var2,expand = c(0,0))
  }
  gRP <- gRP + rptheme + ggplot2::coord_equal(dim(RM)[1]/dim(RM)[2]) #coord_fixed(expand = FALSE)

  gy1 <- gg_plotHolder()
  gy2 <- gg_plotHolder()

  xdims <- ""
  ydims <- ""

  if(!is.null(attr(RM,"emDims1.name"))|nchar(xlabel)>0){
    xdims <- ifelse(nchar(xlabel)>0, xlabel, attr(RM,"emDims1.name"))
  }
  if(!is.null(attr(RM,"emDims2.name"))|nchar(ylabel)>0){
    ydims <- ifelse(nchar(ylabel)>0, ylabel, attr(RM,"emDims2.name"))
    if(AUTO){
      ydims <- xdims
    }
  }

  if(plotDimensions){

    if(!is.null(attr(RM,"emDims1"))){

      gRP <- gRP + ylab("") + xlab("")

      y1 <- data.frame(t1=attr(RM,"emDims1"))
      y2 <- data.frame(t2=attr(RM,"emDims2"))

      # Y1

      colnames(y1) <- paste0("X",1:NCOL(y1))
      y1$tm  <- 1:NROW(y1)
      y1$tmna <- 0
      y1$tmna[is.na(y1[,1])] <- y1$tm[is.na(y1[,1])]
      y1 <- tidyr::gather(y1,key="Dimension",value = "Value", -c("tm","tmna"))
      y1$Value <-  elascer(y1$Value)

      # Y2
      colnames(y2) <- paste0("Y",1:NCOL(y2))
      y2$tm <- 1:NROW(y2)
      y2$tmna <- 0
      y2$tmna[is.na(y2[,1])] <- y2$tm[is.na(y2[,1])]
      y2 <- tidyr::gather(y2,key="Dimension",value = "Value", -c("tm","tmna"))
      y2$Value <-  elascer(y2$Value)


    } else {
      warning("No atrribute with dimensions!")
      if(unthresholded){
        y1 <- data.frame(Value=diag(RM),x=1:length(diag(RM)), Dimension = rep("LOS",length(diag(RM))))
        y2 <- data.frame(Value=diag(RM),x=1:length(diag(RM)), Dimension = rep("LOS",length(diag(RM))))
        xdims <- paste("LOS",xdims)
        ydims <- paste("LOS",ydims)
      } else {
        plotDimensions <- FALSE
      }
    }
  } else {
    gRP <- gRP + ylab(ydims) + xlab(xdims)
  }

  if(plotDimensions){

    gy1 <- ggplot2::ggplot(y1, ggplot2::aes_(y=~Value, x= ~tm,  group= ~Dimension)) +
      ggplot2::geom_line(aes_(colour=~Dimension), show.legend = FALSE) +
      ggplot2::xlab(xdims) + ggplot2::ylab(" ") +
      ggplot2::geom_vline(ggplot2::aes_(xintercept = ~tmna), colour = scales::muted("slategray4"),alpha=.1, size=.5) +
      ggplot2::scale_color_grey() +
      ggplot2::scale_x_continuous(expand = c(0,0)) +
      ggplot2::scale_y_continuous(expand = c(0,0)) +
      ggplot2::theme(panel.background = element_blank(),
                     panel.grid.major  = element_blank(),
                     panel.grid.minor  = element_blank(),
                     legend.background = element_blank(),
                     legend.key = element_blank(),
                     panel.border = element_blank(),
                     axis.text = element_blank(),
                     axis.line = element_blank(),
                     axis.ticks = element_blank(),
                     axis.title.x =element_text(colour = "black",angle = 0, vjust = +3),
                     axis.title.y =element_blank(),
                     plot.margin = margin(0,0,0,0, unit = "pt")) +
      ggplot2::coord_cartesian(expand = FALSE)  # +  coord_fixed(1/10)

    gy2 <- ggplot2::ggplot(y2, ggplot2::aes_(y=~Value, x=~tm, group=~Dimension)) +
      ggplot2::geom_line(ggplot2::aes_(colour=~Dimension), show.legend = FALSE) +
      ggplot2::ylab(" ") + ggplot2::xlab(ydims) +
      ggplot2::geom_vline(ggplot2::aes_(xintercept = ~tmna), colour = scales::muted("slategray4"),alpha=.1, size=.5) +
      ggplot2::scale_color_grey() +
      ggplot2::scale_x_continuous(expand = c(0,0)) +
      ggplot2::theme(panel.background = element_blank(),
                     panel.grid.major  = element_blank(),
                     panel.grid.minor  = element_blank(),
                     legend.background = element_blank(),
                     legend.key = element_blank(),
                     panel.border = element_blank(),
                     axis.text = element_blank(),
                     axis.line = element_blank(),
                     axis.ticks = element_blank(),
                     axis.title.x =element_blank(),
                     axis.title.y =element_text(colour = "black",angle = 90, vjust = -2),
                     plot.margin = margin(0,0,0,0, unit = "pt")) +
      ggplot2::coord_flip(expand = FALSE) +
      ggplot2::scale_y_reverse(expand = c(0,0))


    if(plotDimensionLegend){

      # y1 <- data.frame(t1=attr(RM,"emDims1"))
      # y2 <- data.frame(t2=attr(RM,"emDims2"))

      y1l <- data.frame(Value     = rep(0,length(unique(y1$Dimension)),each=2),
                        tm        = rep(0,length(unique(y1$Dimension)),each=2),
                        Dimension =  rep(unique(y1$Dimension),each=2))

      gdl <- ggplot2::ggplot(y1l, ggplot2::aes_(y=~Value, x=~tm, colour=~Dimension)) +
        ggplot2::geom_line(show.legend = TRUE, size = 2) + ggplot2::scale_color_grey() +
        scale_x_continuous() +
        theme_void() + theme(legend.text = element_text(size=rel(2)),
                             legend.title = element_text(size=rel(2)),
                             legend.position = c(.5,.5))


    }

  } # plotdimensions


  if(plotMeasures){

    rpOUT    <- round(rpOUT,3)
    if(is.na(rpOUT$emRad)){
      rpOUT$Radius <- round(radiusValue,3)
    } else {
      rpOUT$Radius <- rpOUT$emRad
    }


    rpOUTdat <- rpOUT %>%
      dplyr::select(dplyr::one_of(c("Radius","RP_N","RR","DET","MEAN_dl","ENT_dl","LAM_vl","TT_vl","ENT_vl"))) %>%
      tidyr::gather(key="measure",value="value") %>%
      dplyr::mutate(x=rep(0,9),y=9:1)

    #rpOUTdat <- cbind(rpOUTdat,rpOUTdat)
    rpOUTdat$label <-  paste0(rpOUTdat$measure,":\n",rpOUTdat$value)

    gA <-ggplot2::ggplot(rpOUTdat,ggplot2::aes_(x=~x,y=~y)) +
      ggplot2::geom_text(ggplot2::aes_(label=~label), family="mono", hjust="left", vjust="center", size=3, parse = FALSE) +
      #scale_x_continuous(limits = c(0,.3)) +
      ggplot2::theme_void() +
      ggplot2::theme(plot.margin = margin(0,5,0,5, unit = "pt"))

    #geom="text", label = paste("Radius:",rpOUT$Radius,"\nRec points:",rpOUT$RT,"\nRR",rpOUT$RR,"\nDET:",rpOUT$DET,"\nMEAN_dl:",rpOUT$MEAN_dl,"\nENT_dl:",rpOUT$ENT_dl,"\nLAM_vl:",rpOUT$LAM_vl, "\nTT_vl:",rpOUT$TT_vl,"\nENTR_vl:",rpOUT$ENT_vl)) + theme_minimal() + theme(text = element_text(family = "mono"))
    # ,"\nLAM_hl:",rpOUT$LAM_vl, "| TT_hl:",rpOUT$TT_vl,"| ENTR_hl:",rpOUT$ENT_hl))
  }

  # Build Graph using gtable

  g <- ggplot2::ggplotGrob(gRP)

  if(plotDimensions){
    gry2<-ggplot2::ggplotGrob(gy2)
    gry1<-ggplot2::ggplotGrob(gy1)
    if(plotDimensionLegend){
      grDimL <- ggplot2::ggplotGrob(gdl)
    }
  }

  if(unthresholded&plotRadiusRRbar){
    grDist <- ggplot2::ggplotGrob(gDist)
  }

  if(plotMeasures){
    grA <- ggplot2::ggplotGrob(gA)
  }

  # cat(paste("\n\nplotMeasures =", plotMeasures,"\n"))
  # cat(paste("plotDimensions =", plotDimensions,"\n"))
  # cat(paste("plotDimensionLegend =", plotDimensionLegend,"\n"))
  # cat(paste("unthresholded =", unthresholded,"\n"))
  # cat(paste("plotRadiusRRbar =", plotRadiusRRbar,"\n"))

  # Build list for the gtable
  Lmat <- list()
  if(plotMeasures){Lmat[[1]] <- grA} else {Lmat[[1]] <- NA}
  if(plotMeasures){Lmat[[2]] <- grid::nullGrob()} else {Lmat[[2]] <- NA}
  if(plotDimensions){Lmat[[3]] <- gry2} else {Lmat[[3]] <- NA}
  if(plotDimensionLegend){Lmat[[4]] <- grDimL} else {if(plotDimensions){Lmat[[4]] <- grid::nullGrob()} else {Lmat[[4]] <- NA}}
  Lmat[[5]] <- g
  if(plotDimensions){Lmat[[6]] <- gry1} else {Lmat[[6]] <- NA}
  if((unthresholded&plotRadiusRRbar)|chromatic){
    if(chromatic){
      Lmat[[7]] <- grLegend
    } else {
      Lmat[[7]] <- grDist
    }
  } else {
    Lmat[[7]] <- NA
  }
  if(plotDimensions&(unthresholded&plotRadiusRRbar|chromatic)){Lmat[[8]] <- grid::nullGrob()} else {Lmat[[8]] <- NA}

  # Remove unused fields
  Lmat[is.na(Lmat)] <- NULL

  w <- c(.35,.25, 1,.5,.35)[c(plotMeasures,(plotMeasures|plotDimensions),TRUE,plotRadiusRRbar,chromatic)]
  h <- c(1,.25)[c(TRUE,plotDimensions)]

  widths  <- unit(w,"null")  #ifelse(plotMeasures, unit(c(.35,.25, 1,.5),  unit(c(.25, 1), "null")))
  heights <- unit(h,"null") # ifelse(plotMeasures, unit(c(1), "null"), unit(c(1,.25), "null"))

  if(plotDimensions){
    mrows <- 2
  } else {
    mrows <- 1
  }

  mat <- matrix(Lmat,nrow = mrows)
  gt  <- gtable::gtable_matrix("rp", mat, widths = widths, heights =  heights, respect = TRUE)

  #  gtable::gtable_matrix("bi_rp", mat, widths = unit(c(1), "null"), heights =  unit(c(1), "null"), respect = TRUE)

  # if(plotDimensions&plotDimensionLegend&!plotMeasures&unthresholded&!plotRadiusRRbar){
  #   mat <- matrix(list(gry2, grDimL, g, gry1),nrow = 2)
  #   gt  <- gtable::gtable_matrix("di_rp_dim", mat, widths = unit(c(.25, 1), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&plotDimensionLegend&!plotMeasures&unthresholded&plotRadiusRRbar){
  #   mat <- matrix(list(gry2, grDimL, g, gry1, grDist, grid::nullGrob()),nrow = 2)
  #   gt  <- gtable::gtable_matrix("di_rp_dim", mat, widths = unit(c(.25, 1,.5), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&!plotDimensionLegend&!plotMeasures&unthresholded&plotRadiusRRbar){
  #   mat <- matrix(list(gry2, grid::nullGrob(),g, gry1, grDist, grid::nullGrob()),nrow = 2)
  #   gt  <- gtable::gtable_matrix("di_rp_dim", mat, widths = unit(c(.25, 1,.5), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&!plotMeasures&!plotDimensionLegend&unthresholded&!plotRadiusRRbar){
  #   mat <- matrix(list(gry2, grid::nullGrob(),g, gry1),nrow = 2)
  #   gt  <- gtable::gtable_matrix("di_rp_dim", mat, widths = unit(c(.25, 1), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&!plotDimensionLegend&!plotMeasures&!unthresholded){
  #   mat <- matrix(list(gry2, grid::nullGrob(),g, gry1),nrow = 2)
  #   gt  <- gtable::gtable_matrix("bi_rp_dim", mat, widths = unit(c(.25, 1), "null"), heights =  unit(c(1, .25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&!plotDimensionLegend&plotMeasures&unthresholded&plotRadiusRRbar){
  #   mat <- matrix(list(grA, grid::nullGrob(), gry2, grid::nullGrob(),g, gry1, grDist, grid::nullGrob()),nrow = 2)
  #   gt  <- gtable::gtable_matrix("di_rp_dim_meas", mat, widths = unit(c(.35,.25, 1,.5), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&!plotDimensionLegend&plotMeasures&unthresholded&!plotRadiusRRbar){
  #   mat <- matrix(list(grA, grid::nullGrob(), gry2, grid::nullGrob(),g, gry1),nrow = 2)
  #   gt  <- gtable::gtable_matrix("di_rp_dim_meas", mat, widths = unit(c(.35,.25, 1), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(plotDimensions&!plotDimensionLegend&plotMeasures&!unthresholded){
  #   mat <- matrix(list(grA, grid::nullGrob(), gry2, grid::nullGrob(),g, gry1),nrow = 2)
  #   gt  <- gtable::gtable_matrix("bi_rp_dim_meas", mat, widths = unit(c(.35,.25, 1), "null"), heights =  unit(c(1,.25), "null"),respect = TRUE)
  # }
  #
  # if(!plotDimensions&!plotDimensionLegend&plotMeasures&unthresholded&plotRadiusRRbar){
  #   mat <- matrix(list(grA, g, grDist),nrow = 1)
  #   gt  <- gtable::gtable_matrix("di_rp_meas", mat, widths = unit(c(.35, 1,.5), "null"), heights =  unit(c(1), "null"),respect = TRUE)
  # }
  #
  # if(!plotDimensions&!plotDimensionLegend&plotMeasures&unthresholded&!plotRadiusRRbar){
  #   mat <- matrix(list(grA, g),nrow = 1)
  #   gt  <- gtable::gtable_matrix("di_rp_meas", mat, widths = unit(c(.35, 1), "null"), heights =  unit(c(1), "null"),respect = TRUE)
  # }
  #
  # if(!plotDimensions&!plotDimensionLegend&plotMeasures&!unthresholded){
  #   mat <- matrix(list(grA, g),nrow = 1)
  #   gt  <- gtable::gtable_matrix("bi_rp_meas", mat, widths = unit(c(.35, 1), "null"), heights =  unit(c(1), "null"),respect = TRUE)
  # }
  #
  # if(!plotDimensions&!plotDimensionLegend&!plotDimensionLegend&!plotMeasures&unthresholded&plotRadiusRRbar){
  #   mat <- matrix(list(g, grDist),nrow = 1)
  #   gt  <- gtable::gtable_matrix("di_rp", mat, widths = unit(c(1,.5), "null"), heights =  unit(c(1), "null"),respect = TRUE)
  # }
  #
  # if(!plotDimensions&!plotDimensionLegend&!plotMeasures&unthresholded&!plotRadiusRRbar){
  #   mat <- matrix(list(g),nrow = 1)
  #   gt  <- gtable::gtable_matrix("di_rp", mat, widths = unit(c(1), "null"), heights =  unit(c(1), "null"),respect = TRUE)
  # }
  #
  # if(!plotDimensions&!plotDimensionLegend&!plotMeasures&!unthresholded){
  #   mat <- matrix(list(g),nrow = 1)
  #   gt  <- gtable::gtable_matrix("bi_rp", mat, widths = unit(c(1), "null"), heights =  unit(c(1), "null"),respect = TRUE)
  # }

  if(reduced){
    title <- paste(title,"(coursegrained matrix)")
  }

  if(nchar(title)>0){
    grT <- ggplot2::ggplot(data.frame(x=1,y=1)) +
      ggplot2::geom_text(ggplot2::aes_(x=~x,y=~y), label=title) +
      theme_void() +
      theme(plot.margin = margin(0,0,0,0, unit = "pt"))
    gt  <- gtable::gtable_add_rows(x = gt, heights =  unit(c(.1),"null"), pos=0)
    l <- sum(c(plotDimensions,plotMeasures))+1
    gt  <- gtable::gtable_add_grob(x = gt, grobs = ggplot2::ggplotGrob(grT), name = "Title", t=1, l=l)
  }

  g <- gtable::gtable_add_padding(gt, unit(5, "pt"))


  if(!returnOnlyObject){
    if(useGtable){
      grid::grid.newpage()
      grid::grid.draw(g)
    } else {
      # graphics::plot.new()
      graphics::plot(g)
    }
  }
  return(invisible(g))
}


#' rp_size
#'
#' Calculate the maximum possible number of recurrent points in a recurrence matrix.
#'
#' This function can take into account the presence of a `theiler` window, that is the points in the window will be excluded from the calculation. For example, some scholars will exclude the main diagonal from the calculation of the recurrence rate.
#'
#' @param RM A Matrix object
#' @param dims Two element vector representing the dimensions of Matrix `RM`. If `dims` is provided, the Matrix does not have to be passed as an argument (default = `NA`)
#' @param AUTO Is the Matrix an Auto Recurrence Matrix? If so, the length of the diagonal will be subtracted from the matrix size, pass `FALSE` to prevent this behaviour. If `NULL` (default) `AUTO` will take on the value of `isSymmetric(RM)`.
#' @param theiler Should a Theiler window be applied?
#'
#' @return Matrix size for computation of recurrence measures.
#' @export
#'
#' @family Distance matrix operations (recurrence plot)
#'
#' @examples
#' # Create a 10 by 10 matrix
#' library(Matrix)
#' m <- Matrix(rnorm(10),10,10)
#'
#' rp_size(RM = m, AUTO = TRUE, theiler = 0)  # Subtract diagonal
#' rp_size(RM = m, AUTO = FALSE,theiler = 0)  # Do not subtract diagonal
#' rp_size(RM = m, AUTO = NULL, theiler = 0)  # Matrix is symmetrical, AUTO is set to TRUE
#' rp_size(RM = m, AUTO = NULL, theiler = 1)  # Subtract a Theiler window of 1 around and including the diagonal
#'
#' # Calculate without a matrix
#' rp_size(dims = c(10,10), AUTO = TRUE, TRUE,0)
#' rp_size(dims = c(10,10), AUTO = FALSE,FALSE,0)
#'
rp_size <- function(RM = NULL, dims = NULL, AUTO = NULL, theiler = NULL){

  if(is.null(AUTO)){
    if(!is.null(RM)){
      AUTO <- Matrix::isSymmetric(RM)
    } else {
      stop("Please provide a value for AUTO (TRUE/FALSE).")
    }
  }

  if(is.null(RM)&is.null(dims)){
    stop("Need either a Matrix (RM) object or its dimensions (dims)!")
  }

  if(!is.null(dims)&is.null(RM)){
    if(length(dims)!=2){
      stop("dims must be a vector of length 2.")
    }
  }

  if(is.null(theiler)){
    if(!is.null(RM)){
     if(is.null(attributes(RM)$theiler)){
      if(AUTO){
        theiler <- 1
      } else {
        theiler <- 0
      }
    } else {
      theiler <- attributes(RM)$theiler
    }
   }
  }

  if(is.null(dims)){
    if(!is.null(RM)){
      dims <- dim(RM)
      RMdiag <- length(Matrix::diag(RM))
    }
  } else {
    RMdiag <- max(dims, na.rm = TRUE)
  }

  R_N <- cumprod(dims)[2]
  minDiag <- 0

  if(is.na(theiler%00%NA)){
    if(AUTO){
      theiler <- 1
    } else {
      theiler <- 0
    }
  }

  if(length(theiler)==1){
    if(theiler==1){
      minDiag <- RMdiag
    }
    if(theiler>1){
      minDiag <- sum(rep(RMdiag,length(seq(-theiler,theiler))) - abs(seq(-theiler,theiler)))
    }
  }

  if(length(theiler)==2){
    minDiag <- sum(rep(RMdiag,length(seq(min(theiler),max(theiler)))) - abs(seq(min(theiler),max(theiler))))
  }

  if(length(theiler)>2){
    minDiag <- sum(rep(RMdiag,length(theiler)) - abs(theiler))
  }

  return(list(rp_size_total = R_N, rp_size_theiler = R_N - minDiag))
  #cumprod(dim(mat))[2] - ifelse((includeDiag&theiler==0),length(Matrix::diag(mat)),ifelse(theiler>0,Matrix::nnzero(Matrix::band(mat,-theiler,theiler)),0)))
}


#' Empty results vector
#'
#' @return an empty rp_measures
#' @keywords internal
#' @export
#'
rp_empty <- function(){
  data.frame(
    emRad   = NA,
    RP_N     = NA,
    RR       = NA,
    DET      = NA,
    MEAN_dl  = NA,
    MAX_dl   = NA,
    ENT_dl   = NA,
    ENTrel_dl= NA,
    REP_av  = NA,
    CoV_dl   = NA,
    DIV_dl   = NA,
    SING_dl  = NA,
    N_dl     = NA,
    ANI      = NA,
    LAM_vl   = NA,
    TT_vl    = NA,
    MAX_vl   = NA,
    ENT_vl   = NA,
    ENTrel_vl= NA,
    CoV_vl   = NA,
    REP_vl   = NA,
    DIV_vl   = NA,
    SING_vl  = NA,
    N_vl     = NA,
    LAM_hl   = NA,
    TT_hl    = NA,
    MAX_hl   = NA,
    ENT_hl   = NA,
    ENTrel_hl= NA,
    CoV_hl   = NA,
    REP_hl   = NA,
    DIV_hl   = NA,
    SING_hl  = NA,
    N_hl     = NA)
}

#' rp_calc
#'
#' @inheritParams rp_measures
#'
#' @return CRQA measures and matrices of line distributions (if requested)
#' @export
#' @keywords internal
#'
rp_calc <- function(RM,
                    emRad = NULL,
                    DLmin = 2,
                    VLmin = 2,
                    HLmin = 2,
                    DLmax = length(Matrix::diag(RM)),
                    VLmax = length(Matrix::diag(RM)),
                    HLmax = length(Matrix::diag(RM)),
                    theiler   = NA,
                    AUTO      = NULL,
                    chromatic = FALSE,
                    anisotropyHV = FALSE,
                    asymmetryUL = FALSE,
                    returnUL = FALSE,
                    recurrenceTimes = FALSE,
                    matrices  = FALSE){



  RM <- rp_checkfix(RM, checkAUTO = TRUE, fixAUTO = TRUE)

  if(is.null(AUTO)){
    AUTO <- attributes(RM)$AUTO
  }

  if(is.na(theiler)){
    if(is.na(attributes(RM)$theiler%00%NA)){
      RM <- setTheiler(RM, theiler = NA)
    }
    theiler <- attributes(RM)$theiler
  }

  recmatsize <- rp_size(RM = RM, AUTO = AUTO, theiler = theiler)

  if(!is.null(attributes(RM)$emRad)){
    emRad <- attributes(RM)$emRad
  } else {
    emRad <- NA
  }

  if(!all(as.vector(RM)==0|as.vector(RM)==1)){
    if(!chromatic){
      stop("Need a thresholded recurrence matrix")
    } else {
      stop("Chromatic not yet implemented")
    }
  }

  # Total nr. recurrent points
  RP_N <- Matrix::nnzero(RM, na.counted = FALSE)

  #Proportion recurrence / Recurrence Rate
  RR <- RP_N/recmatsize$rp_size_theiler

  if(length(RR)==0){RR<-0}

  if(RR==1){
    warning("Everything is recurring!\nReturning empty vector")
    return(rp_empty())
  }

  lineMeasures_global <- rp_calc_lineMeasures(RM = RM,
                                              RP_N = RP_N,
                                              DLmin = DLmin, DLmax = DLmax,
                                              VLmin = VLmin, VLmax = VLmax,
                                              HLmin = HLmin, HLmax = HLmax,
                                              theiler = theiler,
                                              recurrenceTimes = recurrenceTimes,
                                              AUTO = AUTO,
                                              chromatic = chromatic,
                                              matrices = matrices)

  if(matrices){
    lineMatrices_global <- lineMeasures_global$crqaMatrices
    lineMeasures_global <- lineMeasures_global$crqaMeasures
  }

  # Singularities
  SING <- rp_lineDist(RM,
                      DLmin = 1, DLmax = DLmax,
                      VLmin = 1, VLmax = VLmax,
                      HLmin = 1, HLmax = HLmax,
                      theiler = theiler, AUTO = AUTO)

  # table(SING_N$verticals.dist)
  # table(SING_N$horizontals.dist)
  SING_N <- table(SING$diagonals.dist)[1]
  SING_rate <- SING_N / RP_N

  # H/V Anisotropy ratios
  if(anisotropyHV){
    out_hv_ani <- data.frame(
      Nlines_ani   = ((lineMeasures_global$N_hlp - lineMeasures_global$N_vlp)  / (lineMeasures_global$N_hlp  + lineMeasures_global$N_vlp))%00%NA,
      LAM_ani  = (lineMeasures_global$LAM_hl - lineMeasures_global$LAM_vl) / (lineMeasures_global$LAM_hl + lineMeasures_global$LAM_vl),
      MEAN_ani = (lineMeasures_global$TT_hl  - lineMeasures_global$TT_vl)  / (lineMeasures_global$TT_hl  + lineMeasures_global$TT_vl),
      MAX_ani  = (lineMeasures_global$MAX_hl - lineMeasures_global$MAX_vl) / (lineMeasures_global$MAX_hl + lineMeasures_global$MAX_vl),
      ENT_ani  = (lineMeasures_global$ENT_hl - lineMeasures_global$ENT_vl) / (lineMeasures_global$ENT_hl + lineMeasures_global$ENT_vl)
    )
  } else {
    out_hv_ani <- NA
  }


  # U/L Anisotropy ratios
  if(asymmetryUL){
    RMu <- RM
    RMu[lower.tri(RMu)] <- 0

    recmatsize_u <- rp_size(RM = RMu, AUTO = AUTO, theiler = theiler)

    lineMeasures_upper <- rp_calc_lineMeasures(RM = RMu,
                                               RP_N = RP_N,
                                               DLmin = DLmin, DLmax = DLmax,
                                               VLmin = VLmin, VLmax = VLmax,
                                               HLmin = HLmin, HLmax = HLmax,
                                               theiler = theiler,
                                               recurrenceTimes = recurrenceTimes,
                                               AUTO = AUTO,
                                               chromatic = chromatic,
                                               matrices = FALSE)

    # Total nr. recurrent points in upper
    lineMeasures_upper$RP_N <- Matrix::nnzero(RMu, na.counted = FALSE)

    #Proportion recurrence / Recurrence Rate un upper
    lineMeasures_upper$RR <- lineMeasures_upper$RP_N / recmatsize_u$rp_size_theiler

    if(length(lineMeasures_upper$RR)==0){lineMeasures_upper$RR<-0}


    SING_u <- rp_lineDist(RMu,
                          DLmin = 1, DLmax = DLmax,
                          VLmin = 1, VLmax = VLmax,
                          HLmin = 1, HLmax = HLmax,
                          theiler = theiler, AUTO = AUTO)

    lineMeasures_upper$SING_N  <-  table(SING_u$diagonals.dist)[1]

    rm(RMu, recmatsize_u, SING_u)

    RMl <- RM
    RMl[upper.tri(RMl)] <- 0

    recmatsize_l <- rp_size(RM = RMl, AUTO = AUTO, theiler = theiler)

    lineMeasures_lower <- rp_calc_lineMeasures(RM = RMl,
                                               RP_N = RP_N,
                                               DLmin = DLmin, DLmax = DLmax,
                                               VLmin = VLmin, VLmax = VLmax,
                                               HLmin = HLmin, HLmax = HLmax,
                                               theiler = theiler,
                                               AUTO = AUTO,
                                               recurrenceTimes = recurrenceTimes,
                                               chromatic = chromatic,
                                               matrices = FALSE)


    # Total nr. recurrent points in lower
    lineMeasures_lower$RP_N <- Matrix::nnzero(RMl, na.counted = FALSE)

    #Proportion recurrence / Recurrence Rate in lower
    lineMeasures_lower$RR <- lineMeasures_lower$RP_N / recmatsize_l$rp_size_theiler

    if(length(lineMeasures_lower$RR)==0){lineMeasures_lower$RR<-0}

    SING_l <- rp_lineDist(RMl,
                          DLmin = 1, DLmax = DLmax,
                          VLmin = 1, VLmax = VLmax,
                          HLmin = 1, HLmax = HLmax,
                          theiler = theiler, AUTO = AUTO)

    lineMeasures_lower$SING_N  <-  table(SING_l$diagonals.dist)[1]
    rm(RMl,SING_l)

    # RATIOs
    out_ul_asym <- data.frame(
      Npoints_asym = (lineMeasures_upper$RP_N   - lineMeasures_lower$RP_N)/ (lineMeasures_upper$RP_N  + lineMeasures_lower$RP_N),
      NDlines_asym = ((lineMeasures_upper$N_dl - lineMeasures_lower$N_dl) / (lineMeasures_upper$N_dl  + lineMeasures_lower$N_dl))%00%NA,
      NHlines_asym = ((lineMeasures_upper$N_hl - lineMeasures_lower$N_hl) / (lineMeasures_upper$N_hl  + lineMeasures_lower$N_hl))%00%NA,
      NVlines_asym = ((lineMeasures_upper$N_vl - lineMeasures_lower$N_vl) / (lineMeasures_upper$N_vl  + lineMeasures_lower$N_vl))%00%NA,
      NDpoints_asym = ((lineMeasures_upper$N_dlp - lineMeasures_lower$N_dlp)/ (lineMeasures_upper$N_dlp  + lineMeasures_lower$N_dlp))%00%NA,
      NHpoints_asym = ((lineMeasures_upper$N_hlp - lineMeasures_lower$N_hlp)/ (lineMeasures_upper$N_hlp  + lineMeasures_lower$N_hlp))%00%NA,
      NVpoints_asym = ((lineMeasures_upper$N_vlp - lineMeasures_lower$N_vlp)/ (lineMeasures_upper$N_vlp  + lineMeasures_lower$N_vlp))%00%NA,
      RR_asym      = ((lineMeasures_upper$RR    - lineMeasures_lower$RR)    / (lineMeasures_upper$RR     + lineMeasures_lower$RR))%00%NA,
      SING_N_asym  = (lineMeasures_upper$SING_N - lineMeasures_lower$SING_N)/ (lineMeasures_upper$SING_N + lineMeasures_lower$SING_N),
      DIV_asym     = (lineMeasures_upper$DIV_dl - lineMeasures_lower$DIV_dl)/ (lineMeasures_upper$DIV_dl + lineMeasures_lower$DIV_dl),
      REP_asym     = (lineMeasures_upper$REP_av - lineMeasures_lower$REP_av)/ (lineMeasures_upper$REP_av + lineMeasures_lower$REP_av),
      DET_asym     = (lineMeasures_upper$DET    - lineMeasures_lower$DET)   / (lineMeasures_upper$DET    + lineMeasures_lower$DET),
      LAM_hl_asym  = (lineMeasures_upper$LAM_hl - lineMeasures_lower$LAM_hl)/ (lineMeasures_upper$LAM_hl + lineMeasures_lower$LAM_hl),
      LAM_vl_asym  = (lineMeasures_upper$LAM_vl - lineMeasures_lower$LAM_vl)/ (lineMeasures_upper$LAM_vl + lineMeasures_lower$LAM_vl),
      MEAN_dl_asym = (lineMeasures_upper$MEAN_dl- lineMeasures_lower$MEAN_dl)/ (lineMeasures_upper$MEAN_dl+ lineMeasures_lower$MEAN_dl),
      MEAN_hl_asym = (lineMeasures_upper$TT_hl  - lineMeasures_lower$TT_hl) / (lineMeasures_upper$TT_hl  + lineMeasures_lower$TT_hl),
      MEAN_vl_asym = (lineMeasures_upper$TT_vl  - lineMeasures_lower$TT_vl) / (lineMeasures_upper$TT_vl  + lineMeasures_lower$TT_vl),
      MAX_dl_asym  = (lineMeasures_upper$MAX_dl - lineMeasures_lower$MAX_dl)/ (lineMeasures_upper$MAX_dl + lineMeasures_lower$MAX_dl),
      MAX_hl_asym  = (lineMeasures_upper$MAX_hl - lineMeasures_lower$MAX_hl)/ (lineMeasures_upper$MAX_hl + lineMeasures_lower$MAX_hl),
      MAX_vl_asym  = (lineMeasures_upper$MAX_vl - lineMeasures_lower$MAX_vl)/ (lineMeasures_upper$MAX_vl + lineMeasures_lower$MAX_vl),
      ENT_dl_asym  = (lineMeasures_upper$ENT_dl - lineMeasures_lower$ENT_dl)/ (lineMeasures_upper$ENT_dl + lineMeasures_lower$ENT_dl),
      ENT_hl_asym  = (lineMeasures_upper$ENT_hl - lineMeasures_lower$ENT_hl)/ (lineMeasures_upper$ENT_hl + lineMeasures_lower$ENT_hl),
      ENT_vl_asym  = (lineMeasures_upper$ENT_vl - lineMeasures_lower$ENT_vl)/ (lineMeasures_upper$ENT_vl + lineMeasures_lower$ENT_vl)
    )


  }

  #Output
  out_global <- data.frame(emRad      = emRad,
                           RP_max     = recmatsize$rp_size_theiler,
                           RR         = RR,
                           SING_N     = SING_N,
                           SING_rate  = SING_rate)


  if(asymmetryUL){
    if(returnUL){
    out_UL <- data.frame(upper.tri = lineMeasures_upper,
                         lower.tri = lineMeasures_lower,
                         ratios.ul = out_ul_asym)
    } else {
    out_UL <- data.frame(ratios.ul = out_ul_asym)
    }
  } else {
    out_UL <- NA
  }

  if(anisotropyHV){
    out_HL <- data.frame(ratios.hv = out_hv_ani)
  } else {
    out_HL <- NA
  }

  out <- as.data.frame(list(out_global, lineMeasures_global, out_UL, out_HL)[c(TRUE, TRUE, asymmetryUL, anisotropyHV)])

  if(matrices){
    return(list(crqaMeasures = out,
                crqaMatrices = lineMatrices_global)
    )
  } else {
    return(out)
  }
}


rp_calc_lineMeasures <- function(RM,
                                 RP_N,
                                 DLmin = 2,
                                 VLmin = 2,
                                 HLmin = 2,
                                 DLmax = length(Matrix::diag(RM)),
                                 VLmax = length(Matrix::diag(RM)),
                                 HLmax = length(Matrix::diag(RM)),
                                 d         = NULL,
                                 theiler   = NA,
                                 recurrenceTimes    = FALSE,
                                 AUTO      = NULL,
                                 chromatic = FALSE,
                                 matrices  = FALSE){

  #Get line segments
  # if(Matrix::nnzero(RM)>0)
  lineSegments <- rp_lineDist(RM,
                              DLmin = DLmin, DLmax = DLmax,
                              VLmin = VLmin, VLmax = VLmax,
                              HLmin = HLmin, HLmax = HLmax,
                              d = d, theiler = theiler,
                              recurrenceTimes = recurrenceTimes, AUTO = AUTO,
                              chromatic = chromatic,
                              matrices = matrices)

  dlines <- lineSegments$diagonals.dist%00%0
  vlines <- lineSegments$verticals.dist%00%0
  hlines <- lineSegments$horizontals.dist%00%0

  #Frequency tables of line lengths
  freq_dl <- table(dlines)
  freq_vl <- table(vlines)
  freq_hl <- table(hlines)
  freq_hv <- table(c(hlines,vlines))

  freqvec_dl <- as.numeric(names(freq_dl))
  freqvec_vl <- as.numeric(names(freq_vl))
  freqvec_hl <- as.numeric(names(freq_hl))
  freqvec_hv <- as.numeric(names(freq_hv))


  # Number of lines
  N_dl <- sum(freq_dl, na.rm = TRUE)%00%0
  N_vl <- sum(freq_vl, na.rm = TRUE)%00%0
  N_hl <- sum(freq_hl, na.rm = TRUE)%00%0
  N_hv <- sum(freq_hv, na.rm = TRUE)%00%0

  #Number of recurrent points on diagonal, vertical and horizontal lines
  N_dlp <- sum(freqvec_dl*freq_dl, na.rm = TRUE)
  N_vlp <- sum(freqvec_vl*freq_vl, na.rm = TRUE)
  N_hlp <- sum(freqvec_hl*freq_hl, na.rm = TRUE)
  N_hvp <- sum(freqvec_hv*freq_hv, na.rm = TRUE)

  #Determinism / Horizontal and Vertical Laminarity
  DET    <- N_dlp/RP_N
  LAM_vl <- N_vlp/RP_N
  LAM_hl <- N_hlp/RP_N
  LAM_hv <- N_hvp/(RP_N*2)

  #Array of probabilities that a certain line length will occur (all >1)
  P_dl <- freq_dl/N_dl
  P_vl <- freq_vl/N_vl
  P_hl <- freq_hl/N_hl
  P_hv <- freq_hv/N_hv

  #Entropy of line length distributions
  ENT_dl <- -1 * sum(P_dl * log(P_dl))
  ENT_vl <- -1 * sum(P_vl * log(P_vl))
  ENT_hl <- -1 * sum(P_hl * log(P_hl))
  ENT_hv <- -1 * sum(P_hv * log(P_hv))

  #Relative Entropy (Entropy / Max entropy)
  ENTrel_dl = ENT_dl/(-1 * log(1/DLmax))
  ENTrel_vl = ENT_vl/(-1 * log(1/VLmax))
  ENTrel_hl = ENT_hl/(-1 * log(1/HLmax))
  ENTrel_hv = ENT_hv/(-1 * log(1/max(c(HLmax,VLmax), na.rm = TRUE)))

  #Meanline
  MEAN_dl = mean(dlines, na.rm = TRUE)%00%0
  MEAN_vl = mean(vlines, na.rm = TRUE)%00%0
  MEAN_hl = mean(hlines, na.rm = TRUE)%00%0
  MEAN_hv = mean(c(hlines,vlines), na.rm = TRUE)%00%0

  #Maxline
  MAX_dl = max(freqvec_dl, na.rm = TRUE)%00%0
  MAX_vl = max(freqvec_vl, na.rm = TRUE)%00%0
  MAX_hl = max(freqvec_hl, na.rm = TRUE)%00%0
  MAX_hv = max(freqvec_hv, na.rm = TRUE)%00%0

  # REPetetiveness
  REP_av  <- (N_hlp+N_vlp) / N_dlp
  REP_hl  <-  N_hlp/N_dlp
  REP_vl  <-  N_vlp/N_dlp

  #Coefficient of determination
  CoV_dl = stats::sd(dlines)/mean(dlines)
  CoV_vl = stats::sd(vlines)/mean(vlines)

  CoV_hl = stats::sd(hlines)/mean(hlines)
  CoV_hv = stats::sd(c(hlines,vlines))/mean(c(hlines,vlines))

  #Divergence
  DIV_dl = 1/MAX_dl
  DIV_vl = 1/MAX_vl
  DIV_hl = 1/MAX_hl
  DIV_hv = 1/MAX_hv

  out <- data.frame(
    RP_N      = RP_N,
    DIV_dl    = DIV_dl,
    REP_av    = REP_av,
    N_dl      = N_dl,
    N_dlp     = N_dlp,
    DET       = DET,
    MEAN_dl   = MEAN_dl,
    MAX_dl    = MAX_dl,
    ENT_dl    = ENT_dl,
    ENTrel_dl = ENTrel_dl,
    CoV_dl    = CoV_dl,
    N_vl      = N_vl,
    N_vlp     = N_vlp,
    LAM_vl    = LAM_vl,
    TT_vl     = MEAN_vl,
    MAX_vl    = MAX_vl,
    ENT_vl    = ENT_vl,
    ENTrel_vl = ENTrel_vl,
    CoV_vl    = CoV_vl,
    REP_vl    = REP_vl,
    N_hlp     = N_hlp,
    N_hl      = N_hl,
    LAM_hl    = LAM_hl,
    TT_hl     = MEAN_hl,
    MAX_hl    = MAX_hl,
    ENT_hl    = ENT_hl,
    ENTrel_hl = ENTrel_hl,
    CoV_hl    = CoV_hl,
    REP_hl    = REP_hl,
    N_hvp     = N_hvp,
    N_hv      = N_hv,
    LAM_hv    = LAM_hv,
    TT_hv     = MEAN_hv,
    MAX_hv    = MAX_hv,
    ENT_hv    = ENT_hv,
    ENTrel_hv = ENTrel_hv,
    CoV_hv    = CoV_hv)

  if(matrices){
    return(list(
      crqaMeasures = out,
      crqaMatrices = list(RM     = RM,
                          dlines = dlines,
                          vlines = vlines,
                          hlines = hlines,
                          freq_dl = freq_dl,
                          freq_vl = freq_vl,
                          freq_hl = freq_hl))
    )
  } else {
    return(out)
  }
}


rp_calc_par <- function(X,
                        Y,
                        emRad = NULL,
                        DLmin = 2,
                        VLmin = 2,
                        HLmin = 2,
                        DLmax = NROW(X),
                        VLmax = NROW(X),
                        HLmax = NROW(X),
                        theiler   = NA,
                        AUTO      = NULL,
                        chromatic = FALSE,
                        anisotropyHV = FALSE,
                        asymmetryUL = FALSE,
                        returnUL = FALSE,
                        recurrenceTimes = FALSE,
                        matrices  = FALSE){


}


#' Prepare matrix
#'
#' @param RP Recurrence plot
#' @param emRad Radiuc
#' @param DLmin Minimal diagonal line length
#' @param VLmin Minimal vertical line length
#' @param HLmin Minimal horizontal line length
#' @param DLmax Maximal diagonal line length
#' @param VLmax Maximal vertical line length
#' @param HLmax Maximal horizontal line length
#' @param AUTO Is this an AUTO RQA?
#' @param chromatic Chromatic RQA?
#' @param matrices Return Matrices?
#' @param doHalf Analyse half of the matrix?
#'
#' @return A prepped matrix
#' @keywords internal
#'
#' @export
#'
rp_prep <- function(RP,
                    emRad = NULL,
                    DLmin = 2,
                    VLmin = 2,
                    HLmin = 2,
                    DLmax = length(Matrix::diag(RP)),
                    VLmax = length(Matrix::diag(RP)),
                    HLmax = length(Matrix::diag(RP)),
                    AUTO      = FALSE,
                    chromatic = FALSE,
                    matrices  = FALSE,
                    doHalf    = FALSE){

  out<-rp_calc(RP,
               emRad = emRad,
               DLmin = DLmin,
               VLmin = VLmin,
               HLmin = HLmin,
               DLmax = DLmax,
               VLmax = VLmax,
               HLmax = HLmax,
               AUTO  = AUTO,
               chromatic = chromatic,
               matrices  = matrices)

  if(doHalf){
    if(!AUTO){
      outLo <- rp_calc(Matrix::tril(RP,-1),
                       emRad = emRad,
                       DLmin = DLmin,
                       VLmin = VLmin,
                       HLmin = HLmin,
                       DLmax = DLmax,
                       VLmax = VLmax,
                       HLmax = HLmax,
                       AUTO  = AUTO,
                       chromatic = chromatic,
                       matrices  = matrices)

      outUp <- rp_calc(Matrix::triu(RP,-1),
                       emRad= emRad,
                       DLmin = DLmin,
                       VLmin = VLmin,
                       HLmin = HLmin,
                       DLmax = DLmax,
                       VLmax = VLmax,
                       HLmax = HLmax,
                       AUTO  = AUTO,
                       chromatic = chromatic,
                       matrices  = matrices)
      out <- cbind.data.frame(full  = out,
                              lower = outLo,
                              upper = outUp)
    } else {
      warning("doHalf = TRUE and AUTO = TRUE. Results would be the same for lower and upper triangle!")
      out<- cbind.data.frame(full  = out,
                             lower = rp_empty(),
                             upper = rp_empty())
    }
  }
  return(out)
}


#' False Nearest Neighbours
#'
#' Search for FNN to get an optimal Embedding Dimension using by using [nonlinearTseries::findAllNeighbours()] in a loop.
#'
#' @inheritParams est_parameters
#' @param radius Size of the neighbourhood: Every point smaller than the radius will be considered a near neighbour, see [nonlinearTseries::findAllNeighbours()] (default = `sd(y)/10`).
#' @param number.boxes Integer representing number of boxes to to speed up neighbour search, if `NULL` an optimal number will be chosen [nonlinearTseries::findAllNeighbours()] (default = `NULL`).
#'
#' @return FNN curve
#' @export
#'
fnn <- function(y, emLag = 1, maxDim = 10, radius = sd(y)/10, number.boxes = NULL){

  out <- matrix(NA, nrow=maxDim, ncol = 1)

  # for(emL in 1:length(emLag)){
  for(emD in 1: maxDim){

    yy <- nonlinearTseries::buildTakens(time.series = y, embedding.dim = emD, time.lag = emLag)
    nn <- nonlinearTseries::findAllNeighbours(as.matrix(yy[,1:emD]), radius = radius, number.boxes = number.boxes)
    out[emD,1] <- sum(lengths(nn))

  }
  # }
  return(out)
}