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.