Saturday, December 24, 2011

R array to RcppArmadillo cube (updated)

A quick post here on how to convert a 3D array from R to the RcppArmadillo cube.

Rcpp does not support arrays yet (only vectors and matrices). However, RcppArmadillo has a proper 3D array class: arma::cube. The conversion can be done via Rcpp::NumericVector. The only aspect of this exercise that surprised me was the fact that Rcpp::NumericVector stores dimensions of an array, but does not allow one to extract this info. Hence, if you don't know in advance the dimensions of your array, the only solution I see is to pass the vector of dimensions to Rcpp together with the array. Not too ugly, but a bit awkward. Thanks to the commenter Eli below, who pointed out that dimensions of the array can be extracted via the .attr() method -- much prettier now. My code to illustrate this concept:

require(inline)
require(RcppArmadillo)
 
src <- '
  using namespace Rcpp;
 
  NumericVector vecArray(myArray);
  IntegerVector arrayDims = vecArray.attr("dim");
 
  arma::cube cubeArray(vecArray.begin(), arrayDims[0], arrayDims[1], arrayDims[2], false);
 
  //change one element in the array/cube
  cubeArray(0,0,0) = 518;  
 
  return(wrap(cubeArray));  
'
 
readCube = cxxfunction(signature(myArray="numeric"),body=src, plugin="RcppArmadillo")
Created by Pretty R at inside-R.org

A numerical test:

> set.seed(345)
> testArray = array(rnorm(18), dim=c(3,3,2))
> print(testArray)
, , 1
 
           [,1]        [,2]      [,3]
[1,] -0.7849082 -0.29059656 -0.927724
[2,] -0.2795144 -0.06753159  1.710771
[3,] -0.1614579 -0.63352041  1.654769
 
, , 2
 
          [,1]       [,2]       [,3]
[1,]  1.810483 -0.8496292 -1.4029422
[2,]  1.866772  0.3184496  0.5682982
[3,] -1.399833  0.9035913  1.0457561
 
> readCube(testArray)
, , 1
 
            [,1]        [,2]      [,3]
[1,] 518.0000000 -0.29059656 -0.927724
[2,]  -0.2795144 -0.06753159  1.710771
[3,]  -0.1614579 -0.63352041  1.654769
 
, , 2
 
          [,1]       [,2]       [,3]
[1,]  1.810483 -0.8496292 -1.4029422
[2,]  1.866772  0.3184496  0.5682982
[3,] -1.399833  0.9035913  1.0457561
Created by Pretty R at inside-R.org

3 comments:

  1. Thanks Vladimir, it is very useful! I'll try and do something similar with a complex cube (or complex matrices, if it does not work out).
    Robert

    ReplyDelete
  2. Thanks for the post. fyi, vecArray.attr("dim") will get the dims attribute without passing it in.
    Eli

    ReplyDelete
  3. Thanks a bunch for the tip, Eli! Post is updated.

    ReplyDelete