TexShop engines for Knitr w/ and w/o pstricks

Two years ago, I've posted a link that explains how to set up a painless compilation of Sweave documents in TexShop: http://cameron.bracken.bz/sweave-for-texshop

Since I've switched to Knitr last year, I wanted to share a couple of tricks of using a TexShop+Knitr combination.

First, setting up a Knitr engine in TexShop is analogous to the recipe from the above link: create a file Knitr.engine in your ~/Library/TexShop/Engines/ directory with the following content:

#!/bin/bash

export PATH=$PATH:/usr/texbin:/usr/local/bin
Rscript -e "library(knitr); knit('$1')"

Now if you (re-)open TexShop, you should have Knitr in your list of engines located next to the Typeset button on the top. This worked well for me until today, when I needed to combine Knitr code with some old LaTeX code of mine that used "pstricks" package. Since pdflatex does not work with "pstricks" (there something called PDFTricks, but I haven't tried it), I needed to create another engine in my ~/Library/TexShop/Engines/ directory that I called KnitrDviPs:

Rscript -e "library(knitr); knit('$1')"
latex "${1%.*}"
bibtex "${1%.*}"
latex "${1%.*}"
latex "${1%.*}"
dvips -t letter -o "${1%.*}.ps" "${1%.*}"
pstopdf "${1%.*}.ps"

This engine goes through a less efficient "dvi->ps->pdf" route. However, this still doesn't do the trick just yet, because by default Knitr does not produce postscript files of all embedded figures. So you have to add a "postscript" device to your Knitr document using Knitr chunk options. For example, somewhere near the beginning of your Knitr document you can put the following command:

opts_chunk$set(fig.path='figures/', fig.align='center', dev=c('pdf', 'postscript'))

Now you have two Knitr engines: one that goes through pdflatex and another that goes through "dvi->ps->pdf".

Fundraising for victims of May 2nd, 2014 events in Odessa, Ukraine (updated on May 19, 2014)

Many of you have seen the news about violent clashes that occurred in Odessa, Ukraine on May 2nd, 2014. During these clashes, my very close childhood friend was shot from a machine gun and now he is fighting for his life in one of the hospitals in Odessa. I don’t know all the details yet, but it appears that my friend was wounded by a pro-russian activist during the initial stage of the unrest, when heavily armed pro-russian activists attacked a peaceful demonstration for the unity of Ukraine. The crossed out information hasn't been confirmed. I still don't know all the details, but if you feel that you were mislead by the crossed out sentence, send me a note by email.

My friend’s family has enough financial resources to give him medical care, but my friend’s father says that he sees many families of victims of the clashes that don’t have money to pay for medical expenses. The Ukrainian government promised to help all victims, but their help may arrive too late, if it ever arrives. Since I have a unique opportunity to reach out to the victims’ families in need through my relatives and friends in Odessa, I’ve decided to try to raise funds to provide medical care to the wounded in clashes on May 2nd, 2014. The funds will be distributed to the wounded without
any regard to their political convictions.

If you would like to help me in this fundraising campaign, please donate to my PayPal account using a link on my web page:

http://www.stat.washington.edu/vminin/

-Vladimir

Update 1 (05/06/2014): funds raised so far - ~$1,500; wholehearted thanks to everyone, you are amazing! I've made the first transfer of $1,000 to Odessa today.

Update 2 (05/06/2014): I was on the phone with Odessa today. Hospitals are in dire need of many pieces of medical equipment to treat hundreds of wounded. There is a local fundraising effort going on to buy this equipment and we were offered to join this effort. This would be a slight deviation of the initial plan of distributing funds to the wounded directly. However, I was told that there is where the help is needed the most at the moment, so I support this suggestion. Still, I will ask my contacts in Odessa to give money directly to the families of wounded that need assistance first, and use the rest for medical equipment. Important: if you are not comfortable with this plan or have any hesitation about how your donation will be spent, please contact me by email so I can give you more information and/or refund you.

Update 3 (05/09/2014): funds raised - ~$2,000; again, thank you!

Update 4 (05/09/2014): I am happy to report that I got in touch with a local volunteer in Odessa, whom I personally know and who is involved in distributing financial assistance to the wounded on May 2nd in Odessa. So we are back to the original plan -- providing money directly to the families of wounded!

Update 5 (05/09/2014): Since we have reached a plateau in the total funds raised, I've decided to take a pause, to stop soliciting/receiving donations, and to concentrate on distributing the money we have collected so far. If I have time and energy, I may try for a second round of fundraising. I will be updating this post with details of how the donations are being spent in Odessa.

Update 6 (05/19/2014): I've got an update on how the first transfer (~1,000) to Odessa was spent! First, medical bills, amounting to ~$485 were payed for a person in need who was hurt on May 2, 2014 in Odessa and remains in critical condition. Next, our local volunteer used ~$460 to purchase painkillers for the wounded that still need medical attention. The small remainder of the transferred sum was used to help offset the cost of medical treatment for children of three military personnel relocated from Crimea after the Russian invasion.

Guest post: Where did the Chapman-Kolmogorov name come from?

This is a guest post by Peter Guttorp. Peter and I are working on updating Peter's book Stochastic Modeling of Scientific Data. To provide the reader with productive distractions, we are sprinkling the book with "historical boxes." Also, let's be honest, it keeps us entertained as well, because we are both math history junkies. Well, Peter does real research in the history of mathematics (e.g., here) , so he is a real historian, but I am definitely a junkie. While working on one of the background chapters, we put a historical box for Chapman and Kolmogorov next to the first appearance of the Chapman-Kolmogorov equation, but then wondered why we did this, because Peter knew that neither Chapman nor Kolmogorov were responsible for its discovery. This is, of course, not surprising in light of the Stigler's law, but then who gave this equation its modern name? Peter started digging and this is what he has found so far:

The origin of the Chapman-Kolmogorov equation dates at least back to Bachelier’s 1900 dissertation on the theory of speculation (but it may be earlier than that). Bachelier did not give it a name.  

In physics it is sometimes called Smoluchovsky’s equation, after Smoluchovsky’s 1906 rediscovery of Brownian motion (preceded by Bachelier and by Einstein (1905)). Smoluchovksy’s work was extended by Chapman (1928) to the case of nonstationary transition kernels. It was independently developed by Kolmogorov in his fundamental 1933 Mathematische Annalen paper. 

But where did the name come from? The first uses of the Chapman-Kolmogorov name are in Feller’s 1940 Transactions of the AMS paper, where he says that it is “known as the equation of Chapman and Kolmogoroff” . The 1940 thesis of Lundberg also uses this term, referring to Frechet (1938) and Kolmorogov (1931). Going to Frechet’s book, he does not give the equation a name, but designates it “as a mnemotechnical tool” with the letter I. 

In 1936, Feller writes in his Mathematische Annalen paper “Finally, from the rule of combination of probabilities, what Chapman and von Smoluchowski called the fundamental equality follows”. 

So what happened between 1936 and 1940? Did the name get invented in the corridors of the University of Stockholm (where Feller was a lecturer, Lundberg a PhD student, and Cramer the professor)?

If you know anything about the origin of the name Chapman-Kolmogorov equation, please get in touch with Peter and/or me. We are very curious!

Passing Armadillo matrices by value and by reference

Rcpp and RcppArmadillo make it relatively easy to integrate R code and  C++ code. In fact the resulting C++ code looks so much as R code that it is easy to forget some basic rules by which a C++ coder must play. This is especially true for someone like me, who wrote relatively little C++ code in the past.

Here is an example of a typical calculation that may get one in trouble. Suppose you have a function that takes two matrices from R, converts them into Armadillo matrices and passes these functions to other functions to do some repetitive calculations. For example, imagine passing transition and emission probability matrices to execute the forward-backward algorithm on multiple realizations of the same hidden Markov model.

Below you see my C++ code, saved in the file "matrix_code.cpp" that has two sets of functions: a) doStuff1(), twoMatrices1(); b) doStuff2(), twoMatrices2(). The functions twoMatrices1() and twoMatrices2() are exported to R via Rcpp attributes. The twoMatrices1() calls doStuff1() 100 times, while twoMatrices2() calls doStuff2() 100 times. The only difference between the two sets of functions is in that doStuff1() passes two Armadillo matrices by value, which requires copying these matrices. In contrast, doStuff2() passes Armadillo matrices by reference without calling copy-constructor under the hood.

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>  

