Using RcppArmadillo for a matrix population model

Standard

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”)

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”),

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!

Midwest Mathematical Biology Conference

Standard

Two weekends ago, I was lucky enough to attend the Midwest Mathematical Biology Conference.

This meeting was hosted by the University of Wisconsin-La Crosse and was the first time for the event. Overall, the meeting was great. I enjoyed the great keynote speakers. Meetings like this help inspire me because I enjoy seeing how other people solve problems. I also enjoy the networking and collaborations formed between researchers.

The keynote address that stood out most to me as the one by Carlos Castillo-Chavez. Carlos provided a great overview on epidemiology modeling and to use multiple scales of modeling together. He summarized a project that used an Individual-based Model to examine how disease spread at the city-level and then scaled up the results to a national level. These results were also compared traditional SIR models. The research does a great job illuminating where people catch infection disease during outbreaks. Surprisingly, most people get sick at home. This happens because people have the most contact with others when they are at home. However, when school is not in session, outbreaks are less likely to occur.

Besides the other great talks, my postdoc mentor had a great observation. He realized the USGS has very few mathematicians and all of them are doing geology (or at least non-biology) research. We wished we had more mathematical support. Furthermore, the Midwest Mathematical Biology Meeting expanded his view of biomath. Even if I gained nothing else from the meeting, having my postdoc mentor appreciate math more made it all worth while! However, I did gain more. So, the meeting was great opportunity!