R/RotMat.R

In ODRF: Oblique Decision Random Forest for Classification and Regression

#' Random Rotation Matrix
#'
#' Generate rotation matrices by different distributions, and it comes from the library \code{rerf}.
#'
#' @param dimX The number of dimensions.
#' @param randDist The probability distribution of the random projection direction, including "Binary": \eqn{B\{-1,1\}} binomial distribution (default),
#' "Norm":\eqn{N(0,1)} normal distribution, "Uniform": \eqn{U(-1,1)} uniform distribution.
#' @param numProj The number of projection directions (default ceiling(sqrt(\code{dimX}))).
#' @param dimProj Number of variables to be projected, default dimProj="Rand": random from 1 to \code{dimX}.
#' @param sparsity A real number in \eqn{(0,1)} that specifies the distribution of non-zero elements in the random matrix.
#' When sparsity="pois" means that non-zero elements are generated by the p(\code{lambda}) Poisson distribution.
#' @param prob A probability in \eqn{(0,1)} used for sampling from \eqn{{-1,1}} where \code{prob = 0} will only sample -1 and \code{prob = 1} will only sample 1.
#' @param lambda Parameter of the Poisson distribution (default 1).
#' @param catLabel A category labels of class \code{list} in predictors. (default NULL, for details see Examples of \code{\link{ODT}})
#' @param ... Used to handle superfluous arguments passed in using paramList.
#'
#' @return A random matrix to use in running \code{\link{ODT}}.
#' \itemize{
#' \item{Variable: Variables to be projected.}
#' \item{Number: Number of projections.}
#' \item{Coefficient: Coefficients of the projection matrix.}
#' }
#'
#' @references Tomita, T. M., Browne, J., Shen, C., Chung, J., Patsolic, J. L., Falk, B., ... & Vogelstein, J. T. (2020). Sparse projection oblique randomer forests. Journal of machine learning research, 21(104).
#' @keywords rotation
#'
#' @seealso \code{\link{RotMatPPO}} \code{\link{RotMatRF}} \code{\link{RotMatMake}}
#'
#' @examples
#' set.seed(1)
#' paramList <- list(dimX = 8, numProj = 3, sparsity = 0.25, prob = 0.5)
#' (RotMat <- do.call(RotMatRand, paramList))
#' paramList <- list(dimX = 8, numProj = 3, sparsity = "pois")
#' (RotMat <- do.call(RotMatRand, paramList))
#' paramList <- list(dimX = 8, randDist = "Norm", dimProj = 5)
#' (RotMat <- do.call(RotMatRand, paramList))
#'
#' @export
RotMatRand <- function(dimX, randDist = "Binary", numProj = ceiling(sqrt(dimX)),
                       dimProj = "Rand", sparsity = ifelse(dimX >= 10, 3 / dimX, 1 / dimX),
                       prob = 0.5, lambda = 1, catLabel = NULL, ...) {
  dimX <- dimX - (!is.null(catLabel)) * (length(unlist(catLabel)) - length(catLabel))

  if (dimProj != "Rand") {
    if (dimProj > dimX) {
      stop("ERROR: parameter dimProj is greater than the number of dimensions dimX.")
    }

    nnz <- dimProj * numProj
    nz.rows <- integer(nnz)
    nz.cols <- integer(nnz)
    start.idx <- 1L
    for (i in seq.int(numProj)) {
      end.idx <- start.idx + dimProj - 1L
      nz.rows[start.idx:end.idx] <- sample.int(dimProj, dimProj, replace = FALSE)
      nz.cols[start.idx:end.idx] <- i
      start.idx <- end.idx + 1L
    }
  } else if (sparsity == "pois") {
    nnzPerCol <- stats::rpois(numProj, lambda)
    while (!any(nnzPerCol)) {
      nnzPerCol <- stats::rpois(numProj, lambda)
    }

    nnzPerCol[nnzPerCol > dimX] <- dimX
    nnz <- sum(nnzPerCol)
    nz.rows <- integer(nnz)
    nz.cols <- integer(nnz)
    start.idx <- 1L
    for (i in seq.int(numProj)) {
      if (nnzPerCol[i] != 0L) {
        end.idx <- start.idx + nnzPerCol[i] - 1L
        nz.rows[start.idx:end.idx] <- sample.int(dimX, nnzPerCol[i], replace = FALSE)
        nz.cols[start.idx:end.idx] <- i
        start.idx <- end.idx + 1L
      }
    }
  } else {
    nnzs <- round(dimX * numProj * sparsity)
    ind <- sort(sample.int((dimX * numProj), nnzs, replace = FALSE))
    nz.rows <- ((ind - 1L) %% dimX) + 1L
    nz.cols <- floor((ind - 1L) / dimX) + 1L
  }


  if (randDist == "Binary") {
    proj <- sample(c(-1L, 1L), length(nz.rows), replace = TRUE, prob = c(
      prob, 1 - prob
    ))
  }
  if (randDist == "Norm") {
    proj <- rnorm(length(nz.rows))
  }
  if (randDist == "Uniform") {
    proj <- runif(length(nz.rows), -1, 1)
  }

  randomMatrix <- cbind(nz.rows, nz.cols, proj)

  ## Determine if categorical variables need to be taken d
  ## into consideration
  if (!is.null(catLabel)) {
    ind <- randomMatrix[, 1]
    catVar <- which(ind <= length(catLabel))
    # catVar=which(ind<=length(catLabel))
    randomMatrix[-catVar, 1L] <- randomMatrix[-catVar, 1L] + length(unlist(catLabel)) - length(catLabel)

    catMap <- 1
    for (xj in unique(ind[catVar])) {
      isj <- which(ind == xj)
      randomMatrix[isj, 1L] <- sample(catMap:(catMap + length(catLabel[[xj]]) - 1), length(isj), replace = TRUE)
      catMap <- catMap + length(catLabel[[xj]])
    }
  }

  colnames(randomMatrix) <- c("Variable", "Number", "Coefficient")
  return(randomMatrix)
}


