Merge pull request numpy#13163 from mattip/randomgen

ENH: randomgen

This merges randomgen into numpy, which was originally developed at https://github.com/bashtage/randomgen and provides a new and improved API for random number generation with much new and improved functionality.

.circleci/config.yml

@@ -20,7 +20,7 @@ jobs:
           command: |
             python3 -m venv venv
             . venv/bin/activate
-            pip install cython sphinx==1.8.5 matplotlib
+            pip install cython sphinx==1.8.5 matplotlib ipython
             sudo apt-get update
             sudo apt-get install -y graphviz texlive-fonts-recommended texlive-latex-recommended texlive-latex-extra texlive-generic-extra latexmk texlive-xetex
 

.gitignore

@@ -168,12 +168,16 @@ numpy/core/src/umath/test_rational.c
 numpy/core/src/umath/umath_tests.c
 numpy/distutils/__config__.py
 numpy/linalg/umath_linalg.c
-doc/source/reference/generated
+doc/source/**/generated/
 benchmarks/results
 benchmarks/html
 benchmarks/env
 benchmarks/numpy
 # cythonized files
 cythonize.dat
-numpy/random/mtrand/mtrand.c
-numpy/random/mtrand/randint_helpers.pxi
+numpy/random/_mtrand/_mtrand.c
+numpy/random/*.c
+numpy/random/legacy/*.c
+numpy/random/_mtrand/randint_helpers.pxi
+numpy/random/bounded_integers.pyx
+numpy/random/bounded_integers.pxd

MANIFEST.in

@@ -9,7 +9,7 @@ include pytest.ini
 include *.txt
 include README.md
 include site.cfg.example
-recursive-include numpy/random/mtrand *.pyx *.pxd
+recursive-include numpy/random *.pyx *.pxd *.pyx.in *.pxd.in
 # Add build support that should go in sdist, but not go in bdist/be installed
 # Note that sub-directories that don't have __init__ are apparently not
 # included by 'recursive-include', so list those separately

benchmarks/benchmarks/bench_random.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 
+from numpy.random import RandomState, Generator
 
 class Random(Benchmark):
     params = ['normal', 'uniform', 'weibull 1', 'binomial 10 0.5',
@@ -69,14 +70,102 @@ def time_randint_slow(self, name):
 class Permutation(Benchmark):
     def setup(self):
         self.n = 10000
-        self.a_1d = np.random.random_sample(self.n)
-        self.a_2d = np.random.random_sample((self.n, 2))
-    
+        self.a_1d = np.random.random(self.n)
+        self.a_2d = np.random.random((self.n, 2))
+
     def time_permutation_1d(self):
         np.random.permutation(self.a_1d)
 
     def time_permutation_2d(self):
-        np.random.permutation(self.a_2d)        
+        np.random.permutation(self.a_2d)
 
     def time_permutation_int(self):
         np.random.permutation(self.n)
+
+nom_size = 100000
+
+class RNG(Benchmark):
+    param_names = ['rng']
+    params = ['DSFMT', 'PCG64', 'PCG32', 'MT19937', 'Xoshiro256',
+              'Xoshiro512', 'Philox', 'ThreeFry', 'numpy']
+
+    def setup(self, bitgen):
+        if bitgen == 'numpy':
+            self.rg = np.random.RandomState()
+        else:
+            self.rg = Generator(getattr(np.random, bitgen)())
+        self.rg.random()
+        self.int32info = np.iinfo(np.int32)
+        self.uint32info = np.iinfo(np.uint32)
+        self.uint64info = np.iinfo(np.uint64)
+
+    def time_raw(self, bitgen):
+        if bitgen == 'numpy':
+            self.rg.random_integers(self.int32info.max, size=nom_size)
+        else:
+            self.rg.integers(self.int32info.max, size=nom_size, endpoint=True)
+
+    def time_32bit(self, bitgen):
+        min, max = self.uint32info.min, self.uint32info.max
+        if bitgen == 'numpy':
+            self.rg.randint(min, max + 1, nom_size, dtype=np.uint32)
+        else:
+            self.rg.integers(min, max + 1, nom_size, dtype=np.uint32)
+
+    def time_64bit(self, bitgen):
+        min, max = self.uint64info.min, self.uint64info.max
+        if bitgen == 'numpy':
+            self.rg.randint(min, max + 1, nom_size, dtype=np.uint64)
+        else:
+            self.rg.integers(min, max + 1, nom_size, dtype=np.uint64)
+
+    def time_normal_zig(self, bitgen):
+        self.rg.standard_normal(nom_size)
+
+class Bounded(Benchmark):
+    u8 = np.uint8
+    u16 = np.uint16
+    u32 = np.uint32
+    u64 = np.uint64
+    param_names = ['rng', 'dt_max']
+    params = [['DSFMT', 'PCG64', 'PCG32', 'MT19937','Xoshiro256',
+               'Xoshiro512', 'Philox', 'ThreeFry', 'numpy'],
+              [[u8,    95],
+               [u8,    64],  # Worst case for legacy
+               [u8,   127],  # Best case for legacy
+               [u16,   95],
+               [u16, 1024],  # Worst case for legacy
+               [u16, 1535],   # Typ. avg. case for legacy
+               [u16, 2047],  # Best case for legacy
+               [u32, 1024],  # Worst case for legacy
+               [u32, 1535],   # Typ. avg. case for legacy
+               [u32, 2047],  # Best case for legacy
+               [u64,   95],
+               [u64, 1024],  # Worst case for legacy
+               [u64, 1535],   # Typ. avg. case for legacy
+               [u64, 2047],  # Best case for legacy
+             ]]
+
+    def setup(self, bitgen, args):
+        if bitgen == 'numpy':
+            self.rg = np.random.RandomState()
+        else:
+            self.rg = Generator(getattr(np.random, bitgen)())
+        self.rg.random()
+
+    def time_bounded(self, bitgen, args):
+            """
+            Timer for 8-bit bounded values.
+
+            Parameters (packed as args)
+            ----------
+            dt : {uint8, uint16, uint32, unit64}
+                output dtype
+            max : int
+                Upper bound for range. Lower is always 0.  Must be <= 2**bits.
+            """
+            dt, max = args
+            if bitgen == 'numpy':
+                self.rg.randint(0, max + 1, nom_size, dtype=dt)
+            else:
+                self.rg.integers(0, max + 1, nom_size, dtype=dt)

doc/Makefile

@@ -46,7 +46,8 @@ help:
 	@echo "  upload USERNAME=... RELEASE=... to upload built docs to docs.scipy.org"
 
 clean:
-	-rm -rf build/* source/reference/generated
+	-rm -rf build/* 
+	find . -name generated -type d -prune -exec rm -rf "{}" ";"
 
 version-check:
 ifeq "$(GITVER)" "Unknown"

doc/neps/nep-0010-new-iterator-ufunc.rst

@@ -1877,8 +1877,8 @@ the new iterator.
 Here is one of the original functions, for reference, and some
 random image data.::
 
-    In [5]: rand1 = np.random.random_sample(1080*1920*4).astype(np.float32)
-    In [6]: rand2 = np.random.random_sample(1080*1920*4).astype(np.float32)
+    In [5]: rand1 = np.random.random(1080*1920*4).astype(np.float32)
+    In [6]: rand2 = np.random.random(1080*1920*4).astype(np.float32)
     In [7]: image1 = rand1.reshape(1080,1920,4).swapaxes(0,1)
     In [8]: image2 = rand2.reshape(1080,1920,4).swapaxes(0,1)
 

doc/neps/nep-0019-rng-policy.rst

@@ -143,11 +143,11 @@ We will be more strict about a select subset of methods on these BitGenerator
 objects.  They MUST guarantee stream-compatibility for a specified set
 of methods which are chosen to make it easier to compose them to build other
 distributions and which are needed to abstract over the implementation details
-of the variety of core PRNG algorithms.  Namely,
+of the variety of BitGenerator algorithms.  Namely,
 
     * ``.bytes()``
-    * ``integers`` (formerly ``.random_uintegers()``)
-    * ``random`` (formerly ``.random_sample()``)
+    * ``integers()`` (formerly ``.random_integers()``)
+    * ``random()`` (formerly ``.random_sample()``)
 
 The distributions class (``Generator``) SHOULD have all of the same
 distribution methods as ``RandomState`` with close-enough function signatures

doc/source/about.rst

@@ -8,7 +8,7 @@ needed for scientific computing with Python. This package contains:
 - sophisticated :ref:`(broadcasting) functions <ufuncs>`
 - basic :ref:`linear algebra functions <routines.linalg>`
 - basic :ref:`Fourier transforms <routines.fft>`
-- sophisticated :ref:`random number capabilities <routines.random>`
+- sophisticated :ref:`random number capabilities <numpyrandom>`
 - tools for integrating Fortran code
 - tools for integrating C/C++ code
 

doc/source/conf.py

@@ -19,11 +19,19 @@
 
 sys.path.insert(0, os.path.abspath('../sphinxext'))
 
-extensions = ['sphinx.ext.autodoc', 'numpydoc',
-              'sphinx.ext.intersphinx', 'sphinx.ext.coverage',
-              'sphinx.ext.doctest', 'sphinx.ext.autosummary',
-              'sphinx.ext.graphviz', 'sphinx.ext.ifconfig',
-              'matplotlib.sphinxext.plot_directive']
+extensions = [
+    'sphinx.ext.autodoc',
+    'numpydoc',
+    'sphinx.ext.intersphinx',
+    'sphinx.ext.coverage',
+    'sphinx.ext.doctest',
+    'sphinx.ext.autosummary',
+    'sphinx.ext.graphviz',
+    'sphinx.ext.ifconfig',
+    'matplotlib.sphinxext.plot_directive',
+    'IPython.sphinxext.ipython_console_highlighting',
+    'IPython.sphinxext.ipython_directive',
+]
 
 if sphinx.__version__ >= "1.4":
     extensions.append('sphinx.ext.imgmath')
@@ -234,7 +242,7 @@ def setup(app):
 # -----------------------------------------------------------------------------
 
 import glob
-autosummary_generate = glob.glob("reference/*.rst")
+autosummary_generate = True
 
 # -----------------------------------------------------------------------------
 # Coverage checker
@@ -355,3 +363,4 @@ def linkcode_resolve(domain, info):
     else:
         return "https://github.com/numpy/numpy/blob/v%s/numpy/%s%s" % (
            numpy.__version__, fn, linespec)
+

doc/source/reference/random/bit_generators/dsfmt.rst

@@ -0,0 +1,36 @@
+Double SIMD Mersenne Twister (dSFMT) 
+------------------------------------
+
+.. module:: numpy.random.dsfmt
+
+.. currentmodule:: numpy.random.dsfmt
+
+
+.. autoclass:: DSFMT
+	:exclude-members:
+
+Seeding and State
+=================
+
+.. autosummary::
+   :toctree: generated/
+
+   ~DSFMT.seed
+   ~DSFMT.state
+
+Parallel generation
+===================
+.. autosummary::
+   :toctree: generated/
+
+   ~DSFMT.jumped
+
+Extending
+=========
+.. autosummary::
+   :toctree: generated/
+
+   ~DSFMT.cffi
+   ~DSFMT.ctypes
+
+

doc/source/reference/random/bit_generators/index.rst

@@ -0,0 +1,29 @@
+.. _bit_generator:
+
+Bit Generators
+--------------
+
+.. currentmodule:: numpy.random
+
+The random values produced by :class:`~Generator`
+orignate in a BitGenerator.  The BitGenerators do not directly provide
+random numbers and only contains methods used for seeding, getting or
+setting the state, jumping or advancing the state, and for accessing
+low-level wrappers for consumption by code that can efficiently
+access the functions provided, e.g., `numba <https://numba.pydata.org>`_.
+
+Stable RNGs
+===========
+
+.. toctree::
+   :maxdepth: 1
+
+   DSFMT <dsfmt>
+   MT19937 <mt19937>
+   PCG32 <pcg32>
+   PCG64 <pcg64>
+   Philox <philox>
+   ThreeFry <threefry>
+   Xoshiro256** <xoshiro256>
+   Xoshiro512** <xoshiro512>
+

doc/source/reference/random/bit_generators/mt19937.rst

@@ -0,0 +1,35 @@
+Mersenne Twister (MT19937) 
+--------------------------
+
+.. module:: numpy.random.mt19937
+
+.. currentmodule:: numpy.random.mt19937
+
+.. autoclass:: MT19937
+	:exclude-members:
+
+Seeding and State
+=================
+
+.. autosummary::
+   :toctree: generated/
+
+   ~MT19937.seed
+   ~MT19937.state
+
+Parallel generation
+===================
+.. autosummary::
+   :toctree: generated/
+
+   ~MT19937.jumped
+
+Extending
+=========
+.. autosummary::
+   :toctree: generated/
+
+   ~MT19937.cffi
+   ~MT19937.ctypes
+
+

doc/source/reference/random/bit_generators/pcg32.rst

@@ -0,0 +1,34 @@
+Parallel Congruent Generator (32-bit, PCG32)
+--------------------------------------------
+
+.. module:: numpy.random.pcg32
+
+.. currentmodule:: numpy.random.pcg32
+
+.. autoclass:: PCG32
+	:exclude-members:
+
+Seeding and State
+=================
+
+.. autosummary::
+   :toctree: generated/
+
+   ~PCG32.seed
+   ~PCG32.state
+
+Parallel generation
+===================
+.. autosummary::
+   :toctree: generated/
+
+   ~PCG32.advance
+   ~PCG32.jumped
+
+Extending
+=========
+.. autosummary::
+   :toctree: generated/
+
+   ~PCG32.cffi
+   ~PCG32.ctypes