using namespace Rcpp; 

// First passing Armadillo matrices by value
void doStuff1(arma::mat c, arma::mat d){
  double x = c(1,1) + d(1,1);
} 

// [[Rcpp::export]]
void twoMatrices1(NumericMatrix A, NumericMatrix B) {
    arma::mat M1(A.begin(), A.nrow(), A.ncol(), false);
    arma::mat M2(B.begin(), B.nrow(), B.ncol(), false);
 
    for (int i = 0; i < 100; i++){
      doStuff1(M1, M2);
    }
} 

// Now passign Armadillo matrices by reference
void doStuff2(const arma::mat& c, const arma::mat& d){
  double x = c(1,1) + d(1,1);
}

// [[Rcpp::export]]
void twoMatrices2(NumericMatrix A, NumericMatrix B) {
    arma::mat M1(A.begin(), A.nrow(), A.ncol(), false);
    arma::mat M2(B.begin(), B.nrow(), B.ncol(), false);
 
    for (int i = 0; i < 100; i++){
      doStuff2(M1, M2);
    }
}

Now I am going to call twoMatrices1() and twoMatrices2() from R, using rbenchmark package to compare their speed. I am going to use fairly large matrices (100x100), so the cost of moving them around is significant.

library(RcppArmadillo)
library(rbenchmark)
 
sourceCpp("matrix_code.cpp")
 
x = matrix(1,100,100)
y = matrix(1,100,100)
 
res <- benchmark(twoMatrices1(x, y), twoMatrices2(x, y), 
                 columns = c("test", "replications", "elapsed", "relative"), 
                 order="relative", replications=1000)
print(res)
 
                test replications elapsed relative
2 twoMatrices2(x, y)         1000   0.006    1.000
1 twoMatrices1(x, y)         1000   0.829  138.167

Created by Pretty R at inside-R.org

As you can see, passing two 100x100 Armadillo matrices by value is more than 100 times slower than passing these matrices by reference. Perhaps not a very surprising outcome, but an important one to keep in mind when porting R code to C++.

Dissertation Acknowledgments

