[Rcpp-commits] r3706 - in pkg/RcppGSL/vignettes: . RcppGSL
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 22 00:18:37 CEST 2012
Author: edd
Date: 2012-07-22 00:18:37 +0200 (Sun, 22 Jul 2012)
New Revision: 3706
Added:
pkg/RcppGSL/vignettes/RcppGSL-intro.Rnw
pkg/RcppGSL/vignettes/RcppGSL/RcppGSL-intro.Rnw
Removed:
pkg/RcppGSL/vignettes/RcppGSL.Rnw
pkg/RcppGSL/vignettes/RcppGSL/RcppGSL.Rnw
Modified:
pkg/RcppGSL/vignettes/Makefile
Log:
o remame main vignette to RcppGSL-intro.pdf as ref.manual is also called RcppGSL.pdf
o updated Makefile accordingly
Modified: pkg/RcppGSL/vignettes/Makefile
===================================================================
--- pkg/RcppGSL/vignettes/Makefile 2012-07-21 22:01:50 UTC (rev 3705)
+++ pkg/RcppGSL/vignettes/Makefile 2012-07-21 22:18:37 UTC (rev 3706)
@@ -8,12 +8,8 @@
## on Ubuntu so for now Dirk will insist on pdflatex and this helps.
whoami=$(shell whoami)
-all: RcppGSL-unitTests.pdf RcppGSL.pdf
+all: RcppGSL-unitTests.pdf RcppGSL-intro.pdf
-oldclean:
- rm -f RcppGSL.Rnw
- cp RcppGSL/RcppGSL-fake.Rnw RcppGSL.Rnw
-
clean:
rm -f RcppGSL*.log RcppGSL*.aux RcppGSL*.out RcppGSL*.tex *.blg *.bbl *~
rm -rf auto/
@@ -35,20 +31,20 @@
$(RSCRIPT) --vanilla -e "tools::texi2dvi( 'RcppGSL-unitTests.tex', pdf = TRUE, clean = TRUE )"
rm -fr RcppGSL-unitTests.tex
-RcppGSL.pdf: RcppGSL/RcppGSL.Rnw
- cp -f RcppGSL/RcppGSL.Rnw .
- $(RSCRIPT) --vanilla -e "require(highlight); driver <- HighlightWeaveLatex(boxes = TRUE); Sweave( 'RcppGSL.Rnw', driver = driver ); "
+RcppGSL-intro.pdf: RcppGSL/RcppGSL-intro.Rnw
+ cp -f RcppGSL/RcppGSL-intro.Rnw .
+ $(RSCRIPT) --vanilla -e "require(highlight); driver <- HighlightWeaveLatex(boxes = TRUE); Sweave( 'RcppGSL-intro.Rnw', driver = driver ); "
ifneq (,$(findstring edd,$(whoami)))
- pdflatex RcppGSL
+ pdflatex RcppGSL-intro
else
- $(RSCRIPT) -e "tools::texi2dvi( 'RcppGSL.tex', pdf = TRUE, clean = FALSE )"
+ $(RSCRIPT) -e "tools::texi2dvi( 'RcppGSL-intro.tex', pdf = TRUE, clean = FALSE )"
endif
- bibtex RcppGSL
+ bibtex RcppGSL-intro
ifneq (,$(findstring edd,$(whoami)))
- pdflatex RcppGSL
- pdflatex RcppGSL
+ pdflatex RcppGSL-intro
+ pdflatex RcppGSL-intro
else
- $(RSCRIPT) -e "tools::texi2dvi( 'RcppGSL.tex', pdf = TRUE, clean = TRUE )"
+ $(RSCRIPT) -e "tools::texi2dvi( 'RcppGSL-intro.tex', pdf = TRUE, clean = TRUE )"
endif
- rm -fr RcppGSL.tex RcppGSL.bbl RcppGSL.blg RcppGSL.aux RcppGSL.out RcppGSL.log
- cp RcppGSL/RcppGSL-fake.Rnw RcppGSL.Rnw
+ rm -fr RcppGSL-intro.tex RcppGSL-intro.bbl RcppGSL-intro.blg RcppGSL-intro.aux RcppGSL-intro.out RcppGSL-intro.log
+ cp RcppGSL/RcppGSL-fake.Rnw RcppGSL-intro.Rnw
Copied: pkg/RcppGSL/vignettes/RcppGSL/RcppGSL-intro.Rnw (from rev 3705, pkg/RcppGSL/vignettes/RcppGSL/RcppGSL.Rnw)
===================================================================
--- pkg/RcppGSL/vignettes/RcppGSL/RcppGSL-intro.Rnw (rev 0)
+++ pkg/RcppGSL/vignettes/RcppGSL/RcppGSL-intro.Rnw 2012-07-21 22:18:37 UTC (rev 3706)
@@ -0,0 +1,759 @@
+\documentclass[10pt]{article}
+%\VignetteIndexEntry{RcppGSL}
+
+\usepackage{url,color}
+\usepackage[authoryear,round,longnamesfirst]{natbib}
+\usepackage[colorlinks]{hyperref}
+\definecolor{link}{rgb}{0,0,0.3} %% next few lines courtesy of RJournal.sty
+\hypersetup{
+ colorlinks,%
+ citecolor=link,%
+ filecolor=link,%
+ linkcolor=link,%
+ urlcolor=link
+}
+\usepackage{vmargin}
+\setmargrb{0.75in}{0.75in}{0.75in}{0.75in}
+\usepackage{booktabs} % fancier \hrule
+
+\newcommand{\proglang}[1]{\textsf{#1}}
+\newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}}
+
+<<echo=FALSE,print=FALSE>>=
+library( "RcppGSL" )
+options("width"=65)
+rcppgsl.version <- packageDescription( "RcppGSL" )$Version
+prettyDate <- format(Sys.Date(), "%B %e, %Y")
+require( inline )
+require( RcppGSL )
+@
+% closing $ needed here
+
+\author{Dirk Eddelbuettel \and Romain Fran\c{c}ois}
+\title{\pkg{RcppGSL}: Easier \pkg{GSL} use from \proglang{R} via \pkg{Rcpp}}
+\date{Version \Sexpr{rcppgsl.version} as of \Sexpr{prettyDate}}
+
+\begin{document}
+\maketitle
+
+\abstract{
+ \shortcites{GSL} %% so that we get Galassi et al instead of all names
+ \noindent
+ The GNU Scientific Library, or \pkg{GSL}, is a collection of numerical
+ routines for scientific computing \citep{GSL}. It is particularly useful for
+ \proglang{C} and \proglang{C++} programs as it provides a standard
+ \proglang{C} interface to a wide range of mathematical routines such as
+ special functions, permutations, combinations, fast fourier transforms,
+ eigensystems, random numbers, quadrature, random distributions,
+ quasi-random sequences, Monte Carlo integration, N-tuples, differential
+ equations, simulated annealing, numerical differentiation, interpolation,
+ series acceleration, Chebyshev approximations, root-finding, discrete
+ Hankel transforms physical constants, basis splines and wavelets. There
+ are over 1000 functions in total with an extensive test suite.
+
+ The \pkg{RcppGSL} package provides an easy-to-use interface between
+ \pkg{GSL} data structures and \proglang{R} using concepts from \pkg{Rcpp}
+ \citep{CRAN:Rcpp} which is itself a package that eases the interfaces
+ between \proglang{R} and C++.
+}
+
+\section{Introduction}
+
+The GNU Scientific Library, or \pkg{GSL}, is a collection of numerical
+routines for scientific computing \citep{GSL}. It is a rigourously developed
+and tested library providing support for a wide range of scientific or
+numerical tasks. Among the topics covered in the GSL are
+%% from the GSL manual
+complex numbers, roots of polynomials,
+special functions, vector and matrix data structures,
+permutations, combinations, sorting, BLAS support,
+linear algebra, fast fourier transforms, eigensystems,
+random numbers, quadrature, random distributions, quasi-random sequences,
+Monte Carlo integration, N-tuples,
+differential equations, simulated annealing,
+numerical differentiation, interpolation,
+series acceleration, Chebyshev approximations,
+root-finding, discrete Hankel transforms
+least-squares fitting, minimization,
+physical constants, basis splines and wavelets.
+
+Support for \proglang{C} programming with the GSL is readily available: the GSL itself is written in \proglang{C}
+and provides a \proglang{C}-language Application Programming Interface
+(API). Access from \proglang{C++} is therefore possible, albeit not at the
+abstraction level that can be offered by dedicated \proglang{C++}
+implementations.\footnote{Several \proglang{C++} wrappers for the GSL have
+ been written over the years yet none reached a state of completion
+ comparable to the GSL itself. Three such wrapping library are
+ \url{http://cholm.home.cern.ch/cholm/misc/gslmm/},
+ \url{http://gslwrap.sourceforge.net/} and
+ \url{http://code.google.com/p/gslcpp/}.}
+
+The GSL is somewhat unique among numerical libraries. Its combination of
+broad coverage of scientific topics, serious implementation effort and the
+use of a FLOSS license have lead to a fairly wide usage of the library. As a
+concrete example, we can consider the the CRAN repository network for the
+\proglang{R} language and environment \citep{R:Main}. CRAN contains over a
+dozen packages interfacing the GSL: \pkg{copula}, \pkg{dynamo}, \pkg{gsl},
+\pkg{gstat}, \pkg{magnets}, \pkg{mvabund}, \pkg{QRMlib}, \pkg{RBrownie},
+\pkg{RDieHarder}, \pkg{RHmm}, \pkg{segclust}, \pkg{surveillance}, and
+\pkg{topicmodels}. This is a clear indication that the GSL is popular among
+programmers using either the \proglang{C} or \proglang{C++} language for
+solving problems applied science.
+
+At the same time, the \pkg{Rcpp} package \citep{CRAN:Rcpp} offers a
+higher-level abstraction between \proglang{R} and underlying \proglang{C++}
+(or \proglang{C}) code. \pkg{Rcpp} permits \proglang{R} objects like
+vectors, matrices, lists, functions, environments, $\ldots$, to be
+manipulated directly at the \proglang{C++} level, alleviates the needs for
+complicated and error-prone parameter passing and memory allocation. It also
+permits compact vectorised expressions similar to what can be written in
+\proglang{R} directly at the \proglang{C++} level.
+
+The \pkg{RcppGSL} package discussed here aims the help close the gap. It
+tries to offer access to GSL functions, in particular via the vector and
+matrix data structures used throughout the GSL, while staying closer to the
+`whole object model' familar to the \proglang{R} programmer.
+
+The rest of paper is organised as follows. The next section shows a
+motivating example of a fast linear model fit routine using \pkg{GSL} functions.
+The following section discusses support for \pkg{GSL} vector types, which is
+followed by a section on matrices.
+
+\section{Motivation: FastLm}
+
+Fitting linear models is a key building block of analysing data and
+modeling. \proglang{R} has a very complete and feature-rich function in
+\texttt{lm()} which can provide a model fit as we a number of diagnostic
+measure, either directly or via the corresponding \texttt{summary()} method
+for linear model fits. The \texttt{lm.fit()} function also provides a faster
+alternative (which is however recommend only for for advanced users) which
+provides estimates only and fewer statistics for inference. This sometimes
+leads users request a routine which is both fast and featureful enough. The
+\texttt{fastLm} routine shown here provides such an implementation. It uses
+the \pkg{GSL} for the least-squares fitting functions and therefore provides a nice
+example for \pkg{GSL} integration with \proglang{R}.
+
+<<lang=cpp,size=small>>=
+#include <RcppGSL.h>
+#include <gsl/gsl_multifit.h>
+#include <cmath>
+
+extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
+
+ try {
+ RcppGSL::vector<double> y = ys; // create gsl data structures from SEXP
+ RcppGSL::matrix<double> X = Xs;
+
+ int n = X.nrow(), k = X.ncol();
+ double chisq;
+
+ RcppGSL::vector<double> coef(k); // to hold the coefficient vector
+ RcppGSL::matrix<double> cov(k,k); // and the covariance matrix
+
+ // the actual fit requires working memory we allocate and free
+ gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (n, k);
+ gsl_multifit_linear (X, y, coef, cov, &chisq, work);
+ gsl_multifit_linear_free (work);
+
+ // extract the diagonal as a vector view
+ gsl_vector_view diag = gsl_matrix_diagonal(cov) ;
+
+ // currently there is not a more direct interface in Rcpp::NumericVector
+ // that takes advantage of wrap, so we have to do it in two steps
+ Rcpp::NumericVector std_err ; std_err = diag;
+ std::transform( std_err.begin(), std_err.end(), std_err.begin(), sqrt);
+
+ Rcpp::List res = Rcpp::List::create(Rcpp::Named("coefficients") = coef,
+ Rcpp::Named("stderr") = std_err,
+ Rcpp::Named("df") = n - k);
+
+ // free all the GSL vectors and matrices -- as these are really C data structures
+ // we cannot take advantage of automatic memory management
+ coef.free(); cov.free(); y.free(); X.free();
+
+ return res; // return the result list to R
+
+ } catch( std::exception &ex ) {
+ forward_exception_to_r( ex );
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)" );
+ }
+ return R_NilValue; // -Wall
+}
+@
+
+We first initialize a \textsl{RcppGSL} vector and matrix, each templated to
+the standard numeric type \texttt{double} (and the GSL supports other types
+ranging from lower precision floating point to signed and unsigned integers
+as well as complex numbers). We the reserve another vector and matrix to
+hold the resulting coefficient estimates as well as the estimate of the
+covariance matrix. Next, we allocate workspace using a GSL routine, fit the
+linear model and free the workspace. The next step involves extracting the
+diagonal element from the covariance matrix. We then employ a so-called
+iterator---a common \proglang{C++} idiom from the Standard Template Library
+(STL)---to iterate over the vector of diagonal and transforming it by
+applying the square root function to compute our standard error of the
+estimate. Finally we create a named list with the return value before we free
+temporary memory allocation (a step that has to be done because the
+underlying objects are really \proglang{C} objects conforming to the
+\pkg{GSL} interface and hence without the automatic memory management we
+could have with \proglang{C++} vector or matrix structures as used through
+the \pkg{Rcpp} package) and return the result to \proglang{R}.
+
+We should note that \pkg{RcppArmadillo} \citep{CRAN:RcppArmadillo} implements
+a matching \texttt{fastLm} function using the Armadillo library by
+\cite{Sanderson:2010:Armadillo}, and can do so with more compact code due to
+\proglang{C++} features.
+
+\section{Vectors}
+
+This section details the different vector represenations, starting with their
+definition inside the \pkg{GSL}. We then discuss our layering before showing
+how the two types map. A discussion of read-only `vector view' classes
+concludes the section.
+
+\subsection{\pkg{GSL} Vectors}
+
+\pkg{GSL} defines various vector types to manipulate one-dimensionnal
+data, similar to \proglang{R} arrays. For example the \verb|gsl_vector| and \verb|gsl_vector_int|
+structs are defined as:
+
+<<lang=cpp,size=small>>=
+typedef struct{
+ size_t size;
+ size_t stride;
+ double * data;
+ gsl_block * block;
+ int owner;
+} gsl_vector;
+
+typedef struct {
+ size_t size;
+ size_t stride;
+ int * data;
+ gsl_block_int * block;
+ int owner;
+}
+gsl_vector_int;
+@
+
+A typical use of the \verb|gsl_vector| struct is given below:
+
+<<lang=cpp,size=small>>=
+int i;
+gsl_vector * v = gsl_vector_alloc (3); // allocate a gsl_vector of size 3
+
+for (i = 0; i < 3; i++) { // fill the vector
+ gsl_vector_set (v, i, 1.23 + i);
+}
+
+double sum = 0.0 ; // access elements
+for (i = 0; i < 3; i++) {
+ sum += gsl_vector_set( v, i ) ;
+}
+
+gsl_vector_free (v); // free the memory
+@
+
+\subsection{RcppGSL::vector}
+
+\pkg{RcppGSL} defines the template \texttt{RcppGSL::vector<T>} to manipulate
+\verb|gsl_vector| pointers taking advantage of C++ templates. Using this
+template type, the previous example now becomes:
+
+<<lang=cpp,size=small>>=
+int i;
+RcppGSL::vector<double> v(3); // allocate a gsl_vector of size 3
+
+for (i = 0; i < 3; i++) { // fill the vector
+ v[i] = 1.23 + i ;
+}
+
+double sum = 0.0 ; // access elements
+for (i = 0; i < 3; i++) {
+ sum += v[i] ;
+}
+
+v.free() ; // free the memory
+@
+
+The class \texttt{RcppGSL::vector<double>} is a smart pointer, that can be used
+anywhere where a raw pointer \verb|gsl_vector| can be used, such as the
+\verb|gsl_vector_set| and \verb|gsl_vector_get| functions above.
+
+Beyond the convenience of a nicer syntax for allocation and release of
+memory, the \texttt{RcppGSL::vector} template faciliates interchange of
+\pkg{GSL} vectors with \pkg{Rcpp} objects, and hence \pkg{R} objects. The
+following example defines a \texttt{.Call} compatible function called
+\verb|sum_gsl_vector_int| that operates on a \verb|gsl_vector_int| through
+the \texttt{RcppGSL::vector<int>} template specialization:
+
+<<lang=cpp,size=small>>=
+RCPP_FUNCTION_1( int, sum_gsl_vector_int, RcppGSL::vector<int> vec){
+ int res = std::accumulate( vec.begin(), vec.end(), 0 ) ;
+ vec.free() ; // we need to free vec after use
+ return res ;
+}
+@
+
+The function can then simply be called from \proglang{R} :
+
+<<results=hide,echo=FALSE>>=
+fx <- cxxfunction(
+ list( sum_gsl_vector_int = signature( vec_ = "integer" ) ),
+ c( sum_gsl_vector_int = '
+ RcppGSL::vector<int> vec = as< RcppGSL::vector<int> >(vec_) ;
+ int res = std::accumulate( vec.begin(), vec.end(), 0 ) ;
+ vec.free() ; // we need to free vec after use
+ return wrap( res ) ;
+' ), plugin = "RcppGSL" )
+<<size=small>>=
+.Call( "sum_gsl_vector_int", 1:10 )
+@
+
+A second example shows a simple function that grabs elements of an
+R list as \verb|gsl_vector| objects using implicit conversion mechanisms
+of \pkg{Rcpp}
+
+<<lang=cpp,size=small>>=
+RCPP_FUNCTION_1( double, gsl_vector_sum_2, Rcpp::List data ){
+ // grab "x" as a gsl_vector through the RcppGSL::vector<double> class
+ RcppGSL::vector<double> x = data["x"] ;
+
+ // grab "y" as a gsl_vector through the RcppGSL::vector<int> class
+ RcppGSL::vector<int> y = data["y"] ;
+ double res = 0.0 ;
+ for( size_t i=0; i< x->size; i++){
+ res += x[i] * y[i] ;
+ }
+
+ // as usual with GSL, we need to explicitely free the memory
+ x.free() ;
+ y.free() ;
+
+ // return the result
+ return res ;
+}
+@
+
+called from \proglang{R} :
+
+<<results=hide, echo=FALSE>>=
+fx2 <- cxxfunction( list( gsl_vector_sum_2 = signature( data_ = "list" ) ),
+c( gsl_vector_sum_2 = '
+ List data(data_) ;
+
+ // grab "x" as a gsl_vector through the RcppGSL::vector<double> class
+ RcppGSL::vector<double> x = data["x"] ;
+
+ // grab "y" as a gsl_vector through the RcppGSL::vector<int> class
+ RcppGSL::vector<int> y = data["y"] ;
+ double res = 0.0 ;
+ for( size_t i=0; i< x->size; i++){
+ res += x[i] * y[i] ;
+ }
+
+ // as usual with GSL, we need to explicitely free the memory
+ x.free() ;
+ y.free() ;
+
+ // return the result
+ return wrap(res);
+' ), plugin = "RcppGSL" )
+@
+<<>>=
+data <- list( x = seq(0,1,length=10), y = 1:10 )
+.Call( "gsl_vector_sum_2", data )
+@
+
+\subsection{Mapping}
+
+Table~\ref{tab:mappingVectors} shows the mapping between types defined by the
+\pkg{GSL} and their corresponding types in the \pkg{RcppGSL} package.
+
+\begin{table}[htb]
+ \centering
+ \begin{small}
+ \begin{tabular}{ll}
+ \toprule
+ gsl vector & RcppGSL \\
+ \midrule
+ \texttt{gsl\_vector} & \texttt{RcppGSL::vector<double>} \\
+ \texttt{gsl\_vector\_int} & \texttt{RcppGSL::vector<int>} \\
+ \texttt{gsl\_vector\_float} & \texttt{RcppGSL::vector<float>} \\
+ \texttt{gsl\_vector\_long} & \texttt{RcppGSL::vector<long>} \\
+ \texttt{gsl\_vector\_char} & \texttt{RcppGSL::vector<char>} \\
+ \texttt{gsl\_vector\_complex} & \texttt{RcppGSL::vector<gsl\_complex>} \\
+ \texttt{gsl\_vector\_complex\_float} & \texttt{RcppGSL::vector<gsl\_complex\_float>} \\
+ \texttt{gsl\_vector\_complex\_long\_double} & \texttt{RcppGSL::vector<gsl\_complex\_long\_double>} \\
+ \texttt{gsl\_vector\_long\_double} & \texttt{RcppGSL::vector<long double>} \\
+ \texttt{gsl\_vector\_short} & \texttt{RcppGSL::vector<short>} \\
+ \texttt{gsl\_vector\_uchar} & \texttt{RcppGSL::vector<unsigned char>} \\
+ \texttt{gsl\_vector\_uint} & \texttt{RcppGSL::vector<unsigned int>} \\
+ \texttt{gsl\_vector\_ushort} & \texttt{RcppGSL::vector<insigned short>} \\
+ \texttt{gsl\_vector\_ulong} & \texttt{RcppGSL::vector<unsigned long>} \\
+ \bottomrule
+ \end{tabular}
+ \end{small}
+ \caption{Correspondance between \pkg{GSL} vector types and templates defined in \pkg{RcppGSL}.}
+ \label{tab:mappingVectors}
+\end{table}
+
+\subsection{ Vector Views}
+
+Several \pkg{GSL} algorithms return \pkg{GSL} vector views as their result
+type. \pkg{RcppGSL} defines the template class \texttt{RcppGSL::vector\_view}
+to handle vector views using \proglang{C++} syntax.
+
+<<lang=cpp,size=small>>=
+extern "C" SEXP test_gsl_vector_view(){
+ int n = 10 ;
+ RcppGSL::vector<double> v(n) ;
+ for( int i=0 ; i<n; i++){
+ v[i] = i ;
+ }
+ RcppGSL::vector_view<double> v_even = gsl_vector_subvector_with_stride(v,0,2,n/2);
+ RcppGSL::vector_view<double> v_odd = gsl_vector_subvector_with_stride(v,1,2,n/2);
+
+ List res = List::create(
+ _["even"] = v_even,
+ _["odd" ] = v_odd
+ ) ;
+ v.free() ; // we only need to free v, the views do not own data
+ return res ;
+}
+@
+
+As with vectors, \proglang{C++} objects of type
+\texttt{RcppGSL::vector\_view} can be converted implicitly to their
+associated \pkg{GSL} view type. Table~\ref{tab:mappingVectorViews} displays
+the pairwise correspondance so that the \proglang{C++} objects can be passed
+to compatible \pkg{GSL} algorithms.
+
+\begin{table}[htb]
+ \centering
+ \begin{small}
+ \begin{tabular}{ll}
+ \toprule
+ gsl vector views & RcppGSL \\
+ \midrule
+ \texttt{gsl\_vector\_view} & \texttt{RcppGSL::vector\_view<double>} \\
+ \texttt{gsl\_vector\_view\_int} & \texttt{RcppGSL::vector\_view<int>} \\
+ \texttt{gsl\_vector\_view\_float} & \texttt{RcppGSL::vector\_view<float>} \\
+ \texttt{gsl\_vector\_view\_long} & \texttt{RcppGSL::vector\_view<long>} \\
+ \texttt{gsl\_vector\_view\_char} & \texttt{RcppGSL::vector\_view<char>} \\
+ \texttt{gsl\_vector\_view\_complex} & \texttt{RcppGSL::vector\_view<gsl\_complex>} \\
+ \texttt{gsl\_vector\_view\_complex\_float} & \texttt{RcppGSL::vector\_view<gsl\_complex\_float>} \\
+ \texttt{gsl\_vector\_view\_complex\_long\_double} & \texttt{RcppGSL::vector\_view<gsl\_complex\_long\_double>} \\
+ \texttt{gsl\_vector\_view\_long\_double} & \texttt{RcppGSL::vector\_view<long double>} \\
+ \texttt{gsl\_vector\_view\_short} & \texttt{RcppGSL::vector\_view<short>} \\
+ \texttt{gsl\_vector\_view\_uchar} & \texttt{RcppGSL::vector\_view<unsigned char>} \\
+ \texttt{gsl\_vector\_view\_uint} & \texttt{RcppGSL::vector\_view<unsigned int>} \\
+ \texttt{gsl\_vector\_view\_ushort} & \texttt{RcppGSL::vector\_view<insigned short>} \\
+ \texttt{gsl\_vector\_view\_ulong} & \texttt{RcppGSL::vector\_view<unsigned long>} \\
+ \bottomrule
+ \end{tabular}
+ \end{small}
+ \caption{Correspondance between \pkg{GSL} vector view types and templates defined
+ in \pkg{RcppGSL}.}
+ \label{tab:mappingVectorViews}
+\end{table}
+
+The vector view class also contains a conversion operator to automatically
+transform the data of the view object to a \pkg{GSL} vector object. This
+enables use of vector views where \pkg{GSL} would expect a vector.
+
+
+\section{Matrices}
+
+The \pkg{GSL} also defines a set of matrix data types : \texttt{gsl\_matrix},
+\texttt{gsl\_matrix\_int} etc ... for which \pkg{RcppGSL} defines a corresponding
+convenience \proglang{C++} wrapper generated by the \texttt{RcppGSL::matrix}
+template.
+
+\subsection{Creating matrices}
+
+The \texttt{RcppGSL::matrix} template exposes three constructors.
+
+<<lang=cpp,size=small>>=
+// convert an R matrix to a GSL matrix
+matrix( SEXP x) throw(::Rcpp::not_compatible)
+
+// encapsulate a GSL matrix pointer
+matrix( gsl_matrix* x)
+
+// create a new matrix with the given number of rows and columns
+matrix( int nrow, int ncol)
+@
+
+\subsection{Implicit conversion}
+
+\texttt{RcppGSL::matrix} defines implicit conversion to a pointer to
+the associated \pkg{GSL} matrix type, as well as dereferencing operators, making
+the class \texttt{RcppGSL::matrix} look and feel like a pointer to a GSL
+matrix type.
+
+<<lang=cpp,size=small>>=
+gsltype* data ;
+operator gsltype*(){ return data ; }
+gsltype* operator->() const { return data; }
+gsltype& operator*() const { return *data; }
+@
+
+\subsection{Indexing}
+
+Indexing of \pkg{GSL} matrices is usually the task of the functions
+\texttt{gsl\_matrix\_get}, \texttt{gsl\_matrix\_int\_get}, ... and
+\texttt{gsl\_matrix\_set}, \texttt{gsl\_matrix\_int\_set}, ...
+
+\pkg{RcppGSL} takes advantage of both operator overloading and templates
+to make indexing a \pkg{GSL} matrix much more convenient.
+
+<<lang=cpp,size=small>>=
+RcppGSL::matrix<int> mat(10,10); // create a matrix of size 10x10
+
+for( int i=0; i<10: i++) { // fill the diagonal
+ mat(i,i) = i ;
+}
+@
+
+\subsection{Methods}
+
+The \texttt{RcppGSL::matrix} type also defines the following member functions:
+\begin{quote}
+ \begin{itemize}
+ \item[\texttt{nrow}] extracts the number of rows
+ \item[\texttt{ncol}] extract the number of columns
+ \item[\texttt{size}] extracts the number of elements
+ \item[\texttt{free}] releases the memory
+ \end{itemize}
+\end{quote}
+
+\subsection{Matrix views}
+
+Similar to the vector views discussed above, the \pkg{RcppGSL} also provides
+an implicit conversion operator which returns the underlying matrix stored in
+the matrix view class.
+
+\section{Using \pkg{RcppGSL} in your package}
+
+The \pkg{RcppGSL} package contains a complete example providing a single
+function \texttt{colNorm} which computes a norm for each column of a
+supplied matrix. This example adapts a matrix example from the GSL manual that has
+been chose merely as a means to showing how to set up a package to use
+\pkg{RcppGSL}.
+
+Needless to say, we could compute such a matrix norm easily in \proglang{R}
+using existing facilities. One such possibility is a simple
+\verb|apply(M, 2, function(x) sqrt(sum(x^2)))| as shown on the corresponding
+help page in the example package inside \pkg{RcppGSL}. One point in favour of
+using the \pkg{GSL} code is that it employs a BLAS function so on
+sufficiently large matrices, and with suitable BLAS libraries installed, this
+variant could be faster due to the optimised code in high-performance BLAS
+libraries and/or the inherent parallelism a multi-core BLAS variant which
+compute compute the vector norm in parallel. On all `reasonable' matrix
+sizes, however, the performance difference should be neglible.
+
+\subsection{The \texttt{configure} script}
+
+\subsubsection{Using autoconf}
+
+Using \pkg{RcppGSL} means employing both the \pkg{GSL} and \proglang{R}. We
+may need to find the location of the \pkg{GSL} headers and library, and this
+done easily from a \texttt{configure} source script which \texttt{autoconf}
+generates from a \texttt{configure.in} source file such as the following:
+
+<<lang=sh,size=small>>=
+AC_INIT([RcppGSLExample], 0.1.0)
+
+## Use gsl-config to find arguments for compiler and linker flags
+##
+## Check for non-standard programs: gsl-config(1)
+AC_PATH_PROG([GSL_CONFIG], [gsl-config])
+## If gsl-config was found, let's use it
+if test "${GSL_CONFIG}" != ""; then
+ # Use gsl-config for header and linker arguments (without BLAS which we get from R)
+ GSL_CFLAGS=`${GSL_CONFIG} --cflags`
+ GSL_LIBS=`${GSL_CONFIG} --libs-without-cblas`
+else
+ AC_MSG_ERROR([gsl-config not found, is GSL installed?])
+fi
+
+## Use Rscript to query Rcpp for compiler and linker flags
+## link flag providing libary as well as path to library, and optionally rpath
+RCPP_LDFLAGS=`${R_HOME}/bin/Rscript -e 'Rcpp:::LdFlags()'`
+
+# Now substitute these variables in src/Makevars.in to create src/Makevars
+AC_SUBST(GSL_CFLAGS)
+AC_SUBST(GSL_LIBS)
+AC_SUBST(RCPP_LDFLAGS)
+
+AC_OUTPUT(src/Makevars)
+@
+
+Such a source \texttt{configure.in} gets converted into a script
+\texttt{configure} by invoking the \texttt{autoconf} program.
+
+\subsubsection{Using functions provided by RcppGSL}
+
+\pkg{RcppGSL} provides R functions that allows one to retrieve the same
+information. Therefore the configure script can also be written as:
+
+<<lang=bash,size=small>>=
+#!/bin/sh
+
+GSL_CFLAGS=`${R_HOME}/bin/Rscript -e "RcppGSL:::CFlags()"`
+GSL_LIBS=`${R_HOME}/bin/Rscript -e "RcppGSL:::LdFlags()"`
+RCPP_LDFLAGS=`${R_HOME}/bin/Rscript -e "Rcpp:::LdFlags()"`
+
+sed -e "s|@GSL_LIBS@|${GSL_LIBS}|" \
+ -e "s|@GSL_CFLAGS@|${GSL_CFLAGS}|" \
+ -e "s|@RCPP_LDFLAGS@|${RCPP_LDFLAGS}|" \
+ src/Makevars.in > src/Makevars
+@
+
+Similarly, the configure.win for windows can be written as:
+
+<<lang=bash,size=small>>=
+GSL_CFLAGS=`${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe -e "RcppGSL:::CFlags()"`
+GSL_LIBS=`${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe -e "RcppGSL:::LdFlags()"`
+RCPP_LDFLAGS=`${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe -e "Rcpp:::LdFlags()"`
+
+sed -e "s|@GSL_LIBS@|${GSL_LIBS}|" \
+ -e "s|@GSL_CFLAGS@|${GSL_CFLAGS}|" \
+ -e "s|@RCPP_LDFLAGS@|${RCPP_LDFLAGS}|" \
+ src/Makevars.in > src/Makevars.win
+@
+
+\subsection{The \texttt{src} directory}
+
+The \proglang{C++} source file takes the matrix supplied from \proglang{R}
+and applies the \pkg{GSL} function to each column.
+
+<<lang=cpp,size=small>>=
+#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
+}
+@
+
+The \proglang{Makevars.in} file governs the compilation and uses the values
+supplied by \texttt{configure} during build-time:
+
+<<lang=sh,size=small>>=
+
+# set by configure
+GSL_CFLAGS = @GSL_CFLAGS@
+GSL_LIBS = @GSL_LIBS@
+RCPP_LDFLAGS = @RCPP_LDFLAGS@
+
+# combine with standard arguments for R
+PKG_CPPFLAGS = $(GSL_CFLAGS)
+PKG_LIBS = $(GSL_LIBS) $(RCPP_LDFLAGS)
+@
+
+The variables surrounded by \@ will be filled by \texttt{configure} during
+package build-time.
+
+\subsection{The \texttt{R} directory}
+
+The \proglang{R} source is very simply: a single matrix is passed to \proglang{C++}:
+
+<<size=small>>=
+colNorm <- function(M) {
+ stopifnot(is.matrix(M))
+ res <- .Call("colNorm", M, package="RcppGSLExample")
+}
+@
+
+\section{Using \pkg{RcppGSL} with \pkg{inline}}
+
+The \pkg{inline} package \citep{CRAN:inline} is very helpful for prototyping code in
+\proglang{C}, \proglang{C++} or \proglang{Fortran} as it takes care of code
+compilation, linking and dynamic loading directly from \proglang{R}. It is
+being used extensively by \pkg{Rcpp}, for example in the numerous unit tests.
+
+The example below shows how \pkg{inline} can be deployed with
+\pkg{RcppGSL}. We implement the same column norm example, but this time as an
+\proglang{R} script which is compiled, linked and loaded on-the-fly. Compared
+to standard use of \pkg{inline}, we have to make sure to add a short section
+declaring which header files from \pkg{GSL} we need to use; the \pkg{RcppGSL}
+then communicates with \pkg{inline} to tell it about the location and names
+of libraries used to build code against \pkg{GSL}.
+
+<<size=small>>=
+require(inline)
+
+inctxt='
+ #include <gsl/gsl_matrix.h>
+ #include <gsl/gsl_blas.h>
+'
+
+bodytxt='
+ 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
+'
+
+foo <- cxxfunction(
+ signature(sM="numeric"),
+ body=bodytxt, inc=inctxt, plugin="RcppGSL")
+
+## see Section 8.4.13 of the GSL manual: create M as a sum of two outer products
+M <- outer(sin(0:9), rep(1,10), "*") + outer(rep(1, 10), cos(0:9), "*")
+foo(M)
+@
+
+The \texttt{RcppGSL} inline plugin supports creation of a package skeleton
+based on the inline function.
+
+<<eval=FALSE>>=
+package.skeleton( "mypackage", foo )
+@
+
+\section{Summary}
+
+The GNU Scientific Library (GSL) by \citet{GSL} offers a very comprehensive
+collection of rigorously developed and tested functions for applied
+scientific computing under a common Open Source license. This has lead to
+widespread deployment of \pkg{GSL} among a number of disciplines.
+
+Using the automatic wrapping and converters offered by the \pkg{RcppGSL}
+package presented here, \proglang{R} users and programmers can now deploy
+algorithmns provided by the \pkg{GSL} with greater ease.
+
+\bibliographystyle{plainnat}
+\bibliography{RcppGSL}
+
+\end{document}
+
Deleted: pkg/RcppGSL/vignettes/RcppGSL/RcppGSL.Rnw
===================================================================
--- pkg/RcppGSL/vignettes/RcppGSL/RcppGSL.Rnw 2012-07-21 22:01:50 UTC (rev 3705)
+++ pkg/RcppGSL/vignettes/RcppGSL/RcppGSL.Rnw 2012-07-21 22:18:37 UTC (rev 3706)
@@ -1,759 +0,0 @@
-\documentclass[10pt]{article}
-%\VignetteIndexEntry{RcppGSL}
-
-\usepackage{url,color}
-\usepackage[authoryear,round,longnamesfirst]{natbib}
-\usepackage[colorlinks]{hyperref}
-\definecolor{link}{rgb}{0,0,0.3} %% next few lines courtesy of RJournal.sty
-\hypersetup{
- colorlinks,%
- citecolor=link,%
- filecolor=link,%
- linkcolor=link,%
- urlcolor=link
-}
-\usepackage{vmargin}
-\setmargrb{0.75in}{0.75in}{0.75in}{0.75in}
-\usepackage{booktabs} % fancier \hrule
-
-\newcommand{\proglang}[1]{\textsf{#1}}
-\newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}}
-
-<<echo=FALSE,print=FALSE>>=
-library( "RcppGSL" )
-options("width"=65)
-rcppgsl.version <- packageDescription( "RcppGSL" )$Version
-prettyDate <- format(Sys.Date(), "%B %e, %Y")
-require( inline )
-require( RcppGSL )
-@
-% closing $ needed here
-
-\author{Dirk Eddelbuettel \and Romain Fran\c{c}ois}
-\title{\pkg{RcppGSL}: Easier \pkg{GSL} use from \proglang{R} via \pkg{Rcpp}}
-\date{Version \Sexpr{rcppgsl.version} as of \Sexpr{prettyDate}}
-
-\begin{document}
-\maketitle
-
-\abstract{
- \shortcites{GSL} %% so that we get Galassi et al instead of all names
- \noindent
- The GNU Scientific Library, or \pkg{GSL}, is a collection of numerical
- routines for scientific computing \citep{GSL}. It is particularly useful for
- \proglang{C} and \proglang{C++} programs as it provides a standard
- \proglang{C} interface to a wide range of mathematical routines such as
- special functions, permutations, combinations, fast fourier transforms,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rcpp -r 3706
More information about the Rcpp-commits
mailing list