REF: Remove binomial_t from prng

Remvoe binomial_t from each prng_state and use a single implementation
at the level of a RandomGenerator

core_prng/common.pxd

@@ -21,6 +21,9 @@ cdef enum ConstraintType:
 
 ctypedef ConstraintType constraint_type
 
+cdef int check_constraint(double val, object name, constraint_type cons) except -1
+cdef int check_array_constraint(np.ndarray val, object name, constraint_type cons) except -1
+
 cdef extern from "src/aligned_malloc/aligned_malloc.h":
     cdef void *PyArray_realloc_aligned(void *p, size_t n);
     cdef void *PyArray_malloc_aligned(size_t n);

core_prng/distributions.pxd

@@ -36,7 +36,6 @@ cdef extern from "src/distributions/distributions.h":
         double gauss
         int has_gauss_f
         float gauss_f
-        binomial_t *binomial
 
     ctypedef prng prng_t
 
@@ -93,7 +92,7 @@ cdef extern from "src/distributions/distributions.h":
 
     long random_poisson(prng_t *prng_state, double lam) nogil
     long random_negative_binomial(prng_t *prng_state, double n, double p) nogil
-    long random_binomial(prng_t *prng_state, double p, long n) nogil
+    long random_binomial(prng_t *prng_state, double p, long n, binomial_t *binomial) nogil
     long random_logseries(prng_t *prng_state, double p) nogil
     long random_geometric_search(prng_t *prng_state, double p) nogil
     long random_geometric_inversion(prng_t *prng_state, double p) nogil

core_prng/dsfmt.pyx

@@ -6,7 +6,7 @@ import numpy as np
 cimport numpy as np
 
 from common cimport *
-from distributions cimport prng_t, binomial_t
+from distributions cimport prng_t
 from core_prng.entropy import random_entropy
 import core_prng.pickle
 cimport entropy
@@ -88,7 +88,6 @@ cdef class DSFMT:
         self.rng_state.buffered_uniforms = <double *>PyArray_calloc_aligned(DSFMT_N64, sizeof(double))
         self.rng_state.buffer_loc = DSFMT_N64
         self._prng = <prng_t *>malloc(sizeof(prng_t))
-        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed)
         self._prng.state = <void *>self.rng_state
         self._prng.next_uint64 = &dsfmt_uint64
@@ -114,7 +113,6 @@ cdef class DSFMT:
         PyArray_free_aligned(self.rng_state.state)
         PyArray_free_aligned(self.rng_state.buffered_uniforms)
         free(self.rng_state)
-        free(self._prng.binomial)
         free(self._prng)
 
     cdef _reset_state_variables(self):

core_prng/generator.pyx

@@ -3,19 +3,17 @@
 import operator 
 import warnings
 
-import numpy as np
-cimport numpy as np
-cimport cython
 from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
 from cpython cimport Py_INCREF, PyComplex_RealAsDouble, PyComplex_ImagAsDouble, PyComplex_FromDoubles
 from common cimport *
 from distributions cimport *
 from bounded_integers cimport *
 from libc cimport string
 from libc.stdlib cimport malloc, free
-
 cimport numpy as np
 import numpy as np
+cimport cython
+
 from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
 
 from common cimport *
@@ -58,6 +56,7 @@ cdef class RandomGenerator:
     """
     cdef public object __core_prng
     cdef prng_t *_prng
+    cdef binomial_t *_binomial
     cdef object lock
     poisson_lam_max = POISSON_LAM_MAX
 
@@ -71,11 +70,15 @@ cdef class RandomGenerator:
         if not PyCapsule_IsValid(capsule, anon_name):
             raise ValueError("Invalid prng. The prng must be instantized.")
         self._prng = <prng_t *> PyCapsule_GetPointer(capsule, anon_name)
+        self._binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.lock = Lock()
         with self.lock:
             self._prng.has_gauss = 0
             self._prng.has_gauss_f = 0
 
+    def __dealloc__(self):
+        free(self._binomial)
+
     def __repr__(self):
         return self.__str__() + ' at 0x{:X}'.format(id(self))
 
@@ -3244,10 +3247,64 @@ cdef class RandomGenerator:
         # answer = 0.38885, or 38%.
 
         """
