[Rcpp-devel] bizarre behavior of RcppMatrix
Romain Francois
romain at r-enthusiasts.com
Tue May 25 09:42:53 CEST 2010
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
>
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/cork4b : highlight 0.1-8
|- http://bit.ly/bklUXt : RcppArmadillo 0.2.1
`- http://bit.ly/936ck2 : Rcpp 0.8.0
More information about the Rcpp-devel
mailing list