Revisiting Wrapping a boost random uniform generator in a class.

In a previous post we commented how we had wrapped a boost random generator in a class. This is the code of the class we wrote:

#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_ */

Well, it turned out that this version is wrong if you want to use more than one random generator in your simulation.
We found it out the hard way as we added more tests before adding new features to the class. We didn’t get what we were expecting. After printing out 50 or so random numbers generated using a given seed, we realized that the problem was that the RandomNumberGenerator objects we were using in different tests were not independent.

At work I started a C++ reading club some time ago. This last week we have been reading the Thinking in C++ chapter about name control. There we found out why we were getting a wrong behaviour:

The static objects inside a function (in our code the operator()) are initialized only once, i.e. the first time the function is called. However destructors for static objects are called when main() exits or when the Standard C library function exit() is called.

Therefore, no matter how many new RandomNumberGenerator objects with different seeds we tried to create, they all had the same distribution object than the first one and so the random number sequence, instead of starting again from the begginning for each test as we expected, contained the random numbers that followed the ones we got in the previous test.

Since we want to use more than one random number generator in our simulation, we cannot use the code showed before. We are now using a code that is very similar to the first version showed in a previous post, but avoiding dynamic memory (thanks to another thing we've learned in Thinking in C++ chapter 14).

This is the new code (with some new features) which is passing all 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(seed), distribution(generator) {};

    virtual ~RandomNumberGenerator(){};

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

    double operator() (double lowerLimit, double upperLimit) {
      double width = upperLimit - lowerLimit;
      return width * (*this)() + lowerLimit;
    }

  private:
    boost::mt19937 generator;
    boost::uniform_01<boost::mt19937> distribution;
};
#endif /* RANDOMNUMBERGENERATOR_H_ */

These are the tests:

#include <CppUTest/TestHarness.h>
#include "../Math/RandomNumberGenerator.h"
#include "Constants.h"
#include <iostream>
#include <iomanip>

using namespace std;

static double randomSequence[] = 
 {0.8701241324,0.6960659768,0.5822769315,0.43324903,0.2788389467,
  0.09882423049,0.185911238,0.1808245534,0.4111001296,0.1506180582,
  0.117375552,0.2591784317,0.6849687463,0.6209091851,0.4376110642,
  0.1751907461,0.5562293304,0.2715550633,0.3670803171,0.7612549367,
  0.4023657234,0.9465064413,0.1130407061,0.01872990816,0.447030843,
  0.5319397049,0.5854451165,0.4068064995,0.1619850995,0.6079116813,
  0.5207187904,0.6593447439,0.3260511269,0.3620603553,0.6991862394,
  0.3057702733,0.3663945452,0.9071284649,0.8363745112,0.2565871561,
  0.4813429436,0.2315147296,0.5165022933,0.9432465574,0.3830481281,
  0.4775455245,0.9975408914,0.6317734621,0.5142444882,0.8496764714,
  0.5590532762,0.3113280749,0.03444976453,0.6056593289,0.7199300299,
  0.4580037531,0.4210035447,0.639783914,0.4369351231,0.8001209884,
  0.2817007573,0.1494443719,0.9002743375,0.230675949,0.6696122757,
  0.6084528021,0.4560687505,0.5023682739,0.2898043399,0.7198582094,
  0.5258189649,0.6155499958,0.5592420595,0.8967758908,0.7452838295,
  0.6167110542,0.8283462557,0.1519206658,0.8236944464,0.9035853394,
  0.07714032172,0.4713783977,0.6448620677,0.9572030935,0.3092575835,
  0.5109028067,0.5242537172,0.6404849922,0.9580923351,0.7580416379,
  0.8832009526,0.6661859956,0.2954318929,0.07275792654,0.512375993,
  0.05308085051,0.08870241721,0.9696914584,0.6417167294,0.6661955242};

TEST_GROUP(TestRandomNumberGenerator) {
};

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

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeed) {
  RandomNumberGenerator gen(25);
  for (int i=0; i<100; ++i) {
    DOUBLES_EQUAL(randomSequence[i], gen(), PRECISION);
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeedAndIterval_0_1) {
  RandomNumberGenerator gen(25);
  for (int i=0; i<100; ++i) {
    DOUBLES_EQUAL(randomSequence[i], gen(0, 1), PRECISION);
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeedAndIterval_0_5) {
  RandomNumberGenerator gen(25);
  for (int i=0; i<100; ++i) {
    DOUBLES_EQUAL(randomSequence[i] * 5., gen(0, 5), PRECISION);
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeedAndIterval_1_5) {
  RandomNumberGenerator gen(25);
  for (int i=0; i<100; ++i) {
    DOUBLES_EQUAL(randomSequence[i] * 4. + 1., gen(1, 5), PRECISION);
  }
}

TEST(TestRandomNumberGenerator, 
        GetRandomNumberSequenceWithGivenSeedAndIterval_1_2) {
  RandomNumberGenerator gen(25);
  for (int i=0; i<100; ++i) {
    DOUBLES_EQUAL(randomSequence[i] * 1. + 1., gen(1, 2), PRECISION);
  }
}

As we said in the previous post about this topic, we still have so much to learn...