extending.pyx 2.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677
  1. #cython: language_level=3
  2. from libc.stdint cimport uint32_t
  3. from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
  4. import numpy as np
  5. cimport numpy as np
  6. cimport cython
  7. from numpy.random cimport bitgen_t
  8. from numpy.random import PCG64
  9. np.import_array()
  10. @cython.boundscheck(False)
  11. @cython.wraparound(False)
  12. def uniform_mean(Py_ssize_t n):
  13. cdef Py_ssize_t i
  14. cdef bitgen_t *rng
  15. cdef const char *capsule_name = "BitGenerator"
  16. cdef double[::1] random_values
  17. cdef np.ndarray randoms
  18. x = PCG64()
  19. capsule = x.capsule
  20. if not PyCapsule_IsValid(capsule, capsule_name):
  21. raise ValueError("Invalid pointer to anon_func_state")
  22. rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name)
  23. random_values = np.empty(n)
  24. # Best practice is to acquire the lock whenever generating random values.
  25. # This prevents other threads from modifying the state. Acquiring the lock
  26. # is only necessary if the GIL is also released, as in this example.
  27. with x.lock, nogil:
  28. for i in range(n):
  29. random_values[i] = rng.next_double(rng.state)
  30. randoms = np.asarray(random_values)
  31. return randoms.mean()
  32. # This function is declared nogil so it can be used without the GIL below
  33. cdef uint32_t bounded_uint(uint32_t lb, uint32_t ub, bitgen_t *rng) nogil:
  34. cdef uint32_t mask, delta, val
  35. mask = delta = ub - lb
  36. mask |= mask >> 1
  37. mask |= mask >> 2
  38. mask |= mask >> 4
  39. mask |= mask >> 8
  40. mask |= mask >> 16
  41. val = rng.next_uint32(rng.state) & mask
  42. while val > delta:
  43. val = rng.next_uint32(rng.state) & mask
  44. return lb + val
  45. @cython.boundscheck(False)
  46. @cython.wraparound(False)
  47. def bounded_uints(uint32_t lb, uint32_t ub, Py_ssize_t n):
  48. cdef Py_ssize_t i
  49. cdef bitgen_t *rng
  50. cdef uint32_t[::1] out
  51. cdef const char *capsule_name = "BitGenerator"
  52. x = PCG64()
  53. out = np.empty(n, dtype=np.uint32)
  54. capsule = x.capsule
  55. if not PyCapsule_IsValid(capsule, capsule_name):
  56. raise ValueError("Invalid pointer to anon_func_state")
  57. rng = <bitgen_t *>PyCapsule_GetPointer(capsule, capsule_name)
  58. with x.lock, nogil:
  59. for i in range(n):
  60. out[i] = bounded_uint(lb, ub, rng)
  61. return np.asarray(out)