[Rcpp-devel] Wrapping uBlas Matrices into Rcpp Matrices

[Douglas Bates]

On Tue, Jun 7, 2011 at 11:15 AM, Douglas Bates <bates at stat.wisc.edu> wrote:
> On Mon, Jun 6, 2011 at 6:55 AM, Romain Francois
> <romain at r-enthusiasts.com> wrote:
>> Le 02/06/11 14:43, Cedric Ginestet a écrit :
>>>
>>> Hi again,
>>>
>>> I have tried to do the same for Matrices. Here my naive attempt:
>>> ////////////////////////////////////////////////////////////////
>>> template <typename T>
>>> Rcpp::Matrix< Rcpp::traits::r_sexptype_traits<T>::rtype >
>>> ublas2rcpp( const matrix<T>& x ){
>>> return Rcpp::Matrix< Rcpp::traits::r_sexptype_traits<T>::rtype >(
>>> x.begin(), x.end()
>>> );
>>> }
>>> //////////////////////////////////////////////////////////////
>>> Obviously that doesn't work, and I get the following error message:
>>> templatedFunction.h:63:5: error: ‘const class
>>> boost::numeric::ublas::matrix<int>’ has no member named ‘begin’
>>> templatedFunction.h:63:5: error: ‘const class
>>> boost::numeric::ublas::matrix<int>’ has no member named ‘end’
>>>
>>> I suppose that I either need to 'vectorized' the matrices or to run
>>> through both set of row and column indices. What is the best way to do so?
>>>
>>> Best wishes,
>>> Cedric
>>
>> Again untested, but you might like this constructor from Rcpp::Matrix:
>>
>>        template <typename Iterator>
>>        Matrix( const int& nrows_, const int& ncols, Iterator start ) :
>>            VECTOR( start, start + (nrows_*ncols) ),
>>            nrows(nrows_)
>>        {
>>                VECTOR::attr( "dim" ) = Dimension( nrows, ncols ) ;
>>        }
>>
>>
>> So you'd use it something like this:
>>
>>
>> template <typename T>
>> Rcpp::Matrix< Rcpp::traits::r_sexptype_traits<T>::rtype >
>> ublas2rcpp( const matrix<T>& x ){
>>    return Rcpp::Matrix< Rcpp::traits::r_sexptype_traits<T>::rtype >(
>>        x.size1(), x.size2(), x.begin1()
>>    );
>> }
>>
>> This is untested by guessing what would do the functions found in the
>> documentation of uBlas:
>> http://www.boost.org/doc/libs/1_42_0/libs/numeric/ublas/doc/matrix.htm
>
> I haven't looked much at the dense matrix structures in ublas but the
> sparse matrix templated classes have many arguments with default
> values.  The contents can be represented in different types of vector
> structures for which the default is the ublas::unbounded_array<T>
> structure.  However, you can use a std::vector<T> instead - they just
> feel that their unbounded_array template, which is slightly less
> general than a std::vector template, will be more efficient.
>
> So the first step is either to wrap the different ublas vector
> structures or to use typedef's to specialize the matrix templates to
> use column major dense representations with storage through
> std::vector objects.

If you check the types overview,
http://www.boost.org/doc/libs/1_46_1/libs/numeric/ublas/doc/types_overview.htm,
you will see that you need to specify column-major ordering (default
is row-major) to be able to map a ublas matrix to an R matrix
structure.