#' Create a Projection Matrix: Random Forest (RF)
#'
#' Create a projection matrix with coefficient 1 and 0 such that the ODRF (ODT) has the same partition variables as the Random Forest (CART).
#'
#' @param dimX The number of dimensions.
#' @param numProj The number of projection directions (default ceiling(sqrt(\code{dimX}))).
#' @param catLabel A category labels of class \code{list} in predictors. (default NULL, for details see Examples of \code{\link{ODT}})
#' @param ... Used to handle superfluous arguments passed in using paramList.
#'
#' @return A random matrix to use in running \code{\link{ODT}}.
#' \itemize{
#' \item{Variable: Variables to be projected.}
#' \item{Number: Number of projections.}
#' \item{Coefficient: Coefficients of the projection matrix.}
#' }
#'
#' @keywords rotation
#'
#' @seealso \code{\link{RotMatPPO}} \code{\link{RotMatRand}} \code{\link{RotMatMake}}
#'
#' @examples
#' paramList <- list(dimX = 8, numProj = 3, catLabel = NULL)
#' set.seed(2)
#' (RotMat <- do.call(RotMatRF, paramList))
#'
#' @export
RotMatRF <- function(dimX, numProj, catLabel = NULL, ...) {
  dimX <- dimX - (!is.null(catLabel)) * (length(unlist(catLabel)) - length(catLabel))
  if (numProj > dimX) {
    stop("ERROR: parameter numProj is greater than the number of dimensions dimX.")
  }
  randomMatrix <- cbind(sample.int(dimX, numProj, replace = FALSE), 1:numProj, rep(1L, numProj))
  if (!is.null(catLabel)) {
    ind <- randomMatrix[, 1]
    catVar <- which(ind <= length(catLabel))
    # catVar=which(ind<=length(catLabel))
    randomMatrix[-catVar, 1L] <- randomMatrix[-catVar, 1L] + length(unlist(catLabel)) - length(catLabel)

    catMap <- 1
    for (xj in unique(ind[catVar])) {
      isj <- which(ind == xj)
      randomMatrix[isj, 1L] <- sample(catMap:(catMap + length(catLabel[[xj]]) - 1), length(isj), replace = TRUE)
      catMap <- catMap + length(catLabel[[xj]])
    }
  }

  colnames(randomMatrix) <- c("Variable", "Number", "Coefficient")
  return(randomMatrix)
}


