[Rcpp-devel] Strange behaviro using Rcpp sugar with NumericMatrix

[Dirk Eddelbuettel]

On 15 March 2017 at 14:17, Roland Denis wrote:
| Hello,
| 
| while comparing the performance between pure R and different other 
| syntaxes that mix R and Rcpp, I found out a strange behavior when using 
| Rcpp sugar on NumericMatrix. As far as I understand Rcpp (I'm a newbie 
| in R and Rcpp), syntactic sugar is designed for NumericVector, but since 
| NumericMatrix seems to be seen like an augmented NumericVector, I've 
| supposed this works on NumericMatrix too (and, in fact, it compiles).
| 
| Here is a minimal example :
| 
| library(Rcpp)
| cppFunction("
|      NumericMatrix scaleSugar( NumericMatrix A, double k )
|      {
|          return k*A;
|      }
| ")
| 
| A <- diag(3)
| B = scaleSugar(A, 2)
| 
| Displaying A and B returns :
|  > A
|       [,1] [,2] [,3]
| [1,]    2    0    0
| [2,]    0    2    0
| [3,]    0    0    2
|  > B
|       [,1] [,2] [,3]
| [1,]    2    0    0
| [2,]    0    2    0
| [3,]    0    0    2
| 
| As you see, the matrix A seems to be multiplied in-place and returned.
| 
| Creating the returned matrix before the multiplication doesn't change 
| the result :
| library(Rcpp)
| cppFunction("
|      NumericMatrix scaleSugar( NumericMatrix A, double k )
|      {
|          NumericMatrix R( A.nrow(), A.ncol());
|          R = k*A;
|          return R;
|      }
| ")
| 
| A <- diag(3)
| B = scaleSugar(A, 2)
| 
|  > A
|       [,1] [,2] [,3]
| [1,]    2    0    0
| [2,]    0    2    0
| [3,]    0    0    2
| 
| However, copying A into R before the scaling works :
| library(Rcpp)
| cppFunction("
|      NumericMatrix scaleSugar( NumericMatrix A, double k )
|      {
|          NumericMatrix R( A.nrow(), A.ncol(), A.begin());
|          return k*R;
|      }
| ")
| 
| A <- diag(3)
| B = scaleSugar(A, 2)
| 
|  > A
|       [,1] [,2] [,3]
| [1,]    1    0    0
| [2,]    0    1    0
| [3,]    0    0    1
| 
| 
| and, *more surprisingly*, reversing the multiplication operand's order 
| works too :
| library(Rcpp)
| cppFunction("
|      NumericMatrix scaleSugar( NumericMatrix A, double k )
|      {
|          return A*k;
|      }
| ")
| 
| A <- diag(3)
| B = scaleSugar(A, 2)
| 
|  > A
|       [,1] [,2] [,3]
| [1,]    1    0    0
| [2,]    0    1    0
| [3,]    0    0    1
| 
| Are Rcpp sugars supposed to work with NumericMatrix ?

That could be a bug.  But please note that we said for many years that the
matrix _calculations_ are underdeveloped / incomplete (see eg the Rcpp FAQ
suggesting to do all this with RcppArmadillo and Armadillo matrices).  So if
I were you I'd take a good look at RcppArmadillo.

Dirk



-- 
http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org