>>
>> Romain
>>
>>> On 01/06/11 14:14, Romain Francois wrote:
>>>>
>>>> Le 01/06/11 14:28, Cedric Ginestet a écrit :
>>>>>
>>>>> Dear Romain,
>>>>>
>>>>> Thank you very much for your help. I tried what you suggested by
>>>>> including the following templated function in templatedFunction.h, as
>>>>> follows:
>>>>> template <typename T>
>>>>> Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >
>>>>> ublas2rcpp( const vector<T>& x ){
>>>>> return Rcpp::Vector< r_sexptype_traits<T>::rtype >(
>>>>> x.begin(), x.end()
>>>>> ) ;
>>>>> }
>>>>> In addition, I have tested the function using in subgraph.cpp:
>>>>> Rcpp::Vector<int> xY = ublas2rcpp(Y);
>>>>>
>>>>> And I got the following error messages:
>>>>> templatedFunction.h: In function
>>>>> ‘Rcpp::Vector<Rcpp::traits::r_sexptype_traits<T>::rtype>
>>>>> ublas2rcpp(const boost::numeric::ublas::vector<T>&)’:
>>>>> templatedFunction.h:50:26: error: ‘r_sexptype_traits’ was not declared
>>>>> in this scope
>>>>> templatedFunction.h:50:45: error: template argument 1 is invalid
>>>>> subgraph.cpp: In function ‘SEXPREC* cxx_Mask2Graph(SEXPREC*, SEXPREC*,
>>>>> SEXPREC*, SEXPREC*)’:
>>>>> subgraph.cpp:32:19: error: type/value mismatch at argument 1 in template
>>>>> parameter list for ‘template<int RTYPE> class Rcpp::Vector’
>>>>> subgraph.cpp:32:19: error: expected a constant of type ‘int’, got ‘int’
>>>>> subgraph.cpp:32:24: error: invalid type in declaration before ‘=’ token
>>>>> subgraph.cpp:32:38: error: invalid conversion from ‘SEXPREC*’ to ‘int’
>>>>> subgraph.cpp:34:8: error: invalid conversion from ‘int’ to ‘SEXPREC*’
>>>>> ...
>>>>
>>>> Sure. This was a typo/thinko: go with something like this :
>>>>
>>>> template <typename T>
>>>> Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >
>>>> ublas2rcpp( const vector<T>& x ){
>>>> return Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >(
>>>> x.begin(), x.end()
>>>> ) ;
>>>> }
>>>>
>>>> and Rcpp::Vector<int> makes no sense, you probably want IntegerVector,
>>>> or (the same class):
>>>>
>>>> Rcpp::Vector< r_sexptype_traits<int>::rtype >
>>>>
>>>> Rcpp::Vector is templated by the SEXP type.
>>>>
>>>>> Also, as an aside, I was wondering what I should use instead of
>>>>> push_back for Rcpp Vectors. Do I necessarily have to specify the size of
>>>>> the vector before I assign its elements to specific values?
>>>>
>>>> That is much better yes. ublas probably gives a way to access the size
>>>> of the vector.
>>>>
>>>>> Thanks a lot,
>>>>> Cedric
>>>>>
>>>>>
>>>>> On 01/06/11 11:44, Romain Francois wrote:
>>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I've not used uBlas, but what you are trying to do is quite similar to
>>>>>> what we do in RcppArmadillo.
>>>>>>
>>>>>> You can probably manage to guess the output type from the input type,
>>>>>> so you only have to parameterise your template on the input type.
>>>>>> something like (untested) :
>>>>>>
>>>>>> template <typename T>
>>>>>> Rcpp::Vector< Rcpp::traits::r_sexptype_traits<T>::rtype >
>>>>>> ublas2rcpp( const vector<T>& x ){
>>>>>> return Rcpp::Vector< r_sexptype_traits<T>::rtype >(
>>>>>> x.begin(), x.end()
>>>>>> ) ;
>>>>>> }
>>>>>>
>>>>>> This way you don't have to specify template parameter when you call
>>>>>> ublas2rcpp because the compiler is smart enough.
>>>>>>
>>>>>> Nicer than this would be to implement wrap and as for ublas vectors,
>>>>>> the way to go is documented in the Rcpp-extended vignettes, with
>>>>>> examples implementations in RcppArmadillo and RcppGSL.
>>>>>>
>>>>>> As a side note, you don't want to use push_back on Rcpp types, because
>>>>>> it creates a new vector each time, so this is HUGE memory waste.
>>>>>>
>>>>>> Now, this could get much smarter as ublas has vector expressions,
>>>>>> similar to armadillo, so I suppose someone could write something like
>>>>>> RcppUBlas with nice goodies. This is not me, at least not now ;-)
>>>>>>
>>>>>> Romain
>>>>>>
>>>>>>
>>>>>> Le 01/06/11 12:24, Cedric Ginestet a écrit :
>>>>>>>
>>>>>>> Dear Rcpp experts,
>>>>>>>
>>>>>>> I have started to use the uBlas library, and I am trying to ensure
>>>>>>> that
>>>>>>> I can pass from uBlas vectors to Rcpp vectors relatively easily. So
>>>>>>> far,
>>>>>>> I have tried the following templated function:
>>>>>>>
>>>>>>>
>>>>>>> ///////////////////////////////////////////////////////////////////////
>>>>>>>
>>>>>>> using namespace Rcpp;
>>>>>>> using namespace boost::numeric::ublas;
>>>>>>>
>>>>>>> template <class T1, class T2>
>>>>>>> T1 ublas2rcpp(T2& uVec){
>>>>>>> T1 rcppVec;
>>>>>>> for(int i=0; i<uVec.size(); ++i) rcppVec.push_back(uVec(i));
>>>>>>> return rcppVec;
>>>>>>> }//ublas2rcpp.
>>>>>>>
>>>>>>> ////////
>>>>>>> SEXP foo(vector<int> &y){
>>>>>>> ...
>>>>>>> IntegerVector rcppY=ublas2rcpp<IntegerVector,vector<int> >(y);
>>>>>>> ....
>>>>>>> return rcppY;
>>>>>>> }
>>>>>>>
>>>>>>> ///////////////////////////////////////////////////////////////////////
>>>>>>>
>>>>>>>
>>>>>>> I have got two questions:
>>>>>>> a. This templated function doesn't work. It compiles fine, but hangs
>>>>>>> when I try to run it in R. What do you think is faulty in the codes
>>>>>>> above?
>>>>>>> b. Is there a better way to wrap uBlas vectors into Rcpp ones? Is that
>>>>>>> something that you are planning to implement (or have already
>>>>>>> implemented) within the Rcpp suite?
>>>>>>>
>>>>>>> Thank you very much for your help,
>