[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