Seedable RNG, benchmark suite, and discrete CDF optimisation - #3
Open
ghosteau wants to merge 15 commits into
Open
Seedable RNG, benchmark suite, and discrete CDF optimisation#3ghosteau wants to merge 15 commits into
ghosteau wants to merge 15 commits into
Conversation
The C++ RNG tests compared sample mean/variance against theory using round-number absolute tolerances. Several landed near 2 sigma of the estimator's own standard error, so they failed at roughly the rate a 2 sigma bound fails -- ~7.7% of full-suite runs, measured. Audited all 21 RNG assertions against their standard errors and resized the eight that sat below 7 sigma. Beyond the three cases named in the issue, this also covers discrete_uniform's mean check (2.93 sigma, on the line directly above the flagged variance check and sharing its literal), uniform's variance check (3.73), normal's variance check (3.95), and both negative_binomial checks (~4.3). Each site now carries the standard error and resulting sigma in a comment so the values are not "tidied" back into round numbers later. Measured on MSVC/Windows Release, 500 runs per test plus 100 full ctest sweeps: 0 failures, against a 26/2000 baseline on the same machine. Also drops --repeat until-pass:3 from python-distro.yml, which was the standing workaround. Closes #2 Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every *_sample() function declared its own engine inline:
thread_local std::mt19937 rng{std::random_device{}()};
so the stream was unreachable from outside and reseeded from OS entropy on
every process start. That made sampling results irreproducible for users and
forced the RNG tests to be sized around sampling noise.
Hoists the engine into one shared accessor (src/math/rng.cpp) and rewires all
twelve samplers to draw from it. Adds fastdist::math::seed_rng /
seed_rng_from_entropy, exposed to Python as fastdist.seed() and
fastdist.seed_from_entropy().
Reproducibility is a feature in its own right for a distributions library --
the same argument numpy.random.seed exists for -- and the tests are the first
consumer.
Two limits are documented at the API rather than left to be discovered:
* seeding is per-thread, since the engine is thread_local
* std::mt19937 is specified bit-for-bit by the standard but the distribution
adaptors layered on it are not, so a seed reproduces a run on one platform
and toolchain, not across libstdc++, libc++ and MSVC
With the streams pinned, every RNG block now seeds before sampling, so each is
deterministic: it either always passes or always fails on a given toolchain,
never intermittently. That in turn allows tolerances to come down from ~7-70
sigma to a uniform ~5 sigma. Sensitivity to a systematic bias improves by up
to 12x at the loosest sites (beta variance was 71 sigma, poisson variance 42,
binomial mean 52).
The audit also picked up gamma and chi_square, whose sampling blocks a header
comment had hidden from the earlier pass; both were already above 7 sigma.
Adds tests/cpp/test_rng.cpp and tests/python/test_rng.py covering the contract
itself: same seed replays exactly, distinct seeds diverge, re-seeding rewinds,
and entropy reseeding escapes the fixed stream. The distinct-seeds case guards
against seed_rng ignoring its argument, which would leave every seeded test in
the suite passing while asserting nothing.
test_package_api.py is updated for the widened public API.
Verified on MSVC/Windows Release: 200 full ctest sweeps, 0 failures.
Also reformats tests/cpp/test_app.cpp, which had drifted from .clang-format;
pre-commit corrects it on touch.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The library claims to be fast but had no way to demonstrate it, and no way to
notice if a change made it slower. Adds benchmarks/ and BENCHMARKS.md.
benchmarks/run.py measures fastdist against SciPy, saving each run as JSON
under benchmarks/results/ tagged with version, commit, CPU and toolchain.
compare.py diffs two reports and exits non-zero on a >5% regression, so it can
gate a change. table.py renders a report as Markdown, so entries in the log are
generated from recorded data rather than transcribed.
Methodology choices worth noting:
* The minimum round is reported, not the mean. Noise on a shared machine can
only add time, so the minimum estimates the work; the mean would estimate
how busy the machine was. The min-to-median gap is recorded separately as
noise_pct so a noisy run is visible rather than silently averaged in.
* Every case with a baseline checks both implementations agree numerically
before either is timed. A speedup on a wrong answer is not a speedup, and
the recorded diffs (~1e-16) are what make the numbers quotable.
* Cases are grouped by how fair the comparison is. Only "batch" is a real
throughput comparison; "scalar" measures per-call overhead and flatters
whichever library has the thinner binding layer, so the log says not to
quote it.
The v0.1.0 baseline is recorded, including the results that do not flatter the
library:
* batch PDF/CDF beats SciPy by 3-15x, widest on the cheapest distributions
(uniform 10-15x, bernoulli 11-27x) where SciPy's per-call overhead
dominates, narrowing as arrays grow and real arithmetic takes over
* poisson_cdf is 6.6x *slower* than SciPy. poisson_cdf_scalar sums the PMF
from 0 to k and each term recomputes log(lambda), an lgamma and an exp --
roughly twenty transcendental calls per element where a recurrence needs
none. A real defect, and the clearest optimisation target in the library.
* sampling is 30-100x slower than numpy, because fastdist crosses the
Python/C++ boundary once per variate while numpy fills an array per call.
Structural, not tuning: it needs a batch sampling entry point that does not
exist yet.
scipy is a benchmark-only dependency and is not required to build or use the
library.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two classes of avoidable work, both found by the benchmark suite rather than by
inspection.
Discrete CDFs summed PMF terms
------------------------------
poisson_cdf_scalar, binomial_cdf_scalar and negative_binomial_cdf_scalar each
summed P(0)..P(k) by calling their PMF once per term, and every one of those
calls re-derived the whole term from scratch -- a log, an lgamma and an exp for
Poisson; three lgammas, two logs and an exp for Binomial; two tgammas and two
pows for Negative Binomial.
Consecutive terms are related by a ratio, so the sums now run as recurrences
and need one exp in total:
Poisson P(i) = P(i-1) * lambda / i
Binomial P(k) = P(k-1) * ((n-k+1)/k) * (p/(1-p))
NegativeBinomial P(i) = P(i-1) * ((i+r-1)/i) * (1-p)
Each recurrence starts from a first term that underflows for extreme parameters
(exp(-lambda), (1-p)^n, p^r), which would collapse the sum to zero where the
true CDF is O(1). Each therefore keeps the original per-term evaluation as a
fallback, guarded on the first term's magnitude in log space. Binomial also
special-cases p == 1, whose ratio would divide by zero.
Measured: poisson_cdf over 100k elements goes from 43.47ms to 1.32ms, a 33x
improvement, and from 6.6x slower than SciPy to 4.9x faster. Agreement with
SciPy improves slightly too (3.3e-16 to 2.2e-16) -- fewer operations, less
accumulated rounding.
Fixes a documented overflow while there
---------------------------------------
negative_binomial_pmf_scalar formed C(k+r-1, k) from raw tgammas, which
overflow a double once k+r-1 > 170: inf at k=170, nan beyond, even though the
true PMF at r=3, k=200 is 1.6e-57. Evaluating in log space removes the ceiling
and drops two pows from the hot path. Eight tests in test_negative_binomial.py
were marked xfail(strict) against this defect; they now pass, so the markers
are removed and the cases stay as regression tests.
Batch paths recomputed loop-invariant work
------------------------------------------
Every *_batch function called its scalar counterpart per element, so parameter
validation and any term depending only on the parameters ran once per element
instead of once per array.
Validation is now hoisted (invalid parameters fill the output with NaN and
return), along with the genuinely invariant terms. Where a formula is shared,
scalar and batch both route through one inlined core rather than duplicating
it, and the arithmetic is arranged exactly as before so results stay
bit-identical.
Measured over 100k elements:
normal_logpdf 393.5us -> 77.1us (5.1x) log(sigma) per element
uniform_pdf 129.0us -> 92.7us (1.39x) the density is constant
normal_cdf 711.4us -> 529.3us (1.34x)
uniform_cdf 131.2us -> 102.1us (1.29x)
poisson_pmf 4.72ms -> 4.01ms (1.18x) log(lambda) per element
normal_pdf 462.8us -> 414.5us (1.12x)
exponential_pdf 345.0us -> 308.6us (1.12x)
Applied to the distributions the suite covers. The same pattern is worth
extending to the remaining ones; it is mechanical, but unmeasured work is not
worth claiming.
Also adds the parameter validation beta_sample was missing -- every other
sampler has it, and std::gamma_distribution is undefined rather than
error-valued for a non-positive shape.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Full run at 73ba7d6 against the v0.1.0 baseline on the same machine in the same session: 29 cases improved, none regressed. poisson_cdf -97%, normal_logpdf -80%, and the remaining continuous batch paths 7-29%. The reading notes the two gaps that did not move -- poisson_pmf, whose per-element lgamma is not loop-invariant, and sampling, which needs a batch entry point rather than tuning. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
README
------
The build section described `python3 setup.py bdist_wheel` against a
hand-configured CMake tree. `pip install .` has done the whole job since
pyproject.toml started declaring the build dependencies, so that is the
documented path now, with `python -m build` for a wheel.
Adds the direct CMake instructions for C++ work, including the multi-config
note for the Visual Studio generator and -- because this is the first thing a
fresh checkout hits -- what to do when CMake cannot find pybind11, which it
resolves from the active interpreter rather than a submodule.
Corrects the supported Python range in the build_all flag docs: it said 3.12
through 3.14, while build_all.ps1, setup.py's python_requires and the
classifiers all say 3.10.
New sections for Reproducible Sampling (with both caveats stated up front: the
seed is per-thread, and it reproduces a run on one platform rather than across
toolchains), Benchmarks, and Testing.
The Testing section documents the tolerance policy for sampling tests -- seed
first, then size the tolerance from the estimator's standard error rather than
picking a round number -- and says why, since round numbers are exactly how the
suite acquired a 7.7% flake rate.
The long-term plans list is grouped rather than flat, and the entries that came
out of the benchmark run are specific about what and why: batch sampling entry
points as the largest gap, poisson_pmf's non-hoistable lgamma, and extending
the invariant hoisting to the distributions not yet covered. Completed items
are marked rather than silently dropped.
Import guards
-------------
All thirteen distribution modules wrapped their imports in
except ImportError:
raise ImportError("Internal Error: C++ core (_fastdist) not found...")
which discards the original exception and asserts a cause it has not checked.
Several modules also imported `config` inside the guard, and config imports
numpy -- so a missing numpy reported itself as a missing C++ core. That is not
hypothetical; it cost time during this work.
The guard now covers only the extension import, chains with `from exc` so the
real cause survives, and says what to do about it. `config` is imported outside
the guard, where a numpy failure raises its own accurate error.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
beta_cdf_scalar disagreed with the truth everywhere, not just at the extremes:
0.0015 against 0.1143 for Beta(2,5) at x=0.1, 1.128 for Beta(0.5,0.5) at
x=0.25, and -147 for Beta(0.01,0.01) at x=0.5. A CDF above 1 or below 0 is not
a rounding artifact.
Two independent defects in the series:
* the coefficient ratio was inverted -- (a+n-1)/(a+b+n-1) where the
hypergeometric series for 2F1(a+b, 1; a+1; x) needs (a+b+n-1)/(a+n)
* it normalised by Gamma(a+1) instead of B(a,b)
Replaced with the standard modified-Lentz continued fraction (Numerical
Recipes 6.4) rather than a corrected series. The series converges slowly when
a+b is large -- a corrected version still lost 2e-9 at Beta(100,100) against
MAX_ITER=100 -- while the continued fraction converges quickly across the whole
parameter range. It is also the algorithm the suite's own conftest reference
already implements, so the test and the backend now agree by construction
rather than by coincidence.
The leading factor is formed in log space: B(a,b) overflows for large
parameters and x^a underflows for large a, even where their combination is an
ordinary number. The reflection I_x(a,b) = 1 - I_{1-x}(b,a) subtracts nearby
quantities in the upper tail, so the result is clamped to [0,1] -- a CDF
reporting 1 + 1e-16 breaks callers that treat it as a probability.
Validated against scipy.special.betainc over a 1200-point grid spanning a and b
from 0.01 to 500 and x from 1e-9 to 1-1e-9: worst absolute error 1.1e-12, and
nothing outside [0,1].
Nine tests were marked xfail(strict) against this defect, including the bounded
and monotonic property checks for alpha=beta=0.5. They now pass, so the markers
are removed and the cases stay as regression tests, joined by a direct
comparison against the conftest reference and by the (0.01,0.01) and (100,100)
parameter pairs that the old code failed hardest on.
An audit of every CDF in the library over 50 parameter sets and their full
support now reports no value above 1, below 0, non-monotonic, or NaN.
Also makes test_version.py skip rather than fail when fastdist is not installed.
It compares installed distribution metadata against the compiled constant,
which says nothing about the code when there is no installed distribution to
read -- it failed for anyone running the suite from a source checkout, which
trains people to ignore failures. The message now also names the stale
egg-info case, which reports whatever version it was last generated at.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two issues with the seeding as first written. The parameter was uint32_t, so seed(-1) and seed(2**40) were rejected with pybind11's raw argument-type TypeError. Seeding from hash(x) or a signed counter is ordinary usage and produces both. The binding now takes a signed 64-bit value and reinterprets it, so the whole int64 range is accepted. The value also went straight to mt19937::seed(uint32_t), which derives all 624 state words from a single word by a simple recurrence. std::seed_seq exists to do that mixing properly, and passing both halves of the value through it also lets the full 64 bits contribute rather than only the low 32. Verified rather than assumed. Over 20,000 consecutive seeds, taking the first draw from each: lag-1 autocorrelation -0.0033 (0.46 SE), KS against U(0,1) p=0.09, and a mean absolute gap between adjacent seeds' first draws of 0.33337 against the 1/3 expected of independent uniforms. Within-stream quality is unchanged: 200k draws pass KS for both uniform and normal. Adds tests for the widened range and a guard on adjacent-seed decorrelation, which would catch a regression to naive single-word seeding. Re-verified determinism after the change: 150 full ctest sweeps, 0 failures. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
gamma_cdf_scalar and chi_square_cdf_scalar returned probabilities above 1.0 --
Gamma(1.5, 1.0).cdf_scalar(2.5) gave 1.000498004 and ChiSquare(3).cdf(7.5) gave
1.000258 -- with absolute errors reaching 0.26. Two independent causes.
The Lentz coefficient was written
const double an = -i * (i - a);
with `i` an unsigned loop index. The unary minus applies to the unsigned value
and wraps to 2^32 - i, so at i = 1, a = 1.5 the coefficient was -2147483647.5
instead of 0.5 and the whole continued fraction was garbage. Converting to
double before the negation fixes it. A replica of the loop with C++ unsigned
semantics reproduces the old output bit-for-bit on every tested point, which is
what confirmed the diagnosis rather than merely fitting it.
Separately MAX_ITER was 100. Near x = alpha the gamma series needs roughly
sqrt(2 * alpha * ln(1/EPS)) terms -- 227 at alpha = 1000, 683 at alpha = 10000
-- and the loop simply stopped early and returned the truncated sum. That cost
9e-4 at alpha = 1000 and 0.16 at alpha = 10000, silently. Raised to 1000, which
covers alpha to roughly 20000; the loops exit on convergence, so ordinary calls
pay nothing. The header records where the ceiling still bites and that a larger
one is not the answer beyond that point.
Validated against scipy over 784 gamma points (alpha 0.01 to 15000), 132
chi-square points (k 0.5 to 10000) and 343 beta points: worst absolute error
3.1e-12, nothing outside [0, 1].
This was the last of the CDFs that could exceed 1.0. Thirty tests were marked
xfail(strict) against these defects, including the bounded and monotonic
property checks; they now pass, so the markers are removed and the cases stay
as regression tests, extended to the large-shape regime that the MAX_ITER
ceiling was hiding.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
All five were documented by strict-xfail tests; the markers are now removed and
the cases stay as regression tests.
Validation checks outside their None guards (Beta, Binomial)
------------------------------------------------------------
_validate_params takes each parameter as optional, but the range check sat
outside the `is not None` guard:
if alpha is not None:
if not isinstance(alpha, (int, float)):
raise TypeError(...)
if alpha <= 0: # <- outside the guard
raise ValueError(...)
So validating only `beta` compared None <= 0 and raised TypeError for every
assignment, valid or not, making the property unusable. Binomial had the same
shape on `n`, which made its `p` setter unusable. Gamma and NegativeBinomial
already had it right and were the reference for the fix.
Recursive setter (DiscreteUniform)
----------------------------------
The b setter assigned `self.b = value` rather than `self._b`, re-entering
itself until RecursionError.
Type-changing setter (DiscreteUniform)
--------------------------------------
The a setter stored float(value) although a is an integer parameter that
__init__ stores via int(), so assigning through the setter silently changed the
attribute's type.
Setters that skipped the cross-bound check (Uniform, DiscreteUniform)
----------------------------------------------------------------------
Each setter validated only the bound being assigned, so the a < b check -- which
only runs when both are present -- never fired. Uniform(1.0, 3.0) could be
driven to a = 10.0, b = -10.0, a state the constructor rejects outright, after
which pdf, cdf, mean, variance and sample all returned nan instead of raising.
The setters' own docstrings promised this re-validation. Every setter now passes
the opposite bound alongside the new value, so an assignment that would invalid-
ate the pair is rejected before anything is mutated.
Suite is 1779 passing with 2 xfails remaining, down from 42 at the start of this
branch.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Both were annotated Union[Real, Sequence[Real]] while accepting only one of the two. The tests documenting them said "either the annotation or the implementation is wrong" without picking; these are the picks. law_of_total_probability forwarded its arguments straight to a vector-only binding, so a scalar failed inside pybind11 with an argument-type error. A scalar is the one-element partition P(B) = P(B|A) P(A) -- well defined, already advertised by the signature, and two lines to support -- so scalars are now promoted. Mismatched sequence lengths now raise a named ValueError instead of whatever the binding did with them. sigmoid went the other way. Its body calls float() on the input, and the array path already exists as sigmoid_cpu; that scalar/batch split is the same one the distribution classes use, and returning an ndarray from a function annotated -> float would be worse than not accepting one. So the annotation was the wrong half. sigmoid is now Real-only and rejects a sequence with a message naming sigmoid_cpu, rather than letting float() fail with numpy's "only 0-dimensional arrays can be converted" -- which says nothing about what to do instead. The suite now has no xfails at all: 1782 passing, from 42 documented known bugs at the start of this branch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The suite now times the *_cuda entry points against this library's own *_cpu
path rather than against SciPy, so the reported number is the speedup the GPU
backend buys over the CPU one -- the decision a caller actually faces. GPU
timings include the host-to-device copy and the copy back, because a caller
cannot avoid those. measure() checks the two paths agree numerically before
timing either, which is the part worth having: a kernel that is fast and wrong
is the failure mode.
The cases are generated behind a hasattr(core, "normal_pdf_cuda") guard, so they
are skipped entirely on a CPU-only build and the existing suite is unaffected.
Honest caveat: this code has not been executed. The CUDA backend does not build
on this machine, and I did not want to modify the toolchain to force it.
Every one of the 26 .cu files fails with
nvcc error : 'cudafe++' died with status 0xC0000409
which is cudafe++ crashing on standard library headers from an MSVC newer than
the toolkit supports -- CUDA 12.4 tops out at MSVC 19.39 (VS 17.9) and the
installed toolset is 19.42 (VS 17.12). -allow-unsupported-compiler silences the
version check but not the incompatibility. The fix is either a newer CUDA or
the v14.39 MSVC toolset installed alongside, both of which are the maintainer's
call rather than something to do to their machine unasked.
The README now carries that diagnosis, since the error message names neither
the cause nor the fix. It also records the related trap that cost time here:
the CUDA installer places its MSBuild integration in whichever Visual Studio
instance it finds, so with both Build Tools and a full VS install present the
Visual Studio generator can report "No CUDA toolset found" while nvcc sits on
PATH. Ninja bypasses that integration entirely.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Nothing here changes behaviour; these are comments that actively mislead. All twelve src/math/*.cpp files opened with "Function declarations for ...". They contain the definitions -- the declarations are in the headers, which say so correctly. src/bindings/uniform.cpp and src/bindings/negative_binomial.cpp both claimed to bind normal.cpp and bernoulli.cpp respectively, from copy-paste. Two comments in poisson.cpp had been hard-wrapped at roughly twenty columns, leaving "// Poisson is defined / // on non-negative / // integers" and a five-line fragment of the log-PMF formula. Reflowed, and the log-space one now says why it is in log space rather than only that it is. The gamma CDF section header said "CDF using series / continued fraction", which names the two branches without saying which applies where or why there are two. It now states the crossover and the overflow reason for log space -- the things a reader would otherwise have to derive from the dispatch below it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
gamma, chi_square and beta CDFs had no coverage in the suite, which is how their defects survived: the benchmark checks agreement with SciPy before timing, so a case here would have flagged the beta series and the incomplete gamma immediately. They have no *_cpu batch entry point, so they go in the scalar group. They are also the most expensive routines in the library, so a regression matters more there than anywhere else. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Two runs of an untouched normal_cdf differed by 8% and compare.py called it a regression. Re-measuring put the true figure below both, so the machine was slow for the whole run, not the code. Fixed at both ends. The scalar group now times 21 rounds rather than 7: those loops are dominated by per-call Python overhead, which the interpreter varies far more than it varies compiled work, so the minimum needs more samples to settle. And compare.py will not flag a change smaller than the two runs' combined noise_pct, printing it as "within noise" instead of counting it. The docstring is honest that this is not a complete guard: noise_pct measures spread within a run, so it cannot see a run that was uniformly slow. A flagged case is a prompt to re-run, not a verdict. With both in place the phantom disappeared -- normal_cdf came back at -8.7%, its original level -- and the run has no regressions. BENCHMARKS.md gains an entry for the correctness work whose headline is not a speedup: three CDFs were returning wrong answers, two of them probabilities above 1.0. Worth recording in the performance log because the benchmark suite is what surfaced the first of them, by checking agreement with SciPy before timing anything. A wrong result cannot quietly post a good number. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
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
1 participant
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.
Branches off
repo-organization. Six commits, each independently revertible.What is here
1. Fix the flaky RNG tests (closes #2) —
4f9a5dfThe tolerances were round numbers rather than multiples of the estimator's standard error, so several sat near 2 sigma and failed at about the rate a 2 sigma bound fails. I re-derived every standard error independently and reproduced the issue's numbers exactly; my predicted suite failure rate was 7.87% against the 7.7% measured over 800 runs.
Audited all 21 RNG assertions, not just the three named. Three more sat below 7 sigma, including
test_discrete_uniform.cpp:103— the mean check on the line directly above the flagged variance check, sharing its literal, failing ~1 run in 300. Fixing only the three listed would have left the suite flaky at a lower, harder-to-diagnose rate.Drops
--repeat until-pass:3from CI.2. Seedable shared RNG —
7671393Every sampler declared its own
thread_localengine inline, so the stream was unreachable and reseeded from OS entropy each run. Hoisted into one shared accessor; all twelve samplers rewired. Addsfastdist.seed()/seed_from_entropy()and the C++ equivalents.Two limits are documented at the API rather than left to be discovered: seeding is per-thread, and
std::mt19937is standardised bit-for-bit while the distribution adaptors on top of it are not — so a seed reproduces a run on one toolchain, not across libstdc++/libc++/MSVC.With streams pinned, every RNG block seeds before sampling and is deterministic, which let tolerances come down from ~7–70 sigma to a uniform ~5 sigma. Sensitivity to systematic bias improves up to 12x at the loosest sites.
3. Benchmark suite and evidence log —
7482c06benchmarks/measures against SciPy and saves JSON per run;compare.pygates regressions at 5%;table.pygenerates the Markdown soBENCHMARKS.mdentries are never transcribed by hand. Reports the minimum round rather than the mean, and records how quiet the machine was separately. Every case verifies both implementations agree numerically before either is timed.The v0.1.0 baseline is recorded including the unflattering results.
4. Discrete CDF recurrences and batch invariant hoisting —
73ba7d6Both problems were found by the benchmark, not by inspection.
The three discrete CDFs summed PMF terms, re-deriving each from scratch. They now use recurrences, with log-space fallbacks guarded on the first term's magnitude for the parameter ranges where it underflows. Batch paths hoist validation and loop-invariant terms; scalar and batch share one inlined core so they cannot drift, and the arithmetic is ordered as before so results stay bit-identical.
29 cases improved, none regressed:
poisson_cdfnormal_logpdflog(sigma)per element)uniform_pdfnormal_cdfAlso fixes a documented overflow:
negative_binomial_pmf_scalarbuiltC(k+r-1, k)from rawtgammas, going inf at k=170 and nan beyond. Eight tests werexfail(strict)against it; they now pass, so the markers are gone and the cases stay as regression tests.5–6. Log entry and documentation —
8db7afb,3c568c5README build instructions were stale (
setup.py bdist_wheelagainst a hand-configured tree;pip install .does it all). Adds direct CMake instructions including what to do when it cannot find pybind11 — the first thing a fresh checkout hits. Corrects the supported Python range, which said 3.12–3.14 while every other file says 3.10.New sections for reproducible sampling, benchmarks, and the tolerance policy for sampling tests.
All thirteen distribution modules swallowed the original exception in their import guard and asserted a cause they had not checked — a missing numpy reported itself as a missing C++ core. Now chained with
from exc.Verification
Not done
extras\visual_studio_integration\MSBuildExtensions, but those files are not in either VS instance'sBuildCustomizations, and VS Community has no C++ workload at all. Changing CUDA code I cannot compile or test seemed worse than leaving it.🤖 Generated with Claude Code