[Rcpp-devel] matrix indexing

[Romain Francois]

On 01/23/2010 09:23 PM, Dirk Eddelbuettel wrote:
>
> On 23 January 2010 at 21:05, Romain Francois wrote:
> | Hello,
> |
> | Motivated by Christophe's question on R-help (cc'ed), I'd like to add
> | facilities for matrix (and maybe arrays) indexing into the new Rcpp api.
> |
> | The classic api already ships matrix indexing through RcppMatrix and
> | RcppMatrixView classes. I'm not a big fan of using round brackets for
> | indexing, but apparently one cannot pass more than one parameter to []
> | so ... I guess I'll have to sacrifice look and feel for a minute.
> |
> |
> | Anyway I was thinking we can maybe add NumericVector::operator()( int,
> | int) (other vectors as well), so that we can do :
> |
> | NumericVector m(x) ;
> | double x11 = m(0,0) ;
> |
> |
> | I was also thinking maybe making some sort of indexer class to take care
> | of conversion of 2 indices into 1d index; something like this pseudo code :
> |
> | NumericVector m(x) ;
> | NumericVector::MatrixIndexer indexer(m) ;
> | double x11 = m[ indexer(0,0) ] ;
> |
> | But this might be even more ugly than operator().
> |
> |
> | Given than a REALSXP can change its mind about being a matrix, I'm not
> | sure I want to subclass NumericVector into NumericMatrix, etc. ...
> |
> | Ideas ?
>
> It's a really important topic, and I really don't know yet what I want.
>
> Recall another r-help questions of not too long ago when the almost-always-
> complaining-yet-never-contributing Ivo Welch went on about needing a
> faster lm().
>
> I once wrote one using GNU gsl for the ordinary least squares part, but one
> can do much_ better with proper C++ classes. As a grad student, I published a
> review paper in the J of Applied Econometrics about Newmat (a C++ library I
> used a lot at the time), and things have come a long way since.
>
> Whit (also CC'ed) and I independently landed at using the clever armadillo
> (templated C++) project at (IIRC) http://arma.sf.net, and there is also the
> Eigen2 project (also templated C++, initially used only inside KDE).
>
> Whit essentially answered Ivo with a fast.lm project on github [ the boy is
> lost in the wilderness and has yet not found the goodness that R-Forge is ;-)
> -- see http://github.com/armstrtw/fast.lm ]
>
> Long story short: given that we pass these objects down to C++, we should
> figure out a fast / light / convenient / elegant / [your term here] / ... way
> to make use of them with other libraries.
>
> Dirk

As I see it, there are two distinct issues:

- convenience api to extract data from a matrix, what should it look 
like: m(i,j) or m[ something(i,j) ] or m[i][j] or whatever else ....

- how to pass data back and forth between Rcpp classes and the rest of 
the world. for this we have as and wrap so I guess the only things 
needed here are :

* converting a SEXP to some c++ class

as< armadillo::Matrix > (SEXP)
as< Newmat::Matrix >(SEXP)

(implicit conversion from RObject to SEXP takes care of the rest), so 
that we can write :

IntegerVector dims = { 4, 4 } ;
NumericVector x(16) ;
x.attr( "dim" ) = dims ;

Newmat::Matrix foo = as< Newmat::Matrix >(x) ;

* The other way is wrap's business, and so would come by overloading 
wrap, as in:

RObject wrap( Newmat::Matrix );
RObject wrap( armadillo::Matrix ) ;

(note that I have never used armadilla, and I only used newmat briefly a 
looooong time ago, so these class names are just made up for the purpose 
of the email)...

Anyway, just to say I don't think Rcpp has to ship these converters, but 
it would maybe fit nicely in another package depending on Rcpp... dunno.

Romain

-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://tr.im/KfKn : Rcpp 0.7.2
|- http://tr.im/JOlc : External pointers with Rcpp
`- http://tr.im/JFqa : R Journal, Volume 1/2, December 2009