mercredi 25 janvier 2017

generating random numbers with mt19937_64 for multithread monte carlo simulation

I am running a monte carlo simulation on multiple threads. I have a Draw class which handles random number generation. It uses mt19937-64 as the generator of std::uniform_real_distribution.

template<typename T, typename R>
class Draw {
    T _dist;
    typedef decltype(_dist.min()) result_type;

    public:
    // Draw : T R -> Draw
    //! GIVEN
    //!   1. dist - A distribution.
    //!   2. gen - A random number generator.
    Draw(T dist, R gen) : _dist(dist), _gen(gen) {}

    virtual ~Draw() = default;

    // () : -> Value
    //! GIVEN
    //! RETURNS
    //! value drawn from the distribution, which can be any result type
    //!   supported by the distribution.
    result_type operator()() const { return _draw(); }

    // seed : NonNegInt -> Void
    //! GIVEN
    //!   1. seed - A random number generator (RNG) seed.
    //! EFFECT
    //!   Seeds the RNG with the given seed.
    void seed(unsigned long seed) { _gen.seed(seed);}

    private:
    R _gen;

    // draw : -> Value
    // GIVEN:
    // RETURNS: A value drawn from the distribution, which can be any result
    //          type supported by the distribution.
    std::function<result_type()> _draw = bind(_dist,_gen);
};

// standard uniform distribution ~ Unif(a=0, b=1)
class DrawUnif :
public Draw<std::uniform_real_distribution<double>,std::mt19937_64>
{
    typedef std::mt19937_64 genny;
    typedef std::chrono::system_clock clk;
    typedef std::uniform_real_distribution<double> dist;

    public:
    DrawUnif() :
      Draw(dist{0,1}, genny(clk::now().time_since_epoch().count())) {}

    //! GIVEN
    //! -----
    //!   1. seed - A seed for the random number generator.
    DrawUnif(unsigned long seed) : Draw(dist{0,1}, genny(seed)) {}

     virtual ~DrawUnif() = default;
};

Each thread has access to the following shared pointer

  typedef std::shared_ptr<DrawUnif> DrawUnifShrPtr;
  DrawUnifShrPtr    _unif;

Which is initialized by

_unif(DrawUnifShrPtr(new DrawUnif { seed }));

Each thread has several functions which frequently call draw(*_unif) to generate a random number. The results seem correct but I am wondering if two threads call

draw(*_unif) 

at the same time, what would happen?

Then I allocated a new shared_pointer with different seed to each cohort:

 _unif(std::make_shared<DrawUnif>(*c._unif));
 _unif->seed(a_new_seed);

Now the results look wiered! Each thread is getting exact same random!! numbers and provide exact same result. To summarize:

1- what would happen if multiple treads call the draw at the same time?

2- Why different seeds get exact same results.

Thank you




Aucun commentaire:

Enregistrer un commentaire