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

[Douglas Bates]

In profiling some of the code for generalized linear mixed models I
discovered that a considerable amount of the time was being taken up
in evaluation of some of the functions in the glm family.  Many of
these functions in the family are candidates for a std::transform
application

> str(poisson())
List of 12
 $ family    : chr "poisson"
 $ link      : chr "log"
 $ linkfun   :function (mu)
 $ linkinv   :function (eta)
 $ variance  :function (mu)
 $ dev.resids:function (y, mu, wt)
 $ aic       :function (y, n, mu, wt, dev)
 $ mu.eta    :function (eta)
 $ initialize:  expression({   ...
 $ validmu   :function (mu)
 $ valideta  :function (eta)
 $ simulate  :function (object, nsim)
 - attr(*, "class")= chr "family"

In most cases the time is taken in linkinv, variance and mu.eta.

I do have one implementation of glm families in C++ in the lme4a
package from the lme4 project on R-forge but that implementation is
based on a complex set of concrete classes using virtual functions and
doesn't take advantage of containers like std::vector or
Rcpp::NumericVector.

I think a more idiomatic C++ implementation would use containers and
the std::transform algorithm.  The characteristics like the linkfun,
linkinv and mu.eta functions are associated with the link name ("log"
in the example above).

Would it be feasible to associate the scalar transformation function,
in this case it would be the log function from the math library but in
other cases it could be like

double inverseLink(double x) {return 1/x;}

with the name of the link.  I would think of using a std::map for that
but I don't know what the class or type of the "value" slot would be.
Some exploration indicates that it may be a struct that is derived
from the std::unary_function<double, double> but right now my head is
beginning to hurt reading documentation.  Could someone give me an
example of how I would use, say, the exp or the log function in a
std::transform?

The second level would be to associate the scalar transformation
function with the name, say in a std::map but first I need to make
std::transform work with a scalar function.