NumPy - Advancing Random BitGenerator

Viewed 191

I'm trying to work on a script that, regardless of the count of execution threads/processes I'm using, I'll be able to generate the same dataset of random data deterministically. As the dataset is for consumption by a relational database by means of a csv file(s), row order is not important.

The logic I have been working on is to split the work among n workers, and advance the BitGenerator of the Generator object accordingly. For example, if I want to generate 1000 rows, below will output the equivalent dataset:

  • 1 process generating all 1000 rows.
  • 2 processes (p1 and p2):
    • p1 generating rows from 1-500;
    • p2 generating rows from 501-1000, with the BitGenerator advanced by 500.

Python standard random.Random() object does not allow to select the BitGenerator, and the default is MT19937 which can't be advanced anyway.

I have been thus looking at NumPy objects Generator and BitGenerators and they are very helpful, however, I run into below issue. First however, this is my env:

import sys
import numpy as np

print(sys.version)
print(np.__version__)
3.9.9 (v3.9.9:ccb0e6a345, Nov 15 2021, 13:29:20) 
[Clang 6.0 (clang-600.0.57)]
1.22.0

I have been doing these tests with PCG64, PCG64DXSM and Philox, but the result is similar. What I see is that I can't correctly predict by how much I have to advance the BitGenerator to make sure I get the expected result. See example below:

bitgen = np.random.PCG64(seed=0)
rng = np.random.Generator(bitgen)

for i in range(6):
    print(i, rng.random())

print("-------------------------")

# repeat, advance the bitgen by 1
bitgen = np.random.PCG64(seed=0).advance(1)
rng = np.random.Generator(bitgen)

for i in range(5):
    print(i, rng.random())

See how rows 1-5 of the first group match with 0-4 of the second group. They correctly advanced by 1.

0 0.6369616873214543

1 0.2697867137638703
2 0.04097352393619469
3 0.016527635528529094
4 0.8132702392002724
5 0.9127555772777217
-------------------------
0 0.2697867137638703
1 0.04097352393619469
2 0.016527635528529094
3 0.8132702392002724
4 0.9127555772777217

However, other random functions don't behave the same. Below the output of the above script substituting rng.random() with rng.integers(1000) and rng.bytes(1).

Notice in both cases how the values start matching after 2 rows, not one.

Output using integers(1000):

0 850
1 636

2 511
3 269
4 307
5 40
-------------------------
0 511
1 269
2 307
3 40
4 75

Output using bytes(1):

0 b'_'
1 b'\xcf'

2 b'!'
3 b'\xf8'
4 b'\xb8'
5 b'\xa9'
-------------------------
0 b'!'
1 b'\xf8'
2 b'\xb8'
3 b'\xa9'
4 b'\x1d'

Tricky outputs are also seen when the output is an array. Function rng.choice(5, 3, p=[0.1, 0, 0.3, 0.6, 0]) correctly shifts the number by 1, but the array result per se is not what I wanted: in below, the number start matching after 1 value, but the arrays are all wrong!

0 [3 2 0]
1 [0 3 3]
2 [3 3 3]
3 [3 3 0]
4 [3 0 3]
5 [2 3 3]
-------------------------
0 [2 0 0]
1 [3 3 3]
2 [3 3 3]
3 [3 0 3]
4 [0 3 2]

Given the complexity, it seems there is no easy solution for this?

Parallel execution to reproduce the same dataset is possible IF AND ONLY IF then number of threads/processes remains the same, as per docs here and here, but this was not the intention.

1 Answers

Given that you're parallelising over rows, I'd suggest using SeedSequence.spawn_key to keep determinism within rows. For example, doing:

from numpy.random import Generator, PCG64, SeedSequence

seed = 42

for row in range(7):
  g = Generator(PCG64(SeedSequence(seed, spawn_key=[row])))
  print(g.random())

This way you don't need to worry about knowing how many "random values" are consumed by each generator invocation.

Maybe a more verbose example would help:

def make_row(row_num):
  g = Generator(PCG64(SeedSequence(seed, spawn_key=[row_num])))
  return g.random()

# one process
for i in range(1000):
  print(make_row(i))

# multiple processes:
with multiprocessing.Pool() as pool:
  pool.map(make_row, range(1000))
Related