stan-dev / stan

Stan development repository. The master branch contains the current release. The develop branch contains the latest stable development. See the Developer Process Wiki for details.
https://mc-stan.org
BSD 3-Clause "New" or "Revised" License
2.6k stars 370 forks source link

Add 1 to random seed #3286

Closed SteveBronder closed 6 months ago

SteveBronder commented 6 months ago

Submission Checklist

Summary

Fixes #3285 by just adding 1 to the given seed

How to Verify

python3 ./runTests.py ./src/test/unit/services/util/create_rng_test.cpp

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company): Simon's Foundation

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

WardBrian commented 6 months ago

I've read the boost docs and couple papers on MIXMAX. The actual RNG state is 128 bits, and it does require at least 1 nonzero bit in the seed to function.

Luckily for us, the boost seed taken in is a uint64_t, and our seeds are only uint32s, so we can create the seed like

uint64_t seed_boost = static_cast<uint64_t>(seed) + static_cast<uint64_t>(chain) + 1;

and guarantee that there is always a nonzero bit, without the issue of wrapping around leading to another zero appearing at the top of the range.

This will 1) require updating all the tests which use fixed seeds and 2) be rather difficult if not impossible to test on its own, but I think both of those are fine

nhuurre commented 6 months ago

Luckily for us, the boost seed taken in is a uint64_t, and our seeds are only uint32s, so we can create the seed like

boost::random::mixmax also has a constructor that takes four uint32s. Why not use that?

inline rng_t create_rng(unsigned int seed, unsigned int chain) {
  rng_t rng(1, 1, seed, chain);
  return rng;
}

And if you're going to create a 64-bit seed, I think concatenating the 32-bit parts is better than summing

uint64_t seed_boost = (static_cast<uint64_t>(seed) << 31) + static_cast<uint64_t>(chain) + 1;

EDIT: with one bit of overlap to avoid wrap-around on +1

WardBrian commented 6 months ago

@nhuurre I agree after reading more.

A feature of mixmax which was described by the boost authors is that seed and seed+1 are independent in the case of calling the 4-argument version with (0,0, seed<<32, seed).

If that remains true for values other than 0 for the first two arguments (e.g. (0,1, seed<<32, seed)), then I think we should just call the 4 argument constructor with (0, 1, user_seed, chain_id). I'm inquiring now. If there is some magic in the first two arguments being zero, we can do the shift-by-31 and add 1 you describe.

SteveBronder commented 6 months ago

Going to close this as it sounds like there are better answers