Hello all,<br><br>I've recently started looking into the speed gains from using Rcpp and have began looking into the arma library too. Ive written a small loop to fill in a matrix and done so using the arma::mat class and the Rcpp::NumericMatrix class. In this case, the Rcpp::NumericMatrix class was slightly faster. Is this usually the case, or have I written quite inefficient code in some parts? It's my understanding that the advantage of the arma::mat is that is has more methods associated with it, is this correct?<br>
<br>I've also added an R function to do this and as per usual, the speed gain of simply using Rcpp is incredible!<br><br>Any input massively appreciated<br><br>Jeff<br><br>library(RcppArmadillo)<br>library(inline)<br>
library(rbenchmark)<br><br># start with R function....<br><br>fillMiddleTheta <- function(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem) {<br> alphaCounter <- alphaCounterStart<br>
for (i in (thetaStartIndex + 1L):(thetaEndIndex - 1L)) {<br> addition <- exp(alphas[alphaCounter] + alphas[alphaCounter + 1L] * timerem)<br> thetaMatrix[ , i] <- thetaMatrix[ , i - 1L] + addition<br>
alphaCounter <- alphaCounter + 2L<br> }<br> return(thetaMatrix)<br>}<br><br># cpp function....<br><br>src <- '<br> Rcpp::NumericMatrix thetaMatrix(sthetaMatrix);<br> <br> Rcpp::NumericVector alphas(salphas);<br>
Rcpp::NumericVector timerem(stimerem);<br> <br> int thetaStartIndex = as<int>(sthetaStartIndex);<br> int thetaEndIndex = as<int>(sthetaEndIndex);<br> int alphaCounterStart = as<int>(salphaCounterStart);<br>
<br> int alphaCounter = alphaCounterStart;<br> <br> for (int i = thetaStartIndex + 1; i <= thetaEndIndex - 1; i++) {<br> <br> thetaMatrix(_, i - 1) = thetaMatrix(_, i - 2) + exp(alphas(alphaCounter - 1) + alphas(alphaCounter) * timerem);<br>
alphaCounter = alphaCounter + 2;<br> <br> }<br> <br> return wrap(thetaMatrix);<br> '<br><br>cppFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix", sthetaStartIndex = "integer", <br>
sthetaEndIndex = "integer", salphas = "numericVector", salphaCounterStart = "integer", stimerem = "numericVector"), <br> src, "Rcpp")<br><br># cpp function with armadillo....<br>
<br>armaSrc <- '<br> arma::mat thetaMatrix = Rcpp::as<arma::mat>(sthetaMatrix);<br> <br> arma::colvec alphas = Rcpp::as<arma::colvec>(salphas);<br> arma::colvec timerem = Rcpp::as<arma::colvec>(stimerem);<br>
<br> int thetaStartIndex = as<int>(sthetaStartIndex);<br> int thetaEndIndex = as<int>(sthetaEndIndex);<br> int alphaCounterStart = as<int>(salphaCounterStart);<br> <br>
int alphaCounter = alphaCounterStart;<br> <br> for (int i = thetaStartIndex + 1; i <= thetaEndIndex - 1; i++) {<br> <br> thetaMatrix.col(i - 1) = thetaMatrix.col(i - 2) + exp(alphas(alphaCounter - 1) + alphas(alphaCounter) * timerem);<br>
alphaCounter = alphaCounter + 2;<br> <br> }<br> <br> return wrap(thetaMatrix);<br> '<br><br>cppArmaFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix", sthetaStartIndex = "integer", <br>
sthetaEndIndex = "integer", salphas = "numericVector", salphaCounterStart = "integer", stimerem = "numericVector"), <br> armaSrc, "RcppArmadillo")<br>
<br># test<br>n <- 1e6<br><br>thetaMatrix <- matrix(0, ncol = 10, nrow = 10)<br>thetaMatrix[, 1L] <- 1<br>thetaMatrix[, 10L] <- 100<br>thetaStartIndex <- 1L<br>thetaEndIndex <- 10L<br>alphas <- seq(0, 0.06, length = 16)<br>
alphaCounterStart <- 1L<br>timerem <- rnorm(10, mean = 40, sd = 2)<br><br>r <- fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem)<br>cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem)<br>
cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem)<br>all.equal(r, cpp, cppArma)<br><br>benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br> cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
<br> columns = c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),<br> order = "relative", <br> replications = n)<br>
<br>Here is some of my output;<br><br>> n <- 1e6<br>> r <- fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem)<br>> cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem)<br>
> cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem)<br>> all.equal(r, cpp, cppArma)<br>[1] TRUE<br><br>> benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
+ cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>+ cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
+ <br>+ columns = c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),<br>+ order = "relative",<br>+ replications = n)<br>
> benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>+ cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
+ cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem),<br>+ <br>+ columns = c("test", "replications", "elapsed", "relative", "user.self", "sys.self"),<br>
+ order = "relative",<br>+ replications = n)<br><br> test replications elapsed relative user.self sys.self<br>
2 cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem) 1000000 12.22 1.000000 12.17 0.01<br>3 cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem) 1000000 15.07 1.233224 15.01 0.02<br>
1 fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas, alphaCounterStart, timerem) 1000000 86.62 7.088380 86.02 0.10<br>