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