expokit_wrapalldgexpv_tvals: Run EXPOKIT's dgexpv on one or more t-values

View source: R/wrappers.r

expokit_wrapalldgexpv_tvalsR Documentation

Run EXPOKIT's dgexpv on one or more t-values


The function runs EXPOKIT's dgexpv function on a Q matrix and one or more time values. If Qmat is NULL (default), a default matrix is input.


  expokit_wrapalldgexpv_tvals(Qmat = NULL, tvals = c(2.1),
    inputprobs_for_fast = NULL, transpose_needed = TRUE,
    transform_to_coo_TF = TRUE, coo_n = NULL,
    force_list_if_1_tval = FALSE, check_for_0_rows = TRUE)



an input Q transition matrix


one or more time values to exponentiate by (doesn't have to literally be a time value, obviously)


If NULL (default), the full probability matrix (Pmat) is returned. However, the full speed of EXPOKIT on sparse matrices will be exploited if inputprobs_for_fast=c(starting probabilities). In this case these starting probabilities are input to myDMEXPV directly, as v, and w, the output probabilities, are returned.


If TRUE (default), matrix will be transposed (apparently EXPOKIT needs the input matrix to be transposed compared to normal)


Should the matrix be tranposed to COO? COO format is required for EXPOKIT's sparse-matrix functions (like dmexpv and unlike the padm-related functions. Default TRUE; if FALSE, user must put a COO-formated matrix in Qmat. Supplying the coo matrix is probably faster for repeated calculations on large matrices.


If a COO matrix is input, coo_n specified the order (# rows, equals # columns) of the matrix.


Default FALSE, but set to TRUE if you want a single matrix to be returned inside a list


If TRUE or a numeric value, the input Qmat is checked for all-zero rows, since these will crash the FORTRAN wrapalldmexpv function. A small nonzero value set to check_for_0_rows or the default (0.0000000000001) is input to off-diagonal cells in the row (and the diagonal value is normalized), which should fix the problem.


NOTE: DGEXPV vs. DMEXPV. According to the EXPOKIT documentation, DGEXPV should be faster than DMEXPV, however DMEXPV runs an accuracy check appropriate for Markov chains, which is not done in DGEXPV.


tmpoutmat the output matrix, if 1 t-value is input; list_of_matrices_output, if more than 1 t-value is input; to get a single output matrix in a list, set force_list_if_1_tval=TRUE


Nicholas J. Matzke nickmatzke.ncse@gmail.com and Drew Schmidt schmidt@math.utk.edu

See Also



# Example:
# Make a square instantaneous rate matrix (Q matrix)
# This matrix is taken from Peter Foster's (2001) "The Idiot's Guide
# to the Zen of Likelihood in a Nutshell in Seven Days for Dummies,
# Unleashed" at:
# \url{http://www.bioinf.org/molsys/data/idiots.pdf}
# The Q matrix includes the stationary base freqencies, which Pmat
# converges to as t becomes large.
Qmat = matrix(c(-1.218, 0.504, 0.336, 0.378, 0.126, -0.882, 0.252, 0.504, 0.168,
0.504, -1.05, 0.378, 0.126, 0.672, 0.252, -1.05), nrow=4, byrow=TRUE)

# Make a series of t values
tvals = c(0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 2, 5, 14)

# Exponentiate each with EXPOKIT's dgexpv (should be fast for large sparse matrices)
for (t in tvals)
	Pmat = expokit_dgexpv_Qmat(Qmat=Qmat, t=t, transpose_needed=TRUE)
	cat("\n\nTime=", t, "\n", sep="")

# DMEXPV and DGEXPV are designed for large, sparse Q matrices (sparse = lots of zeros).
# DMEXPV is specifically designed for Markov chains and so may be slower, but more accurate.

# DMEXPV, single t-value

# DGEXPV, single t-value
expokit_wrapalldgexpv_tvals(Qmat=Qmat, tvals=tvals[1], transpose_needed=TRUE)
expokit_wrapalldgexpv_tvals(Qmat=Qmat, tvals=2)

# These functions runs the for-loop itself (sadly, we could not get mapply() to work
# on a function that calls dmexpv/dgexpv), returning a list of probability matrices.

# DGEXPV functions
list_of_P_matrices_dgexpv = expokit_wrapalldgexpv_tvals(Qmat=Qmat,
tvals=tvals, transpose_needed=TRUE)

nmatzke/rexpokit documentation built on Nov. 28, 2023, 9:35 p.m.