[Rcpp-devel] bizarre behavior of RcppMatrix
Dirk Eddelbuettel
edd at debian.org
Tue May 25 04:38:04 CEST 2010
Hi Erica,
Thanks for your interest in Rcpp, and for posting on the list!
On 24 May 2010 at 20:58, Yi-Hsin Erica Tsai wrote:
| Hi,
|
| I've just started using Rcpp with the hopes of speeding up some R loops
| through replacement with C++ code. I've been struggling with
| understanding the behavior of the NumericMatrix datatype and more
| generally how to pass matrices between R and C++.
|
| I seem to be getting bizarre behavior when I pass default parameters to
| my C++ functions. For example, take this C++ function and the
| corresponding R wrapper function:
|
| ______________________________________
|
| //testRcppMatrix.cpp
|
| #include <Rcpp.h>
|
| RcppExport SEXP plusonediag_RcppMatrix(SEXP matrix)
| {
| Rcpp::NumericMatrix orig(matrix);
| Rcpp::NumericMatrix plus(matrix);
I think I once tried something like that and learned that you can get into
trouble with two matrices copied that way -- there is only one underlying
SEXP and pointer, so if you need two distinct ones you can't use it this way
but may need so called 'deep copies' (ie distinct memory).
| int n = orig.rows();
| int k = orig.cols();
|
| for (int i=0; i<n; i++)
| {
| if(i <= k)
| {
| plus(i,i) = plus(i,i) + 1;
| }
| }
|
| Rcpp::Pairlist res(Rcpp::Named( "original", orig),
| Rcpp::Named( "plus", plus)
| );
|
| return res;
| }
| ______________________________________
|
| #testRcpp.r
|
| plusonediag_RcppMatrix_wrapper <- function(mat=matrix(seq(1,9)^2, ncol=3))
| {
| ## Make the call...
| val <- .Call("plusonediag_RcppMatrix",
| mat)
| return(val)
| }
|
| ______________________________________
|
| If I pass actual values to the plusonediag_RcppMatrix_wrapper, then I
What do you mean by 'actual values' ?
| get correct output; however, if I try to use the default parameter (mat)
| or pass a matrix equal to the default matrix, then the original matrix
| appears to be overwritten. See the session output below.
Yes it would be overwritten as there is only one matrix here.
| Anyone have any idea what I'm doing wrong? I'd appreciate any help or
| insight into this problem.
|
| I've attached the test files if that's helpful.
|
| Thanks and best wishes,
| --Erica Tsai
|
| > # main.r
| >
| > library(Rcpp)
| Loading required package: inline
| > dyn.load("testRcppMatrix.so")
| > source("testRcppMatrix.r")
| >
| > foo <- matrix(1:9, 3, 3)
| > bar <- foo ^ 2
| >
| > plusonediag_RcppMatrix_wrapper #default mat is same as bar
| function(mat=matrix(seq(1,9)^2, ncol=3))
| {
| ## Make the call...
| val <- .Call("plusonediag_RcppMatrix",
| mat)
| return(val)
| }
| > plusonediag_RcppMatrix_wrapper() #does not produce correct output!!
|
| $original
| [,1] [,2] [,3]
| [1,] 2 16 49
| [2,] 4 26 64
| [3,] 9 36 82
|
| $plus
| [,1] [,2] [,3]
| [1,] 2 16 49
| [2,] 4 26 64
| [3,] 9 36 82
|
| > plusonediag_RcppMatrix_wrapper(foo) #produces correct output
| $original
| [,1] [,2] [,3]
| [1,] 1 4 7
| [2,] 2 5 8
| [3,] 3 6 9
|
| $plus
| [,1] [,2] [,3]
| [1,] 2 4 7
| [2,] 2 6 8
| [3,] 3 6 10
|
| > plusonediag_RcppMatrix_wrapper(bar) #does not produce correct output!!
| $original
| [,1] [,2] [,3]
| [1,] 2 16 49
| [2,] 4 26 64
| [3,] 9 36 82
|
| $plus
| [,1] [,2] [,3]
| [1,] 2 16 49
| [2,] 4 26 64
| [3,] 9 36 82
|
| > plusonediag_RcppMatrix_wrapper(matrix((1:9), 3, 3)) #produces
| correct output
| $original
| [,1] [,2] [,3]
| [1,] 1 4 7
| [2,] 2 5 8
| [3,] 3 6 9
|
| $plus
| [,1] [,2] [,3]
| [1,] 2 4 7
| [2,] 2 6 8
| [3,] 3 6 10
|
| > plusonediag_RcppMatrix_wrapper(matrix((1:9)^2, 3, 3)) #does not
| produce correct output!!
| $original
| [,1] [,2] [,3]
| [1,] 2 16 49
| [2,] 4 26 64
| [3,] 9 36 82
|
| $plus
| [,1] [,2] [,3]
| [1,] 2 16 49
| [2,] 4 26 64
| [3,] 9 36 82
|
| >
Some of these are indeed goofy!!
I wonder if this may also be related to a bug that Doug Bates reported a few
days ago and which Romain has fixed in SVN -- it had to do with evaluations
of environment related to function calls.
When I do this (using current SVN sources) I get:
> src <- '
+ Rcpp::NumericMatrix orig(matrix);
+ Rcpp::NumericMatrix plus(matrix);
+
+ int n = orig.rows();
+ int k = orig.cols();
+
+ for (int i=0; i<n; i++) {
+ if(i <= k) {
+ plus(i,i) = plus(i,i) + 1;
+ }
+ }
+
+ Rcpp::List res = Rcpp::List::create(Rcpp::Named( "original", orig),
+ Rcpp::Named( "plus", plus));
+ return res;
+ '
>
> fun <- cppfunction(signature(matrix="numeric"), src)
> fun( matrix(1:9, ncol=3))
$original
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
$plus
[,1] [,2] [,3]
[1,] 2 4 7
[2,] 2 6 8
[3,] 3 6 10
> wrapper <- function(mat=matrix(seq(1,9)^2, ncol=3)) { fun(mat) }
> wrapper()
$original
[,1] [,2] [,3]
[1,] 2 16 49
[2,] 4 26 64
[3,] 9 36 82
$plus
[,1] [,2] [,3]
[1,] 2 16 49
[2,] 4 26 64
[3,] 9 36 82
> wrapper(matrix(1:9,ncol=3))
$original
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
$plus
[,1] [,2] [,3]
[1,] 2 4 7
[2,] 2 6 8
[3,] 3 6 10
>
So it 'works' in the first case but not with the wrapper. Odd.
I can't quite make up my mind whether that is really a bug in Rcpp or just my
mind playing tricks with lexical scoping. I mean -- if you don't want your
wrapper to have default arg, can you just not use one?
Dirk
| --
| Yi-Hsin Erica Tsai, Ph.D.
| NSF Postdoctoral Fellow
| Department of Biological Sciences (Carstens lab)
| Louisiana State University
|
| http://www.duke.edu/~yet2/
| etsai at lsu.edu
|
| ----------------------------------------------------------------------
| # main.r
|
| library(Rcpp)
| dyn.load("testRcppMatrix.so")
| source("testRcppMatrix.r")
|
| foo <- matrix(1:9, 3, 3)
| bar <- foo ^ 2
|
| foo
| bar
|
| #"original" matrix is changed
| plusonediag_RcppMatrix_wrapper #default mat is same as bar
| plusonediag_RcppMatrix_wrapper() #does not produce correct output!!
| plusonediag_RcppMatrix_wrapper(foo) #produces correct output
| plusonediag_RcppMatrix_wrapper(bar) #does not produce correct output!!
| plusonediag_RcppMatrix_wrapper(matrix((1:9), 3, 3)) #produces correct output
| plusonediag_RcppMatrix_wrapper(matrix((1:9)^2, 3, 3)) #does not produce correct output!!
|
|
| ----------------------------------------------------------------------
| //testRcppMatrix.cpp
|
| #include <Rcpp.h>
|
| RcppExport SEXP plusonediag_RcppMatrix(SEXP matrix)
| {
| Rcpp::NumericMatrix orig(matrix);
| Rcpp::NumericMatrix plus(matrix);
|
| int n = orig.rows();
| int k = orig.cols();
|
| for (int i=0; i<n; i++)
| {
| if(i <= k)
| {
| plus(i,i) = plus(i,i) + 1;
| }
| }
|
| Rcpp::Pairlist res(Rcpp::Named( "original", orig),
| Rcpp::Named( "plus", plus)
| );
|
| return res;
| }
|
| /*
| RcppExport SEXP plusminusonediag_RcppMatrix(SEXP matrix)
| {
| Rcpp::NumericMatrix orig(matrix);
| Rcpp::NumericMatrix plus(matrix);
| Rcpp::NumericMatrix minus(matrix);
|
| int n = orig.rows();
| int k = orig.cols();
|
| for (int i=0; i<n; i++)
| {
| if(i <= k)
| {
| plus(i,i) = plus(i,i) + 1;
| minus(i,i) = minus(i,i) - 1;
|
| }
| }
|
| Rcpp::Pairlist res(Rcpp::Named( "original", orig),
| Rcpp::Named( "plus", plus),
| Rcpp::Named( "minus", plus)
| );
|
| return res;
| }
| */
|
| ----------------------------------------------------------------------
| #testRcpp.r
|
| plusonediag_RcppMatrix_wrapper <- function(mat=matrix(seq(1,9)^2, ncol=3))
| {
| ## Make the call...
| val <- .Call("plusonediag_RcppMatrix",
| mat)
| return(val)
| }
|
| plusminusonediag_RcppMatrix_wrapper <- function(mat=matrix(seq(1,9)^2, ncol=3))
| {
| ## Make the call...
| val <- .Call("plusminusonediag_RcppMatrix",
| mat)
| return(val)
| }
|
|
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
--
Regards, Dirk
More information about the Rcpp-devel
mailing list