[Rcpp-devel] Organization of C++ classes for glm families

[Douglas Bates]

On Tue, Mar 16, 2010 at 3:24 PM, Dirk Eddelbuettel <edd at debian.org> wrote:

> That's an interesting problem.  AFAICT folks prefer iterators and transform()
> and its ilk because it allows you to swap vector() for, say, list() with minimal
> fuss.  It is less clear that you pick up speed that way.  Wouldn't the cost
> of log() dominate the looping, whether it's old school C style or via iterators?

I was not sufficiently explicit.  The large amounts of time spent in
evaluating the inverse link and things like that are in the R
versions.   They use R constructions like

> poisson()$linkinv
function (eta)
pmax(exp(eta), .Machine$double.eps)
<environment: 0x1c435e0>

and it is (or, at least, was) the calls to pmax, pmin, etc. that
gobbled up the time.

Several years ago I wrote C code using .Call for the usual link
functions in the binomial family, just to speed up such cases.  For
example, the default link for the binomial is

> binomial()$linkinv
function (eta)
.Call("logit_linkinv", eta, PACKAGE = "stats")
<environment: 0x1c4b9e0>

Using Rcpp I could certainly extract the R function from the family
and set up the call but for the heavily used families I would prefer
to set up a way to recognize the link name and use compiled code.

In any case I now have code that seems to work

otherd <- '
'
bod <- '
  Rcpp::NumericVector xv(x);
  Rcpp::NumericVector y(xv.size());

  std::map<std::string, double(*)(double)> lst;
  lst["log"] = &log;
  std::transform(xv.begin(), xv.end(), y.begin(), *(lst["log"]));
  return y;
'
ff <- cfunction(signature(x = "numeric"), bod, Rcpp = TRUE, verbose = FALSE,
                otherdefs = otherd, includes = "#include <cmath>")