#' ROPES algorithm
#' 
#' @aliases ropes
#'
#' @description 
#' This function runs the ROPES algorithm.
#' It is used inside \code{\link{executing_ropes}}.
#' 
#' @usage 
#' ropes(b, w, lambda2Ages = 1, lambda1Ages = 0, lambda0Ages = 0.01, lambda0Inds = 0.1, 
#'       trace = 0, UStart = NULL, VStart = NULL, method = "cg", maxIter = 100000, 
#'       smoothMethod = splineMatrix)
#'  
#' @param b Matrix b from \code{\link{b_w}}.  
#' @param w Matrix w from \code{\link{b_w}}.
#' @param lambda2Ages Tuning parameter (see page 22 of the paper 
#' referenced below).
#' @param lambda1Ages Tuning parameter (page 22).
#' @param lambda0Ages Tuning parameter (page 22.)
#' @param lambda0Inds Additional tuning parameter.
#' @param trace Tracing information about the progress of the
#' optimization, see \code{\link{optim}}.
#' @param UStart Initial value for matrix U.
#' @param VStart Initial value for matrix V.
#' @param method Optimization method to be used. Options 'cg' (default) 
#' or 'bfgs'.
#' @param maxIter Maximum number of iterations.
#' @param smoothMethod Method to compute the average of all curves.         
#' 
#' @return 
#' A list with the following elements:
#' \emph{Z}: Matrix of the same dimension as \code{b}, close to \code{b} at
#' the points where observations are available. Z=UV^T
#' \emph{U}: Matrix with the 'scores'.
#' \emph{V}: Matrix with the 'features'.
#' \emph{Avg}: Average of all curves.
#' \emph{Canonical}: Canonical ROPES solution.
#' 
#' @note 
#' 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{executing_ropes}}
#' 
#' @export

# If we don't remove the observations with only one measurement, we get these warnings:
#In predict.lm(fit, newdata = data.frame(nNaInds = seq_along(v))) :
#  prediction from a rank-deficient fit may be misleading
options(warn = -1)

ropes <- function(b, w, lambda2Ages = 1, lambda1Ages = 0, lambda0Ages = 0.01, 
                  lambda0Inds = 0.1, trace = 0, UStart = NULL, VStart = NULL, 
                  method = "cg", maxIter = 100000, smoothMethod = splineMatrix
                  ){
  Y = b
  W = w
  
  #print(lambda2Ages)
  #print(lambda1Ages)
  #print(lambda0Ages)
  
  b[w == 0] = NA
  avgs = smoothMethod(b)
  for (j in 1:ncol(Y)) {
    Y[,j] = Y[,j] - avgs[1,j] # We prepare the data by subtracting the 
                              # average of all curves (22 of 32 in the paper).
  }
  Y[is.na(Y)] = 0
  
  if (lambda0Inds <= 0) A0 = NULL else
    A0 = weightedNDiffMatrix(nrow(Y), 0, lambda0Inds)
  
  if (lambda0Ages <= 0) B0 = NULL else
    B0 = weightedNDiffMatrix(ncol(Y), 0, lambda0Ages)
  
  if (lambda1Ages <= 0) B1 = NULL else
    B1 = weightedNDiffMatrix(ncol(Y), 1, lambda1Ages)
  
  if (lambda2Ages <= 0) B2 = NULL else
    B2 = weightedNDiffMatrix(ncol(Y), 2, lambda2Ages)
  
  A = A0
  B = rbind(B2, B1, B0)
  
  if (is.null(UStart) || is.null(VStart)) {
    for (i in 1:nrow(b)) {
      #print(i)
      v = b[i,]
      nNaInds = seq_along(v)[!is.na(v)]
      nNaVals = v[!is.na(v)]
      fit = lm(nNaVals ~ nNaInds)
      #print(fit)
      #if(is.na(fit$coefficients[2])){ # This would be needed if we don't remove players with 
      # only one measurement.
      #  #fit$coefficients[2] <- 0 # This still gives an error in predict
      #  fcast <- rep(fit$coefficients[1], ncol(b))
      #  names(fcast) <- 1:ncol(b)
      #}
      fcast = predict(fit, newdata = data.frame(nNaInds = seq_along(v)))
      b[i,] = fcast
    }
    
    #Calculation of the starting values for U and V
    svdY = svd(b)
    sqrtD = diag(sqrt(svdY$d))
    UStart = svdY$u %*% sqrtD
    VStart = svdY$v %*% sqrtD
  }
  
  if (method == "bfgs") {
    fcast = bfgs(W = W, Y = Y, A = A, B = B, UStart = UStart, VStart = VStart, 
                 maxIter = maxIter, factr = 10, trace = trace, lmm = 10)
  }
  else if (method == "cg") {
    fcast = cg(W = W, Y = Y, A = A, B = B, UStart = UStart, VStart = VStart, 
               maxIter = maxIter, type = 3, reltol = 1E-8, trace = trace)
  }
  else {
    stop(paste("\nMethod name '", method, "' is incorrect"))
  }
  
  result = fcast$U %*% t(fcast$V)
  
  #Canonical
  zStar = A %*% result %*% t(B)
  udvStar = svd(zStar)
  canonical = list(P = ginv(A) %*% udvStar$u, D = udvStar$d, 
                   Q = ginv(B) %*% udvStar$v)
  
  for (j in 1:ncol(result)) {
    result[,j] = result[,j] + avgs[1,j]
  }
  return(list(Z = result, U = fcast$U, V = fcast$V, 
              Avg = avgs, Canonical = canonical))
}