c++rcpparrayfire

Alternating Error: "Invalid dimension for argument 0"

In converting the example below to a gfor loop. I encountered an error of the type "Invalid dimension for argument 0", the full error message below. However, the error occurs, then the function runs, then the same error. This pattern repeats. I am confused and am wondering if this error is in someway system dependent.

Full error message:

 Error in random_shuffle(theta, 5, 1) : 
 ArrayFire Exception (Invalid input size:203):
In function af_err af_assign_seq(af_array *, const af_array, const unsigned int, const af_seq *, const af_array)
In file src/api/c/assign.cpp:168
Invalid dimension for argument 0
Expected: (outDims.ndims() >= inDims.ndims())

A second problem, the seed fails to change with the input parameter, when using the gfor loop.

#include "RcppArrayFire.h"

using namespace Rcpp;
using namespace RcppArrayFire;

// [[Rcpp::export]]
af::array random_shuffle(const RcppArrayFire::typed_array<f64> theta, int counts, int seed){

  const int theta_size = theta.dims()[0];
  af::array out(counts, theta_size, f64);
  af::array seed_seq = af::seq(seed, seed+counts);

// for(int f = 0; f < counts; f++){
gfor ( af::seq f, counts-1 ){

  af::randomEngine engine;
  engine.setSeed(af::sum<double>(seed_seq(f)));

  af::array index_shuffle(1, u16);

  af::array temp_rand(1, f64);
  af::array temp_end(1, f64);

  af::array shuffled = theta;

// implementation of the Knuth-Fisher-Yates shuffle algo
  for(int i = theta_size-1; i > 1; i --){

    index_shuffle = af::round(af::randu(1, u16, engine)/(65536/(i+1)));
    temp_rand = shuffled(index_shuffle);
    temp_end = shuffled(i);
    shuffled(index_shuffle) = temp_end;
    shuffled(i) = temp_rand;
    }

  out(f, af::span) = shuffled;
  }
  return out;
}


/*** R
theta <- 10:20
random_shuffle(theta, 5, 1)
random_shuffle(theta, 5, 2)
*/

Updated with Ralf Stunber's solution, but 'shuffled' samples in Column space.

    // [[Rcpp::export]]
af::array random_shuffle2(const RcppArrayFire::typed_array<f64> theta, int counts, int seed) {
  int len = theta.dims(0);
  af::setSeed(seed);
  af::array tmp = af::randu(len, counts, 1);
  af::array val, idx;
  af::sort(val, idx, tmp, 1);
  af::array shuffled = theta(idx);
  return af::moddims(shuffled, len, counts);
}
/*** R
random_shuffle2(theta, 5, 1)
*/

Here is a picture of output, sampling with replacement:enter image description here

In the second part, of 50 repetitions, the samples moves towards anergodic outcome.

Solution

  • Why do you want to use multiple RNG engines in parallel? There is really no need for this. In general, it should be sufficient to use only the global RNG engine. It should also be sufficient to set the seed of this engine only once. You can do this from R with RcppArrayFire::arrayfire_set_seed. Besides, random number generation within a gfor loop does not work as one might expect, c.f. http://arrayfire.org/docs/page_gfor.htm.

    Anyway, I am not an expert in writing efficient GPU algorithms, which is why I like using the methods implemented in libraries like ArrayFire. Unfortunately ArrayFire does not have a shuffle algorithm, but the corresponding issue has a nice implementation, which can be generalized to your case with multiple shuffles:

    // [[Rcpp::depends(RcppArrayFire)]]
#include "RcppArrayFire.h"

// [[Rcpp::export]]
af::array random_shuffle(const RcppArrayFire::typed_array<f64> theta, int counts, int seed) {
    int len = theta.dims(0);
    af::setSeed(seed);
    af::array tmp = af::randu(counts, len, 1);
    af::array val, idx;
    af::sort(val, idx, tmp, 1);
    af::array shuffled = theta(idx);
    return af::moddims(shuffled, counts, len);
}

    BTW, depending on the later usage it might make more sense to arrange the different samples in columns instead of rows, since both R and AF use column major layout.