Dimension swap from cube to matrix

[mboeck11]

I encountered a bug in RcppArmadillo when saving a row of a cube to a matrix. If one of the remaining dimensions is just of dimension 1, then dimensions are swapped. I'll illustrate the bug with an example

cppFunction("arma::mat test_row(arma::cube Rarray) {
arma::mat test = Rarray.row(0);
return test;
}", depends="RcppArmadillo")
Rarray1 <- array(NA, dim=c(100,1,2))
dim(test_row(Rarray1))
> 1 2
Rarray2 <- array(NA, dim=c(100,2,1))
dim(test_row(Rarray2))
> 1 2

When pulling out the first row and saving it to a matrix, dimensions are not kept in the second example, but swapped. I'll tested also whether this happens also when pulling out the first column or slice, but there everything seems fine. I'll just attach the code

cppFunction("arma::mat test_col(arma::cube Rarray) {
arma::mat test = Rarray.col(0);
return test;
}", depends="RcppArmadillo")
Rarray3 <- array(NA, dim=c(1,100,2))
dim(test_col(Rarray3))
> 1 2
Rarray4 <- array(NA, dim=c(2,100,1))
dim(test_col(Rarray4))
> 2 1

cppFunction("arma::mat test_slice(arma::cube Rarray) {
arma::mat test = Rarray.slice(0);
return test;
}", depends="RcppArmadillo")

Rarray5 <- array(NA, dim=c(1,2,100))
dim(test_slice(Rarray5))
> 1 2
Rarray6 <- array(NA, dim=c(2,1,100))
dim(test_slice(Rarray6))
> 2 1

[coatless]

Cubes act as a 3d object with (W x H x D). Meanwhile, a matrix is a 2D object with (W x H). The subset operation on a cube is going across the depth dimension for the row/column.

With a debug statement and smaller sample size, we have:

#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::mat test_row(arma::cube Rarray) {
  arma::mat test = Rarray.row(0);
  
  Rcpp::Rcout << "test:" << std::endl << test << std::endl;
  
  return test;
}

View of the data:

Rarray1 <- array(1:10, dim=c(5,1,2))
# , , 1
#
#      [,1]
# [1,]    1
# [2,]    2
# [3,]    3
# [4,]    4
# [5,]    5
#
#, , 2
#
#      [,1]
# [1,]    6
# [2,]    7
# [3,]    8
# [4,]    9
# [5,]   10

Calling the function shows the first and second entry value in the slice being retrieved an placed in sequential order.

dim(test_row(Rarray1))
# test:
#    1.0000   6.0000
# 
# [1] 1 2

I think you want:

#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::mat test_slice_row(arma::cube Rarray) {
  arma::mat test = Rarray.slice(0).row(0);
                                         // ^^^^ Retrieve first matrix in cube
  Rcpp::Rcout << "test:" << std::endl << test << std::endl;
  
  return test;
}

dim(test_slice_row(Rarray1))
# test:
#   1.0000
# 
# [1] 1 1

[eddelbuettel]

Thanks, useful discussion!

Maybe it would help further if either (or preferably both) of the example used an actual data set. Maybe even as simple as array(1:27, dim=c(3,3,3)) ?

[duckmayr]

This is a really interesting issue!

I found it slightly easier to understand, at least for myself (I can be a little thick), the issue @coatless explains by also adding a line printing the whole array, and using the function on two differently dimensioned arrays:

#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::mat test_row(arma::cube Rarray) {
    Rcpp::Rcout << "Rarray:" << std::endl << Rarray << std::endl;
    arma::mat test = Rarray.row(0);
    Rcpp::Rcout << "test:" << std::endl << test << std::endl;
    return test;
}

Then I think the issue with OP's original code becomes quite clear; it takes the first row of every slice!

Rcpp::sourceCpp("so.cpp")

Rarray <- array(1:6, dim = c(3, 1, 2))
test_row(Rarray)
Rarray:
[cube slice 0]
   1.0000
   2.0000
   3.0000

[cube slice 1]
   4.0000
   5.0000
   6.0000


test:
   1.0000   4.0000

     [,1] [,2]
[1,]    1    4

Rarray <- array(1:6, dim = c(3, 2, 1))