-        return disc(&random_binomial, self._prng, size, self.lock, 1, 1,
-                    p, 'p', CONS_BOUNDED_0_1_NOTNAN,
-                    n, 'n', CONS_NON_NEGATIVE,
-                    0.0, '', CONS_NONE)
+
+        # Uses a custom implementation since self._binomial is required
+        cdef double _dp = 0
+        cdef long _in = 0
+        cdef bint is_scalar = True
+        cdef np.npy_intp i, cnt
+        cdef np.ndarray randoms
+        cdef np.int_t *randoms_data
+        cdef np.broadcast it
+
+        p_arr = <np.ndarray>np.PyArray_FROM_OTF(p, np.NPY_DOUBLE, np.NPY_ALIGNED)
+        is_scalar = is_scalar and np.PyArray_NDIM(p_arr) == 0
+        n_arr = <np.ndarray>np.PyArray_FROM_OTF(n, np.NPY_LONG, np.NPY_ALIGNED)
+        is_scalar = is_scalar and np.PyArray_NDIM(n_arr) == 0
+
+        if not is_scalar:
+            check_array_constraint(p_arr, 'p', CONS_BOUNDED_0_1_NOTNAN)
+            check_array_constraint(n_arr, 'n', CONS_NON_NEGATIVE)
+            if size is not None:
+                randoms = <np.ndarray>np.empty(size, np.int)
+            else:
+                it = np.PyArray_MultiIterNew2(p_arr, n_arr)
+                randoms = <np.ndarray>np.empty(it.shape, np.int)
+
+            randoms_data = <np.int_t *>np.PyArray_DATA(randoms)
+            cnt = np.PyArray_SIZE(randoms)
+
+            it = np.PyArray_MultiIterNew3(randoms, p_arr, n_arr)
+            with self.lock, nogil:
+                for i in range(cnt):
+                    _dp = (<double*>np.PyArray_MultiIter_DATA(it, 1))[0]
+                    _in = (<long*>np.PyArray_MultiIter_DATA(it, 2))[0]
+                    (<long*>np.PyArray_MultiIter_DATA(it, 0))[0] = random_binomial(self._prng, _dp, _in, self._binomial)
+
+                    np.PyArray_MultiIter_NEXT(it)
+
+            return randoms
+
+        _dp = PyFloat_AsDouble(p)
+        _in = PyInt_AsLong(n)
+        check_constraint(_dp, 'p', CONS_BOUNDED_0_1_NOTNAN)
+        check_constraint(<double>_in, 'n', CONS_NON_NEGATIVE)
+
+        if size is None:
+            with self.lock:
+                return random_binomial(self._prng, _dp, _in, self._binomial)
+
+        randoms = <np.ndarray>np.empty(size, np.int)
+        cnt = np.PyArray_SIZE(randoms)
+        randoms_data =  <np.int_t *>np.PyArray_DATA(randoms)
+
+        with self.lock, nogil:
+            for i in range(cnt):
+                randoms_data[i] = random_binomial(self._prng, _dp, _in,
+                                                  self._binomial)
+
+        return randoms
+
 
     def negative_binomial(self, n, p, size=None):
         """
@@ -4000,7 +4057,8 @@ cdef class RandomGenerator:
                 Sum = 1.0
                 dn = n
                 for j in range(d-1):
-                    mnix[i+j] = random_binomial(self._prng, pix[j]/Sum, dn)
+                    mnix[i+j] = random_binomial(self._prng, pix[j]/Sum, dn,
+                                                self._binomial)
                     dn = dn - mnix[i+j]
                     if dn <= 0:
                         break

core_prng/mt19937.pyx

@@ -7,7 +7,7 @@ import numpy as np
 cimport numpy as np
 
 from common cimport *
-from distributions cimport prng_t, binomial_t
+from distributions cimport prng_t
 import core_prng.pickle
 from core_prng.entropy import random_entropy
 
@@ -61,7 +61,6 @@ cdef class MT19937:
     def __init__(self, seed=None):
         self.rng_state = <mt19937_state *>malloc(sizeof(mt19937_state))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
-        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed)
 
         self._prng.state = <void *>self.rng_state
@@ -75,7 +74,6 @@ cdef class MT19937:
 
     def __dealloc__(self):
         free(self.rng_state)
-        free(self._prng.binomial)
         free(self._prng)
 
     cdef _reset_state_variables(self):

core_prng/pcg32.pyx

@@ -5,7 +5,7 @@ import numpy as np
 cimport numpy as np
 
 from common cimport *
-from distributions cimport prng_t, binomial_t
+from distributions cimport prng_t
 from core_prng.entropy import random_entropy
 import core_prng.pickle
 cimport entropy
@@ -68,7 +68,6 @@ cdef class PCG32:
         self.rng_state = <pcg32_state *>malloc(sizeof(pcg32_state))
         self.rng_state.pcg_state = <pcg32_random_t *>malloc(sizeof(pcg32_random_t))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
-        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed, inc)
 
         self._prng.state = <void *>self.rng_state
@@ -94,7 +93,6 @@ cdef class PCG32:
 
     def __dealloc__(self):
         free(self.rng_state)
-        free(self._prng.binomial)
         free(self._prng)
 
     cdef _reset_state_variables(self):

core_prng/pcg64.pyx

@@ -5,7 +5,7 @@ import numpy as np
 cimport numpy as np
 
 from common cimport *
-from distributions cimport prng_t, binomial_t
+from distributions cimport prng_t
 from core_prng.entropy import random_entropy
 import core_prng.pickle
 cimport entropy
@@ -83,7 +83,6 @@ cdef class PCG64:
         self.rng_state = <pcg64_state *>malloc(sizeof(pcg64_state))
         self.rng_state.pcg_state = <pcg64_random_t *>malloc(sizeof(pcg64_random_t))
         self._prng = <prng_t *>malloc(sizeof(prng_t))
-        self._prng.binomial = <binomial_t *>malloc(sizeof(binomial_t))
         self.seed(seed, inc)
 
         self._prng.state = <void *>self.rng_state
@@ -109,7 +108,6 @@ cdef class PCG64:
 
     def __dealloc__(self):
         free(self.rng_state)
-        free(self._prng.binomial)
         free(self._prng)
 
     cdef _reset_state_variables(self):

core_prng/philox.pyx

@@ -4,7 +4,7 @@ from cpython.pycapsule cimport PyCapsule_New
 import numpy as np
 
 from common cimport *
-from distributions cimport prng_t, binomial_t
+from distributions cimport prng_t
 from core_prng.entropy import random_entropy
 import core_prng.pickle
 cimport entropy
@@ -76,7 +76,6 @@ cdef class Philox:
         self.rng_state.key = <philox4x64_key_t *> malloc(
             sizeof(philox4x64_key_t))
         self._prng = <prng_t *> malloc(sizeof(prng_t))
-        self._prng.binomial = <binomial_t *> malloc(sizeof(binomial_t))
         self.seed(seed, counter, key)
 
         self._prng.state = <void *> self.rng_state
@@ -104,7 +103,6 @@ cdef class Philox:
         free(self.rng_state.ctr)
         free(self.rng_state.key)
         free(self.rng_state)
-        free(self._prng.binomial)
         free(self._prng)
 
     cdef _reset_state_variables(self):