Skip to content

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

4 messages · Søren Højsgaard, Dirk Eddelbuettel, Kevin Ushey

Original
Søren Højsgaard
# 
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?? 

If not, is the "code" 25 safe to use in an R-package??

Are there more elegant ways of handling the dispatch ??

The source file is attached. Thanks in advance.

S?ren


-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: dispatch.cpp
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20131208/6845e470/attachment.ksh>
Dirk Eddelbuettel
# 
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 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-devel
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
Kevin Ushey
# 
I think the safest bet would be to handle dispatch at the R level, which
would then call your internal C++ implementations. At least, it would be
easier / cleaner. Otherwise, you're right, you need to perform an 'ad-hoc'
dispatch based on the internal type / class and such.

-Kevin
On Sun, Dec 8, 2013 at 2:08 PM, S?ren H?jsgaard <sorenh at math.aau.dk> wrote:

            

      
      
      
    

        

      
      
    

        
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20131208/594996b1/attachment.html>