[Rcpp-devel] Dispatch on sparse and dense matrices...

[Søren Højsgaard]

Thanks, Dirk! 

Follow up: That also means (I guess) that to be on the safe side I should write a check-function that will check that the INTSXP and REALSXP really is a matrix, i.e. that it has a dim attribute; right? Or are there facilities for that?

Cheers
Søren


-----Original Message-----
From: Dirk Eddelbuettel [mailto:edd at debian.org] 
Sent: 8. december 2013 22:59
To: Søren Højsgaard
Cc: rcpp-devel at lists.r-forge.r-project.org
Subject: Re: [Rcpp-devel] Dispatch on sparse and dense matrices...


Hi Søren,

On 8 December 2013 at 21:47, Søren Højsgaard wrote:
| Dear all,
| 
| I want to do some computation on a matrix, irrespective of whether it is a numeric matrix, an integer matrix or a sparse matrix (a dgCMatrix which can be handled with RcppEigen). For simplicity, I want to compute the sum of the elements. 
| 
| To do so I use a template
| 
| template <typename TT>
| SEXP do_compute( SEXP XX_ ){
|   const TT X(as<TT>(XX_));
|   return wrap(X.sum()); // some computation on X...
| };
| 
| I then dispatch depending on input type with
| 
| // [[Rcpp::export]]
| SEXP compute ( SEXP XX_ ){
|   int type = TYPEOF(XX_) ;
|   Rf_PrintValue(wrap(type));
|   switch( type ){
|   case 13 : return do_compute<MapMati>(XX_); // matrix - integer 
|   case 14 : return do_compute<MapMatd>(XX_); // matrix - double
|   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
|   }
|   return R_NilValue ;
| }
| 
| and the numbers 13, 14, 25 are disclosed to by simply printing the type. I know that for integers and doubles I can do something "nicer":
| 
| // [[Rcpp::export]]
| SEXP compute2 ( SEXP XX_ ){
|   int type = TYPEOF(XX_) ;
|   Rf_PrintValue(wrap(type));
|   switch( type ){
|   case INTSXP : return do_compute<MapMati>(XX_); // matrix - integer 
|   case REALSXP : return do_compute<MapMatd>(XX_); // matrix - double
|   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
|   }
|   return R_NilValue ;
| }
| 
| Questions:
| 
| Is there a similar keyword for sparse matrices?? 

I don't think so as sparse matrices do not have a first-class representation in a SEXP. So no TYPEOF(...) that uniquely identifies them.
 
| If not, is the "code" 25 safe to use in an R-package??

No, it simply says 'S4'. You can get that for lots of things which are not Sparse matrices. From Rinternals:

  #define S4SXP       25    /* S4, non-vector */

So you could replace '25' with 'S4SXP' but that does not protect you from
S4 types which are not sparse matrices.

| Are there more elegant ways of handling the dispatch ??

Not sure. The TYPEOF() idiom used to be the standard; I think a few nicer 
C++ wrappers for is<> were added more recently.  Those will still 
C++ dispatch on
base types, and sparse matrices are not base types.  So you may have to go via S4SXP _and_ checking that the required slots exist.  Which you could then encode is a isSpMat() helper ...

Dirk

| The source file is attached. Thanks in advance.
| 
| Søren
| 
| 
| 
| ----------------------------------------------------------------------
| #include <RcppEigen.h>
| //[[Rcpp::depends(RcppEigen)]]
| 
| using namespace Rcpp;
| 
| typedef Eigen::Map<Eigen::MatrixXd> MapMatd; typedef 
| Eigen::Map<Eigen::MatrixXi> MapMati; typedef 
| Eigen::MappedSparseMatrix<double> MSpMat;
| 
| // [[Rcpp::export]]
| SEXP getType ( SEXP XX_ ){
|   int type = TYPEOF(XX_) ;
|   return wrap(type);
| }
| 
| template <typename TT>
| SEXP do_compute( SEXP XX_ ){
|   const TT X(as<TT>(XX_));
|   return wrap(X.sum()); // some computation on X...
| };
| 
| // [[Rcpp::export]]
| SEXP compute ( SEXP XX_ ){
|   int type = TYPEOF(XX_) ;
|   Rf_PrintValue(wrap(type));
|   switch( type ){
|   case 13 : return do_compute<MapMati>(XX_); // matrix - integer 
|   case 14 : return do_compute<MapMatd>(XX_); // matrix - double
|   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
|   }
|   return R_NilValue ;
| }
| 
| // [[Rcpp::export]]
| SEXP compute2 ( SEXP XX_ ){
|   int type = TYPEOF(XX_) ;
|   Rf_PrintValue(wrap(type));
|   switch( type ){
|   case INTSXP : return do_compute<MapMati>(XX_); // matrix - integer 
|   case REALSXP : return do_compute<MapMatd>(XX_); // matrix - double
|   case 25 : return do_compute<MSpMat>(XX_); // dgCMatrix
|   }
|   return R_NilValue ;
| }
| 
| /*** R
| 
| library(Matrix)
| m <- matrix(1:9, nrow=3)
| M <- as(m, "dgCMatrix")
| 
| getType( m ); getType( M ); getType( 1 ); getType( 1L )
| 
| compute( m )   # integer
| compute( 1*m ) # double
| compute( M )   # sparse dgCMatrix
| */
| 
| 
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-deve
| l
--
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com