#################################################################################
#' Create a Projection Matrix: RotMatPPO
#'
#' Create a projection matrix using projection pursuit optimization (\code{\link{PPO}}).
#'
#' @param X An n by d numeric matrix (preferable) or data frame.
#' @param y A response vector of length n.
#' @param model Model for projection pursuit (for details see \code{\link{PPO}}).
#' @param split One of three criteria, 'gini': gini impurity index (classification), 'entropy': information gain (classification, default)
#' or 'mse': mean square error (regression).
#' @param weights A vector of length same as \code{data} that are positive weights. (default NULL)
#' @param dimProj Number of variables to be projected, \code{dimProj}=min(ceiling(n^0.4),ceiling(ncol(X)*2/3)) (default) or dimProj="Rand": random from 1 to ncol(X).
#' @param numProj The number of projection directions, when \code{dimProj}="Rand" default
#' \code{numProj} = sample(ceiling(ncol(\code{X})/3),1) otherwise default \code{numProj}=ceiling(ncol(\code{X})/\code{dimProj}).
#' @param catLabel A category labels of class \code{list} in predictors. (default NULL, for details see Examples of \code{\link{ODT}})
#' @param ... Used to handle superfluous arguments passed in using paramList.
#'
#' @return A random matrix to use in running \code{\link{ODT}}.
#' \itemize{
#' \item{Variable: Variables to be projected.}
#' \item{Number: Number of projections.}
#' \item{Coefficient: Coefficients of the projection matrix.}
#' }
#'
#' @keywords rotation
#'
#' @seealso \code{\link{RotMatMake}} \code{\link{RotMatRand}} \code{\link{RotMatRF}} \code{\link{PPO}}
#'
#' @examples
#' set.seed(220828)
#' X <- matrix(rnorm(1000), 100, 10)
#' y <- (rnorm(100) > 0) + 0
#' (RotMat <- RotMatPPO(X, y))
#' (RotMat <- RotMatPPO(X, y, dimProj = "Rand"))
#' (RotMat <- RotMatPPO(X, y, dimProj = 6, numProj = 4))
#'
#' # classification
#' data(seeds)
#' (PP <- RotMatPPO(seeds[, 1:7], seeds[, 8], model = "Log", split = "entropy"))
#' (PP <- RotMatPPO(seeds[, 1:7], seeds[, 8], model = "PPR", split = "entropy"))
#' (PP <- RotMatPPO(seeds[, 1:7], seeds[, 8], model = "LDA", split = "entropy"))
#'
#' # regression
#' data(body_fat)
#' (PP <- RotMatPPO(body_fat[, 2:15], body_fat[, 1], model = "Log", split = "mse"))
#' (PP <- RotMatPPO(body_fat[, 2:15], body_fat[, 1], model = "Rand", split = "mse"))
#' (PP <- RotMatPPO(body_fat[, 2:15], body_fat[, 1], model = "PPR", split = "mse"))
#'
#' @importFrom stats rnorm lm ppr
#' @importFrom nnet nnet
#' @export
RotMatPPO <- function(X, y, model = "PPR", split = "entropy", weights = NULL, dimProj = min(ceiling(length(y)^0.4), ceiling(ncol(X) * 2 / 3)),
                      numProj = ifelse(dimProj == "Rand", sample(floor(ncol(X) / 3), 1), ceiling(ncol(X) / dimProj)), catLabel = NULL, ...) {
  # numProj = ifelse(dimProj == "Rand", max(5, sample(floor(ncol(X) / 3), 1)), max(5, ceiling(ncol(X) / dimProj)))
  if (dimProj != "Rand") {
    if (dimProj > ncol(X)) {
      stop("ERROR: parameter dimProj is greater than the number of dimensions.")
    }
  }
  # cancor glm
  X <- as.matrix(X)
  p <- ncol(X)
  n <- length(y)
  weights <- if (is.null(weights)) rep(1, n)
  # d=min(100, max(5, ceiling(p/q))) d q
  p0 <- p - (!is.null(catLabel)) * (length(unlist(catLabel)) - length(catLabel))

  if ((numProj == 1) && (p0 < 10)) {
    sparseM <- NULL
    d2 <- 0
  } else {
    d2 <- ceiling(sqrt(p0))
    ind <- sample(p0, d2, replace = FALSE)
    sparseM <- cbind(ind, 1:d2, rep(1, d2))
  }

  # ceiling(sqrt(p))#d#round(p/3)#
  # if((n>Q)&(n>10)&(p>1)){
  sparseM1 <- NULL
  if (n > 10) {
    numProj <- min(p0, numProj)
    d1 <- numProj
    indp <- p0
    if (d1 == 1) {
      sparseM1 <- cbind(seq(indp), 1, 1)
    } else {
      if (dimProj == "Rand") {
        # if(is.null(numProj)){max(5,sample(floor(p0/3),1))}
        spd <- sample(p0, d1)
        indp <- sum(spd)
        ind <- unlist(sapply(spd, function(pd) sample(p0, pd)))
        sparseM1 <- cbind(ind, d2 + rep(1:d1, spd), rep(1, indp))
      } else {
        # if(is.null(dimProj)){dimProj =min(ceiling(n^0.4),ceiling(p0*2/3))}
        # if(is.null(numProj)){numProj=max(5, ceiling(p0/dimProj))}
        # indp <- p0
        # sparseM=NULL
        # for (k in 1:ceiling(2*q*d/p)) {
        ind <- sample(1:indp, replace = FALSE)
        s <- c(0, floor(quantile(1:indp, (1:d1) / d1)))
        sparseM1 <- cbind(ind, d2 + rep(1:d1, s[-1] - s[-d1 - 1]), rep(1, indp))
      }
    }
  }

  sparseM <- rbind(sparseM, sparseM1)
  if (is.null(sparseM)) sparseM <- t(c(sample(p0, 1), 1, 1))

  if (!is.null(catLabel)) {
    ind <- sparseM[, 1]
    catVar <- which(ind <= length(catLabel))
    sparseM[-catVar, 1L] <- sparseM[-catVar, 1L] + length(unlist(catLabel)) - length(catLabel)

    catMap <- 1
    for (xj in unique(ind[catVar])) {
      isj <- which(ind == xj)
      sparseM[isj, 1L] <- sample(catMap:(catMap + length(catLabel[[xj]]) - 1), length(isj), replace = TRUE)
      catMap <- catMap + length(catLabel[[xj]])
    }
  }

  ##########################
  if (n > 10 && nrow(sparseM1) > 1) {
    Yi <- c(y)
    indC <- 0L
    if (split != "mse") {
      y <- as.factor(y)
      indC <- levels(y)
      if (length(indC) > 2) {
        Yi <- (matrix(y, n, length(indC)) == matrix(indC, n, length(indC), byrow = TRUE)) + 0
      } else {
        Yi <- as.integer(y)
      }
    }

    sparseM1 <- sparseM[d2 + (1:indp), ]
    sparseM <- sparseM[-(d2 + (1:indp)), ]
    indTab <- table(sparseM1[, 2])
    ind <- as.numeric(names(indTab)[which(indTab > 1)])
    jx <- which(sparseM1[, 2] %in% ind)
    d11 <- min(max(c(indTab, d1)), length(unique(sparseM1[jx, 1])))
    # d11=min(c(length(unique(sparseM[jx,1])),ifelse(split=='r',Inf,Inf)))#
    # d11=max(max(indTab),min(length(unique(sparseM1[jx,1])),d1))
    if (length(unique(sparseM1[jx, 2])) == 1) d11 <- 0

    sparseM1 <- rbind(sparseM1, matrix(d1 + d2 + 1, d11, 3))
    for (ni in unique(sparseM1[, 2L])) {
      lrows <- which(sparseM1[, 2L] == ni)
      j.I <- sparseM1[lrows, 1L]
      if (ni == (d1 + d2 + 1)) {
        # jx=which(sparseM1[, 3L]!=1L)
        ix <- order(abs(sparseM1[jx, 3L]), decreasing = TRUE)
        j.I <- unique(sparseM1[jx[ix], 1])[1:d11]
        lrows <- indp + (1:d11)
        sparseM1[lrows, 1L] <- j.I
      }
      Xi <- X[, j.I, drop = FALSE]
      pi <- length(j.I)

      # S = eigen(cov(Xi))                # using PCA to handle the colinearity
      # Ii = 1:min(which((cumsum(S$values)+1e-4)/sum(S$values,1e-4) > 0.99))
      # Xi = Xi%*%S$vectors[,Ii, drop = FALSE]
      if (pi > 1L) {
        if (model == "PPR") {
          PP <- try(ppr(Xi, Yi, weights = weights, nterms = 1, bass = 1)$alpha, silent = TRUE) # sm.method="spline",sm.method="gcvspline"
        } else if (model == "Log") {
          if ((n > 5 * pi) && (pi < 20)) {
            PP <- try(ppr(Xi, Yi, weights = weights, nterms = 1, bass = 1)$alpha, silent = TRUE)
          } else {
            # PP = try(nnet(Xi, y, size=1,trace=FALSE)$wts[2:(1+pi)], silent = TRUE)
            # PP = nnet(y~Xi,weights=weights, size=1,linout=TRUE,trace=FALSE)$wts[2:(1+pi)]#
            PP <- nnet(Xi, Yi, weights = weights, size = 1, linout = TRUE, trace = FALSE)$wts[2:(1 + pi)] #
          }
        } else {
          PP <- PPO(Xi, y, model, split)
        }

        if (inherits(PP, "try-error")) {
          LM <- lm(Yi ~ ., data = data.frame(Xi), weights)
          if (length(indC) > 2) {
            theta.i <- as.matrix(LM$coefficients)[-1, , drop = FALSE]
            for (j in seq(ncol(theta.i))) {
              theta.i[is.na(theta.i[, j]), j] <- 0.0
            }
            theta.i <- eigen(theta.i %*% t(theta.i))$vectors[, 1]
          } else {
            theta.i <- LM$coefficients[-1]
            theta.i[is.na(theta.i)] <- 0.0
          }
        } else {
          theta.i <- PP
        }

        # theta.i=S$vectors[,Ii]%*%theta.i
      } else {
        theta.i <- rep(1, pi)
      }

      theta.i <- theta.i / sqrt(sum(theta.i^2))
      sparseM1[lrows, 3L] <- c(theta.i)
    }

    # sparseM1[,2]=sparseM1[,2]+d2
    # sparseM[d2+(1:p0),]=sparseM1
    sparseM <- rbind(sparseM, sparseM1)
  }

  sparseM[is.na(sparseM[, 3]), 3] <- 1

  colnames(sparseM) <- c("Variable", "Number", "Coefficient")
  return(sparseM)
}

