[Rcpp-commits] r1226 - in pkg/RcppGSL: inst/unitTests src

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

Author: romain
Date: 2010-05-13 14:53:04 +0200 (Thu, 13 May 2010)
New Revision: 1226

Added:
   pkg/RcppGSL/inst/unitTests/runit.fastLm.R
Modified:
   pkg/RcppGSL/src/fastLm.cpp
Log:
update fastLm to use RcppGSL classes

Added: pkg/RcppGSL/inst/unitTests/runit.fastLm.R
===================================================================
--- pkg/RcppGSL/inst/unitTests/runit.fastLm.R	                        (rev 0)
+++ pkg/RcppGSL/inst/unitTests/runit.fastLm.R	2010-05-13 12:53:04 UTC (rev 1226)
@@ -0,0 +1,34 @@
+#!/usr/bin/r -t
+#                        Emacs make this -*- mode: R; tab-width: 4 -*-
+#
+# Copyright (C) 2010	Romain Francois and Dirk Eddelbuettel
+#
+# This file is part of RcppGSL.
+#
+# RcppGSL 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.
+#
+# RcppGSL 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 RcppGSL.  If not, see <http://www.gnu.org/licenses/>.
+
+test.fastLm <- function() {
+    data(trees)
+    flm <- .Call("fastLm",
+                 log(trees$Volume),
+                 cbind(rep(1,31), log(trees$Girth)),
+                 PACKAGE="RcppGSL")
+    fit <- lm(log(Volume) ~ log(Girth), data=trees)
+
+    checkEquals( as.numeric(flm$coef), as.numeric(coef(fit)),
+                msg="fastLm.coef")
+    checkEquals( as.numeric(flm$stderr), as.numeric(coef(summary(fit))[,2]),
+                msg="fastLm.stderr")
+}
+

Modified: pkg/RcppGSL/src/fastLm.cpp
===================================================================
--- pkg/RcppGSL/src/fastLm.cpp	2010-05-13 12:21:22 UTC (rev 1225)
+++ pkg/RcppGSL/src/fastLm.cpp	2010-05-13 12:53:04 UTC (rev 1226)
@@ -22,42 +22,40 @@
 #include <RcppGSL.h>
 
 #include <gsl/gsl_multifit.h>
+#include <cmath>
 
+RCPP_FUNCTION_2( Rcpp::List, fastLm, SEXP ys, SEXP Xs ){
+	
+	RcppGSL::vector<double> y = ys ;
+	RcppGSL::matrix<double> X = Xs ; 
 
-extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
+	using Rcpp::_ ;
 
-    RcppGSL::vector<double> y(ys);		// creates RcppGSL vector from SEXP
-    RcppGSL::matrix<double> X(Xs);		// creates RcppGSL matrix from SEXP
-
-    //Rcpp::NumericVector Yr(ys);			
-    //Rcpp::NumericMatrix Xr(Xs);			// creates Rcpp matrix from SEXP
-
-    //int i, j, n = Xr.nrow(), k = Xr.ncol();
-    int i, j, n=4,k=2;  // FIXME -- need coef extraction
+    int n=X->size1,k=X->size2;
     double chisq;
 
-    //gsl_matrix *X = gsl_matrix_alloc (n, k);
-    //gsl_vector *y = gsl_vector_alloc (n);
-    gsl_vector *c = gsl_vector_alloc (k);   // FIXME trivial ctor
-    gsl_matrix *cov = gsl_matrix_alloc (k, k);
-    //for (i = 0; i < n; i++) {
-    //    for (j = 0; j < k; j++)
-    //        gsl_matrix_set (X, i, j, Xr(i,j));
-    //    gsl_vector_set (y, i, Yr(i));
-    //}
-
+    RcppGSL::vector<double> c(k) ;
+    RcppGSL::matrix<double> cov(k,k) ;
+    
     gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k);
     gsl_multifit_linear (X, y, c, cov, &chisq, work);
     gsl_multifit_linear_free (work);
-
-	// FIXME need free()
-
-	//Rcpp::NumericVector coefr = Rcpp::wrap(c);
-    //RcppGSL::vector<double> coefr(5);
-    return Rcpp::List::create( Rcpp::Named( "coef", c),
-                               Rcpp::Named( "stderr", c));
-    //return Rcpp::List::create( Rcpp::_["coef"] = c,
-    //                           Rcpp::_["stderr"], c);
-    return Xs;
+         
+    gsl_vector_view diag = gsl_matrix_diagonal(cov) ;
+          
+    Rcpp::NumericVector stderr ; stderr = diag ;
+    std::transform( stderr.begin(), stderr.end(), stderr.begin(), sqrt ) ;
+    
+    Rcpp::List res = Rcpp::List::create( 
+    	_["coef"] = c, 
+    	_["stderr"] = stderr
+    	) ;
+    
+    c.free() ;
+    cov.free();
+    y.free();
+    X.free();
+	
+    return res ;
 }