[Rcpp-devel] bizarre behavior of RcppMatrix

Tue May 25 17:53:29 CEST 2010

Hi Dirk and Romain,

Thanks for your help.  I see now that I was treating the NumericMatrix 
declarations as copy constructors, but rather should be thinking of them 
as wrapper objects, and I need to be careful to only assign one Rcpp 
variable to each SEXP.  Also, a good catch of the difference between the 
integer and numeric matrix; I would definitely not have spotted that for 
a long time.

Thanks again for the helpful replies.  It's sometimes difficult to be a 
newbie learning a new package, especially when my c++ is so rusty!  :)

Best,
--Erica

On 5/25/2010 2:42 AM, Romain Francois wrote:
> Le 25/05/10 04:38, Dirk Eddelbuettel a écrit :
>>
>>
>> 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).
>
> It depends.
>
> If "matrix" is a numeric matrix, there is only one object. If it is an
> integer matrix, there are 3 objects.
>
> I agree though that wrapping up one SEXP into two distinct Rcpp objects
> is looking for trouble.
>
>> | 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.
>
> Unrelated.
>
>> 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.
>
> All of the examples work, when you have the correct definition of "work"
> ;-)
>
>> 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
>>
>
>