#' Prediction intervals
#' 
#' @aliases forec_intervals_ropes
#'
#' @description 
#' This function returns the predictions intervals with the ROPES
#' algorithm, see Section 5 of the paper referenced below.
#' 
#' @usage 
#' forec_intervals_ropes(res_ropes, nMCIter, ...)
#' 
#' @param res_ropes List returned by \code{\link{executing_ropes}}.
#' @param nMCIter Number of Monte-Carlo iterations.
#' @param ... Further parameters for \code{\link{ropes}}.
#' 
#' @return 
#' Matrix with the prediction intervals.
#' 
#' @note 
#' The code of this function was kindly provided by 
#' Alexander Dokumentov and Rob Hyndman.
#'
#' @references 
#' Dokumentov, A., Hyndman, R. J., 2016. Low-dimensional decomposition, 
#' smoothing and forecasting of sparse functional data, 
#' \url{http://robjhyndman.com/papers/ROPES.pdf}. Working paper.
#'  
#' @author 
#' Alexander Dokumentov and Rob Hyndman
#' 
#' @seealso 
#' \code{\link{ropes}}, \code{\link{executing_ropes}}
#' 
#' @examples 
#' # See the code.
#'
#' @importFrom parallel detectCores
#' @importFrom foreach foreach
#' @importFrom doMC registerDoMC
#' @importFrom bigmemory big.matrix describe attach.big.matrix
#' @importFrom synchronicity boost.mutex attach.mutex lock unlock
#'
#' @export

forec_intervals_ropes <- function(res_ropes, nMCIter, ...){
  z <- res_ropes$result$Z
  z1 <- z # This is to define m = array(0.0, c(nMCIter, dim(z))) 

  b <- res_ropes$b_12
  w <- res_ropes$w_12
  resid <- b - z
  stdev <- sd(resid, na.rm = TRUE)

  nCores <- detectCores()
  registerDoMC(nCores)  # Number of CPU cores

  # CALCULATION OF THE PREDICTION INTERVALS:
  file.remove("log.txt")
  nMCIter <- nMCIter
  counter <- big.matrix(1, 1, type = 'double', init = nMCIter)
  counterDescr <- describe(counter)
  mutex = boost.mutex()
  mutexDescr = describe(mutex)
  result <- foreach(j = 1:nCores) %dopar%
  {
    l = list()
    cr <- attach.big.matrix(counterDescr)
    mx <- attach.mutex(mutexDescr)
    for (i in 1:nMCIter) {
      lock(mx)
      doWork = cr[1,1]
      if (doWork > 0) {
        cr[1,1] <- doWork - 1
      }
      unlock(mx)
      if (doWork <= 0) break
      
      sink("log.txt", append = TRUE)
      cat(paste("\n", Sys.time(), " -> Thread ", j,", Iteration ", i, ", Work ", doWork, sep = ""))
      
      db = b # This is to define the matrix \Delta in step 1. for 
             # prediction intervals (page 21 of 32 of the reference above).
      fnNa = !is.na(db)       
      nnNa = sum(!is.na(db)) 
      
      db[fnNa] = db[fnNa] + rnorm(nnNa, mean = 0, sd = stdev) # This is the set of "distorted" solutions. 
      # Y_{\Delta} = Y + \Delta 
      # (rnorm is because both the matrices \Delta_k and \Upsilon_k 
      # have elements i.i.d. N(0,\sigma^2) ) and we find a set of 
      # solutions Z(Y_{\Delta}) with the bone function:
      dResult = ropes(db, w, ...)
      
      l[[i]] <- list(thread = j, interation = i, z = dResult$Z)
    }
    l
  }

  # THIS CODE GENERATES A log.txt FILE
  # THE FOLLOWING WARNING APPEARS THE FIRST TIME I RUN THE CODE BECAUSE 
  # I DIDN'T HAVE ANY log.txt FILE
  #Warning message:
  #  In file.remove("log.txt") :
  #  cannot remove file 'log.txt', reason 'No such file or directory'
  
  m <- array(0.0, c(nMCIter, dim(z1))) 
  k <- 1
  for (i in 1:length(result)) {
    l = result[[i]]
    for (j in 1:length(l)) {
      m[k,,] = l[[j]]$z # each m is the element Z_{ij}, 
      # specifically this is Z(Y_{\Delta}) in page 21 of 32 for 
      # prediction intervals.
      k = k + 1
    }
  }
  
  nNorm <- 20
  probs <- c(0.025, 0.2, 0.5, 0.8, 0.975)

  q <- array(0.0, c(length(probs), dim(z)))
  dimnames(q) <- list(as.character(probs), rownames(b), colnames(b))

  for (k in 1:(dim(m))[2]) {    
    for (l in 1:(dim(m))[3]) {   
      v = c()
      for (r in 1:dim(m)[1]) {   
        v = c(v, rnorm(nNorm, mean = m[r,k,l], sd = stdev)) # stdev = sd(resid, na.rm=T)
        # Specifically, this is the set of "distorted" solutions Z(Y_{\Delta})
        # k and l correspond to the subindices i and j, respectively.
        # Remember z = result$Z
      }
      q[,k,l] = quantile(v, probs = probs)
    }
  }
  # q are the prediction intervals.
 
  return(q)
}  