[Rcpp-devel] Rcpp: Distinguishing between input types to function call
Søren Højsgaard
sorenh at math.aau.dk
Wed Feb 20 00:15:17 CET 2013
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 this 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
`- http://bit.ly/SweN1Z : SuperStorm Sandy
blog: http://romainfrancois.blog.free.fr
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible
_______________________________________________
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
More information about the Rcpp-devel
mailing list