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 my understanding that the advantage of the arma::mat is that is
has more methods associated with it, is this correct?
I've also added an R function to do this and as per usual, the speed gain of
simply using Rcpp is incredible!
Any input massively appreciated
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
