Merge pull request #3 from bashtage/move-binomial

REF: Remove binomial_t from prng

_randomgen/README.md

@@ -7,7 +7,14 @@ Experimental Core Pseudo Random Number Generator interface for future
 NumPy RandomState evolution.
 
 This is a library and generic interface for alternative random 
-generators in Python and NumPy. 
+generators in Python and NumPy.
+
+
+### Compatibility Warning
+Core PRNG no longer supports Box-Muller normal variates and so it not
+100% compatible with NumPy (or randomstate). Box-Muller normals are slow
+to generate and all functions which previously relied on Box-Muller
+normals now use the faster Ziggurat implementation.
 
 ## Features
 
@@ -42,14 +49,14 @@ y = rnd.standard_gamma(5.5, 10000, method='zig')
 
     * Uniforms (`random_sample`)
     * Exponentials (`standard_exponential`, both Inverse CDF and Ziggurat)
-    * Normals (`standard_normal`, both Box-Muller and Ziggurat)
-    * Standard Gammas (via `standard_gamma`, both Inverse CDF and Ziggurat)
+    * Normals (`standard_normal`)
+    * Standard Gammas (via `standard_gamma`)
 
   **WARNING**: The 32-bit generators are **experimental** and subject 
   to change.
 
   **Note**: There are _no_ plans to extend the alternative precision 
-  generation to all random number types.
+  generation to all distributions.
 
 * Support for filling existing arrays using `out` keyword argument. Currently
   supported in (both 32- and 64-bit outputs)
@@ -61,7 +68,7 @@ y = rnd.standard_gamma(5.5, 10000, method='zig')
 
 ## Included Pseudo Random Number Generators
 
-This modules includes a number of alternative random 
+This module includes a number of alternative random
 number generators in addition to the MT19937 that is included in NumPy. 
 The RNGs include:
 
