[Rcpp-devel] Rcpp: Distinguishing between input types to function call

Søren Højsgaard sorenh at math.aau.dk
Wed Feb 20 10:05:05 CET 2013


Thanks a lot! 

I thought that Rf_inherits was the Rcpp version of inherits (which gives what I expect).

> XX1 <- letters[1:4] # character
> XX2 <- 1:4          # integer
> XX3 <- (1:4)+.5     # numeric
> inherits(XX1, "character")
[1] TRUE
> inherits(XX1, "numeric")
[1] FALSE
> inherits(XX1, "integer")
[1] FALSE
> inherits(XX2, "character")
[1] FALSE
> inherits(XX2, "numeric")
[1] FALSE
> inherits(XX2, "integer")
[1] TRUE
> inherits(XX3, "character")
[1] FALSE
> inherits(XX3, "numeric")
[1] TRUE
> inherits(XX3, "integer")
[1] FALSE

Best regards
Søren

-----Original Message-----
From: Romain Francois [mailto:romain at r-enthusiasts.com] 
Sent: 20. februar 2013 09:53
To: Søren Højsgaard
Cc: rcpp-devel at lists.r-forge.r-project.org
Subject: Re: [Rcpp-devel] Rcpp: Distinguishing between input types to function call

Hello,

Here is a shorter version of your code. The key idea was to use TYPEOF instead of Rf_inherits which uses the class attribute (simple vectors don't have them).

#include <Rcpp.h>
using namespace Rcpp;

template <int RTYPE>
SEXP allpairsXtemplate_( SEXP XX_ ){
   Vector<RTYPE> X(XX_);
   Matrix<RTYPE> ans(2, X.size()*(X.size()-1)/2);
   int col=0;
   for (int ii=0; ii<X.size(); ii++){
     for (int jj=ii+1; jj<X.size(); jj++){
       ans(0,col) = X(ii);
       ans(1,col++) = X(jj);
     }
   }
   return ans ;
};

// [[Rcpp::export]]
SEXP allpairsX_ ( SEXP XX_ ){
     int type = TYPEOF(XX_) ;
     switch( type ){
         case INTSXP : return allpairsXtemplate_<INTSXP> ( XX_ ) ;
         case REALSXP: return allpairsXtemplate_<REALSXP>( XX_ ) ;
         case STRSXP : return allpairsXtemplate_<STRSXP> ( XX_ ) ;
     }
     return R_NilValue ;
}


/*** R

XX1 <- letters[1:4] # character
XX2 <- 1:4          # integer
XX3 <- (1:4)+.5     # numeric

allpairsX_( XX1 )
allpairsX_( XX2 )
allpairsX_( XX3 )

***/

Also, I'm templating allpairsXtemplate_ on the R type rather than the actual classes, because NumericVector = Vector<REALSXP>, etc ...


About your code, with e.g. TT = NumericVector, you don't need as in :

TT X = as<TT>(XX_);

because NumericVector already has a SEXP constructor, that is why I do: 
Vector<RTYPE> X(XX_);



Same for return(wrap(ans));  you don't need to call wrap here because 
ans can convert itself to SEXP.



Another way to write this using Rcpp's builtin dispatch mechanism is to 
use RCPP_RETURN_VECTOR. For example :

#include <Rcpp.h>
using namespace Rcpp;

template <typename T>
SEXP allpairsXtemplate_( const T& X){
   const int RTYPE = T::r_type::value ;
   Matrix<RTYPE> ans(2, X.size()*(X.size()-1)/2);
   int col=0;
   for (int ii=0; ii<X.size(); ii++){
     for (int jj=ii+1; jj<X.size(); jj++){
       ans(0,col) = X(ii);
       ans(1,col++) = X(jj);
     }
   }
   return ans ;
};

// [[Rcpp::export]]
SEXP allpairsX_ ( SEXP XX_ ){
     RCPP_RETURN_VECTOR( allpairsXtemplate_, XX_ ) ;
     return R_NilValue ; // never used
}


So we call one of the generated overloads of allpairsXtemplate_ which 
takes a Vector as input. From this vector, we can deduce the RTYPE (at 
compile time):

const int RTYPE = T::r_type::value ;

use it to get the correct Matrix type : Matrix<RTYPE>.



Yet another way, probably the one I would use:

template <int RTYPE>
Matrix<RTYPE> allpairsXtemplate_( const Vector<RTYPE>& X){
   Matrix<RTYPE> ans(2, X.size()*(X.size()-1)/2);
   int col=0;
   for (int ii=0; ii<X.size(); ii++){
     for (int jj=ii+1; jj<X.size(); jj++){
       ans(0,col) = X(ii);
       ans(1,col++) = X(jj);
     }
   }
   return ans ;
};

This works because RCPP_RETURN_VECTOR will cast to the appropriate 
Vector type.

And knowing the RTYPE at first lets us use it on the output signture.


RCPP_RETURN_VECTOR is defined in dispatch.h (macro haters beware):

#define ___RCPP_HANDLE_CASE___( ___RTYPE___ , ___FUN___ , ___OBJECT___ , 
___RCPPTYPE___ )	\
	case ___RTYPE___ :																	\
		return ___FUN___( ::Rcpp::___RCPPTYPE___< ___RTYPE___ >( 
___OBJECT___ ) ) ;	
		
