Hey Dirk,
Thanks for the informative and speedy reply! Having a look through the
archives now....
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.
Thanks again,
Jeff
On Fri, Sep 16, 2011 at 3:34 PM, Dirk Eddelbuettel <e...@debian.org> wrote:
>
> Hi Jeffrey,
>
> And welcome! Thanks for a detailed post with examples and benchmarks!
>
> On 16 September 2011 at 14:48, Jeffrey Pollock wrote:
> | Hello all,
> |
> | 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
>
> This may not be worth worrying about. See the list archives just two or
> three weeks back: Darren Cook noticed that even within Rcpp, the times
> between accessing vector (not matrix) elements using () and [] for indexing
> differed far more (which may well be indicative of us doing something
> silly,
> there may be a cache lookup gone wrong...)
>
> For all real applications, _actual work_ is bound to dominate. So as the
> difference is small... Also, check the Arma docs. Conrad is very generous
> is his API, there are often several different ways to do similar things
>
> | my understanding that the advantage of the arma::mat is that is has more
> | methods associated with it, is this correct?
>
> Armadillo is a full-blown linear algebra library. If you want to do 'math'
> with your matrix, it is a good choice to deploy arma matrices
>
> Rcpp::NumericMatrix is a very convenient container to get matrix data back
> and from R, and its vector variant also goes easily to STL types etc so
> that
> you can (relatively) easily exchange data with other C++ libraries etc.
>
> So it depends.
>
> | I've also added an R function to do this and as per usual, the speed gain
> of
> | simply using Rcpp is incredible!
>
> Yes, compiled loops often win rather handily.
>
> | Any input massively appreciated
>
> Hope this helps, Dirk
>
> | Jeff
> |
> | library(RcppArmadillo)
> | library(inline)
> | library(rbenchmark)
> |
> | # start with R function....
> |
> | fillMiddleTheta <- function(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem) {
> | alphaCounter <- alphaCounterStart
> | for (i in (thetaStartIndex + 1L):(thetaEndIndex - 1L)) {
> | addition <- exp(alphas[alphaCounter] + alphas[alphaCounter + 1L]
> | * timerem)
> | thetaMatrix[ , i] <- thetaMatrix[ , i - 1L] + addition
> | alphaCounter <- alphaCounter + 2L
> | }
> | return(thetaMatrix)
> | }
> |
> | # cpp function....
> |
> | src <- '
> | Rcpp::NumericMatrix thetaMatrix(sthetaMatrix);
> |
> | Rcpp::NumericVector alphas(salphas);
> | Rcpp::NumericVector timerem(stimerem);
> |
> | int thetaStartIndex = as<int>(sthetaStartIndex);
> | int thetaEndIndex = as<int>(sthetaEndIndex);
> | int alphaCounterStart = as<int>(salphaCounterStart);
> |
> | int alphaCounter = alphaCounterStart;
> |
> | for (int i = thetaStartIndex + 1; i <= thetaEndIndex - 1; i++) {
> |
> | thetaMatrix(_, i - 1) = thetaMatrix(_, i - 2) + exp(alphas
> | (alphaCounter - 1) + alphas(alphaCounter) * timerem);
> | alphaCounter = alphaCounter + 2;
> |
> | }
> |
> | return wrap(thetaMatrix);
> | '
> |
> | cppFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix",
> | sthetaStartIndex = "integer",
> | sthetaEndIndex = "integer", salphas =
> | "numericVector", salphaCounterStart = "integer", stimerem =
> "numericVector"),
> | src, "Rcpp")
> |
> | # cpp function with armadillo....
> |
> | armaSrc <- '
> | arma::mat thetaMatrix = Rcpp::as<arma::mat>(sthetaMatrix);
> |
> | arma::colvec alphas = Rcpp::as<arma::colvec>(salphas);
> | arma::colvec timerem = Rcpp::as<arma::colvec>(stimerem);
> |
> | int thetaStartIndex = as<int>(sthetaStartIndex);
> | int thetaEndIndex = as<int>(sthetaEndIndex);
> | int alphaCounterStart = as<int>(salphaCounterStart);
> |
> | int alphaCounter = alphaCounterStart;
> |
> | for (int i = thetaStartIndex + 1; i <= thetaEndIndex - 1; i++) {
> |
> | thetaMatrix.col(i - 1) = thetaMatrix.col(i - 2) + exp(alphas
> | (alphaCounter - 1) + alphas(alphaCounter) * timerem);
> | alphaCounter = alphaCounter + 2;
> |
> | }
> |
> | return wrap(thetaMatrix);
> | '
> |
> | cppArmaFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix",
> | sthetaStartIndex = "integer",
> | sthetaEndIndex = "integer", salphas =
> | "numericVector", salphaCounterStart = "integer", stimerem =
> "numericVector"),
> | armaSrc, "RcppArmadillo")
> |
> | # test
> | n <- 1e6
> |
> | thetaMatrix <- matrix(0, ncol = 10, nrow = 10)
> | thetaMatrix[, 1L] <- 1
> | thetaMatrix[, 10L] <- 100
> | thetaStartIndex <- 1L
> | thetaEndIndex <- 10L
> | alphas <- seq(0, 0.06, length = 16)
> | alphaCounterStart <- 1L
> | timerem <- rnorm(10, mean = 40, sd = 2)
> |
> | r <- fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,
> | alphaCounterStart, timerem)
> | cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> alphas,
> | alphaCounterStart, timerem)
> | cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
> thetaEndIndex,
> | alphas, alphaCounterStart, timerem)
> | all.equal(r, cpp, cppArma)
> |
> | benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> alphas,
> | alphaCounterStart, timerem),
> | cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem),
> | cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
> | thetaEndIndex, alphas, alphaCounterStart, timerem),
> |
> | columns = c("test", "replications", "elapsed", "relative",
> | "user.self", "sys.self"),
> | order = "relative",
> | replications = n)
> |
> | Here is some of my output;
> |
> | > n <- 1e6
> | > r <- fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> alphas,
> | alphaCounterStart, timerem)
> | > cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem)
> | > cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
> | thetaEndIndex, alphas, alphaCounterStart, timerem)
> | > all.equal(r, cpp, cppArma)
> | [1] TRUE
> |
> | > benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem),
> | + cppFillMiddleTheta(thetaMatrix, thetaStartIndex,
> | thetaEndIndex, alphas, alphaCounterStart, timerem),
> | + cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
> | thetaEndIndex, alphas, alphaCounterStart, timerem),
> | +
> | + columns = c("test", "replications", "elapsed", "relative",
> | "user.self", "sys.self"),
> | + order = "relative",
> | + replications = n)
> | > benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem),
> | + cppFillMiddleTheta(thetaMatrix, thetaStartIndex,
> | thetaEndIndex, alphas, alphaCounterStart, timerem),
> | + cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
> | thetaEndIndex, alphas, alphaCounterStart, timerem),
> | +
> | + columns = c("test", "replications", "elapsed", "relative",
> | "user.self", "sys.self"),
> | + order = "relative",
> | + replications = n)
> |
> |
>
> | test replications elapsed relative user.self sys.self
> | 2 cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem) 1000000 12.22
> | 1.000000 12.17 0.01
> | 3 cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> alphas,
> | alphaCounterStart, timerem) 1000000 15.07 1.233224
> | 15.01 0.02
> | 1 fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
> | alphas, alphaCounterStart, timerem) 1000000 86.62
> | 7.088380 86.02 0.10
> |
> | ----------------------------------------------------------------------
> | _______________________________________________
> | Rcpp-devel mailing list
> | Rcpp-devel@lists.r-forge.r-project.org
> | https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
> --
> New Rcpp master class for R and C++ integration is scheduled for
> San Francisco (Oct 8), more details / reg.info available at
>
> http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php
>
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
