[Rcpp-devel] Random matrix row and column permutations in RcppEigen

[Kouros Owzar]

Hi Dirk

On Sat, 1 Jun 2013, Dirk Eddelbuettel wrote:

>
> Hi Kouros,
>
> On 1 June 2013 at 16:13, owzar001 at duke.edu 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
>