bitgenerator被设计成可以使用高性能Python的标准工具numba和Cython进行扩展。这个 Generator 对象也可以与用户提供的位生成器一起使用,只要这些位生成器导出一小部分所需的函数。
Generator
Numba可以与CTypes或CFFI一起使用。两个小比特集的接口函数都通过迭代导出当前所有比特集。
这个例子展示了如何使用一个纯Python实现来生成gaussian样本,然后编译这个实现。随机数由 ctypes.next_double .
ctypes.next_double
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()) 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(f'{t1:.2f} secs for {n} PCG64 (Numba/PCG64) gaussian randoms') t2 = timeit(numpycall, number=1000) print(f'{t2:.2f} secs for {n} PCG64 (NumPy/PCG64) gaussian randoms')
CTypes和CFFI都允许在将distributions.c文件编译成 DLL 或 so . 下面的示例展示了如何使用更复杂的分布 examples 下面部分。
DLL
so
Cython可用于打开 PyCapsule 由位发生器提供。本例使用 PCG64 还有上面的例子。并提供了通常的数组对齐检查——使用Cython来消除代码的高性能限制。
PyCapsule
PCG64
#!/usr/bin/env python3 #cython: language_level=3 """ 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 libc.stdint cimport uint16_t, uint64_t from numpy.random cimport bitgen_t from numpy.random import PCG64 from numpy.random.c_distributions cimport ( random_standard_uniform_fill, random_standard_uniform_fill_f) @cython.boundscheck(False) @cython.wraparound(False) 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" cdef double[::1] random_values 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, dtype='float64') with x.lock, nogil: for i in range(n): # Call the function random_values[i] = rng.next_double(rng.state) randoms = np.asarray(random_values) return randoms
还可以使用 bitgen_t 结构。
bitgen_t
@cython.boundscheck(False) @cython.wraparound(False) def uint10_uniforms(Py_ssize_t n): """Uniform 10 bit integers stored as 16-bit unsigned integers""" cdef Py_ssize_t i cdef bitgen_t *rng cdef const char *capsule_name = "BitGenerator" cdef uint16_t[::1] random_values cdef int bits_remaining cdef int width = 10 cdef uint64_t buff, mask = 0x3FF x = PCG64() capsule = x.capsule if not PyCapsule_IsValid(capsule, capsule_name): raise ValueError("Invalid pointer to anon_func_state") rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name) random_values = np.empty(n, dtype='uint16') # Best practice is to release GIL and acquire the lock bits_remaining = 0 with x.lock, nogil: for i in range(n): if bits_remaining < width: buff = rng.next_uint64(rng.state) random_values[i] = buff & mask buff >>= width randoms = np.asarray(random_values) return randoms
Cython可以用来直接访问 numpy/random/c_distributions.pxd . 这需要与 npyrandom 类库位于 numpy/random/lib .
numpy/random/c_distributions.pxd
npyrandom
numpy/random/lib
def uniforms_ex(bit_generator, Py_ssize_t n, dtype=np.float64): """ Create an array of `n` uniformly distributed doubles via a "fill" function. A 'real' distribution would want to process the values into some non-uniform distribution Parameters ---------- bit_generator: BitGenerator instance n: int Output vector length dtype: {str, dtype}, optional Desired dtype, either 'd' (or 'float64') or 'f' (or 'float32'). The default dtype value is 'd' """ cdef Py_ssize_t i cdef bitgen_t *rng cdef const char *capsule_name = "BitGenerator" cdef np.ndarray randoms capsule = bit_generator.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) _dtype = np.dtype(dtype) randoms = np.empty(n, dtype=_dtype) if _dtype == np.float32: with bit_generator.lock: random_standard_uniform_fill_f(rng, n, <float*>np.PyArray_DATA(randoms)) elif _dtype == np.float64: with bit_generator.lock: random_standard_uniform_fill(rng, n, <double*>np.PyArray_DATA(randoms)) else: raise TypeError('Unsupported dtype %r for random' % _dtype) return randoms
见 延伸 numpy.random 西顿大道 有关这些示例的完整列表和 setup.py 构建c扩展模块。
setup.py
CFFI可用于直接访问 include/numpy/random/distributions.h . 需要对头文件进行一些“按摩”:
include/numpy/random/distributions.h
""" Use cffi to access any of the underlying C functions from distributions.h """ import os import numpy as np import cffi from .parse import parse_distributions_h ffi = cffi.FFI() inc_dir = os.path.join(np.get_include(), 'numpy') # Basic numpy types ffi.cdef(''' typedef intptr_t npy_intp; typedef unsigned char npy_bool; ''') parse_distributions_h(ffi, inc_dir)
一旦标头被 ffi.cdef ,可以直接从 _generator 共享对象,使用 BitGenerator.cffi 接口。
ffi.cdef
_generator
BitGenerator.cffi
# Compare the distributions.h random_standard_normal_fill to # Generator.standard_random bit_gen = np.random.PCG64() rng = np.random.Generator(bit_gen) state = bit_gen.state interface = rng.bit_generator.cffi n = 100 vals_cffi = ffi.new('double[%d]' % n) lib.random_standard_normal_fill(interface.bit_generator, n, vals_cffi) # reset the state bit_gen.state = state vals = rng.standard_normal(n) for i in range(n): assert vals[i] == vals_cffi[i]
Generator 可与用户提供的 BitGenerator 编写一个新的位生成器最简单的方法是检查一个现有位生成器的pyx文件。必须提供的关键结构是 capsule 其中包含 PyCapsule 类型的结构指针 bitgen_t ,
BitGenerator
capsule
typedef struct bitgen { void *state; uint64_t (*next_uint64)(void *st); uint32_t (*next_uint32)(void *st); double (*next_double)(void *st); uint64_t (*next_raw)(void *st); } bitgen_t;
它提供了5个指针。第一个是指向位生成器使用的数据结构的不透明指针。下三个是函数指针,返回下一个64位和32位无符号整数、下一个随机双精度和下一个原始值。最后一个函数用于测试,因此如果不需要,可以设置为下一个64位无符号整数函数。内部功能 Generator 如中所示使用此结构
bitgen_state->next_uint64(bitgen_state->state)