[Rcpp-commits] r831 - in pkg/RcppArmadillo: . R inst src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Mar 4 04:38:18 CET 2010


Author: edd
Date: 2010-03-04 04:38:14 +0100 (Thu, 04 Mar 2010)
New Revision: 831

Added:
   pkg/RcppArmadillo/src/fastLm.cpp
Modified:
   pkg/RcppArmadillo/NEWS
   pkg/RcppArmadillo/R/fastLm.R
   pkg/RcppArmadillo/inst/ChangeLog
   pkg/RcppArmadillo/src/RcppArmadillo.cpp
Log:
fastLm now in its own file
with vector ctor also 0.9.0-style
added some would-be-nice TODOs in fastLm.R


Modified: pkg/RcppArmadillo/NEWS
===================================================================
--- pkg/RcppArmadillo/NEWS	2010-03-03 15:53:34 UTC (rev 830)
+++ pkg/RcppArmadillo/NEWS	2010-03-04 03:38:14 UTC (rev 831)
@@ -1,5 +1,10 @@
 0.1-0	(under development)
 
+    o	the fastLm() implementation of a bare-bones lm() fit
+	(using Armadillo's solve() function) provides an example
+	of how efficient code can be written compactly using 
+	the combination of Rcpp, RcppAramadillo and Armadillo
+
     o	support for Rcpp implicit wrap of these types :
     	Mat<T>, Col<T>, Row<T>, Cube<T> where T is one of : 
     	int, unsigned int, double, float

Modified: pkg/RcppArmadillo/R/fastLm.R
===================================================================
--- pkg/RcppArmadillo/R/fastLm.R	2010-03-03 15:53:34 UTC (rev 830)
+++ pkg/RcppArmadillo/R/fastLm.R	2010-03-04 03:38:14 UTC (rev 831)
@@ -25,3 +25,18 @@
 
     res <- .Call("fastLm", y, X, package="RcppArmadillo")
 }
+
+## What would be nice here:
+##
+##  fastLm <- function(x, ...) UseMethod("fastLm")
+##
+##  fastLm.formula <- ...
+##
+##  fastLm.default <- ...
+##
+##  print.fastLm <- function(x, ...)
+##
+##  summary.fastLm <- function(object, ...
+##
+##  print.summary.fastLm <- ...
+

Modified: pkg/RcppArmadillo/inst/ChangeLog
===================================================================
--- pkg/RcppArmadillo/inst/ChangeLog	2010-03-03 15:53:34 UTC (rev 830)
+++ pkg/RcppArmadillo/inst/ChangeLog	2010-03-04 03:38:14 UTC (rev 831)
@@ -1,14 +1,18 @@
+2010-03-03  Dirk Eddelbuettel  <edd at debian.org>
+
+	* src/fastLm.cpp: Moved into its own file, some more polish
+
 2010-03-03  Romain Francois <romain at r-enthusiasts.com>
 
-	* inst/include/RcppArmadillo.h: avoid an extra memory copy when 
+	* inst/include/RcppArmadillo.h: avoid an extra memory copy when
 	possible (i.e. in wrap( eGlue) and wrap( eOp ) when the elem_type
 	is int or double).
 
-	* inst/include/RcppArmadilloDefines.in.h added the SCALAR macro 
+	* inst/include/RcppArmadilloDefines.in.h added the SCALAR macro
 	to take care of the 0.9.0 api change
 
 	* src/RcppArmadillo.cpp : fix a runtime error uisng the SCALAR macro
-	
+
 2010-03-02  Dirk Eddelbuettel  <edd at dexter>
 
 	* src/RcppArmadillo.cpp: Added bare-bones 'fastLm' function