Jonathan Eisen posted his dissertation acknowledgments. I think this is a great idea. Maybe it is just me, but when I read someone's dissertation, not necessarily as a committee member, I really enjoy reading the acknowledgements. For me, these personal comments place the work in the context of "real life" and personal relationships. 

Here are my acknowledgments:

First and foremost, I would like to thank my mom and dad for creating many opportunities for my intellectual growth. Feeling my parents' constant love and support makes my scientific journey so much more meaningful! I am very thankful to my middle and high school mathematics teachers, Maya Radionovna Matihina and Svetlana Efimovna Gaydamak, for introducing me to the beautiful world of mathematical reasoning and abstract thinking.  I am very grateful to my undergraduate advisor, Galina Evgenievna Samkova, for patiently listening to my naive ideas and taking seriously my modest attempts to produce new mathematical results. Special thanks goes to Stephen Krone, who hooked me onto mathematical biology by regularly taking me to the bioinformatics seminars at the University of Idaho. My first scientific collaborators, Zaid Abdo and Paul Joyce, became my very good friends and influenced my scientific career more than they probably realize.    

UCLA Department of Biomathematics has been my second home for the last four years. I would like to thank the departmental staff, David Tomita, Martha Riemer, Wendy Hagar, and Fred Hughes, for taking care of administrative hurdles so well. During my first years at UCLA, I received a lot of advice and guidance from my fellow Biomathematics students. I am very grateful to Robert Rovetti for being a wonderful colleague and a good friend. Robert thoroughly read and constructively criticized Chapters 3 and 4 of this dissertation. Intellectual discussions with Benjamin Redelings, John O'Brien, and Alex Alekseyenko were certainly worth inhaling cigarette smoke and drinking really bad coffee from vending machines. 

All members of my dissertation committee played important roles in my graduate school development. I enjoyed collaborating with Christina Kitchen, who is very passionate about both statistics and biology. I am grateful to Janet Sinsheimer for giving me many opportunities to assist her in graduate and undergraduate level teaching. Elliot Landaw provided me with excellent academic guidance at the beginning of my studies and was extremely supportive throughout all these years. Ken Lange has frequently given me always timely and wise pieces of advice. Chapter 6 would not be possible without
his enthusiasm about Markov chain induced counting processes and his encouragement to pursue this line of research. 

The chair of my dissertation committee, Marc Suchard, is the greatest advisor one could wish for! I am very grateful to Marc for introducing me to interesting problems and for not revealing solutions to these problems, giving me an opportunity to become an independent researcher. It is very difficult to summarize everything I have learned from Marc in one paragraph. Perhaps the most important lesson that Marc taught me is how not to be afraid: not to be afraid to attack challenging problems, not to be afraid to seek rigorous solutions for these problems, and not to be afraid to go somewhere no one has gone before.    

This dissertation would never be finished without support of my wife, Yuliya, who patiently tolerated my "writing weekends'' and many other strange habits of a finishing Ph.D. student.  

Guest post: Functional Programming in R with an HMM example

Below is a guest post by Chris Glazner, who works on his Ph.D. in Statistics with Elizabeth Thompson here at UW. I especially recommend checking out Chris' slick example of using Reduce for HMM forward-backward algorithm.

Chris writes:

I spent some time over the summer looking at ways of using functional programming to implement the statistical methods in my dissertation research. Functional programming is a style which treats functions as the basic building block of computations. The flow of a program is expressed as the composition of functions, rather than a sequence of steps which manipulate variables. Programming this way tends to feel a little more like doing math, so I've found it to be a good fit for statistical programming. If you can describe your program as a function between sets (say, from a data vector to a confidence interval) you've already figured out the overall structure of a functional program to implement it. 

R is a functional language, meaning that functions can be given names, created anonymously, and passed to other functions. Functions whose inputs or outputs are themselves functions are called "higher order functions". A number of higher order functions are included in the R base library, including the popular apply family: apply, sapply, lapply, and mapply. These functions take a function and a list or vector and apply the function to every element of the list or vector. They replace for loops that look like this:  

for (i in 1:length(x)) f(x[i])  

Created by Pretty R at inside-R.org

where f is any suitable function. Replacing a loop like this with a function call simplifies code and makes it easier to understand.
Other higher order functions express other common patterns. Reduce captures the behavior of an accumulating loop like: 

y <- startValue
for (i in 1:length(x)) y <- g(y,x[i]) 

Created by Pretty R at inside-R.org 

This pattern may not seem as general as apply, but it encompasses many common calculations. If g is `+` and startValue is 0, then the call to Reduce is equivalent to sum(x). If g is `*` and startValue is 1, we have product(x). Reduce is known by many names in other languages and is a powerful abstraction for expressing computations. 

A nice application of Reduce arises in the Hidden Markov Models (HMMs). An HMM is an unobserved Markov chain along with data which are generated independently according to the hidden Markov state at each step. The goal is usually to infer something about the hidden process from the data. The standard method for calculating the marginal probabilities of the hidden states at each step is known as the forward-backward algorithm, which is explained in a nice tutorial by LR Rabiner. We can use Reduce to express this algorithm concisely. 

As an example we will implement the simple Occasionally Dishonest Casino model described in a classic Biologial Sequence Analysis book by Durbin, Eddy, Krogh, and Mitchison. This casino has a fair die and a biased die which rolls three 50% of the time. We observe repeated rolls of a single die, but we never know which one. The casino switches the dice according to a two-state Markov chain. We want to infer which die is in use at a given time. 

