Wed, 01 Dec 2010

RcppGSL 0.1.0

Earlier in the year, Romain and I did a bunch of initial work on a wrapper from R to the GNU GSL by way of our Rcpp package for seamless R and C++ integration. But other work kept us busy and this fell a little to the side.

We have now found some time to finish this work for a first release, together with a nicely detailed eleven page package vignette. As of today, the package is now a CRAN package, and Romain already posted a nice announcement on his blog and on the rcpp-devel list.

So what does RcppGSL do? I gave the package its own webpage here as well and listed these points as key features of RcppGSL:

• templated vector and matrix classes: these are similar to Rcpp's own vector and matrix classes, but really are just smart pointers around the C structure expected by the library
• this means you can transfer data from R to your GSL-using programs in pretty much the same way you would in other C++ programs using Rcpp---by relying on the Rcpp::as() and Rcpp::wrap() converterrs
• at the C++ level you can use these GSL vectors in a more C++-alike way (using eg foo[i] to access an element at index i)
• yet at the same time you can just pass these vector and matrix objects to the GSL functions expecting its C objects: thanks to some cleverness in these classes they pass the right object on (see the example below)
• we also provide the lightweight views for vectors and matrices as the GSL API uses these in many places.

Also provided is a simple example which is a simple implementation of a column norm (which we could easily compute directly in R, but we are simply re-using an example from Section 8.4.14 of the GSL manual):

#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

extern "C" SEXP colNorm(SEXP sM) {

  try {

        RcppGSL::matrix<double> M = sM;     // create gsl data structures from SEXP
        int k = M.ncol();
        Rcpp::NumericVector n(k);           // to store results

        for (int j = 0; j < k; j++) {
            RcppGSL::vector_view<double> colview = gsl_matrix_column (M, j);
            n[j] = gsl_blas_dnrm2(colview);
        }
        M.free() ;
        return n;                           // return vector

  } catch( std::exception &ex ) {
        forward_exception_to_r( ex );

  } catch(...) {
        ::Rf_error( "c++ exception (unknown reason)" );
  }
  return R_NilValue; // -Wall
}

This example function is implemented in an example package contained in the RcppGSL package itself -- so that users have a complete stanza to use in their packages. This will then build a user package on Linux, OS X and Windows provided the GSL is installed (and on Windows you have to do all the extra steps of defining an environment variable pointing to and of course install Rtools to build in the first place---Linux and OS X are so much easier for development).

Another complete example is in the package itself and provides a faster (compiled) alternative to the standard lm() function in R; this example is the continuation of the same example I had in several versions of my Intro to HPC with R tutorials and in the Rcpp package itself as an early example.

We will try to touch base with CRAN package authors using both GSL and Rcpp to see how this can help them. The API in our package may well be incomplete, but we are always happy to try to respond to requests for additional features brought to our attention, preferably via the rcpp-devel list.

More information is on the RcppGSL page. Questions, comments etc should go to the rcpp-devel mailing list off the R-Forge page.

/code/rcpp | permanent link