Hi Dirk
On Sat, 1 Jun 2013, Dirk Eddelbuettel wrote:
Hi Kouros,
On 1 June 2013 at 16:13, [email protected] wrote:
| The problem is that I cannot (or do not know how to) use the R random seeding
| mechanism to make the permutation resampling process
| reproducible. I would appreciate any advice on this matter (e.g., replacing
| std::random_shuffle with a suitable Rcpp function or to be pointed to
| existing solutions ).
No direct solution for you but
a) you can probably seed the standard library RNGs as well (which may just
need a helper function)
b) you can look at the very nice sample() function contributed to
RcppArmadillo and see if you can cook something similar up for
RcppEigen.
There are posts on the Rcpp Gallery about b).
Dirk
Thanks for the prompt and helpful response. I adopted the solution you had
outlined in http://gallery.rcpp.org/articles/stl-random-shuffle/
to resolve this issue along the lines of a) as shown below.
I agree that the sample() function contributed to RcppArmadillo is
very nice.
Thanks again for your help.
Take care, Kouros
# Adopted from
http://stackoverflow.com/questions/15858569/randomly-permute-rows-columns-of-a-matrix-with-eigen
# This function returns a random column and a row random permutation replicate of a
# given matrix
permmatrix<-'
List permmatrix(NumericMatrix Xr){
RNGScope scope;
const Eigen::Map<Eigen::MatrixXd> X(as<Eigen::Map<Eigen::MatrixXd> >(Xr));
int nr = X.rows();
int nc = X.cols();
// Permute Columns
Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic> permc(nc);
permc.setIdentity();
std::random_shuffle(permc.indices().data(),permc.indices().data()+permc.indices().size(),randWrapper);
Eigen::MatrixXd Xpc = X * permc;
// Permute Rows
Eigen::PermutationMatrix<Eigen::Dynamic,Eigen::Dynamic> permr(nr);
permr.setIdentity();
std::random_shuffle(permr.indices().data(),
permr.indices().data()+permr.indices().size(),randWrapper);
Eigen::MatrixXd Xpr = permr * X;
return(List::create(Named("X")=X,Named("Xpc")=Xpc,Named("Xpr")=Xpr));
}'
# Adopted from http://gallery.rcpp.org/articles/stl-random-shuffle/
randWrapper<-'inline int randWrapper(const int n) { return
floor(unif_rand()*n); }'
library(RcppEigen)
permMatrix<-cppFunction(permmatrix,depends="RcppEigen",include=randWrapper)
permMatrix(matrix(1:24,4,6))
$X
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 5 9 13 17 21
[2,] 2 6 10 14 18 22
[3,] 3 7 11 15 19 23
[4,] 4 8 12 16 20 24
$Xpc
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 5 17 13 9 1 21
[2,] 6 18 14 10 2 22
[3,] 7 19 15 11 3 23
[4,] 8 20 16 12 4 24
$Xpr
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3 7 11 15 19 23
[2,] 2 6 10 14 18 22
[3,] 1 5 9 13 17 21
[4,] 4 8 12 16 20 24
set.seed(12123)
permMatrix(matrix(1:24,4,6))$Xpc
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 17 21 13 9 5
[2,] 2 18 22 14 10 6
[3,] 3 19 23 15 11 7
[4,] 4 20 24 16 12 8
set.seed(12123)
permMatrix(matrix(1:24,4,6))$Xpc
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 17 21 13 9 5
[2,] 2 18 22 14 10 6
[3,] 3 19 23 15 11 7
[4,] 4 20 24 16 12 8
print(sessionInfo(),locale=FALSE)
R version 3.0.1 (2013-05-16)
Platform: x86_64-pc-linux-gnu (64-bit)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] RcppEigen_0.3.1.2.1 Matrix_1.0-12 lattice_0.20-15
[4] Rcpp_0.10.3
loaded via a namespace (and not attached):
[1] grid_3.0.1 tools_3.0.1
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
