R/util.R

In keyuan/bitphyloR: bitphylogeny in R

#' Convert beta random variables to edges in stick-breaking process
#' @param sticks A beta random variable vector
#' @return A vector of edges
#' @export
SticksToEdges <- function(sticks){
  return(1.0 - cumprod(1.0 - sticks))
}


#' Trim leading or trailing zeros in a vector
#' @param x A numerical vector
#' @param trim A string with "f" representing trim from front and "b" to
#' trim from back. Default is "fb", trim zeros from both front and back of the
#' array.
#' @return The result of trimming the input
#' @export
TrimZeros <- function(x, trim = "fb") {
  nonZeroIds <- which(x != 0)
  if (length(nonZeroIds) == 0) {
    return(c())
  }

  if (trim == "fb") {
    return(x[min(nonZeroIds):max(nonZeroIds)])
  } else if (trim == "f") {
    return(x[min(nonZeroIds):length(x)])
  } else {
    return(x[1:max(nonZeroIds)])
  }
}

#' Sample from a user supplied distribution with slice sampler
#' @param initX initial parameter
#' @param logprob a user supplied log-distribution function
#' @param sigma step size for slice construction
#' @param setpOut use setp out procedure or not
#' @param maxStepsOut maximum step out number
#' @param compwise use component wise update or not
#' @param verbose disply steps out&in
#' @param ... additional args for logprob
#' @return a sample of parameter
#' @export
SliceSampler <- function(initX,
                         logprob,
                         sigma = 1.0,
                         stepOut = T,
                         maxStepsOut = 100,
                         compwise = F, verbose = F, ...) {

  DirectionSlice <- function(direction, initX) {

    DirLogProb <- function(z) {
      return(logprob(direction*z + initX, ...))
    }

    llhS <- log(runif(1)) + DirLogProb(0)
    upper <- sigma * runif(1)
    lower <- upper - sigma
    lowerStepsOut <- 0
    upperStepsOut <- 0
    if (stepOut) {

      while (DirLogProb(lower) > llhS && lowerStepsOut < maxStepsOut) {
        lowerStepsOut <- lowerStepsOut + 1
        lower <- lower - sigma
      }
      while (DirLogProb(upper) > llhS && upperStepsOut < maxStepsOut) {
        upperStepsOut <- upperStepsOut + 1
        upper <- upper + sigma
      }
    }

    stepsIn <- 0
    while(T) {
      stepsIn <- stepsIn + 1
      newZ <- (upper - lower) * runif(1) + lower
      llhNew <- DirLogProb(newZ)
      if (is.nan(llhNew)) {
        cat("%f, %f, %f", newZ, initX, direction*newZ + initX, logprob(initX, ...))
        stop("Slice sampler got a NaN")
      }
      if (llhNew > llhS) {
        break
      } else if (newZ < 0) {
        lower <- newZ
      } else if (newZ > 0) {
        upper <- newZ
      } else {
        stop("Slice sampler shrank to zero")
      }
    }

    if (verbose) {
      cat("Steps Out:", lowerStepsOut, upperStepsOut, " Steps In:", stepsIn)
    }
    return(direction*newZ + initX)
  }

  dims <- length(initX)
  if (compwise) {
    ordering <- sample(dims)
    newX <- Reduce(
      rbind,
      Map(
        function(d) {
          direction = array(0, dim = dims)
          direction[d] = 1
          return(direction*DirectionSlice(direction, initX))
        },
        ordering),
      c())
    return(colSums(newX))
  } else {
    direction <- rnorm(dims)
    direction <- direction / sqrt(sum(direction^2))
    return(DirectionSlice(direction, initX))
  }
}

#' Bound samples generated by rbeta
#' @param ... rbeta parameters
#' @return bounded beta random variable
#' @export
BoundBeta <- function(...){
  return(
    (1-.Machine$double.eps)*(rbeta(...) - 0.5) + 0.5
    )
}


dbb <- function(x, N, u, v) {
  beta(x+u, N-x+v)/beta(u,v)*choose(N,x)
}

pbb <- function(q, N, u, v) {
  sapply(q, function(xx) sum(dbb(0:xx, N, u, v)))
}

qbb <- function(p, N, u, v) {
  pp <- cumsum(dbb(0:N, N, u, v))
  sapply(p, function(x) sum(pp < x))
}

rbb <- function(n, N, u, v) {
  p <- rbeta(n, u, v)
  rbinom(n, N, p)
}



repmat <- function(x, n, m){
  x <- as.matrix(x)
  return(kronecker(matrix(1,n,m),x))
}

ConvertFunctionNameToVariableName <- function(f) {
  res <- unlist(strsplit(gsub("(.)([[:upper:]])", "\\1 \\2", f), " "))
  res[2] <- paste(tolower(substring(res[2], 1, 1)), substring(res[2], 2), sep = "" )
  return(paste(res[-1], collapse = "" ))
}

lexcmp <- function(vec1, vec2) {
  index <- which.min(vec1 == vec2)  # find the first diff
  sign(vec1[index] - vec2[index])   # assumes numeric
}

ComparePath <- function(x, y) {
  if (is.null(x) && is.null(y)) {
    return(0)
  } else if (is.null(x)) {
    return(1)
  } else if (is.null(y)) {
    return(-1)
  }

  n1 <- length(x)
  n2 <- length(y)
  if (n1 > n2) {
    yy <- rep(0, n1)
    yy[1:n2] <- y
    return(lexcmp(yy, x))
  } else if (n1 <n2) {
    xx <- rep(0, n2)
    xx[1:n1] <- x
    return(lexcmp(y, xx))
  } else {
    return(lexcmp(y, x))
  }
}