[Rcpp-commits] r3089 - pkg/RcppEigen/src

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

Author: dmbates
Date: 2011-06-15 21:14:44 +0200 (Wed, 15 Jun 2011)
New Revision: 3089

Modified:
   pkg/RcppEigen/src/fastLm.cpp
Log:
Initial (not very good) implementation of "fastLmSVD", which isn't really that fast.


Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp	2011-06-15 19:13:45 UTC (rev 3088)
+++ pkg/RcppEigen/src/fastLm.cpp	2011-06-15 19:14:44 UTC (rev 3089)
@@ -181,3 +181,36 @@
     }
     return R_NilValue; // -Wall
 }
+
+extern "C" SEXP fastLmSVD(SEXP Xs, SEXP ys) {
+    using namespace Eigen;
+    using namespace Rcpp;
+    try {
+	NumericMatrix   X(Xs);
+	NumericVector   y(ys);
+	size_t          n = X.nrow(), p = X.ncol();
+	ptrdiff_t      df = n - p;
+	if ((size_t)y.size() != n)
+	    throw std::invalid_argument("size mismatch");
+	
+	MatrixXd        A = Map<MatrixXd>(X.begin(), n, p); // shares storage
+	VectorXd        b = Map<VectorXd>(y.begin(), n);
+
+	VectorXd     coef = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);
+//	double         s2 = (b - A*coef).squaredNorm()/df;
+
+	// ArrayXd        se = (Ch.solve(MatrixXd::Identity(p, p)).diagonal().array() * s2).sqrt();
+	// NumericVector Rse(p);	// should define a wrap method for ArrayXd, ArrayXXd, etc.
+	// std::copy(se.data(), se.data() + p, Rse.begin());
+			    
+	return List::create(_["coefficients"] = coef,
+			    // _["se"]           = Rse,
+			    _["df.residual"]  = df);
+    } catch( std::exception &ex ) {
+	forward_exception_to_r( ex );
+    } catch(...) { 
+	::Rf_error( "c++ exception (unknown reason)" ); 
+    }
+    return R_NilValue; // -Wall
+}
+