Modified: pkg/RcppArmadillo/src/RcppArmadillo.cpp
===================================================================
--- pkg/RcppArmadillo/src/RcppArmadillo.cpp	2010-03-03 15:53:34 UTC (rev 830)
+++ pkg/RcppArmadillo/src/RcppArmadillo.cpp	2010-03-04 03:38:14 UTC (rev 831)
@@ -171,44 +171,14 @@
 }
 
 extern "C" SEXP RcppArmadillo_wrap_Op(){
+
     arma::mat m1 = eye<mat>( 3, 3 ) ;
 	
     List res ;
     res["mat+mat"] = - m1 ;
     return res ;
-	
 }
 
-extern "C" SEXP fastLm(SEXP sy, SEXP sX) {
-
-    Rcpp::NumericVector yr(sy);
-    Rcpp::NumericMatrix Xr(sX);
-    int n = Xr.nrow(), k = Xr.ncol() ;
-
-#if ARMA_VERSION_GE_090
-    arma::mat X(Xr.begin(), n, k, false);   // reuses memory and avoids extra copy
-#else
-    arma::mat X(Xr.begin(), n, k);
-#endif
-
-    arma::colvec y(yr.begin(), yr.size());
-
-    arma::colvec coef = solve(X, y);            // fit model y ~ X
-
-    arma::colvec resid = y - X*coef; 
-
-#if ARMA_VERSION_GE_090
-    double sig2 = arma::as_scalar( trans(resid)*resid/(n-k) );
-#else
-    double sig2 = ( trans(resid)*resid/(n-k) );
-#endif
- 
-    arma::colvec stderrest = sqrt( sig2 * diagvec( arma::inv(arma::trans(X)*X)) );
-    Rcpp::Pairlist res(Rcpp::Named( "coef", coef),
-                       Rcpp::Named( "stderr", stderrest));
-    return res;
-}
-
 extern "C" SEXP armadillo_version(SEXP single_){
     struct arma::arma_version av;
     bool single = Rcpp::as<bool>( single_) ;

Added: pkg/RcppArmadillo/src/fastLm.cpp
===================================================================
--- pkg/RcppArmadillo/src/fastLm.cpp	                        (rev 0)
+++ pkg/RcppArmadillo/src/fastLm.cpp	2010-03-04 03:38:14 UTC (rev 831)
@@ -0,0 +1,54 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
+//
+// fastLm.cpp: Rcpp/Armadillo glue example of a simple lm() alternative
+//
+// Copyright (C)  2010 Dirk Eddelbuettel and Romain Francois
+//
+// This file is part of RcppArmadillo.
+//
+// RcppArmadillo is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// RcppArmadillo is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with RcppArmadillo.  If not, see <http://www.gnu.org/licenses/>.
+
+#include <RcppArmadillo.h>
+
+extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
+
+    Rcpp::NumericVector yr(ys);			// creates Rcpp vector from SEXP
+    Rcpp::NumericMatrix Xr(Xs);			// creates Rcpp matrix from SEXP
+    int n = Xr.nrow(), k = Xr.ncol();
+
+#if ARMA_VERSION_GE_090
+    arma::mat X(Xr.begin(), n, k, false);   	// reuses memory and avoids extra copy
+    arma::colvec y(yr.begin(), yr.size(), false);
+#else
+    arma::mat X(Xr.begin(), n, k);
+    arma::colvec y(yr.begin(), yr.size());
+#endif
+
+    arma::colvec coef = solve(X, y);            // fit model y ~ X
+
+    arma::colvec resid = y - X*coef; 		// residuals
+
+#if ARMA_VERSION_GE_090
+    double sig2 = arma::as_scalar( trans(resid)*resid/(n-k) );
+#else
+    double sig2 = ( trans(resid)*resid/(n-k) );
+#endif
+    						// std.error of estimate 
+    arma::colvec stderrest = sqrt( sig2 * diagvec( arma::inv(arma::trans(X)*X)) );
+
+    Rcpp::Pairlist res(Rcpp::Named( "coef", coef),
+                       Rcpp::Named( "stderr", stderrest));
+    return res;
+}
+



More information about the Rcpp-commits mailing list