#define ___RCPP_RETURN___( __FUN__, __SEXP__ , __RCPPTYPE__ )						\
	SEXP __TMP__ = __SEXP__ ;															\
	switch( TYPEOF( __TMP__ ) ){										 				\
		___RCPP_HANDLE_CASE___( INTSXP  , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( REALSXP , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( RAWSXP  , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( LGLSXP  , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( CPLXSXP , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( STRSXP  , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( VECSXP  , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
		___RCPP_HANDLE_CASE___( EXPRSXP , __FUN__ , __TMP__ , __RCPPTYPE__ )			\
	default:																			\
		throw std::range_error( "not a vector" ) ;									\
	}

#define RCPP_RETURN_VECTOR( _FUN_, _SEXP_ )  ___RCPP_RETURN___( _FUN_, 
_SEXP_ , Vector )
#define RCPP_RETURN_MATRIX( _FUN_, _SEXP_ )  ___RCPP_RETURN___( _FUN_, 
_SEXP_ , Matrix )



Romain

Le 20/02/13 00:15, Søren Højsgaard a écrit :
> Dear all
>
> I have tried to follow Romains suggestion (thanks) below to obtain all pairs of elements of a vector, for various input types; i.e.
>
>   XX1 <- letters[1:4] # character
>   XX2 <- 1:4          # integer
>   XX3 <- (1:4)+.5     # numeric
>   combn(XX1, 2)
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,] "a"  "a"  "a"  "b"  "b"  "c"
> [2,] "b"  "c"  "d"  "c"  "d"  "d"
>   combn(XX2, 2)
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]    1    1    1    2    2    3
> [2,]    2    3    4    3    4    4
>   combn(XX3, 2)
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]  1.5  1.5  1.5  2.5  2.5  3.5
> [2,]  2.5  3.5  4.5  3.5  4.5  4.5
>
> My take on this is as follows:
> ------------------------------
> #include <Rcpp.h>
>
> #ifndef BEGIN_RCPP
> #define BEGIN_RCPP
> #endif
>
> #ifndef END_RCPP
> #define END_RCPP
> #endif
>
> using namespace Rcpp;
>
> // [[Rcpp::export]]
> template <typename TT, typename UU>
> SEXP allpairsXtemplate_( SEXP XX_ ){
>    TT X = as<TT>(XX_);
>    UU ans(2, X.size()*(X.size()-1)/2);
>    int col=0;
>    for (int ii=0; ii<X.size(); ii++){
>      for (int jj=ii+1; jj<X.size(); jj++){
>        ans(0,col) = X(ii);
>        ans(1,col++) = X(jj);
>      }
>    }
>    return(wrap(ans));
> };
>
> // [[Rcpp::export]]
> RcppExport SEXP allpairsX_char ( SEXP XX_ ){
>    return allpairsXtemplate_<CharacterVector, CharacterMatrix>(XX_);
> }
>
> // [[Rcpp::export]]
> RcppExport SEXP allpairsX_int ( SEXP XX_ ){
>    return allpairsXtemplate_<IntegerVector, IntegerMatrix>(XX_);
> }
>
> // [[Rcpp::export]]
> RcppExport SEXP allpairsX_num ( SEXP XX_ ){
>    return allpairsXtemplate_<NumericVector, NumericMatrix>(XX_);
> }
>
> // [[Rcpp::export]]
> RcppExport SEXP allpairsX_ ( SEXP XX_ ){
>    if( Rf_inherits( XX_, "character" ) ){
>      Rcout << "character\n";
>      return allpairsXtemplate_<CharacterVector, CharacterMatrix>(XX_);
>    }
>    if (Rf_inherits( XX_, "integer" ) ){
>      Rcout << "integer\n";
>      return allpairsXtemplate_<IntegerVector, IntegerMatrix>(XX_);
>    }
>    if (Rf_inherits( XX_, "numeric" ) ){
>      Rcout << "numeric\n";
>      return allpairsXtemplate_<NumericVector, NumericMatrix>(XX_);
>    }
>    return R_NilValue;
> }
>
> ------------------------------
>
> I correctly get:
>
> dyn.load("template.dll")
>   .Call("allpairsX_char", XX1)
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,] "a"  "a"  "a"  "b"  "b"  "c"
> [2,] "b"  "c"  "d"  "c"  "d"  "d"
>   .Call("allpairsX_int",  XX2)
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]    1    1    1    2    2    3
> [2,]    2    3    4    3    4    4
>   .Call("allpairsX_num",  XX3)
>       [,1] [,2] [,3] [,4] [,5] [,6]
> [1,]  1.5  1.5  1.5  2.5  2.5  3.5
> [2,]  2.5  3.5  4.5  3.5  4.5  4.5
>   dyn.unload("template.dll")
>
> However the function allpairsX_ fails:
>
> dyn.load("template.dll")
>   .Call("allpairsX_", XX1)
> NULL
>   .Call("allpairsX_", XX2)
> NULL
>   .Call("allpairsX_", XX3)
> NULL
>   dyn.unload("template.dll")
>
> Now for the questions:
>
> 1) From various tests it seems that Rf_inherits does not work - or perhaps I have misunderstood its usage. Any experiences with that? Any other suggestions on how to dispatch on the input type?
>
> 2) I have never used templates before. Is the approach above what "one would normally do"?
>
> 3) Using sourceCpp I get the following:
>
> sourceCpp("template.cpp")
> g++ -m64 -I"C:/programs/R/current/include" -DNDEBUG     -I"C:/programs/R/current/library/Rcpp/include" -I"d:/RCompile/CRANpkg/extralibs64/local/include"     -O2 -Wall  -mtune=core2 -c template.cpp -o template.o template.cpp: In function 'SEXPREC* sourceCpp_55702_allpairsXtemplate_(SEXP)': template.cpp:67:5: error: a template declaration cannot appear at block scope template.cpp:68:5: error: expected ';' before 'return' template.cpp:66:10: warning: unused variable 'XX_' [-Wunused-variable] template.cpp: In function 'SEXPREC* sourceCpp_55702_allpairsX_char(SEXP)': template.cpp:76:5: error: expected unqualified-id before string constant template.cpp:77:23: error: '__result' was not declared in this scope template.cpp:75:10: warning: unused variable 'XX_' [-Wunused-variable] template.cpp: In function 'SEXPREC* sourceCpp_55702_allpairsX_int(SEXP)': template.cpp:85:5: error: expected unqualified-id before string constant template.cpp:86:23: error: '__result' was not declared in t
 his scope
 template.cpp:84:10: warning: unused variable 'XX_' [-Wunused-variable] template.cpp: In function 'SEXPREC* sourceCpp_55702_allpairsX_num(SEXP)': template.cpp:94:5: error: expected unqualified-id before string constant template.cpp:95:23: error: '__result' was not declared in this scope template.cpp:93:10: warning: unused variable 'XX_' [-Wunused-variable] template.cpp: In function 'SEXPREC* sourceCpp_55702_allpairsX_(SEXP)': template.cpp:103:5: error: expected unqualified-id before string constant template.cpp:104:23: error: '__result' was not declared in this scope template.cpp:102:10: warning: unused variable 'XX_' [-Wunused-variable] make: *** [template.o] Error 1
