Repository navigation
Stop passing degenerate parameters to the MKL discrete generators - #195
Open
vlad-perevezentsev wants to merge 10 commits into
Open
vlad-perevezentsev wants to merge 10 commits into
vlad-perevezentsev wants to merge 10 commits into
Conversation
viRngPoisson requires lam > 0, so lam == 0 made it fail with "Parameter 5 was incorrect on entry to viRngPoisson" and leave the output buffer untouched, which was then handed back as random values. The assert on the MKL status code cannot catch that: meson.build defines NDEBUG, so every assert in mkl_distributions.cpp is a no-op in shipped builds. The scalar path hit this directly, and so did array-valued lam with method="PTPE", which calls the scalar kernel per element; only the default POISNORM array path was correct, since viRngPoissonV accepts a zero lambda. The distribution is degenerate at lam == 0, so the kernels now fill zeros without touching the stream, which is what the array path and NumPy already returned. Streams for a positive lam are unaffected. NaN lam passed validation as well, because `lam < 0` is false for NaN, and MKL then returned arbitrary values. It now raises ValueError with NumPy's wording. The isnan check runs first, as the comparison ufuncs warn on NaN input.
viRngNegbinomial requires 0 < p < 1, and both endpoints reached it. p == 1 only reported "INTERNAL ERROR: Condition -6", and p == 0 passed a check that tested `p < 0`, so in both cases the untouched output buffer was handed back as random values, on the scalar and on the array path alike. The assert on the MKL status code is a no-op, as meson.build defines NDEBUG. p == 1 means no trial ever fails, so the number of failures is zero and the kernel now fills zeros without touching the stream, as NumPy returns. p == 0 has no finite answer, NumPy's own Generator rejects it, and it now raises ValueError, which makes the documented range of p (0, 1]. NaN was unchecked for both n and p, so it is rejected together with them.
Neither validated NaN: `p <= 0.0`, `p > 1.0` and `p >= 1.0` are all false for NaN, so it reached the kernels. geometric then wrote nothing and handed back uninitialized memory as random values, since viRngGeometric is only called for 0 < p < 1; logseries silently returned ones. NumPy raises ValueError for both. The checks are merged into one condition per path, with isnan first, as the comparison ufuncs warn on NaN input.
vlad-perevezentsev
requested review from
antonwolfy,
jharlow-intel,
ndgrigorian and
xaleryb
as code owners
October 16, 2026 12:21
vlad-perevezentsev
force-pushed
the
fix/degenerate-params-uninitialized-output
branch
from
October 16, 2026 12:22
83356fb to
ed246ad
Compare
antonwolfy
approved these changes
Oct 16, 2026
antonwolfy
left a comment
antonwolfy
left a comment
Collaborator
There was a problem hiding this comment.
Thanks, this is a nice catch! The scalar and array checks now agree with each other and with NumPy at the boundaries, and the C guards follow the existing p == 1 handling in irk_geometric_vec, so valid-parameter streams are untouched. A couple of small things inline.
- Not for this PR, but while you're in this area:
binomialwith arraypstill lets NaN through toviRngBinomial, since only the scalar path checksisnan.noncentral_chisquare(0.5, 0)callsviRngPoisson(PTPE, ..., lambda=0)directly inirk_noncentral_chisquare_vec, bypassing the new guard. Worth a follow-up issue?
Per review: with lam = 5 for the middle element, the assertion on the zero elements would also hold if viRngPoissonV had rejected the whole call and left the zero-filled output buffer untouched, which is exactly the failure mode this test is about. A lam of 1000 cannot produce a zero, so the element now proves the call went through. That viRngPoissonV accepts a zero lambda was checked directly, not assumed: poisson([0, 1000, 0]) returns [0, 1005, 0] with nothing on stderr, while the scalar and PTPE paths printed "Parameter 5 was incorrect on entry to viRngPoisson" before this branch.
Per review: gh-193 is the unrelated set_state buffer-length issue. The entries were written before the PR existed and kept the guessed number.
antonwolfy
reviewed
Oct 16, 2026
Per review: both RandomState and Generator document p in [0, 1) and return 1 there, so rejecting it was a visible difference when mkl_random is patched into numpy.random. The sampler needs no special case, as noted in the review: with p == 0, r = log(1 - p) is zero and the acceptance test V >= p always holds, so irk_logseries_vec already fills ones. Checked: logseries(0.0, 5), logseries(np.zeros(5)) and a mixed array all return ones where p is zero, matching NumPy. The docstring prose said 0 < p < 1 and now matches the code. negative_binomial keeps rejecting p == 0: MKL has no such distribution there and the number of failures has no finite answer.
Same point as the review on gh-196: assert_array_equal ignores dtypes by default, so the int32 expectations were not actually checked and an int64 or float result would have passed. strict=True also pins the shape and needs NumPy 1.24, below the 1.26.4 minimum.
antonwolfy
reviewed
Oct 16, 2026
Per review: vec_Poisson_array allocates the output with PyArray_SimpleNew, and the scalar path with np.empty, so a rejected call leaves whatever was in memory rather than zeros. The check itself is unaffected.
Per review, older than this PR but in the lines being touched: the Parameters block described loc and scale, which the function does not take, and never mentioned p, while Returns claimed the values lie in [0, n]. The support is the positive integers. gh-165 rewrites this block as well, so that PR keeps the final wording; the range there needs to read [0, 1) after this one.
antonwolfy
approved these changes
Oct 16, 2026
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
MKL's discrete generators reject degenerate parameter values that NumPy accepts.
mkl_randompassed them through and only checked the status code with anassert, whichmeson.builddisables via-DNDEBUG, so theuntouched output buffer was handed back as random values.
Fixed
poissonwith a scalarlam == 0, and with array-valuedlamandmethod="PTPE", returned uninitialized memory; it now returns zeros, as the default array path and NumPy already did.negative_binomialwithp == 1returned uninitialized memory; it now returns zeros, as NumPy does.geometricwithpbeingNaNreturned uninitialized memory.Changed
negative_binomialnow rejectsp == 0, which MKL does not support, witha
ValueError; the documented range ofpbecomes(0, 1].poisson,negative_binomial,geometricandlogseriesnow rejectNaNparameters with aValueError, as NumPy does.Notes
MKLRandomStatein the API reference and align docstrings with NumPy #165.geometricis shifted by one relative to NumPy for everyp(MKL counts failures, NumPy counts trials), which is a breaking change and gets its own PR.