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))); }

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 >

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

## No comments:

## Post a Comment