R/learning_EM.R

In franciscorichter/emphasis: Evolutionary Modeling on Phylogenetic Applications with Simulations and Importance Sampling

#' function to perform one step of the E-M algorithm
#' @param brts vector of branching times
#' @param init_pars vector of initial parameter files
#' @param sample_size number of samples
#' @param maxN maximum number of failed trees
#' @param soc number of lineages at the root/crown (1/2)
#' @param max_missing maximum number of species missing
#' @param max_lambda maximum speciation rate
#' @param lower_bound vector of lower bound values for optimization, should
#' be equal in length to the vector of init_pars
#' @param upper_bound vector of upper bound values for optimization, should
#' be equal in length to the vector of init_pars
#' @param xtol_rel relative tolerance for optimization
#' @param num_threads number of threads used.
#' @return a list with the following components: 
#' \itemize{
#'  \item{trees}{list of trees}
#'  \item{rejected}{number of rejected trees}
#'  \item{rejected_overruns}{number of trees rejected due to too large size}
#'  \item{rejected_lambda}{number of trees rejected due to lambda errors}
#'  \item{rejected_zero_weights}{number of trees rejected to zero weight}
#'  \item{time_elapsed}{time used}
#'  \item{weights}{vector of weights}
#'  \item{fhat}{vector of fhat values}
#'  \item{logf}{vector of logf values}
#'  \item{logg}{vector of logg values}
#' }
#' @export
e_cpp <- function(brts, init_pars, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads=1) {
  .Call('_emphasis_rcpp_mce', PACKAGE = 'emphasis', brts, init_pars, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads)
}

augmentPD_old <- function(phylo, pars, maxN, max_missing, lower_bound, upper_bound) {
    num_threads=1
    sample_size = 1
    init_pars = pars
    soc = 2 
    max_lambda = max_missing*100
    xtol_rel = 0.00001
    brts = ape::branching.times(phylo)
    # brts = sort(max(brts) - brts)
    # brts = brts[-1]
  .Call('_emphasis_rcpp_mce', PACKAGE = 'emphasis', brts, init_pars, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads)
}


#' Monte Carlo Grid Calculation
#'
#' This function performs a Monte Carlo Expectation-Maximization (MCEM) grid calculation
#' for given parameters, using Rcpp for improved performance. It is designed to work
#' with the 'emphasis' package.
#'
#' @param pars_R Numeric matrix of parameters for the MCEM algorithm.
#' @param brts Numeric vector of branching times.
#' @param sample_size Integer, the number of samples to draw in the Monte Carlo simulation.
#' @param maxN Integer, the maximum number of iterations for the EM algorithm.
#' @param soc Numeric, scale parameter for optimization control. Defaults to 2.
#' @param max_missing Numeric, the maximum allowed missing information.
#' @param max_lambda Numeric, the maximum allowed lambda value.
#' @param lower_bound Numeric vector, the lower bounds for each parameter.
#' @param upper_bound Numeric vector, the upper bounds for each parameter.
#' @param xtol_rel Numeric, the relative tolerance for optimization convergence. Defaults to 0.00001.
#' @param num_threads Integer, the number of threads to use for parallel computation.
#'
#' @return Returns the result of the MCEM grid calculation as a numeric matrix.
#' @export
#' @examples
#' # Define parameters for the MCEM grid calculation
#' pars_R <- matrix(runif(20), nrow = 5) # Example parameter matrix
#' brts <- runif(10, 0, 1) # Example branching times
#' sample_size <- 100
#' maxN <- 1000
#' lower_bound <- rep(0, ncol(pars_R))
#' upper_bound <- rep(1, ncol(pars_R))
#' result <- mcemGridCalculation(pars_R, brts, sample_size, maxN, 2, 1e4, 1e4, lower_bound, upper_bound, 0.00001, 1)
#' print(result)
mcGrid <- function(pars_R, brts, sample_size, maxN, soc=2, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel=0.00001, num_threads) {
  .Call('_emphasis_rcpp_mce_grid', PACKAGE = 'emphasis', pars_R, brts, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads)
}



rcpp_mce_grid <- function(pars_R, brts, sample_size, maxN, soc=2, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel=0.00001, num_threads) {
  .Call('_emphasis_rcpp_mce_grid', PACKAGE = 'emphasis', pars_R, brts, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads)
}

rcpp_mce_grid_factorial <- function(pars_R, brts, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads) {
  .Call('_emphasis_rcpp_mce_grid_factorial', PACKAGE = 'emphasis', pars_R, brts, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads)
}

#' function to perform one step of the E-M algorithm
#' @param brts vector of branching times
#' @param init_pars vector of initial parameter files
#' @param sample_size number of samples
#' @param maxN maximum number of failed trees
#' @param soc number of lineages at the root/crown (1/2)
#' @param max_missing maximum number of species missing
#' @param max_lambda maximum speciation rate
#' @param lower_bound vector of lower bound values for optimization, should
#' be equal in length to the vector of init_pars
#' @param upper_bound vector of upper bound values for optimization, should
#' be equal in length to the vector of init_pars
#' @param xtol_rel relative tolerance for optimization
#' @param num_threads number of threads used.
#' @param copy_trees if set to true, the trees generated are returned as well
#' @param rconditional R function that evaluates the GAM function.
#' @return a list with the following components: 
#' \itemize{
#'  \item{trees}{list of trees}
#'  \item{rejected}{number of rejected trees}
#'  \item{rejected_overruns}{number of trees rejected due to too large size}
#'  \item{rejected_lambda}{number of trees rejected due to lambda errors}
#'  \item{rejected_zero_weights}{number of trees rejected to zero weight}
#'  \item{estimates}{vector of estimates}
#'  \item{nlopt}{nlopt status}
#'  \item{fhat}{vector of fhat values}
#'  \item{time}{time elapsed}
#'  \item{weights}{vector of weights}
#' }
#' @export
em_cpp <- function(brts, pars, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads, copy_trees, rconditional = NULL) {
  .Call('_emphasis_rcpp_mcem', PACKAGE = 'emphasis', brts, pars, sample_size, maxN, soc, max_missing, max_lambda, lower_bound, upper_bound, xtol_rel, num_threads, copy_trees, rconditional)
}

#' function to perform one step of the E-M algorithm
#' @param e_step result of e_step function, as a list
#' @param init_pars vector of initial parameter values
#' @param plugin string indicating plugin used, currently available: 'rpd1' and
#' 'rpd5c'
#' @param lower_bound vector of lower bound values for optimization, should
#' be equal in length to the vector of init_pars
#' @param upper_bound vector of upper bound values for optimization, should
#' be equal in length to the vector of init_pars 
#' @param xtol_rel relative tolerance for optimization
#' @param num_threads number of threads used.
#' @param rconditional R function that evaluates the GAM function.
#' @return list with the following entries:
#' \itemize{
#'  \item{estimates}{vector of estimates}
#'  \item{nlopt}{nlopt status}
#'  \item{time}{used computation time} 
#' }
#' @export
m_cpp <- function(e_step, init_pars, plugin, lower_bound, upper_bound, xtol_rel, num_threads, rconditional = NULL) {
  .Call('_emphasis_rcpp_mcm', PACKAGE = 'emphasis', e_step, init_pars, plugin, lower_bound, upper_bound, xtol_rel, num_threads, rconditional)
}