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);
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
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.
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
>
--
Yi-Hsin Erica Tsai, Ph.D.
NSF Postdoctoral Fellow
Department of Biological Sciences (Carstens lab)
Louisiana State University
http://www.duke.edu/~yet2/
[email protected]
# 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
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
