The handling of random numbers

[quaquel]

Currently, the seeded random number generator resides in the model class. However, this means that all other classes that need random number generation need a reference to the model class. In particular, for AgentSet, and if we add some CellCollection, this becomes quite burdensome because any class that generates them also needs a reference to the model.

One possible solution would be to turn the seeded random number generator into a singleton. Then, in the __init__, or wherever you need it, you instantiate this class and get the singleton instance back.

It is even possible to expand upon this. A common technique to reduce variance in stochastic models is to use common random numbers. At the moment this is not easy to do in MESA unless you want to write everything yourself. With a Random class that tracks its instances based on the seed, this becomes easier to do.

Something like the following API would work well for this:

from mesa import Random


random = Random() # returns the generic random number generator as per the current model.random

my_random = Random(42) # returns a unique new random number generator seeded with 42

# somewhere else in the code
# these use cases need a bit more thought
another_random = Random(42, new=False) # returns my_random
yet_more_random = Random(42, new=True) # returns a new seeded generator

Implementation is pretty straightforward and can even be done in a backward-compatible manner.

[rht]

If you reuse the seed within a simulation run, then the random values are repeated: another_random and yet_more_random will start with the same exact same numbers. This turns the numbers generated to be less random.

Another problem, even when you have fixed the first problem, is that mapping a random state to a single seed number doesn't apply to NumPy's default RNG, PCG64.

Tangent point: I looked at Julia's RNG just now. They had switched from MersenneTwister to Xoshiro256++. I'm not sure how it compares with PCG64. A HN discussion about PCG and alternatives.

[rht]

I might have misread your code.

another_random = Random(42, new=False) # returns my_random
yet_more_random = Random(42, new=True) # returns a new seeded generator

Do next(another_random) and next(another_random) give the same number?

Another question: when you do from mesa import Random, initializes an RNG object, and reuse the same object, is this a global RNG that is not tied to any particular model object? If so, how do you ensure thread safety when doing computation in parallel within a single simulation run?

[quaquel]

I am trying to flash out a possible API here, so this is very much WIP. But a few thoughts.

1. No idea about within-model parallelization. Is that even in scope for MESA? If you could parallelize using threats (so not Python). The same seed object would be returned. It's a singleton, after all. You still get race conditions and other nasty stuff, which would cause a model to become non-deterministic. This is why most discrete simulation packages I have worked with avoid thread-based parallelization of single runs. Rather, you parallelize the replications, which is much easier (also in python).
2. On the code, if new=True, you would return a new seeded RNG starting afresh, while new=False would return the already seeded RNG and whatever state it is at. So, assuming this second RNG has already been called several times, you would get different results. Hope this is clear.

3. If you reuse the seed within a simulation run, then the random values are repeated: another_random and yet_more_random will start with the same exact same numbers. This turns the numbers generated to be less random.

Which is exactly the point of using common random numbers. For example, you might use common random numbers to ensure that your agents' initial location and state are always the same. And then use a set of different seeds only to explore the dynamic behavior. Or, in discrete event simulation, you might use it to ensure that the arrival process is always the same.
4. Are the details of RNG generation within the scope for MESA, or are we just using random and numpy.random.Generator? I would simply not worry about it.

[rht]

Rather, you parallelize the replications, which is much easier (also in python).

This makes sense and what I have done so far as well. Also, for larger scale model, it would be better to use Agents.jl, etc. As such, my thread-safety concern doesn't need be considered in the design anymore.

Which is exactly the point of using common random numbers. For example, you might use common random numbers to ensure that your agents' initial location and state are always the same. And then use a set of different seeds only to explore the dynamic behavior. Or, in discrete event simulation, you might use it to ensure that the arrival process is always the same.

Did you mean that the simulation ought to be deterministic (always the same initial condition and outcome, for a given seed)? I don't see how this is related to variance reduction. It's about reproducibility instead.

Are the details of RNG generation within the scope for MESA, or are we just using random and numpy.random.Generator?

It's definitely possible that some users would use their own RNG. That said, even if we are just using random and numpy.random.Generator, the latter's can't be stored as a single integer like CPython's random does.

[rht]

I see this statement to be the main motivation of this discussion (with the variance reduction being a consequence of not accessing the RNG via the model instance):

However, this means that all other classes that need random number generation need a reference to the model class. In particular, for AgentSet, and if we add some CellCollection, this becomes quite burdensome because any class that generates them also needs a reference to the model.

From the agent's perspective, it needs to access the model to get both its RNG and clock anyway. There could be other model attributes beyond the RNG and the clock that the space/cell may need to access to. This is regardless of whether there is a solution to access the RNG without via the model instance.

[quaquel]

The critical point is the need for AgentSet and CellCollection (see https://github.com/projectmesa/mesa/tree/experimental/cell-space) to access the model's RNG. Using a Singleton is a simple and convenient solution for this problem. It is also much easier than passing a reference to the model to them, because this requires whatever code is instantiating them needs to have it as well.

But, with just a little bit of thought it is easy also to use this to support common random numbers, which for people that are familiar with them is a nice benefit.

Likewise, with a little bit of care, it becomes trivial to have both numpy RNG and python RNG covered by this as well.