The model is fully specified by the Markov transition matrix, the emission distribution of the data conditional on the hidden state, and the initial distribution. Since this is R, we will specify these using vectors and matrices, but we could use functions if we wanted.

tMat <- t(matrix(c(0.95, 0.05 # fair -> biased
                 , 0.10, 0.90 # biased -> fair
                 ),2,2))
 
eMat <- cbind( rep(1/6,6)            # fair
              ,c(.1,.1,.5,.1,.1,.1)) # biased
initDist <- c(0.5,0.5)

Created by Pretty R at inside-R.org 

Now we need two functions to execute the forward and backward passes. These can be concisely expressed in matrix notation: 

fwd <- function(distr,obs) normalize(( t( tMat ) %*% distr ) * eMat[obs,] )
backwd <- function(obs,distr) normalize( tMat %*% ( distr * eMat[obs,] ))
normalize <- function(x) x/sum(x)

Created by Pretty R at inside-R.org 

The fwd function takes two inputs, a distribution on the hidden state space and the observation at the current node, and returns another distribution. Reduce uses the output of one call as the distr argument to the next call, with the obs argument taken from the list of data observations. We set the accum option to Reduce so that it collects the intermediate results from each call. Assuming the data are in a vector called dieRolls, we write 

dieRolls <- c(3,3,3,4,4,6,2,3,4)
forwardProbs <- Reduce(fwd,dieRolls,init=initDist,acc=T)
backwardsProbs <- Reduce(backwd,dieRolls,init=normalize(c(1,1)),acc=T,right=T)
marginalProbs <- mapply(`*`,forwardProbs,backwardsProbs)

Created by Pretty R at inside-R.org 

Since the backwards algorithm starts at the end of the data vector, we use the right argument (and switch the order of the arguments to backwd). The initial values for both passes are as described in the Rabiner tutorial. The higher order function mapply is used to combine the resulting lists with pointwise multiplication. 

The Rabiner tutorial presents this algorithm in terms of loops over arrays. This style is easy to translate into computer code, but it doesn't provide any intuition as to what is happening; I have a hard time connecting the abstract statistical ideas to the implementation. In contrast, Reduce expresses the recursive nature of the dynamic programming problem. Coding the algorithm becomes part of the research process, rather than a chore. Of course, some people can grasp the loops and arrays just fine, but if you are like me and have spent more time taking math and statistics classes than programming, you might want to check out functional programming.

RcppArmadillo and R-forge

Jane Lange and I have recently finished a paper on continuous-time hidden Markov models and published an accompanying R package cthmm on R-forge. Jane's package cthmm relies heavily on Rcpp and RcppArmadillo packages. Somewhat surprisingly, the last dependency caused compilation problems on R-forge. Specifically, the package compiled only on Mac OS X, but not on Linux or Windows. Since Jane and I use Macs, we didn't have any problems compiling the package locally on our computers.

It turns out that the root of the problem was in using Armadillo's .inv() function to invert complex matrices. Armadillo library outsources some of its linear algebra calculations to lapack, which should be no problem since R also relies on lapack (who doesn't, right?). But for some reason r-project.org systems have only an abbreviated version of the lapack library, missing some functions for complex matrix calculations. This problem is discussed here.

Following suggestions from the above link, my solution was to download missing fortran routines (zgetri.f, ztrti2.f, ztrtri.f) from http://www.netlib.org/lapack/ and add them to /src directory of our package. A bit ugly, but worked and the package compiled on R-forge on all architectures.

P.S. As of July 22, 2012, our cthmm failed to compile on R-forge again, but we believe that this has to do with dependent packages (looks like survival package was not available on Windows). We hope to resolve this problem shortly.

Update on the P.S.: It looks like the problem with the survival package was resolved by the R-forge compilers. Our cthmm compiles now.

Blog name change

As I find myself working on problems outside of evolutionary biology more often, I decided that a blog name change would be in place. This doesn't mean that I am planning to stop working on phylogenetics and population genetics though.

TexShop+Sweave

A useful link for TexShop and Sweave lovers:

http://cameron.bracken.bz/sweave-for-texshop

Packaging and exposing Rcpp functions

In my previous posts, I used the inline package to test my Rcpp functions. This strategy becomes very limited as your software project grows even modestly beyond a couple of functions. So a natural solution is to create an R package with your C++ functions and R code. Creating R packages that depend on the Rcpp package is very well documented in the corresponding Rcpp vignette.

I've created a package phyloRcppExamples. It is hosted on R-forge, so you can browse the source code online or download the source via svn. I hope that through this package I will be able to illustrate more elaborate examples of using Rcpp for phylogenetic computations. I am not sure if this package will ever make its way to CRAN, but the good people at R-forge provide binaries for the most recent version of R (not tested yet on phyloRcppExamples).

Now, since we are not using inline anymore, every time I add a new C++ function to my package, I need to expose it to R (e.g. via .Call()). This usually requires writing a small C++ wrapper, which becomes cumbersome. Rcpp modules come to rescue. Below is an illustration of this concept.

In the /src directory of the phyloRcppExamples packages, I have a file two_state_model.cpp. One of the functions in this file calculates transition probabilities of a two-state continuous-time Markov chain:

//Computes finite time transition probabilities using analytic formulae.

