Extracting a column with NA's from a bigmemory

问题:

I'm trying to create a function that extracts a column from a big.matrix object in Rcpp (so that it can be analyzed in cpp before bringing the results to R), but I can't figure out how to get it to recognise NA's (they are now presented as -2147483648 - as shown in my minimal example below). It would be even better if I could access the function GetMatrixCols (src/bigmemory.cpp) straight from Rcpp, but I've yet to discover a way to do that.

#include <Rcpp.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(BH, bigmemory)]]
#include <bigmemory/MatrixAccessor.hpp>
#include <bigmemory/isna.hpp>
using namespace Rcpp;

//Logic for extracting column from a Big Matrix object
template <typename T>
NumericVector GetColumn_logic(XPtr<BigMatrix> pMat,  MatrixAccessor<T> mat,   int cn) {
  NumericVector nv(pMat->nrow());
  for(int i = 0; i < pMat->nrow(); i++) {
    if(isna(mat[cn][i])) {
      nv[i] = NA_INTEGER;
    } else {
      nv[i] = mat[cn][i];
    }
  }
  return nv;
}

//' Extract Column from a Big Matrix.
//' 
//' @param pBigMat A bigmemory object address.
//' @param colNum Column Number to extract. Indexing starts from zero.
//' @export
// [[Rcpp::export]]
NumericVector GetColumn(SEXP pBigMat, int colNum) {
  XPtr<BigMatrix> xpMat(pBigMat);

  switch(xpMat->matrix_type()) {
    case 1: return GetColumn_logic(xpMat, MatrixAccessor<char>(*xpMat), colNum);
    case 2: return GetColumn_logic(xpMat, MatrixAccessor<short>(*xpMat), colNum);
    case 4: return GetColumn_logic(xpMat, MatrixAccessor<int>(*xpMat), colNum);
    case 6: return GetColumn_logic(xpMat, MatrixAccessor<float>(*xpMat), colNum);
    case 8: return GetColumn_logic(xpMat, MatrixAccessor<double>(*xpMat), colNum);
    default: throw Rcpp::exception("Unknown type detected for big.matrix object!");
  }
}

/*** R
bm <- bigmemory::as.big.matrix(as.matrix(reshape2::melt(matrix(c(1:4,NA,6:20),4,5))))
bigmemory:::CGetType(bm@address)
bigmemory:::GetCols.bm(bm, 3)
GetColumn(bm@address, 2)
*/

回答1:

That's a great one! Stay with me for a moment:

tl;dr: It works once fixed:

R> sourceCpp("/tmp/bigmemEx.cpp")

R> bm <- bigmemory::as.big.matrix(as.matrix(reshape2::melt(matrix(c(1:4,NA,6:20),4,5))))

R> bigmemory:::CGetType(bm@address)
[1] 4

R> bigmemory:::GetCols.bm(bm, 3)
 [1]  1  2  3  4 NA  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

R> GetColumn(bm@address, 2)
 [1]  1  2  3  4 NA  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
R> 

The trouble starts at the inside. When you create your matrix as

matrix(c(1:4,NA,6:20),4,5)

what do you get? Integer!

