mardi 29 janvier 2019

using gsl random number generator in armadillo

I'm writing some code in C++ for scientific computing purposes. The linear algebra library I chose to use is Armadillo. My problem is that Armadillo uses either the simple rand() function or a Mersenne Twister rng, while I need better quality numbers, for example ranlux. GSL provides an implementation of ranulx, and Armadillo allows for an external rng by defining the preprocessor macro ARMA_RNG_ALT and writing an appropriate header file that we will call arma_rng_alt.hpp.

Now, I'm no C++ wizard but I managed to do it by more or less copying and pasting from the existing headers arma_rng_cxx98.hpp and arma_rng_cxx11.hpp. My solution might not be the prettiest but I think it works alright. This is an extract of my header:

#if defined(ARMA_RNG_ALT)
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#define RNG_TYPE gsl_rng_ranlxd1

class arma_rng_alt
{
    public:

    typedef unsigned long int seed_type;

    inline void set_seed(const seed_type val);

    arma_inline int    randi_val();
    arma_inline double randu_val();
    arma_inline double randn_val();


    private:

    gsl_rng* engine;
};

inline
void
arma_rng_alt::set_seed(const arma_rng_alt::seed_type val)
{
    engine = gsl_rng_alloc(RNG_TYPE);
    gsl_rng_set(engine, val);
}

arma_inline
int
arma_rng_alt::randi_val()
{
    return gsl_rng_get(engine);
}

arma_inline
double
arma_rng_alt::randu_val()
{
    return gsl_rng_uniform(engine);
}

arma_inline
double
arma_rng_alt::randn_val()
{
    return gsl_ran_gaussian_ziggurat(engine, 1.);
}

#endif

Now, this would be ok if this external rng was treated by armadillo as a class without a private variable (like in the arma_rng_cxx98.hpp header, that uses rand() and doesn't need to allocate an engine). Instead I need to allocate the engine and call it every time, therefore I need to instantiate the class otherwise I get a "cannot call member function without object" error. So I actually needed to modify a header from the Armadillo library itself, namely arma_rng.hpp. To make things clearer, this is how a typical chunk of the original arma_rng.hpp to generate random integers looks like:

template<typename eT>
struct arma_rng::randi
  {
  arma_inline
  operator eT ()
    {
    #if   defined(ARMA_RNG_ALT)
      {
      return eT( arma_rng_alt::randi_val() );
      }
    #elif defined(ARMA_USE_EXTERN_CXX11_RNG)
      {
      return eT( arma_rng_cxx11_instance.randi_val() );
      }
    #else
      {
      return eT( arma_rng_cxx98::randi_val() );
      }
    #endif
    }

As you can see, if cxx98 or alt is used, the function randi_val is simply called, while if cxx11 is used, randi_val is called on an actual instance of the calss arma_rng_cxx11 (a thread local variable declared earlier in the header). What I need is something similar to cxx11, so this is how I modified the library header arma_rng.hpp for this specific example:

template<typename eT>
struct arma_rng::randi
  {
  arma_inline
  operator eT ()
    {
    #if   defined(ARMA_RNG_ALT)
      {
      return eT( arma_rng_alt_instance.randi_val() );
      }
    #elif defined(ARMA_USE_EXTERN_CXX11_RNG)
      {
      return eT( arma_rng_cxx11_instance.randi_val() );
      }
    #else
      {
      return eT( arma_rng_cxx98::randi_val() );
      }
    #endif
    }

after declaring my own thread local instance of the class arma_rng_alt I've created.

Works like a charm, the catch is that now I need to run my code on a HPC cluster and, while I can create my own arma_rng_alt.hpp, I can't modify the library header arma_rng.hpp itself. So I need a way to make the original header work (so without instantiating the class). Do you have any idea on how to do that? My C++ expertise is frankly stretched to the limit here.

Thanks




Aucun commentaire:

Enregistrer un commentaire