I’ve had a busy winter and spring with training for the Birkie and then recovering from the Birkie, hence the few posts. One of the things I’ve been doing is teaching myself the Rcpp package in R. This package lets people easily use C++ code within R. In this post, I demonstrate how I used Rcpp and specifically the RcppArmadillo package to create a population model.

Matrix models are popular in ecology. These models are series of difference equations (i.e., discrete time). I was interested in coding a simple example for a two life-stage species with the projection matrix R.

To code this up in R, I use the following code

popModelR A = matrix(c(0, 2, 0.25, 0.5), nrow = 2, byrow = TRUE),

P0 = c(0, 50)){

P = matrix(0, nrow = dim(A)[1], ncol = nYears + 1)

P[, 1] = P0

for( t in 1:nYears){

P[ , t + 1] = A %*% P[, t]

}

return(P)

}nYears = 10

A = matrix(c(0, 2, 0.25, 0.5),

nrow = 2, byrow = TRUE)

P0 = matrix(c(0, 50), nrow = 2)

P0

popModelROut <- popModelR(nYears = 10,

A = A,

P0 = P0)

Obviously, this simple model runs fast, but how would one code this with Rcpp? In order to get matrix operators, I needed to use the RcppArmadillo package, so my code looks like this:

library(“inline”)

library(“Rcpp”)

library(“RcppArmadillo”)src1 <- ‘

int nYearsX = Rcpp::as<int>(nYears);

arma::mat P0X = Rcpp::as<arma::mat>(P0);

arma::mat AX = Rcpp::as<arma::mat>(A);

arma::mat PX(AX.n_cols, nYearsX + 1);

PX.col(0) = P0X;for(int t = 0; t < nYearsX; t++) {

PX.col(t + 1) = AX * PX.col(t);

}return Rcpp::wrap(PX);

‘popModelRcpp <- cxxfunction(signature(nYears = “integer”,

A = “matrix”,

P0 = “matrix”),

src1, plugin = “RcppArmadillo”)popModelRcpp(nYears, A, P0)

Now, to compare the two functions, I use the benchmark package and run the model for 100 simulated years:

library(rbenchmark)

nYears = 100

res popModelR(nYears, A, P0),

columns = c(“test”, “replications”, “elapsed”,

“relative”, “user.self”, “sys.self”),

order = “relative”,

replications = 1000)

print(res)

test | replications | elapsed | relative | user.self | sys.self | |
---|---|---|---|---|---|---|

1 | popModelRcpp(nYears, A, P0) | 1000 | 0.02 | 1.00 | 0.00 | 0.00 |

2 | popModelR(nYears, A, P0) | 1000 | 0.53 | 29.61 | 0.64 | 0.00 |

The Rcpp code is almost 30 times quicker than the base code in R!