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

[Douglas Bates]

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

OK, it was a dumb mistake.  I misread the documentation.  If I change
the code to

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

it's happy.  Now on to that question about std::map<std::string,
std::unary_function<double, double>>