diff options
Diffstat (limited to 'numpy/random/examples')
| -rw-r--r-- | numpy/random/examples/__init__.py | 0 | ||||
| -rw-r--r-- | numpy/random/examples/cython/extending.pyx | 2 | ||||
| -rw-r--r-- | numpy/random/examples/cython/extending_distributions.pyx | 34 | ||||
| -rw-r--r-- | numpy/random/examples/cython/setup.py | 17 | ||||
| -rw-r--r-- | numpy/random/examples/numba/extending.py | 97 |
5 files changed, 83 insertions, 67 deletions
diff --git a/numpy/random/examples/__init__.py b/numpy/random/examples/__init__.py new file mode 100644 index 000000000..e69de29bb --- /dev/null +++ b/numpy/random/examples/__init__.py diff --git a/numpy/random/examples/cython/extending.pyx b/numpy/random/examples/cython/extending.pyx index a6a4ba4bf..2a866648d 100644 --- a/numpy/random/examples/cython/extending.pyx +++ b/numpy/random/examples/cython/extending.pyx @@ -8,7 +8,7 @@ import numpy as np cimport numpy as np cimport cython -from numpy.random.common cimport bitgen_t +from numpy.random._bit_generator cimport bitgen_t from numpy.random import PCG64 np.import_array() diff --git a/numpy/random/examples/cython/extending_distributions.pyx b/numpy/random/examples/cython/extending_distributions.pyx index 3cefec97e..d17da45c1 100644 --- a/numpy/random/examples/cython/extending_distributions.pyx +++ b/numpy/random/examples/cython/extending_distributions.pyx @@ -1,21 +1,23 @@ #!/usr/bin/env python #cython: language_level=3 """ -This file shows how the distributions that are accessed through -distributions.pxd can be used Cython code. +This file shows how the to use a BitGenerator to create a distribution. """ import numpy as np cimport numpy as np cimport cython from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer -from numpy.random.common cimport * -from numpy.random.distributions cimport random_gauss_zig +from numpy.random._bit_generator cimport bitgen_t from numpy.random import PCG64 @cython.boundscheck(False) @cython.wraparound(False) -def normals_zig(Py_ssize_t n): +def uniforms(Py_ssize_t n): + """ Create an array of `n` uniformly distributed doubles. + A 'real' distribution would want to process the values into + some non-uniform distribution + """ cdef Py_ssize_t i cdef bitgen_t *rng cdef const char *capsule_name = "BitGenerator" @@ -23,21 +25,25 @@ def normals_zig(Py_ssize_t n): x = PCG64() capsule = x.capsule + # Optional check that the capsule if from a BitGenerator if not PyCapsule_IsValid(capsule, capsule_name): raise ValueError("Invalid pointer to anon_func_state") + # Cast the pointer rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name) - random_values = np.empty(n) - # Best practice is to release GIL and acquire the lock + random_values = np.empty(n, dtype='float64') with x.lock, nogil: for i in range(n): - random_values[i] = random_gauss_zig(rng) + # Call the function + random_values[i] = rng.next_double(rng.state) randoms = np.asarray(random_values) - return randoms + return randoms + +# cython example 2 @cython.boundscheck(False) @cython.wraparound(False) -def uniforms(Py_ssize_t n): +def uint16_uniforms(Py_ssize_t n): cdef Py_ssize_t i cdef bitgen_t *rng cdef const char *capsule_name = "BitGenerator" @@ -45,15 +51,13 @@ def uniforms(Py_ssize_t n): x = PCG64() capsule = x.capsule - # Optional check that the capsule if from a BitGenerator if not PyCapsule_IsValid(capsule, capsule_name): raise ValueError("Invalid pointer to anon_func_state") - # Cast the pointer rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name) - random_values = np.empty(n) + random_values = np.empty(n, dtype='uint32') + # Best practice is to release GIL and acquire the lock with x.lock, nogil: for i in range(n): - # Call the function - random_values[i] = rng.next_double(rng.state) + random_values[i] = rng.next_uint32(rng.state) randoms = np.asarray(random_values) return randoms diff --git a/numpy/random/examples/cython/setup.py b/numpy/random/examples/cython/setup.py index 69f057ed5..315527a2d 100644 --- a/numpy/random/examples/cython/setup.py +++ b/numpy/random/examples/cython/setup.py @@ -9,15 +9,20 @@ import numpy as np from distutils.core import setup from Cython.Build import cythonize from setuptools.extension import Extension -from os.path import join +from os.path import join, abspath, dirname + +curpath = abspath(dirname(__file__)) extending = Extension("extending", - sources=['extending.pyx'], - include_dirs=[np.get_include()]) + sources=[join(curpath, 'extending.pyx')], + include_dirs=[ + np.get_include(), + join(curpath, '..', '..') + ], + ) distributions = Extension("extending_distributions", - sources=['extending_distributions.pyx', - join('..', '..', 'src', - 'distributions', 'distributions.c')], + sources=[join(curpath, 'extending_distributions.pyx'), + ], include_dirs=[np.get_include()]) extensions = [extending, distributions] diff --git a/numpy/random/examples/numba/extending.py b/numpy/random/examples/numba/extending.py index d41c2d76f..0d240596b 100644 --- a/numpy/random/examples/numba/extending.py +++ b/numpy/random/examples/numba/extending.py @@ -1,14 +1,57 @@ -import datetime as dt - import numpy as np import numba as nb from numpy.random import PCG64 +from timeit import timeit + +bit_gen = PCG64() +next_d = bit_gen.cffi.next_double +state_addr = bit_gen.cffi.state_address + +def normals(n, state): + out = np.empty(n) + for i in range((n + 1) // 2): + x1 = 2.0 * next_d(state) - 1.0 + x2 = 2.0 * next_d(state) - 1.0 + r2 = x1 * x1 + x2 * x2 + while r2 >= 1.0 or r2 == 0.0: + x1 = 2.0 * next_d(state) - 1.0 + x2 = 2.0 * next_d(state) - 1.0 + r2 = x1 * x1 + x2 * x2 + f = np.sqrt(-2.0 * np.log(r2) / r2) + out[2 * i] = f * x1 + if 2 * i + 1 < n: + out[2 * i + 1] = f * x2 + return out + +# Compile using Numba +normalsj = nb.jit(normals, nopython=True) +# Must use state address not state with numba +n = 10000 + +def numbacall(): + return normalsj(n, state_addr) + +rg = np.random.Generator(PCG64()) -x = PCG64() -f = x.ctypes.next_uint32 -s = x.ctypes.state +def numpycall(): + return rg.normal(size=n) +# Check that the functions work +r1 = numbacall() +r2 = numpycall() +assert r1.shape == (n,) +assert r1.shape == r2.shape + +t1 = timeit(numbacall, number=1000) +print('{:.2f} secs for {} PCG64 (Numba/PCG64) gaussian randoms'.format(t1, n)) +t2 = timeit(numpycall, number=1000) +print('{:.2f} secs for {} PCG64 (NumPy/PCG64) gaussian randoms'.format(t2, n)) + +# example 2 + +next_u32 = bit_gen.ctypes.next_uint32 +ctypes_state = bit_gen.ctypes.state @nb.jit(nopython=True) def bounded_uint(lb, ub, state): @@ -19,14 +62,14 @@ def bounded_uint(lb, ub, state): mask |= mask >> 8 mask |= mask >> 16 - val = f(state) & mask + val = next_u32(state) & mask while val > delta: - val = f(state) & mask + val = next_u32(state) & mask return lb + val -print(bounded_uint(323, 2394691, s.value)) +print(bounded_uint(323, 2394691, ctypes_state.value)) @nb.jit(nopython=True) @@ -36,42 +79,6 @@ def bounded_uints(lb, ub, n, state): out[i] = bounded_uint(lb, ub, state) -bounded_uints(323, 2394691, 10000000, s.value) - -g = x.cffi.next_double -cffi_state = x.cffi.state -state_addr = x.cffi.state_address - - -def normals(n, state): - out = np.empty(n) - for i in range((n + 1) // 2): - x1 = 2.0 * g(state) - 1.0 - x2 = 2.0 * g(state) - 1.0 - r2 = x1 * x1 + x2 * x2 - while r2 >= 1.0 or r2 == 0.0: - x1 = 2.0 * g(state) - 1.0 - x2 = 2.0 * g(state) - 1.0 - r2 = x1 * x1 + x2 * x2 - f = np.sqrt(-2.0 * np.log(r2) / r2) - out[2 * i] = f * x1 - if 2 * i + 1 < n: - out[2 * i + 1] = f * x2 - return out +bounded_uints(323, 2394691, 10000000, ctypes_state.value) -print(normals(10, cffi_state).var()) -# Warm up -normalsj = nb.jit(normals, nopython=True) -normalsj(1, state_addr) - -start = dt.datetime.now() -normalsj(1000000, state_addr) -ms = 1000 * (dt.datetime.now() - start).total_seconds() -print('1,000,000 Polar-transform (numba/PCG64) randoms in ' - '{ms:0.1f}ms'.format(ms=ms)) - -start = dt.datetime.now() -np.random.standard_normal(1000000) -ms = 1000 * (dt.datetime.now() - start).total_seconds() -print('1,000,000 Polar-transform (NumPy) randoms in {ms:0.1f}ms'.format(ms=ms)) |