arma::mat trans_prob(double lambda_01, double lambda_10, double time){

  double total_rate = lambda_01 + lambda_10;               

  arma::mat prob_matrix = arma::mat(2,2);
  
  prob_matrix(0,0) = (lambda_10 + lambda_01*exp(-total_rate*time))/total_rate;
  prob_matrix(0,1) = (lambda_01 - lambda_01*exp(-total_rate*time))/total_rate;
  prob_matrix(1,0) = (lambda_10 - lambda_10*exp(-total_rate*time))/total_rate;
  prob_matrix(1,1) = (lambda_01 + lambda_10*exp(-total_rate*time))/total_rate;

  return prob_matrix;
}

One way to expose this function R is to write a wrapper function

RcppExport SEXP trans_prob_to_R(SEXP lambda_01, SEXP lambda_10, SEXP time){
  return wrap(trans_prob(as<double>(lambda_01), as<double>(lambda_10), as<double>(time)));
}

Another method uses Rcpp modules:

RCPP_MODULE(two_state){ 
  function( "trans_prob", &trans_prob );
}

Here is how the two methods are used within R:

> require(phyloRcppExamples)
> 
> forward.rate = 2.3
> backward.rate = 3.1
> elapsed.time = 0.03
> 
> ## via wrapper
> .Call("trans_prob_to_R", forward.rate, backward.rate, elapsed.time, PACKAGE="phyloRcppExamples")
           [,1]       [,2]
[1,] 0.93629903 0.06370097
[2,] 0.08585783 0.91414217
> 
> ## same as above, but using Rcpp modules
> 
> two.state.Rcpp = Module("two_state", PACKAGE="phyloRcppExamples")
> 
> two.state.Rcpp$trans_prob(forward.rate, backward.rate, elapsed.time)
           [,1]       [,2]
[1,] 0.93629903 0.06370097
[2,] 0.08585783 0.91414217
> 

Created by Pretty R at inside-R.org

Keep in mind that you can add multiple functions to the same module.

Phylogenetic likelihood with Rcpp+RcppArmadillo

In this post, I want to show how to implement Felsenstein's algorithm for calculating phylogenetic likelihood with the help of Rcpp and RcppArmadillo. I am going to use ape's representation of a phylogenetic tree. If you are not familiar with this data structure, look at this document. My implementation is very similar to the likelihood calculator of ape's ace function.

Problem formulation: given a rooted binary phylogenetic tree with branch length and a continuous-time Markov chain model of discrete character evolution (e.g. DNA, amino acids, morphological traits), compute the probability of observing certain characters at the tips of the tree. I am going to stop short of actually computing this probability. Instead, I will compute what's called partial likelihoods for each node in the tree. The actual likelihood is a weighted sum of the partial likelihoods at the root of the phylogeny.

C++ code with inline compilation:

require(inline)
require(RcppArmadillo)
 
src <- '
  /* Arguments: 
        tE: edge matrix of the ape object phylo
        nIN: number of internal nodes (not necessary, but convinient to get it from phylo)
        tS: integer vector of tip states (-1=missing value)
        nS: state space size (e.g. 2 for a binary trait)
        pM: array of probability matrices for each edge of the tree
 
     Two important assumptions:
        1. edges in the edge matrix and probability matrices are in the "pruningwise" order; 

        see ?reorder.phylo for more details

        2. tip state vector is ordered according to the tip numbering in the edge matrix
  */
 
  using namespace Rcpp;
 
  IntegerMatrix treeEdges(tE);
 
  // get the number of edges
  int numEdges = treeEdges.nrow();
  int numIntNodes = as<int>(nIN);
  IntegerVector tipStates(tS);
  int numStates = as<int>(nS);
  NumericVector vecProbMat(pM);
  arma::cube cubeProbMat(vecProbMat.begin(), numStates, numStates, numEdges, false);
 
  // get the number of tips in the tree
  int numTips = tipStates.size();
 
  // prepare a matrix for storing regular (backward) partial likelihoods
  arma::mat partialLike = arma::zeros<arma::mat>(numTips + numIntNodes, numStates);
 
  for (int i=0; i < numTips; i++){
    if (tipStates[i] == -1){// -1 denotes a missing value
      partialLike.row(i) = arma::ones<arma::rowvec>(numStates);
    }else{
      partialLike(i, tipStates[i]) = 1.0;
    }
  }
 
  // compute regular partial likelihoods for all internal nodes
  for (int i=0; i < numEdges; i+=2){      
    // parent=treeEdges(i,0) or treeEdges(i+1,0); treeEdges indices should be shifted by one
    partialLike.row(treeEdges(i,0)-1) = (partialLike.row(treeEdges(i,1)-1)*cubeProbMat.slice(i).t())%(partialLike.row(treeEdges(i+1,1)-1)*cubeProbMat.slice(i+1).t());            
  }
 
  return wrap(partialLike);
'  
 
 
partLike = cxxfunction(signature(tE="integer", nIN="integer", tS="integer", nS="integer", pM="numeric"),body=src, plugin="RcppArmadillo")

Created by Pretty R at inside-R.org

I am going to test the code with a simple two state model, for which transition probabilities are available in closed form:

two.state.trans.prob = function(forward.rate, backward.rate, elapsed.time){
  total.rate = forward.rate + backward.rate               
 
  return((matrix(c(rep(backward.rate,2),rep(forward.rate,2)),2,2) +
    matrix(c(forward.rate, -backward.rate, -forward.rate, backward.rate),2,2)*
    exp(-total.rate*elapsed.time))/total.rate)    
}

Created by Pretty R at inside-R.org

R wrapper function that takes care of the prerequisites:

two.state.part.like = function(my.tree, my.data, forward.rate, backward.rate){
 
  ## reorder the edges in the "pruningwise" order
  my.tree = reorder(my.tree, order = "pr")
 
  if (!("phylo" %in% class(my.tree)))
    stop("Error: object \"my.tree\" is not of class \"phylo\"")
 
  if (is.null(my.tree$edge.length))
    stop("Error: tree \" my.tree\" must have branch lengths.")
 
  ## reorder data on tips to match the order of the my.tree phylo object
  if (!is.null(names(my.data))) {
    if(!any(is.na(match(names(my.data), my.tree$tip.label)))){
      my.data = my.data[my.tree$tip.label]
    }else{
      warning('the names of argument "my.data" and the names of the tip labels
did not match: the former were ignored in the analysis.')
    }
  }
 
  ## prepare transition probability matrices (this of course can and should be done in C++ as well)
  prob.array = array(0, dim=c(2,2,length(my.tree$edge.length)))            
  for (i in 1:length(my.tree$edge.length)){
    prob.array[,,i] = two.state.trans.prob(forward.rate, backward.rate, my.tree$edge.length[i])          
  }            
 
  return(partLike(my.tree$edge, my.tree$Nnode, my.data, 2, prob.array))
}

Created by Pretty R at inside-R.org

Testing: 

First, let's simulate a tree with tip states. This can be done by several packages in R. I am going to use diversitree package:

require(diversitree)
 
set.seed(34344)
test.tree = tree.bisse(c(0.1, 0.1, 0.03, 0.03, 0.01, 0.07), x0=0, max.taxa = 50)

Created by Pretty R at inside-R.org

This produces the following tree:

Now, my function returns partial likelihoods for all nodes in the tree. I am going to extract only partial likelihoods corresponding to the root (recall that a weighted sum of these numbers gives the actual probability of observing the data). I am going to do the same using diversitree's likelhood calculator:

> test.like = two.state.part.like(test.tree, test.tree$tip.state, 0.02, 0.05)
> test.like[51,]
[1] 1.564918e-08 1.195940e-08
> 
> mk2.lik = make.mk2(test.tree, test.tree$tip.state)
> exp(mk2.lik(c(0.02,0.05),root=ROOT.BOTH))
[1] 1.564918e-08 1.195940e-08

Created by Pretty R at inside-R.org

The above agreement is encouraging.

R array to RcppArmadillo cube (updated)

A quick post here on how to convert a 3D array from R to the RcppArmadillo cube.

Rcpp does not support arrays yet (only vectors and matrices). However, RcppArmadillo has a proper 3D array class: arma::cube. The conversion can be done via Rcpp::NumericVector. The only aspect of this exercise that surprised me was the fact that Rcpp::NumericVector stores dimensions of an array, but does not allow one to extract this info. Hence, if you don't know in advance the dimensions of your array, the only solution I see is to pass the vector of dimensions to Rcpp together with the array. Not too ugly, but a bit awkward. Thanks to the commenter Eli below, who pointed out that dimensions of the array can be extracted via the .attr() method -- much prettier now. My code to illustrate this concept:

require(inline)
require(RcppArmadillo)
 
src <- '
  using namespace Rcpp;
 
  NumericVector vecArray(myArray);
  IntegerVector arrayDims = vecArray.attr("dim");

 
  arma::cube cubeArray(vecArray.begin(), arrayDims[0], arrayDims[1], arrayDims[2], false);
 
  //change one element in the array/cube
  cubeArray(0,0,0) = 518;  
 
  return(wrap(cubeArray));  
'
 
readCube = cxxfunction(signature(myArray="numeric"),body=src, plugin="RcppArmadillo")

Created by Pretty R at inside-R.org

A numerical test:

> set.seed(345)
> testArray = array(rnorm(18), dim=c(3,3,2))
> print(testArray)
, , 1
 
           [,1]        [,2]      [,3]
[1,] -0.7849082 -0.29059656 -0.927724
[2,] -0.2795144 -0.06753159  1.710771
[3,] -0.1614579 -0.63352041  1.654769
 
, , 2
 
          [,1]       [,2]       [,3]
[1,]  1.810483 -0.8496292 -1.4029422
[2,]  1.866772  0.3184496  0.5682982
[3,] -1.399833  0.9035913  1.0457561
 
> readCube(testArray)
, , 1
 
            [,1]        [,2]      [,3]
[1,] 518.0000000 -0.29059656 -0.927724
[2,]  -0.2795144 -0.06753159  1.710771
[3,]  -0.1614579 -0.63352041  1.654769
 
, , 2
 
          [,1]       [,2]       [,3]
[1,]  1.810483 -0.8496292 -1.4029422
[2,]  1.866772  0.3184496  0.5682982
[3,] -1.399833  0.9035913  1.0457561

Created by Pretty R at inside-R.org

Uniformization-based CTMC sampling using Rcpp+RcppArmadillo

Suppose you have a continuous-time Markov chain on a finite state space and you want to draw realizations of this chain conditional on starting and ending states of the chain during some finite time interval. For a good review of available algorithms to accomplish this task see (Hobolth and Stone, 2009). Notice that supplementary materials to the the above paper contain R code, written by Asger Hobolth, that implements all algorithms described in the main text.

Here, I am going to implement a uniformization method (Fearnhead and Sherlock, 2006; Lartillot, 2006; Mateiu and Rannala, 2006) using Rcpp and RcppArmadillo. In essence, I just Rcpp-ed Asger's code. I want to point out that my motivation is not necessarily to speed up Asger's R implementation. Instead, I wanted to have a standalone C++ code that I can use in my packages internally. For portability, I am going to use the inline package, so that one can copy and paste my code into an R console. 

I am fairly new to Rcpp and especially RcppArmadillo, so here are a couple of lessons that I've learned during this coding exercise:

1. Rcpp sugar doesn't have "sample()" implemented, so I had to write my own. Hence, the awkward "sampleOnce" function below.
2. The only way to convert between "colvec" and "rovec" in Armadillo is to use a transpose. Thanks to Jane Lange for pointing this possibility.
3. If you want to grow a vector, using Rcpp NumericVector class is not computationally efficient. According to Rcpp developers it is better to create a std::vector instance, use .push_back() to grow it and then create NumericVector from std:vector.

So here is the C++ code to use with inline:

inc <-
'int sampleOnce(arma::colvec weights, double rUnif) {
 
  double total = sum(weights);
  double cumProb=0;
 
  int i;
 
  for (i = 0; i < weights.n_elem; i++) {
    cumProb += weights(i) / total;
    if (rUnif < cumProb) break;
  }
  return(i);
}
'
 