> Error in sourceCpp("template.cpp") :
>    Error 1 occurred building shared library.
>
> It seems that the error occurs because of the template. Am I doing something wrong or is it just not possible to use sourceCpp when templates are involved.
>
> 4) Not a question, but an observation: On windows the above error message comes as one long line which means that I must manually scroll to the end of the line ("to the far right"). Slightly annoying. For comparison, cxxfunction() produces error in more readable form: One line per error. It would be nice if sourceCpp did the same thing.
>
> Thanks in advance - and thanks for making Rcpp available.
>
> Best regards
> Søren
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> -----Original Message-----
> From: rcpp-devel-bounces at lists.r-forge.r-project.org [mailto:rcpp-devel-bounces at lists.r-forge.r-project.org] On Behalf Of Romain Francois
> Sent: 3. december 2012 23:18
> To: rcpp-devel at lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] Rcpp: Distinguishing between input types to function call
>
>
> Hello,
>
> I have not tested this, but I think you are looking for templates. So you'd put generic code in this template function:
>
> template <typename T>
> SEXP topoSort( SEXP XX_ ){
>     const T X = Rcpp::as<T>(XX_) ;
>
>     ...
> }
>
>
> and two other functions to instantiate the template:
>
> RcppExport SEXP C_topoSort_st ( SEXP XX_ ){
>       return topoSort< Eigen::Map<Eigen::MatrixXi> >( XX_ ) ; }
>
> RcppExport SEXP C_topoSort_sp ( SEXP XX_ ){
>       return topoSort< Eigen::MappedSparseMatrix<double> >( XX_ ) ; }
>
> or Perhaps you would have one instead of the two, something like this:
>
> RcppExport SEXP topoSort_facade ( SEXP XX_ ){
>
> 	if( Rf_inherits( XX_, "dgCMatrix" ) ){
> 		return topoSort< Eigen::MappedSparseMatrix<double> >( XX_ ) ;
> 	} else {
> 		return topoSort< Eigen::Map<Eigen::MatrixXi> >( XX_ ) ;
> 	}
> }
>
>
>
>
> Le 03/12/12 22:58, Søren Højsgaard a écrit :
>> Dear list,
>>
>> I represent a directed acyclic graph (DAG) as an adjacency matrix. This can be either a "standard matrix" in R or as a sparse matrix (dgCMatrix from the matrix package). I have implemented a topological sort function for DAGs for these two representations (using the RcppEigen package):
>>
>> // standard matrix
>> RcppExport SEXP C_topoSort_st ( SEXP XX_ ){
>>     typedef Eigen::Map<Eigen::MatrixXi> MapMati;
>>     const MapMati X(Rcpp::as<MapMati>(XX_));
>>     //typedef Eigen::MappedSparseMatrix<double> MSpMat;
>>     //const MSpMat   X(as<MSpMat>(XX_));
>>     .... some code
>> }
>>
>> // sparse matrix
>> RcppExport SEXP C_topoSort_sp ( SEXP XX_ ){
>>     //   typedef Eigen::Map<Eigen::MatrixXi> MapMati;
>>     //   const MapMati X(Rcpp::as<MapMati>(XX_));
>>     typedef Eigen::MappedSparseMatrix<double> MSpMat;
>>     const MSpMat   X(as<MSpMat>(XX_));
>>     .... some code
>> }
>>
>> Notice: The functions only differ with respect to the first four lines.
>>
>> Question: Is there any way in which I can "reduce" these two functions to only one which then checks the "type" of XX_ at the entry and then creates the appropriate "type" of X?
>
> Templates.
>
>> Question: Is it correct (haven't tried, just guessing from what I've read) that I can not directly store 'some code' in an inline function (because the correct type of X would need to be known?
>
> I think I understand what you mean, and that you are fine.
>
>> Apologies for trivial C++ questions - I am working on learning it...
>
> Those are good kind of questions to ask yourself. I hope this will give you enough motivation to find out more about C++ templates.
>
> Romain
>
>> The functions are listed below.
>>
>> Best regards
>> Søren
>>
>> ----------------------------------
>>
>> # include <RcppEigen.h>
>> # include <Rcpp.h>
>>
>> #ifndef BEGIN_RCPP
>> #define BEGIN_RCPP
>> #endif
>>
>> #ifndef END_RCPP
>> #define END_RCPP
>> #endif
>>
>> using namespace Rcpp;
>>
>> // standard matrix
>> RcppExport SEXP C_topoSort_st ( SEXP XX_ ){
>>     typedef Eigen::Map<Eigen::MatrixXi> MapMati;
>>     const MapMati X(Rcpp::as<MapMati>(XX_));
>>     //typedef Eigen::MappedSparseMatrix<double> MSpMat;
>>     //const MSpMat   X(as<MSpMat>(XX_));
>>     int ii, jj, kk=0, count=0, ll=0, flagsum=0;
>>     int ncX(X.rows());
>>     Eigen::VectorXi indegree(ncX);
>>     Eigen::VectorXi flag(ncX);
>>     Eigen::VectorXi ans(ncX);
>>
>>     for (ii = 0; ii < ncX; ii++) {
>>       indegree[ii] = 0; flag[ii] = 0; ans[ii] = 0;
>>     }
>>     for (jj = 0; jj < ncX; jj++)
>>       for (ii = 0; ii < ncX; ii++)
>>         indegree[jj] = indegree[jj] +  X.coeff(ii,jj);
>>
>>     /*   Rcout<<"indegree: ";for (ii=0;ii<ncX;ii++) Rcout << indegree[ii]<<" " ; Rcout << std::endl;*/
>>     /*   Rcout<<"flag    : ";for (ii=0;ii<ncX;ii++) Rcout << flag[ii]<<" " ; Rcout << std::endl;*/
>>     while (count < ncX){
>>       /* Rcout << "count=" << count << std::endl;*/
>>       for (kk = 0; kk < ncX; kk++){
>>         /* Rcout <<" kk="<<kk<<" indeg="<<indegree[kk]<<" flag="<<flag[kk] << std::endl;*/
>>         if ((indegree[kk] == 0) && (flag[kk] == 0)){
>> 	/*Rcout << "   no incomming:" << kk << std::endl;*/
>> 	ans[ll++] = kk+1;
>> 	flag[kk]  = 1;
>> 	flagsum++;
>> 	for (jj = 0; jj < ncX; jj++){
>> 	  /*  Rcout <<"kk,jj="<<kk<<","<<jj<<" entry=" << X.coeff(kk,jj) << std::endl;*/
>> 	  if (X.coeff(kk,jj) == 1){
>> 	    indegree[jj]--;
>> 	    /* Rcout <<" updating indegree at entry="<<jj<<std::endl;*/
>> 	  }
>> 	}
>>         }
>>         /* Rcout<<"indegree: ";for (ii=0;ii<ncX;ii++) Rcout << indegree[ii]<<" " ; Rcout << std::endl;	*/
>>       }
>>       if (flagsum==ncX)
>>         break;
>>       count++;
>>       /* Rcout<<"flag    : ";for (ii=0;ii<ncX;ii++) Rcout << flag[ii]<<" " ; Rcout << std::endl;	*/
>>     }
>>     if (flagsum<ncX)
>>       ans[0] = -1;
>>     return(wrap(ans));
>> }
>>
>> // sparse matrix
>> RcppExport SEXP C_topoSort_sp ( SEXP XX_ ){
>>     //   typedef Eigen::Map<Eigen::MatrixXi> MapMati;
>>     //   const MapMati X(Rcpp::as<MapMati>(XX_));
>>     typedef Eigen::MappedSparseMatrix<double> MSpMat;
>>     const MSpMat   X(as<MSpMat>(XX_));
>>     int ii, jj, kk=0, count=0, ll=0, flagsum=0;
>>     int ncX(X.rows());
>>     Eigen::VectorXi indegree(ncX);
>>     Eigen::VectorXi flag(ncX);
>>     Eigen::VectorXi ans(ncX);
>>
>>     for (ii = 0; ii < ncX; ii++) {
>>       indegree[ii] = 0; flag[ii] = 0; ans[ii] = 0;
>>     }
>>     for (jj = 0; jj < ncX; jj++)
>>       for (ii = 0; ii < ncX; ii++)
>>         indegree[jj] = indegree[jj] +  X.coeff(ii,jj);
>>
>>     /*   Rcout<<"indegree: ";for (ii=0;ii<ncX;ii++) Rcout << indegree[ii]<<" " ; Rcout << std::endl;*/
>>     /*   Rcout<<"flag    : ";for (ii=0;ii<ncX;ii++) Rcout << flag[ii]<<" " ; Rcout << std::endl;*/
>>     while (count < ncX){
>>       /* Rcout << "count=" << count << std::endl;*/
>>       for (kk = 0; kk < ncX; kk++){
>>         /* Rcout <<" kk="<<kk<<" indeg="<<indegree[kk]<<" flag="<<flag[kk] << std::endl;*/
>>         if ((indegree[kk] == 0) && (flag[kk] == 0)){
>> 	/*Rcout << "   no incomming:" << kk << std::endl;*/
>> 	ans[ll++] = kk+1;
>> 	flag[kk]  = 1;
>> 	flagsum++;
>> 	for (jj = 0; jj < ncX; jj++){
>> 	  /*  Rcout <<"kk,jj="<<kk<<","<<jj<<" entry=" << X.coeff(kk,jj) << std::endl;*/
>> 	  if (X.coeff(kk,jj) == 1){
>> 	    indegree[jj]--;
>> 	    /* Rcout <<" updating indegree at entry="<<jj<<std::endl;*/
>> 	  }
>> 	}
>>         }
>>         /* Rcout<<"indegree: ";for (ii=0;ii<ncX;ii++) Rcout << indegree[ii]<<" " ; Rcout << std::endl;	*/
>>       }
>>       if (flagsum==ncX)
>>         break;
>>       count++;
>>       /* Rcout<<"flag    : ";for (ii=0;ii<ncX;ii++) Rcout << flag[ii]<<" " ; Rcout << std::endl;	*/
>>     }
>>     if (flagsum<ncX)
>>       ans[0] = -1;
>>     return(wrap(ans));
>> }


-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30

R Graph Gallery: http://gallery.r-enthusiasts.com

blog:            http://romainfrancois.blog.free.fr
|- http://bit.ly/14LJhmm : bibtex 0.3-5
`- http://bit.ly/RE6sYH  : OOP with Rcpp modules



More information about the Rcpp-devel mailing list