test_row(Rarray)
Rarray:
[cube slice 0]
   1.0000   4.0000
   2.0000   5.0000
   3.0000   6.0000


test:
   1.0000   4.0000

     [,1] [,2]
[1,]    1    4

I realize this is what @coatless already said, but the exercise helped me fully understand the issue, so I put it up here in case it helps anyone else understand it also.

[eddelbuettel]

I had been meaning to take another look at this as I think the basic bug report is onto something. If we just pull one dimension we seem be creating a cube from a cube (as the arma-only print shows), but then return a matrix in one case, and vectors in the two others. That seems ... hardly ideal. May need some more though.

Code

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

#include <RcppArmadillo.h>

using namespace arma;

// [[Rcpp::export]]
SEXP getRow(arma::cube cb, int i) {
  cb.print("Cube:");
  auto ro = cb.row(i);
  ro.print("Row:");
  return Rcpp::wrap(ro);
}

// [[Rcpp::export]]
SEXP getCol(arma::cube cb, int i) {
  cb.print("Cube:");
  auto co = cb.col(i);
  co.print("Col:");
  return Rcpp::wrap(co);
}

// [[Rcpp::export]]
SEXP getSlice(arma::cube cb, int i) {
  cb.print("Cube:");
  auto sl = cb.slice(i);
  return Rcpp::wrap(sl);
}

Output

> Rcpp::sourceCpp("/tmp/armacube.cpp")
> cb <- array(1:8, dim=c(2,2,2))
> getRow(cb, 0)                  # this gets us a vector -- should it return a matrix?
Cube:
[cube slice: 0]
   1.0000   3.0000
   2.0000   4.0000

[cube slice: 1]
   5.0000   7.0000
   6.0000   8.0000
Row:
[cube slice: 0]
   1.0000   3.0000

[cube slice: 1]
   5.0000   7.0000
[1] 1 3 5 7
> 
> 
> getCol(cb, 0)                 # same
Cube:
[cube slice: 0]
   1.0000   3.0000
   2.0000   4.0000

[cube slice: 1]
   5.0000   7.0000
   6.0000   8.0000
Col:
[cube slice: 0]
   1.0000
   2.0000

[cube slice: 1]
   5.0000
   6.0000
[1] 1 2 5 6
> 
> getSlice(cb, 0)    # returns a matrix
Cube:
[cube slice: 0]
   1.0000   3.0000
   2.0000   4.0000

[cube slice: 1]
   5.0000   7.0000
   6.0000   8.0000
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> 

That does not seem all that consistent.

[cgiachalis]

You can modify getRow() and getCol() to get a sliced array :

CODE

#include <RcppArmadillo.h>
using namespace arma;

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


// [[Rcpp::export]]
arma::cube getRow2(arma::cube cb, int i) {
  cb.print("Cube:");
  auto ro = cb.row(i);
  ro.print("Row:");
  return ro;
}

// [[Rcpp::export]]
arma::cube  getCol2(arma::cube cb, int i) {
  cb.print("Cube:");
  auto co = cb.col(i);
  co.print("Col:");
  return co;
}

OUTPUT

Rcpp::sourceCpp("cube/armacube2.cpp")
cb <- array(1:8, dim=c(2,2,2))
getCol2(cb, 0) 
Cube:
[cube slice: 0]
   1.0000   3.0000
   2.0000   4.0000

[cube slice: 1]
   5.0000   7.0000
   6.0000   8.0000
Col:
[cube slice: 0]
   1.0000
   2.0000

[cube slice: 1]
   5.0000
   6.0000
, , 1

     [,1]
[1,]    1
[2,]    2

, , 2

     [,1]
[1,]    5
[2,]    6

getRow2(cb, 0) 
Cube:
[cube slice: 0]
   1.0000   3.0000
   2.0000   4.0000

[cube slice: 1]
   5.0000   7.0000
   6.0000   8.0000
Row:
[cube slice: 0]
   1.0000   3.0000

[cube slice: 1]
   5.0000   7.0000
, , 1

     [,1] [,2]
[1,]    1    3

, , 2

     [,1] [,2]
[1,]    5    7

[eddelbuettel]

Yes, that is good. My think may have been to assume the type would change. Your code does not do that -- a row or col select from a cube is still a cube. If that is in fact what one wants then may have a case for a small correction that should make the behavior more type-consistent.