R> matrix(c(1:4,NA,6:20),4,5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1   NA    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
R> class(matrix(c(1:4,NA,6:20),4,5))
[1] "matrix"
R> typeof(matrix(c(1:4,NA,6:20),4,5))
[1] "integer"
R> 

Not a problem per se, but a problem once you remember that the IEEE 754standard has NaN defined for floating point only (correct if I'm wrong).

The other issue is that you reflexively used NumericVector in your, but operate on integers. Now R has NaN, and even NA, for floating point and integer, but 'normal libraries' outside of R do not. And a bigmemory by design represents things outside of R, you're stuck.

The fix is simple enough: use IntegerVector (or equivalently convert your integer data on input). Below is my altered version of your code.

// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

#include <Rcpp.h>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(BH, bigmemory)]]

#include <bigmemory/MatrixAccessor.hpp>
#include <bigmemory/isna.hpp>

using namespace Rcpp;

//Logic for extracting column from a Big Matrix object
template <typename T>
IntegerVector GetColumn_logic(XPtr<BigMatrix> pMat,  MatrixAccessor<T> mat,   int cn) {
    IntegerVector nv(pMat->nrow());
    for(int i = 0; i < pMat->nrow(); i++) {
        if(isna(mat[cn][i])) {
            nv[i] = NA_INTEGER;
        } else {
            nv[i] = mat[cn][i];
        }
    }
    return nv;
}

//' Extract Column from a Big Matrix.
//' 
//' @param pBigMat A bigmemory object address.
//' @param colNum Column Number to extract. Indexing starts from zero.
//' @export
// [[Rcpp::export]]
IntegerVector GetColumn(SEXP pBigMat, int colNum) {
    XPtr<BigMatrix> xpMat(pBigMat);

    switch(xpMat->matrix_type()) {
    case 1: return GetColumn_logic(xpMat, MatrixAccessor<char>(*xpMat), colNum);
    case 2: return GetColumn_logic(xpMat, MatrixAccessor<short>(*xpMat), colNum);
    case 4: return GetColumn_logic(xpMat, MatrixAccessor<int>(*xpMat), colNum);
    case 6: return GetColumn_logic(xpMat, MatrixAccessor<float>(*xpMat), colNum);
    case 8: return GetColumn_logic(xpMat, MatrixAccessor<double>(*xpMat), colNum);
    default: throw Rcpp::exception("Unknown type detected for big.matrix object!");
    }
}

/*** R
bm <- bigmemory::as.big.matrix(as.matrix(reshape2::melt(matrix(c(1:4,NA,6:20),4,5))))
bigmemory:::CGetType(bm@address)
bigmemory:::GetCols.bm(bm, 3)
GetColumn(bm@address, 2)
*/

回答2:

Accessing a column of a big.matrix in Rcpp is not difficult, you can for example get an std vector, an Armadillo vector or an Eigen vector with the following code (there may exist cleaner code):

// [[Rcpp::depends(RcppEigen, RcppArmadillo, bigmemory, BH)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
#include <bigmemory/BigMatrix.h>
#include <bigmemory/MatrixAccessor.hpp>

using namespace Rcpp;
using namespace arma;
using namespace Eigen;
using namespace std;

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]]
ListOf<IntegerVector> AccessVector(SEXP pBigMat, int j) {
  XPtr<BigMatrix> xpMat(pBigMat);
  MatrixAccessor<int> macc(*xpMat);

  int n = xpMat->nrow();

  // Bigmemory
  cout << "Bigmemory:"; 
  for (int i = 0; i < n; i++) {
    cout << macc[j][i] << ' ';
  }
  cout << endl;    

  // STD VECTOR
  vector<int> stdvec(macc[j], macc[j] + n); 

  // ARMA VECTOR
  Row<int> armavec(macc[j], n); // Replace Row by Col if you want

  // EIGEN VECTOR
  VectorXi eigenvec(n);
  memcpy(&(eigenvec(0)), macc[j], n * sizeof(int));

  return(List::create(_["Std vector"] = stdvec, 
                      _["Arma vector"] = armavec,
                      _["Eigen vector"] = eigenvec));
}

AccessVector(bm@address, 2) gets you:

Bigmemory:1 2 3 4 -2147483648 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 
$`Std vector`
 [1]  1  2  3  4 NA  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

$`Arma vector`
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,]    1    2    3    4   NA    6    7    8    9    10    11    12    13    14    15
     [,16] [,17] [,18] [,19] [,20]
[1,]    16    17    18    19    20

$`Eigen vector`
 [1]  1  2  3  4 NA  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

You can see that C doesn't know about NAs but when returning to R, you keep them.

So, it depends on what operations you want to do in Rcpp on the columns. I think if you use directly Eigen or Armadillo operations, it should be OK, but you will certainly get lots of NAs in your result.

Maybe it would be clearer if you say what are these operations you want to do.