Deterministically reproducible monte carlo simulation using Rayon?

[masonk]

I'm trying to understand if it's possible to write a monte carlo simulation in Rayon that produces deterministic output. I currently think it's not possible, but I'd like your authoritative answer.

rust-random/rand#399 (comment) (related thread)

I can control the seeds that are used by each Rayon batch using map_init, but:

1. I can't directly control the number of batches that will be spawned (is this guaranteed to == # threads in the current threadpool? Because I think I have seen examples where there are actually many more batches than threads).

2. I can't stop threads from stealing work nondeterministically, and when that happens some items get processed using a different seed.

[adamreichold]

Personally, I have written such simulations by forking the PRNG during initialisation for each entity/agent/individual so that during parallel loops, each iteration would use a separate PRNG state and the result was independent of the order/grouping in which they were processed. Of course, this implies using a reasonably "small" PRNG like PCG64 and choosing a procedure to derive these individual states from a single seed, e.g. using a CSRNG to initialise the PRNG, jumping ahead or utilising PCG64's stream parameter.

There were some subtleties around interactions between entities: They can be synchronized using mutexes/atomics but rather all interactors need to be collected during the first parallel phase and then the interactions themselves must be simulated in an order-independent way in a second parallel pass to achieve determinism. In the end, these complications (and their computational overhead) meant that we dropped work-stealing task-based parallelism and split our simulation along its spatial coordinates with one thread par region instead, but YMMV.

[adamreichold]

(Maybe just to complete the picture: Work-stealing does work though, i.e. load balancing could be rather bad with the spatial work distribution if the number of interactions per region differed significantly due to the simulation setup (heterogeneous landscapes in our case.))

[masonk]

You're right. It's definitely possible to deterministically seed if I am willing to set up the seed on a per-item basis. However, that can be expensive. Even calling rand::thread_local() from the hot function adds 25% to my wall time, and that's pretty cheap. since it's just an atomic check that the ThreadRng is set.

In my particular use-case, since it's a pretty homogeneous workload, I'd happily opt out of workstealing on it, but for many simulations that wouldn't be a good solution.

[adamreichold]

Even calling rand::thread_local() from the hot function adds 25% to my wall time, and that's pretty cheap.

I think the general purpose thread-local PRNG is relatively expensive compared to a cheap small PRNG like PCG64 from rand_pcg with per-item state. And the seeding will happen only during initialisation.

[masonk]

In the Rand issue tracker we hit on a pattern which is mutlithreaded, deterministic, and very fast. I've got a PR out how to show the pattern in the Rand crate.
rust-random/rand#1236