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

[Douglas Bates]

On Tue, Mar 16, 2010 at 2:00 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> 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.

To follow up on my own posting, the things I have tried are various versions on

cpp <- '
  Rcpp::NumericVector xv(x);
  Rcpp::NumericVector y(xv.size());
  std::transform(xv.begin(), y.begin(), ll);
  return y;
'
ff <- cfunction(signature(x = "numeric"), cpp, Rcpp = TRUE, verbose = FALSE,
                otherdefs = "double ll(const double xx){return log(xx);}",
                includes = "#include <cmath>")

which always produces an error of the form

:14: error: no matching function for call to ‘transform(double*,
double*, double (&)(double))’

I'm at the point one often gets to with a new language that my code
will not compile and I have no idea of how to fix it.