Wednesday, February 15, 2012

Wrapping a boost random uniform generator in a class.

@irrati0nal and I decided to use Boost.Random as the random generator library in our pet project.
We wanted to wrap it in a class so we could avoid using it directly in our code. This way it would be easier to change to a different random library if we decided to do it. It was our first experience using the boost libraries. Since we don’t understand templates very well yet, we suffered a bit to get our random generator class working.
We checked the boost documentation example, Bojan Nikolic's examples and some answers in Stack Overflow like this one, but we still didn’t know how to do it. So we decided to hack a solution based on the Stack Overflow answer and Nikolic's examples.
After struggling for a while we finally got to a clumsy solution that compiled and passed the tests:
#ifndef RANDOMNUMBERGENERATOR_H_
#define RANDOMNUMBERGENERATOR_H_

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_01.hpp>

class RandomNumberGenerator {
  public:
    RandomNumberGenerator(int seed) {
     generator = new boost::mt19937(seed);
     distribution = 
          new boost::uniform_01<boost::mt19937>(*generator);
    };

    virtual ~RandomNumberGenerator() {
     delete generator;
     delete distribution;
    };

    double operator() () {
     return (*distribution)();
    }

  private:
    boost::mt19937 *generator;
    boost::uniform_01<boost::mt19937> *distribution;
};
#endif /* RANDOMNUMBERGENERATOR_H_ */
#include <CppUTest/TestHarness.h>
#include "../Math/RandomNumberGenerator.h"
#include <iostream>
#include <iomanip>

using namespace std;

const double PRECISION = 1.0e-8;
static const int seed = 25;
static RandomNumberGenerator gen(seed);

TEST_GROUP(TestRandomNumberGenerator) {
};

TEST(TestRandomNumberGenerator, 
        GetRandomNumbersBetween01) {
  for (int i=0;i<100;++i) {
    double r = gen();
    CHECK( (r < 1.0) and (r >= 0.0) );
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeed) {
  double a[] = {0.8701241323724389,
  0.6960659767501056,
  0.5822769314981997,
  0.4332490300294012,
  0.2788389467168599};

  for (int i=0;i<5;++i) {
   DOUBLES_EQUAL(a[i], gen(), PRECISION);
  }
}
Then we decided to try if we could refine the code to use templates and avoid using dynamic memory, and this was the result:
#ifndef RANDOMNUMBERGENERATOR_H_
#define RANDOMNUMBERGENERATOR_H_

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_01.hpp>

template <int seed>
class RandomNumberGenerator {
  public:
    RandomNumberGenerator() {};
    virtual ~RandomNumberGenerator() {};

    double operator() () {
     static boost::mt19937 generator(seed);
     static boost::uniform_01<boost::mt19937> 
                               distribution(generator);
     return distribution();
    }
};
#endif /* RANDOMNUMBERGENERATOR_H_ */
#include <CppUTest/TestHarness.h>
#include "../Math/RandomNumberGenerator.h"
#include <iostream>
#include <iomanip>

using namespace std;

const double PRECISION = 1.0e-8;
static const int seed = 25;
static RandomNumberGenerator<seed> gen;

TEST_GROUP(TestRandomNumberGenerator) {
};

TEST(TestRandomNumberGenerator, 
        GetRandomNumbersBetween01) {
  for (int i=0;i<100;++i) {
    double r = gen();
    CHECK( (r < 1.0) and (r >= 0.0) );
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeed) {
  double a[] = {0.8701241323724389,
  0.6960659767501056,
  0.5822769314981997,
  0.4332490300294012,
  0.2788389467168599};

  for (int i=0;i<5;++i) {
   DOUBLES_EQUAL(a[i], gen(), PRECISION);
  }
}
We were really excited with our little success so we tried to push the template solution a bit further and got to this:
#ifndef RANDOMNUMBERGENERATOR_H_
#define RANDOMNUMBERGENERATOR_H_

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_01.hpp>

template <int seed, class rng, class distrib>
class RandomNumberGenerator {
  public:
    RandomNumberGenerator() {};
    virtual ~RandomNumberGenerator() {};

    double operator() () {
     static rng generator(seed);
     static distrib distribution(generator);
     return distribution();
    }
};
#endif /* RANDOMNUMBERGENERATOR_H_ */
#include <CppUTest/TestHarness.h>
#include "../Math/RandomNumberGenerator.h"
#include <iostream>
#include <iomanip>

using namespace std;

const double PRECISION = 1.0e-8;
static const int seed = 25;
static RandomNumberGenerator<seed, 
                     boost::mt19937, 
                     boost::uniform_01<boost::mt19937> > gen;

TEST_GROUP(TestRandomNumberGenerator) {
};

TEST(TestRandomNumberGenerator, 
        GetRandomNumbersBetween01) {
  for (int i=0;i<100;++i) {
    double r = gen();
    CHECK( (r < 1.0) and (r >= 0.0) );
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeed) {
  double a[] = {0.8701241323724389,
  0.6960659767501056,
  0.5822769314981997,
  0.4332490300294012,
  0.2788389467168599};

  for (int i=0;i<5;++i) {
   DOUBLES_EQUAL(a[i], gen(), PRECISION);
  }
}
It also worked but we realized that we didn’t understand very well what was going on, and that the declaration of the generator had become much more complicated.
Since we were not going to be using different random generators or statistical distributions, we decided to keep it simple and stick to the second solution which is easier to understand (for us) and also makes the generator creation easier.
Besides having a lot of fun, we’ve realized how little C++ we know and that we are eager to learn more about templates and the boost library.

PS: By the way we are using a Mersenne Twister: "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura, ACM Transactions on Modeling and Computer Simulation: Special Issue on Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30. 
I talked about another version of this generator in a previous post.

Update
There is an update of this post: Revisiting Wrapping a boost random uniform generator in a class.

No comments:

Post a Comment