##################################################################################
#' Create rotation matrix used to determine the linear combination of features.
#'
#' Create any projection matrix with a self-defined projection matrix function and projection optimization model function
#'
#' @param X An n by d numeric matrix (preferable) or data frame.
#' @param y A response vector of length n.
#' @param RotMatFun A self-defined projection matrix function name, which can also be \code{\link{RotMatRand}} and \code{\link{RotMatPPO}}. Note that \code{(,...)} is necessary.
#' @param PPFun A self-defined projection function name, which can also be \code{\link{PPO}}. Note that \code{(,...)} is necessary.
#' @param FunDir The path to the \code{function} of the user-defined \code{NodeRotateFun} (default current Workspace).
#' @param paramList List of parameters used by the functions \code{RotMatFun} and \code{PPFun}. If left unchanged, default values will be used, for details see \code{\link[ODRF]{defaults}}.
#' @param ... Used to handle superfluous arguments passed in using paramList.
#'
#' @details There are two ways for the user to define a projection direction function. The first way is to connect two custom functions with the function \code{RotMatMake()}.
#' Specifically, \code{RotMatFun()} is defined to determine the variables to be projected, the projection dimensions and the number of projections (the first two columns of the rotation matrix).
#' \code{PPFun()} is defined to determine the projection coefficients (the third column of the rotation matrix). After that let the argument \code{RotMatFun="RotMatMake"},
#' and the argument \code{paramList} must contain the parameters \code{RotMatFun} and \code{PPFun}. The second way is to define a function directly,
#' and just let the argument \code{RotMatFun} be the name of the defined function and let the argument \code{paramList} be the arguments list used in the defined function.
#'
#' @return A random matrix to use in running \code{\link{ODT}}.
#' \itemize{
#' \item{Variable: Variables to be projected.}
#' \item{Number: Number of projections.}
#' \item{Coefficient: Coefficients of the projection matrix.}
#' }
#'
#' @seealso \code{\link{RotMatPPO}} \code{\link{RotMatRand}} \code{\link{RotMatRF}}
#'
#' @examples
#' set.seed(220828)
#' X <- matrix(rnorm(1000), 100, 10)
#' y <- (rnorm(100) > 0) + 0
#' (RotMat <- RotMatMake(X, y, "RotMatRand", "PPO"))
#' library(nnet)
#' (RotMat <- RotMatMake(X, y, "RotMatPPO", "PPO", paramList = list(model = "Log")))
#'
#' ## Define projection matrix function makeRotMat and projection pursuit function makePP.##
#' ##  Note that '...' is necessary.
#' makeRotMat <- function(dimX, dimProj, numProj, ...) {
#'   RotMat <- matrix(1, dimProj * numProj, 3)
#'   for (np in seq(numProj)) {
#'     RotMat[(dimProj * (np - 1) + 1):(dimProj * np), 1] <-
#'       sample(1:dimX, dimProj, replace = FALSE)
#'     RotMat[(dimProj * (np - 1) + 1):(dimProj * np), 2] <- np
#'   }
#'   return(RotMat)
#' }
#'
#' makePP <- function(dimProj, prob, ...) {
#'   pp <- sample(c(1L, -1L), dimProj, replace = TRUE, prob = c(prob, 1 - prob))
#'   return(pp)
#' }
#' \donttest{
#' RotMat <- RotMatMake(
#'   RotMatFun = "makeRotMat", PPFun = "makePP",
#'   paramList = list(dimX = 8, dimProj = 5, numProj = 4, prob = 0.5)
#' )
#' head(RotMat)
#' #>      Variable Number Coefficient
#' #> [1,]        6      1           1
#' #> [2,]        8      1           1
#' #> [3,]        1      1          -1
#' #> [4,]        4      1          -1
#' #> [5,]        5      1          -1
#' #> [6,]        6      2           1
#' }
#' \donttest{
#' # train ODT with defined projection matrix function
#' tree <- ODT(X, y,
#'   split = "entropy", NodeRotateFun = "makeRotMat",
#'   paramList = list(dimX = ncol(X), dimProj = 5, numProj = 4)
#' )
#' # train ODT with defined projection matrix function and projection optimization model function
#' tree <- ODT(X, y,
#'   split = "entropy", NodeRotateFun = "RotMatMake", paramList =
#'     list(
#'       RotMatFun = "makeRotMat", PPFun = "makePP",
#'       dimX = ncol(X), dimProj = 5, numProj = 4, prob = 0.5
#'     )
#' )
#' }
#' @keywords rotation
#' @export
RotMatMake <- function(X = NULL, y = NULL, RotMatFun = "RotMatPPO", PPFun = "PPO", FunDir = getwd(), paramList = NULL, ...) {
  p <- 1
  if (!is.null(X)) {
    X <- as.matrix(X)
    p <- ncol(X)
  } else if (!is.null(paramList$dimX)) {
    p <- paramList$dimX
  }

  if (RotMatFun == "RotMatPPO") {
    if (is.null(paramList$dimProj)) {
      paramList$dimProj <- min(ceiling(length(y)^0.4), ceiling(p * 2 / 3))
    }
    if (is.null(paramList$numProj)) {
      paramList$numProj <- ifelse(paramList$dimProj == "Rand", sample(floor(p / 3), 1), ceiling(p / paramList$dimProj))
    }
  }

  paramList <- defaults(paramList, dimX = p)

  if ((!RotMatFun %in% ls("package:ODRF", pattern = "RotMat")) && (!RotMatFun %in% ls(envir = .GlobalEnv))) {
    source(paste0(FunDir, "/", RotMatFun, ".R"))
  }
  if ((PPFun != "PPO") && (!RotMatFun %in% ls(envir = .GlobalEnv))) {
    source(paste0(FunDir, "/", PPFun, ".R"))
  }

  RotMatFun <- match.fun(RotMatFun, descend = TRUE)
  PPFun <- match.fun(PPFun, descend = TRUE)

  paramList$X <- X
  paramList$y <- y
  sparseM <- do.call(RotMatFun, paramList)
  for (ni in unique(sparseM[, 2L])) {
    lrows <- which(sparseM[, 2L] == ni)
    paramList$X <- X[, sparseM[lrows, 1L], drop = FALSE]
    sparseM[lrows, 3L] <- do.call(PPFun, paramList)
  }

  paramList$X <- NULL
  paramList$y <- NULL

  colnames(sparseM) <- c("Variable", "Number", "Coefficient")
  return(sparseM)
}


#' Samples a p x p uniformly random rotation matrix
#'
#' Samples a p x p uniformly random rotation matrix via QR decomposition
#' of a matrix with elements sampled iid from a standard normal distribution.
#'
#' @param p The columns of an n by p numeric matrix or data frame.
#'
#' @return A p x p uniformly random rotation matrix.
#'
#' @seealso \code{\link{RotMatPPO}} \code{\link{RotMatRand}} \code{\link{RotMatRF}} \code{\link{RotMatMake}}
#'
#' @examples
#' set.seed(220828)
#' (RandRot(10))
#'
#' @keywords tree rotation
#' @importFrom stats rnorm
#' @export
RandRot <- function(p) {
  RotMat <- qr.Q(qr(matrix(stats::rnorm(p^2), p, p)))
  if (det(RotMat) < 0) {
    RotMat[, 1] <- -RotMat[, 1]
  }
  return(RotMat)
}

# roxygen2::roxygenise()