Update benchmark script

- Use `microbenchmark` instead of `rbenchmark`

  Report relative and absolute timings

- Use `fastLmPure` functions instead of calling
  C++ directly

Author: michaelweylandt

ChangeLog

@@ -1,7 +1,8 @@
 2018-05-30  Michael Weylandt  <michael.weylandt@gmail.com>
 
-	* inst/examples/lmBenchmark.R: Update call to RcppArmadillo's version
-	of fastLm in benchmark script
+	* inst/examples/lmBenchmark.R: Update benchmark script to use
+	microbenchmark and to use exposed fastLm functions from Rcpp
+	packages rather than invoking .Call directly
 
 2018-05-25  Ralf Stubner  <ralf.stubner@daqana.com>
 

inst/examples/lmBenchmark.R

@@ -5,7 +5,7 @@
 ## This file is part of RcppEigen.
 
 require("stats", character=TRUE, quietly=TRUE)
-require("rbenchmark", character=TRUE, quietly=TRUE)
+require("microbenchmark", character=TRUE, quietly=TRUE)
 require("RcppEigen", character=TRUE, quietly=TRUE)
 
 ## define different versions of lm
@@ -15,45 +15,61 @@ exprs <- list()
 ## handle rank-deficient cases.
 
                                         # default version used in lm()
-exprs$lm.fit <- expression(stats::lm.fit(mm, y))
+exprs["lm.fit"] <- alist(stats::lm.fit(mm, y))
                                         # versions from RcppEigen
 ## column-pivoted QR decomposition - similar to lm.fit
-exprs$PivQR <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 0L, PACKAGE="RcppEigen"))
+exprs["PivQR"] <- alist(RcppEigen::fastLmPure(mm, y, 0L))
 ## LDLt Cholesky decomposition with rank detection
-exprs$LDLt <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 2L, PACKAGE="RcppEigen"))
+exprs["LDLt"] <- alist(RcppEigen::fastLmPure(mm, y, 2L))
 ## SVD using the Lapack subroutine dgesdd and Eigen support
-exprs$GESDD <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 6L, PACKAGE="RcppEigen"))
+exprs["GESDD"] <- alist(RcppEigen::fastLmPure(mm, y, 6L))
 ## SVD (the JacobiSVD class from Eigen)
-exprs$SVD <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 4L, PACKAGE="RcppEigen"))
+exprs["SVD"] <- alist(RcppEigen::fastLmPure(mm, y, 4L))
 ## eigenvalues and eigenvectors of X'X
-exprs$SymmEig <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 5L, PACKAGE="RcppEigen"))
+exprs["SymmEig"] <- alist(RcppEigen::fastLmPure(mm, y, 5L))
 
 ## Non-rank-revealing decompositions.  These work fine except when
 ## they don't.
 
 ## Unpivoted  QR decomposition
-exprs$QR <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 1L, PACKAGE="RcppEigen"))
+exprs["QR"] <- alist(RcppEigen::fastLmPure(mm, y, 1L))
 ## LLt Cholesky decomposition
-exprs$LLt <- expression(.Call("RcppEigen_fastLm_Impl", mm, y, 3L, PACKAGE="RcppEigen"))
+exprs["LLt"] <- alist(RcppEigen::fastLmPure(mm, y, 3L))
 
 if (suppressMessages(require("RcppArmadillo", character=TRUE, quietly=TRUE))) {
-    exprs$arma <- expression(.Call("_RcppArmadillo_fastLm_impl", mm, y, PACKAGE="RcppArmadillo"))
+    exprs["arma"] <- alist(RcppArmadillo::fastLmPure(mm, y))
 }
 
 if (suppressMessages(require("RcppGSL", character=TRUE, quietly=TRUE))) {
-    exprs$GSL <- expression(.Call("RcppGSL_fastLm", mm, y, PACKAGE="RcppGSL"))
+    exprs["GSL"] <- alist(RcppGSL::fastLmPure(mm, y))
 }
 
 do_bench <- function(n=100000L, p=40L, nrep=20L, suppressSVD=(n > 100000L)) {
     mm <- cbind(1, matrix(rnorm(n * (p - 1L)), nc=p-1L))
     y <- rnorm(n)
     if (suppressSVD) exprs <- exprs[!names(exprs) %in% c("SVD", "GSL")]
     cat("lm benchmark for n = ", n, " and p = ", p, ": nrep = ", nrep, "\n", sep='')
-    do.call(benchmark, c(exprs,
-                         list(order="relative",
-                              columns = c("test", "relative",
-                              "elapsed", "user.self", "sys.self"),
-                              replications = nrep)))
+    mb <- microbenchmark(list=exprs, times = nrep)
+    
+    op <- options(microbenchmark.unit="relative")
+    on.exit(options(op))
+    
+    mb_relative <- summary(mb)
+    levels(mb_relative$expr) <- names(exprs)
+    
+    options(microbenchmark.unit=NULL)
+    mb_absolute <- summary(mb)
+    levels(mb_absolute$expr) <- names(exprs)
+
+    mb_combined <- merge(mb_relative[, c("expr", "median")], 
+                         mb_absolute[, c("expr", "median")], 
+                         by="expr")
+    
+    colnames(mb_combined) <- c("Method", 
+                               "Relative", 
+                               paste0("Elapsed (", attr(mb_absolute, "unit"), ")"))
+    
+    mb_combined[order(mb_combined$Relative),]
 }
 
 print(do_bench())