[Rcpp-devel] Returning a named (and classed?) list from an RcppExport SEXP function

Douglas Bates bates at stat.wisc.edu
Fri Mar 12 18:23:52 CET 2010

On Fri, Mar 12, 2010 at 10:47 AM, Dirk Eddelbuettel <edd at debian.org> 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.  I will play
around a bit with expressions like

Rcpp::List res(Rcpp::PairList(Rcpp::Named(...), ...))

I would like to have the capability of attaching an attribute named
"class" to the VECSXP before returning it.

> | 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.

> | 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.  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.

Thanks for the response.

> Hope that answer the first batch of questions.  Keep'em coming!
> Dirk
> --
>  Registration is open for the 2nd International conference R / Finance 2010
>  See http://www.RinFinance.com for details, and see you in Chicago in April!

More information about the Rcpp-devel mailing list