diff --git a/CHANGELOG.md b/CHANGELOG.md index b05a04e..155597b 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/mkl_random/mklrand.pyx b/mkl_random/mklrand.pyx index 7f55da6..f8d09c1 100644 --- a/mkl_random/mklrand.pyx +++ b/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 ) diff --git a/mkl_random/src/mkl_distributions.cpp b/mkl_random/src/mkl_distributions.cpp index 81c6f0e..52fcdc7 100644 --- a/mkl_random/src/mkl_distributions.cpp +++ b/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); diff --git a/mkl_random/tests/test_random.py b/mkl_random/tests/test_random.py index 13eb55b..be50a25 100644 --- a/mkl_random/tests/test_random.py +++ b/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):