<div dir="ltr">I have incorporated the feedback from this thread, and have significantly altered the proposal. I think this version will be more palatable to everyone.<div><br></div><div>  <a href="https://github.com/numpy/numpy/pull/11356" target="_blank">https://github.com/numpy/numpy/pull/11356</a><br></div><div>  <a href="https://github.com/rkern/numpy/blob/nep/rng-clarification/doc/neps/nep-0019-rng-policy.rst" target="_blank">https://github.com/rkern/numpy/blob/nep/rng-clarification/doc/neps/nep-0019-rng-policy.rst</a></div><div><br></div><div>I'm pretty sure that Kevin Sheppard's prototype already implements the broad strokes of my proposal (seriously, he thinks of everything; I'm just playing catch up), so I don't think there is any technical risk. I think it's just a matter of the fine details of shoving this into numpy.random per se rather than a third party package.</div><div><br></div><div>  <a href="https://bashtage.github.io/randomgen/devel/legacy.html" target="_blank">https://bashtage.github.io/randomgen/devel/legacy.html</a><br></div><div><div><br></div><div>---</div><div><br></div><div><div>==============================</div><div>Random Number Generator Policy</div><div>==============================</div><div><br></div><div>:Author: Robert Kern <<a href="mailto:robert.kern@gmail.com">robert.kern@gmail.com</a>></div><div>:Status: Draft</div><div>:Type: Standards Track</div><div>:Created: 2018-05-24</div><div><br></div><div><br></div><div>Abstract</div><div>--------</div><div><br></div><div>For the past decade, NumPy has had a strict backwards compatibility policy for</div><div>the number stream of all of its random number distributions.  Unlike other</div><div>numerical components in ``numpy``, which are usually allowed to return</div><div>different when results when they are modified if they remain correct, we have</div><div>obligated the random number distributions to always produce the exact same</div><div>numbers in every version.  The objective of our stream-compatibility guarantee</div><div>was to provide exact reproducibility for simulations across numpy versions in</div><div>order to promote reproducible research.  However, this policy has made it very</div><div>difficult to enhance any of the distributions with faster or more accurate</div><div>algorithms.  After a decade of experience and improvements in the surrounding</div><div>ecosystem of scientific software, we believe that there are now better ways to</div><div>achieve these objectives.  We propose relaxing our strict stream-compatibility</div><div>policy to remove the obstacles that are in the way of accepting contributions</div><div>to our random number generation capabilities.</div><div><br></div><div><br></div><div>The Status Quo</div><div>--------------</div><div><br></div><div>Our current policy, in full:</div><div><br></div><div>    A fixed seed and a fixed series of calls to ``RandomState`` methods using the</div><div>    same parameters will always produce the same results up to roundoff error</div><div>    except when the values were incorrect.  Incorrect values will be fixed and</div><div>    the NumPy version in which the fix was made will be noted in the relevant</div><div>    docstring.  Extension of existing parameter ranges and the addition of new</div><div>    parameters is allowed as long the previous behavior remains unchanged.</div><div><br></div><div>This policy was first instated in Nov 2008 (in essence; the full set of weasel</div><div>words grew over time) in response to a user wanting to be sure that the</div><div>simulations that formed the basis of their scientific publication could be</div><div>reproduced years later, exactly, with whatever version of ``numpy`` that was</div><div>current at the time.  We were keen to support reproducible research, and it was</div><div>still early in the life of ``numpy.random``.  We had not seen much cause to</div><div>change the distribution methods all that much.</div><div><br></div><div>We also had not thought very thoroughly about the limits of what we really</div><div>could promise (and by “we” in this section, we really mean Robert Kern, let’s</div><div>be honest).  Despite all of the weasel words, our policy overpromises</div><div>compatibility.  The same version of ``numpy`` built on different platforms, or</div><div>just in a different way could cause changes in the stream, with varying degrees</div><div>of rarity.  The biggest is that the ``.multivariate_normal()`` method relies on</div><div>``numpy.linalg`` functions.  Even on the same platform, if one links ``numpy``</div><div>with a different LAPACK, ``.multivariate_normal()`` may well return completely</div><div>different results.  More rarely, building on a different OS or CPU can cause</div><div>differences in the stream.  We use C ``long`` integers internally for integer</div><div>distribution (it seemed like a good idea at the time), and those can vary in</div><div>size depending on the platform.  Distribution methods can overflow their</div><div>internal C ``longs`` at different breakpoints depending on the platform and</div><div>cause all of the random variate draws that follow to be different.</div><div><br></div><div>And even if all of that is controlled, our policy still does not provide exact</div><div>guarantees across versions.  We still do apply bug fixes when correctness is at</div><div>stake.  And even if we didn’t do that, any nontrivial program does more than</div><div>just draw random numbers.  They do computations on those numbers, transform</div><div>those with numerical algorithms from the rest of ``numpy``, which is not</div><div>subject to so strict a policy.  Trying to maintain stream-compatibility for our</div><div>random number distributions does not help reproducible research for these</div><div>reasons.</div><div><br></div><div>The standard practice now for bit-for-bit reproducible research is to pin all</div><div>of the versions of code of your software stack, possibly down to the OS itself.</div><div>The landscape for accomplishing this is much easier today than it was in 2008.</div><div>We now have ``pip``.  We now have virtual machines.  Those who need to</div><div>reproduce simulations exactly now can (and ought to) do so by using the exact</div><div>same version of ``numpy``.  We do not need to maintain stream-compatibility</div><div>across ``numpy`` versions to help them.</div><div><br></div><div>Our stream-compatibility guarantee has hindered our ability to make</div><div>improvements to ``numpy.random``.  Several first-time contributors have</div><div>submitted PRs to improve the distributions, usually by implementing a faster,</div><div>or more accurate algorithm than the one that is currently there.</div><div>Unfortunately, most of them would have required breaking the stream to do so.</div><div>Blocked by our policy, and our inability to work around that policy, many of</div><div>those contributors simply walked away.</div><div><br></div><div><br></div><div>Implementation</div><div>--------------</div><div><br></div><div>Work on a proposed new PRNG subsystem is already underway in the randomgen_</div><div>project.  The specifics of the new design are out of scope for this NEP and up</div><div>for much discussion, but we will discuss general policies that will guide the</div><div>evolution of whatever code is adopted.  We will also outline just a few of the</div><div>requirements that such a new system must have to support the policy proposed in</div><div>this NEP.</div><div><br></div><div>First, we will maintain API source compatibility just as we do with the rest of</div><div>``numpy``.  If we *must* make a breaking change, we will only do so with an</div><div>appropriate deprecation period and warnings.</div><div><br></div><div>Second, breaking stream-compatibility in order to introduce new features or</div><div>improve performance will be *allowed* with *caution*.  Such changes will be</div><div>considered features, and as such will be no faster than the standard release</div><div>cadence of features (i.e. on ``X.Y`` releases, never ``X.Y.Z``).  Slowness will</div><div>not be considered a bug for this purpose.  Correctness bug fixes that break</div><div>stream-compatibility can happen on bugfix releases, per usual, but developers</div><div>should consider if they can wait until the next feature release.  We encourage</div><div>developers to strongly weight user’s pain from the break in</div><div>stream-compatibility against the improvements.  One example of a worthwhile</div><div>improvement would be to change algorithms for a significant increase in</div><div>performance, for example, moving from the `Box-Muller transform</div><div><<a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform">https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform</a>>`_ method of</div><div>Gaussian variate generation to the faster `Ziggurat algorithm</div><div><<a href="https://en.wikipedia.org/wiki/Ziggurat_algorithm">https://en.wikipedia.org/wiki/Ziggurat_algorithm</a>>`_.  An example of a</div><div>discouraged improvement would be tweaking the Ziggurat tables just a little bit</div><div>for a small performance improvement.</div><div><br></div><div>Any new design for the RNG subsystem will provide a choice of different core</div><div>uniform PRNG algorithms.  A promising design choice is to make these core</div><div>uniform PRNGs their own lightweight objects with a minimal set of methods</div><div>(randomgen_ calls them “basic RNGs”).  The broader set of non-uniform</div><div>distributions will be its own class that holds a reference to one of these core</div><div>uniform PRNG objects and simply delegates to the core uniform PRNG object when</div><div>it needs uniform random numbers.  To borrow an example from randomgen_, the</div><div>class ``MT19937`` is a basic RNG that implements the classic Mersenne Twister</div><div>algorithm.  The class ``RandomGenerator`` wraps around the basic RNG to provide</div><div>all of the non-uniform distribution methods::</div><div><br></div><div>    # This is not the only way to instantiate this object.</div><div>    # This is just handy for demonstrating the delegation.</div><div>    >>> brng = MT19937(seed)</div><div>    >>> rg = RandomGenerator(brng)</div><div>    >>> x = rg.standard_normal(10)</div><div><br></div><div>We will be more strict about a select subset of methods on these basic RNG</div><div>objects.  They MUST guarantee stream-compatibility for a specified set</div><div>of methods which are chosen to make it easier to compose them to build other</div><div>distributions and which are needed to abstract over the implementation details</div><div>of the variety of core PRNG algorithms.  Namely,</div><div><br></div><div>    * ``.bytes()``</div><div>    * ``.random_uintegers()``</div><div>    * ``.random_sample()``</div><div><br></div><div>The distributions class (``RandomGenerator``) SHOULD have all of the same</div><div>distribution methods as ``RandomState`` with close-enough function signatures</div><div>such that almost all code that currently works with ``RandomState`` instances</div><div>will work with ``RandomGenerator`` instances (ignoring the precise stream</div><div>values).  Some variance will be allowed for integer distributions: in order to</div><div>avoid some of the cross-platform problems described above, these SHOULD be</div><div>rewritten to work with ``uint64`` numbers on all platforms.</div><div><br></div><div>.. _randomgen: <a href="https://github.com/bashtage/randomgen">https://github.com/bashtage/randomgen</a></div><div><br></div><div><br></div><div>Supporting Unit Tests</div><div>:::::::::::::::::::::</div><div><br></div><div>Because we did make a strong stream-compatibility guarantee early in numpy’s</div><div>life, reliance on stream-compatibility has grown beyond reproducible</div><div>simulations.  One use case that remains for stream-compatibility across numpy</div><div>versions is to use pseudorandom streams to generate test data in unit tests.</div><div>With care, many of the cross-platform instabilities can be avoided in the</div><div>context of small unit tests.</div><div><br></div><div>The new PRNG subsystem MUST provide a second, legacy distributions class that</div><div>uses the same implementations of the distribution methods as the current</div><div>version of ``numpy.random.RandomState``.  The methods of this class will keep</div><div>the same strict stream-compatibility guarantees.  It is intended that this</div><div>class will no longer be modified, except to keep it working when numpy</div><div>internals change.  All new development should go into the primary distributions</div><div>class.  The purpose of ``RandomState`` will be documented as providing certain</div><div>fixed functionality for backwards compatibility and stable numbers for the</div><div>limited purpose of unit testing, and not making whole programs reproducible</div><div>across numpy versions.</div><div><br></div><div>This legacy distributions class MUST be accessible under the name</div><div>``numpy.random.RandomState`` for backwards compatibility.  All current ways of</div><div>instantiating ``numpy.random.RandomState`` with a given state should</div><div>instantiate the Mersenne Twister basic RNG with the same state.  The legacy</div><div>distributions class MUST be capable of accepting other basic RNGs.  The purpose</div><div>here is to ensure that one can write a program with a consistent basic RNG</div><div>state with a mixture of libraries that may or may not have upgraded from</div><div>``RandomState``.  Instances of the legacy distributions class MUST respond</div><div>``True`` to ``isinstance(rg, numpy.random.RandomState)`` because there is</div><div>current utility code that relies on that check.  Similarly, old pickles of</div><div>``numpy.random.RandomState`` instances MUST unpickle correctly.</div><div><br></div><div><br></div><div>``numpy.random.*``</div><div>::::::::::::::::::</div><div><br></div><div>The preferred best practice for getting reproducible pseudorandom numbers is to</div><div>instantiate a generator object with a seed and pass it around.  The implicit</div><div>global ``RandomState`` behind the ``numpy.random.*`` convenience functions can</div><div>cause problems, especially when threads or other forms of concurrency are</div><div>involved.  Global state is always problematic.  We categorically recommend</div><div>avoiding using the convenience functions when reproducibility is involved.</div><div><br></div><div>That said, people do use them and use ``numpy.random.seed()`` to control the</div><div>state underneath them.  It can be hard to categorize and count API usages</div><div>consistently and usefully, but a very common usage is in unit tests where many</div><div>of the problems of global state are less likely.</div><div><br></div><div>The initial release of the new PRNG subsystem MUST leave these convenience</div><div>functions as aliases to the methods on a global ``RandomState`` that is</div><div>initialized with a Mersenne Twister basic RNG object.  A call to</div><div>``numpy.random.seed()`` will be forwarded to that basic RNG object.  In order</div><div>to allow certain workarounds, it MUST be possible to replace the basic RNG</div><div>underneath the global ``RandomState`` with any other basic RNG object (we leave</div><div>the precise API details up to the new subsystem).  Calling ``numpy.random.seed()``</div><div>thereafter SHOULD just pass the given seed to the current basic RNG object and</div><div>not attempt to reset the basic RNG to the Mersenne Twister.  The global</div><div>``RandomState`` instance MUST be accessible by the name</div><div>``numpy.random.mtrand._rand``: Robert Kern long ago promised ``scikit-learn``</div><div>that this name would be stable.  Whoops.</div><div><br></div><div>The set of ``numpy.random.*`` convenience functions SHALL remain the same as</div><div>they currently are.  They SHALL be aliases to the ``RandomState`` methods and</div><div>not the new less-stable distributions class (``RandomGenerator``, in the</div><div>examples above). Users who want to get the fastest, best distributions can</div><div>follow best practices and instantiate generator objects explicitly.</div><div><br></div><div>After we have experience with the new PRNG subsystem, we can and should revisit</div><div>these issues in future NEPs.</div><div><br></div><div><br></div><div>Alternatives</div><div>------------</div><div><br></div><div>Versioning</div><div>::::::::::</div><div><br></div><div>For a long time, we considered that the way to allow algorithmic improvements</div><div>while maintaining the stream was to apply some form of versioning.  That is,</div><div>every time we make a stream change in one of the distributions, we increment</div><div>some version number somewhere.  ``numpy.random`` would keep all past versions</div><div>of the code, and there would be a way to get the old versions.</div><div><br></div><div>We will not be doing this.  If one needs to get the exact bit-for-bit results</div><div>from a given version of ``numpy``, whether one uses random numbers or not, one</div><div>should use the exact version of ``numpy``.</div><div><br></div><div>Proposals of how to do RNG versioning varied widely, and we will not</div><div>exhaustively list them here.  We spent years going back and forth on these</div><div>designs and were not able to find one that sufficed.  Let that time lost, and</div><div>more importantly, the contributors that we lost while we dithered, serve as</div><div>evidence against the notion.</div><div><br></div><div>Concretely, adding in versioning makes maintenance of ``numpy.random``</div><div>difficult.  Necessarily, we would be keeping lots of versions of the same code</div><div>around.  Adding a new algorithm safely would still be quite hard.</div><div><br></div><div>But most importantly, versioning is fundamentally difficult to *use* correctly.</div><div>We want to make it easy and straightforward to get the latest, fastest, best</div><div>versions of the distribution algorithms; otherwise, what's the point?  The way</div><div>to make that easy is to make the latest the default.  But the default will</div><div>necessarily change from release to release, so the user’s code would need to be</div><div>altered anyway to specify the specific version that one wants to replicate.</div><div><br></div><div>Adding in versioning to maintain stream-compatibility would still only provide</div><div>the same level of stream-compatibility that we currently do, with all of the</div><div>limitations described earlier.  Given that the standard practice for such needs</div><div>is to pin the release of ``numpy`` as a whole, versioning ``RandomState`` alone</div><div>is superfluous.</div><div><br></div><div><br></div><div>``StableRandom``</div><div>::::::::::::::::</div><div><br></div><div>A previous version of this NEP proposed to leave ``RandomState`` completely</div><div>alone for a deprecation period and build the new subsystem alongside with new</div><div>names.  To satisfy the unit testing use case, it proposed introducing a small</div><div>distributions class nominally called ``StableRandom``. It would have provided</div><div>a small subset of distribution methods that were considered most useful in unit</div><div>testing, but not the full set such that it would be too likely to be used</div><div>outside of the testing context.</div><div><br></div><div>During discussion about this proposal, it became apparent that there was no</div><div>satisfactory subset.  At least some projects used a fairly broad selection of</div><div>the ``RandomState`` methods in unit tests.</div><div><br></div><div>Downstream project owners would have been forced to modify their code to</div><div>accomodate the new PRNG subsystem.  Some modifications might be simply</div><div>mechanical, but the bulk of the work would have been tedious churn for no</div><div>positive improvement to the downstream project, just avoiding being broken.</div><div><br></div><div>Furthermore, under this old proposal, we would have had a quite lengthy</div><div>deprecation period where ``RandomState`` existed alongside the new system of</div><div>basic RNGs and distribution classes. Leaving the implementation of</div><div>``RandomState`` fixed meant that it could not use the new basic RNG state</div><div>objects.  Developing programs that use a mixture of libraries that have and</div><div>have not upgraded would require managing two sets of PRNG states.  This would</div><div>notionally have been time-limited, but we intended the deprecation to be very</div><div>long.</div><div><br></div><div>The current proposal solves all of these problems.  All current usages of</div><div>``RandomState`` will continue to work in perpetuity, though some may be</div><div>discouraged through documentation.  Unit tests can continue to use the full</div><div>complement of ``RandomState`` methods.  Mixed ``RandomState/RandomGenerator``</div><div>code can safely share the common basic RNG state.  Unmodified ``RandomState``</div><div>code can make use of the new features of alternative basic RNGs like settable</div><div>streams.</div><div><br></div><div><br></div><div>Discussion</div><div>----------</div><div><br></div><div>- `NEP discussion <<a href="https://mail.python.org/pipermail/numpy-discussion/2018-June/078126.html">https://mail.python.org/pipermail/numpy-discussion/2018-June/078126.html</a>>`_</div><div>- `Earlier discussion <<a href="https://mail.python.org/pipermail/numpy-discussion/2018-January/077608.html">https://mail.python.org/pipermail/numpy-discussion/2018-January/077608.html</a>>`_</div><div><br></div><div><br></div><div>Copyright</div><div>---------</div><div><br></div><div>This document has been placed in the public domain.</div></div><div><br></div><div><br></div>-- <br><div dir="ltr" class="gmail-m_-889815295177569256gmail_signature">Robert Kern</div></div></div>