Stop passing degenerate parameters to the MKL discrete generators

CHANGELOG.md

@@ -29,8 +29,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Sped up `randint` with `array_like` bounds by generating words in cache-sized chunks and by picking the Lemire rejection test by range width [gh-173](https://github.com/IntelPython/mkl_random/pull/173)
 * Reduced the memory use and sped up integer generation with `WH`, `MCG31`, `R250` and `MRG32K3A`, which lack `viRngUniformBits` support [gh-178](https://github.com/IntelPython/mkl_random/pull/178)
 * Replaced the `mklrand_py_helper.h` shim with `cpython.bytes` cimports and dropped the no-op `__future__` import, since the project requires Python 3.10+ [gh-192](https://github.com/IntelPython/mkl_random/pull/192)
+* `negative_binomial` now rejects `p == 0`, which MKL does not support, with a `ValueError`; the documented range of `p` is `(0, 1]` [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* `poisson`, `negative_binomial`, `geometric` and `logseries` now reject `NaN` parameters with a `ValueError`, as NumPy does, instead of passing them to MKL [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 
 ### Fixed
+* Fixed `poisson` with a scalar `lam == 0` returning uninitialized memory as random values, since MKL requires `lam > 0`. The same applied to array-valued `lam` with `method="PTPE"`. Both now return zeros, as the default array path and NumPy already did [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* Fixed `negative_binomial` with `p == 1` returning uninitialized memory as random values, since MKL requires `p < 1`; it now returns zeros, as NumPy does [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* Fixed `geometric` with `p` being `NaN` returning uninitialized memory as random values [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
+* Fixed `logseries` rejecting `p == 0`, which NumPy documents as valid; it now returns ones, which its sampler already produced without a special case [gh-195](https://github.com/IntelPython/mkl_random/pull/195)
 * Fixed `multinormal_cholesky` and the other functions taking `method` routing integer-like ids such as `np.int64(0)`, `np.int64(2)` or `False` to `BoxMuller`; NumPy integers now select the method they name [gh-189](https://github.com/IntelPython/mkl_random/pull/189)
 * Fixed `multinormal_cholesky` raising `ZeroDivisionError` for an empty `mean`; it now raises `ValueError` [gh-189](https://github.com/IntelPython/mkl_random/pull/189)
 * Fixed an out-of-range integer `brng` indexing `brng_list` past its end, which returned uninitialized memory as random values; it now warns and falls back to `MT19937` [gh-177](https://github.com/IntelPython/mkl_random/pull/177)

mkl_random/mklrand.pyx

@@ -43,6 +43,7 @@ from cpython.bytes cimport (
     PyBytes_Check,
     PyBytes_FromStringAndSize,
 )
+from libc.math cimport isnan
 from libc.stdint cimport int64_t
 from libc.string cimport memcpy, memset
 
@@ -5871,14 +5872,14 @@ cdef class _MKLRandomState:
 
         Samples are drawn from a negative binomial distribution with specified
         parameters, `n` trials and `p` probability of success where `n` is an
-        integer > 0 and `p` is in the interval [0, 1].
+        integer > 0 and `p` is in the interval (0, 1].
 
         Parameters
         ----------
         n : int
             Parameter, > 0.
         p : float
-            Parameter, >= 0 and <=1.
+            Parameter, > 0 and <= 1.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  Default is None, in which case a
@@ -5949,12 +5950,10 @@ cdef class _MKLRandomState:
             fp = PyFloat_AsDouble(p)
             fn = PyFloat_AsDouble(n)
 
-            if fn <= 0:
-                raise ValueError("n <= 0")
-            if fp < 0:
-                raise ValueError("p < 0")
-            elif fp > 1:
-                raise ValueError("p > 1")
+            if isnan(fn) or fn <= 0:
+                raise ValueError("n <= 0 or n is NaN")
+            if isnan(fp) or fp <= 0 or fp > 1:
+                raise ValueError("p <= 0, p > 1 or p is NaN")
             return vec_discdd_array_sc(
                 self.internal_state,
                 irk_negbinomial_vec,
@@ -5964,12 +5963,14 @@ cdef class _MKLRandomState:
                 self.lock
             )
 
-        if np.any(np.less_equal(n, 0)):
-            raise ValueError("n <= 0")
-        if np.any(np.less(p, 0)):
-            raise ValueError("p < 0")
-        if np.any(np.greater(p, 1)):
-            raise ValueError("p > 1")
+        if np.any(np.isnan(on)) or np.any(np.less_equal(on, 0)):
+            raise ValueError("n <= 0 or n is NaN")
+        if (
+            np.any(np.isnan(op))
+            or np.any(np.less_equal(op, 0))
+            or np.any(np.greater(op, 1))
+        ):
+            raise ValueError("p <= 0, p > 1 or p is NaN")
         return vec_discdd_array(
             self.internal_state, irk_negbinomial_vec, size, on, op, self.lock
         )
@@ -6052,9 +6053,9 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(olam) == 0:
             flam = PyFloat_AsDouble(lam)
 
-            if lam < 0:
-                raise ValueError("lam < 0")
-            if lam > self._poisson_lam_max:
+            if isnan(flam) or flam < 0:
+                raise ValueError("lam < 0 or lam is NaN")
+            if flam > self._poisson_lam_max:
                 raise ValueError("lam value too large")
             method = choose_method(
                 method, [POISNORM, PTPE], _method_alias_dict_poisson
@@ -6076,8 +6077,8 @@ cdef class _MKLRandomState:
                     self.lock
                 )
 
-        if np.any(np.less(olam, 0)):
-            raise ValueError("lam < 0")
+        if np.any(np.isnan(olam)) or np.any(np.less(olam, 0)):
+            raise ValueError("lam < 0 or lam is NaN")
         if np.any(np.greater(olam, self._poisson_lam_max)):
             raise ValueError("lam value too large.")
         method = choose_method(
@@ -6251,18 +6252,18 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(op) == 0:
             fp = PyFloat_AsDouble(p)
 
-            if fp <= 0.0:
-                raise ValueError("p <= 0.0")
-            if fp > 1.0:
-                raise ValueError("p > 1.0")
+            if isnan(fp) or fp <= 0.0 or fp > 1.0:
+                raise ValueError("p <= 0, p > 1 or p is NaN")
             return vec_discd_array_sc(
                 self.internal_state, irk_geometric_vec, size, fp, self.lock
             )
 
-        if np.any(np.less_equal(op, 0.0)):
-            raise ValueError("p <= 0.0")
-        if np.any(np.greater(op, 1.0)):
-            raise ValueError("p > 1.0")
+        if (
+            np.any(np.isnan(op))
+            or np.any(np.less_equal(op, 0.0))
+            or np.any(np.greater(op, 1.0))
+        ):
+            raise ValueError("p <= 0, p > 1 or p is NaN")
         return vec_discd_array(
             self.internal_state, irk_geometric_vec, size, op, self.lock
         )
@@ -6435,14 +6436,12 @@ cdef class _MKLRandomState:
         Draw samples from a logarithmic series distribution.
 
         Samples are drawn from a log series distribution with specified
-        shape parameter, 0 < ``p`` < 1.
+        shape parameter, 0 <= ``p`` < 1.
 
         Parameters
         ----------
-        loc : float
-
-        scale : float > 0.
-
+        p : float or array_like of floats
+            Shape parameter, 0 <= ``p`` < 1.
         size : int or tuple of ints, optional
             Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
             ``m * n * k`` samples are drawn.  Default is None, in which case a
@@ -6451,7 +6450,7 @@ cdef class _MKLRandomState:
         Returns
         -------
         samples : ndarray or scalar
-                  where the values are all integers in  [0, n].
+                  where the values are all positive integers.
 
         See Also
         --------
@@ -6513,18 +6512,18 @@ cdef class _MKLRandomState:
         if cnp.PyArray_NDIM(op) == 0:
             fp = PyFloat_AsDouble(p)
 
-            if fp <= 0.0:
-                raise ValueError("p <= 0.0")
-            if fp >= 1.0:
-                raise ValueError("p >= 1.0")
+            if isnan(fp) or fp < 0.0 or fp >= 1.0:
+                raise ValueError("p < 0, p >= 1 or p is NaN")
             return vec_discd_array_sc(
                 self.internal_state, irk_logseries_vec, size, fp, self.lock
             )
 
-        if np.any(np.less_equal(op, 0.0)):
-            raise ValueError("p <= 0.0")
-        if np.any(np.greater_equal(op, 1.0)):
-            raise ValueError("p >= 1.0")
+        if (
+            np.any(np.isnan(op))
+            or np.any(np.less(op, 0.0))
+            or np.any(np.greater_equal(op, 1.0))
+        ):
+            raise ValueError("p < 0, p >= 1 or p is NaN")
         return vec_discd_array(
             self.internal_state, irk_logseries_vec, size, op, self.lock
         )

mkl_random/src/mkl_distributions.cpp

@@ -1348,6 +1348,14 @@ void irk_negbinomial_vec(irk_state *state,
     if (len < 1)
         return;
 
+    /* MKL requires p < 1, while p == 1 means every trial succeeds */
+    if (p == 1.0) {
+        npy_intp i;
+        for (i = 0; i < len; ++i)
+            res[i] = 0;
+        return;
+    }
+
     while (len > MKL_INT_MAX) {
         err = viRngNegbinomial(VSL_RNG_METHOD_NEGBINOMIAL_NBAR, state->stream,
                                MKL_INT_MAX, res, a, p);
@@ -1399,6 +1407,14 @@ void irk_poisson_vec_PTPE(irk_state *state,
     if (len < 1)
         return;
 
+    /* MKL requires lambda > 0, while lambda == 0 always gives zero events */
+    if (lambda == 0.0) {
+        npy_intp i;
+        for (i = 0; i < len; ++i)
+            res[i] = 0;
+        return;
+    }
+
     while (len > MKL_INT_MAX) {
         err = viRngPoisson(VSL_RNG_METHOD_POISSON_PTPE, state->stream,
                            MKL_INT_MAX, res, lambda);
@@ -1423,6 +1439,14 @@ void irk_poisson_vec_POISNORM(irk_state *state,
     if (len < 1)
         return;
 
+    /* MKL requires lambda > 0, while lambda == 0 always gives zero events */
+    if (lambda == 0.0) {
+        npy_intp i;
+        for (i = 0; i < len; ++i)
+            res[i] = 0;
+        return;
+    }
+
     while (len > MKL_INT_MAX) {
         err = viRngPoisson(VSL_RNG_METHOD_POISSON_POISNORM, state->stream,
                            MKL_INT_MAX, res, lambda);

mkl_random/tests/test_random.py

@@ -926,6 +926,12 @@ def test_randomdsit_geometric(randomdist):
     np.testing.assert_array_equal(actual, desired)
 
 
+def test_geometric_exceptions():
+    for p in (0.0, -0.5, 1.5, np.nan):
+        assert_raises(ValueError, rnd.geometric, p)
+        assert_raises(ValueError, rnd.geometric, [p] * 5)
+
+
 def test_randomdist_gumbel(randomdist):
     rnd.seed(randomdist.seed, randomdist.brng)
     actual = rnd.gumbel(loc=0.123456789, scale=2.0, size=(3, 2))
@@ -1025,6 +1031,21 @@ def test_randomdist_logseries(randomdist):
     np.testing.assert_array_equal(actual, desired)
 
 
+def test_logseries_exceptions():
+    for p in (-0.5, 1.0, np.nan):
+        assert_raises(ValueError, rnd.logseries, p)
+        assert_raises(ValueError, rnd.logseries, [p] * 5)
+
+
+def test_logseries_zero_p():
+    # NumPy documents p in [0, 1) and the sampler needs no special case:
+    # with p == 0 the acceptance test V >= p always holds
+    ones = np.ones(5, dtype="int32")
+    assert rnd.logseries(0.0) == 1
+    np.testing.assert_array_equal(rnd.logseries(0.0, 5), ones, strict=True)
+    np.testing.assert_array_equal(rnd.logseries(np.zeros(5)), ones, strict=True)
+
+
 def test_randomdist_multinomial(randomdist):
     rs = rnd.MKLRandomState(randomdist.seed, brng=randomdist.brng)
     actual = rs.multinomial(20, [1 / 6.0] * 6, size=(3, 2))
@@ -1169,6 +1190,27 @@ def test_randomdist_negative_binomial(randomdist):
     np.testing.assert_array_equal(actual, desired)
 
 
+def test_negative_binomial_p_is_one():
+    # MKL requires p < 1, while p == 1 means no trial ever fails
+    zeros = np.zeros(5, dtype="int32")
+    assert rnd.negative_binomial(5, 1.0) == 0
+    np.testing.assert_array_equal(
+        rnd.negative_binomial(5, 1.0, 5), zeros, strict=True
+    )
+    np.testing.assert_array_equal(
+        rnd.negative_binomial(np.full(5, 5.0), 1.0), zeros, strict=True
+    )
+
+
+def test_negative_binomial_exceptions():
+    for p in (0.0, -0.5, 1.5, np.nan):
+        assert_raises(ValueError, rnd.negative_binomial, 5, p)
+        assert_raises(ValueError, rnd.negative_binomial, 5, [p] * 5)
+    for n in (0, -1, np.nan):
+        assert_raises(ValueError, rnd.negative_binomial, n, 0.5)
+        assert_raises(ValueError, rnd.negative_binomial, [n] * 5, 0.5)
+
+
 def test_randomdist_noncentral_chisquare(randomdist):
     rnd.seed(randomdist.seed, brng=randomdist.brng)
     actual = rnd.noncentral_chisquare(df=5, nonc=5, size=(3, 2))
@@ -1282,6 +1324,26 @@ def test_poisson_exceptions():
     assert_raises(ValueError, rnd.poisson, [lamneg] * 10)
     assert_raises(ValueError, rnd.poisson, lambig)
     assert_raises(ValueError, rnd.poisson, [lambig] * 10)
+    assert_raises(ValueError, rnd.poisson, np.nan)
+    assert_raises(ValueError, rnd.poisson, [np.nan] * 10)
+
+
+@pytest.mark.parametrize("method", ["POISNORM", "PTPE"])
+def test_poisson_zero_lam(method):
+    # MKL requires lam > 0, while lam == 0 always gives zero events
+    zeros = np.zeros(5, dtype="int32")
+    assert rnd.poisson(0.0, method=method) == 0
+    np.testing.assert_array_equal(
+        rnd.poisson(0.0, 5, method=method), zeros, strict=True
+    )
+    np.testing.assert_array_equal(
+        rnd.poisson(np.zeros(5), method=method), zeros, strict=True
+    )
+    # A lam too large to draw a zero, so that a call rejected as a whole,
+    # which would leave the uninitialized buffer alone, fails here
+    mixed = rnd.poisson(np.array([0.0, 1000.0, 0.0]), method=method)
+    np.testing.assert_array_equal(mixed[::2], zeros[:2], strict=True)
+    assert mixed[1] > 0
 
 
 def test_randomdist_power(randomdist):