[Rcpp-devel] Wrapping uBlas Matrices into Rcpp Matrices

[Douglas Bates]

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.

>
> 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,