@@ -71,23 +78,22 @@ The RNGs include:
   SSE2-aware version of the MT19937 generator that is especially fast at 
   generating doubles
 * [xoroshiro128+](http://xoroshiro.di.unimi.it/) and
-  [xorshift1024*](http://xorshift.di.unimi.it/)
+  [xorshift1024*φ](http://xorshift.di.unimi.it/)
 * [PCG64](http:w//www.pcg-random.org/)
-* ThreeFry and Philox implementationf from [Random123](https://www.deshawrsearch.com/resources_random123.html)
+* ThreeFry and Philox from [Random123](https://www.deshawrsearch.com/resources_random123.html)
 ## Differences from `numpy.random.RandomState`
 
 ### New Features
 * `standard_normal`, `normal`, `randn` and `multivariate_normal` all 
-  support an additional `method` keyword argument which can be `bm` or
-  `zig` where `bm` corresponds to the current method using the Box-Muller
-  transformation and `zig` uses the much faster (100%+) Ziggurat method.
-* `standard_exponential` and `standard_gamma` both support an additional
+  use the much faster (100%+) Ziggurat method.
+* `standard_gamma` and `gamma` both use the much faster Ziggurat method.
+* `standard_exponential` `exponential` both support an additional
   `method` keyword argument which can be `inv` or
   `zig` where `inv` corresponds to the current method using the inverse
   CDF and `zig` uses the much faster (100%+) Ziggurat method.
 * Core random number generators can produce either single precision
   (`np.float32`) or double precision (`np.float64`, the default) using
-  an the optional keyword argument `dtype`
+  the optional keyword argument `dtype`
 * Core random number generators can fill existing arrays using the
   `out` keyword argument
 
@@ -126,7 +132,7 @@ The version matched the latest version of NumPy where
 ## Documentation
 
 An occasionally updated build of the documentation is available on
-[my github pages](http://bashtage.github.io/core-prng/).
+[my GitHub pages](http://bashtage.github.io/core-prng/).
 
 ## Plans
 This module is essentially complete.  There are a few rough edges that 
@@ -139,8 +145,8 @@ need to be smoothed.
 Building requires:
 
   * Python (2.7, 3.4, 3.5, 3.6)
-  * NumPy (1.9, 1.10, 1.11, 1.12)
-  * Cython (0.22, **not** 0.23, 0.24, 0.25)
+  * NumPy (1.10, 1.11, 1.12, 1.13, 1.14)
+  * Cython (0.25+)
   * tempita (0.5+), if not provided by Cython
 
 Testing requires pytest (3.0+).
@@ -151,12 +157,12 @@ versions.
 ## Development and Testing
 
 All development has been on 64-bit Linux, and it is regularly tested on 
-Travis-CI (Linux) and Appveyor (Windows). The library is occasionally 
-tested on Linux 32-bit, OSX 10.13, Free BSD 11.1.
+Travis-CI (Linux/OSX) and Appveyor (Windows). The library is occasionally
+tested on Linux 32-bit and Free BSD 11.1.
 
 Basic tests are in place for all RNGs. The MT19937 is tested against 
 NumPy's implementation for identical results. It also passes NumPy's 
-test suite.
+test suite where still relevant.
 
 ## Installing
 
@@ -206,37 +212,37 @@ NumPy's mt19937.
 ```
 Speed-up relative to NumPy (Uniform Doubles)
 ************************************************************
-MT19937          22.9%
-PCG64           109.6%
-Philox           -6.2%
-ThreeFry        -16.6%
-Xoroshiro128    161.0%
-Xorshift1024    119.9%
+DSFMT           137.1%
+MT19937          21.0%
+PCG32           101.2%
+PCG64           110.7%
+Philox           -2.7%
+ThreeFry        -11.4%
+ThreeFry32      -62.3%
+Xoroshiro128    181.4%
+Xorshift1024    141.8%
 
 Speed-up relative to NumPy (64-bit unsigned integers)
 ************************************************************
-MT19937           6.2%
-PCG64            88.2%
-Philox          -23.0%
-ThreeFry        -26.5%
-Xoroshiro128    142.4%
-Xorshift1024    107.5%
-
-Speed-up relative to NumPy (Standard normals (Box-Muller))
-************************************************************
-MT19937          17.7%
-PCG64            35.6%
-Philox          -26.2%
-ThreeFry        -16.9%
-Xoroshiro128     57.9%
-Xorshift1024     40.9%
-
-Speed-up relative to NumPy (Standard normals (Ziggurat))
+DSFMT            24.8%
+MT19937          15.0%
+PCG32            92.6%
+PCG64            99.0%
+Philox          -20.4%
+ThreeFry        -21.7%
+ThreeFry32      -64.4%
+Xoroshiro128    164.2%
+Xorshift1024    120.8%
+
+Speed-up relative to NumPy (Standard normals)
 ************************************************************
-MT19937         107.9%
-PCG64           149.6%
-Philox           11.1%
-ThreeFry         78.8%
-Xoroshiro128    224.7%
-Xorshift1024    158.6%
-```
+DSFMT           299.4%
+MT19937         271.2%
+PCG32           364.5%
+PCG64           364.2%
+Philox          256.9%
+ThreeFry        236.0%
+ThreeFry32       97.0%
+Xoroshiro128    477.4%
+Xorshift1024    360.7%
+```

_randomgen/benchmark.py

@@ -103,16 +103,9 @@ def timer_64bit():
     run_timer(dist, command, command_numpy, SETUP, '64-bit unsigned integers')
 
 
-def timer_normal():
-    dist = 'standard_normal'
-    command = 'rg.standard_normal(1000000, method="bm")'
-    command_numpy = 'rg.standard_normal(1000000)'
-    run_timer(dist, command, command_numpy, SETUP, 'Box-Muller normals')
-
-
 def timer_normal_zig():
     dist = 'standard_normal'
-    command = 'rg.standard_normal(1000000, method="zig")'
+    command = 'rg.standard_normal(1000000)'
     command_numpy = 'rg.standard_normal(1000000)'
     run_timer(dist, command, command_numpy, SETUP,
               'Standard normals (Ziggurat)')
@@ -129,5 +122,4 @@ def timer_normal_zig():
         timer_raw()
         timer_32bit()
         timer_64bit()
-        timer_normal()
         timer_normal_zig()

_randomgen/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);

_randomgen/core_prng/distributions.pxd

@@ -32,7 +32,6 @@ cdef extern from "src/distributions/distributions.h":
         uint32_t (*next_uint32)(void *st) nogil
         double (*next_double)(void *st) nogil
         uint64_t (*next_raw)(void *st) nogil
-        binomial_t *binomial
 
     ctypedef prng prng_t
 
@@ -85,7 +84,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

_randomgen/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):

_randomgen/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,8 +70,12 @@ 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()
 
+    def __dealloc__(self):
+        free(self._binomial)
+
     def __repr__(self):
         return self.__str__() + ' at 0x{:X}'.format(id(self))
 
@@ -3160,10 +3163,64 @@ cdef class RandomGenerator:
         >>> sum(np.random.binomial(9, 0.1, 20000) == 0)/20000.
         # 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):
         """
@@ -3913,7 +3970,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