Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Index element from list in Rcpp

Tags:

r

rcpp

Suppose I have a List in Rcpp, here called x containing matrices. I can extract one of the elements using x[0] or something. However, how do I extract a specific element of that matrix? My first thought was x[0](0,0) but that does not seem to work. I tried using * signs but also doesn't work.

Here is some example code that prints the matrix (shows matrix can easily be extracted):

library("Rcpp")

cppFunction(
includes = ' 
NumericMatrix RandMat(int nrow, int ncol)
 {
  int N = nrow * ncol;
  NumericMatrix Res(nrow,ncol);
  NumericVector Rands  = runif(N);
   for (int i = 0; i < N; i++) 
  {
    Res[i] = Rands[i];
  }
  return(Res);
 }',

code = '
void foo()
{
  List x;
  x[0] = RandMat(3,3);
  Rf_PrintValue(wrap( x[0] )); // Prints first matrix in list.
}
')


foo()

How could I change the line Rf_PrintValue(wrap( x[0] )); here to print the the element in the first row and column? In the code I want to use it for I need to extract this element to do computations.

like image 736
Sacha Epskamp Avatar asked Jul 31 '13 14:07

Sacha Epskamp


2 Answers

Quick ones:

  1. Compound expression in C++ can bite at times; the template magic gets in the way. So just assign from the List object to a whatever the element is, eg a NumericMatrix.

  2. Then pick from the NumericMatrix as you see fit. We have row, col, element, ... access.

  3. Printing can be easier using Rcpp::Rcout << anElement but note that we currently cannot print entire matrices or vectors -- but the int or double types are fine.

Edit:

Here is a sample implementation.

#include <Rcpp.h>

// [[Rcpp::export]]
double sacha(Rcpp::List L) {
    double sum = 0;
    for (int i=0; i<L.size(); i++) {
        Rcpp::NumericMatrix M = L[i];
        double topleft = M(0,0);
        sum += topleft;
        Rcpp::Rcout << "Element is " << topleft << std::endl;
    }
    return sum;    
}

/*** R
set.seed(42)
L <- list(matrix(rnorm(9),3), matrix(1:9,3), matrix(sqrt(1:4),2))
sacha(L) # fix typo   
*/

And its result:

R> Rcpp::sourceCpp('/tmp/sacha.cpp')

R> set.seed(42)

R> L <- list(matrix(rnorm(9),3), matrix(1:9,3), matrix(sqrt(1:4),2))

R> sacha(L)
Element is 1.37096
Element is 1
Element is 1
[1] 3.37096
R>
like image 87
Dirk Eddelbuettel Avatar answered Oct 17 '22 10:10

Dirk Eddelbuettel


You have to be explicit at some point. The List class has no idea about the types of elements it contains, it does not know it is a list of matrices.

Dirk has shown you what we usually do, fetch the element as a NumericMatrix and process the matrix.

Here is an alternative that assumes that all elements of your list have the same structure, using a new class template: ListOf with enough glue to make the user code seamless. This just moves to a different place the explicitness.

#include <Rcpp.h>
using namespace Rcpp ;

template <typename WHAT>
class ListOf : public List {
public:
    template <typename T>
    ListOf( const T& x) : List(x){}

    WHAT operator[](int i){ return as<WHAT>( ( (List*)this)->operator[]( i) ) ; }

} ;

// [[Rcpp::export]]
double sacha( ListOf<NumericMatrix> x){
    double sum = 0.0 ; 
    for( int i=0; i<x.size(); i++){
        sum += x[i](0,0) ;    
    }
    return sum ;
}

/*** R
    L <- list(matrix(rnorm(9),3), matrix(1:9,3), matrix(sqrt(1:4),2))
    sacha( L )
*/

When I sourceCpp this file, I get:

> L <- list(matrix(rnorm(9), 3), matrix(1:9, 3), matrix(sqrt(1:4), 2))    
> sacha(L)
[1] 1.087057
like image 6
Romain Francois Avatar answered Oct 17 '22 10:10

Romain Francois