Vfrdtyky

The input of my function consists of two matrices mat1 and mat2 and the number of permutations B. Both mat1 and mat2 have m columns, but different number of rows.

The function first permutes the rows of both matrices (while maintaining column information). It then performs some operation that compares columns of permuted versions of mat1 and mat2.

The following is an example of my function permute_data. The comparison function CompareMatCols() outputs a vector of length m.

Question
What is the best way to initialize my output list object? I've seen several posts indicating the limitations of push_back. Both B and m will be on the order of ~10000, so an efficient way would be ideal.

#include <Rcpp.h>

#include <math.h>

//#include <random> //for std::shuffle



using namespace std;

using namespace Rcpp;



// [[Rcpp::export]]

NumericVector ColMax(NumericMatrix X) {

  NumericVector out = no_init(X.cols());

  for(int j = 0; j < X.cols(); ++j) {

      double omax = X(0,j);

    for(int i = 0; i < X.rows(); ++i){

      omax = std::max(X(i,j),omax);

    }

    out[j] = omax;

  }

  return out;

}



// [[Rcpp::export]]

NumericVector vecmin(NumericVector vec1, NumericVector vec2) {

  int n = vec1.size();

  if(n != vec2.size()) return 0;

  else {

    NumericVector out = no_init(n);

    for(int i = 0; i < n; i++) {

      out[i] = std::min(vec1[i], vec2[i]);

    }

    return out;

  }

}



// [[Rcpp::export]]

List permute_data(NumericMatrix mat1,NumericMatrix mat2,int B) {



  List out(B); // How to initialize this???, Will be large ~10000 elements

  int N1 = mat1.rows();

  int N2 = mat2.rows();

  int m = mat1.cols(); //Will be large ~10000 elements



  // Row labels to be permuted

  IntegerVector permindx = seq(0,N1+N2-1);

  NumericMatrix M1 = no_init_matrix(N1,m);

  NumericMatrix M2 = no_init_matrix(N2,m);



  for(int b = 0; b<B; ++b){

    // Permute the N1+N2 rows

   /*std::random_device rng;

   std::mt19937 urng(rng()); //uniform rng

   std::shuffle(permindx.begin(),permindx.end(),urng);*/

    permindx = sample(permindx,N1+N2); //Use Rcpp's function to work with R's RNG

    for(int j=0; j<m; ++j){

        // Pick out first N1 elements of permindx

        for(int i=0; i<N1; ++i){

          if(permindx[i]>=N1){ //Check that shuffled index is in bounds

            M1(i,j) = mat2(permindx[i],j);

          } else{

            M1(i,j) = mat1(permindx[i],j);

          }

        }

        // Pick out last N2 elements of permindx

        for(int k=0; k<N2; ++k){

          if(permindx[k+N1]<N1){ //Check that shuffled index is in bounds

            M2(k,j) = mat1(permindx[k+N1],j);

          } else{

            M2(k,j) = mat2(permindx[k+N1],j);

          }

        }

    }

    out[b] = vecmin(ColMax(M1),ColMax(M2)); //a vector of length m

  }

  return(out);

}



/***R

set.seed(1)

X = matrix(rnorm(3*5),ncol=5)

Y = matrix(rnorm(5*5),ncol=5)

B = 5

res = permute_data(X,Y,B)

*/

Edit: Added lines to address @duckmayr's point.
Edit 2: Per Dirk's suggestion, I included a minimally complete verifiable example.

In the end, the solution I went with is the one seen above:

// [[Rcpp::export]]

List permute_data(NumericMatrix mat1,NumericMatrix mat2,int B) {



  List out(B); // Will be large ~5000 elements

  int N1 = mat1.rows();

  int N2 = mat2.rows();

  int m = mat1.cols(); //Will be large ~10000 elements



  // Row labels to be permuted

  IntegerVector permindx = seq(0,N1+N2-1);

  NumericMatrix M1 = no_init_matrix(N1,m);

  NumericMatrix M2 = no_init_matrix(N2,m);



  for(int b = 0; b<B; ++b){

    // Permute the N1+N2 rows

    permindx = sample(permindx,N1+N2); //Use Rcpp's function to work with R's RNG

    for(int j=0; j<m; ++j){

        // Pick out first N1 elements of permindx

        for(int i=0; i<N1; ++i){

          if(permindx[i]>=N1){ //Check that shuffled index is in bounds

            M1(i,j) = mat2(permindx[i],j);

          } else{

            M1(i,j) = mat1(permindx[i],j);

          }

        }

        // Pick out last N2 elements of permindx

        for(int k=0; k<N2; ++k){

          if(permindx[k+N1]<N1){ //Check that shuffled index is in bounds

            M2(k,j) = mat1(permindx[k+N1],j);

          } else{

            M2(k,j) = mat2(permindx[k+N1],j);

          }

        }

    }

    out[b] = vecmin(ColMax(M1),ColMax(M2)); //a vector of length m

  }

  return(out);

}