Hey Dirk,<br><br>Thanks for the informative and speedy reply! Having a look through the archives now....<br><br>I should hopefully have time to post up quite a few benchmarks of work I'm doing at the moment/as I learn more of Rcpp.<br>
<br>Thanks again,<br><br>Jeff<br><br><div class="gmail_quote">On Fri, Sep 16, 2011 at 3:34 PM, Dirk Eddelbuettel <span dir="ltr"><<a href="mailto:edd@debian.org">edd@debian.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<br>
Hi Jeffrey,<br>
<br>
And welcome! Thanks for a detailed post with examples and benchmarks!<br>
<div class="im"><br>
On 16 September 2011 at 14:48, Jeffrey Pollock wrote:<br>
| Hello all,<br>
|<br>
| I've recently started looking into the speed gains from using Rcpp and have<br>
| began looking into the arma library too. Ive written a small loop to fill in a<br>
| matrix and done so using the arma::mat class and the Rcpp::NumericMatrix class.<br>
| In this case, the Rcpp::NumericMatrix class was slightly faster. Is this<br>
| usually the case, or have I written quite inefficient code in some parts? It's<br>
<br>
</div>This may not be worth worrying about. See the list archives just two or<br>
three weeks back: Darren Cook noticed that even within Rcpp, the times<br>
between accessing vector (not matrix) elements using () and [] for indexing<br>
differed far more (which may well be indicative of us doing something silly,<br>
there may be a cache lookup gone wrong...)<br>
<br>
For all real applications, _actual work_ is bound to dominate. So as the<br>
difference is small... Also, check the Arma docs. Conrad is very generous<br>
is his API, there are often several different ways to do similar things<br>
<div class="im"><br>
| my understanding that the advantage of the arma::mat is that is has more<br>
| methods associated with it, is this correct?<br>
<br>
</div>Armadillo is a full-blown linear algebra library. If you want to do 'math'<br>
with your matrix, it is a good choice to deploy arma matrices<br>
<br>
Rcpp::NumericMatrix is a very convenient container to get matrix data back<br>
and from R, and its vector variant also goes easily to STL types etc so that<br>
you can (relatively) easily exchange data with other C++ libraries etc.<br>
<br>
So it depends.<br>
<div class="im"><br>
| I've also added an R function to do this and as per usual, the speed gain of<br>
| simply using Rcpp is incredible!<br>
<br>
</div>Yes, compiled loops often win rather handily.<br>
<br>
| Any input massively appreciated<br>
<br>
Hope this helps, Dirk<br>
<div><div></div><div class="h5"><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,<br>
| alphas, alphaCounterStart, timerem) {<br>
| alphaCounter <- alphaCounterStart<br>
| for (i in (thetaStartIndex + 1L):(thetaEndIndex - 1L)) {<br>
| addition <- exp(alphas[alphaCounter] + alphas[alphaCounter + 1L]<br>
| * 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<br>
| (alphaCounter - 1) + alphas(alphaCounter) * timerem);<br>
| alphaCounter = alphaCounter + 2;<br>
| <br>
| }<br>
| <br>
| return wrap(thetaMatrix);<br>
| '<br>
|<br>
| cppFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix",<br>
| sthetaStartIndex = "integer",<br>
| sthetaEndIndex = "integer", salphas =<br>
| "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<br>
| (alphaCounter - 1) + alphas(alphaCounter) * timerem);<br>
| alphaCounter = alphaCounter + 2;<br>
| <br>
| }<br>
| <br>
| return wrap(thetaMatrix);<br>
| '<br>
|<br>
| cppArmaFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix",<br>
| sthetaStartIndex = "integer",<br>
| sthetaEndIndex = "integer", salphas =<br>
| "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,<br>
| alphaCounterStart, timerem)<br>
| cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,<br>
| alphaCounterStart, timerem)<br>
| cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem)<br>
| all.equal(r, cpp, cppArma)<br>
|<br>
| benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,<br>
| alphaCounterStart, timerem),<br>
| cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem),<br>
| cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,<br>
| thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
| <br>
| columns = c("test", "replications", "elapsed", "relative",<br>
| "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,<br>
| alphaCounterStart, timerem)<br>
| > cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem)<br>
| > cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,<br>
| thetaEndIndex, alphas, alphaCounterStart, timerem)<br>
| > all.equal(r, cpp, cppArma)<br>
| [1] TRUE<br>
|<br>
| > benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem),<br>
| + cppFillMiddleTheta(thetaMatrix, thetaStartIndex,<br>
| thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
| + cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,<br>
| thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
| + <br>
| + columns = c("test", "replications", "elapsed", "relative",<br>
| "user.self", "sys.self"),<br>
| + order = "relative",<br>
| + replications = n)<br>
| > benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem),<br>
| + cppFillMiddleTheta(thetaMatrix, thetaStartIndex,<br>
| thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
| + cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,<br>
| thetaEndIndex, alphas, alphaCounterStart, timerem),<br>
| + <br>
| + columns = c("test", "replications", "elapsed", "relative",<br>
| "user.self", "sys.self"),<br>
| + order = "relative",<br>
| + replications = n)<br>
|<br>
| <br>
| test replications elapsed relative user.self sys.self<br>
| 2 cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem) 1000000 12.22<br>
| 1.000000 12.17 0.01<br>
| 3 cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,<br>
| alphaCounterStart, timerem) 1000000 15.07 1.233224 <br>
| 15.01 0.02<br>
| 1 fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,<br>
| alphas, alphaCounterStart, timerem) 1000000 86.62<br>
| 7.088380 86.02 0.10<br>
|<br>
</div></div>| ----------------------------------------------------------------------<br>
| _______________________________________________<br>
| Rcpp-devel mailing list<br>
| <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
| <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
<font color="#888888">--<br>
New Rcpp master class for R and C++ integration is scheduled for<br>
San Francisco (Oct 8), more details / <a href="http://reg.info" target="_blank">reg.info</a> available at<br>
<a href="http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php" target="_blank">http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php</a><br>
</font></blockquote></div><br>