[Rcpp-devel] Organization of C++ classes for glm families
Romain Francois
romain.francois at dbmail.com
Tue Mar 16 21:33:08 CET 2010
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
More information about the Rcpp-devel
mailing list