Efficiently sample a single row of an array at random

Question

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 <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
#include <set>
using namespace Rcpp;

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

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

IntegerVector ind_index = RcppArmadillo::sample(x, k, false);
std::set<int> distinct_container;

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

return distinct_container.size();
}

// [[Rcpp::export]]
arma::Cube<int> accumulate(const arma::Cube<int>& pop,
                     const IntegerVector& specs,
                     int perms,
                     int K) {

int num_specs = specs.size();
arma::Cube<int> res(perms, num_specs, K);

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

for (i = 0; i < K; i++) {
    for (k = 0; k < num_specs; k++) {
        for (j = 0; j < perms; j++) {
            pop_ptr = &(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 <- N

## Create an ID for each tag ##
tags <- 1:h

## Assign individuals (N) ##
specs <- 1:num.specs

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

  pop <- array(dim = c(perms, num.specs, K))

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

## Perform accumulation ##
HAC.mat <- accumulate(pop, specs, perms, K)

## Example
K <- 1
N <- 100
h <- 5
probs <- rep(1/h, h)
perms <- 100

Answer 1

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<int> 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 < K; ++i) {
    for (int k = 0;  k < num_specs;  ++k) {
        for (int j = 0;  j < perms;  ++j) {
            const int *const pop_ptr = &(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& 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<int> distinct_container;
    for (int i: ind_index) {
        distinct_container.insert(pop_ptr[i]);
    }

    return distinct_container.size();
}

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

    Rcpp::IntegerVector specs_C(specs - 1);

    for (int i = 0;  i < K;  ++i) {
        for (int k = 0;  k < num_specs; ++k) {
            for (int j = 0;  j < perms; ++j) {
                const auto& sampled = pop(sample_one(perms), 0, sample_one(K));
                res(j, k, i) = sample_n_distinct(specs_C, k + 1, &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& 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<int> distinct(sampled.begin(), sampled.end());

    return distinct.size();
}

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

    for (int i = 0;  i < K;  ++i) {
        for (int k = 0;  k < num_specs; ++k) {
            for (int j = 0;  j < perms; ++j) {
                const auto& 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?