Use Philox PRNG in the Reaction Ensemble

[jngrad]

The Reaction Ensemble (RE) is one of the last consumers of the Mersenne Twister (MT) pseudo-random number generator (PRNG). The MT19937 generator is known to produce correlated data in the first 1 million generated values (Figure 4 in Panneton et al. 2006, DOI:10.1145/1132973.1132974), which we discard in ESPResSo (5ac12d2). The Philox PRNG doesn't suffer from this correlation.

ESPResSo already provides methods to generate a single random value drawn from a uniform or Gaussian distribution:

uint64_t counter = 0;
uint32_t const seed = 42;
int const key_uniform = 1;
int const key_gaussian = 2;
double random_uniform_value  = Random::noise_uniform<RNGSalt::RE,  1>(counter++, seed, key_uniform)
double random_gaussian_value = Random::noise_gaussian<RNGSalt::RE, 1>(counter++, seed, key_gaussian)[0];

The salt RNGSalt::RE needs to be added to the enum class (src/core/random.hpp#L40-L61) beforehand. The key decorrelates sequences with the same seed. The mean and standard deviation of the distribution may differ between MT19937 and Philox.

Once the core has been updated, the samples, tutorial and tests may need adjustments. Many of them rely on a fixed random seed, and with the new PRNG the random sequence will change. Here is a list:

testsuite/python/reaction_ensemble.py
testsuite/python/constant_pH.py
testsuite/python/widom_insertion.py
testsuite/python/wang_landau_reaction_ensemble.py
testsuite/scripts/samples/test_wang_landau_reaction_ensemble.py
testsuite/scripts/samples/test_reaction_ensemble_complex_reaction.py
testsuite/scripts/samples/test_reaction_ensemble.py
testsuite/scripts/samples/test_widom_insertion.py
testsuite/scripts/samples/test_grand_canonical.py
testsuite/scripts/tutorials/test_constant_pH__ideal.py
testsuite/scripts/tutorials/test_constant_pH__interactions.py

The tests in testsuite/python are statistical tests whose seed and parameters were cherry-picked to reduce their runtime without sacrificing their statistical power; if their runtime increases with the new PRNG, it will be necessary to adjust their parameters again.

[jngrad]

There is also another issue with the MT19937 generator: it's not multi-platform. In particular, particle ids are selected at random using std::uniform_int_distribution<int>, whose random sequence depends on the seed and the compiler version. Same problem with the uniformly-distributed and normally-distributed real numbers generators.

The C++ standard says the following about random number generators and their probability distributions in N4713 (link):

29.6.8 Random number distribution class templates
29.6.8.1 In general
³ The algorithms for producing each of the specified distributions are implementation-defined.

ESPResSo compiled with two different compilers, or with two different releases of the same compiler, will not produce the same simulation results for the same reaction ensemble script. Below is a MWE to reproduce the issue; it can be executed online with two releases of Clang to produce different number sequences. This would not happen with Philox, which is multi-platform.

The MD thermostats were changed from MT19937 to Philox a couple of years ago to solve these issues.

MWE

#include <random>
#include <iostream>

namespace ReactionMethods {
class ReactionAlgorithm {
public:
  ReactionAlgorithm(int seed)
      : m_generator(std::mt19937(seed)),
        m_normal_distribution(0.0, 1.0), m_uniform_real_distribution(0.0, 1.0) {
    m_generator.discard(1'000'000);
  }

  int i_random(int maxint) {
    std::uniform_int_distribution<int> uniform_int_dist(0, maxint - 1);
    return uniform_int_dist(m_generator);
  }

  double uniform_random() {
    return m_uniform_real_distribution(m_generator);
  }

  double normal_random() {
    return m_normal_distribution(m_generator);
  }

private:
  std::mt19937 m_generator;
  std::normal_distribution<double> m_normal_distribution;
  std::uniform_real_distribution<double> m_uniform_real_distribution;
};
} // namespace ReactionMethods

int main() {
    int seed = 42;
    auto reaction = ReactionMethods::ReactionAlgorithm(seed);

    std::cout << "random uniform integers: ";
    for (int i = 0; i < 10; i++) {
        std::cout << reaction.i_random(10) << " ";
    }
    std::cout << "\n";

    std::cout << "random uniform reals: ";
    for (int i = 0; i < 4; i++) {
        std::cout << reaction.uniform_random() << " ";
    }
    std::cout << "\n";

    std::cout << "random normal reals: ";
    for (int i = 0; i < 4; i++) {
        std::cout << reaction.normal_random() << " ";
    }
    std::cout << "\n";
    return 0;
}

[jngrad]

In addition, the internal state of the generator and distributions cannot be serialized, which prevents checkpointing of reaction simulations.

[RudolfWeeber]

@jngrad how far is this still needed after your refactoring?

[jngrad]

My refactoring didn't touch the RNG, so that we would have a fair comparison between the C++ and Python implementations. As such, the instructions in the original post still stand. The work should be done on top of the reaction method refactoring PR, since the reaction methods script interface policy changed from LOCAL to GLOBAL, and therefore affects how objects are serialized. One important thing is that SingleReaction objects now have to be serialized using the Python functions __getstate__() and __setstate__(), because GLOBAL policy objects cannot be deserialized recursively.

[jngrad]

• SingleReaction no longer have custom checkpointing code
• a new python test is needed to prove that reactions are decorrelated from one another
  • since a system with Philox instead of MersenneTwister can be reset, once could run a simulation with two concurrent Widom reactions with the same parameters and verify they give the same results with the same seed + salt, and different results with different seed + salt
• it is probably best if we use a different global salt for each reaction method type, and let the user take responsibility for ensuring two concurrent methods of the same type (i.e. same salt) are not used with the same seed
  • e.g. RNGSalt::ReactionEnsemble, RNGSalt::WindomInsertion, etc.

[jngrad]

Offline discussion: we are moving towards a full Python re-implementation of the reaction methods. The RNG will become available at the Python level, and we'll have to choose one that has a serializable state to allow checkpointing. Closing.