Following the
Rcpp 0.10.0 release, I had written about
simulating pi easily by using the wonderful new
Rcpp Attributes feature.
Now with
Rcpp 0.10.1 released a good week ago, it is time to look
at how
Rcpp Attributes can help with external libraries. As this posts aims to show, it is a breeze!
One key aspect is the use of the
plugins for the
inline package. They provide
something akin to a callback mechanism so that compilation and linking steps can be informed about header and library locations and names.
We are going to illustrate this with an example from
GNU Scientific Library (GSL). The example I picked
uses B-Spline estimation from the GSL. This is a little redundant as R has its own spline routines and package, but serves well as a simple
illustration---and, by reproducing an existing example, followed an established path. So we will look at
Section 39.7 of the GSL manual which has a
complete example as a standalone C program, generating both the data and the fit via cubic B-splines.
We can decompose this two parts: data generation, and fitting. We will provide one function each, and then use both from R. These two
function will follow the
aforementioned example from Section 39.7
somewhat closely.
We start with the first function to generate the data.
// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
const int N = 200; // number of data points to fit
const int NCOEFFS = 12; // number of fit coefficients */
const int NBREAK = (NCOEFFS - 2); // nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
// [[Rcpp::export]]
Rcpp::List genData()
const size_t n = N;
size_t i;
double dy;
gsl_rng *r;
RcppGSL::vector<double> w(n), x(n), y(n);
gsl_rng_env_setup();
r = gsl_rng_alloc(gsl_rng_default);
//printf("#m=0,S=0\n");
/* this is the data to be fitted */
for (i = 0; i < n; ++i)
double sigma;
double xi = (15.0 / (N - 1)) * i;
double yi = cos(xi) * exp(-0.1 * xi);
sigma = 0.1 * yi;
dy = gsl_ran_gaussian(r, sigma);
yi += dy;
gsl_vector_set(x, i, xi);
gsl_vector_set(y, i, yi);
gsl_vector_set(w, i, 1.0 / (sigma * sigma));
//printf("%f %f\n", xi, yi);
Rcpp::DataFrame res = Rcpp::DataFrame::create(Rcpp::Named("x") = x,
Rcpp::Named("y") = y,
Rcpp::Named("w") = w);
x.free();
y.free();
w.free();
gsl_rng_free(r);
return(res);
We include a few header files, define (in what is common for C programs) a few constants and then define a single function
genData() which returns and
Rcpp::List as a list object to R. A primary importance here are the two attributes:
one to declare a dependence on the
RcppGSL package, and one to declare the
export of the data generator function.
That is all it takes! The plugin of the
RcppGSL will provide information about the headers and library, and
Rcpp
Attributes will do the rest.
The core of the function is fairly self-explanatory, and closely follows the
original example. Space gets
allocated, the RNG is setup and a simple functional form generates some data plus noise (see below). In the original, the data is written to
the standard output; here we return it to R as three columns in a
data.frame object familiar to R users. We then free the GSL
vectors; this manual step is needed as they are implemented as C vectors which do not have a destructor.
Next, we can turn the fitting function.
// [[Rcpp::export]]
Rcpp::List fitData(Rcpp::DataFrame ds)
const size_t ncoeffs = NCOEFFS;
const size_t nbreak = NBREAK;
const size_t n = N;
size_t i, j;
Rcpp::DataFrame D(ds); // construct the data.frame object
RcppGSL::vector<double> y = D["y"]; // access columns by name,
RcppGSL::vector<double> x = D["x"]; // assigning to GSL vectors
RcppGSL::vector<double> w = D["w"];
gsl_bspline_workspace *bw;
gsl_vector *B;
gsl_vector *c;
gsl_matrix *X, *cov;
gsl_multifit_linear_workspace *mw;
double chisq, Rsq, dof, tss;
bw = gsl_bspline_alloc(4, nbreak); // allocate a cubic bspline workspace (k = 4)
B = gsl_vector_alloc(ncoeffs);
X = gsl_matrix_alloc(n, ncoeffs);
c = gsl_vector_alloc(ncoeffs);
cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
mw = gsl_multifit_linear_alloc(n, ncoeffs);
gsl_bspline_knots_uniform(0.0, 15.0, bw); // use uniform breakpoints on [0, 15]
for (i = 0; i < n; ++i) // construct the fit matrix X
double xi = gsl_vector_get(x, i);
gsl_bspline_eval(xi, B, bw); // compute B_j(xi) for all j
for (j = 0; j < ncoeffs; ++j) // fill in row i of X
double Bj = gsl_vector_get(B, j);
gsl_matrix_set(X, i, j, Bj);
gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw); // do the fit
dof = n - ncoeffs;
tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
Rsq = 1.0 - chisq / tss;
Rcpp::NumericVector FX(151), FY(151); // output the smoothed curve
double xi, yi, yerr;
for (xi = 0.0, i=0; xi < 15.0; xi += 0.1, i++)
gsl_bspline_eval(xi, B, bw);
gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
FX[i] = xi;
FY[i] = yi;
Rcpp::List res =
Rcpp::List::create(Rcpp::Named("X")=FX,
Rcpp::Named("Y")=FY,
Rcpp::Named("chisqdof")=Rcpp::wrap(chisq/dof),
Rcpp::Named("rsq")=Rcpp::wrap(Rsq));
gsl_bspline_free(bw);
gsl_vector_free(B);
gsl_matrix_free(X);
gsl_vector_free(c);
gsl_matrix_free(cov);
gsl_multifit_linear_free(mw);
y.free();
x.free();
w.free();
return(res);
The second function closely follows the second part of the GSL example and, given the input data, fits the output data. Data structures are
setup, the spline basis is created, data is fit and then the fit is evaluated at a number of points. These two vectors are returned along
with two goodness of fit measures.
We only need to load the
Rcpp package and source a file containing the two
snippets shown above, and we are ready to deploy this:
library(Rcpp)
sourceCpp("bSpline.cpp") # compile two functions
dat <- genData() # generate the data
fit <- fitData(dat) # fit the model, returns matrix and gof measures
And with that, we generate a chart such as
via a simple four lines, or as much as it took to create the C++ functions, generate the data and fit it!
op <- par(mar=c(3,3,1,1))
plot(dat[,"x"], dat[,"y"], pch=19, col="#00000044")
lines(fit[[1]], fit[[2]], col="orange", lwd=2)
par(op)
The RcppArmadillo and RcppEigen package support plugin use in the same way. Add an attribute to export a function, and an attribute for the
depends -- and you're done. Extending R with (potentially much faster) C++ code has never been easier, and opens a whole new set of doors.
A long time ago, I spent many many hours playing Frontier (Elite 2) on
the family Atari ST. It was nice-looking (oh, those sunrises on
distant planets!), playable, varied, had ample room for progression,
very different possible roles, a huge universe of which nobody could
ever hope of exploring more than a tiny fraction, and it was basically
the first game I encountered with no set goals. You are a spaceship
pilot in a galaxy colonised left and right, and it's up to you to
decide what you want the game to be. Pirate, headhunter, mercenary,
trader, taxi, explorer, there's no end goal to achieve except what you
fix for yourself. That was, for me as a youngster, a mind-opening
experience, which I never felt as strongly since (with the possible
exception of the Creatures game, but I barely tried it).
Then came modernity, the PC, and the First Encounters sequel
(Elite 3). It added moderate complexity to the gameplay, some sort of
political evolution among the factions and a kind of plot so the
galaxy was no longer quite so static, but the end results were still a
bit disappointing. It felt like Frontier with a few
not-exactly-overwhelming textures stuck onto the ships, a loose plot
arc, and a third faction, the Alliance of Independent Systems, beyond
the Federation and the Empire. Not that exciting, although the name
of the Alliance's capital star system stayed in my mind.
So, like many others, I came to wait for Elite 4. That wait started
in 1995, more than seventeen years ago. During that time, Duke Nukem
Forever was announced, promised, developed, postponed, abandoned,
found again, developed again, and it was even actually released in
2011, while Elite 4 went into the same kind of hell only without the
actual release (yet).
At some point I started trying various substitutes to help pass time.
There was Parsec which was a nice-looking space combat simulator,
multiplayer and all, but never felt like finished. Development slowed
down to a thin trickle, and it seems like it only recently started
again as Open Parsec. Hopefully it won't
wither off again, but since it's only focused on combat, it's not
exactly the same spirit as Elite.
There's Oolite, but it looks and feels very
much like a lightly-modernised Frontier. Better textures, sure, but
the UI still looks like it got taken from the times where games
mentioned they required an EGA/VGA video card.
There's Vega Strike. It's very promising
too, and I did try it at some point, but it was crashing too often for
my tastes. Development seems to have slowed down, but not stopped;
maybe it's time I gave it another go.
I just recently found out about
Pioneer, and I admit I'm excited.
Really excited. It feels like like Frontier felt at the time. It's
not finished yet (the on-planet cities in particular look a bit
weird), but from the alpha I tried and the videos on the website it's
really impressive. I'm going to keep an eye on that, knowing I could
spend quite some time in there even before it's quite finished, if
only to fly a spaceship in the canyons of Europa with Jupiter rising
in the distance. I'll let you know if I find a black monolith.
But the main event triggering this blog post is that apparently
Elite 4 is not dead (yet). Much better, it's under heavy development,
under the Elite Dangerous name, and the current state seems to be
rather good already. And Frontier Developments, the company behind
it, opened a
Kickstarter campaign
to fund the project. The screenshots and the videos look impressive,
and the interviews of the lead developer imply that the gameplay will
be awesome. The galaxy will apparently be really dynamic, and the
player's actions could really influence its geopolitical
(galactopolitical?) evolution, the profitability of trade routes, the
prevalence and repartition of space pirates in various sectors, and so
on. Combine with a multi-player mode, and this could become the
greatest game ever made (for my personal taste at least). So I'm very
much hoping the funding campaign reaches its target.
The release date is set to March 2014. It's going to be a long wait,
fraught with incertitudes.
But then again, Duke Nukem Forever was eventually released, wasn't it?
[No, I never played the original Elite. I'm not that old.]
Yeah, wait, what?
No, really!
I took a look through 
Have you ever heard someone talking extensively about Cloud Computing or generally Clouds? and have you ever noticed the fact many people (even the ones who present themselves as experts) don t really understand what a Cloud is at all? That happened to me multiple times and one of the most common misunderstandings is many see the Cloud as something being on the internet. Many companies add a little logo representing a cloud on their frontpage and without a single change on their infrastructure (but surely with a price increment) they start calling their products as being on the Cloud. Given the lack of knowledge about this specific topic people tend to buy the product presented as being on the Cloud without understanding what they really bought.
But what Cloud Computing really means? it took several years and more than fifteen drafts to the National Institute of Standards and Technology s (NIST) to find a definition. The final accepted proposal:
Before we begin, everyone should read up on hashtables and what open
addressing / closed hashing is. The context is lines 111 190 of 
So I finally got the first laser-cut version of the 
Some brief software updates:
Here is a bit of simple (yet subtly-flawed) sql, which produces different answers on different database engines: 
I've recently replaced my laptop's HDD as the S.M.A.R.T. data
suggested a potential failure. I took the opportunity to
try out beta4 of 