src <- '
  using namespace Rcpp;
 
  NumericMatrix rM(rm);
  int numStates = rM.nrow();
  arma::mat rateMatrix(rM.begin(), numStates, numStates, false); // reuses memory and avoids extra copy
 
  int startState = as<int>(bS);
  int endState = as<int>(eS);
  double elapsedTime = as<double>(eT);
  double transProb = as<double>(tP);  
 
  RNGScope scope;
 
  // Set the rate of the dominating Poisson process
 
  double poissonRate = -1.0*arma::min(rateMatrix.diag());
  arma::mat dominTransProb(numStates,numStates);
  dominTransProb.eye();
  arma::cube powerTransProb = arma::cube(numStates, numStates,1);
  powerTransProb.slice(0) = dominTransProb;
  dominTransProb = dominTransProb + rateMatrix/poissonRate;
 
 
  // simulate the number of jumps, including fictitious ones
 
  double rU = as<double>(runif(1));
  double cum = 0;  
 
  if (startState==endState){
    cum = ::Rf_dpois(0,poissonRate*elapsedTime,0)/transProb;
  }
 
  bool notExceed = true;
 
  if (cum > rU){ 
    notExceed = false; 
  }
 
  int numJumps = 0;
  double nextProb;  
 
  while(notExceed){
 
    numJumps++;
    powerTransProb.insert_slices(numJumps,1,false);
    powerTransProb.slice(numJumps) = powerTransProb.slice(numJumps-1)*dominTransProb;
    nextProb = ::Rf_dpois(numJumps,poissonRate*elapsedTime,0)*powerTransProb(startState,endState,numJumps)/transProb;
    cum += nextProb;
 
    if (cum > rU){
      notExceed = false;
    }
  }
 
  List simPath = List::create(Named("states") = NULL,
       Named("times") = NULL);
 
  // if numJumps = 0: done
  if (numJumps == 0 || ((numJumps == 1) && (startState == endState))){
 
    simPath["states"] = NumericVector::create(startState, endState);
    simPath["times"] = NumericVector::create(0, elapsedTime);
 
  }else{ // if one true jump: done
    if ((numJumps == 1) && (startState != endState)){
 
      simPath["states"] = NumericVector::create(startState,endState,endState);
      simPath["times"] = NumericVector::create(0,elapsedTime*(as<double>(runif(1))),elapsedTime);
 
    }else{ // Case (nJmp >= 2) 
 
      // Simulate jumping times
      NumericVector dominJumpTimes = elapsedTime*runif(numJumps);
      std::sort(dominJumpTimes.begin(),dominJumpTimes.end());
      NumericVector dominStates(numJumps+1);
      dominStates[0] = startState;
      dominStates[numJumps] = endState;
 
      // Simulate states of the dominating Markov chain
      for (int i = 1; i < numJumps; i++){
        dominStates[i] = sampleOnce(powerTransProb.slice(1).row(dominStates[i-1]).t()%powerTransProb.slice(numJumps-i).col(endState), as<double>(runif(1)));
      }
 
      // Remove virtual substitutions
      std::vector<int> trueStates(1);
      std::vector<double> trueTimes(1);
      trueStates[0] = startState;
      trueTimes[0] = 0.0;
 
      for (int i = 1; i < dominStates.size(); i++){
        if (dominStates[i-1] != dominStates[i]){
          trueStates.push_back(dominStates[i]);
          trueTimes.push_back(dominJumpTimes[i-1]);
        }
      }
 
      trueStates.push_back(endState);
      trueTimes.push_back(elapsedTime);
 
      simPath["states"] = NumericVector(trueStates.begin(), trueStates.end());
      simPath["times"] = NumericVector(trueTimes.begin(), trueTimes.end());
    }
  }
 
 return simPath;
