[Rcpp-devel] Rcpp: Distinguishing between input types to function call
Romain Francois
romain at r-enthusiasts.com
Thu Feb 21 08:11:18 CET 2013
Hello,
Soren and I had the same idea, so we'll jointly write it for the Rcpp
gallery.
Romain
Le 20/02/13 17:55, Kevin Ushey a écrit :
> Hi Romain,
>
> Just wanted to give you a thanks for putting together this answer. It's
> a really great example of more generic programming with Rcpp by
> templating over the RTYPEs.
>
> You should consider submitting some version of this to the Rcpp gallery;
> I think it has a lot of use cases.
>
> -Kevin
>
> On Wed, Feb 20, 2013 at 12:53 AM, Romain Francois
> <romain at r-enthusiasts.com <mailto:romain at r-enthusiasts.com>> wrote:
>
> 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
> <tel:2.5%20%203.5%20%204.5%20%203.5%20%204.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
> <tel:2.5%20%203.5%20%204.5%20%203.5%20%204.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>
> [mailto: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
> <mailto: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