R/troppca.linsp2poly.R

In HoujieWang/Rtropical: Data Analysis Tools over Space of Phylogenetic Trees Using Tropical Geometry

#' Tropical Principal Component Analysis by Polytope Converted from Linear Space
#'
#' Approximate the principal component as a tropical polytope converted from tropical linear space for a given data matrix
#' via MCMC and return the results as an object of class \code{troppca}.
#'
#' @importFrom parallel parLapply
#' @importFrom parallel makeCluster
#' @importFrom parallel stopCluster
#' @importFrom Rfast rowMaxs
#' @importFrom Rfast colMins
#' @importFrom stats runif
#'
#' @param x a data matrix, of size n x e, with each row an observation vector.
#' e is the dimension of the tropical space
#' @param pcs a numeric value indicating the order of principal component. (default: 2)
#' @param iteration a list with arguments controlling the iteration of the algorithm.
#' \describe{
#' \item{exhaust}{a logical variable indicating if to iterate over all possible combinations of the linear space
#' based on the given data matrix \code{x}. If FALSE, please input a number of iteration for \code{niter}.
#' If TRUE, please enter 0 for \code{niter} and this function will iterate over all possible combinations of linear space.
#' This could be time consuming when \code{x} is large. (default: FALSE)}
#' \item{niter}{a numeric variable indicating the number of iterations. (default: 100)}
#' }
#' @param ncores a numeric value indicating the number of threads utilized for multi-cored CPUs. (default: 2)
#'
#' @return A list of S3 class \code{"troppca"}, including:
#' \item{pc}{The principal component as a tropical linear space}
#' \item{obj}{The tropical PCA objective, the sum of tropical distance from each point to the projection.}
#' \item{projection}{The projections of all data points.}
#' \item{type}{The geometry of principal component.}
#'
#'
#' @examples
#' \donttest{
#' library(Rfast)
#' n <- 50
#' e <- 50
#' s <- 5
#' x <- rbind(
#'   rmvnorm(n, mu = c(5, -5, rep(0, e - 2)), sigma = diag(s, e)),
#'   rmvnorm(n, mu = c(-5, 5, rep(0, e - 2)), sigma = diag(s, e))
#' )
#' troppca_fit <- troppca.linsp2poly(x)
#' }
#'
#' @export
#' @export troppca.linsp2poly
troppca.linsp2poly = function(x, pcs = 2, iteration = list(), ncores = 2){
  con <- list(
    exhaust = FALSE,
    niter = 100
  )
  con[names(iteration)] <- iteration

  exhaust <- con$exhaust
  niter <- con$niter
  pcs <- pcs + 1
  if (exhaust) {
    warning("Iterating all possible PC choices enabled, this could take long for high order PCs.")
    all_choices <- comboGeneral(nrow(x), pcs)
    all_choices <- lapply(1:nrow(all_choices), function(i) all_choices[i, ])
  } else {

    all_choices <- lapply(1: niter, function(i){sample(1: nrow(x), pcs, replace = F)})
  }
  x_list <- lapply(seq_len(nrow(x)), function(i) x[i, ])
  cl <- makeCluster(ncores)
  all_objs <- unlist(parLapply(cl, all_choices, function(ind){
    V <- x[ind, ]
    proj = do.call("rbind", lapply(x_list, tropproj.poly, t(linsp_to_poly(V))))
    temp <- x - proj
    sum(rowMaxs(temp, T) - rowMins(temp, T))
  }))
  stopCluster(cl)

  best_choice <- all_choices[[which.min(all_objs)]]
  pc <- x[best_choice, ]
  rownames(pc) <- paste("pc", 1: pcs, sep = "")
  proj_points <- do.call("rbind", lapply(x_list, tropproj.poly, t(linsp_to_poly(pc))))
  troppca.out <- list("pc" = pc,
                     "obj" = min(all_objs),
                     "projection" = proj_points,
                     "type" = "linear space")
  class(troppca.out) <- "troppca"
  troppca.out
}
# troppca.linsp2poly <- function(x, pcs = 2, nsample = 1000, ncores = 2) {
#   pcs <- pcs + 1
#   n <- nrow(x)
#   cl <- makeCluster(ncores)
#   x_list <- lapply(seq_len(n), function(i) x[i, ])
#   troppca_objs <- vector(mode = "numeric", nsample)
#   samples <- matrix(NA, nrow = nsample, ncol = pcs)
#   samples[1, ] <- sample(1:n, pcs)
#   troppca_objs[1] <- troppca.obj2(x[samples[1, ], ], x_list, cl)
#
#   t <- 1
#   while (t < nsample) {
#     # Find a new proposal by changing a randomly selected vertex of the current polytope
#     current_choice <- samples[t, ]
#     current_obj <- troppca_objs[t]
#
#     change_ind <- sample(pcs, 1)
#     out_change <- sample(c(1:n)[-current_choice], 1)
#     new_choice <- c(current_choice[-change_ind], out_change)
#     new_obj <- troppca.obj2(x[new_choice, ], x_list, cl)
#
#     # Compute the probability we accept the new PCA base
#     p <- min(1, current_obj / new_obj)
#
#     if (sample(c(0, 1), 1, prob = c(1 - p, p)) == 1) {
#       samples[(t + 1), ] <- new_choice
#       troppca_objs[(t + 1)] <- new_obj
#       t <- t + 1
#     }
#   }
#   min_index <- which(troppca_objs == min(troppca_objs))[1]
#   best_obj <- troppca_objs[min_index]
#   pc <- x[samples[min_index, ], ]
#   pc2 <- linsp_to_poly(pc)
#   proj_points <- do.call("rbind", parLapply(cl, x_list, tropproj.poly, tconv = t(pc2)))
#   stopCluster(cl)
#   rownames(pc) <- paste("pc", 1:pcs, sep = "")
#   troppca.out <- list(
#     "pc" = pc,
#     "obj" = troppca_objs[min_index],
#     "projection" = proj_points,
#     "type" = "linear space"
#   )
#   class(troppca.out) <- "troppca"
#   troppca.out
# }