CODEWITHSUNDEEP
rrcpp
Barely_sufficient
Barely_sufficient

Reputation: 41

Matrix operations in Rcpp

I'm very new to Rcpp and I would like to manipulate a matrix in Rcpp. I have a 3D matrix in R and I would like to transfer it to a function in Rcpp. I can't seem to do so!

 include <Rcpp.h>
   using namespace Rcpp;



  // [[Rcpp::export]]
  NumericMatrix plusfive(NumericMatrix v) {
    int cols = v.cols();
    int rows = v.rows();
  
    for(int i = 0; i < rows; i++)
    {
      for(int j = 0; j < cols; j++)
      {
        v[i][j] = v[i][j] + 5;
      }
    }
  
    return(v);
  }

  /*** R
  plusfive()
  */

I'm getting the error -

invalid types 'Rcpp::traits::storage_type<14>::type {aka double}[int]' for array subscriptyour text

Upvotes: 2

Views: 73

Answers (1)

Ada Lovelace
Ada Lovelace

Reputation: 1265

Now, matrices are 'by design / by construction' always two-dimensional. Internally, they are represented as a single vector with a dimension attribution of length two -- enforcing just rows and columns. So a 'three-dimensional matrix' cannot exist, and that means you cannot use NumericMatrix.

But Array objects are supported in whichever dimension you want to use, including three. There are examples in the package i.e. this is adapted from one of the tests

Rcpp::cppFunction("IntegerVector x(){return IntegerVector(Dimension(2,3,4));}")

But what you probably really want is to look at RcppArmadillo and its cube class which has (way) more support for three-dimensional arrays right out of the box.

Here is an improved version of your function (with corrected indexing, that creates a new matrix to return to not alter incoming matrix) alongside a simpler variant as the + operator is overloaded for the Matrix type.

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericMatrix plusfive(Rcpp::NumericMatrix v) {
    int cols = v.cols();
    int rows = v.rows();

    Rcpp::NumericMatrix w(rows, cols);
    for(int i = 0; i < rows; i++) {
        for(int j = 0; j < cols; j++) {
            w(i,j) = v(i,j) + 5;
        }
    }
    return(w);
}

// [[Rcpp::export]]
Rcpp::NumericMatrix plusfive2(Rcpp::NumericMatrix v) {
    return v + 5;
}

/*** R
M <- matrix(1.0*(1:4), 2, 2)
plusfive(M)
plusfive2(M)
*/

Upvotes: 5

Related Questions