On Fri, Mar 12, 2010 at 2:23 PM, Romain Francois <[email protected]> wrote: > Le 12/03/10 18:23, Douglas Bates a écrit : >> >> On Fri, Mar 12, 2010 at 10:47 AM, Dirk Eddelbuettel<[email protected]> >> wrote: >>> >>> Hi Doug, >>> >>> Welcome! We're honoured to have you here. And glad you picked the list. >>> >>> On 12 March 2010 at 10:14, Douglas Bates wrote: >>> | I have looked at your (i.e. Dirk and Roman's) submitted Rjournal >>> | article on Rcpp and decided that I should take the plunge. >>> | >>> | To get my feet wet I am starting with a relatively simple example >>> | interfacing to (God help us) 1960's style Fortran written by Mike >>> | Powell and incorporated in the minqa package. >>> | >>> | The value to be returned is a named list of heterogeneous values. I >>> | can see how to construct such a list in the classic Rcpp API using a >>> | RcppResultSet object but I don't know the best way to accomplish this >>> | in the more modern Rcpp namespace. Can someone point me in the right >>> | direction? Also, is there a convenient idiom for attaching an S3 >>> | class name to a list to be returned? >>> >>> Yes. See eg the examples I put on my blog, eg inside of RcppArmadillo we >>> simply do (including actual fastLm source here modulu two ifdefs dealing >>> with >>> older versions -- which means you need Armadillo 0.9.0 or later) >>> >>> #include<RcppArmadillo.h> >>> >>> extern "C" SEXP fastLm(SEXP ys, SEXP Xs) { >>> >>> Rcpp::NumericVector yr(ys); // creates Rcpp vector >>> from SEXP >>> Rcpp::NumericMatrix Xr(Xs); // creates Rcpp matrix >>> from SEXP >>> int n = Xr.nrow(), k = Xr.ncol(); >>> >>> arma::mat X(Xr.begin(), n, k, false); // reuses memory and >>> avoids extra copy >>> arma::colvec y(yr.begin(), yr.size(), false); >>> >>> arma::colvec coef = solve(X, y); // fit model y ~ X >>> >>> arma::colvec resid = y - X*coef; // residuals >>> >>> double sig2 = arma::as_scalar( trans(resid)*resid/(n-k) ); >>> // std.error of estimate >>> arma::colvec stderrest = sqrt( sig2 * diagvec( >>> arma::inv(arma::trans(X)*X)) ); >>> >>> Rcpp::Pairlist res(Rcpp::Named( "coefficients", coef), >>> Rcpp::Named( "stderr", stderrest)); >>> return res; >>> } >>> >>> [ By the way, I still think I should try to implement all four of your >>> "how >>> to do lm" solutions from the older R News / R Journal piece in Armadillo >>> one >>> day :-) And see what the timing impacts would be. ] >>> >>> That works "on the fly" with Pairlist and named. We also have Rcpp::List >>> examples, but AFAICR you need to 'predeclare' how long your list is and >>> the >>> set each element. The 200+ unit tests will have examples. Better docs >>> will >>> appear one day... >> >> Nice approach. I presume that the Pairlist object is converted to an >> VECSXP in the process of being packaged for return. > > No, they are shipped to R as a LISTSXP.
Hmm. That may cause problems in the long run. IIRC, having LISTSXPs, other than those related to language objects, loose in R is frowned upon. I'll check before saying anything definite. >> I will play >> around a bit with expressions like >> >> Rcpp::List res(Rcpp::PairList(Rcpp::Named(...), ...)) > > This should work since the List( SEXP ) constructor uses coercion. Let us > know otherwise. I tried it and was terribly frustrated that my example didn't work and I couldn't decide why. After much exploration I found that it had nothing to do with that construction and was in fact due to my using 6 SEXPs in a Pairlist constructor. Apparently my version is not compiled with HAS_VARIADIC_TEMPLATES set and that means the constructors only go up to 5. I see now that because I am using g++ 4.4.1 I should probably set the compiler flag -std=c++0x I have tried it again initializing a list from the Pairlist object and that works. Do you happen to know if that conversion determines the length of the result before allocating it? I tried to read the code for internal::ListInitialization but I am not sufficiently fluent in C++ to decide what is going on. >> I would like to have the capability of attaching an attribute named >> "class" to the VECSXP before returning it. > > res.attr( "class" ) = "foo" ; > >>> | An unrelated question -- has anyone looked at accessing S4 classed >>> | objects in Rcpp? I have had considerable experience with S4 classes >>> | and objects and could keep such a project in mind as I begin to >>> | explore Rcpp. >>> >>> I'd be happy to go there (time permitting) but very gently as I am still >>> much >>> more firmly in S3 land. Romain may be more eager :) >> >> One of the nice properties of S4 classed objects is that you can >> define validity checks for classes and bypass some of the code that >> does validity checking of arguments. > > For now the only things we have in the api to deal with S4'ness is these > methods : "isS4", "slot", "hasSlot". > > Not sure right now how we could formalize S4 objects in the api, but this > would definitely be a fine addition. > >>> | Yet another question, if i may. Does coercion take place in >>> | constructing an object like >>> | >>> | Rcpp::NumericVector aa(a); >>> | >>> | where a is an SEXP, say an SEXP argument to a C++ function, where I >>> | accidentally passed an integer vector instead of a double vector? If >>> >>> Very much so. Again, see my blog and the RcppVectorExample code. Int >>> in, >>> Double out. All magic hidden. >>> >>> And of course you can just try it: >>> >>> R> library(RcppExamples) >>> Loading required package: Rcpp >>> R> RcppVectorExample(c(1,4,9), api="new") >>> $result >>> [1] 1 2 3 >>> >>> $original >>> [1] 1 4 9 >>> >>> R> RcppVectorExample(c(2,5,10), api="new") >>> $result >>> [1] 1.414 2.236 3.162 >>> >>> $original >>> [1] 2 5 10 >>> >>> R> >> >> Perhaps you are assuming that c(1,4,9) generates an integer vector in >> this example. In fact, c(1,4,9) is a numeric (in the sense of double >> precision) vector. You need to write c(1L, 4L, 9L) to ensure you have >> an integer vector or coerce it with as.integer(). >> >>> At this level, behaviour is the same with 'classic' API but we know of at >>> least one case where 'classic' doesn't cast a Matrix right and 'new' >>> does. >>> >>> | so, what happens if there is no coercion? Are the contents of the >>> | SEXP copied to newly allocated storage or not? >>> >>> It depends. In the new API we do fewer copies. >>> >>> | The reason that I ask is because there are some places in the lme4 >>> | code where I want to make sure that I have a pointer to the contents >>> | of the original vector because I am want to change the contents >>> | without copying. >>> >>> I think we do that. If you look at the vignette and the timing example -- >>> but >>> using pointer accessors (eg double *x = foo.begin()) and pointer ops (eg >>> x++) >>> we get essentially the same speed as pure and ugly C code. Because we do >>> away with the copy, and the operator[], for purest speed. >> >> Yes, I like that idea of being able to use foo.begin() to access the >> pointer. Part of my question related to whether I could count on >> getting the pointer to the original contents from R. > > Yes. The objects of the new api have the SEXP as their data member and > methods applied on objects are relayed back to the SEXP in terms of calls t > the R API. We use R memory management all the way. > > There are some times though where memory is reallocated, because there is no > other way. for example if you add values to a vector (push_back, push_front, > etc ...) > > In some way I >> >> would need to check whether the R object had been coerced, and hence >> allocated new storage for its contents, or not. That's not an >> immediate problem, though. > > There is something that might be helpful (or not). In the file RcppCommon.h, > we have the RCPP_RETURN_VECTOR macro, this is an initial attempt to avoid > writing the same switch( TYPEOF(x), case REALSXP: ... code over and over > again. There is an example use in RcppCommon.cpp > > >> Thanks for the response. >> >>> Hope that answer the first batch of questions. Keep'em coming! >>> >>> Dirk > > -- > 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 > > _______________________________________________ Rcpp-devel mailing list [email protected] https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
