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

[Romain Francois]

Le 16/03/10 21:01, Douglas Bates a écrit :
>
> 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.

I'll reply to the ful posting later, but for now the problem is that you 
are one argument short. transform needs 4 arguments.

I often visit cplusplus.com for a reference on STL containers and 
algorithms, for example: http://cplusplus.com/reference/algorithm/transform/

Also, Dirk recommended the books from Scott Meyers : Effective C++, More 
Effective C++ and Effective STL which are very good (although not 
introductory).


-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://tr.im/OIXN : raster images and RImageJ
|- http://tr.im/OcQe : Rcpp 0.7.7
`- http://tr.im/O1wO : highlight 0.1-5