Initialize a variable with different type based on a switch statement

rcpp

Solution

Thanks to the pointers from Kevin Ushey on separating the dispatching and the function logic. The template code required is sufficiently different to Kevin's suggestions to warrant its own answer.

To write a function that generally works on all types of `big.matrix`, the pattern is as follows:

#include <Rcpp.h>
using namespace Rcpp;

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

// Generic Implementation of GetDiag:
template <typename T>
NumericVector GetDiag_impl(XPtr<BigMatrix> xpMat, MatrixAccessor<T> mat) {
  NumericVector diag(xpMat->nrow()); // Assume matrix is square
  for (unsigned int i = 0; i < xpMat->nrow(); i++) { 
    diag[i] = mat[i][i];
  }
  return diag;
}

// Dispatch code
// [[Rcpp::export]]
NumericVector GetDiag(SEXP pMat) {
  XPtr<BigMatrix> xpMat(pMat);
  switch(xpMat->matrix_type()) {
    case 1:
      return GetDiag_impl(xpMat, MatrixAccessor<char>(*xpMat));
    case 2:
      return GetDiag_impl(xpMat, MatrixAccessor<short>(*xpMat));
    case 4:
      return GetDiag_impl(xpMat, MatrixAccessor<int>(*xpMat));
    case 8:
      return GetDiag_impl(xpMat, MatrixAccessor<double>(*xpMat));
    default:
      // This should be impossible to reach, but shuts up the compiler
      throw Rcpp::exception("Unknown type of big.matrix detected! Aborting.");
  }
}

You do not need to template the return type of `GetDiag_impl`: all four `big.matrix` types are stored as numeric in R (See This Answer for a discussion on 'char' big.matrix objects).

Problem

I am developing some functions in `Rcpp` that operate on `big.matrix` objects from the `bigmemory` package. These objects are passed to `Rcpp` as `SEXP` objects which i then have to cast to an `XPtr<BigMatrix>`, and then to a `MatrixAccessor` object to access elements of the matrix. For example, if I want to implement a function that get's the diagonal: ``` #include <Rcpp.h> using namespace Rcpp; // [[Rcpp::depends(BH, bigmemory) #include <bigmemory/MatrixAccessor.hpp> #include <numeric> // [[Rcpp::export]] NumericVector GetDiag(SEXP pmat) { XPtr<BigMatrix> xpMat(pMat); // allows you to access attributes MatrixAccessor<double> mat(*xpMat); // allows you to access matrix elements NumericVector diag(xpMat->nrow()); // Assume matrix is square for (int i = 0; i < xpMat->nrow(); i++) { diag[i] = mat[i][i]; } return diag; } ``` This function works beautifully, provided the `big.matrix` object in R is filled with doubles. However, if you call this function on an integer matrix (e.g. `diag(as.big.matrix(matrix(1:9, 3))@address)`), You get garbage as a result, because the `MatrixAccessor` has been initialized as `<double>`. Internally, `big.matrix` objects can come in four types: ``` void typeof(SEXP pMat) { XPtr<BigMatrix> xpMat(pMat); int type = xpMat->matrix_type(); type == 1 // char type == 2 // short type == 4 // int type == 8 // double } ``` Since all we're doing is accessing the elements of the matrix, the `diag` function should be able to handle each of these types. But for now, since our function signature is `NumericVector`, I'll ignore character matrices. To handle this, I figured I could just throw in a switch statement, initializing the corresponding `mat` with the appropriate type at runtime: ``` // [[Rcpp::export]] NumericVector GetDiag(SEXP pmat) { XPtr<BigMatrix> xpMat(pMat); // Determine the typeof(pmat), and initialize accordingly: switch(xpMat->matrix_type()) { case == 1: { // Function expects to return a NumericVector. throw; } case == 2: { MatrixAccessor<short> mat(*xpMat); break; } case == 4: { MatrixAccessor<int> mat(*xpMat); break; } case == 8: { MatrixAccessor<double> mat(*xpMat); } } MatrixAccessor<double> mat(*xpMat); // allows you to access matrix elements NumericVector diag(xpMat->nrow()); // Assume matrix is square for (int i = 0; i < xpMat->nrow(); i++) { diag[i] = mat[i][i]; } return diag; } ``` However, this results in compiler errors, because I'm redefining `mat` after it has been declared already in the first `case`. The only way I can see to do this, is to write three different `diag` functions, one for each type, whose code is the same with the exception of the initialization of `mat`. Is there a better way?

Original source

Related problems