[Rcpp-devel] 3D array proof of concept

[Sparapani, Rodney]

Hi Gang:

For the last several years, I have wanted a way to create 3D arrays on the 
C++ side and return them to R.  For many of my projects, we are doing 
MCMC and we DON'T NEED LINEAR ALGEBRA so we are not using Eigen nor
Armadillo.  Of course, we could use Eigen tensors (but they can't be 
returned to R as of yet) or Armadillo cubes (which apparently can be).  
So, after thinking about this way too long, I have created the following 
proof of concept code which can easily be manipulated to create 3D arrays
and populate them.  Here is a doxygen representation straight from
a rejected Rcpp Gallery entry for your personal enjoyment.

/**
 * @title Proof of concept: creating and returning three-dimensional array
 * @author Rodney Sparapani
 * @license GPL (>= 2)
 * @tags array
 * @summary An example of how to access the entries of an array 
 *   since generic (i, j, k)-like operators don't exist.
 */

#include <Rcpp.h>
// [[Rcpp::export]]
using namespace Rcpp;
IntegerVector get3d(IntegerVector x, IntegerVector args) {
  IntegerVector rc(1), dim=x.attr("dim");
  rc[0]=R_NaN;
  const size_t K=args.size();
  if(K!=dim.size()) return rc;
  size_t i;
  if(K==1) {
    i=args[0]-1;
    if(i>=0 && i<dim[0]) rc[0]=x[i];
  }
  else if(K==2) {
    i=(args[1]-1)*dim[0]+args[0]-1;
    if(i>=0 && i<(dim[0]*dim[1])) rc[0]=x[i];
  }
  else if(K==3) {
    i=((args[2]-1)*dim[1]+(args[1]-1))*dim[0]+args[0]-1;
    if(i>=0 && i<(dim[0]*dim[1]*dim[2])) rc[0]=x[i];
  }
  return rc;
}

/**
 * Occasionally, a 3-dimensional array (or higher) is needed.  For those
 * familiar with Armadillo cubes or Eigen tensors this is rather easy on
 * the C++ side; however, returning this object to R might be difficult.
 * And, for many of us, we are not familiar with these in any case.
 * Furthermore, there is always a hesitancy to add yet another dependency 
 * by having been burned too many times in the past.  These are not
 * evil home-made types: they are present in R and Rcpp by design.
 * Constructing, such arrays with Rcpp is trivial.  However, reading 
 * and writing to their storage locations is not currently supported.
 * This function returns the entry corresponding to "args" from "x" which 
 * is either a 1, 2 or 3-dimensional array.  Since the `(i, j, k)` operator
 * doesn't exist, we resort to "args" which is an integer vector.
 */

/*** R
library(Rcpp)

b = array(1:8, dim=8)
c = array(1:8, dim=c(2, 4))
a <- array(1:24, dim=c(2, 3, 4))

get3d(b, 3)
b[3]
for(k in 1:8) print(b[k]-get3d(b, k))

get3d(c, c(2, 1))
c[2, 1]
for(i in 1:2)
        for(k in 1:4) print(c[i, k]-get3d(c, c(i, k)))

get3d(a, c(2, 1, 1))
a[2, 1, 1]

for(i in 1:2)
    for(j in 1:3)
        for(k in 1:4) 
            print(a[i, j, k]-get3d(a, c(i, j, k)))
*/

-- 

Rodney Sparapani, Associate Professor of Biostatistics
Chair-elect ISBA Section on Biostatistics and Pharmaceutical Statistics
Institute for Health and Equity, Division of Biostatistics
Medical College of Wisconsin, Milwaukee Campus