Efficiently sample a single row of an array at random I have the following C++ code snippet that randomly samples a single row in an array called pop with num_specs columns and perms rows. In addition, K = 1. The triply-nested for loop uses a pointer for referencing.

Some of the below syntax (such as IntegerVector) is from Rcpp, an R package to integrate C++ code with R code.

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppProgress)]]
#define ARMA_DONT_PRINT_OPENMP_WARNING
#include &lt;RcppArmadillo.h&gt;
#include &lt;RcppArmadilloExtensions/sample.h&gt;
#include &lt;set&gt;
using namespace Rcpp;

int sample_one(int n) {
  return n * unif_rand();
} 

int sample_n_distinct(const IntegerVector&amp; x, 
                      int k,
                      const int * pop_ptr) {

IntegerVector ind_index = RcppArmadillo::sample(x, k, false);
std::set&lt;int&gt; distinct_container;

for (int i = 0; i &lt; k; i++) {
    distinct_container.insert(pop_ptr[ind_index[i]]);
}

return distinct_container.size();
}

// [[Rcpp::export]]
arma::Cube&lt;int&gt; accumulate(const arma::Cube&lt;int&gt;&amp; pop,
                     const IntegerVector&amp; specs,
                     int perms,
                     int K) {

int num_specs = specs.size();
arma::Cube&lt;int&gt; res(perms, num_specs, K);

IntegerVector specs_C = specs - 1;
const int * pop_ptr;
int i, j, k;

for (i = 0; i &lt; K; i++) {
    for (k = 0; k &lt; num_specs; k++) {
        for (j = 0; j &lt; perms; j++) {
            pop_ptr = &amp;(pop(sample_one(perms), 0, sample_one(K)));
            res(j, k, i) = sample_n_distinct(specs_C, k + 1, pop_ptr);
        }
    }
}
return res;
}


While loops in compiled languages aren't bad, it is possible to write slow code. 

I'm not a native C++ programmer, so I don't know all the tricks of the trade. 

Is there a way to reduce the number of levels in the triply-nested for loop above, possibly by employing modular arithmetic in order to see a gain in speed for large input values?

The R code is below:

## Set up container to hold the identity of each individual from each permutation ##
num.specs &lt;- N

## Create an ID for each tag ##
tags &lt;- 1:h

## Assign individuals (N) ##
specs &lt;- 1:num.specs

## Generate permutations. Assume each permutation has N individuals, and sample those 
# individuals' tags from the probabilities ##
  gen.perms &lt;- function() {
      sample(tags, size = num.specs, replace = TRUE, prob = probs)
  }

  pop &lt;- array(dim = c(perms, num.specs, K))

  for (i in 1:K) {
    pop[,, i] &lt;- replicate(perms, gen.perms())
  }

## Perform accumulation ##
HAC.mat &lt;- accumulate(pop, specs, perms, K)

## Example
K &lt;- 1
N &lt;- 100
h &lt;- 5
probs &lt;- rep(1/h, h)
perms &lt;- 100
I can only comment on the C++ code, as I've never used R nor Rcpp.

For an ignorant reader such as myself, it would have helped to omit the using namespace Rcpp;, so that I could see which names come from that library (I did quite a bit of external reading to even understand the code).

Assuming I'm right that sample() returns an array of length k, then we can use range-based for to iterate through its values (since IntegerVector has suitable begin() and end()):

// sample k values from x, without replacement
const Rcpp::IntegerVector ind_index = Rcpp::RcppArmadillo::sample(x, k, false);

// how many different values do they index in pop_ptr?
std::set&lt;int&gt; distinct_container;
for (int i: ind_index) {
    distinct_container.insert(pop_ptr[i]);
}


In accumulate() (which name I dislike, because it sounds like std::accumulate()), we can reduce the scope of i, j, k and pop_ptr quite simply:

for (int i = 0; i &lt; K; ++i) {
    for (int k = 0;  k &lt; num_specs;  ++k) {
        for (int j = 0;  j &lt; perms;  ++j) {
            const int *const pop_ptr = &amp;(pop(sample_one(perms), 0, sample_one(K)));
            res(j, k, i) = sample_n_distinct(specs_C, k + 1, pop_ptr);
        }
    }
}


Also, I'd be inclined to show that the creation of specs_C is more than just a conversion from specs - 1 - that's actually a constructor argument specifying the length of the new vector:

IntegerVector specs_C(specs - 1);


I don't see any code that ever modifies this vector; if it's supposed to remain full of zeros, then it can be declared const.



With those changes, and an automated re-indent, we get something like this for those two functions (I've compiled, but not tested this):

int sample_one(int n)
{
    return n * R::unif_rand();
}

int sample_n_distinct(const Rcpp::IntegerVector&amp; x,
                      int k,
                      const int * pop_ptr)
{
    // sample k values from x, without replacement
    const auto ind_index = Rcpp::RcppArmadillo::sample(x, k, false);

    // how many different values do they index in pop_ptr?
    std::set&lt;int&gt; distinct_container;
    for (int i: ind_index) {
        distinct_container.insert(pop_ptr[i]);
    }

    return distinct_container.size();
}

// [[Rcpp::export]]
arma::Cube&lt;int&gt; accumulate(const arma::Cube&lt;int&gt;&amp; pop,
                           const Rcpp::IntegerVector&amp; specs,
                           int perms,
                           int K)
{
    auto const num_specs = specs.size();
    arma::Cube&lt;int&gt; res(perms, num_specs, K);

    Rcpp::IntegerVector specs_C(specs - 1);

    for (int i = 0;  i &lt; K;  ++i) {
        for (int k = 0;  k &lt; num_specs; ++k) {
            for (int j = 0;  j &lt; perms; ++j) {
                const auto&amp; sampled = pop(sample_one(perms), 0, sample_one(K));
                res(j, k, i) = sample_n_distinct(specs_C, k + 1, &amp;sampled);
            }
        }
    }
    return res;
}




Sadly, I don't have the domain knowledge to suggest more meaningful improvements other than these fairly mechanical changes that improve readability.

In particular, I don't understand why we have to pass specs_C to accumulate() instead of sampling directly from pop_ptr.

I'm also very wary of taking the address of the element returned by pop() - that seems to assume a memory layout that could vary between platforms, for instance.  I'd expect to use one of the subcube view operations to construct the population:

int sample_n_distinct(const Rcpp::IntegerVector&amp; population, int k)
{
    // sample k values from population, without replacement
    const auto sampled = Rcpp::RcppArmadillo::sample(population, k, false);

    // how many different values do have?
    std::set&lt;int&gt; distinct(sampled.begin(), sampled.end());

    return distinct.size();
}

// [[Rcpp::export]]
arma::Cube&lt;int&gt; accumulate(const arma::Cube&lt;int&gt;&amp; pop,
                           const Rcpp::IntegerVector&amp; specs,
                           int perms,
                           int K)
{
    auto const num_specs = specs.size();
    arma::Cube&lt;int&gt; res(perms, num_specs, K);

    for (int i = 0;  i &lt; K;  ++i) {
        for (int k = 0;  k &lt; num_specs; ++k) {
            for (int j = 0;  j &lt; perms; ++j) {
                const auto&amp; sampled = pop.tube(sample_one(perms), sample_one(K));
                const auto sampled_vector
                    = Rcpp::IntegerVector(sampled.begin(), sampled.end());
                res(j, k, i) = sample_n_distinct(sampled_vector, k + 1);
            }
        }
    }
    return res;
}


I may be wrong on exactly what subcube view is required, so pop.tube() may need to be a different call, but you get the idea.  See how much simpler life is when you can use iterators instead of grubbing around with pointers?