Re: [boost] [random] new threefry random engine

Subject: Re: [boost] [random] new threefry random engine
From: Thijs van den Berg (thijs_at_[hidden])
Date: 2014-05-06 15:51:12

> 
> I think you'll like the new example code better.  From the new README:
> 
>    typedef threefry<4, uint32_t> Prf;
>    counter_based_engine<Prf, 32> cbeng(seed);
>    for(size_t i=0; i<Natoms; ++i){
>        float boltzmannfactor = sqrt(kT/atoms[i].mass);
>        normal_distribution nd(0., rmsvelocity);
>        Prf::domain_type base = {atoms[i].id, timestep, THERMALIZE_CTXT};
>        cbeng.restart(base);
>        atoms[i].vx = boltzmannfactor*nd(cbeng);
>        atoms[i].vy = boltzmannfactor*nd(cbeng);
>        atoms[i].vz = boltzmannfactor*nd(cbeng);
>    }
> 
> You can seed the counter_based_engine directly with an arithmetic seed, or,
> if you prefer,
> you can use the the much larger key_type, or you can pre-initialize a Prf
> and use that.

I tried to examine what happend when you draw more samples inside the loop, like this

        for (size_t j=0; j<100; ++j) {
            atoms[i].vx += mbd(cbeng);
            atoms[i].vy += mbd(cbeng);
            atoms[i].vz += mbd(cbeng);
        }

and that caused an abort trap. I expect that the cbeng ran out of samples and that is not the behaviour  you want. To be compatible you’ll need to allow the distribution mdb() can draw any amount of samples. Mbd() can be a probability distribution that consumes many random samples from the engine. E.g. I’ve just finished a distribution sampler that consumes 4 per draw. The same issues arises if you have many more variables to set than vx,vy,vz, which make the usage error-prone. One solution for this particular example is that you move the domain_type variable inside the 2nd loop and make it a function of both i,j from the loops. However that still isn’t a nice design because it can be very inefficient. Every time time you call restart() it will generate 4x64 bits of random data on it’s first call, ..but you only consume 3x32 bits. Putting the burden of optimising this on the user is not a good idea. 

The version I’ve finished here https://github.com/sitmo/threefry is a standard random engine concept and doesn’t need a domain_type, key_type,  restart() and doesn’t have a  low limit on the number of samples that users need to be aware of. Its intuitive to use like any other random engine -as can be seen in the next example-. I *think* it also uses less memory and is also faster because I’ve optimised the engine for it’s purpose. 

    for(size_t i=0; i<Natoms; ++i){
        float rmsvelocity = sqrt(kT/atoms[i].mass);
        normal_distribution<float> mbd(0., rmsvelocity);
        threefry4x64_13 eng(i);
        atoms[i].vx += mbd(eng);
        atoms[i].vy += mbd(eng);
        atoms[i].vz += mbd(eng);
    }