[Rcpp-commits] r3625 - in pkg/RcppArmadillo: . inst inst/examples

[noreply at r-forge.r-project.org]

Author: edd
Date: 2012-05-26 20:12:56 +0200 (Sat, 26 May 2012)
New Revision: 3625

Modified:
   pkg/RcppArmadillo/ChangeLog
   pkg/RcppArmadillo/inst/NEWS
   pkg/RcppArmadillo/inst/examples/fastLm.r
Log:
small enhancement to fastLm also comparing 'one vs two' casts


Modified: pkg/RcppArmadillo/ChangeLog
===================================================================
--- pkg/RcppArmadillo/ChangeLog	2012-05-25 18:31:21 UTC (rev 3624)
+++ pkg/RcppArmadillo/ChangeLog	2012-05-26 18:12:56 UTC (rev 3625)
@@ -1,3 +1,8 @@
+2012-05-26  Dirk Eddelbuettel  <edd at debian.org>
+
+	* inst/examples/fastLm.r: Also examine the 'two versus one' casts
+	impact when converting SEXP to Armadillo types
+
 2012-05-21  Dirk Eddelbuettel  <edd at debian.org>
 
 	* DESCRIPTION: Release 0.3.2.0

Modified: pkg/RcppArmadillo/inst/NEWS
===================================================================
--- pkg/RcppArmadillo/inst/NEWS	2012-05-25 18:31:21 UTC (rev 3624)
+++ pkg/RcppArmadillo/inst/NEWS	2012-05-26 18:12:56 UTC (rev 3625)
@@ -1,3 +1,7 @@
+0.3.x.y  2012-xx-yy
+
+    o   Small enhancement to fastLm 
+
 0.3.2.0  2012-05-21
 
     o   Upgraded to Armadillo release 3.2.0 "Creamfields"

Modified: pkg/RcppArmadillo/inst/examples/fastLm.r
===================================================================
--- pkg/RcppArmadillo/inst/examples/fastLm.r	2012-05-25 18:31:21 UTC (rev 3624)
+++ pkg/RcppArmadillo/inst/examples/fastLm.r	2012-05-26 18:12:56 UTC (rev 3625)
@@ -45,10 +45,34 @@
                               Rcpp::Named("df")          =df);
 '
 
-fLm <- cxxfunction(signature(Xs="numeric", ys="numeric"),
-                   src, plugin="RcppArmadillo")
+fLmTwoCasts <- cxxfunction(signature(Xs="numeric", ys="numeric"),
+                           src, plugin="RcppArmadillo")
 
-fastLmPure2 <- function(X, y) {
+src <- '
+    arma::mat X = Rcpp::as<arma::mat>(Xs);
+    arma::colvec y = Rcpp::as<arma::colvec>(ys);
+    int df = X.n_rows - X.n_cols;
+
+    // fit model y ~ X, extract residuals
+    arma::colvec coef = arma::solve(X, y);
+    arma::colvec res  = y - X*coef;
+
+    double s2 = std::inner_product(res.begin(), res.end(),
+                                   res.begin(), 0.0)/df;
+    // std.errors of coefficients
+    arma::colvec sderr = arma::sqrt(s2 *
+       arma::diagvec(arma::pinv(arma::trans(X)*X)));
+
+    return Rcpp::List::create(Rcpp::Named("coefficients")=coef,
+                              Rcpp::Named("stderr")      =sderr,
+                              Rcpp::Named("df")          =df);
+'
+
+fLmOneCast <- cxxfunction(signature(Xs="numeric", ys="numeric"),
+                          src, plugin="RcppArmadillo")
+
+
+fastLmPureDotCall <- function(X, y) {
     .Call("fastLm", X, y, package = "RcppArmadillo")
 }
 
@@ -57,15 +81,16 @@
 X <- cbind(1, log(trees$Girth))
 frm <- formula(log(Volume) ~ log(Girth))
 
-res <- benchmark(fLm(X, y),             	# inline'd above
+res <- benchmark(fLmOneCast(X, y),             	# inline'd above
+                 fLmTwoCasts(X, y),            	# inline'd above
                  fastLmPure(X, y),              # similar, but with 2 error checks
-                 fastLmPure2(X, y),             # now with the 2 error checks
+                 fastLmPureDotCall(X, y),       # now without the 2 error checks
                  fastLm(frm, data=trees),       # using model matrix
                  lm.fit(X, y),                  # R's fast function, no stderr
                  lm(frm, data=trees),           # R's standard function
                  columns = c("test", "replications", "relative",
                              "elapsed", "user.self", "sys.self"),
                  order="relative",
-                 replications=2500)
+                 replications=5000)
 
 print(res)