1Parallel Random Number Generation
2=================================
3
4There are three strategies implemented that can be used to produce
5repeatable pseudo-random numbers across multiple processes (local
6or distributed).
7
8.. currentmodule:: numpy.random
9
10.. _seedsequence-spawn:
11
12`~SeedSequence` spawning
13------------------------
14
15`~SeedSequence` `implements an algorithm`_ to process a user-provided seed,
16typically as an integer of some size, and to convert it into an initial state for
17a `~BitGenerator`. It uses hashing techniques to ensure that low-quality seeds
18are turned into high quality initial states (at least, with very high
19probability).
20
21For example, `MT19937` has a state consisting of 624
22`uint32` integers. A naive way to take a 32-bit integer seed would be to just set
23the last element of the state to the 32-bit seed and leave the rest 0s. This is
24a valid state for `MT19937`, but not a good one. The Mersenne Twister
25algorithm `suffers if there are too many 0s`_. Similarly, two adjacent 32-bit
26integer seeds (i.e. ``12345`` and ``12346``) would produce very similar
27streams.
28
29`~SeedSequence` avoids these problems by using successions of integer hashes
30with good `avalanche properties`_ to ensure that flipping any bit in the input
31input has about a 50% chance of flipping any bit in the output. Two input seeds
32that are very close to each other will produce initial states that are very far
33from each other (with very high probability). It is also constructed in such
34a way that you can provide arbitrary-sized integers or lists of integers.
35`~SeedSequence` will take all of the bits that you provide and mix them
36together to produce however many bits the consuming `~BitGenerator` needs to
37initialize itself.
38
39These properties together mean that we can safely mix together the usual
40user-provided seed with simple incrementing counters to get `~BitGenerator`
41states that are (to very high probability) independent of each other. We can
42wrap this together into an API that is easy to use and difficult to misuse.
43
44.. code-block:: python
45
46  from numpy.random import SeedSequence, default_rng
47
48  ss = SeedSequence(12345)
49
50  # Spawn off 10 child SeedSequences to pass to child processes.
51  child_seeds = ss.spawn(10)
52  streams = [default_rng(s) for s in child_seeds]
53
54.. end_block
55
56Child `~SeedSequence` objects can also spawn to make grandchildren, and so on.
57Each `~SeedSequence` has its position in the tree of spawned `~SeedSequence`
58objects mixed in with the user-provided seed to generate independent (with very
59high probability) streams.
60
61.. code-block:: python
62
63  grandchildren = child_seeds[0].spawn(4)
64  grand_streams = [default_rng(s) for s in grandchildren]
65
66.. end_block
67
68This feature lets you make local decisions about when and how to split up
69streams without coordination between processes. You do not have to preallocate
70space to avoid overlapping or request streams from a common global service. This
71general "tree-hashing" scheme is `not unique to numpy`_ but not yet widespread.
72Python has increasingly-flexible mechanisms for parallelization available, and
73this scheme fits in very well with that kind of use.
74
75Using this scheme, an upper bound on the probability of a collision can be
76estimated if one knows the number of streams that you derive. `~SeedSequence`
77hashes its inputs, both the seed and the spawn-tree-path, down to a 128-bit
78pool by default. The probability that there is a collision in
79that pool, pessimistically-estimated ([1]_), will be about :math:`n^2*2^{-128}` where
80`n` is the number of streams spawned. If a program uses an aggressive million
81streams, about :math:`2^{20}`, then the probability that at least one pair of
82them are identical is about :math:`2^{-88}`, which is in solidly-ignorable
83territory ([2]_).
84
85.. [1] The algorithm is carefully designed to eliminate a number of possible
86       ways to collide. For example, if one only does one level of spawning, it
87       is guaranteed that all states will be unique. But it's easier to
88       estimate the naive upper bound on a napkin and take comfort knowing
89       that the probability is actually lower.
90
91.. [2] In this calculation, we can ignore the amount of numbers drawn from each
92       stream. Each of the PRNGs we provide has some extra protection built in
93       that avoids overlaps if the `~SeedSequence` pools differ in the
94       slightest bit. `PCG64` has :math:`2^{127}` separate cycles
95       determined by the seed in addition to the position in the
96       :math:`2^{128}` long period for each cycle, so one has to both get on or
97       near the same cycle *and* seed a nearby position in the cycle.
98       `Philox` has completely independent cycles determined by the seed.
99       `SFC64` incorporates a 64-bit counter so every unique seed is at
100       least :math:`2^{64}` iterations away from any other seed. And
101       finally, `MT19937` has just an unimaginably huge period. Getting
102       a collision internal to `SeedSequence` is the way a failure would be
103       observed.
104
105.. _`implements an algorithm`: http://www.pcg-random.org/posts/developing-a-seed_seq-alternative.html
106.. _`suffers if there are too many 0s`: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
107.. _`avalanche properties`: https://en.wikipedia.org/wiki/Avalanche_effect
108.. _`not unique to numpy`: https://www.iro.umontreal.ca/~lecuyer/myftp/papers/parallel-rng-imacs.pdf
109
110
111.. _independent-streams:
112
113Independent Streams
114-------------------
115
116`Philox` is a counter-based RNG based which generates values by
117encrypting an incrementing counter using weak cryptographic primitives. The
118seed determines the key that is used for the encryption. Unique keys create
119unique, independent streams. `Philox` lets you bypass the
120seeding algorithm to directly set the 128-bit key. Similar, but different, keys
121will still create independent streams.
122
123.. code-block:: python
124
125  import secrets
126  from numpy.random import Philox
127
128  # 128-bit number as a seed
129  root_seed = secrets.getrandbits(128)
130  streams = [Philox(key=root_seed + stream_id) for stream_id in range(10)]
131
132.. end_block
133
134This scheme does require that you avoid reusing stream IDs. This may require
135coordination between the parallel processes.
136
137
138.. _parallel-jumped:
139
140Jumping the BitGenerator state
141------------------------------
142
143``jumped`` advances the state of the BitGenerator *as-if* a large number of
144random numbers have been drawn, and returns a new instance with this state.
145The specific number of draws varies by BitGenerator, and ranges from
146:math:`2^{64}` to :math:`2^{128}`.  Additionally, the *as-if* draws also depend
147on the size of the default random number produced by the specific BitGenerator.
148The BitGenerators that support ``jumped``, along with the period of the
149BitGenerator, the size of the jump and the bits in the default unsigned random
150are listed below.
151
152+-----------------+-------------------------+-------------------------+-------------------------+
153| BitGenerator    | Period                  |  Jump Size              | Bits                    |
154+=================+=========================+=========================+=========================+
155| MT19937         | :math:`2^{19937}`       | :math:`2^{128}`         | 32                      |
156+-----------------+-------------------------+-------------------------+-------------------------+
157| PCG64           | :math:`2^{128}`         | :math:`~2^{127}` ([3]_) | 64                      |
158+-----------------+-------------------------+-------------------------+-------------------------+
159| Philox          | :math:`2^{256}`         | :math:`2^{128}`         | 64                      |
160+-----------------+-------------------------+-------------------------+-------------------------+
161
162.. [3] The jump size is :math:`(\phi-1)*2^{128}` where :math:`\phi` is the
163       golden ratio. As the jumps wrap around the period, the actual distances
164       between neighboring streams will slowly grow smaller than the jump size,
165       but using the golden ratio this way is a classic method of constructing
166       a low-discrepancy sequence that spreads out the states around the period
167       optimally. You will not be able to jump enough to make those distances
168       small enough to overlap in your lifetime.
169
170``jumped`` can be used to produce long blocks which should be long enough to not
171overlap.
172
173.. code-block:: python
174
175  import secrets
176  from numpy.random import PCG64
177
178  seed = secrets.getrandbits(128)
179  blocked_rng = []
180  rng = PCG64(seed)
181  for i in range(10):
182      blocked_rng.append(rng.jumped(i))
183
184.. end_block
185
186When using ``jumped``, one does have to take care not to jump to a stream that
187was already used. In the above example, one could not later use
188``blocked_rng[0].jumped()`` as it would overlap with ``blocked_rng[1]``. Like
189with the independent streams, if the main process here wants to split off 10
190more streams by jumping, then it needs to start with ``range(10, 20)``,
191otherwise it would recreate the same streams. On the other hand, if you
192carefully construct the streams, then you are guaranteed to have streams that
193do not overlap.
194