'

Created by Pretty R at inside-R.org

Now, we compile it with inline:

require(inline)
require(RcppArmadillo)
 
 
unifSample = cxxfunction(signature(rm="numeric", bS="numeric", eS="numeric", eT="numeric", tP="numeric"), includes=inc,body=src, plugin="RcppArmadillo")

Created by Pretty R at inside-R.org

The arguments of unifSample are the CTMC rate matrix (rm), beginning state (bS), ending state (es), time interval length (eT), and transition probability corresponding to the beginning and ending state (tP). The last argument is not strictly required, but I am operating under the assumption that by the time I will want to call this function, I will have transition probabilities pre-computed.

I am going to use markovjumps package to create a simple nucleotide substitution model and simulate from it:

require(markovjumps)
 
my.start = 1
my.end = 2
elapsed.time = 10.0
 
 
my.stat = c(0.1,0.3,0.25,0.35)
 
## define a K80 markov chain (HKY reduces to K80 when stationary distribution is uniform)
my.model = hky.mc(c(0.1,0.3),my.stat,F)
 
## compute eigen-decomposition of the generator
my.model.eigen = as.eigen(my.model)
 
## get transition probabilities
my.prob = matexp(my.model.eigen,elapsed.time)
 
## simulate CTMC path
unifSample(my.model$rate.matrix, my.start-1, my.end-1, elapsed.time, my.prob[my.start,my.end])

Created by Pretty R at inside-R.org

The output looks something like this:

> unifSample(my.model$rate.matrix, my.start-1, my.end-1, elapsed.time, my.prob[my.start,my.end])
$states
[1] 0 3 1 1
 
$times
[1]  0.000000  4.430341  9.295083 10.000000

Created by Pretty R at inside-R.org

Notice that in R my state space is {1,2,3,4}, while in C++ it is {0,1,2,3}. To see more transitions, we can increase the time interval:

> elapsed.time = 100.0
> my.prob = matexp(my.model.eigen,elapsed.time)
> unifSample(my.model$rate.matrix, my.start-1, my.end-1, elapsed.time, my.prob[my.start,my.end])
$states
 [1] 0 3 1 3 0 3 1 3 0 1 3 1 2 3 1 1
 
$times
 [1]   0.000000   4.299258   5.135858   7.616836  14.274938  19.912332  27.649715  35.538612
 [9]  52.687339  53.404439  83.587054  86.185838  88.686555  91.954041  98.642818 100.000000

Created by Pretty R at inside-R.org

To test the simulator, I am going to use analytic formulae for CTMC path statistics, implemented in markovjumps. Namely, I am going to look at the total number of transitions and for total time spent in state 1 (2 in the R state space):

> sim.size = 10000
> num.jumps = array(dim=c(4,4,sim.size))
> dwell.in.one = array(dim=c(4,4,sim.size))
> 
> set.seed(1111)
> for (i in 1:sim.size){
+   for (my.start in 1:4){
+     for (my.end in 1:4){
+       temp1 = unifSample(my.model$rate.matrix, my.start-1, my.end-1, elapsed.time, my.prob[my.start,my.end])
+   
+       num.jumps[my.start,my.end,i] = length(temp1$states)-2
+       dwell.in.one[my.start,my.end,i] = sum(diff(temp1$times)[temp1$states[1:(length(temp1$states)-1)]==1])
+     }
+   }
+ }
> 
> apply(num.jumps,MARGIN=c(1,2),FUN=mean)
        [,1]    [,2]    [,3]    [,4]
[1,] 17.0541 16.9925 16.8286 16.8386
[2,] 16.9071 16.8813 16.7615 16.7366
[3,] 16.9179 16.7869 16.6631 16.6248
[4,] 16.8299 16.7299 16.5572 16.5783
> cond.mean.markov.jumps(my.model.eigen, matrix(1,4,4)-diag(rep(1,4)),elapsed.time)
         [,1]     [,2]     [,3]     [,4]
[1,] 17.06970 16.97879 16.88114 16.82559
[2,] 16.97879 16.88788 16.79024 16.73468
[3,] 16.88114 16.79024 16.69259 16.63704
[4,] 16.82559 16.73468 16.63704 16.58148
> 
> apply(dwell.in.one, MARGIN=c(1,2),FUN=mean)
         [,1]     [,2]     [,3]     [,4]
[1,] 26.20608 30.72730 27.08764 27.11346
[2,] 30.64288 35.35969 31.76749 31.54262
[3,] 27.00934 31.44582 28.03639 27.87922
[4,] 27.16154 31.64142 28.04596 28.07696
> cond.mean.markov.rewards(my.model.eigen, c(0,1,0,0), elapsed.time)
         [,1]     [,2]     [,3]     [,4]
[1,] 26.18182 30.72727 27.09091 27.09091
[2,] 30.72727 35.27273 31.63636 31.63636
[3,] 27.09091 31.63636 28.00000 28.00000
[4,] 27.09091 31.63636 28.00000 28.00000
> 

Created by Pretty R at inside-R.org

Seems to work. I noticed that I can break the program by jacking up the time interval length. I think the problem is in that Poisson density returns non-numbers in some cases and confuses Rcpp. I am not sure if this is worth fixing though, because in practice such parameter configuration shouldn't occur...