Search code examples
c++rcpprcpparmadillo

Get each row of an arma::mat matrix as arma::vec in for loop


I am using RcppArmadillo to create a function using stochastic simulation. I have trouble pulling out each row of a arma::mat as an arma::vec.

Below is a simplified example of my problem. I have used R nomenclature to illustrate what I am trying to achieve.

I believe there should be a fairly simple way of achieving this in C++, but I'm afraid I haven't figured it out yet. Any help would be much appreciated.

#include <RcppArmadillo.h>   
using namespace Rcpp;

    
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
List function(List params) {
  arma::mat c= params["c"];
  arma::mat I= params["I"];
  
  
  for (int istep = 0; istep < (I.n_elem); istep++) {
    arma::vec loopedrows = I[istep,] //Here I have used the R indexing method, but this does not work in C++

    double product= accu(c*loopedrows)

    arma:vec newvec = stochastic_simulation(product) 

    I[istep+1,] =  newvec // store the output of the in matrix I, again the nomenclature is in R.
  }  

  
  return wrap(I);
};

Solution

  • Here is a simple (and very pedestrian, going step by step in the loop) answer for you.

    Code

    #include <RcppArmadillo.h>
    
    // [[Rcpp::depends(RcppArmadillo)]]
    
    // [[Rcpp::export]]
    arma::mat rowwiseAdd(arma::mat A, arma::mat B) {
        if (A.n_rows != B.n_rows || A.n_cols != B.n_cols)
            Rcpp::stop("Matrices must conform.");
    
        arma::mat C(A.n_rows, A.n_cols);
    
        for (size_t i=0; i < A.n_rows; i++) {
            arma::rowvec a = A.row(i);
            arma::rowvec b = B.row(i);
            arma::rowvec c = a + b;
            C.row(i) = c;
        }
        return C;
    }
    
    /*** R
    A <- matrix(1:9, 3, 3)
    B <- matrix(9:1, 3, 3)
    rowwiseAdd(A, B)
    */
    

    Output

    > Rcpp::sourceCpp("~/git/stackoverflow/70251105/answer.cpp")
    
    > A <- matrix(1:9, 3, 3)
    
    > B <- matrix(9:1, 3, 3)
    
    > rowwiseAdd(A, B)
         [,1] [,2] [,3]
    [1,]   10   10   10
    [2,]   10   10   10
    [3,]   10   10   10
    >