From 4f9a5df93a24c5be1075ec3f74462b3454313eb7 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:03:29 -0400 Subject: [PATCH 01/15] Size RNG test tolerances from estimator standard error 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 --- .github/workflows/python-distro.yml | 7 +++---- tests/cpp/test_discrete_uniform.cpp | 8 ++++++-- tests/cpp/test_geometric.cpp | 5 ++++- tests/cpp/test_negative_binomial.cpp | 7 +++++-- tests/cpp/test_normal.cpp | 4 +++- tests/cpp/test_uniform.cpp | 7 +++++-- 6 files changed, 26 insertions(+), 12 deletions(-) diff --git a/.github/workflows/python-distro.yml b/.github/workflows/python-distro.yml index bbc0609..3f04c39 100644 --- a/.github/workflows/python-distro.yml +++ b/.github/workflows/python-distro.yml @@ -46,11 +46,10 @@ jobs: - name: Build C++ tests run: cmake --build build --target fastdist_tests --parallel - # Current C++ tests have flakes (geometric, discrete uniform, and uniform) - # Tracking: ghosteau/fastdist#2 - # Remove the --repeat once that issue is closed. + # RNG tolerances are sized from the estimator standard error (>=7 sigma), + # so a failure here is a real regression rather than a flake. - name: Run C++ tests - run: ctest --test-dir build --output-on-failure --repeat until-pass:3 + run: ctest --test-dir build --output-on-failure # Build isolation is left on so this also validates that # build-system.requires in pyproject.toml is complete. diff --git a/tests/cpp/test_discrete_uniform.cpp b/tests/cpp/test_discrete_uniform.cpp index 9471ead..3d41440 100644 --- a/tests/cpp/test_discrete_uniform.cpp +++ b/tests/cpp/test_discrete_uniform.cpp @@ -100,8 +100,12 @@ void test_discrete_uniform() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::discrete_uniform_mean(a, b)) < 1e-2); - assert(std::abs(var - fastdist::math::discrete_uniform_variance(a, b)) < 1e-2); + // Tolerances are sized from the estimator standard error at N, not round + // numbers: SE(mean) = 3.4e-3 and SE(var) = 5.0e-3 here, so these sit at + // ~8.8 and ~10 sigma. Anything near 2 sigma flakes a few percent of runs. + // See ghosteau/fastdist#2. + assert(std::abs(mean - fastdist::math::discrete_uniform_mean(a, b)) < 0.03); + assert(std::abs(var - fastdist::math::discrete_uniform_variance(a, b)) < 0.05); } std::cout << "Discrete uniform tests passed.\n"; diff --git a/tests/cpp/test_geometric.cpp b/tests/cpp/test_geometric.cpp index 8d76a0b..5e65808 100644 --- a/tests/cpp/test_geometric.cpp +++ b/tests/cpp/test_geometric.cpp @@ -114,7 +114,10 @@ void test_geometric() { const double var = sumsq / N - mean * mean; assert(std::abs(mean - fastdist::math::geometric_mean(p)) < 0.15); - assert(std::abs(var - fastdist::math::geometric_variance(p)) < 0.15); + // SE(var) = 0.068 at N, so 0.5 is ~7.3 sigma. The mean assertion above + // shares the 0.15 literal but has SE 6.9e-3, which is ~22 sigma already. + // See ghosteau/fastdist#2. + assert(std::abs(var - fastdist::math::geometric_variance(p)) < 0.5); } std::cout << "Geometric distribution tests passed.\n"; diff --git a/tests/cpp/test_negative_binomial.cpp b/tests/cpp/test_negative_binomial.cpp index 9f22f5c..84536ac 100644 --- a/tests/cpp/test_negative_binomial.cpp +++ b/tests/cpp/test_negative_binomial.cpp @@ -101,9 +101,12 @@ void test_negative_binomial() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::negative_binomial_mean(r, p)) < 5e-2); + // SE(mean) = 0.011 and SE(var) = 0.117 at N, so these are ~9 and ~8.5 + // sigma; the previous values sat near 4.3 sigma. + // See ghosteau/fastdist#2. + assert(std::abs(mean - fastdist::math::negative_binomial_mean(r, p)) < 0.1); - assert(std::abs(var - fastdist::math::negative_binomial_variance(r, p)) < 5e-1); + assert(std::abs(var - fastdist::math::negative_binomial_variance(r, p)) < 1.0); } std::cout << "Negative Binomial distribution tests passed!\n"; diff --git a/tests/cpp/test_normal.cpp b/tests/cpp/test_normal.cpp index ab65ffc..37e5522 100644 --- a/tests/cpp/test_normal.cpp +++ b/tests/cpp/test_normal.cpp @@ -118,7 +118,9 @@ void test_normal() { const double var = sumsq / N - mean * mean; assert(std::abs(mean - mu) < 5e-2); - assert(std::abs(var - sigma * sigma) < 5e-2); + // SE(var) = 0.013 at N, so 0.1 is ~7.9 sigma; 5e-2 was 3.95 sigma. + // See ghosteau/fastdist#2. + assert(std::abs(var - sigma * sigma) < 0.1); } // ------------------------- diff --git a/tests/cpp/test_uniform.cpp b/tests/cpp/test_uniform.cpp index 7d4b3c1..dbc9b18 100644 --- a/tests/cpp/test_uniform.cpp +++ b/tests/cpp/test_uniform.cpp @@ -97,8 +97,11 @@ void test_uniform() { double expected_mean = 0.5 * (a + b); double expected_var = (b - a) * (b - a) / 12.0; - assert(std::abs(mean - expected_mean) < 5e-3); - assert(std::abs(var - expected_var) < 5e-3); + // SE(mean) = 1.7e-3 and SE(var) = 1.3e-3 at N, so 0.02 is ~11.5 and ~15 + // sigma. The previous 5e-3 put the mean check at 2.9 sigma. + // See ghosteau/fastdist#2. + assert(std::abs(mean - expected_mean) < 0.02); + assert(std::abs(var - expected_var) < 0.02); } std::cout << "Continuous uniform distribution tests passed!\n"; From 7671393bc4524020842382b5fe320156de64b616 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:11:14 -0400 Subject: [PATCH 02/15] Add a seedable shared RNG and make the sampling tests deterministic 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 --- CMakeLists.txt | 5 +- include/fastdist/math/rng.h | 39 +++++++++++ python/fastdist/__init__.py | 4 +- src/bindings/bindings.cpp | 2 + src/bindings/rng.cpp | 26 ++++++++ src/math/bernoulli.cpp | 4 +- src/math/beta.cpp | 6 +- src/math/binomial.cpp | 4 +- src/math/discrete_uniform.cpp | 4 +- src/math/exponential.cpp | 4 +- src/math/gamma.cpp | 4 +- src/math/geometric.cpp | 4 +- src/math/negative_binomial.cpp | 4 +- src/math/normal.cpp | 7 +- src/math/poisson.cpp | 6 +- src/math/rng.cpp | 21 ++++++ src/math/uniform.cpp | 6 +- tests/cpp/test_app.cpp | 34 +++++----- tests/cpp/test_bernoulli.cpp | 10 ++- tests/cpp/test_beta.cpp | 10 ++- tests/cpp/test_binomial.cpp | 10 ++- tests/cpp/test_chi_square.cpp | 10 ++- tests/cpp/test_discrete_uniform.cpp | 14 ++-- tests/cpp/test_exponential.cpp | 10 ++- tests/cpp/test_gamma.cpp | 10 ++- tests/cpp/test_geometric.cpp | 13 ++-- tests/cpp/test_negative_binomial.cpp | 13 ++-- tests/cpp/test_normal.cpp | 16 +++-- tests/cpp/test_poisson.cpp | 10 ++- tests/cpp/test_rng.cpp | 96 ++++++++++++++++++++++++++++ tests/cpp/test_uniform.cpp | 13 ++-- tests/python/test_package_api.py | 15 ++++- tests/python/test_rng.py | 84 ++++++++++++++++++++++++ 33 files changed, 431 insertions(+), 87 deletions(-) create mode 100644 include/fastdist/math/rng.h create mode 100644 src/bindings/rng.cpp create mode 100644 src/math/rng.cpp create mode 100644 tests/cpp/test_rng.cpp create mode 100644 tests/python/test_rng.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 00e23a9..69dcfce 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -74,6 +74,7 @@ pybind11_add_module(_fastdist src/bindings/beta.cpp src/bindings/chi_square.cpp src/bindings/utils.cpp + src/bindings/rng.cpp src/wrappers/normal_wrapper.cpp src/wrappers/poisson_wrapper.cpp @@ -102,6 +103,7 @@ add_library(fastdist_core STATIC src/api/chi_square.cpp src/api/utils.cpp + src/math/rng.cpp src/math/normal.cpp src/math/exponential.cpp src/math/poisson.cpp @@ -243,6 +245,7 @@ add_executable(fastdist_tests tests/cpp/test_beta.cpp tests/cpp/test_chi_square.cpp tests/cpp/test_utils.cpp + tests/cpp/test_rng.cpp tests/cpp/test_app.cpp ) @@ -268,7 +271,7 @@ target_compile_options(fastdist_tests PRIVATE # works. Each runs as its own process. set(FASTDIST_TEST_CASES bernoulli binomial negative_binomial discrete_uniform exponential - geometric normal poisson uniform beta gamma chi_square utils + geometric normal poisson uniform beta gamma chi_square utils rng ) foreach (case IN LISTS FASTDIST_TEST_CASES) diff --git a/include/fastdist/math/rng.h b/include/fastdist/math/rng.h new file mode 100644 index 0000000..a3efafa --- /dev/null +++ b/include/fastdist/math/rng.h @@ -0,0 +1,39 @@ +// Header file for the shared pseudo-random number generator +#ifndef RNG_H +#define RNG_H + +#include +#include + +namespace fastdist::math { + // The engine every *_sample() function draws from. + // + // The engine is thread_local: each thread owns an independent stream, so + // concurrent sampling needs no locking and threads never interleave draws. + // The flip side is that seeding is also per-thread -- see seed_rng below. + std::mt19937& rng(); + + // Pins the calling thread's stream to a fixed sequence. The same seed + // reproduces the same draws on every run, which makes sampling-based tests + // deterministic and lets callers reproduce a result exactly. + // + // Two caveats worth knowing: + // + // 1. This seeds the calling thread only. A thread that has not been seeded + // keeps its own entropy-initialised stream. + // + // 2. std::mt19937 is specified bit-for-bit by the standard, but the + // distribution adaptors layered on top of it (std::normal_distribution + // and friends) are not. Identical engine state therefore yields + // different samples across libstdc++, libc++ and MSVC. A seed makes a + // run reproducible on one platform and toolchain, not across all of + // them. + void seed_rng(std::uint32_t value); + + // Returns the calling thread's stream to non-deterministic behaviour by + // drawing a fresh seed from the OS entropy source. This is the state every + // thread starts in, so it is only needed to undo a previous seed_rng call. + void seed_rng_from_entropy(); +} // namespace fastdist::math + +#endif // RNG_H diff --git a/python/fastdist/__init__.py b/python/fastdist/__init__.py index b943f97..3c7e6b7 100644 --- a/python/fastdist/__init__.py +++ b/python/fastdist/__init__.py @@ -2,6 +2,7 @@ from importlib.metadata import PackageNotFoundError, version as _pkg_version from . import _fastdist +from ._fastdist import seed, seed_from_entropy from .distributions import ( Bernoulli, Beta, Binomial, ChiSquare, DiscreteUniform, Exponential, Gamma, Geometric, NegativeBinomial, Normal, @@ -17,4 +18,5 @@ "Bernoulli", "Beta", "Binomial", "ChiSquare", "DiscreteUniform", "Exponential", "Gamma", "Geometric", "NegativeBinomial", "Normal", - "Poisson", "Uniform", "Utils"] \ No newline at end of file + "Poisson", "Uniform", "Utils", + "seed", "seed_from_entropy"] \ No newline at end of file diff --git a/src/bindings/bindings.cpp b/src/bindings/bindings.cpp index 136ff88..2d033a9 100644 --- a/src/bindings/bindings.cpp +++ b/src/bindings/bindings.cpp @@ -19,6 +19,7 @@ void bind_gamma(py::module_ &m); void bind_beta(py::module_ &m); void bind_chi_square(py::module_ &m); void bind_utils(py::module_ &m); +void bind_rng(py::module_ &m); // Pybind Module @@ -36,6 +37,7 @@ PYBIND11_MODULE(_fastdist, m) { bind_beta(m); bind_chi_square(m); bind_utils(m); + bind_rng(m); m.attr("__version__") = FASTDIST_VERSION_STRING; } diff --git a/src/bindings/rng.cpp b/src/bindings/rng.cpp new file mode 100644 index 0000000..315aa36 --- /dev/null +++ b/src/bindings/rng.cpp @@ -0,0 +1,26 @@ +// pybind11 bindings for /src/math/rng.cpp +#include "fastdist/math/rng.h" +#include "pybind11/pybind11.h" + +namespace py = pybind11; + +void bind_rng(py::module_ &m) { + m.def("seed", &fastdist::math::seed_rng, py::arg("value"), + R"pbdoc(Seed the sampling engine so draws are reproducible. + +Pins the calling thread's random stream to a fixed sequence: the same seed +replays the same draws on every run. + +Two caveats. The seed applies to the calling thread only, so a worker thread +that has not been seeded keeps its own entropy-initialised stream. And while +the underlying Mersenne Twister engine is specified bit-for-bit by the C++ +standard, the distribution adaptors built on it are not -- the same seed +therefore yields different samples on Linux, macOS and Windows. A seed makes a +run reproducible on one platform and toolchain, not across all of them.)pbdoc"); + + m.def("seed_from_entropy", &fastdist::math::seed_rng_from_entropy, + R"pbdoc(Return the sampling engine to non-deterministic behaviour. + +Draws a fresh seed from the OS entropy source. This is the state every thread +starts in, so it is only needed to undo a previous seed() call.)pbdoc"); +} diff --git a/src/math/bernoulli.cpp b/src/math/bernoulli.cpp index 1dc27aa..863400b 100644 --- a/src/math/bernoulli.cpp +++ b/src/math/bernoulli.cpp @@ -1,6 +1,7 @@ // Function declarations for Bernoulli distribution functions #include #include +#include #include #include @@ -73,10 +74,9 @@ namespace fastdist::math { return -1; // signal invalid input } - thread_local std::mt19937 rng{std::random_device{}()}; std::bernoulli_distribution dist(p); - return dist(rng) ? 1 : 0; + return dist(rng()) ? 1 : 0; } // Batch Functions diff --git a/src/math/beta.cpp b/src/math/beta.cpp index e004f8c..3e7f486 100644 --- a/src/math/beta.cpp +++ b/src/math/beta.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include @@ -63,11 +64,10 @@ namespace fastdist::math { // RNG // ------------------------- double beta_sample(const double alpha, const double beta) { - thread_local std::mt19937 rng{std::random_device{}()}; std::gamma_distribution ga(alpha, 1.0); std::gamma_distribution gb(beta, 1.0); - double a = ga(rng); - double b = gb(rng); + double a = ga(rng()); + double b = gb(rng()); return a / (a + b); } diff --git a/src/math/binomial.cpp b/src/math/binomial.cpp index ac9c9e1..100693e 100644 --- a/src/math/binomial.cpp +++ b/src/math/binomial.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include @@ -89,10 +90,9 @@ namespace fastdist::math { return -1; // signal invalid input } - thread_local std::mt19937 rng{std::random_device{}()}; std::binomial_distribution dist(n, p); - return dist(rng); + return dist(rng()); } } // namespace fastdist::math diff --git a/src/math/discrete_uniform.cpp b/src/math/discrete_uniform.cpp index 7425211..56f037c 100644 --- a/src/math/discrete_uniform.cpp +++ b/src/math/discrete_uniform.cpp @@ -1,6 +1,7 @@ // Function declarations for discrete uniform distribution functions #include #include +#include #include #include @@ -86,10 +87,9 @@ namespace fastdist::math { return std::numeric_limits::min(); } - thread_local std::mt19937 rng{std::random_device{}()}; std::uniform_int_distribution dist(a, b); - return dist(rng); + return dist(rng()); } } // namespace fastdist::math diff --git a/src/math/exponential.cpp b/src/math/exponential.cpp index 5f8b924..aba4e0d 100644 --- a/src/math/exponential.cpp +++ b/src/math/exponential.cpp @@ -1,6 +1,7 @@ // Function declarations for exponential distribution functions #include "fastdist/math/exponential.h" #include +#include #include #include @@ -70,10 +71,9 @@ namespace fastdist::math { return std::numeric_limits::quiet_NaN(); } - thread_local std::mt19937 rng{std::random_device{}()}; std::exponential_distribution dist(lambda); - return dist(rng); + return dist(rng()); } // Batch Functions diff --git a/src/math/gamma.cpp b/src/math/gamma.cpp index 87e0f83..8740abe 100644 --- a/src/math/gamma.cpp +++ b/src/math/gamma.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include @@ -84,9 +85,8 @@ namespace fastdist::math { if (!std::isfinite(alpha) || !std::isfinite(theta) || alpha <= 0.0 || theta <= 0.0) return std::numeric_limits::quiet_NaN(); - thread_local std::mt19937 rng{std::random_device{}()}; std::gamma_distribution dist(alpha, theta); - return dist(rng); + return dist(rng()); } // ------------------------- diff --git a/src/math/geometric.cpp b/src/math/geometric.cpp index 879b6ec..7eba2f6 100644 --- a/src/math/geometric.cpp +++ b/src/math/geometric.cpp @@ -1,6 +1,7 @@ // Function declarations for geometric distribution functions #include #include +#include #include #include @@ -90,11 +91,10 @@ namespace fastdist::math { return -1; } - thread_local std::mt19937 rng{std::random_device{}()}; std::geometric_distribution dist(p); // std::geometric_distribution returns #failures before first success - return dist(rng) + 1; + return dist(rng()) + 1; } } // namespace fastdist::math diff --git a/src/math/negative_binomial.cpp b/src/math/negative_binomial.cpp index b553044..3398e10 100644 --- a/src/math/negative_binomial.cpp +++ b/src/math/negative_binomial.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include @@ -108,9 +109,8 @@ namespace fastdist::math { return -1; // invalid input } - thread_local std::mt19937 rng{std::random_device{}()}; std::negative_binomial_distribution dist(r, p); - return dist(rng); + return dist(rng()); } } // namespace fastdist::math diff --git a/src/math/normal.cpp b/src/math/normal.cpp index e30c277..da84fff 100644 --- a/src/math/normal.cpp +++ b/src/math/normal.cpp @@ -1,6 +1,7 @@ // Function declarations for normal distribution functions #include #include +#include #include #include #include @@ -74,18 +75,16 @@ namespace fastdist::math { if (!std::isfinite(mu) || !std::isfinite(sigma) || sigma <= 0.0) { return std::numeric_limits::quiet_NaN(); } - thread_local std::mt19937 rng{std::random_device{}()}; std::normal_distribution dist(mu, sigma); - return dist(rng); + return dist(rng()); } double normal_log_sample(const double mu, const double sigma) { if (!std::isfinite(mu) || !std::isfinite(sigma) || sigma <= 0.0) { return std::numeric_limits::quiet_NaN(); } - thread_local std::mt19937 rng{std::random_device{}()}; std::normal_distribution dist(mu, sigma); - return std::exp(dist(rng)); + return std::exp(dist(rng())); } double z_score(const double x, const double mu, const double sigma) { return (x - mu) / sigma; } diff --git a/src/math/poisson.cpp b/src/math/poisson.cpp index c9891a4..57ca3ad 100644 --- a/src/math/poisson.cpp +++ b/src/math/poisson.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include #include @@ -86,14 +87,13 @@ namespace fastdist::math { return lambda * (std::exp(t) - 1.0); } - // RNG via std::poisson_distribution + // X ~ Poisson(lambda) int poisson_sample(const double lambda) { if (!std::isfinite(lambda) || lambda <= 0.0) { return -1; // signal invalid input } - thread_local std::mt19937 rng{std::random_device{}()}; std::poisson_distribution dist(lambda); - return dist(rng); + return dist(rng()); } // Batch Functions diff --git a/src/math/rng.cpp b/src/math/rng.cpp new file mode 100644 index 0000000..97441c7 --- /dev/null +++ b/src/math/rng.cpp @@ -0,0 +1,21 @@ +// Function definitions for the shared pseudo-random number generator +#include +#include +#include + +namespace fastdist::math { + + std::mt19937& rng() { + // Function-local rather than namespace-scope so the engine is created on + // first use. A namespace-scope thread_local would be constructed on every + // thread that touches the library, including threads that never sample, + // and would pay a std::random_device read to do it. + thread_local std::mt19937 engine{std::random_device{}()}; + return engine; + } + + void seed_rng(const std::uint32_t value) { rng().seed(value); } + + void seed_rng_from_entropy() { rng().seed(std::random_device{}()); } + +} // namespace fastdist::math diff --git a/src/math/uniform.cpp b/src/math/uniform.cpp index 28e2d53..cfee811 100644 --- a/src/math/uniform.cpp +++ b/src/math/uniform.cpp @@ -1,5 +1,6 @@ // Function declarations for continuous uniform distribution functions #include +#include #include #include #include @@ -83,15 +84,14 @@ namespace fastdist::math { return std::log(mgf); } - // RNG: simple thread-local uniform_real_distribution + // X ~ Uniform(a, b) double uniform_sample(const double a, const double b) { if (!std::isfinite(a) || !std::isfinite(b) || a >= b) { return std::numeric_limits::quiet_NaN(); } - thread_local std::mt19937 rng{std::random_device{}()}; std::uniform_real_distribution dist(a, b); - return dist(rng); + return dist(rng()); } // Batch Functions diff --git a/tests/cpp/test_app.cpp b/tests/cpp/test_app.cpp index 45996a0..d243c89 100644 --- a/tests/cpp/test_app.cpp +++ b/tests/cpp/test_app.cpp @@ -20,28 +20,30 @@ void test_normal(); void test_poisson(); void test_uniform(); void test_utils(); +void test_rng(); // Run one case by name, or all of them when given no argument. Keeping the // list ordered (rather than a map) preserves the original execution order. int main(int argc, char** argv) { const std::vector> tests{ - {"bernoulli", test_bernoulli}, - {"binomial", test_binomial}, - {"negative_binomial", test_negative_binomial}, - {"discrete_uniform", test_discrete_uniform}, - {"exponential", test_exponential}, - {"geometric", test_geometric}, - {"normal", test_normal}, - {"poisson", test_poisson}, - {"uniform", test_uniform}, - {"beta", test_beta}, - {"gamma", test_gamma}, - {"chi_square", test_chi_square}, - {"utils", test_utils}, + {"bernoulli", test_bernoulli}, + {"binomial", test_binomial}, + {"negative_binomial", test_negative_binomial}, + {"discrete_uniform", test_discrete_uniform}, + {"exponential", test_exponential}, + {"geometric", test_geometric}, + {"normal", test_normal}, + {"poisson", test_poisson}, + {"uniform", test_uniform}, + {"beta", test_beta}, + {"gamma", test_gamma}, + {"chi_square", test_chi_square}, + {"utils", test_utils}, + {"rng", test_rng}, }; const auto lookup = [&tests](const std::string& name) -> void (*)() { - for (const auto& [n, fn] : tests) { + for (const auto& [n, fn]: tests) { if (n == name) { return fn; } @@ -65,7 +67,7 @@ int main(int argc, char** argv) { selected.push_back(fn); } - for (const auto fn : selected) { + for (const auto fn: selected) { fn(); } return 0; @@ -73,7 +75,7 @@ int main(int argc, char** argv) { std::cout << "Starting fastdist C++ tests...\n"; - for (const auto& [name, fn] : tests) { + for (const auto& [name, fn]: tests) { fn(); } diff --git a/tests/cpp/test_bernoulli.cpp b/tests/cpp/test_bernoulli.cpp index 2bf8e96..87296d6 100644 --- a/tests/cpp/test_bernoulli.cpp +++ b/tests/cpp/test_bernoulli.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include void test_bernoulli() { @@ -39,6 +40,11 @@ void test_bernoulli() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 9.2e-4, SE(var) = 3.7e-4 at this N). + fastdist::math::seed_rng(1001); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -52,8 +58,8 @@ void test_bernoulli() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::bernoulli_mean(p)) < 1e-2); - assert(std::abs(var - fastdist::math::bernoulli_variance(p)) < 1e-2); + assert(std::abs(mean - fastdist::math::bernoulli_mean(p)) < 0.005); + assert(std::abs(var - fastdist::math::bernoulli_variance(p)) < 0.002); } std::cout << "Bernoulli tests passed\n"; diff --git a/tests/cpp/test_beta.cpp b/tests/cpp/test_beta.cpp index fe79db7..3efe5f5 100644 --- a/tests/cpp/test_beta.cpp +++ b/tests/cpp/test_beta.cpp @@ -1,6 +1,7 @@ // Unit tests for Beta distribution #include #include +#include #include #include @@ -52,6 +53,11 @@ void test_beta() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 3.2e-4, SE(var) = 7.0e-5 at this N). + fastdist::math::seed_rng(1002); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -70,8 +76,8 @@ void test_beta() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::beta_mean(alpha, beta)) < 5e-3); - assert(std::abs(var - fastdist::math::beta_variance(alpha, beta)) < 5e-3); + assert(std::abs(mean - fastdist::math::beta_mean(alpha, beta)) < 0.002); + assert(std::abs(var - fastdist::math::beta_variance(alpha, beta)) < 4e-4); } std::cout << "Beta distribution tests passed!\n"; diff --git a/tests/cpp/test_binomial.cpp b/tests/cpp/test_binomial.cpp index b3d06ed..1e5c9ac 100644 --- a/tests/cpp/test_binomial.cpp +++ b/tests/cpp/test_binomial.cpp @@ -1,6 +1,7 @@ // Unit tests for Binomial distribution #include #include +#include #include #include @@ -68,6 +69,11 @@ void test_binomial() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 2.9e-3, SE(var) = 5.8e-3 at this N). + fastdist::math::seed_rng(1003); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -85,8 +91,8 @@ void test_binomial() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::binomial_mean(n, p)) < 0.15); - assert(std::abs(var - fastdist::math::binomial_variance(n, p)) < 0.15); + assert(std::abs(mean - fastdist::math::binomial_mean(n, p)) < 0.015); + assert(std::abs(var - fastdist::math::binomial_variance(n, p)) < 0.03); } std::cout << "Binomial tests passed.\n"; diff --git a/tests/cpp/test_chi_square.cpp b/tests/cpp/test_chi_square.cpp index c35068c..a6ed2c5 100644 --- a/tests/cpp/test_chi_square.cpp +++ b/tests/cpp/test_chi_square.cpp @@ -1,6 +1,7 @@ // Unit tests for Chi-square distribution #include #include +#include #include #include @@ -79,6 +80,11 @@ void test_chi_square() { // RNG sanity check // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 6.9e-3, SE(var) = 4.8e-2 at this N). + fastdist::math::seed_rng(1004); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -95,9 +101,9 @@ void test_chi_square() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::chi_square_mean(k)) < 5e-2); + assert(std::abs(mean - fastdist::math::chi_square_mean(k)) < 0.04); - assert(std::abs(var - fastdist::math::chi_square_variance(k)) < 5e-1); + assert(std::abs(var - fastdist::math::chi_square_variance(k)) < 0.25); } std::cout << "Chi-square distribution tests passed!\n"; diff --git a/tests/cpp/test_discrete_uniform.cpp b/tests/cpp/test_discrete_uniform.cpp index 3d41440..5ae7114 100644 --- a/tests/cpp/test_discrete_uniform.cpp +++ b/tests/cpp/test_discrete_uniform.cpp @@ -1,6 +1,7 @@ // Unit tests for Discrete Uniform distribution #include #include +#include #include #include @@ -82,6 +83,11 @@ void test_discrete_uniform() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 3.4e-3, SE(var) = 5.0e-3 at this N). + fastdist::math::seed_rng(1005); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -100,12 +106,8 @@ void test_discrete_uniform() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - // Tolerances are sized from the estimator standard error at N, not round - // numbers: SE(mean) = 3.4e-3 and SE(var) = 5.0e-3 here, so these sit at - // ~8.8 and ~10 sigma. Anything near 2 sigma flakes a few percent of runs. - // See ghosteau/fastdist#2. - assert(std::abs(mean - fastdist::math::discrete_uniform_mean(a, b)) < 0.03); - assert(std::abs(var - fastdist::math::discrete_uniform_variance(a, b)) < 0.05); + assert(std::abs(mean - fastdist::math::discrete_uniform_mean(a, b)) < 0.02); + assert(std::abs(var - fastdist::math::discrete_uniform_variance(a, b)) < 0.025); } std::cout << "Discrete uniform tests passed.\n"; diff --git a/tests/cpp/test_exponential.cpp b/tests/cpp/test_exponential.cpp index 6582611..3d8857b 100644 --- a/tests/cpp/test_exponential.cpp +++ b/tests/cpp/test_exponential.cpp @@ -1,6 +1,7 @@ // Unit tests for Exponential distribution #include #include +#include #include #include @@ -95,6 +96,11 @@ void test_exponential() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 1.0e-3, SE(var) = 1.4e-3 at this N). + fastdist::math::seed_rng(1006); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -112,8 +118,8 @@ void test_exponential() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::exponential_mean(lambda)) < 1e-2); - assert(std::abs(var - fastdist::math::exponential_variance(lambda)) < 1e-2); + assert(std::abs(mean - fastdist::math::exponential_mean(lambda)) < 0.005); + assert(std::abs(var - fastdist::math::exponential_variance(lambda)) < 0.0075); } std::cout << "Exponential distribution tests passed.\n"; diff --git a/tests/cpp/test_gamma.cpp b/tests/cpp/test_gamma.cpp index ff4492f..06ef7fb 100644 --- a/tests/cpp/test_gamma.cpp +++ b/tests/cpp/test_gamma.cpp @@ -1,6 +1,7 @@ // Unit tests for Gamma distribution #include #include +#include #include #include @@ -79,6 +80,11 @@ void test_gamma() { // RNG sanity check // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 6.9e-3, SE(var) = 4.8e-2 at this N). + fastdist::math::seed_rng(1007); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -96,9 +102,9 @@ void test_gamma() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::gamma_mean(alpha, theta)) < 5e-2); + assert(std::abs(mean - fastdist::math::gamma_mean(alpha, theta)) < 0.04); - assert(std::abs(var - fastdist::math::gamma_variance(alpha, theta)) < 5e-1); + assert(std::abs(var - fastdist::math::gamma_variance(alpha, theta)) < 0.25); } std::cout << "Gamma distribution tests passed!\n"; diff --git a/tests/cpp/test_geometric.cpp b/tests/cpp/test_geometric.cpp index 5e65808..a3adcc0 100644 --- a/tests/cpp/test_geometric.cpp +++ b/tests/cpp/test_geometric.cpp @@ -1,6 +1,7 @@ // Unit tests for Geometric distribution #include #include +#include #include #include @@ -96,6 +97,11 @@ void test_geometric() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 6.9e-3, SE(var) = 6.8e-2 at this N). + fastdist::math::seed_rng(1008); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -113,11 +119,8 @@ void test_geometric() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - fastdist::math::geometric_mean(p)) < 0.15); - // SE(var) = 0.068 at N, so 0.5 is ~7.3 sigma. The mean assertion above - // shares the 0.15 literal but has SE 6.9e-3, which is ~22 sigma already. - // See ghosteau/fastdist#2. - assert(std::abs(var - fastdist::math::geometric_variance(p)) < 0.5); + assert(std::abs(mean - fastdist::math::geometric_mean(p)) < 0.04); + assert(std::abs(var - fastdist::math::geometric_variance(p)) < 0.4); } std::cout << "Geometric distribution tests passed.\n"; diff --git a/tests/cpp/test_negative_binomial.cpp b/tests/cpp/test_negative_binomial.cpp index 84536ac..6f7e732 100644 --- a/tests/cpp/test_negative_binomial.cpp +++ b/tests/cpp/test_negative_binomial.cpp @@ -1,6 +1,7 @@ // Unit tests for Negative Binomial distribution #include #include +#include #include #include @@ -84,6 +85,11 @@ void test_negative_binomial() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 1.1e-2, SE(var) = 1.2e-1 at this N). + fastdist::math::seed_rng(1009); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -101,12 +107,9 @@ void test_negative_binomial() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - // SE(mean) = 0.011 and SE(var) = 0.117 at N, so these are ~9 and ~8.5 - // sigma; the previous values sat near 4.3 sigma. - // See ghosteau/fastdist#2. - assert(std::abs(mean - fastdist::math::negative_binomial_mean(r, p)) < 0.1); + assert(std::abs(mean - fastdist::math::negative_binomial_mean(r, p)) < 0.075); - assert(std::abs(var - fastdist::math::negative_binomial_variance(r, p)) < 1.0); + assert(std::abs(var - fastdist::math::negative_binomial_variance(r, p)) < 0.75); } std::cout << "Negative Binomial distribution tests passed!\n"; diff --git a/tests/cpp/test_normal.cpp b/tests/cpp/test_normal.cpp index 37e5522..9ba682f 100644 --- a/tests/cpp/test_normal.cpp +++ b/tests/cpp/test_normal.cpp @@ -1,6 +1,7 @@ // Unit tests for Normal distribution #include #include +#include #include #include @@ -103,6 +104,11 @@ void test_normal() { // Normal RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 4.5e-3, SE(var) = 1.3e-2 at this N). + fastdist::math::seed_rng(1012); + constexpr int N = 200000; double sum = 0.0; double sumsq = 0.0; @@ -117,16 +123,16 @@ void test_normal() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - mu) < 5e-2); - // SE(var) = 0.013 at N, so 0.1 is ~7.9 sigma; 5e-2 was 3.95 sigma. - // See ghosteau/fastdist#2. - assert(std::abs(var - sigma * sigma) < 0.1); + assert(std::abs(mean - mu) < 0.025); + assert(std::abs(var - sigma * sigma) < 0.075); } // ------------------------- // Log-normal RNG tests // ------------------------- { + fastdist::math::seed_rng(1013); + constexpr int N = 250000; double sum = 0.0; @@ -137,7 +143,7 @@ void test_normal() { } const double log_mean = sum / N; - assert(std::abs(log_mean - mu) < 5e-2); + assert(std::abs(log_mean - mu) < 0.02); } std::cout << "Normal distribution tests passed.\n"; diff --git a/tests/cpp/test_poisson.cpp b/tests/cpp/test_poisson.cpp index 58da367..d5503fd 100644 --- a/tests/cpp/test_poisson.cpp +++ b/tests/cpp/test_poisson.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include void test_poisson() { @@ -49,6 +50,11 @@ void test_poisson() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 4.0e-3, SE(var) = 1.2e-2 at this N). + fastdist::math::seed_rng(1010); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -62,8 +68,8 @@ void test_poisson() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - lambda) < 5e-2); - assert(std::abs(var - lambda) < 5e-1); + assert(std::abs(mean - lambda) < 0.02); + assert(std::abs(var - lambda) < 0.075); } std::cout << "Poisson distribution tests passed!\n"; diff --git a/tests/cpp/test_rng.cpp b/tests/cpp/test_rng.cpp new file mode 100644 index 0000000..a0eccd7 --- /dev/null +++ b/tests/cpp/test_rng.cpp @@ -0,0 +1,96 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + // Draws a fixed-length trace from several samplers. Mixing distributions + // matters: they all share one engine, so this also pins down the order in + // which they consume it. + std::vector trace(const std::uint32_t seed, const int n = 64) { + fastdist::math::seed_rng(seed); + + std::vector out; + out.reserve(static_cast(n) * 3); + for (int i = 0; i < n; ++i) { + out.push_back(fastdist::math::uniform_sample(0.0, 1.0)); + out.push_back(fastdist::math::normal_sample(0.0, 1.0)); + out.push_back(static_cast(fastdist::math::bernoulli_sample(0.5))); + } + return out; + } +} // namespace + +void test_rng() { + std::cout << "Running RNG seeding tests...\n"; + + // ------------------------- + // Reproducibility + // ------------------------- + { + // The core contract: same seed, same draws. Compared bitwise rather than + // with a tolerance -- replaying a stream is exact, not approximate. + const std::vector first = trace(12345); + const std::vector second = trace(12345); + + assert(first == second); + } + + // ------------------------- + // Distinct seeds give distinct streams + // ------------------------- + { + // Guards against seed_rng being wired up to something that ignores its + // argument, which would make every "seeded" test vacuous. + const std::vector a = trace(1); + const std::vector b = trace(2); + + assert(a != b); + } + + // ------------------------- + // Re-seeding rewinds the stream + // ------------------------- + { + fastdist::math::seed_rng(777); + const double first_draw = fastdist::math::uniform_sample(0.0, 1.0); + + // Advance the engine well past that point. + for (int i = 0; i < 100; ++i) { + (void) fastdist::math::uniform_sample(0.0, 1.0); + } + + fastdist::math::seed_rng(777); + assert(fastdist::math::uniform_sample(0.0, 1.0) == first_draw); + } + + // ------------------------- + // Entropy reseeding escapes the fixed stream + // ------------------------- + { + // seed_from_entropy has to actually move off the seeded sequence. A + // single draw could collide by chance, so compare a run of them; the + // odds of a full trace repeating are nil. + fastdist::math::seed_rng(999); + const std::vector seeded = trace(999); + + fastdist::math::seed_rng_from_entropy(); + std::vector entropic; + entropic.reserve(seeded.size()); + for (std::size_t i = 0; i < seeded.size(); ++i) { + entropic.push_back(fastdist::math::uniform_sample(0.0, 1.0)); + } + + assert(entropic != seeded); + } + + // Leave the engine non-deterministic so a later test that forgot to seed + // does not silently inherit this one's fixed stream. + fastdist::math::seed_rng_from_entropy(); + + std::cout << "RNG seeding tests passed.\n"; +} diff --git a/tests/cpp/test_uniform.cpp b/tests/cpp/test_uniform.cpp index dbc9b18..b253d0b 100644 --- a/tests/cpp/test_uniform.cpp +++ b/tests/cpp/test_uniform.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -77,6 +78,11 @@ void test_uniform() { // RNG tests // ------------------------- { + // Seeded, so this block is deterministic: it either always passes or + // always fails on a given toolchain, never intermittently. Tolerances + // below are ~5x the estimator standard error (SE(mean) = 1.7e-3, SE(var) = 1.3e-3 at this N). + fastdist::math::seed_rng(1011); + constexpr int N = 250000; double sum = 0.0; double sumsq = 0.0; @@ -97,11 +103,8 @@ void test_uniform() { double expected_mean = 0.5 * (a + b); double expected_var = (b - a) * (b - a) / 12.0; - // SE(mean) = 1.7e-3 and SE(var) = 1.3e-3 at N, so 0.02 is ~11.5 and ~15 - // sigma. The previous 5e-3 put the mean check at 2.9 sigma. - // See ghosteau/fastdist#2. - assert(std::abs(mean - expected_mean) < 0.02); - assert(std::abs(var - expected_var) < 0.02); + assert(std::abs(mean - expected_mean) < 0.01); + assert(std::abs(var - expected_var) < 0.0075); } std::cout << "Continuous uniform distribution tests passed!\n"; diff --git a/tests/python/test_package_api.py b/tests/python/test_package_api.py index e84974d..bccc584 100644 --- a/tests/python/test_package_api.py +++ b/tests/python/test_package_api.py @@ -12,10 +12,21 @@ "Poisson", "Uniform", "Utils", ] +# Module-level functions, as opposed to the distribution classes. These live at +# the package root because they act on the one engine every sampler shares, so +# they belong to no single distribution. +EXPECTED_FUNCTIONS = ["seed", "seed_from_entropy"] + def test_top_level_all_matches_expected(): - # the top level exports every distribution class plus __version__ - assert sorted(fastdist.__all__) == sorted(EXPECTED + ["__version__"]) + # the top level exports every distribution class, the module-level + # functions, and __version__ + assert sorted(fastdist.__all__) == sorted(EXPECTED + EXPECTED_FUNCTIONS + ["__version__"]) + + +@pytest.mark.parametrize("name", EXPECTED_FUNCTIONS) +def test_every_exported_function_is_callable(name): + assert callable(getattr(fastdist, name)), f"{name} is not callable" def test_distributions_all_matches_top_level(): diff --git a/tests/python/test_rng.py b/tests/python/test_rng.py new file mode 100644 index 0000000..4b07ee4 --- /dev/null +++ b/tests/python/test_rng.py @@ -0,0 +1,84 @@ +""" +Tests for the seedable sampling engine exposed as fastdist.seed(). + +Every sampler in the library draws from one shared engine, so these tests mix +distributions deliberately: that pins down both the reproducibility guarantee +and the order in which samplers consume the stream. +""" + +import fastdist +from fastdist import Bernoulli, Normal, Uniform + + +def _trace(seed: int, n: int = 64) -> list: + """Replayable trace across several samplers sharing the one engine.""" + fastdist.seed(seed) + + uniform = Uniform(0.0, 1.0) + normal = Normal(0.0, 1.0) + bernoulli = Bernoulli(0.5) + + out = [] + for _ in range(n): + out.append(uniform.sample()) + out.append(normal.sample()) + out.append(bernoulli.sample()) + return out + + +def test_same_seed_reproduces_stream(): + """The core contract: same seed, same draws. + + Compared exactly rather than with a tolerance -- replaying a seeded stream + is bit-for-bit reproducible, not merely close. + """ + assert _trace(12345) == _trace(12345) + + +def test_distinct_seeds_give_distinct_streams(): + """Guards against seed() ignoring its argument. + + If it did, every seeded test in the suite would still pass while asserting + nothing, so this failure mode is worth pinning explicitly. + """ + assert _trace(1) != _trace(2) + + +def test_reseeding_rewinds_the_stream(): + """Re-seeding mid-stream returns the engine to the same point.""" + uniform = Uniform(0.0, 1.0) + + fastdist.seed(777) + first = uniform.sample() + + for _ in range(100): + uniform.sample() + + fastdist.seed(777) + assert uniform.sample() == first + + +def test_seed_from_entropy_escapes_the_fixed_stream(): + """seed_from_entropy() has to actually leave the seeded sequence. + + A single draw could collide by chance, so a whole trace is compared; the + odds of one repeating are nil. + """ + seeded = _trace(999) + + fastdist.seed_from_entropy() + uniform = Uniform(0.0, 1.0) + entropic = [uniform.sample() for _ in range(len(seeded))] + + assert entropic != seeded + + +def test_seed_is_exported_from_the_package_root(): + """seed() is part of the public API, not an implementation detail.""" + assert "seed" in fastdist.__all__ + assert "seed_from_entropy" in fastdist.__all__ + + +def teardown_module(module): + """Leave the engine non-deterministic for any test module that follows.""" + fastdist.seed_from_entropy() From 7482c0665c2e7a33311ef097ee669dca0fb71456 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:17:28 -0400 Subject: [PATCH 03/15] Add a benchmark suite and performance evidence log 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 --- BENCHMARKS.md | 170 +++++++ benchmarks/compare.py | 136 ++++++ benchmarks/harness.py | 250 ++++++++++ ...0.1.0_20260905T061458+0000_7671393bc4.json | 450 ++++++++++++++++++ benchmarks/run.py | 228 +++++++++ benchmarks/table.py | 88 ++++ 6 files changed, 1322 insertions(+) create mode 100644 BENCHMARKS.md create mode 100644 benchmarks/compare.py create mode 100644 benchmarks/harness.py create mode 100644 benchmarks/results/0.1.0_20260905T061458+0000_7671393bc4.json create mode 100644 benchmarks/run.py create mode 100644 benchmarks/table.py diff --git a/BENCHMARKS.md b/BENCHMARKS.md new file mode 100644 index 0000000..b3f8992 --- /dev/null +++ b/BENCHMARKS.md @@ -0,0 +1,170 @@ +# fastdist benchmarks + +Performance evidence log. One entry per release, plus entries for changes made +specifically to move performance. + +The point of this file is to be checkable. Every number here was produced by +`benchmarks/run.py`, saved as JSON under `benchmarks/results/`, and rendered by +`benchmarks/table.py` rather than typed in by hand. The raw reports are +committed alongside, so any claim below can be traced to the run that produced +it, on a named CPU, at a named commit. + +Results that are unflattering are recorded too. A log that only contains wins +is marketing, not evidence, and it would not catch a regression. + +--- + +## Running the suite + +```bash +pip install scipy # baseline, not needed to build or use fastdist +python benchmarks/run.py +``` + +The run writes `benchmarks/results/__.json` and +prints a summary. `--quick` uses one array size for a fast check; `--no-write` +prints without saving. + +To compare two runs: + +```bash +python benchmarks/compare.py --latest +``` + +`compare.py` exits non-zero if any case regressed by more than 5%, so it can +gate a change. To render a report for this log: + +```bash +python benchmarks/table.py --latest +``` + +--- + +## Method, and what the numbers do not say + +The baseline is SciPy, because that is the realistic alternative for someone +who would otherwise use this library. + +Each case is timed as several independent rounds, and the **minimum** round is +reported. Noise on a shared machine can only ever add time, so the minimum is +the best available estimate of the work itself; the mean would measure how busy +the machine was. Each report also records `noise_pct`, the gap between the +minimum and median, as a check on how quiet the run was. + +Every case with a baseline verifies that both implementations produce the same +numbers before either is timed. The `max abs diff` column carries that +agreement, and it is the reason the speedups can be taken at face value: at +1e-16 the two are computing the same function. + +Three caveats worth stating plainly: + +- **These are single-machine numbers.** Every figure below is one desktop CPU + on Windows. They are directionally useful, not a portable claim. +- **`scalar` is not a throughput number.** It measures the cost of one Python + call into each library. fastdist wins by ~60x there, but that mostly reflects + a thinner binding layer than SciPy's dispatch machinery, not faster math. + Quote the `batch` numbers instead. +- **The comparison is single-threaded on both sides.** Neither library is + parallelising these calls. + +--- + +## v0.1.0 — initial baseline + +First recorded measurement, taken at the point the benchmark suite was added so +that later work has something to be compared against. + + +- **Version** 0.1.0 (`7671393bc4` on `fix/flaky-rng-tolerances`, working tree dirty) +- **Measured** 2026-09-05T06:14:58+00:00 +- **CPU** AMD Ryzen 7 7700 8-Core Processor +- **Platform** Windows-11-10.0.26200-SP0 +- **Toolchain** Python 3.14.2, numpy 2.5.2, scipy 1.18.1 +- **CUDA** not built + +### batch (vs vectorised SciPy) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_pdf` | 1,000 | 5.70 us | 31.52 us (scipy) | **5.53x** | 1.1e-16 | +| `normal_cdf` | 1,000 | 7.76 us | 29.88 us (scipy) | **3.85x** | 2.2e-16 | +| `normal_logpdf` | 1,000 | 5.01 us | 32.29 us (scipy) | **6.44x** | 8.9e-16 | +| `exponential_pdf` | 1,000 | 4.49 us | 29.98 us (scipy) | **6.68x** | 0.0e+00 | +| `exponential_cdf` | 1,000 | 4.40 us | 31.13 us (scipy) | **7.07x** | 8.3e-17 | +| `uniform_pdf` | 1,000 | 2.37 us | 33.17 us (scipy) | **13.98x** | 0.0e+00 | +| `uniform_cdf` | 1,000 | 2.41 us | 31.55 us (scipy) | **13.07x** | 0.0e+00 | +| `poisson_pmf` | 1,000 | 47.22 us | 38.13 us (scipy) | **0.81x** | 2.0e-19 | +| `poisson_cdf` | 1,000 | 425.69 us | 70.40 us (scipy) | **0.17x** | 3.3e-16 | +| `bernoulli_pmf` | 1,000 | 1.92 us | 51.16 us (scipy) | **26.67x** | 2.2e-16 | +| `normal_pdf` | 100,000 | 462.80 us | 1.92 ms (scipy) | **4.16x** | 1.1e-16 | +| `normal_cdf` | 100,000 | 711.40 us | 2.32 ms (scipy) | **3.26x** | 2.2e-16 | +| `normal_logpdf` | 100,000 | 393.50 us | 2.03 ms (scipy) | **5.16x** | 8.9e-16 | +| `exponential_pdf` | 100,000 | 345.00 us | 1.81 ms (scipy) | **5.24x** | 0.0e+00 | +| `exponential_cdf` | 100,000 | 336.00 us | 1.99 ms (scipy) | **5.93x** | 1.1e-16 | +| `uniform_pdf` | 100,000 | 129.00 us | 1.95 ms (scipy) | **15.11x** | 0.0e+00 | +| `uniform_cdf` | 100,000 | 131.20 us | 1.86 ms (scipy) | **14.18x** | 0.0e+00 | +| `poisson_pmf` | 100,000 | 4.72 ms | 3.58 ms (scipy) | **0.76x** | 2.0e-19 | +| `poisson_cdf` | 100,000 | 43.47 ms | 6.45 ms (scipy) | **0.15x** | 3.3e-16 | +| `bernoulli_pmf` | 100,000 | 290.60 us | 3.59 ms (scipy) | **12.37x** | 2.2e-16 | +| `normal_pdf` | 1,000,000 | 5.34 ms | 19.98 ms (scipy) | **3.74x** | 1.1e-16 | +| `normal_cdf` | 1,000,000 | 7.67 ms | 21.87 ms (scipy) | **2.85x** | 2.2e-16 | +| `normal_logpdf` | 1,000,000 | 4.47 ms | 21.47 ms (scipy) | **4.80x** | 8.9e-16 | +| `exponential_pdf` | 1,000,000 | 4.45 ms | 16.78 ms (scipy) | **3.77x** | 0.0e+00 | +| `exponential_cdf` | 1,000,000 | 4.36 ms | 22.64 ms (scipy) | **5.20x** | 1.7e-16 | +| `uniform_pdf` | 1,000,000 | 1.88 ms | 20.21 ms (scipy) | **10.75x** | 0.0e+00 | +| `uniform_cdf` | 1,000,000 | 2.02 ms | 20.22 ms (scipy) | **10.03x** | 0.0e+00 | +| `poisson_pmf` | 1,000,000 | 50.15 ms | 38.98 ms (scipy) | **0.78x** | 2.0e-19 | +| `poisson_cdf` | 1,000,000 | 451.33 ms | 65.91 ms (scipy) | **0.15x** | 3.3e-16 | +| `bernoulli_pmf` | 1,000,000 | 3.66 ms | 40.59 ms (scipy) | **11.08x** | 2.2e-16 | + +### scalar (per-call cost, not throughput) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_pdf` | 20,000 | 8.55 ms | 488.55 ms (scipy) | **57.13x** | 1.1e-16 | +| `normal_cdf` | 20,000 | 7.15 ms | 477.79 ms (scipy) | **66.78x** | 2.2e-16 | + +### sample (vs numpy) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_sample` | 100,000 | 27.29 ms | 871.20 us (numpy) | **0.03x** | - | +| `uniform_sample` | 100,000 | 23.74 ms | 227.30 us (numpy) | **0.01x** | - | +| `normal_sample` | 1,000,000 | 303.10 ms | 10.19 ms (numpy) | **0.03x** | - | +| `uniform_sample` | 1,000,000 | 273.73 ms | 3.30 ms (numpy) | **0.01x** | - | + +### Reading of this baseline + +**Where the library is genuinely fast.** The continuous PDF/CDF batch paths beat +SciPy by 3–15x, and agreement to ~1e-16 confirms both sides compute the same +function. The largest margins are on the cheapest distributions — `uniform` at +10–15x, `bernoulli` at 11–27x — which is what one would expect: when the math +per element is trivial, the fixed overhead SciPy pays per call dominates, and +fastdist has less of it. The margin narrows as arrays grow (`uniform_pdf` falls +from 14x at n=1,000 to 10.75x at n=1,000,000), which is the same effect seen +from the other side — at a million elements the actual arithmetic starts to +dominate the fixed cost. + +**Where it is slower, and why.** `poisson_cdf` is 6.6x *slower* than SciPy, and +`poisson_pmf` about 25% slower. `poisson_cdf_scalar` sums the PMF from 0 to k, +calling `poisson_pmf_scalar` once per term, and each of those recomputes +`log(lambda)`, an `lgamma`, and an `exp`. For the benchmark's counts that is on +the order of twenty transcendental calls per element where a recurrence needs +none. This is a real defect, not a measurement artifact, and it is the clearest +optimisation target in the library. + +**Sampling is 30–100x slower than numpy.** This is a structural gap, not a +tuning problem: fastdist crosses the Python/C++ boundary once per variate, +while numpy fills an entire array per call. Closing it needs a batch sampling +entry point — `normal_sample_batch(n)` returning an array — which does not +exist yet. Until it does, the honest statement is that this library is for +evaluating distribution functions, not for bulk variate generation. + +--- + +## Changes to record here + +Add an entry when a release ships, or when a change is made specifically to +move performance. Each entry should carry the generated table, the commit, and +a short reading of what moved and why. `compare.py` output makes a good basis +for the reading. diff --git a/benchmarks/compare.py b/benchmarks/compare.py new file mode 100644 index 0000000..fd459ca --- /dev/null +++ b/benchmarks/compare.py @@ -0,0 +1,136 @@ +""" +Compare two benchmark reports and print a regression table. + + python benchmarks/compare.py BEFORE.json AFTER.json + python benchmarks/compare.py --latest # newest two reports + +Reads the JSON written by run.py. Cases are matched on (group, case, n), so +adding or removing benchmarks between runs is fine -- unmatched cases are +listed separately rather than silently dropped. + +Reading the output +------------------ +`change` is the change in fastdist's own time: negative is faster. It is the +number to look at when judging a code change. + +`speedup` columns are fastdist against the baseline library in each report. A +change there can come from either side, so a moved speedup with an unmoved +`change` means the baseline moved, not this library. + +The threshold for calling something a regression defaults to 5%, which is +comfortably above the noise on a quiet machine. Check the `noise_pct` field in +the reports before trusting a smaller difference: if either run was noisy, the +comparison is not meaningful at that resolution. +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +RESULTS_DIR = Path(__file__).resolve().parent / "results" + + +def load(path: Path) -> tuple[dict, dict]: + payload = json.loads(path.read_text(encoding="utf-8")) + index = {(r["group"], r["case"], r["n"]): r for r in payload["results"]} + return payload["environment"], index + + +def _fmt_time(seconds: float) -> str: + if seconds >= 1e-3: + return f"{seconds * 1e3:8.2f}ms" + return f"{seconds * 1e6:8.2f}us" + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("before", nargs="?", type=Path) + parser.add_argument("after", nargs="?", type=Path) + parser.add_argument("--latest", action="store_true", help="use the newest two reports") + parser.add_argument("--threshold", type=float, default=5.0, + help="percent change before a case is called out (default 5)") + args = parser.parse_args() + + if args.latest: + reports = sorted(RESULTS_DIR.glob("*.json")) + if len(reports) < 2: + return _fail(f"need two reports in {RESULTS_DIR}, found {len(reports)}") + before_path, after_path = reports[-2], reports[-1] + elif args.before and args.after: + before_path, after_path = args.before, args.after + else: + return _fail("pass two report paths, or --latest") + + before_env, before = load(before_path) + after_env, after = load(after_path) + + print(f"before {before_path.name}") + print(f" {before_env['fastdist_version']} @ {(before_env.get('git_commit') or '')[:10]}" + f" {before_env['timestamp_utc']}") + print(f"after {after_path.name}") + print(f" {after_env['fastdist_version']} @ {(after_env.get('git_commit') or '')[:10]}" + f" {after_env['timestamp_utc']}") + + if before_env.get("processor") != after_env.get("processor"): + print("\n!! different CPUs -- these reports are not comparable") + print(f" before: {before_env.get('processor')}") + print(f" after: {after_env.get('processor')}") + + print() + header = f"{'case':<34} {'before':>11} {'after':>11} {'change':>9} {'speedup':>16}" + print(header) + print("-" * len(header)) + + regressions, improvements = [], [] + + for key in sorted(before.keys() & after.keys()): + b, a = before[key], after[key] + group, case, n = key + + change = (a["fastdist_s"] - b["fastdist_s"]) / b["fastdist_s"] * 100.0 + speed = "" + if b.get("speedup") and a.get("speedup"): + speed = f"{b['speedup']:6.2f}x -> {a['speedup']:6.2f}x" + + label = f"{group}/{case} n={n:,}" + flag = "" + if change > args.threshold: + flag, _ = " REGRESSED", regressions.append((label, change)) + elif change < -args.threshold: + flag, _ = " faster", improvements.append((label, change)) + + print(f"{label:<34} {_fmt_time(b['fastdist_s']):>11} {_fmt_time(a['fastdist_s']):>11} " + f"{change:+8.1f}% {speed:>16}{flag}") + + only_before = before.keys() - after.keys() + only_after = after.keys() - before.keys() + for label, keys in (("only in before", only_before), ("only in after", only_after)): + if keys: + print(f"\n{label}:") + for group, case, n in sorted(keys): + print(f" {group}/{case} n={n:,}") + + print() + if regressions: + print(f"{len(regressions)} regression(s) beyond {args.threshold:g}%:") + for label, change in sorted(regressions, key=lambda t: -t[1]): + print(f" {label} {change:+.1f}%") + if improvements: + print(f"{len(improvements)} improvement(s) beyond {args.threshold:g}%:") + for label, change in sorted(improvements, key=lambda t: t[1]): + print(f" {label} {change:+.1f}%") + if not regressions and not improvements: + print(f"no case moved more than {args.threshold:g}%") + + return 1 if regressions else 0 + + +def _fail(message: str) -> int: + print(f"error: {message}") + return 2 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/benchmarks/harness.py b/benchmarks/harness.py new file mode 100644 index 0000000..ff2faa8 --- /dev/null +++ b/benchmarks/harness.py @@ -0,0 +1,250 @@ +""" +Timing harness for the fastdist benchmark suite. + +Separated from the benchmark definitions so the measurement methodology lives +in one place and can be reviewed on its own. + +Methodology +----------- +Each case is timed as `repeat` independent rounds of `inner` calls. The round +total is divided by `inner` to get a per-call time, and the reported figure is +the **minimum** round rather than the mean. + +Minimum is the right summary here. The quantity of interest is how long the +work takes; every source of noise on a shared machine (scheduler preemption, +frequency scaling, another process touching the cache) can only ever add time, +never remove it. The mean estimates "time under typical interference", which is +a property of the machine that day. The minimum estimates the work itself, +which is the thing a regression would move. + +The spread between the minimum and the median is reported alongside as +`noise_pct`. It is not an error bar on the measurement -- it is a measure of +how quiet the machine was. A large value means the numbers are still usable but +the machine was busy; comparisons across runs with very different noise_pct +deserve suspicion. + +A warmup round runs before timing and is discarded, so first-call costs (lazy +imports, page faults on the output buffer, branch predictor cold start) do not +land in the result. +""" + +from __future__ import annotations + +import json +import platform +import subprocess +import statistics +import sys +import time +from dataclasses import dataclass, asdict, field +from datetime import datetime, timezone +from pathlib import Path +from typing import Callable + +REPO_ROOT = Path(__file__).resolve().parents[1] + + +@dataclass +class Result: + """One measured case.""" + + group: str + case: str + n: int + # Seconds per call, minimum over rounds. + fastdist_s: float + baseline_s: float | None + baseline_name: str | None + # How quiet the machine was, as (median - min) / min, in percent. + fastdist_noise_pct: float + baseline_noise_pct: float | None + # Largest absolute difference between the two implementations' outputs. + # None when there is no baseline to compare against. + max_abs_diff: float | None + + @property + def speedup(self) -> float | None: + if self.baseline_s is None or self.fastdist_s == 0.0: + return None + return self.baseline_s / self.fastdist_s + + +def _time_one(fn: Callable[[], object], repeat: int, inner: int) -> tuple[float, float]: + """Return (min seconds per call, noise percent).""" + # Warmup, discarded. + for _ in range(inner): + fn() + + rounds = [] + for _ in range(repeat): + start = time.perf_counter() + for _ in range(inner): + fn() + rounds.append((time.perf_counter() - start) / inner) + + best = min(rounds) + median = statistics.median(rounds) + noise = ((median - best) / best * 100.0) if best > 0 else 0.0 + return best, noise + + +def measure( + group: str, + case: str, + n: int, + fastdist_fn: Callable[[], object], + baseline_fn: Callable[[], object] | None = None, + baseline_name: str | None = None, + repeat: int = 7, + inner: int = 1, +) -> Result: + """Time one case against an optional baseline and check they agree.""" + fd_s, fd_noise = _time_one(fastdist_fn, repeat, inner) + + base_s = base_noise = max_abs_diff = None + if baseline_fn is not None: + base_s, base_noise = _time_one(baseline_fn, repeat, inner) + + # A speedup only means something if both sides computed the same thing. + # Any case whose outputs disagree is a bug in the benchmark or the + # library, and reporting its timing would be misleading either way. + max_abs_diff = _max_abs_diff(fastdist_fn(), baseline_fn()) + + return Result( + group=group, + case=case, + n=n, + fastdist_s=fd_s, + baseline_s=base_s, + baseline_name=baseline_name, + fastdist_noise_pct=fd_noise, + baseline_noise_pct=base_noise, + max_abs_diff=max_abs_diff, + ) + + +def _max_abs_diff(a, b) -> float: + """Largest absolute elementwise difference, ignoring matching non-finites.""" + import numpy as np + + a = np.asarray(a, dtype=float) + b = np.asarray(b, dtype=float) + if a.shape != b.shape: + return float("inf") + + finite = np.isfinite(a) & np.isfinite(b) + # Disagreeing on *where* the non-finites are is a real mismatch. + if not np.array_equal(np.isfinite(a), np.isfinite(b)): + return float("inf") + if not finite.any(): + return 0.0 + return float(np.max(np.abs(a[finite] - b[finite]))) + + +def environment() -> dict: + """Everything needed to judge whether two runs are comparable.""" + import numpy + + try: + import scipy + + scipy_version = scipy.__version__ + except ImportError: + scipy_version = None + + sys.path.insert(0, str(REPO_ROOT / "python")) + import fastdist + import fastdist._fastdist as core + + return { + "timestamp_utc": datetime.now(timezone.utc).isoformat(timespec="seconds"), + "fastdist_version": core.__version__, + "git_commit": _git("rev-parse", "HEAD"), + "git_branch": _git("rev-parse", "--abbrev-ref", "HEAD"), + "git_dirty": bool(_git("status", "--porcelain")), + "cuda_available": hasattr(core, "normal_pdf_cuda"), + "python": platform.python_version(), + "numpy": numpy.__version__, + "scipy": scipy_version, + "platform": platform.platform(), + "processor": _cpu_name(), + "machine": platform.machine(), + } + + +def _git(*args: str) -> str | None: + try: + out = subprocess.run( + ["git", *args], + cwd=REPO_ROOT, + capture_output=True, + text=True, + check=True, + ) + return out.stdout.strip() + except (OSError, subprocess.CalledProcessError): + return None + + +def _cpu_name() -> str: + """platform.processor() is often empty or useless; try harder on each OS. + + Every probe here is best-effort. A missing CPU name is cosmetic -- it must + never take the benchmark run down with it -- so each branch swallows its + own failures and the function always returns a string. + """ + if sys.platform == "win32": + # wmic was removed in recent Windows 11 builds, so try the registry + # first and only then fall back to the (vaguer) environment variable. + try: + import winreg + + key = winreg.OpenKey( + winreg.HKEY_LOCAL_MACHINE, + r"HARDWARE\DESCRIPTION\System\CentralProcessor\0", + ) + with key: + return str(winreg.QueryValueEx(key, "ProcessorNameString")[0]).strip() + except OSError: + pass + + import os + + return os.environ.get("PROCESSOR_IDENTIFIER") or platform.processor() or "unknown" + + if sys.platform == "darwin": + try: + return subprocess.run( + ["sysctl", "-n", "machdep.cpu.brand_string"], + capture_output=True, + text=True, + check=True, + ).stdout.strip() + except (OSError, subprocess.CalledProcessError): + return platform.processor() or "unknown" + + try: + for line in Path("/proc/cpuinfo").read_text().splitlines(): + if line.startswith("model name"): + return line.split(":", 1)[1].strip() + except OSError: + pass + return platform.processor() or "unknown" + + +def write_report(results: list[Result], env: dict, out_dir: Path) -> Path: + """Write one run to a JSON file named for the version and commit.""" + out_dir.mkdir(parents=True, exist_ok=True) + + commit = (env.get("git_commit") or "nocommit")[:10] + stamp = env["timestamp_utc"].replace(":", "").replace("-", "") + path = out_dir / f"{env['fastdist_version']}_{stamp}_{commit}.json" + + payload = { + "environment": env, + "results": [ + {**asdict(r), "speedup": r.speedup} for r in results + ], + } + path.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8") + return path diff --git a/benchmarks/results/0.1.0_20260905T061458+0000_7671393bc4.json b/benchmarks/results/0.1.0_20260905T061458+0000_7671393bc4.json new file mode 100644 index 0000000..3a6395c --- /dev/null +++ b/benchmarks/results/0.1.0_20260905T061458+0000_7671393bc4.json @@ -0,0 +1,450 @@ +{ + "environment": { + "timestamp_utc": "2026-09-05T06:14:58+00:00", + "fastdist_version": "0.1.0", + "git_commit": "7671393bc4524020842382b5fe320156de64b616", + "git_branch": "fix/flaky-rng-tolerances", + "git_dirty": true, + "cuda_available": false, + "python": "3.14.2", + "numpy": "2.5.2", + "scipy": "1.18.1", + "platform": "Windows-11-10.0.26200-SP0", + "processor": "AMD Ryzen 7 7700 8-Core Processor", + "machine": "AMD64" + }, + "results": [ + { + "group": "batch", + "case": "normal_pdf", + "n": 1000, + "fastdist_s": 5.6979997316375374e-06, + "baseline_s": 3.1522000208497046e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.492808268471156, + "baseline_noise_pct": 0.9073023524743741, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 5.5321168292575535 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000, + "fastdist_s": 7.758000283502042e-06, + "baseline_s": 2.98839999595657e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.12889267469594837, + "baseline_noise_pct": 0.5688667197854018, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.8520235714757867 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000, + "fastdist_s": 5.013999762013555e-06, + "baseline_s": 3.229000023566186e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.3988861834695411, + "baseline_noise_pct": 5.500155226436825, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 6.439968442019757 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000, + "fastdist_s": 4.489999846555293e-06, + "baseline_s": 2.9980000108480452e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 6.77061423085315, + "baseline_noise_pct": 8.972646826203025, + "max_abs_diff": 0.0, + "speedup": 6.677060385977735 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000, + "fastdist_s": 4.400000325404107e-06, + "baseline_s": 3.1128000118769704e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.4545299645740266, + "baseline_noise_pct": 9.348496620924507, + "max_abs_diff": 8.326672684688674e-17, + "speedup": 7.074544958337208 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000, + "fastdist_s": 2.3720000172033907e-06, + "baseline_s": 3.317199996672571e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.674516848197502, + "baseline_noise_pct": 8.543349950206235, + "max_abs_diff": 0.0, + "speedup": 13.984822818777124 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000, + "fastdist_s": 2.414000337012112e-06, + "baseline_s": 3.154599980916828e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.242735798261051, + "baseline_noise_pct": 1.0904719874476938, + "max_abs_diff": 0.0, + "speedup": 13.06793513053681 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000, + "fastdist_s": 4.721800039988011e-05, + "baseline_s": 3.813199989963323e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.532931897961865, + "baseline_noise_pct": 2.5805096198875774, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.8075733740671079 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000, + "fastdist_s": 0.00042568799981381746, + "baseline_s": 7.040199998300522e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.9171036164473829, + "baseline_noise_pct": 2.72435449723879, + "max_abs_diff": 3.3306690738754696e-16, + "speedup": 0.1653840371675895 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000, + "fastdist_s": 1.9180000526830553e-06, + "baseline_s": 5.115600011777133e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.10427607789629564, + "baseline_noise_pct": 4.718898563957289, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 26.671532175512787 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 100000, + "fastdist_s": 0.0004628000024240464, + "baseline_s": 0.0019247000163886696, + "baseline_name": "scipy", + "fastdist_noise_pct": 18.60415187394897, + "baseline_noise_pct": 24.81945063091742, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 4.158815916826938 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 100000, + "fastdist_s": 0.0007114000036381185, + "baseline_s": 0.0023206000041682273, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.6887797411158246, + "baseline_noise_pct": 4.942685984556709, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.2620185441391865 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 100000, + "fastdist_s": 0.0003934999986086041, + "baseline_s": 0.0020284999918658286, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.2198166096263638, + "baseline_noise_pct": 4.426917168030424, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 5.155019057277004 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 100000, + "fastdist_s": 0.0003449999785516411, + "baseline_s": 0.001806199987186119, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.797103772202533, + "baseline_noise_pct": 11.582328802420228, + "max_abs_diff": 0.0, + "speedup": 5.235362607176971 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 100000, + "fastdist_s": 0.00033599999733269215, + "baseline_s": 0.0019929000118281692, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.6666707088653734, + "baseline_noise_pct": 13.558130837496488, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 5.931250082287616 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 100000, + "fastdist_s": 0.0001289999927394092, + "baseline_s": 0.001948899996932596, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.387599872394361, + "baseline_noise_pct": 8.95889942989107, + "max_abs_diff": 0.0, + "speedup": 15.107752764525632 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 100000, + "fastdist_s": 0.00013119998038746417, + "baseline_s": 0.001860599993960932, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.5335407277334036, + "baseline_noise_pct": 8.991722771107007, + "max_abs_diff": 0.0, + "speedup": 14.181404512913385 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 100000, + "fastdist_s": 0.004719599994132295, + "baseline_s": 0.003583799989428371, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.409272252477938, + "baseline_noise_pct": 4.631955271802637, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.7593440109085469 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 100000, + "fastdist_s": 0.043472599994856864, + "baseline_s": 0.006450299988500774, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.7701403121477435, + "baseline_noise_pct": 2.599879475137506, + "max_abs_diff": 3.3306690738754696e-16, + "speedup": 0.14837621833669698 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 100000, + "fastdist_s": 0.0002906000008806586, + "baseline_s": 0.0035945000126957893, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.9291092385132008, + "baseline_noise_pct": 4.721101989472907, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 12.369236069520698 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 1000000, + "fastdist_s": 0.005337300011888146, + "baseline_s": 0.019984699989436194, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.5893237249814374, + "baseline_noise_pct": 1.7533414054507943, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 3.744346382051385 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000000, + "fastdist_s": 0.007669299986446276, + "baseline_s": 0.021866399998543784, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.4174307312383236, + "baseline_noise_pct": 7.199173142639623, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 2.851159823867578 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000000, + "fastdist_s": 0.004473199980566278, + "baseline_s": 0.021465000027092174, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.2384873162387213, + "baseline_noise_pct": 3.8616350615994133, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 4.798578225956007 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000000, + "fastdist_s": 0.004447999992407858, + "baseline_s": 0.01677919999929145, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.4572839127608495, + "baseline_noise_pct": 18.562863529680822, + "max_abs_diff": 0.0, + "speedup": 3.7723021645529005 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000000, + "fastdist_s": 0.00435590001870878, + "baseline_s": 0.022635200002696365, + "baseline_name": "scipy", + "fastdist_noise_pct": 14.398860648416289, + "baseline_noise_pct": 3.1654236948603165, + "max_abs_diff": 1.6653345369377348e-16, + "speedup": 5.1964461777077515 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000000, + "fastdist_s": 0.001879399991594255, + "baseline_s": 0.020212400006130338, + "baseline_name": "scipy", + "fastdist_noise_pct": 12.317761921931861, + "baseline_noise_pct": 8.196453579813765, + "max_abs_diff": 0.0, + "speedup": 10.754709001027818 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000000, + "fastdist_s": 0.002015199977904558, + "baseline_s": 0.02021650000824593, + "baseline_name": "scipy", + "fastdist_noise_pct": 16.191941959174223, + "baseline_noise_pct": 3.426409009372913, + "max_abs_diff": 0.0, + "speedup": 10.032006862796523 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000000, + "fastdist_s": 0.05014619999565184, + "baseline_s": 0.03898069998831488, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.2479510337370519, + "baseline_noise_pct": 7.372366392640059, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.7773410545902757 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000000, + "fastdist_s": 0.4513342999853194, + "baseline_s": 0.06590710001182742, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.1765019924326166, + "baseline_noise_pct": 1.7460940827250089, + "max_abs_diff": 3.3306690738754696e-16, + "speedup": 0.14602723527542932 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000000, + "fastdist_s": 0.003662899980554357, + "baseline_s": 0.040586399991298094, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.861476047858873, + "baseline_noise_pct": 6.957010213170335, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 11.080400831790005 + }, + { + "group": "scalar", + "case": "normal_pdf", + "n": 20000, + "fastdist_s": 0.008551600010832772, + "baseline_s": 0.48854549997486174, + "baseline_name": "scipy", + "fastdist_noise_pct": 11.60367636008294, + "baseline_noise_pct": 7.529718325774552, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 57.129133654052445 + }, + { + "group": "scalar", + "case": "normal_cdf", + "n": 20000, + "fastdist_s": 0.007154999999329448, + "baseline_s": 0.4777868000091985, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.894479460155713, + "baseline_noise_pct": 6.3198899538978965, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 66.77663173360946 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 100000, + "fastdist_s": 0.027288599987514317, + "baseline_s": 0.0008711999980732799, + "baseline_name": "numpy", + "fastdist_noise_pct": 3.617627905760217, + "baseline_noise_pct": 7.64462953482912, + "max_abs_diff": null, + "speedup": 0.031925419349907676 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 100000, + "fastdist_s": 0.02374400000553578, + "baseline_s": 0.0002272999845445156, + "baseline_name": "numpy", + "fastdist_noise_pct": 2.0763140137978535, + "baseline_noise_pct": 0.17599303863985102, + "max_abs_diff": null, + "speedup": 0.00957294409078175 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 1000000, + "fastdist_s": 0.30309840000700206, + "baseline_s": 0.010189299995545298, + "baseline_name": "numpy", + "fastdist_noise_pct": 1.4449762836444306, + "baseline_noise_pct": 0.7861187596350936, + "max_abs_diff": null, + "speedup": 0.03361713554182374 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 1000000, + "fastdist_s": 0.2737265000178013, + "baseline_s": 0.0033044999872799963, + "baseline_name": "numpy", + "fastdist_noise_pct": 5.09954279220077, + "baseline_noise_pct": 5.843547160249549, + "max_abs_diff": null, + "speedup": 0.012072269170376614 + } + ] +} diff --git a/benchmarks/run.py b/benchmarks/run.py new file mode 100644 index 0000000..08824ec --- /dev/null +++ b/benchmarks/run.py @@ -0,0 +1,228 @@ +""" +The fastdist benchmark suite. + +Run from the repo root: + + python benchmarks/run.py # full suite, writes a JSON report + python benchmarks/run.py --quick # fewer sizes, for a fast check + python benchmarks/run.py --no-write # print only, write nothing + +Results land in benchmarks/results/ as one JSON file per run, tagged with the +version, commit and machine. benchmarks/compare.py turns two of those into a +regression table, and BENCHMARKS.md is the curated log across releases. + +What is being compared +---------------------- +The baseline is SciPy, because that is what someone reaching for this library +would otherwise use. Comparisons are grouped by how fair they are: + + batch fastdist's *_cpu entry points against the equivalent vectorised SciPy + call. Both take a numpy array and return one, both do the loop in + compiled code, and neither pays per-element Python overhead. This is + the honest headline comparison. + + scalar fastdist's *_scalar entry points against SciPy called on one value at + a time. Both sides pay Python call overhead per element, so this + measures the cost of a single call rather than throughput. It is + reported because users do write scalar loops, but it flatters + whichever library has the thinner binding layer and should not be + quoted as a throughput number. + + sample Drawing variates. fastdist samples one value per call, while numpy + fills an array in one call, so numpy is expected to win by a wide + margin. It is measured anyway: this is a real gap in the library and + the log should record it rather than quietly omit it. + +Every case with a baseline also checks that the two implementations agree +numerically (see harness.max_abs_diff). A speedup on a wrong answer is not a +speedup. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np + +REPO_ROOT = Path(__file__).resolve().parents[1] +sys.path.insert(0, str(REPO_ROOT / "python")) +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from harness import Result, environment, measure, write_report # noqa: E402 + +try: + from scipy import stats as sps +except ImportError: # pragma: no cover + sys.exit("benchmarks require scipy: pip install scipy") + +import fastdist._fastdist as core # noqa: E402 + +SIZES = (1_000, 100_000, 1_000_000) +QUICK_SIZES = (100_000,) +SCALAR_N = 20_000 + + +# --------------------------------------------------------------------------- +# Batch: fastdist *_cpu vs vectorised scipy +# --------------------------------------------------------------------------- +def batch_cases(sizes): + """(case name, fastdist callable factory, scipy callable factory). + + The trailing 0.0 on every *_cpu call is step_size. The C++ headers declare + a default for it but the pybind11 bindings do not expose one, so it has to + be passed explicitly; 0.0 means "evaluate x as given". + """ + rng = np.random.default_rng(20260905) + + for n in sizes: + x_real = rng.normal(0.0, 1.0, n) + x_pos = np.abs(rng.normal(2.0, 1.0, n)) + 0.05 + x_unit = rng.uniform(0.01, 0.99, n) + k_count = rng.integers(0, 20, n).astype(float) + # bernoulli_pmf_batch takes int32, and poisson takes an int step_size: + # the discrete batch entry points are not typed uniformly with the + # continuous ones. + k_binary = rng.integers(0, 2, n).astype(np.int32) + + yield from [ + ("normal_pdf", n, + lambda x=x_real: core.normal_pdf_cpu(x, 0.0, 1.0, 0.0), + lambda x=x_real: sps.norm.pdf(x, 0.0, 1.0)), + ("normal_cdf", n, + lambda x=x_real: core.normal_cdf_cpu(x, 0.0, 1.0, 0.0), + lambda x=x_real: sps.norm.cdf(x, 0.0, 1.0)), + ("normal_logpdf", n, + lambda x=x_real: core.normal_logpdf_cpu(x, 0.0, 1.0, 0.0), + lambda x=x_real: sps.norm.logpdf(x, 0.0, 1.0)), + ("exponential_pdf", n, + lambda x=x_pos: core.exponential_pdf_cpu(x, 2.0, 0.0), + lambda x=x_pos: sps.expon.pdf(x, scale=0.5)), + ("exponential_cdf", n, + lambda x=x_pos: core.exponential_cdf_cpu(x, 2.0, 0.0), + lambda x=x_pos: sps.expon.cdf(x, scale=0.5)), + ("uniform_pdf", n, + lambda x=x_real: core.uniform_pdf_cpu(x, -3.0, 3.0, 0.0), + lambda x=x_real: sps.uniform.pdf(x, loc=-3.0, scale=6.0)), + ("uniform_cdf", n, + lambda x=x_real: core.uniform_cdf_cpu(x, -3.0, 3.0, 0.0), + lambda x=x_real: sps.uniform.cdf(x, loc=-3.0, scale=6.0)), + ("poisson_pmf", n, + lambda x=k_count: core.poisson_pmf_cpu(x, 4.0, 0), + lambda x=k_count: sps.poisson.pmf(x, 4.0)), + ("poisson_cdf", n, + lambda x=k_count: core.poisson_cdf_cpu(x, 4.0, 0), + lambda x=k_count: sps.poisson.cdf(x, 4.0)), + ("bernoulli_pmf", n, + lambda x=k_binary: core.bernoulli_pmf_cpu(x, 0.3, 0), + lambda x=k_binary: sps.bernoulli.pmf(x, 0.3)), + ] + + +# --------------------------------------------------------------------------- +# Scalar: per-call cost +# --------------------------------------------------------------------------- +def scalar_cases(): + xs = np.random.default_rng(7).normal(0.0, 1.0, SCALAR_N) + + def fd_loop(fn, *args): + return lambda: [fn(float(v), *args) for v in xs] + + def sp_loop(fn, *args, **kw): + return lambda: [float(fn(float(v), *args, **kw)) for v in xs] + + return [ + ("normal_pdf", SCALAR_N, + fd_loop(core.normal_pdf_scalar, 0.0, 1.0), + sp_loop(sps.norm.pdf, 0.0, 1.0)), + ("normal_cdf", SCALAR_N, + fd_loop(core.normal_cdf_scalar, 0.0, 1.0), + sp_loop(sps.norm.cdf, 0.0, 1.0)), + ] + + +# --------------------------------------------------------------------------- +# Sampling +# --------------------------------------------------------------------------- +def sample_cases(sizes): + for n in sizes: + rng = np.random.default_rng(11) + yield ("normal_sample", n, + lambda n=n: [core.normal_sample(0.0, 1.0) for _ in range(n)], + lambda n=n, rng=rng: rng.normal(0.0, 1.0, n)) + yield ("uniform_sample", n, + lambda n=n: [core.uniform_sample(0.0, 1.0) for _ in range(n)], + lambda n=n, rng=rng: rng.uniform(0.0, 1.0, n)) + + +def run(sizes, sample_sizes) -> list[Result]: + results: list[Result] = [] + + for case, n, fd, sp in batch_cases(sizes): + # Big arrays are slow enough that one call per round is plenty; small + # ones need repetition to rise above timer resolution. + inner = 50 if n <= 1_000 else 1 + results.append(measure("batch", case, n, fd, sp, "scipy", inner=inner)) + print(f" batch {case:<18} n={n:<9,} {_fmt(results[-1])}") + + for case, n, fd, sp in scalar_cases(): + results.append(measure("scalar", case, n, fd, sp, "scipy")) + print(f" scalar {case:<18} n={n:<9,} {_fmt(results[-1])}") + + for case, n, fd, np_fn in sample_cases(sample_sizes): + # No correctness check: both draw from their own RNG, so the outputs + # are different random numbers by construction. measure() would flag + # that as a mismatch, so the baseline is timed as a separate case. + fd_result = measure("sample", case, n, fd) + np_result = measure("sample", case, n, np_fn) + fd_result.baseline_s = np_result.fastdist_s + fd_result.baseline_name = "numpy" + fd_result.baseline_noise_pct = np_result.fastdist_noise_pct + results.append(fd_result) + print(f" sample {case:<18} n={n:<9,} {_fmt(fd_result)}") + + return results + + +def _fmt(r: Result) -> str: + speed = r.speedup + verdict = "-" if speed is None else f"{speed:6.2f}x" + diff = "" if r.max_abs_diff is None else f" maxdiff={r.max_abs_diff:.2e}" + return f"fastdist={r.fastdist_s * 1e6:10.2f}us {verdict}{diff}" + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--quick", action="store_true", help="one array size only") + parser.add_argument("--no-write", action="store_true", help="do not save a report") + args = parser.parse_args() + + env = environment() + print("fastdist benchmark suite") + print(f" version {env['fastdist_version']} ({env['git_branch']}@{(env['git_commit'] or '')[:10]}" + f"{', dirty' if env['git_dirty'] else ''})") + print(f" cpu {env['processor']}") + print(f" python {env['python']} numpy {env['numpy']} scipy {env['scipy']}") + print(f" cuda {'available' if env['cuda_available'] else 'not built'}") + print() + + sizes = QUICK_SIZES if args.quick else SIZES + sample_sizes = (100_000,) if args.quick else (100_000, 1_000_000) + results = run(sizes, sample_sizes) + + mismatched = [r for r in results if r.max_abs_diff is not None and r.max_abs_diff > 1e-9] + if mismatched: + print("\nWARNING: outputs disagree with the baseline, timings below are not comparable:") + for r in mismatched: + print(f" {r.group}/{r.case} n={r.n}: max abs diff {r.max_abs_diff:.3e}") + + if not args.no_write: + path = write_report(results, env, REPO_ROOT / "benchmarks" / "results") + print(f"\nwrote {path.relative_to(REPO_ROOT)}") + + return 1 if mismatched else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/benchmarks/table.py b/benchmarks/table.py new file mode 100644 index 0000000..e524981 --- /dev/null +++ b/benchmarks/table.py @@ -0,0 +1,88 @@ +""" +Render one benchmark report as a Markdown table, for pasting into BENCHMARKS.md. + + python benchmarks/table.py benchmarks/results/.json + python benchmarks/table.py --latest + +Exists so entries in the evidence log are generated from the recorded JSON +rather than transcribed by hand -- a typo in a number nobody can check later is +worse than no number at all. +""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path + +RESULTS_DIR = Path(__file__).resolve().parent / "results" + + +def _fmt_time(seconds: float) -> str: + if seconds >= 1.0: + return f"{seconds:.2f} s" + if seconds >= 1e-3: + return f"{seconds * 1e3:.2f} ms" + return f"{seconds * 1e6:.2f} us" + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("report", nargs="?", type=Path) + parser.add_argument("--latest", action="store_true") + parser.add_argument("--group", help="only this group (batch, scalar, sample)") + args = parser.parse_args() + + if args.latest: + reports = sorted(RESULTS_DIR.glob("*.json")) + if not reports: + print(f"error: no reports in {RESULTS_DIR}") + return 2 + path = reports[-1] + elif args.report: + path = args.report + else: + print("error: pass a report path, or --latest") + return 2 + + payload = json.loads(path.read_text(encoding="utf-8")) + env = payload["environment"] + + print(f"") + print() + print(f"- **Version** {env['fastdist_version']} " + f"(`{(env.get('git_commit') or '')[:10]}` on `{env.get('git_branch')}`" + f"{', working tree dirty' if env.get('git_dirty') else ''})") + print(f"- **Measured** {env['timestamp_utc']}") + print(f"- **CPU** {env['processor']}") + print(f"- **Platform** {env['platform']}") + print(f"- **Toolchain** Python {env['python']}, numpy {env['numpy']}, scipy {env['scipy']}") + print(f"- **CUDA** {'available' if env.get('cuda_available') else 'not built'}") + print() + + groups = {} + for r in payload["results"]: + if args.group and r["group"] != args.group: + continue + groups.setdefault(r["group"], []).append(r) + + for group, rows in groups.items(): + print(f"### {group}") + print() + print("| case | n | fastdist | baseline | speedup | max abs diff |") + print("|---|---:|---:|---:|---:|---:|") + for r in rows: + base = _fmt_time(r["baseline_s"]) if r["baseline_s"] is not None else "-" + if r["baseline_name"]: + base += f" ({r['baseline_name']})" + speed = f"{r['speedup']:.2f}x" if r.get("speedup") else "-" + diff = "-" if r["max_abs_diff"] is None else f"{r['max_abs_diff']:.1e}" + print(f"| `{r['case']}` | {r['n']:,} | {_fmt_time(r['fastdist_s'])} " + f"| {base} | **{speed}** | {diff} |") + print() + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From 73ba7d69055948cba21190122d6c409a29c4a903 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:24:13 -0400 Subject: [PATCH 04/15] Remove quadratic work from the discrete CDFs and hoist batch invariants 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 --- src/math/beta.cpp | 7 ++ src/math/binomial.cpp | 29 ++++++++- src/math/exponential.cpp | 25 +++++++- src/math/negative_binomial.cpp | 39 ++++++++++-- src/math/normal.cpp | 88 +++++++++++++++++++++++--- src/math/poisson.cpp | 51 +++++++++++++-- src/math/uniform.cpp | 40 +++++++++++- tests/python/test_negative_binomial.py | 23 ++----- 8 files changed, 260 insertions(+), 42 deletions(-) diff --git a/src/math/beta.cpp b/src/math/beta.cpp index 3e7f486..3f50989 100644 --- a/src/math/beta.cpp +++ b/src/math/beta.cpp @@ -64,6 +64,13 @@ namespace fastdist::math { // RNG // ------------------------- double beta_sample(const double alpha, const double beta) { + // Every other sampler validates its parameters; this one did not, and + // std::gamma_distribution has undefined behaviour for a non-positive + // shape rather than a defined error value. + if (!std::isfinite(alpha) || !std::isfinite(beta) || alpha <= 0.0 || beta <= 0.0) { + return std::numeric_limits::quiet_NaN(); + } + std::gamma_distribution ga(alpha, 1.0); std::gamma_distribution gb(beta, 1.0); double a = ga(rng()); diff --git a/src/math/binomial.cpp b/src/math/binomial.cpp index 100693e..45f8af8 100644 --- a/src/math/binomial.cpp +++ b/src/math/binomial.cpp @@ -36,12 +36,37 @@ namespace fastdist::math { if (x < 0) return 0.0; if (x >= n) return 1.0; + // p == 1 puts all mass at n, and x < n here, so nothing has accumulated + // yet. Handled separately because the ratio below divides by (1 - p). + if (p == 1.0) return 0.0; + + // Consecutive PMF terms satisfy + // P(k) = P(k-1) * ((n - k + 1) / k) * (p / (1 - p)) + // so the sum costs one exp overall instead of three lgammas, two logs + // and an exp per term. + // + // P(0) = (1-p)^n underflows for large n, which would collapse the whole + // recurrence to zero; fall back to per-term log-space evaluation there. + const double log_p0 = static_cast(n) * std::log1p(-p); + constexpr double MIN_RECURRENCE_LOG = -700.0; + + if (log_p0 > MIN_RECURRENCE_LOG) { + const double odds = p / (1.0 - p); + double term = std::exp(log_p0); + double sum = term; + for (int k = 1; k <= x; ++k) { + term *= (static_cast(n - k + 1) / static_cast(k)) * odds; + sum += term; + } + // Summing PMF terms accumulates rounding error, so the total can + // land a few ULP above 1.0 + return std::min(sum, 1.0); + } + double sum = 0.0; for (int k = 0; k <= x; ++k) { sum += binomial_pmf_scalar(k, n, p); } - // Summing PMF terms accumulates rounding error, so the total can - // land a few ULP above 1.0 return std::min(sum, 1.0); } diff --git a/src/math/exponential.cpp b/src/math/exponential.cpp index aba4e0d..87e9dfd 100644 --- a/src/math/exponential.cpp +++ b/src/math/exponential.cpp @@ -1,5 +1,6 @@ // Function declarations for exponential distribution functions #include "fastdist/math/exponential.h" +#include #include #include #include @@ -79,15 +80,35 @@ namespace fastdist::math { // Batch Functions void exponential_pdf_batch(const double* x_data, double* output, const size_t n, const double lambda, const double stepSize) { + if (!std::isfinite(lambda) || lambda <= 0.0) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + for (size_t i = 0; i < n; i++) { - output[i] = exponential_pdf_scalar(x_data[i] + stepSize * static_cast(i), lambda); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + output[i] = (x < 0.0) ? 0.0 : lambda * std::exp(-lambda * x); } } void exponential_cdf_batch(const double* x_data, double* output, const size_t n, const double lambda, const double stepSize) { + if (!std::isfinite(lambda) || lambda <= 0.0) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + for (size_t i = 0; i < n; i++) { - output[i] = exponential_cdf_scalar(x_data[i] + stepSize * static_cast(i), lambda); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + output[i] = (x < 0.0) ? 0.0 : 1.0 - std::exp(-lambda * x); } } diff --git a/src/math/negative_binomial.cpp b/src/math/negative_binomial.cpp index 3398e10..47c2a2b 100644 --- a/src/math/negative_binomial.cpp +++ b/src/math/negative_binomial.cpp @@ -23,9 +23,17 @@ namespace fastdist::math { return 0.0; } - const double comb = std::tgamma(k + r) / (std::tgamma(r) * std::tgamma(k + 1.0)); - - return comb * std::pow(p, r) * std::pow(1.0 - p, k); + // Evaluated in log space. Forming C(k + r - 1, k) from raw factorials + // overflows a double once k + r - 1 > 170 -- inf at k = 170 and nan + // beyond -- even though the coefficient and the resulting PMF are + // comfortably inside range (for r = 3, k = 200 the true PMF is 1.6e-57). + // lgamma keeps the intermediate values small, and folding the two pows + // into the same exponent removes them from the hot path. + const double log_pmf = std::lgamma(static_cast(k) + r) - std::lgamma(static_cast(r)) - + std::lgamma(static_cast(k) + 1.0) + r * std::log(p) + + static_cast(k) * std::log1p(-p); + + return std::exp(log_pmf); } // ------------------------- @@ -40,12 +48,33 @@ namespace fastdist::math { return 0.0; } + // Consecutive PMF terms satisfy + // P(i) = P(i-1) * ((i + r - 1) / i) * (1 - p) + // so the sum costs one exp overall instead of three tgammas and two + // pows per term. + // + // P(0) = p^r underflows for small p with large r, which would collapse + // the recurrence to zero; fall back to per-term evaluation there. + const double log_p0 = static_cast(r) * std::log(p); + constexpr double MIN_RECURRENCE_LOG = -700.0; + + if (log_p0 > MIN_RECURRENCE_LOG) { + const double q = 1.0 - p; + double term = std::exp(log_p0); + double sum = term; + for (int i = 1; i <= k; ++i) { + term *= (static_cast(i + r - 1) / static_cast(i)) * q; + sum += term; + } + // Summing PMF terms accumulates rounding error, so the total can + // land a few ULP above 1.0 + return std::min(sum, 1.0); + } + double sum = 0.0; for (int i = 0; i <= k; ++i) { sum += negative_binomial_pmf_scalar(i, r, p); } - // Summing PMF terms accumulates rounding error, so the total can - // land a few ULP above 1.0 return std::min(sum, 1.0); } diff --git a/src/math/normal.cpp b/src/math/normal.cpp index da84fff..f145b0b 100644 --- a/src/math/normal.cpp +++ b/src/math/normal.cpp @@ -1,4 +1,5 @@ // Function declarations for normal distribution functions +#include #include #include #include @@ -8,13 +9,38 @@ namespace fastdist::math { + namespace { + // The scalar formulas with parameter validation and every loop-invariant + // term lifted into arguments, so the batch paths can compute those once + // instead of once per element. Scalar and batch both route through these, + // so there is still only one copy of each formula. + // + // The arithmetic is arranged exactly as the scalar versions had it -- + // same operations in the same order -- so hoisting does not perturb + // rounding and the results are bit-identical to before. + inline double normal_pdf_core(const double x, const double mu, const double sigma, const double denom) { + const double z = (x - mu) / sigma; + return std::exp(-0.5 * z * z) / denom; + } + + inline double normal_logpdf_core(const double x, const double mu, const double inv_sigma, + const double log_sigma) { + const double z = (x - mu) * inv_sigma; + return -0.5 * z * z - log_sigma - LOG_SQRT_2PI; + } + + inline double normal_cdf_core(const double x, const double mu, const double scale) { + return 0.5 * (1.0 + std::erf((x - mu) / scale)); + } + } // namespace + + double normal_pdf_scalar(const double x, const double mu, const double sigma) { if (!std::isfinite(x) || !std::isfinite(mu) || !std::isfinite(sigma) || sigma <= 0.0) { return std::numeric_limits::quiet_NaN(); } - const double z = (x - mu) / sigma; - return std::exp(-0.5 * z * z) / (sigma * SQRT_2PI); + return normal_pdf_core(x, mu, sigma, sigma * SQRT_2PI); } double normal_logpdf_scalar(const double x, const double mu, const double sigma) { @@ -22,9 +48,7 @@ namespace fastdist::math { return std::numeric_limits::quiet_NaN(); } - const double inv_sigma = 1.0 / sigma; - const double z = (x - mu) * inv_sigma; - return -0.5 * z * z - std::log(sigma) - LOG_SQRT_2PI; + return normal_logpdf_core(x, mu, 1.0 / sigma, std::log(sigma)); } double normal_cdf_scalar(const double x, const double mu, const double sigma) { @@ -32,8 +56,7 @@ namespace fastdist::math { return std::numeric_limits::quiet_NaN(); } - const double z = (x - mu) / (sigma * std::sqrt(2.0)); - return 0.5 * (1.0 + std::erf(z)); + return normal_cdf_core(x, mu, sigma * std::sqrt(2.0)); } double normal_mean(const double mu) { @@ -92,22 +115,67 @@ namespace fastdist::math { // Batch Functions void normal_pdf_batch(const double* x_data, double* output, const size_t n, const double mu, const double sigma, const double stepSize) { + // Parameter validity does not vary across the array, so it is checked + // once here rather than on every element. + if (!std::isfinite(mu) || !std::isfinite(sigma) || sigma <= 0.0) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + + const double denom = sigma * SQRT_2PI; + for (size_t i = 0; i < n; i++) { - output[i] = normal_pdf_scalar(x_data[i] + stepSize * static_cast(i), mu, sigma); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + output[i] = normal_pdf_core(x, mu, sigma, denom); } } void normal_logpdf_batch(const double* x_data, double* output, const size_t n, const double mu, const double sigma, const double stepSize) { + // Parameter validity does not vary across the array, so it is checked + // once here rather than on every element. + if (!std::isfinite(mu) || !std::isfinite(sigma) || sigma <= 0.0) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + + // log(sigma) in particular is a transcendental call that used to run + // once per element for a value that never changes. + const double inv_sigma = 1.0 / sigma; + const double log_sigma = std::log(sigma); + for (size_t i = 0; i < n; i++) { - output[i] = normal_logpdf_scalar(x_data[i] + stepSize * static_cast(i), mu, sigma); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + output[i] = normal_logpdf_core(x, mu, inv_sigma, log_sigma); } } void normal_cdf_batch(const double* x_data, double* output, const size_t n, const double mu, const double sigma, const double stepSize) { + // Parameter validity does not vary across the array, so it is checked + // once here rather than on every element. + if (!std::isfinite(mu) || !std::isfinite(sigma) || sigma <= 0.0) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + + const double scale = sigma * std::sqrt(2.0); + for (size_t i = 0; i < n; i++) { - output[i] = normal_cdf_scalar(x_data[i] + stepSize * static_cast(i), mu, sigma); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + output[i] = normal_cdf_core(x, mu, scale); } } diff --git a/src/math/poisson.cpp b/src/math/poisson.cpp index 57ca3ad..c953d4b 100644 --- a/src/math/poisson.cpp +++ b/src/math/poisson.cpp @@ -40,13 +40,35 @@ namespace fastdist::math { const int ki = static_cast(std::floor(x)); + // Consecutive PMF terms are related by P(i) = P(i-1) * lambda / i, so + // the sum needs one exp in total rather than a log, an lgamma and an + // exp per term. That is the difference between this being the slowest + // path in the library and it being competitive -- see BENCHMARKS.md. + // + // The recurrence has to start from P(0) = exp(-lambda), which underflows + // to zero for large lambda and would collapse the whole sum to zero even + // where the true CDF is O(1). Past that point, fall back to evaluating + // each term in log space, which stays accurate because the exponent + // i*log(lambda) - lambda - lgamma(i+1) remains small near i = lambda. + constexpr double MAX_RECURRENCE_LAMBDA = 700.0; + + if (lambda <= MAX_RECURRENCE_LAMBDA) { + double term = std::exp(-lambda); + double sum = term; + for (int i = 1; i <= ki; ++i) { + term *= lambda / static_cast(i); + sum += term; + } + // Summing PMF terms accumulates rounding error, so the total can + // land a few ULP above 1.0 + return std::min(sum, 1.0); + } + + const double log_lambda = std::log(lambda); double sum = 0.0; for (int i = 0; i <= ki; ++i) { - sum += poisson_pmf_scalar(i, lambda); + sum += std::exp(static_cast(i) * log_lambda - lambda - std::lgamma(i + 1.0)); } - - // Summing PMF terms accumulates rounding error, so the total can - // land a few ULP above 1.0 return std::min(sum, 1.0); } @@ -99,8 +121,27 @@ namespace fastdist::math { // Batch Functions void poisson_pmf_batch(const double* x_data, double* output, const size_t n, const double lambda, const int stepSize) { + // lambda is fixed across the array, so both its validation and log() are + // hoisted; log(lambda) used to be a transcendental call per element for a + // value that never changes. + if (!std::isfinite(lambda) || lambda <= 0.0) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + + const double log_lambda = std::log(lambda); + for (size_t i = 0; i < n; i++) { - output[i] = poisson_pmf_scalar(x_data[i] + stepSize * static_cast(i), lambda); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + if (x < 0.0 || std::floor(x) != x) { + output[i] = 0.0; + continue; + } + output[i] = std::exp(x * log_lambda - lambda - std::lgamma(x + 1.0)); } } void poisson_cdf_batch(const double* x_data, double* output, const size_t n, const double lambda, diff --git a/src/math/uniform.cpp b/src/math/uniform.cpp index cfee811..cc69ef7 100644 --- a/src/math/uniform.cpp +++ b/src/math/uniform.cpp @@ -1,4 +1,5 @@ // Function declarations for continuous uniform distribution functions +#include #include #include #include @@ -97,14 +98,49 @@ namespace fastdist::math { // Batch Functions void uniform_pdf_batch(const double* x_data, double* output, const size_t n, const double a, const double b, const double stepSize) { + // The density is constant across the support, so the whole value -- not + // just the validation -- is loop-invariant. This used to be a division + // per element for a number that never changes. + if (!std::isfinite(a) || !std::isfinite(b) || a >= b) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + + const double density = 1.0 / (b - a); + for (size_t i = 0; i < n; i++) { - output[i] = uniform_pdf_scalar(x_data[i] + stepSize * static_cast(i), a, b); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + output[i] = (x < a || x > b) ? 0.0 : density; } } void uniform_cdf_batch(const double* x_data, double* output, const size_t n, const double a, const double b, const double stepSize) { + if (!std::isfinite(a) || !std::isfinite(b) || a >= b) { + std::fill_n(output, n, std::numeric_limits::quiet_NaN()); + return; + } + + // Kept as a division by the hoisted range rather than a multiply by its + // reciprocal, so results stay bit-identical to the scalar path. + const double range = b - a; + for (size_t i = 0; i < n; i++) { - output[i] = uniform_cdf_scalar(x_data[i] + stepSize * static_cast(i), a, b); + const double x = x_data[i] + stepSize * static_cast(i); + if (!std::isfinite(x)) { + output[i] = std::numeric_limits::quiet_NaN(); + continue; + } + if (x <= a) { + output[i] = 0.0; + } else if (x >= b) { + output[i] = 1.0; + } else { + output[i] = (x - a) / range; + } } } diff --git a/tests/python/test_negative_binomial.py b/tests/python/test_negative_binomial.py index 8279e34..b93f67a 100644 --- a/tests/python/test_negative_binomial.py +++ b/tests/python/test_negative_binomial.py @@ -152,21 +152,16 @@ def test_cdf_approaches_one_in_the_tail(r, p): # --------------------------------------------------------------------------- # Numeric range # -# KNOWN BUG: the PMF evaluates C(k + r - 1, k) using raw factorials, so the -# intermediate (k + r - 1)! overflows a double once k + r - 1 > 170. The result -# is inf at k = 170 and nan from k = 200 onward, and cdf_scalar(200) is nan. +# REGRESSION: the PMF used to evaluate C(k + r - 1, k) from raw factorials, so +# the intermediate (k + r - 1)! overflowed a double once k + r - 1 > 170 -- +# inf at k = 170, nan from k = 200 onward, and cdf_scalar(200) nan with it. # -# The coefficient itself is small -- for r = 3, k = 200 it is C(202, 200) = -# 20301 and the true PMF is 1.58e-57, comfortably inside double range. The -# Binomial implementation handles n = 1000 without difficulty, so this is a -# defect in the negative binomial routine rather than an inherent limit. +# The coefficient itself is small: for r = 3, k = 200 it is C(202, 200) = 20301 +# and the true PMF is 1.58e-57, comfortably inside double range. Evaluating the +# whole PMF in log space via lgamma keeps the intermediates small and removes +# the ceiling entirely. These cases stay to hold that fixed. # --------------------------------------------------------------------------- -OVERFLOW_BUG = "C(k + r - 1, k) computed via raw factorials; overflows for k >= 170" - - -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason=OVERFLOW_BUG) @pytest.mark.parametrize("k", [170, 200, 500, 1000]) def test_pmf_stays_finite_in_the_far_tail(k): value = NegativeBinomial(3, 0.5).pmf_scalar(k) @@ -174,8 +169,6 @@ def test_pmf_stays_finite_in_the_far_tail(k): assert value >= 0.0 -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason=OVERFLOW_BUG) @pytest.mark.parametrize("k", [200, 500]) def test_pmf_matches_closed_form_in_the_far_tail(k): assert NegativeBinomial(3, 0.5).pmf_scalar(k) == pytest.approx( @@ -183,8 +176,6 @@ def test_pmf_matches_closed_form_in_the_far_tail(k): ) -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason=OVERFLOW_BUG) @pytest.mark.parametrize("k", [200, 500]) def test_cdf_stays_finite_in_the_far_tail(k): assert NegativeBinomial(3, 0.5).cdf_scalar(k) == pytest.approx(1.0, abs=1e-9) From 8db7afb78625a7c220acb2b49a9e45049ecf3165 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:25:44 -0400 Subject: [PATCH 05/15] Record the optimisation pass in the benchmark log Full run at 73ba7d6905 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 --- BENCHMARKS.md | 108 +++++ ...0.1.0_20260905T062420+0000_73ba7d6905.json | 450 ++++++++++++++++++ 2 files changed, 558 insertions(+) create mode 100644 benchmarks/results/0.1.0_20260905T062420+0000_73ba7d6905.json diff --git a/BENCHMARKS.md b/BENCHMARKS.md index b3f8992..4448616 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -162,6 +162,114 @@ evaluating distribution functions, not for bulk variate generation. --- +## Unreleased — discrete CDF recurrences and batch invariant hoisting + +Commit `73ba7d6905`, measured against the v0.1.0 baseline above on the same +machine in the same session. 29 cases improved, none regressed. + +Two changes: the three discrete CDFs that summed PMF terms now use recurrences, +and the batch paths hoist parameter validation and loop-invariant terms out of +their loops. See the commit for the derivations and the underflow fallbacks. + +The largest movements, as reported by `compare.py --latest`. `change` is +fastdist's own time, so negative is faster: + +| case | change | +|---|---:| +| `batch/poisson_cdf n=1,000` | -97.7% | +| `batch/poisson_cdf n=100,000` | -97.0% | +| `batch/poisson_cdf n=1,000,000` | -96.9% | +| `batch/normal_logpdf n=100,000` | -80.2% | +| `batch/normal_logpdf n=1,000,000` | -68.3% | +| `batch/normal_logpdf n=1,000` | -62.0% | +| `batch/uniform_pdf n=100,000` | -28.7% | +| `batch/uniform_pdf n=1,000,000` | -26.4% | +| `batch/normal_cdf n=100,000` | -25.7% | +| `batch/uniform_cdf n=1,000,000` | -25.2% | +| `batch/normal_cdf n=1,000,000` | -23.6% | +| `batch/normal_cdf n=1,000` | -22.5% | + + +- **Version** 0.1.0 (`73ba7d6905` on `fix/flaky-rng-tolerances`, working tree dirty) +- **Measured** 2026-09-05T06:24:20+00:00 +- **CPU** AMD Ryzen 7 7700 8-Core Processor +- **Platform** Windows-11-10.0.26200-SP0 +- **Toolchain** Python 3.14.2, numpy 2.5.2, scipy 1.18.1 +- **CUDA** not built + +### batch (vs vectorised SciPy) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_pdf` | 1,000 | 5.24 us | 31.57 us (scipy) | **6.02x** | 1.1e-16 | +| `normal_cdf` | 1,000 | 6.01 us | 29.63 us (scipy) | **4.93x** | 2.2e-16 | +| `normal_logpdf` | 1,000 | 1.91 us | 32.20 us (scipy) | **16.89x** | 8.9e-16 | +| `exponential_pdf` | 1,000 | 4.14 us | 29.31 us (scipy) | **7.08x** | 0.0e+00 | +| `exponential_cdf` | 1,000 | 4.16 us | 29.63 us (scipy) | **7.12x** | 8.3e-17 | +| `uniform_pdf` | 1,000 | 2.02 us | 32.99 us (scipy) | **16.30x** | 0.0e+00 | +| `uniform_cdf` | 1,000 | 2.13 us | 31.45 us (scipy) | **14.78x** | 0.0e+00 | +| `poisson_pmf` | 1,000 | 40.52 us | 37.70 us (scipy) | **0.93x** | 2.0e-19 | +| `poisson_cdf` | 1,000 | 9.95 us | 70.07 us (scipy) | **7.05x** | 2.2e-16 | +| `bernoulli_pmf` | 1,000 | 1.93 us | 48.90 us (scipy) | **25.36x** | 2.2e-16 | +| `normal_pdf` | 100,000 | 413.90 us | 1.44 ms (scipy) | **3.48x** | 1.1e-16 | +| `normal_cdf` | 100,000 | 528.60 us | 2.06 ms (scipy) | **3.89x** | 2.2e-16 | +| `normal_logpdf` | 100,000 | 77.90 us | 1.67 ms (scipy) | **21.45x** | 8.9e-16 | +| `exponential_pdf` | 100,000 | 312.50 us | 1.52 ms (scipy) | **4.87x** | 0.0e+00 | +| `exponential_cdf` | 100,000 | 312.10 us | 1.64 ms (scipy) | **5.27x** | 1.1e-16 | +| `uniform_pdf` | 100,000 | 92.00 us | 1.55 ms (scipy) | **16.79x** | 0.0e+00 | +| `uniform_cdf` | 100,000 | 102.00 us | 1.65 ms (scipy) | **16.18x** | 0.0e+00 | +| `poisson_pmf` | 100,000 | 4.01 ms | 3.16 ms (scipy) | **0.79x** | 2.0e-19 | +| `poisson_cdf` | 100,000 | 1.31 ms | 6.05 ms (scipy) | **4.63x** | 2.2e-16 | +| `bernoulli_pmf` | 100,000 | 290.80 us | 3.48 ms (scipy) | **11.96x** | 2.2e-16 | +| `normal_pdf` | 1,000,000 | 4.73 ms | 19.66 ms (scipy) | **4.16x** | 1.1e-16 | +| `normal_cdf` | 1,000,000 | 5.86 ms | 22.84 ms (scipy) | **3.90x** | 2.2e-16 | +| `normal_logpdf` | 1,000,000 | 1.42 ms | 21.70 ms (scipy) | **15.33x** | 8.9e-16 | +| `exponential_pdf` | 1,000,000 | 3.92 ms | 18.36 ms (scipy) | **4.68x** | 0.0e+00 | +| `exponential_cdf` | 1,000,000 | 3.74 ms | 18.55 ms (scipy) | **4.96x** | 1.7e-16 | +| `uniform_pdf` | 1,000,000 | 1.38 ms | 17.70 ms (scipy) | **12.80x** | 0.0e+00 | +| `uniform_cdf` | 1,000,000 | 1.51 ms | 19.54 ms (scipy) | **12.97x** | 0.0e+00 | +| `poisson_pmf` | 1,000,000 | 42.65 ms | 37.75 ms (scipy) | **0.89x** | 2.0e-19 | +| `poisson_cdf` | 1,000,000 | 14.09 ms | 63.76 ms (scipy) | **4.53x** | 2.2e-16 | +| `bernoulli_pmf` | 1,000,000 | 3.52 ms | 38.78 ms (scipy) | **11.01x** | 2.2e-16 | + +### scalar (per-call cost, not throughput) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_pdf` | 20,000 | 7.12 ms | 451.64 ms (scipy) | **63.40x** | 1.1e-16 | +| `normal_cdf` | 20,000 | 7.35 ms | 439.65 ms (scipy) | **59.79x** | 2.2e-16 | + +### sample (vs numpy) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_sample` | 100,000 | 27.25 ms | 867.30 us (numpy) | **0.03x** | - | +| `uniform_sample` | 100,000 | 23.73 ms | 227.10 us (numpy) | **0.01x** | - | +| `normal_sample` | 1,000,000 | 288.47 ms | 10.02 ms (numpy) | **0.03x** | - | +| `uniform_sample` | 1,000,000 | 253.13 ms | 3.19 ms (numpy) | **0.01x** | - | + +### Reading + +`poisson_cdf` was the headline defect in the baseline and is now the largest +win: 43.47ms to 1.32ms at 100k elements, moving from 6.6x slower than SciPy to +4.9x faster. Agreement with SciPy tightened from 3.3e-16 to 2.2e-16 at the same +time, which is the expected consequence of doing far fewer floating-point +operations to reach the same answer. + +`normal_logpdf` improved 80% purely from hoisting `log(sigma)`, which the batch +path had been recomputing per element for a value fixed across the whole array. +It is now the fastest continuous case in the suite at 22x SciPy. + +The remaining known gaps are unchanged and still worth recording: + +- `poisson_pmf` is 0.81x. The per-element `lgamma` dominates and is not + loop-invariant, so hoisting cannot reach it. Beating SciPy here needs a + different evaluation strategy, not tuning. +- Sampling is still 30-100x slower than numpy. Unchanged, and structural: it + needs a batch sampling entry point that does not exist yet. + +--- + ## Changes to record here Add an entry when a release ships, or when a change is made specifically to diff --git a/benchmarks/results/0.1.0_20260905T062420+0000_73ba7d6905.json b/benchmarks/results/0.1.0_20260905T062420+0000_73ba7d6905.json new file mode 100644 index 0000000..1f22308 --- /dev/null +++ b/benchmarks/results/0.1.0_20260905T062420+0000_73ba7d6905.json @@ -0,0 +1,450 @@ +{ + "environment": { + "timestamp_utc": "2026-09-05T06:24:20+00:00", + "fastdist_version": "0.1.0", + "git_commit": "73ba7d69055948cba21190122d6c409a29c4a903", + "git_branch": "fix/flaky-rng-tolerances", + "git_dirty": true, + "cuda_available": false, + "python": "3.14.2", + "numpy": "2.5.2", + "scipy": "1.18.1", + "platform": "Windows-11-10.0.26200-SP0", + "processor": "AMD Ryzen 7 7700 8-Core Processor", + "machine": "AMD64" + }, + "results": [ + { + "group": "batch", + "case": "normal_pdf", + "n": 1000, + "fastdist_s": 5.243999767117202e-06, + "baseline_s": 3.157400002237409e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 27.8794839408897, + "baseline_noise_pct": 7.607524107998575, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 6.020976625582756 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000, + "fastdist_s": 6.011999794282019e-06, + "baseline_s": 2.9628000338561832e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.7651422461864626, + "baseline_noise_pct": 3.8612102898679925, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.928143937519902 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000, + "fastdist_s": 1.905999961309135e-06, + "baseline_s": 3.219599951989949e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.6295955727972818, + "baseline_noise_pct": 4.683813037916011, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 16.891920342844962 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000, + "fastdist_s": 4.138000076636672e-06, + "baseline_s": 2.930999966338277e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.531647777544893, + "baseline_noise_pct": 3.596043912606601, + "max_abs_diff": 0.0, + "speedup": 7.083131735271901 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000, + "fastdist_s": 4.161999677307904e-06, + "baseline_s": 2.9628000338561832e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.6247044688298693, + "baseline_noise_pct": 1.8225983834309571, + "max_abs_diff": 8.326672684688674e-17, + "speedup": 7.118693569367608 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000, + "fastdist_s": 2.0239996956661345e-06, + "baseline_s": 3.299200034234673e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.49410374136698193, + "baseline_noise_pct": 0.7092629520292216, + "max_abs_diff": 0.0, + "speedup": 16.300397877030544 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000, + "fastdist_s": 2.127999905496836e-06, + "baseline_s": 3.144800022710115e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.6578997755941307, + "baseline_noise_pct": 3.6313913646507094, + "max_abs_diff": 0.0, + "speedup": 14.778196251732826 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000, + "fastdist_s": 4.052400006912649e-05, + "baseline_s": 3.7700000102631745e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.24676778995692783, + "baseline_noise_pct": 0.5888592667993294, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.9303129019426137 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000, + "fastdist_s": 9.946000063791872e-06, + "baseline_s": 7.007200038060546e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.764126243087026, + "baseline_noise_pct": 17.904440319454316, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 7.0452443124045985 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000, + "fastdist_s": 1.92799954675138e-06, + "baseline_s": 4.889799980446696e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.8299122235512824, + "baseline_noise_pct": 4.5400630554029675, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 25.362039055900503 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 100000, + "fastdist_s": 0.00041390000842511654, + "baseline_s": 0.0014389999851118773, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.7489718381244842, + "baseline_noise_pct": 6.587909961282603, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 3.4766850829195466 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 100000, + "fastdist_s": 0.0005286000086925924, + "baseline_s": 0.002058300015050918, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.502460501433774, + "baseline_noise_pct": 3.25997116906677, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.8938705660293005 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 100000, + "fastdist_s": 7.790001109242439e-05, + "baseline_s": 0.0016709000046830624, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.2567039673857815, + "baseline_noise_pct": 6.679033835044771, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 21.449290972508653 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 100000, + "fastdist_s": 0.0003124999930150807, + "baseline_s": 0.0015232000150717795, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.48000366568573566, + "baseline_noise_pct": 1.1948524740376525, + "max_abs_diff": 0.0, + "speedup": 4.874240157177452 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 100000, + "fastdist_s": 0.00031209998996928334, + "baseline_s": 0.0016442999767605215, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.2816409749121431, + "baseline_noise_pct": 0.9122439407595886, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 5.268503779581513 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 100000, + "fastdist_s": 9.200000204145908e-05, + "baseline_s": 0.0015450000064447522, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.10869647742428575, + "baseline_noise_pct": 6.2912610371331725, + "max_abs_diff": 0.0, + "speedup": 16.793477958277762 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 100000, + "fastdist_s": 0.00010199999087490141, + "baseline_s": 0.0016506000247318298, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.19607994195211964, + "baseline_noise_pct": 1.7145252109326794, + "max_abs_diff": 0.0, + "speedup": 16.182354631347167 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 100000, + "fastdist_s": 0.0040145000093616545, + "baseline_s": 0.003160399995977059, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.02740015617255305, + "baseline_noise_pct": 1.629541003106192, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.7872462295695932 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 100000, + "fastdist_s": 0.0013082000077702105, + "baseline_s": 0.0060505000001285225, + "baseline_name": "scipy", + "fastdist_noise_pct": 3.6920951798922386, + "baseline_noise_pct": 0.8825716743936082, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.625057303310544 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 100000, + "fastdist_s": 0.0002908000024035573, + "baseline_s": 0.0034771999926306307, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.0316347024944024, + "baseline_noise_pct": 1.1618548450721544, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 11.957358885455411 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 1000000, + "fastdist_s": 0.004731399996671826, + "baseline_s": 0.01965940001537092, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.5875641031334836, + "baseline_noise_pct": 1.1078668232204452, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 4.155091522424608 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000000, + "fastdist_s": 0.00585980000323616, + "baseline_s": 0.02283909998368472, + "baseline_name": "scipy", + "fastdist_noise_pct": 7.048022139731513, + "baseline_noise_pct": 1.9790622564895242, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.8975903565090095 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000000, + "fastdist_s": 0.0014157999830786139, + "baseline_s": 0.021702900005038828, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.230824994110061, + "baseline_noise_pct": 3.3313520519695023, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 15.329072089580432 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000000, + "fastdist_s": 0.003922700008843094, + "baseline_s": 0.01836290000937879, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.514058878042331, + "baseline_noise_pct": 5.364621011637729, + "max_abs_diff": 0.0, + "speedup": 4.681188968818058 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000000, + "fastdist_s": 0.0037422000023070723, + "baseline_s": 0.018549100001109764, + "baseline_name": "scipy", + "fastdist_noise_pct": 3.2547698158203184, + "baseline_noise_pct": 2.051312392023963, + "max_abs_diff": 1.6653345369377348e-16, + "speedup": 4.956736676199619 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000000, + "fastdist_s": 0.0013827999937348068, + "baseline_s": 0.01770200001192279, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.1405837932073917, + "baseline_noise_pct": 4.28765104071291, + "max_abs_diff": 0.0, + "speedup": 12.801562114642069 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000000, + "fastdist_s": 0.0015073000104166567, + "baseline_s": 0.019544200011296198, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.4082806488935775, + "baseline_noise_pct": 2.691335478509691, + "max_abs_diff": 0.0, + "speedup": 12.96636361456249 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000000, + "fastdist_s": 0.04264639999018982, + "baseline_s": 0.0377493999840226, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.0702802582503077, + "baseline_noise_pct": 1.6071779195555445, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.885172019038097 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000000, + "fastdist_s": 0.014089000003878027, + "baseline_s": 0.06376290001207963, + "baseline_name": "scipy", + "fastdist_noise_pct": 6.059337143188412, + "baseline_noise_pct": 1.2599802964161702, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.525722194231582 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000000, + "fastdist_s": 0.0035237999982200563, + "baseline_s": 0.03878130001248792, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.5693281273613617, + "baseline_noise_pct": 1.84057778511608, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 11.005533807842996 + }, + { + "group": "scalar", + "case": "normal_pdf", + "n": 20000, + "fastdist_s": 0.007124100025976077, + "baseline_s": 0.4516408999916166, + "baseline_name": "scipy", + "fastdist_noise_pct": 3.4221861536973686, + "baseline_noise_pct": 5.455130392311249, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 63.396204200507 + }, + { + "group": "scalar", + "case": "normal_cdf", + "n": 20000, + "fastdist_s": 0.007353400025749579, + "baseline_s": 0.43965400001616217, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.219381383429801, + "baseline_noise_pct": 1.218799321600988, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 59.78921294593183 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 100000, + "fastdist_s": 0.02724950000992976, + "baseline_s": 0.0008673000265844166, + "baseline_name": "numpy", + "fastdist_noise_pct": 2.121506735704241, + "baseline_noise_pct": 1.2567700044449364, + "max_abs_diff": null, + "speedup": 0.03182810790173659 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 100000, + "fastdist_s": 0.02373069999157451, + "baseline_s": 0.00022709998302161694, + "baseline_name": "numpy", + "fastdist_noise_pct": 17.008347842272727, + "baseline_noise_pct": 0.1761480314992921, + "max_abs_diff": null, + "speedup": 0.009569881339456816 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 1000000, + "fastdist_s": 0.28846969999722205, + "baseline_s": 0.01002149999840185, + "baseline_name": "numpy", + "fastdist_noise_pct": 1.3494658138910405, + "baseline_noise_pct": 0.9210198507262023, + "max_abs_diff": null, + "speedup": 0.034740217078252436 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 1000000, + "fastdist_s": 0.2531314000079874, + "baseline_s": 0.0031857999856583774, + "baseline_name": "numpy", + "fastdist_noise_pct": 3.073542036221331, + "baseline_noise_pct": 4.57028163990082, + "max_abs_diff": null, + "speedup": 0.012585558273520597 + } + ] +} From 3c568c5366d2981621a7e541498a68ffa2116706 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:28:19 -0400 Subject: [PATCH 06/15] Document the build, seeding, benchmarks and tolerance policy 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 --- README.md | 201 ++++++++++++++++-- python/fastdist/distributions/bernoulli.py | 12 +- python/fastdist/distributions/beta.py | 9 +- python/fastdist/distributions/binomial.py | 9 +- python/fastdist/distributions/chi_square.py | 9 +- .../distributions/discrete_uniform.py | 9 +- python/fastdist/distributions/exponential.py | 12 +- python/fastdist/distributions/gamma.py | 9 +- python/fastdist/distributions/geometric.py | 9 +- .../distributions/negative_binomial.py | 9 +- python/fastdist/distributions/normal.py | 12 +- python/fastdist/distributions/poisson.py | 12 +- python/fastdist/distributions/uniform.py | 12 +- python/fastdist/distributions/utils.py | 12 +- 14 files changed, 283 insertions(+), 53 deletions(-) diff --git a/README.md b/README.md index 7ca48c2..810786b 100644 --- a/README.md +++ b/README.md @@ -21,28 +21,59 @@ git clone https://github.com/ghosteau/fastdist.git ## Building and Installing `fastdist` -### Building the C++ Extension and Python Wheel +### Installing -1. Build the C++ project using CMake. - - This produces the compiled extension (`.pyd` on Windows) in your CMake build directory (e.g., - `cmake-build-debug`). +From the **project root**: -2. From the **project root**, build the Python wheel: +```bash +pip install . +``` + +That is the whole thing. `setup.py` drives CMake, and `pyproject.toml` declares +the build dependencies (including `pybind11` and `cmake`), so pip provisions +them in an isolated build environment. + +To build a wheel without installing it: + +```bash +pip install build +python -m build --wheel +``` + +The wheel lands in `dist/`. (`python setup.py bdist_wheel` still works but is +deprecated upstream; prefer `python -m build`.) + +### Building the C++ project directly + +Needed when working on the C++ side, running the C++ tests, or using an IDE's +CMake integration: ```bash -python3 setup.py bdist_wheel +cmake -S . -B build -DCMAKE_BUILD_TYPE=Release +cmake --build build --target fastdist_tests --parallel +ctest --test-dir build --output-on-failure ``` -**Important:** +On Windows the Visual Studio generator is multi-config, so `CMAKE_BUILD_TYPE` +is ignored -- pick the configuration at build and test time instead: -- This command must be run from the project root. +```powershell +cmake -S . -B build +cmake --build build --target fastdist_tests --config Release --parallel +ctest --test-dir build -C Release --output-on-failure +``` -3. Install the generated wheel: +**If CMake reports it cannot find pybind11**, it is resolving pybind11 from the +active interpreter and that interpreter has none installed. Either install it +(`pip install pybind11`) or point CMake at one explicitly: ```bash -pip install .\dist\fastdist--cpXXX-cpXXX-win_amd64.whl --force-reinstall +cmake -S . -B build -Dpybind11_DIR="$(python -m pybind11 --cmakedir)" ``` +Pass `-DPython_EXECUTABLE=...` as well if the interpreter you want is not the +first one on `PATH`. + --- ## Creating A Distribution Class @@ -124,6 +155,100 @@ When cleanup is enabled, only the final wheel files will remain. --- +## Reproducible Sampling + +Every `*_sample()` function draws from one shared Mersenne Twister engine. +Seeding it makes a run reproducible: + +```python +import fastdist + +fastdist.seed(12345) +a = [fastdist.Normal(0, 1).sample() for _ in range(5)] + +fastdist.seed(12345) +b = [fastdist.Normal(0, 1).sample() for _ in range(5)] + +assert a == b # exactly equal, not merely close + +fastdist.seed_from_entropy() # back to non-deterministic +``` + +From C++, the same thing lives in `fastdist/math/rng.h`: + +```cpp +#include + +fastdist::math::seed_rng(12345); +fastdist::math::seed_rng_from_entropy(); +``` + +Two limits are worth knowing before relying on this: + +1. **Seeding is per-thread.** The engine is `thread_local`, so a worker thread + that has not been seeded keeps its own entropy-initialised stream. This is + what makes concurrent sampling lock-free; it also means one `seed()` call + does not cover threads you spawn. + +2. **Reproducible per platform, not across them.** `std::mt19937` is specified + bit-for-bit by the C++ standard, but the distribution adaptors built on it + (`std::normal_distribution` and friends) are not. The same seed therefore + produces different samples under libstdc++, libc++ and MSVC. A seed pins a + run on one platform and toolchain, not across all of them. + +--- + +## Benchmarks + +Performance is measured against SciPy and recorded in +[BENCHMARKS.md](BENCHMARKS.md), with the raw JSON for every run kept under +`benchmarks/results/`. + +```bash +pip install scipy # baseline only; not needed to build or use fastdist +python benchmarks/run.py # full suite, writes a report +python benchmarks/run.py --quick # one array size, for a fast check +python benchmarks/compare.py --latest # diff the two newest reports +``` + +`compare.py` exits non-zero if any case regressed by more than 5%, so it can +gate a change. Every case checks that fastdist and the baseline agree +numerically before either is timed -- a speedup on a wrong answer is not a +speedup. + +Add an entry to BENCHMARKS.md when a release ships, or when a change is made +specifically to move performance. Generate the table with +`python benchmarks/table.py --latest` rather than transcribing numbers. + +--- + +## Testing + +```bash +ctest --test-dir build --output-on-failure # C++ +pytest # Python +``` + +### Tolerances in sampling tests + +The RNG test blocks draw a large sample and compare its mean and variance +against theory. Two rules keep those honest: + +1. **Seed first.** Every RNG block calls `seed_rng()` before sampling, so it is + deterministic: it either always passes or always fails on a given toolchain, + never intermittently. + +2. **Size the tolerance from the estimator's standard error**, not from a round + number. Each block's tolerance is roughly 5x the standard error of the + statistic being checked, and the SE is recorded in a comment at the site. + +Round-number tolerances are how this suite acquired a 7.7% flake rate: several +sat near 2 sigma of their estimator's own noise and failed at about the rate a +2 sigma bound fails (ghosteau/fastdist#2). If you change `N` or a distribution +parameter, recompute the standard error and resize the tolerance with it. + +--- + ## Code Formatting and Pre-Commit Hooks This repository enforces consistent formatting using `clang-format`. @@ -220,26 +345,60 @@ Python Bindings: Long-term plans: +Performance (see [BENCHMARKS.md](BENCHMARKS.md) for what is measured today): + +- Add batch sampling entry points (`normal_sample_batch(n)` returning an array). + Sampling is currently 30-100x slower than numpy because every variate crosses + the Python/C++ boundary individually. This is the single largest gap in the + library. +- Speed up `poisson_pmf`, still ~0.8x SciPy. The per-element `lgamma` dominates + and is not loop-invariant, so it needs a different evaluation strategy. +- Extend batch invariant hoisting to the distributions the benchmark suite does + not yet cover (binomial, geometric, beta, gamma, chi-square, discrete uniform, + negative binomial). The pattern is established in `normal.cpp`. +- Precalculate reused values on CPU and send to GPU +- Add a memory-constraint option to CUDA, so a limited GPU can cap its streaming + budget + +Correctness and API: + - Add Hypergeometric Distribution - Make auto_tune() dynamically find the sign flip +- Check for all isfinite values (currently only set up in normal) +- Add specific parameters in all return _core._(x, a, b) → (x=x, a=a, b=b) +- Expose the `step_size` defaults the C++ headers declare through the pybind11 + bindings; callers currently have to pass it explicitly +- Make `step_size` consistently typed -- it is `double` on the continuous batch + functions and `int` on the discrete ones +- Merge validation checks and CUDA availability into a singular function for cleanliness +- Look into the usage of @classmethod and check for redundancies in the Python classes +- Refine Utils class to be more efficient and comprehensive +- Seed the CUDA RNG alongside the CPU engine, so `seed()` covers both backends + +CUDA: + - Add CUDA/Batch extern functions - Set up CI for cuda tests - Add cuda implementation for all classes -- Fix up the python-distro.yml file to be more efficient and comprehensive -- Look into the usage of @classmethod and check for redundancies in the Python classes -- Add specific parameters in all return _core._(x, a, b) → (x=x, a=a, b=b) -- Check for all isfinite values (currently only set up in normal) -- Merge validation checks and CUDA availability into a singular function for cleanliness -- Use size_t instead of int in all cuda files -- Add memory constraint option to cuda where if you have limited gpu memory you can set what your limit for streaming is -- Update all docstrings to match each other and be comprehensive - Create new cuda tests -- Refine Utils class to be more efficient and comprehensive +- Use size_t instead of int in all cuda files - Add batch and cuda functions to the C API +- Benchmark the CUDA backend and record it in BENCHMARKS.md + +Tooling and docs: + +- Fix up the python-distro.yml file to be more efficient and comprehensive +- Update all docstrings to match each other and be comprehensive - Update pynvml to nvidia-ml-py -- Precalculate reused values on CPU and send to GPU - Make a full, comprehensive documentation page -- In the future, try to get the library on pip +- In the future, try to get the library on pip + +Done since v0.1.0: + +- ~~Seedable RNG for reproducible sampling~~ (`fastdist.seed()`) +- ~~Flaky RNG tests~~ (tolerances sized from estimator standard error; #2) +- ~~Performance measurement and evidence log~~ (`benchmarks/`, BENCHMARKS.md) + --- ## Contributors diff --git a/python/fastdist/distributions/bernoulli.py b/python/fastdist/distributions/bernoulli.py index 79dfc51..e121811 100644 --- a/python/fastdist/distributions/bernoulli.py +++ b/python/fastdist/distributions/bernoulli.py @@ -1,9 +1,15 @@ # python/distributions/bernoulli.py try: from fastdist import _fastdist as _core - from fastdist import config -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc + +from fastdist import config import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/beta.py b/python/fastdist/distributions/beta.py index 08b382c..54e9bfa 100644 --- a/python/fastdist/distributions/beta.py +++ b/python/fastdist/distributions/beta.py @@ -1,8 +1,13 @@ # python/distributions/bernoulli.py try: from .. import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/binomial.py b/python/fastdist/distributions/binomial.py index df55b70..f290e94 100644 --- a/python/fastdist/distributions/binomial.py +++ b/python/fastdist/distributions/binomial.py @@ -1,8 +1,13 @@ # python/distributions/binomial.py try: from fastdist import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/chi_square.py b/python/fastdist/distributions/chi_square.py index 9ccfc88..75f4aec 100644 --- a/python/fastdist/distributions/chi_square.py +++ b/python/fastdist/distributions/chi_square.py @@ -1,8 +1,13 @@ # python/distributions/chi_square.py try: from fastdist import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/discrete_uniform.py b/python/fastdist/distributions/discrete_uniform.py index 7d29a43..64a3164 100644 --- a/python/fastdist/distributions/discrete_uniform.py +++ b/python/fastdist/distributions/discrete_uniform.py @@ -1,8 +1,13 @@ # python/distributions/discrete_uniform.py try: from fastdist import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/exponential.py b/python/fastdist/distributions/exponential.py index 025d8c4..6e6ebe0 100644 --- a/python/fastdist/distributions/exponential.py +++ b/python/fastdist/distributions/exponential.py @@ -1,9 +1,15 @@ # python/distributions/exponential.py try: from fastdist import _fastdist as _core - from fastdist import config -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc + +from fastdist import config import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/gamma.py b/python/fastdist/distributions/gamma.py index c6cb086..3656778 100644 --- a/python/fastdist/distributions/gamma.py +++ b/python/fastdist/distributions/gamma.py @@ -2,8 +2,13 @@ try: from .. import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/geometric.py b/python/fastdist/distributions/geometric.py index 843b8c4..8e245c9 100644 --- a/python/fastdist/distributions/geometric.py +++ b/python/fastdist/distributions/geometric.py @@ -1,8 +1,13 @@ # python/distributions/geometric.py try: from fastdist import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/negative_binomial.py b/python/fastdist/distributions/negative_binomial.py index 1aed9b2..abeae9e 100644 --- a/python/fastdist/distributions/negative_binomial.py +++ b/python/fastdist/distributions/negative_binomial.py @@ -1,8 +1,13 @@ # python/distributions/poisson.py try: from fastdist import _fastdist as _core -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc import numpy as np from numbers import Real diff --git a/python/fastdist/distributions/normal.py b/python/fastdist/distributions/normal.py index 0727a52..11f32ba 100644 --- a/python/fastdist/distributions/normal.py +++ b/python/fastdist/distributions/normal.py @@ -1,9 +1,15 @@ # python/distributions/normal.py try: from .. import _fastdist as _core - from fastdist import config -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc + +from .. import config import math from numbers import Real diff --git a/python/fastdist/distributions/poisson.py b/python/fastdist/distributions/poisson.py index 8249f27..9bd48b4 100644 --- a/python/fastdist/distributions/poisson.py +++ b/python/fastdist/distributions/poisson.py @@ -1,9 +1,15 @@ # python/distributions/poisson.py try: from .. import _fastdist as _core - from fastdist import config -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc + +from .. import config from numbers import Real from typing import Sequence, Union diff --git a/python/fastdist/distributions/uniform.py b/python/fastdist/distributions/uniform.py index 86a9050..281fd66 100644 --- a/python/fastdist/distributions/uniform.py +++ b/python/fastdist/distributions/uniform.py @@ -1,9 +1,15 @@ # python/distributions/uniform.py try: from .. import _fastdist as _core - from fastdist import config -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc + +from .. import config from numbers import Real from typing import Sequence, Union diff --git a/python/fastdist/distributions/utils.py b/python/fastdist/distributions/utils.py index c1d4f49..c50812a 100644 --- a/python/fastdist/distributions/utils.py +++ b/python/fastdist/distributions/utils.py @@ -1,9 +1,15 @@ # python/distributions/utils.py try: from fastdist import _fastdist as _core - from fastdist import config -except ImportError: - raise ImportError("Internal Error: C++ core (_fastdist) not found. Check package structure.") +except ImportError as exc: # pragma: no cover - only hit in a broken install + raise ImportError( + "fastdist's compiled extension (_fastdist) could not be imported. " + "Build it with `pip install .` from the repository root; importing " + "the package straight from a source checkout will not work until the " + "extension has been built." + ) from exc + +from fastdist import config import numpy as np from numbers import Real From 0e27b85b306844929db7c96558488a11f17ee710 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:48:09 -0400 Subject: [PATCH 07/15] Fix the regularized incomplete beta, which was wrong at every point 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 --- src/math/beta.cpp | 85 +++++++++++++++++++++++++++++++----- tests/python/test_beta.py | 61 +++++++++++++------------- tests/python/test_version.py | 27 ++++++++++-- 3 files changed, 128 insertions(+), 45 deletions(-) diff --git a/src/math/beta.cpp b/src/math/beta.cpp index 3f50989..330c22e 100644 --- a/src/math/beta.cpp +++ b/src/math/beta.cpp @@ -29,10 +29,26 @@ namespace fastdist::math { } // Forward declarations for internal functions - static double beta_inc_series(double a, double b, double x); + static double beta_continued_fraction(double a, double b, double x); // ------------------------- - // CDF using series + // CDF + // + // The regularized incomplete beta I_x(a,b), evaluated as + // + // I_x(a,b) = x^a (1-x)^b / (a B(a,b)) * CF(a,b,x) + // + // where CF is the standard continued fraction for the function. The + // fraction converges rapidly for x below the transition point + // (a+1)/(a+b+2) and slowly above it, so the reflection + // + // I_x(a,b) = 1 - I_{1-x}(b,a) + // + // maps every input onto the fast side before evaluating. + // + // The leading factor is formed in log space. Computing B(a,b) directly + // overflows for large a or b, and x^a underflows for large a, even when + // their combination is an ordinary number. // ------------------------- double beta_cdf_scalar(const double x, const double alpha, const double beta) { if (!std::isfinite(x) || !std::isfinite(alpha) || !std::isfinite(beta) || alpha <= 0.0 || beta <= 0.0) { @@ -42,11 +58,21 @@ namespace fastdist::math { if (x <= 0.0) return 0.0; if (x >= 1.0) return 1.0; + const double log_prefactor = std::lgamma(alpha + beta) - std::lgamma(alpha) - std::lgamma(beta) + + alpha * std::log(x) + beta * std::log1p(-x); + const double prefactor = std::exp(log_prefactor); + + double result; if (x < (alpha + 1.0) / (alpha + beta + 2.0)) { - return beta_inc_series(alpha, beta, x); + result = prefactor * beta_continued_fraction(alpha, beta, x) / alpha; } else { - return 1.0 - beta_inc_series(beta, alpha, 1.0 - x); + result = 1.0 - prefactor * beta_continued_fraction(beta, alpha, 1.0 - x) / beta; } + + // The reflection above subtracts two nearby quantities in the upper + // tail, so the result can land a few ULP outside [0, 1]. A CDF that + // reports 1 + 1e-16 breaks callers that treat it as a probability. + return std::min(std::max(result, 0.0), 1.0); } // ------------------------- @@ -81,15 +107,50 @@ namespace fastdist::math { // ------------------------- // Internal: incomplete beta series // ------------------------- - static double beta_inc_series(const double a, const double b, const double x) { - double sum = 1.0 / a; - double term = sum; - for (unsigned int n = 1; n <= MAX_ITER; ++n) { - term *= x * (a + n - 1) / (a + b + n - 1); - sum += term; - if (std::fabs(term) < EPS * std::fabs(sum)) break; + // Modified Lentz evaluation of the continued fraction for the incomplete + // beta function (Numerical Recipes 6.4). Each iteration applies two + // coefficients, the even and odd terms of the fraction. + // + // FPMIN guards the standard Lentz failure mode: a denominator that lands + // exactly on zero would otherwise propagate an infinity through the whole + // recurrence. + static double beta_continued_fraction(const double a, const double b, const double x) { + const double qab = a + b; + const double qap = a + 1.0; + const double qam = a - 1.0; + + double c = 1.0; + double d = 1.0 - qab * x / qap; + if (std::fabs(d) < FPMIN) d = FPMIN; + d = 1.0 / d; + double h = d; + + for (unsigned int m = 1; m <= MAX_ITER; ++m) { + const double m2 = 2.0 * m; + + // Even step. + double numerator = m * (b - m) * x / ((qam + m2) * (a + m2)); + d = 1.0 + numerator * d; + if (std::fabs(d) < FPMIN) d = FPMIN; + c = 1.0 + numerator / c; + if (std::fabs(c) < FPMIN) c = FPMIN; + d = 1.0 / d; + h *= d * c; + + // Odd step. + numerator = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); + d = 1.0 + numerator * d; + if (std::fabs(d) < FPMIN) d = FPMIN; + c = 1.0 + numerator / c; + if (std::fabs(c) < FPMIN) c = FPMIN; + d = 1.0 / d; + const double delta = d * c; + h *= delta; + + if (std::fabs(delta - 1.0) < EPS) break; } - return sum * std::pow(x, a) * std::pow(1.0 - x, b) / std::tgamma(a + 1.0); + + return h; } } // namespace fastdist::math diff --git a/tests/python/test_beta.py b/tests/python/test_beta.py index f5e7e38..52c1f36 100644 --- a/tests/python/test_beta.py +++ b/tests/python/test_beta.py @@ -4,7 +4,7 @@ import pytest -from conftest import EXACT, ITERATIVE +from conftest import EXACT, ITERATIVE, regularized_incomplete_beta from fastdist.distributions.beta import Beta @@ -141,48 +141,40 @@ def test_mean_lies_inside_the_support(alpha, beta): # --------------------------------------------------------------------------- # CDF # -# KNOWN BUG: the regularized incomplete beta in src/math/beta.cpp is incorrect. -# It disagrees with numerical integration at every tested point, returns 0.6534 -# for Beta(1, 1) at x=0.5 where the exact answer is 0.5, and returns a negative -# value (-0.5958) for Beta(0.5, 0.5) at x=0.5, which is impossible for a CDF. -# These tests assert the correct behaviour and are marked strict-xfail so they -# turn into failures the moment the backend is fixed and the marker goes stale. +# REGRESSION: the regularized incomplete beta in src/math/beta.cpp used to be +# wrong at every point. It summed a hypergeometric series with the coefficient +# ratio inverted and normalised by Gamma(a+1) instead of B(a,b), which returned +# 0.6534 for Beta(1, 1) at x=0.5 where the exact answer is 0.5, and -0.5958 for +# Beta(0.5, 0.5) at x=0.5 -- impossible for a CDF. +# +# It is now the standard modified-Lentz continued fraction with the reflection +# I_x(a,b) = 1 - I_{1-x}(b,a), validated against scipy.special.betainc over a +# 1200-point grid of (a, b, x) to a worst absolute error of 1.1e-12. These +# cases stay to hold that fixed. # --------------------------------------------------------------------------- -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="beta_cdf_scalar is numerically incorrect") @pytest.mark.parametrize("x", [0.1, 0.3, 0.5, 0.75, 0.9]) def test_cdf_of_uniform_special_case_is_identity(x): """Beta(1, 1) is uniform on [0, 1], so its CDF is F(x) = x exactly.""" assert Beta(1.0, 1.0).cdf_scalar(x) == pytest.approx(x, **ITERATIVE) -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="beta_cdf_scalar is numerically incorrect") def test_cdf_matches_analytic_integral(): """For Beta(2, 3), F(0.5) = 0.6875 by direct integration of 12x(1-x)^2.""" assert Beta(2.0, 3.0).cdf_scalar(0.5) == pytest.approx(0.6875, **ITERATIVE) -# The CDF is wrong everywhere, but it still happens to be bounded and monotonic -# for most parameters. It breaks both properties only for alpha = beta = 0.5, -# where it returns values as low as -0.4273. Only that case is xfailed, so the -# structural guarantees stay enforced for every other parameter pair. +# alpha = beta = 0.5 is the arcsine distribution, whose density diverges at both +# endpoints; it was the case that exposed the old implementation most sharply, +# so it stays in the property list rather than being treated as special. CDF_PROPERTY_PARAMS = [ (2.0, 3.0), (1.0, 1.0), - pytest.param( - 0.5, 0.5, - marks=[ - pytest.mark.known_bug, - pytest.mark.xfail( - strict=True, - reason="beta_cdf_scalar returns negative values for alpha=beta=0.5", - ), - ], - ), + (0.5, 0.5), (5.0, 2.0), (3.0, 3.0), + (0.01, 0.01), + (100.0, 100.0), ] @@ -200,11 +192,20 @@ def test_cdf_is_monotonic(alpha, beta): assert values == sorted(values) -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="beta_cdf_scalar is numerically incorrect") -def test_cdf_is_symmetric_for_symmetric_parameters(): - """For a == b the distribution is symmetric about 0.5, so F(0.5) == 0.5.""" - assert Beta(3.0, 3.0).cdf_scalar(0.5) == pytest.approx(0.5, **ITERATIVE) +@pytest.mark.parametrize("alpha, beta", [(0.5, 0.5), (3.0, 3.0), (0.01, 0.01), (100.0, 100.0)]) +def test_cdf_is_symmetric_for_symmetric_parameters(alpha, beta): + """I_x(a,a) = 1/2 at x = 1/2 for any a, by symmetry of the density.""" + assert Beta(alpha, beta).cdf_scalar(0.5) == pytest.approx(0.5, **ITERATIVE) + + +@pytest.mark.parametrize("alpha, beta", CDF_PROPERTY_PARAMS) +def test_cdf_matches_reference_implementation(alpha, beta): + """Against the independent Lentz reference in conftest, not a stored value.""" + dist = Beta(alpha, beta) + for x in (0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99): + assert dist.cdf_scalar(x) == pytest.approx( + regularized_incomplete_beta(alpha, beta, x), **ITERATIVE + ) # --------------------------------------------------------------------------- diff --git a/tests/python/test_version.py b/tests/python/test_version.py index f899aa1..a1534b9 100644 --- a/tests/python/test_version.py +++ b/tests/python/test_version.py @@ -1,11 +1,18 @@ -import re +import re from pathlib import Path +import pytest + import fastdist import fastdist._fastdist as core CMAKELISTS = Path(__file__).resolve().parents[2] / "CMakeLists.txt" +# What __init__.py reports when importlib.metadata cannot find an installed +# distribution -- i.e. when the suite is running against a source checkout +# rather than an installed wheel. +UNINSTALLED = "0.0.0+unknown" + def _cmake_version() -> str: """The single source of truth: the project() call in CMakeLists.txt.""" @@ -23,5 +30,19 @@ def test_compiled_module_matches_cmake(): def test_python_package_matches_compiled_module(): - """The wheel metadata must match the C++ constant.""" - assert fastdist.__version__ == core.__version__ \ No newline at end of file + """The installed wheel's metadata must match the C++ constant. + + Only meaningful against an installed distribution. Running from a source + checkout there is no metadata to read, and a stale egg-info left over from + an older build reports whatever version it was generated at -- neither says + anything about the code under test, so both are skipped rather than failed. + """ + if fastdist.__version__ == UNINSTALLED: + pytest.skip("fastdist is not installed; no distribution metadata to check") + + assert fastdist.__version__ == core.__version__, ( + f"installed metadata says {fastdist.__version__} but the compiled module " + f"says {core.__version__}. If a stale python/fastdist.egg-info is present, " + f"delete it and rebuild -- importlib.metadata will read it in preference " + f"to nothing, even when it predates the current version." + ) From 3f0a23bd67aac77b511c52cf1d90aebfaaabf6db Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:51:45 -0400 Subject: [PATCH 08/15] Widen seed() to 64 bits and mix it through seed_seq 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 --- include/fastdist/math/rng.h | 6 +++++- src/bindings/rng.cpp | 10 ++++++++-- src/math/rng.cpp | 9 ++++++++- tests/python/test_rng.py | 39 +++++++++++++++++++++++++++++++++++++ 4 files changed, 60 insertions(+), 4 deletions(-) diff --git a/include/fastdist/math/rng.h b/include/fastdist/math/rng.h index a3efafa..da21e13 100644 --- a/include/fastdist/math/rng.h +++ b/include/fastdist/math/rng.h @@ -17,6 +17,10 @@ namespace fastdist::math { // reproduces the same draws on every run, which makes sampling-based tests // deterministic and lets callers reproduce a result exactly. // + // The value is spread across mt19937's full 19937-bit state via seed_seq + // rather than handed to the engine as a single word, so seeds differing by + // one bit give unrelated streams from the first draw. + // // Two caveats worth knowing: // // 1. This seeds the calling thread only. A thread that has not been seeded @@ -28,7 +32,7 @@ namespace fastdist::math { // different samples across libstdc++, libc++ and MSVC. A seed makes a // run reproducible on one platform and toolchain, not across all of // them. - void seed_rng(std::uint32_t value); + void seed_rng(std::uint64_t value); // Returns the calling thread's stream to non-deterministic behaviour by // drawing a fresh seed from the OS entropy source. This is the state every diff --git a/src/bindings/rng.cpp b/src/bindings/rng.cpp index 315aa36..807e86b 100644 --- a/src/bindings/rng.cpp +++ b/src/bindings/rng.cpp @@ -1,12 +1,18 @@ // pybind11 bindings for /src/math/rng.cpp #include "fastdist/math/rng.h" +#include #include "pybind11/pybind11.h" namespace py = pybind11; void bind_rng(py::module_ &m) { - m.def("seed", &fastdist::math::seed_rng, py::arg("value"), - R"pbdoc(Seed the sampling engine so draws are reproducible. + // Taken as a signed 64-bit value and reinterpreted, so the negative and + // large results of expressions like hash(x) are accepted rather than + // rejected by an unsigned overload. + m.def( + "seed", [](const std::int64_t value) { fastdist::math::seed_rng(static_cast(value)); }, + py::arg("value"), + R"pbdoc(Seed the sampling engine so draws are reproducible. Pins the calling thread's random stream to a fixed sequence: the same seed replays the same draws on every run. diff --git a/src/math/rng.cpp b/src/math/rng.cpp index 97441c7..3233c25 100644 --- a/src/math/rng.cpp +++ b/src/math/rng.cpp @@ -14,7 +14,14 @@ namespace fastdist::math { return engine; } - void seed_rng(const std::uint32_t value) { rng().seed(value); } + void seed_rng(const std::uint64_t value) { + // mt19937::seed(uint32_t) derives all 624 state words from one word by a + // simple recurrence. seed_seq exists to do this mixing properly, and it + // also lets the full 64 bits contribute rather than just the low 32. + std::seed_seq sequence{static_cast(value & 0xFFFFFFFFu), + static_cast(value >> 32)}; + rng().seed(sequence); + } void seed_rng_from_entropy() { rng().seed(std::random_device{}()); } diff --git a/tests/python/test_rng.py b/tests/python/test_rng.py index 4b07ee4..2b25697 100644 --- a/tests/python/test_rng.py +++ b/tests/python/test_rng.py @@ -6,6 +6,8 @@ and the order in which samplers consume the stream. """ +import pytest + import fastdist from fastdist import Bernoulli, Normal, Uniform @@ -73,6 +75,43 @@ def test_seed_from_entropy_escapes_the_fixed_stream(): assert entropic != seeded +@pytest.mark.parametrize("seed", [0, 1, -1, 2**32, 2**40, -(2**40), 2**63 - 1, -(2**63)]) +def test_seed_accepts_the_full_signed_64_bit_range(seed): + """Negative and large seeds are accepted, not just uint32. + + Seeding from `hash(x)` or a signed counter is ordinary usage, and both + routinely produce values outside the unsigned 32-bit range. + """ + fastdist.seed(seed) + first = [Uniform(0.0, 1.0).sample() for _ in range(4)] + + fastdist.seed(seed) + assert [Uniform(0.0, 1.0).sample() for _ in range(4)] == first + + +def test_adjacent_seeds_are_not_correlated(): + """Seeding in a loop is common, so adjacent seeds must give unrelated draws. + + mt19937 seeded from a single word mixes slowly; the implementation spreads + the seed across the full state with seed_seq to avoid that. If it regressed + to naive seeding, consecutive seeds would produce clustered first draws and + the mean gap would collapse well below the 1/3 expected of independent + uniforms. + """ + firsts = [] + for seed in range(512): + fastdist.seed(seed) + firsts.append(Uniform(0.0, 1.0).sample()) + + gaps = [abs(b - a) for a, b in zip(firsts, firsts[1:])] + mean_gap = sum(gaps) / len(gaps) + + # 1/3 is the expected absolute difference between two independent U(0,1) + # draws. The window is wide enough not to flake on 512 samples but far + # tighter than anything clustered seeding would produce. + assert 0.25 < mean_gap < 0.42, f"mean gap {mean_gap} suggests correlated seeding" + + def test_seed_is_exported_from_the_package_root(): """seed() is part of the public API, not an implementation detail.""" assert "seed" in fastdist.__all__ From 13fa38d10be504fa889c60474259bce63313b776 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:56:52 -0400 Subject: [PATCH 09/15] Fix the incomplete gamma: unsigned negation and a truncated series 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 --- include/fastdist/config.h | 19 +++- src/math/gamma.cpp | 8 +- tests/python/test_chi_square.py | 54 +++-------- tests/python/test_gamma.py | 165 ++++++++++++++++++++------------ 4 files changed, 140 insertions(+), 106 deletions(-) diff --git a/include/fastdist/config.h b/include/fastdist/config.h index d0fa182..459abe4 100644 --- a/include/fastdist/config.h +++ b/include/fastdist/config.h @@ -3,8 +3,23 @@ #ifndef CONFIG_H #define CONFIG_H -// Macros for iterations in Beta and Gamma calculations -constexpr unsigned int MAX_ITER = 100; +// Iteration ceiling for the Beta and Gamma series and continued fractions. +// +// Every loop using this exits as soon as its term falls below EPS, so the bound +// only matters for parameters that converge slowly, and raising it costs +// ordinary calls nothing. +// +// The gamma series is the binding constraint: near x = alpha it needs roughly +// sqrt(2 * alpha * ln(1/EPS)) terms, about 227 at alpha = 1000 and 683 at +// alpha = 10000. At the previous ceiling of 100 it simply stopped early and +// returned the truncated sum, so Gamma(1000, 0.5).cdf(500) was wrong by 9e-4 +// and Gamma(10000, ...) by 0.16, with no indication anything had gone wrong. +// +// 1000 covers alpha up to roughly 20000. Beyond that the result degrades +// silently again; a shape parameter that large needs a different algorithm +// (a normal approximation, or Temme's uniform asymptotic expansion) rather +// than a larger ceiling. +constexpr unsigned int MAX_ITER = 1000; constexpr double EPS = 1e-12; constexpr double FPMIN = 1e-30; diff --git a/src/math/gamma.cpp b/src/math/gamma.cpp index 8740abe..045cb10 100644 --- a/src/math/gamma.cpp +++ b/src/math/gamma.cpp @@ -115,7 +115,13 @@ namespace fastdist::math { double h = d; for (unsigned int i = 1; i <= MAX_ITER; ++i) { - const double an = -i * (i - a); + // i is converted to double *before* the negation. Written as + // -i * (i - a), the unary minus applies to the unsigned loop + // index and wraps to 2^32 - i, so the first coefficient came out + // as -2147483647.5 instead of 0.5 and the whole fraction was + // wrong -- returning probabilities above 1.0. + const double di = static_cast(i); + const double an = -di * (di - a); b += 2.0; d = an * d + b; if (std::fabs(d) < FPMIN) d = FPMIN; diff --git a/tests/python/test_chi_square.py b/tests/python/test_chi_square.py index fd61a29..72ae818 100644 --- a/tests/python/test_chi_square.py +++ b/tests/python/test_chi_square.py @@ -103,38 +103,19 @@ def test_pdf_of_two_degrees_is_exponential(x): # CDF # --------------------------------------------------------------------------- -# KNOWN BUG: chi_square_cdf_scalar delegates to the same regularized lower -# incomplete gamma as Gamma.cdf_scalar, whose continued-fraction branch is -# incorrect. ChiSquare(3.0).cdf(7.5) returns 1.000498004, a probability greater -# than one. The failing points below were determined empirically against -# conftest.regularized_lower_gamma; k = 2 is correct throughout because that -# case reduces to the exact exponential form. -CF_BUG = ( - "regularized lower incomplete gamma is incorrect in the " - "continued-fraction branch (x/2 >= k/2 + 1)" -) +# REGRESSION: chi_square_cdf_scalar delegates to the same regularized lower +# incomplete gamma as Gamma.cdf_scalar, whose continued-fraction branch used an +# unsigned loop index under a unary minus and was wrong as a result. +# ChiSquare(3.0).cdf(7.5) returned 1.000498004, a probability greater than one. +# k = 2 was correct throughout because that case reduces to the exact +# exponential form. See test_gamma.py for the full diagnosis. +# +# Validated against scipy.special.gammainc over 132 points with k from 0.5 to +# 10000: worst absolute error 3.1e-12, nothing outside [0, 1]. CHI2_X = (0.5, 1.0, 3.0, 7.5, 20.0, 50.0) -_CDF_KNOWN_BAD = { - (1.0, 3.0), (1.0, 7.5), (1.0, 20.0), - (3.0, 7.5), (3.0, 20.0), - (5.0, 7.5), (5.0, 20.0), (5.0, 50.0), - (10.0, 20.0), (10.0, 50.0), -} - - -def _cdf_case(k, x): - marks = [] - if (k, x) in _CDF_KNOWN_BAD: - marks = [ - pytest.mark.known_bug, - pytest.mark.xfail(strict=True, reason=CF_BUG), - ] - return pytest.param(k, x, marks=marks) - - -CDF_CASES = [_cdf_case(k, x) for k in DEGREES for x in CHI2_X] +CDF_CASES = [(k, x) for k in DEGREES for x in CHI2_X] @pytest.mark.parametrize("k, x", CDF_CASES) @@ -150,29 +131,20 @@ def test_cdf_exact_for_two_degrees_of_freedom(x): assert ChiSquare(2.0).cdf(x) == pytest.approx(1 - math.exp(-x / 2), **ITERATIVE) -_CDF_PROPERTY_DEGREES = [ - 1.0, - 2.0, - pytest.param(3.0, marks=[ - pytest.mark.known_bug, - pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF exceeds 1.0"), - ]), - 5.0, - 10.0, -] +_CDF_PROPERTY_DEGREES = [1.0, 2.0, 3.0, 5.0, 10.0, 100.0, 1000.0] @pytest.mark.parametrize("k", _CDF_PROPERTY_DEGREES) def test_cdf_is_bounded(k): dist = ChiSquare(k) - for x in (0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0): + for x in (0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0, 1000.0): assert 0.0 <= dist.cdf(x) <= 1.0 @pytest.mark.parametrize("k", _CDF_PROPERTY_DEGREES) def test_cdf_is_monotonic(k): dist = ChiSquare(k) - values = [dist.cdf(x) for x in (0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0)] + values = [dist.cdf(x) for x in (0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0, 1000.0)] assert values == sorted(values) diff --git a/tests/python/test_gamma.py b/tests/python/test_gamma.py index 9146036..14a75a1 100644 --- a/tests/python/test_gamma.py +++ b/tests/python/test_gamma.py @@ -115,51 +115,109 @@ def test_pdf_is_non_negative(x, alpha, theta): assert Gamma(alpha, theta).pmf_scalar(x) >= 0.0 +@pytest.mark.parametrize("theta", [0.5, 1.0, 3.0]) +@pytest.mark.parametrize("x", [0.25, 1.0, 4.0]) +def test_pdf_with_unit_shape_is_exponential(x, theta): + """Gamma(1, th) is Exponential with rate 1/th.""" + assert Gamma(1.0, theta).pmf_scalar(x) == pytest.approx( + math.exp(-x / theta) / theta, **EXACT + ) + + +@pytest.mark.parametrize("alpha, theta", PARAMS) +@pytest.mark.parametrize("x", [0.1, 2.0, 15.0]) +def test_pdf_is_non_negative(x, alpha, theta): + assert Gamma(alpha, theta).pmf_scalar(x) >= 0.0 + + # --------------------------------------------------------------------------- # CDF # -# KNOWN BUG: the regularized lower incomplete gamma used by gamma_cdf_scalar is -# correct in its series branch (x/theta < alpha+1) but wrong in the continued- -# fraction branch. Absolute errors reach 0.26, and the CDF can exceed 1.0 -- -# Gamma(1.5, 1.0).cdf_scalar(2.5) returns 1.000498004. +# REGRESSION: the regularized lower incomplete gamma was wrong in its +# continued-fraction branch (x/theta >= alpha+1). The Lentz coefficient was +# written `-i * (i - a)` with `i` an unsigned loop index, so the unary minus +# wrapped to 2^32 - i and the first coefficient came out as -2147483647.5 +# instead of 0.5. Absolute errors reached 0.26 and the CDF exceeded 1.0 -- +# Gamma(1.5, 1.0).cdf_scalar(2.5) returned 1.000498004. # -# The failing points below were determined empirically by comparing against -# conftest.regularized_lower_gamma. They are marked strict-xfail so they become -# failures the moment the backend is fixed and the markers go stale. Note that -# alpha = 1 is correct throughout, so the failure set is not simply -# "everything in the continued-fraction branch". +# Separately, MAX_ITER was 100, which silently truncated the series near +# x = alpha for large shapes: alpha = 1000 was wrong by 9e-4 and alpha = 10000 +# by 0.16. +# +# Both are fixed. Validated against scipy.special.gammainc over 784 points +# spanning alpha from 0.01 to 15000: worst absolute error 3.1e-12, nothing +# outside [0, 1]. These cases stay to hold that fixed. # --------------------------------------------------------------------------- -CF_BUG = ( - "regularized lower incomplete gamma is incorrect in the " - "continued-fraction branch (x/theta >= alpha+1)" -) - GAMMA_SHAPES = [(0.5, 1.0), (1.0, 1.0), (1.5, 1.0), (2.0, 3.0), (3.0, 1.0), (5.0, 2.0)] GAMMA_X = (0.2, 0.8, 1.5, 2.5, 4.0, 10.0, 20.0) -_CDF_KNOWN_BAD = { - (0.5, 1.0, 1.5), (0.5, 1.0, 2.5), (0.5, 1.0, 4.0), - (0.5, 1.0, 10.0), (0.5, 1.0, 20.0), - (1.5, 1.0, 2.5), (1.5, 1.0, 4.0), (1.5, 1.0, 10.0), (1.5, 1.0, 20.0), - (2.0, 3.0, 10.0), (2.0, 3.0, 20.0), - (3.0, 1.0, 4.0), (3.0, 1.0, 10.0), (3.0, 1.0, 20.0), - (5.0, 2.0, 20.0), -} +CDF_CASES = [ + pytest.param(alpha, theta, x) + for (alpha, theta) in GAMMA_SHAPES + for x in GAMMA_X +] -def _cdf_case(alpha, theta, x): - marks = [] - if (alpha, theta, x) in _CDF_KNOWN_BAD: - marks = [ - pytest.mark.known_bug, - pytest.mark.xfail(strict=True, reason=CF_BUG), - ] - return pytest.param(alpha, theta, x, marks=marks) +@pytest.mark.parametrize("alpha, theta, x", CDF_CASES) +def test_cdf_matches_reference(alpha, theta, x): + assert Gamma(alpha, theta).cdf_scalar(x) == pytest.approx( + regularized_lower_gamma(alpha, x / theta), **ITERATIVE + ) +@pytest.mark.parametrize("alpha", [200.0, 1000.0, 5000.0]) +def test_cdf_is_accurate_for_large_shape(alpha): + """Near x = alpha the series needs ~sqrt(2*alpha*ln(1/EPS)) terms. + + At the old MAX_ITER of 100 it stopped early and returned the truncated sum + with no indication anything was wrong. + """ + assert Gamma(alpha, 1.0).cdf_scalar(alpha) == pytest.approx( + regularized_lower_gamma(alpha, alpha), **ITERATIVE + ) + + +@pytest.mark.parametrize("theta", [0.5, 1.0, 3.0]) +@pytest.mark.parametrize("x", [0.25, 1.0, 4.0]) +def test_pdf_with_unit_shape_is_exponential(x, theta): + """Gamma(1, th) is Exponential with rate 1/th.""" + assert Gamma(1.0, theta).pmf_scalar(x) == pytest.approx( + math.exp(-x / theta) / theta, **EXACT + ) + + +@pytest.mark.parametrize("alpha, theta", PARAMS) +@pytest.mark.parametrize("x", [0.1, 2.0, 15.0]) +def test_pdf_is_non_negative(x, alpha, theta): + assert Gamma(alpha, theta).pmf_scalar(x) >= 0.0 + + +# --------------------------------------------------------------------------- +# CDF +# +# REGRESSION: the regularized lower incomplete gamma was wrong in its +# continued-fraction branch (x/theta >= alpha+1). The Lentz coefficient was +# written `-i * (i - a)` with `i` an unsigned loop index, so the unary minus +# wrapped to 2^32 - i: the first coefficient came out as -2147483647.5 instead +# of 0.5. Absolute errors reached 0.26 and the CDF exceeded 1.0 -- +# Gamma(1.5, 1.0).cdf_scalar(2.5) returned 1.000498004. alpha = 1 was correct +# throughout because that case reduces to the exact exponential form. +# +# Separately MAX_ITER was 100, which silently truncated the series near +# x = alpha for large shapes: alpha = 1000 was wrong by 9e-4, alpha = 10000 by +# 0.16. +# +# Both are fixed. Validated against scipy.special.gammainc over 784 points with +# alpha from 0.01 to 15000: worst absolute error 3.1e-12, nothing outside +# [0, 1]. These cases stay to hold that fixed. +# --------------------------------------------------------------------------- + +GAMMA_SHAPES = [(0.5, 1.0), (1.0, 1.0), (1.5, 1.0), (2.0, 3.0), (3.0, 1.0), (5.0, 2.0)] +GAMMA_X = (0.2, 0.8, 1.5, 2.5, 4.0, 10.0, 20.0) + CDF_CASES = [ - _cdf_case(alpha, theta, x) + (alpha, theta, x) for (alpha, theta) in GAMMA_SHAPES for x in GAMMA_X ] @@ -172,6 +230,18 @@ def test_cdf_matches_reference(alpha, theta, x): ) +@pytest.mark.parametrize("alpha", [200.0, 1000.0, 5000.0]) +def test_cdf_is_accurate_for_large_shape(alpha): + """Near x = alpha the series needs ~sqrt(2*alpha*ln(1/EPS)) terms. + + At the old MAX_ITER of 100 it stopped early and returned the truncated sum + with no indication anything had gone wrong. + """ + assert Gamma(alpha, 1.0).cdf_scalar(alpha) == pytest.approx( + regularized_lower_gamma(alpha, alpha), **ITERATIVE + ) + + @pytest.mark.parametrize("theta", [0.5, 1.0, 3.0]) @pytest.mark.parametrize("x", [0.25, 1.0, 2.0, 6.0]) def test_cdf_with_unit_shape_is_exponential(x, theta): @@ -181,43 +251,14 @@ def test_cdf_with_unit_shape_is_exponential(x, theta): ) -_BOUNDED_SHAPES = [ - (0.5, 1.0), - (1.0, 1.0), - pytest.param(1.5, 1.0, marks=[ - pytest.mark.known_bug, - pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF exceeds 1.0"), - ]), - pytest.param(2.0, 3.0, marks=[ - pytest.mark.known_bug, - pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF exceeds 1.0"), - ]), - (3.0, 1.0), - (5.0, 2.0), -] - - -@pytest.mark.parametrize("alpha, theta", _BOUNDED_SHAPES) +@pytest.mark.parametrize("alpha, theta", GAMMA_SHAPES) def test_cdf_is_bounded(alpha, theta): dist = Gamma(alpha, theta) for x in (0.1, 0.5, 1.0, 3.0, 8.0, 25.0): assert 0.0 <= dist.cdf_scalar(x) <= 1.0 -_MONOTONIC_SHAPES = [ - (0.5, 1.0), - (1.0, 1.0), - pytest.param(1.5, 1.0, marks=[ - pytest.mark.known_bug, - pytest.mark.xfail(strict=True, reason=CF_BUG + "; CDF is non-monotonic"), - ]), - (2.0, 3.0), - (3.0, 1.0), - (5.0, 2.0), -] - - -@pytest.mark.parametrize("alpha, theta", _MONOTONIC_SHAPES) +@pytest.mark.parametrize("alpha, theta", GAMMA_SHAPES) def test_cdf_is_monotonic(alpha, theta): dist = Gamma(alpha, theta) values = [dist.cdf_scalar(x) for x in (0.1, 0.5, 1.0, 3.0, 8.0, 25.0)] From 7b9e5edc83043d737d318fe6b5353a5848760600 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 02:59:00 -0400 Subject: [PATCH 10/15] Fix five defects in the Python distribution classes 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 --- python/fastdist/distributions/beta.py | 8 ++++---- python/fastdist/distributions/binomial.py | 4 ++-- .../distributions/discrete_uniform.py | 12 +++++++---- python/fastdist/distributions/uniform.py | 8 ++++++-- tests/python/test_beta.py | 20 +++++++------------ tests/python/test_binomial.py | 12 +++-------- tests/python/test_discrete_uniform.py | 15 +++++--------- tests/python/test_uniform.py | 17 +++++----------- 8 files changed, 40 insertions(+), 56 deletions(-) diff --git a/python/fastdist/distributions/beta.py b/python/fastdist/distributions/beta.py index 54e9bfa..7cc0d4d 100644 --- a/python/fastdist/distributions/beta.py +++ b/python/fastdist/distributions/beta.py @@ -33,12 +33,12 @@ def beta(self): @alpha.setter def alpha(self, value): - self._validate_params(alpha=value) + self._validate_params(alpha=value, beta=self._beta) self._alpha = float(value) @beta.setter def beta(self, value): - self._validate_params(beta=value) + self._validate_params(alpha=self._alpha, beta=value) self._beta = float(value) def __repr__(self): @@ -50,8 +50,8 @@ def _validate_params(alpha: Union[int, float] = None, beta: Union[int, float] = if alpha is not None: if not isinstance(alpha, (int, float)): raise TypeError("alpha must be a real number") - if alpha <= 0: - raise ValueError("alpha must be positive") + if alpha <= 0: + raise ValueError("alpha must be positive") if beta is not None: if not isinstance(beta, (int, float)): diff --git a/python/fastdist/distributions/binomial.py b/python/fastdist/distributions/binomial.py index f290e94..152d163 100644 --- a/python/fastdist/distributions/binomial.py +++ b/python/fastdist/distributions/binomial.py @@ -50,8 +50,8 @@ def _validate_params(n: int = None, p: Real = None) -> None: if n is not None: if not isinstance(n, int): raise TypeError("n must be an integer") - if n < 0: - raise ValueError("n must be a non-negative integer") + if n < 0: + raise ValueError("n must be a non-negative integer") if p is not None: if not isinstance(p, (int, float)): diff --git a/python/fastdist/distributions/discrete_uniform.py b/python/fastdist/distributions/discrete_uniform.py index 64a3164..fcf2c6d 100644 --- a/python/fastdist/distributions/discrete_uniform.py +++ b/python/fastdist/distributions/discrete_uniform.py @@ -33,13 +33,17 @@ def b(self): @a.setter def a(self, value): - self._validate_params(a=value) - self._a = float(value) + # Both bounds are passed so the a < b relationship is re-checked against + # the current opposite bound, and int() matches how __init__ stores it -- + # a is an integer parameter, so assigning through the setter must not + # quietly change its type to float. + self._validate_params(a=value, b=self._b) + self._a = int(value) @b.setter def b(self, value): - self._validate_params(b=value) - self.b = value + self._validate_params(a=self._a, b=value) + self._b = int(value) def __repr__(self): return f"DiscreteUniform(a={self.a}, b={self.b})" diff --git a/python/fastdist/distributions/uniform.py b/python/fastdist/distributions/uniform.py index 281fd66..da06827 100644 --- a/python/fastdist/distributions/uniform.py +++ b/python/fastdist/distributions/uniform.py @@ -121,7 +121,11 @@ def a(self, value): 0.2 """ - self._validate_params(a=value) + # The opposite bound is passed too: validating `a` alone skips the + # a < b check entirely, which let Uniform(1.0, 3.0) be driven to + # a = 10.0, b = -10.0 -- a state the constructor rejects outright, and + # from which pdf, cdf, mean, variance and sample all silently return nan. + self._validate_params(a=value, b=self._b) self._a = float(value) @b.setter @@ -149,7 +153,7 @@ def b(self, value): 2.0 """ - self._validate_params(b=value) + self._validate_params(a=self._a, b=value) self._b = float(value) def __repr__(self): diff --git a/tests/python/test_beta.py b/tests/python/test_beta.py index 52c1f36..4b12932 100644 --- a/tests/python/test_beta.py +++ b/tests/python/test_beta.py @@ -67,20 +67,17 @@ def test_alpha_setter_updates_and_validates(): dist.alpha = -1.0 -# KNOWN BUG: the beta setter calls _validate_params(beta=value), leaving alpha -# as None. Because the `if alpha <= 0` check sits outside the `if alpha is not -# None` guard (beta.py line 47), the comparison None <= 0 raises TypeError for -# *every* assignment, valid or not. The beta property is unusable. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="beta setter raises TypeError for any value") +# REGRESSION: the beta setter called _validate_params(beta=value), leaving +# alpha as None, and the `if alpha <= 0` check sat outside the `if alpha is not +# None` guard -- so None <= 0 raised TypeError for every assignment, valid or +# not, and the property was unusable. The check is now inside its guard and the +# setters pass both parameters. def test_beta_setter_updates_value(): dist = Beta(alpha=2.0, beta=3.0) dist.beta = 5.0 assert dist.beta == 5.0 -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="beta setter raises TypeError before validating") def test_beta_setter_rejects_non_positive(): dist = Beta(alpha=2.0, beta=3.0) with pytest.raises(ValueError, match="beta must be positive"): @@ -244,13 +241,10 @@ def test_classmethods_reject_invalid_parameters(method_name, args): # --------------------------------------------------------------------------- # Validation edge case # -# KNOWN BUG: in Beta._validate_params the `if alpha <= 0` check sits outside the -# `if alpha is not None` guard (beta.py line 47), so validating only `beta` -# raises TypeError comparing None to int. +# REGRESSION: the `if alpha <= 0` check sat outside the `if alpha is not None` +# guard, so validating only `beta` raised TypeError comparing None to int. # --------------------------------------------------------------------------- -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="alpha <= 0 check sits outside the None guard") def test_validate_params_accepts_a_single_named_parameter(): Beta._validate_params(beta=3.0) diff --git a/tests/python/test_binomial.py b/tests/python/test_binomial.py index 9daac6f..0863c9b 100644 --- a/tests/python/test_binomial.py +++ b/tests/python/test_binomial.py @@ -85,21 +85,15 @@ def test_n_setter_updates_and_validates(): dist.n = -5 -# KNOWN BUG: the p setter calls _validate_params(p=value), leaving n as None. -# Because the `if n < 0` check sits outside the `if n is not None` guard -# (binomial.py line 47), the comparison None < 0 raises TypeError for *every* -# assignment, valid or not. The p property is unusable. This is the same defect -# pattern as Beta._validate_params. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="p setter raises TypeError for any value") +# REGRESSION: the `if n < 0` check sat outside the `if n is not None` guard, so +# setting p (which leaves n as None) raised TypeError comparing None to int for +# every assignment. Same defect pattern as Beta._validate_params had. def test_p_setter_updates_value(): dist = Binomial(n=10, p=0.3) dist.p = 0.6 assert dist.p == 0.6 -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="p setter raises TypeError before validating") def test_p_setter_rejects_out_of_range(): dist = Binomial(n=10, p=0.3) with pytest.raises(ValueError, match=r"p must be in the interval \[0, 1\]"): diff --git a/tests/python/test_discrete_uniform.py b/tests/python/test_discrete_uniform.py index 8181a6f..73127e0 100644 --- a/tests/python/test_discrete_uniform.py +++ b/tests/python/test_discrete_uniform.py @@ -69,22 +69,17 @@ def test_repr(): assert repr(DiscreteUniform(a=1, b=6)) == "DiscreteUniform(a=1, b=6)" -# KNOWN BUG: the b setter assigns `self.b = value` instead of `self._b = value` -# (discrete_uniform.py), so it re-enters itself and raises RecursionError for -# every assignment. The b property is unusable. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="b setter recurses into itself (self.b = value)") +# REGRESSION: the b setter assigned `self.b = value` instead of `self._b`, so +# it re-entered itself and raised RecursionError for every assignment. def test_b_setter_updates_value(): dist = DiscreteUniform(a=1, b=6) dist.b = 10 assert dist.b == 10 -# KNOWN BUG: the a setter stores float(value) even though a is an integer -# parameter that __init__ stores via int(). Setting a therefore changes the -# attribute's type from int to float. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="a setter stores float(value) instead of int") +# REGRESSION: the a setter stored float(value) even though a is an integer +# parameter that __init__ stores via int(), so assigning through the setter +# silently changed the attribute's type. def test_a_setter_preserves_integer_type(): dist = DiscreteUniform(a=1, b=6) dist.a = 2 diff --git a/tests/python/test_uniform.py b/tests/python/test_uniform.py index 10a0e3e..20acea7 100644 --- a/tests/python/test_uniform.py +++ b/tests/python/test_uniform.py @@ -84,32 +84,25 @@ def test_property_setters_update_values(): assert dist.b == 5.0 -# KNOWN BUG: each setter validates only the bound being assigned, never the -# a < b relationship against the other one. Uniform(1.0, 3.0) can be driven to -# a = 10.0, b = -10.0 -- a state the constructor rejects outright. Once there, -# pdf, cdf, mean, variance and sample all return nan rather than raising, so the -# corruption propagates silently. -CROSS_BOUND_BUG = "setters do not re-validate a < b against the other bound" +# REGRESSION: each setter validated only the bound being assigned, never the +# a < b relationship against the other one. 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 rather than raising, so +# the corruption propagated silently. Both setters now pass the opposite bound. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason=CROSS_BOUND_BUG) def test_a_setter_rejects_value_above_b(): dist = Uniform(a=1.0, b=3.0) with pytest.raises(ValueError, match="a must be less than b"): dist.a = 10.0 -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason=CROSS_BOUND_BUG) def test_b_setter_rejects_value_below_a(): dist = Uniform(a=1.0, b=3.0) with pytest.raises(ValueError, match="a must be less than b"): dist.b = -10.0 -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason=CROSS_BOUND_BUG + "; results become nan") def test_instance_stays_usable_after_setter_assignments(): dist = Uniform(a=1.0, b=3.0) try: From 82342668c537740896832f339c8db2a3be0071fe Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 03:00:25 -0400 Subject: [PATCH 11/15] Resolve the last two annotation/implementation mismatches in Utils 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 --- python/fastdist/distributions/utils.py | 29 +++++++++++++++++++++++- tests/python/test_utils.py | 31 +++++++++++++++----------- 2 files changed, 46 insertions(+), 14 deletions(-) diff --git a/python/fastdist/distributions/utils.py b/python/fastdist/distributions/utils.py index c50812a..515b1e0 100644 --- a/python/fastdist/distributions/utils.py +++ b/python/fastdist/distributions/utils.py @@ -109,10 +109,37 @@ def law_of_total_probability(cls, p_A: Union[Real, Sequence[Real]], p_A_valid = p_A_valid.tolist() if isinstance(p_B_given_A_valid, np.ndarray): p_B_given_A_valid = p_B_given_A_valid.tolist() + + # The binding takes vectors. A scalar is the one-element partition + # P(B) = P(B|A) P(A), which the signature already advertises, so promote + # rather than letting it fail inside pybind11 with an argument-type error. + if isinstance(p_A_valid, Real): + p_A_valid = [p_A_valid] + if isinstance(p_B_given_A_valid, Real): + p_B_given_A_valid = [p_B_given_A_valid] + + if len(p_A_valid) != len(p_B_given_A_valid): + raise ValueError("p_A and p_B_given_A must have the same length") + return _core.law_of_total_probability(p_B_given_A_valid, p_A_valid) @classmethod - def sigmoid(cls, x: Union[Real, Sequence[Real]]) -> float: + def sigmoid(cls, x: Real) -> float: + """Logistic function for a single value. + + Scalar only. The signature used to advertise Sequence[Real] as well, but + the body calls float() on the input so any sequence raised TypeError. + Use sigmoid_cpu for arrays -- the 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. + """ + # _validate_input accepts a sequence even when asked for Real, and + # float() on the resulting array then fails with a numpy message about + # 0-dimensional arrays, which says nothing useful. Reject it here with + # the name of the function that does handle arrays. + if not isinstance(x, Real): + raise TypeError("x must be a real number; use Utils.sigmoid_cpu for arrays") + validated_input = cls._validate_input(_input=x, input_name="x", input_type=Real) return _core.sigmoid(float(validated_input)) diff --git a/tests/python/test_utils.py b/tests/python/test_utils.py index 7f0bcff..f500ce6 100644 --- a/tests/python/test_utils.py +++ b/tests/python/test_utils.py @@ -69,12 +69,11 @@ def test_law_of_total_probability_over_a_partition_is_bounded(): assert 0.0 <= result <= 1.0 -# KNOWN BUG: the signature annotates both arguments as -# Union[Real, Sequence[Real]], but the implementation validates them as -# sequences and forwards them to a vector-only binding, so scalar inputs raise -# TypeError from pybind11. Either the annotation or the implementation is wrong. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="scalar inputs are annotated but unsupported") +# REGRESSION: the signature annotates both arguments as +# Union[Real, Sequence[Real]] but they were forwarded straight to a vector-only +# binding, so scalar inputs raised TypeError from pybind11. A scalar is the +# one-element partition P(B) = P(B|A) P(A), so it is promoted rather than +# rejected and the signature now tells the truth. def test_law_of_total_probability_accepts_scalars(): assert Utils.law_of_total_probability(0.3, 0.2) == pytest.approx(0.06, **EXACT) @@ -133,13 +132,19 @@ def test_logit_returns_nan_outside_the_open_unit_interval(p): assert math.isnan(Utils.logit(p)) -# KNOWN BUG: sigmoid is annotated Union[Real, Sequence[Real]] but its body calls -# float() on the validated input, so any sequence raises TypeError. The array -# path exists separately as sigmoid_cpu. -@pytest.mark.known_bug -@pytest.mark.xfail(strict=True, reason="sequence inputs are annotated but unsupported") -def test_sigmoid_accepts_a_sequence(): - result = Utils.sigmoid([0.0, 2.0]) +# RESOLVED: sigmoid was annotated Union[Real, Sequence[Real]] while its body +# called float() on the input, so a sequence raised an opaque numpy error. The +# annotation was the wrong half: the library's convention is an explicit *_cpu +# entry point for arrays, and returning an ndarray from a function annotated +# -> float would be worse than not accepting one. sigmoid is now scalar-only +# and says so, pointing at sigmoid_cpu. +def test_sigmoid_rejects_a_sequence_and_names_the_array_path(): + with pytest.raises(TypeError, match="sigmoid_cpu"): + Utils.sigmoid([0.0, 2.0]) + + +def test_sigmoid_cpu_handles_what_sigmoid_rejects(): + result = Utils.sigmoid_cpu([0.0, 2.0]) np.testing.assert_allclose(result, [0.5, 1.0 / (1.0 + math.exp(-2.0))], rtol=1e-12) From 6377d3e1957f6c7e7b785660e4dcbb2c4b28f509 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 03:04:06 -0400 Subject: [PATCH 12/15] Wire CUDA into the benchmark suite and document the toolchain constraint 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 --- README.md | 53 +++++++++++++++++++++++++++++++++++++++++++++++ benchmarks/run.py | 45 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 98 insertions(+) diff --git a/README.md b/README.md index 810786b..117051b 100644 --- a/README.md +++ b/README.md @@ -74,6 +74,53 @@ cmake -S . -B build -Dpybind11_DIR="$(python -m pybind11 --cmakedir)" Pass `-DPython_EXECUTABLE=...` as well if the interpreter you want is not the first one on `PATH`. +### Building with CUDA + +```bash +cmake -S . -B build -DFASTDIST_ENABLE_CUDA=ON -DCMAKE_BUILD_TYPE=Release +cmake --build build --parallel +``` + +`Normal.is_cuda_available()` reports whether the extension was built with the +backend; the Python classes dispatch to it automatically above the per-function +thresholds in `fastdist.config`. + +**The CUDA toolkit and the MSVC toolset have to be compatible versions.** This +is the failure most likely to stop you, and it does not announce itself +clearly: `nvcc` aborts with + +``` +nvcc error : 'cudafe++' died with status 0xC0000409 +``` + +on every `.cu` file. That is `cudafe++` crashing on standard library headers +from an MSVC newer than the toolkit supports. `-allow-unsupported-compiler` +silences the version *check* but does not fix the incompatibility. + +Each CUDA release supports host compilers up to a specific MSVC version -- +CUDA 12.4 tops out at MSVC 19.39 (Visual Studio 17.9), so MSVC 19.42 +(VS 17.12) fails. Check with `nvcc --version` and `cl` (in a developer prompt), +then either: + +- upgrade the CUDA toolkit to one that supports your MSVC, or +- install the older MSVC toolset alongside the current one through the Visual + Studio Installer (Individual components -> "MSVC v143 ... build tools + (v14.39)") and point CMake at it with `-T version=14.39`. + +Note also that the CUDA toolkit installs its Visual Studio MSBuild integration +into whichever VS instance it finds. If you have both Build Tools and a full +Visual Studio install, the integration may land in the one without the C++ +workload, and the Visual Studio generator will then report "No CUDA toolset +found" even though `nvcc` is on `PATH`. The Ninja generator does not use that +integration at all: + +```bash +cmake -S . -B build -G Ninja -DFASTDIST_ENABLE_CUDA=ON -DCMAKE_BUILD_TYPE=Release +``` + +Ninja needs the MSVC environment on `PATH`, so run it from a "x64 Native Tools +Command Prompt" (or after `vcvars64.bat`). + --- ## Creating A Distribution Class @@ -211,6 +258,12 @@ python benchmarks/run.py --quick # one array size, for a fast check python benchmarks/compare.py --latest # diff the two newest reports ``` +When the extension is built with `FASTDIST_ENABLE_CUDA=ON`, the suite also +times the `*_cuda` entry points against this library's own `*_cpu` path, so the +reported speedup is the one a caller actually decides on. GPU timings include +the host-to-device copy and the copy back, since a caller cannot avoid those. +The CUDA cases are skipped entirely on a CPU-only build. + `compare.py` exits non-zero if any case regressed by more than 5%, so it can gate a change. Every case checks that fastdist and the baseline agree numerically before either is timed -- a speedup on a wrong answer is not a diff --git a/benchmarks/run.py b/benchmarks/run.py index 08824ec..f271b99 100644 --- a/benchmarks/run.py +++ b/benchmarks/run.py @@ -28,6 +28,13 @@ whichever library has the thinner binding layer and should not be quoted as a throughput number. + cuda The *_cuda entry points against the *_cpu ones on the same input, + so the number is the speedup the GPU backend buys over this library's + own CPU path -- the decision a caller actually faces. Skipped + entirely unless the extension was built with FASTDIST_ENABLE_CUDA. + GPU timings include the host-to-device copy and the copy back, + because a caller cannot avoid those. + sample Drawing variates. fastdist samples one value per call, while numpy fills an array in one call, so numpy is expected to win by a wide margin. It is measured anyway: this is a real gap in the library and @@ -156,6 +163,36 @@ def sample_cases(sizes): lambda n=n, rng=rng: rng.uniform(0.0, 1.0, n)) +# --------------------------------------------------------------------------- +# CUDA: *_cuda against this library's own *_cpu path +# --------------------------------------------------------------------------- +def cuda_cases(sizes): + """Empty unless the extension was built with the CUDA backend.""" + if not hasattr(core, "normal_pdf_cuda"): + return + + rng = np.random.default_rng(31337) + for n in sizes: + x_real = rng.normal(0.0, 1.0, n) + x_pos = np.abs(rng.normal(2.0, 1.0, n)) + 0.05 + + yield ("normal_pdf", n, + lambda x=x_real: core.normal_pdf_cuda(x, 0.0, 1.0, 0.0), + lambda x=x_real: core.normal_pdf_cpu(x, 0.0, 1.0, 0.0)) + yield ("normal_cdf", n, + lambda x=x_real: core.normal_cdf_cuda(x, 0.0, 1.0, 0.0), + lambda x=x_real: core.normal_cdf_cpu(x, 0.0, 1.0, 0.0)) + yield ("normal_logpdf", n, + lambda x=x_real: core.normal_logpdf_cuda(x, 0.0, 1.0), + lambda x=x_real: core.normal_logpdf_cpu(x, 0.0, 1.0, 0.0)) + yield ("exponential_pdf", n, + lambda x=x_pos: core.exponential_pdf_cuda(x, 2.0, 0.0), + lambda x=x_pos: core.exponential_pdf_cpu(x, 2.0, 0.0)) + yield ("uniform_pdf", n, + lambda x=x_real: core.uniform_pdf_cuda(x, -3.0, 3.0, 0.0), + lambda x=x_real: core.uniform_pdf_cpu(x, -3.0, 3.0, 0.0)) + + def run(sizes, sample_sizes) -> list[Result]: results: list[Result] = [] @@ -170,6 +207,14 @@ def run(sizes, sample_sizes) -> list[Result]: results.append(measure("scalar", case, n, fd, sp, "scipy")) print(f" scalar {case:<18} n={n:<9,} {_fmt(results[-1])}") + for case, n, gpu, cpu in cuda_cases(sizes): + # The "fastdist" column is the GPU path and the baseline is the CPU + # path, so `speedup` reads as "how much the GPU buys over the CPU". + # measure() also checks the two agree numerically, which is the part + # worth having: a kernel that is fast and wrong is the failure mode. + results.append(measure("cuda", case, n, gpu, cpu, "fastdist-cpu")) + print(f" cuda {case:<18} n={n:<9,} {_fmt(results[-1])}") + for case, n, fd, np_fn in sample_cases(sample_sizes): # No correctness check: both draw from their own RNG, so the outputs # are different random numbers by construction. measure() would flag From 01c4568cdd669159ac45d5df0981b9d4918a3aa7 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 03:05:39 -0400 Subject: [PATCH 13/15] Correct comments that were wrong, stale, or unreadable 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 --- src/bindings/negative_binomial.cpp | 2 +- src/bindings/uniform.cpp | 2 +- src/math/bernoulli.cpp | 2 +- src/math/beta.cpp | 2 +- src/math/binomial.cpp | 2 +- src/math/chi_square.cpp | 2 +- src/math/discrete_uniform.cpp | 2 +- src/math/exponential.cpp | 2 +- src/math/gamma.cpp | 10 ++++++++-- src/math/geometric.cpp | 2 +- src/math/negative_binomial.cpp | 2 +- src/math/normal.cpp | 2 +- src/math/poisson.cpp | 15 ++++++--------- src/math/uniform.cpp | 2 +- 14 files changed, 26 insertions(+), 23 deletions(-) diff --git a/src/bindings/negative_binomial.cpp b/src/bindings/negative_binomial.cpp index 153557e..58dc864 100644 --- a/src/bindings/negative_binomial.cpp +++ b/src/bindings/negative_binomial.cpp @@ -1,4 +1,4 @@ -// pybind11 bindings for /src/math/bernoulli.cpp +// pybind11 bindings for /src/math/negative_binomial.cpp #include "fastdist/math/negative_binomial.h" #include diff --git a/src/bindings/uniform.cpp b/src/bindings/uniform.cpp index 721efc9..ae38000 100644 --- a/src/bindings/uniform.cpp +++ b/src/bindings/uniform.cpp @@ -1,4 +1,4 @@ -// pybind11 bindings for /src/math/normal.cpp +// pybind11 bindings for /src/math/uniform.cpp #include "fastdist/math/uniform.h" #include "pybind11/pybind11.h" #include "wrappers/uniform_wrapper.h" diff --git a/src/math/bernoulli.cpp b/src/math/bernoulli.cpp index 863400b..a63da08 100644 --- a/src/math/bernoulli.cpp +++ b/src/math/bernoulli.cpp @@ -1,4 +1,4 @@ -// Function declarations for Bernoulli distribution functions +// Function definitions for Bernoulli distribution functions #include #include #include diff --git a/src/math/beta.cpp b/src/math/beta.cpp index 330c22e..c3c6d60 100644 --- a/src/math/beta.cpp +++ b/src/math/beta.cpp @@ -1,4 +1,4 @@ -// Function declarations for beta distribution functions +// Function definitions for beta distribution functions #include #include #include diff --git a/src/math/binomial.cpp b/src/math/binomial.cpp index 45f8af8..d2be9bb 100644 --- a/src/math/binomial.cpp +++ b/src/math/binomial.cpp @@ -1,4 +1,4 @@ -// Function declarations for binomial distribution functions +// Function definitions for binomial distribution functions #include #include #include diff --git a/src/math/chi_square.cpp b/src/math/chi_square.cpp index 8fa4e18..4a4d836 100644 --- a/src/math/chi_square.cpp +++ b/src/math/chi_square.cpp @@ -1,4 +1,4 @@ -// Function declarations for chi-square distribution functions +// Function definitions for chi-square distribution functions #include #include #include diff --git a/src/math/discrete_uniform.cpp b/src/math/discrete_uniform.cpp index 56f037c..4fd9978 100644 --- a/src/math/discrete_uniform.cpp +++ b/src/math/discrete_uniform.cpp @@ -1,4 +1,4 @@ -// Function declarations for discrete uniform distribution functions +// Function definitions for discrete uniform distribution functions #include #include #include diff --git a/src/math/exponential.cpp b/src/math/exponential.cpp index 87e9dfd..8d99866 100644 --- a/src/math/exponential.cpp +++ b/src/math/exponential.cpp @@ -1,4 +1,4 @@ -// Function declarations for exponential distribution functions +// Function definitions for exponential distribution functions #include "fastdist/math/exponential.h" #include #include diff --git a/src/math/gamma.cpp b/src/math/gamma.cpp index 045cb10..d379075 100644 --- a/src/math/gamma.cpp +++ b/src/math/gamma.cpp @@ -1,4 +1,4 @@ -// Function declarations for gamma distribution functions +// Function definitions for gamma distribution functions #include #include #include @@ -27,7 +27,13 @@ namespace fastdist::math { static double gamma_p_cf(double a, double x); // ------------------------- - // CDF using series / continued fraction + // CDF + // + // The regularized lower incomplete gamma P(a, x). The series converges + // quickly below x = a+1 and the continued fraction above it, so the + // dispatch below picks whichever is on its fast side. Both are evaluated + // in log space, since x^a and Gamma(a) each overflow well before their + // ratio does. // ------------------------- double gamma_cdf_scalar(const double x, const double alpha, const double theta) { if (!std::isfinite(x) || !std::isfinite(alpha) || !std::isfinite(theta) || x < 0.0 || alpha <= 0.0 || diff --git a/src/math/geometric.cpp b/src/math/geometric.cpp index 7eba2f6..9435a06 100644 --- a/src/math/geometric.cpp +++ b/src/math/geometric.cpp @@ -1,4 +1,4 @@ -// Function declarations for geometric distribution functions +// Function definitions for geometric distribution functions #include #include #include diff --git a/src/math/negative_binomial.cpp b/src/math/negative_binomial.cpp index 47c2a2b..4cccb47 100644 --- a/src/math/negative_binomial.cpp +++ b/src/math/negative_binomial.cpp @@ -1,4 +1,4 @@ -// Function declarations for negative binomial distribution functions +// Function definitions for negative binomial distribution functions #include #include #include diff --git a/src/math/normal.cpp b/src/math/normal.cpp index f145b0b..850bfbb 100644 --- a/src/math/normal.cpp +++ b/src/math/normal.cpp @@ -1,4 +1,4 @@ -// Function declarations for normal distribution functions +// Function definitions for normal distribution functions #include #include #include diff --git a/src/math/poisson.cpp b/src/math/poisson.cpp index c953d4b..670b4d7 100644 --- a/src/math/poisson.cpp +++ b/src/math/poisson.cpp @@ -1,4 +1,4 @@ -// Function declarations for poisson distribution functions +// Function definitions for poisson distribution functions #include #include #include @@ -12,18 +12,15 @@ namespace fastdist::math { return std::numeric_limits::quiet_NaN(); } - // Poisson is defined - // on non-negative - // integers + // Poisson is defined on non-negative integers if (x < 0.0 || std::floor(x) != x) { return 0.0; } - // log PMF for - // numerical - // stability: log P = - // k * log(lambda) - - // lambda - log(k!) + // Evaluated in log space for numerical stability: + // log P(k) = k log(lambda) - lambda - log(k!) + // Forming lambda^k / k! directly overflows both terms for modest k + // even where their ratio is an ordinary number. const double log_p = x * std::log(lambda) - lambda - std::lgamma(x + 1.0); return std::exp(log_p); diff --git a/src/math/uniform.cpp b/src/math/uniform.cpp index cc69ef7..e18af33 100644 --- a/src/math/uniform.cpp +++ b/src/math/uniform.cpp @@ -1,4 +1,4 @@ -// Function declarations for continuous uniform distribution functions +// Function definitions for continuous uniform distribution functions #include #include #include From 8f0f46fdd0221ac2a772c26e1c508b217d50fdc5 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 03:06:35 -0400 Subject: [PATCH 14/15] Benchmark the three iterative CDFs 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 --- benchmarks/run.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/benchmarks/run.py b/benchmarks/run.py index f271b99..c193647 100644 --- a/benchmarks/run.py +++ b/benchmarks/run.py @@ -139,6 +139,17 @@ def fd_loop(fn, *args): def sp_loop(fn, *args, **kw): return lambda: [float(fn(float(v), *args, **kw)) for v in xs] + # Positive support for the distributions defined on it, and the unit + # interval for beta. + xs_pos = np.abs(np.random.default_rng(8).normal(3.0, 1.5, SCALAR_N)) + 0.01 + xs_unit = np.random.default_rng(9).uniform(0.01, 0.99, SCALAR_N) + + def fd_loop_on(values, fn, *args): + return lambda: [fn(float(v), *args) for v in values] + + def sp_loop_on(values, fn, *args, **kw): + return lambda: [float(fn(float(v), *args, **kw)) for v in values] + return [ ("normal_pdf", SCALAR_N, fd_loop(core.normal_pdf_scalar, 0.0, 1.0), @@ -146,6 +157,18 @@ def sp_loop(fn, *args, **kw): ("normal_cdf", SCALAR_N, fd_loop(core.normal_cdf_scalar, 0.0, 1.0), sp_loop(sps.norm.cdf, 0.0, 1.0)), + # The three iterative CDFs. They have no *_cpu batch path, so scalar is + # the only way to track them -- and they are the most expensive + # routines in the library, so a regression here matters most. + ("gamma_cdf", SCALAR_N, + fd_loop_on(xs_pos, core.gamma_cdf_scalar, 3.0, 2.0), + sp_loop_on(xs_pos, sps.gamma.cdf, 3.0, 0.0, 2.0)), + ("chi_square_cdf", SCALAR_N, + fd_loop_on(xs_pos, core.chi_square_cdf_scalar, 6.0), + sp_loop_on(xs_pos, sps.chi2.cdf, 6.0)), + ("beta_cdf", SCALAR_N, + fd_loop_on(xs_unit, core.beta_cdf_scalar, 2.0, 5.0), + sp_loop_on(xs_unit, sps.beta.cdf, 2.0, 5.0)), ] From 97d4dae5e10a8659f98075a75df7609144581372 Mon Sep 17 00:00:00 2001 From: ghosteau Date: Sat, 5 Sep 2026 03:09:47 -0400 Subject: [PATCH 15/15] Make compare.py noise-aware and record the correctness pass 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 --- BENCHMARKS.md | 99 ++++ benchmarks/compare.py | 34 +- ...0.1.0_20260905T070637+0000_8f0f46fdd0.json | 486 ++++++++++++++++++ ...0.1.0_20260905T070807+0000_8f0f46fdd0.json | 486 ++++++++++++++++++ benchmarks/run.py | 7 +- 5 files changed, 1104 insertions(+), 8 deletions(-) create mode 100644 benchmarks/results/0.1.0_20260905T070637+0000_8f0f46fdd0.json create mode 100644 benchmarks/results/0.1.0_20260905T070807+0000_8f0f46fdd0.json diff --git a/BENCHMARKS.md b/BENCHMARKS.md index 4448616..4eda950 100644 --- a/BENCHMARKS.md +++ b/BENCHMARKS.md @@ -270,6 +270,105 @@ The remaining known gaps are unchanged and still worth recording: --- +## Unreleased — correctness pass + +The headline of this entry is not a speedup. Three of the library's CDFs were +returning wrong answers, two of them probabilities above 1.0, and the benchmark +suite is what surfaced the first of them: it checks agreement with SciPy before +it times anything, so a wrong result cannot quietly post a good number. + +- **beta_cdf** was wrong at every point -- 0.0015 against a true 0.1143 for + Beta(2,5) at x=0.1, and -147 for Beta(0.01,0.01) at x=0.5. The series had an + inverted coefficient ratio and normalised by Gamma(a+1) instead of B(a,b). + Replaced with the modified-Lentz continued fraction. +- **gamma_cdf** and **chi_square_cdf** shared a continued fraction whose Lentz + coefficient was written `-i * (i - a)` with an unsigned loop index, so the + unary minus wrapped to 2^32 - i. Errors reached 0.26 and + Gamma(1.5,1.0).cdf(2.5) returned 1.000498004. +- **MAX_ITER** was 100, silently truncating the gamma series for large shapes + (wrong by 0.16 at alpha = 10000). Now 1000. + +All three now agree with SciPy to ~1e-12 across the parameter ranges recorded in +the commits, and they are benchmarked from here on -- their absence from the +suite is why the defects survived this long. They land at 46-56x SciPy in the +scalar group, which is where their cost can be tracked since they have no +`*_cpu` batch path. + +Performance is otherwise unchanged from the previous entry: no regression +outside measurement noise. + +One methodology change came out of this run. The scalar cases are dominated by +per-call Python overhead, which the interpreter varies far more than it varies +compiled work; at 7 rounds an untouched `normal_cdf` differed by 8% between two +runs and `compare.py` reported it as a regression. The scalar group now uses 21 +rounds, and `compare.py` will not flag a change smaller than the two runs' +combined `noise_pct`. Re-running confirmed the phantom: `normal_cdf` came back +at -8.7%, its original level. + + +- **Version** 0.1.0 (`8f0f46fdd0` on `fix/flaky-rng-tolerances`, working tree dirty) +- **Measured** 2026-09-05T07:08:07+00:00 +- **CPU** AMD Ryzen 7 7700 8-Core Processor +- **Platform** Windows-11-10.0.26200-SP0 +- **Toolchain** Python 3.14.2, numpy 2.5.2, scipy 1.18.1 +- **CUDA** not built + +### batch (vs vectorised SciPy) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_pdf` | 1,000 | 5.22 us | 31.30 us (scipy) | **6.00x** | 1.1e-16 | +| `normal_cdf` | 1,000 | 6.01 us | 30.26 us (scipy) | **5.04x** | 2.2e-16 | +| `normal_logpdf` | 1,000 | 1.88 us | 31.97 us (scipy) | **17.04x** | 8.9e-16 | +| `exponential_pdf` | 1,000 | 4.17 us | 29.55 us (scipy) | **7.09x** | 0.0e+00 | +| `exponential_cdf` | 1,000 | 4.20 us | 30.39 us (scipy) | **7.24x** | 8.3e-17 | +| `uniform_pdf` | 1,000 | 2.04 us | 33.67 us (scipy) | **16.50x** | 0.0e+00 | +| `uniform_cdf` | 1,000 | 2.10 us | 32.46 us (scipy) | **15.44x** | 0.0e+00 | +| `poisson_pmf` | 1,000 | 40.67 us | 41.04 us (scipy) | **1.01x** | 2.0e-19 | +| `poisson_cdf` | 1,000 | 9.94 us | 72.03 us (scipy) | **7.25x** | 2.2e-16 | +| `bernoulli_pmf` | 1,000 | 2.00 us | 49.44 us (scipy) | **24.72x** | 2.2e-16 | +| `normal_pdf` | 100,000 | 414.70 us | 1.44 ms (scipy) | **3.47x** | 1.1e-16 | +| `normal_cdf` | 100,000 | 529.20 us | 2.09 ms (scipy) | **3.95x** | 2.2e-16 | +| `normal_logpdf` | 100,000 | 77.60 us | 1.66 ms (scipy) | **21.45x** | 8.9e-16 | +| `exponential_pdf` | 100,000 | 308.80 us | 1.54 ms (scipy) | **4.99x** | 0.0e+00 | +| `exponential_cdf` | 100,000 | 312.30 us | 1.64 ms (scipy) | **5.26x** | 1.1e-16 | +| `uniform_pdf` | 100,000 | 93.70 us | 1.51 ms (scipy) | **16.07x** | 0.0e+00 | +| `uniform_cdf` | 100,000 | 99.30 us | 1.59 ms (scipy) | **16.04x** | 0.0e+00 | +| `poisson_pmf` | 100,000 | 4.02 ms | 3.23 ms (scipy) | **0.80x** | 2.0e-19 | +| `poisson_cdf` | 100,000 | 1.31 ms | 6.18 ms (scipy) | **4.70x** | 2.2e-16 | +| `bernoulli_pmf` | 100,000 | 299.20 us | 3.54 ms (scipy) | **11.82x** | 2.2e-16 | +| `normal_pdf` | 1,000,000 | 4.82 ms | 20.38 ms (scipy) | **4.23x** | 1.1e-16 | +| `normal_cdf` | 1,000,000 | 5.96 ms | 22.55 ms (scipy) | **3.78x** | 2.2e-16 | +| `normal_logpdf` | 1,000,000 | 1.31 ms | 21.51 ms (scipy) | **16.39x** | 8.9e-16 | +| `exponential_pdf` | 1,000,000 | 3.69 ms | 17.35 ms (scipy) | **4.70x** | 0.0e+00 | +| `exponential_cdf` | 1,000,000 | 3.80 ms | 19.79 ms (scipy) | **5.20x** | 1.7e-16 | +| `uniform_pdf` | 1,000,000 | 1.39 ms | 18.56 ms (scipy) | **13.36x** | 0.0e+00 | +| `uniform_cdf` | 1,000,000 | 1.49 ms | 19.02 ms (scipy) | **12.80x** | 0.0e+00 | +| `poisson_pmf` | 1,000,000 | 42.78 ms | 38.78 ms (scipy) | **0.91x** | 2.0e-19 | +| `poisson_cdf` | 1,000,000 | 14.15 ms | 63.90 ms (scipy) | **4.52x** | 2.2e-16 | +| `bernoulli_pmf` | 1,000,000 | 3.51 ms | 38.80 ms (scipy) | **11.06x** | 2.2e-16 | + +### scalar (per-call cost, not throughput) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_pdf` | 20,000 | 7.22 ms | 441.28 ms (scipy) | **61.08x** | 1.1e-16 | +| `normal_cdf` | 20,000 | 7.24 ms | 425.74 ms (scipy) | **58.77x** | 2.2e-16 | +| `gamma_cdf` | 20,000 | 8.78 ms | 426.74 ms (scipy) | **48.63x** | 1.2e-13 | +| `chi_square_cdf` | 20,000 | 7.91 ms | 445.09 ms (scipy) | **56.25x** | 1.2e-13 | +| `beta_cdf` | 20,000 | 9.48 ms | 467.80 ms (scipy) | **49.36x** | 8.9e-16 | + +### sample (vs numpy) + +| case | n | fastdist | baseline | speedup | max abs diff | +|---|---:|---:|---:|---:|---:| +| `normal_sample` | 100,000 | 28.22 ms | 893.60 us (numpy) | **0.03x** | - | +| `uniform_sample` | 100,000 | 23.78 ms | 226.90 us (numpy) | **0.01x** | - | +| `normal_sample` | 1,000,000 | 295.14 ms | 9.94 ms (numpy) | **0.03x** | - | +| `uniform_sample` | 1,000,000 | 258.23 ms | 3.14 ms (numpy) | **0.01x** | - | + +--- + ## Changes to record here Add an entry when a release ships, or when a change is made specifically to diff --git a/benchmarks/compare.py b/benchmarks/compare.py index fd459ca..9edc68b 100644 --- a/benchmarks/compare.py +++ b/benchmarks/compare.py @@ -17,10 +17,17 @@ change there can come from either side, so a moved speedup with an unmoved `change` means the baseline moved, not this library. -The threshold for calling something a regression defaults to 5%, which is -comfortably above the noise on a quiet machine. Check the `noise_pct` field in -the reports before trusting a smaller difference: if either run was noisy, the -comparison is not meaningful at that resolution. +The threshold for calling something a regression defaults to 5%, but a case is +only flagged when the change also exceeds the two runs' combined `noise_pct` +(the gap between each run's minimum and median round). A change smaller than +the noise says nothing, and the scalar cases -- dominated by Python call +overhead -- routinely swing 10% between runs on an otherwise idle machine. +Changes above the threshold but inside the noise are printed and marked rather +than counted. + +`noise_pct` measures spread *within* a run, so it does not catch a run that was +uniformly slow. Treat a flagged case as a prompt to re-run rather than a +verdict; if it does not reproduce, it was the machine. """ from __future__ import annotations @@ -50,7 +57,9 @@ def main() -> int: parser.add_argument("after", nargs="?", type=Path) parser.add_argument("--latest", action="store_true", help="use the newest two reports") parser.add_argument("--threshold", type=float, default=5.0, - help="percent change before a case is called out (default 5)") + help="minimum percent change before a case is called out " + "(default 5); the run's measured noise raises this " + "further when the machine was busy") args = parser.parse_args() if args.latest: @@ -95,11 +104,22 @@ def main() -> int: speed = f"{b['speedup']:6.2f}x -> {a['speedup']:6.2f}x" label = f"{group}/{case} n={n:,}" + + # A run is only as trustworthy as it was quiet. Each report records + # noise_pct -- the gap between the minimum and median round -- and a + # change smaller than the noise in the two runs combined says nothing. + # Without this the scalar cases, which are dominated by Python call + # overhead and routinely swing 10%, produce phantom regressions. + noise = (b.get("fastdist_noise_pct") or 0.0) + (a.get("fastdist_noise_pct") or 0.0) + limit = max(args.threshold, noise) + flag = "" - if change > args.threshold: + if change > limit: flag, _ = " REGRESSED", regressions.append((label, change)) - elif change < -args.threshold: + elif change < -limit: flag, _ = " faster", improvements.append((label, change)) + elif abs(change) > args.threshold: + flag = f" (within noise, +/-{noise:.0f}%)" print(f"{label:<34} {_fmt_time(b['fastdist_s']):>11} {_fmt_time(a['fastdist_s']):>11} " f"{change:+8.1f}% {speed:>16}{flag}") diff --git a/benchmarks/results/0.1.0_20260905T070637+0000_8f0f46fdd0.json b/benchmarks/results/0.1.0_20260905T070637+0000_8f0f46fdd0.json new file mode 100644 index 0000000..9a64cd9 --- /dev/null +++ b/benchmarks/results/0.1.0_20260905T070637+0000_8f0f46fdd0.json @@ -0,0 +1,486 @@ +{ + "environment": { + "timestamp_utc": "2026-09-05T07:06:37+00:00", + "fastdist_version": "0.1.0", + "git_commit": "8f0f46fdd0221ac2a772c26e1c508b217d50fdc5", + "git_branch": "fix/flaky-rng-tolerances", + "git_dirty": true, + "cuda_available": false, + "python": "3.14.2", + "numpy": "2.5.2", + "scipy": "1.18.1", + "platform": "Windows-11-10.0.26200-SP0", + "processor": "AMD Ryzen 7 7700 8-Core Processor", + "machine": "AMD64" + }, + "results": [ + { + "group": "batch", + "case": "normal_pdf", + "n": 1000, + "fastdist_s": 5.243999767117202e-06, + "baseline_s": 3.0984000186435875e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 9.344017170150746, + "baseline_noise_pct": 7.40382129985198, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 5.908467117165566 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000, + "fastdist_s": 6.031999946571887e-06, + "baseline_s": 2.9643999878317116e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.19894050862418441, + "baseline_noise_pct": 4.034544530745808, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.914456256778389 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000, + "fastdist_s": 1.906000543385744e-06, + "baseline_s": 3.142200002912432e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.20983458767855567, + "baseline_noise_pct": 5.703011443412737, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 16.485829523063785 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000, + "fastdist_s": 4.1919999057427046e-06, + "baseline_s": 2.888199989683926e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.38550051577377, + "baseline_noise_pct": 1.4057204575563182, + "max_abs_diff": 0.0, + "speedup": 6.889790206644144 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000, + "fastdist_s": 4.2560003930702805e-06, + "baseline_s": 2.997600007802248e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.37592900088391296, + "baseline_noise_pct": 0.520414468992185, + "max_abs_diff": 8.326672684688674e-17, + "speedup": 7.043232450549136 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000, + "fastdist_s": 2.0680000307038427e-06, + "baseline_s": 3.320000017993152e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.2243882777428565, + "baseline_noise_pct": 21.493974648039124, + "max_abs_diff": 0.0, + "speedup": 16.0541584559996 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000, + "fastdist_s": 2.1299999207258226e-06, + "baseline_s": 3.1925999792292714e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.09389743208559582, + "baseline_noise_pct": 7.273069703884072, + "max_abs_diff": 0.0, + "speedup": 14.988732854700554 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000, + "fastdist_s": 4.0530000114813446e-05, + "baseline_s": 3.820000041741878e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.08758795237858, + "baseline_noise_pct": 29.450260190527082, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.9425117273428512 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000, + "fastdist_s": 9.954000124707819e-06, + "baseline_s": 7.066399964969605e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.22101835687388482, + "baseline_noise_pct": 12.433488887920745, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 7.099055531885505 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000, + "fastdist_s": 2.004000125452876e-06, + "baseline_s": 4.9125999794341624e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.6985790244466411, + "baseline_noise_pct": 4.547490412223295, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 24.513970418659447 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 100000, + "fastdist_s": 0.00041459998465143144, + "baseline_s": 0.0014866999990772456, + "baseline_name": "scipy", + "fastdist_noise_pct": 7.428847765550876, + "baseline_noise_pct": 4.573888071487637, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 3.5858660253621712 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 100000, + "fastdist_s": 0.0005287000094540417, + "baseline_s": 0.0021744999976363033, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.4917705026215466, + "baseline_noise_pct": 2.120027676812596, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.11291840127218 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 100000, + "fastdist_s": 7.830001413822174e-05, + "baseline_s": 0.0016923999937716872, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.8939669219972881, + "baseline_noise_pct": 3.456630176417811, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 21.614299976806148 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 100000, + "fastdist_s": 0.0003093999985139817, + "baseline_s": 0.0015406000020448118, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.424046198145162, + "baseline_noise_pct": 4.303518758132793, + "max_abs_diff": 0.0, + "speedup": 4.979314833368341 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 100000, + "fastdist_s": 0.00031209998996928334, + "baseline_s": 0.0017124000005424023, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.12816503000743962, + "baseline_noise_pct": 14.786265303475247, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 5.486703157891596 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 100000, + "fastdist_s": 9.379998664371669e-05, + "baseline_s": 0.0016820999735500664, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.10661063506242587, + "baseline_noise_pct": 1.3613952329132646, + "max_abs_diff": 0.0, + "speedup": 17.93283809238947 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 100000, + "fastdist_s": 0.00010030000703409314, + "baseline_s": 0.0016482999781146646, + "baseline_name": "scipy", + "fastdist_noise_pct": 63.70886073745561, + "baseline_noise_pct": 8.699873626197833, + "max_abs_diff": 0.0, + "speedup": 16.433697532587292 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 100000, + "fastdist_s": 0.004086000000825152, + "baseline_s": 0.0032562999986112118, + "baseline_name": "scipy", + "fastdist_noise_pct": 9.476260615926657, + "baseline_noise_pct": 3.2275891478311203, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.7969407728716629 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 100000, + "fastdist_s": 0.0013057999894954264, + "baseline_s": 0.0061944999906700104, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.9572681858043112, + "baseline_noise_pct": 3.013963932834326, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.743835227831196 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 100000, + "fastdist_s": 0.00030079999123699963, + "baseline_s": 0.0035217000113334507, + "baseline_name": "scipy", + "fastdist_noise_pct": 23.636974598683032, + "baseline_noise_pct": 2.2602718765671654, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 11.707779634071569 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 1000000, + "fastdist_s": 0.004805000004125759, + "baseline_s": 0.02025329999742098, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.8054111052936427, + "baseline_noise_pct": 1.6851575231154283, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 4.215046822066746 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000000, + "fastdist_s": 0.005918100010603666, + "baseline_s": 0.022419999993871897, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.3774521259423325, + "baseline_noise_pct": 0.9322034784986444, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.7883780189082983 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000000, + "fastdist_s": 0.0012892999802716076, + "baseline_s": 0.021157800016226247, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.095247382562667, + "baseline_noise_pct": 2.972898821720191, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 16.410300426568753 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000000, + "fastdist_s": 0.003711299999849871, + "baseline_s": 0.01773330001742579, + "baseline_name": "scipy", + "fastdist_noise_pct": 6.808934964564618, + "baseline_noise_pct": 3.8650447441297393, + "max_abs_diff": 0.0, + "speedup": 4.778190935290366 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000000, + "fastdist_s": 0.003899499977706, + "baseline_s": 0.019318000006023794, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.1724586125728145, + "baseline_noise_pct": 1.921523833775782, + "max_abs_diff": 1.6653345369377348e-16, + "speedup": 4.953968487361858 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000000, + "fastdist_s": 0.0014127000176813453, + "baseline_s": 0.01817789999768138, + "baseline_name": "scipy", + "fastdist_noise_pct": 7.397181248629095, + "baseline_noise_pct": 2.419421396144653, + "max_abs_diff": 0.0, + "speedup": 12.86748762664889 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000000, + "fastdist_s": 0.0014718000020366162, + "baseline_s": 0.019254600018030033, + "baseline_name": "scipy", + "fastdist_noise_pct": 3.607826933176652, + "baseline_noise_pct": 3.5300655139699106, + "max_abs_diff": 0.0, + "speedup": 13.082348139275927 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000000, + "fastdist_s": 0.04324629998882301, + "baseline_s": 0.03819900000235066, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.1973278960312919, + "baseline_noise_pct": 2.45582340782798, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.8832894377605292 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000000, + "fastdist_s": 0.014147600013529882, + "baseline_s": 0.0646400999976322, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.3945827472622407, + "baseline_noise_pct": 0.5515152481566593, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.568979893113634 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000000, + "fastdist_s": 0.003634300024714321, + "baseline_s": 0.03918700001668185, + "baseline_name": "scipy", + "fastdist_noise_pct": 9.157196703822052, + "baseline_noise_pct": 5.2499552848056865, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 10.782544025038824 + }, + { + "group": "scalar", + "case": "normal_pdf", + "n": 20000, + "fastdist_s": 0.00728070002514869, + "baseline_s": 0.44155819999286905, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.9253504959495436, + "baseline_noise_pct": 0.6323968137963056, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 60.647767174537776 + }, + { + "group": "scalar", + "case": "normal_cdf", + "n": 20000, + "fastdist_s": 0.007933199987746775, + "baseline_s": 0.41405970000778325, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.482428663221429, + "baseline_noise_pct": 2.2169508381686462, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 52.19327643918207 + }, + { + "group": "scalar", + "case": "gamma_cdf", + "n": 20000, + "fastdist_s": 0.0086775999807287, + "baseline_s": 0.4146003999921959, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.816124302226725, + "baseline_noise_pct": 1.4974177566050075, + "max_abs_diff": 1.177946629127291e-13, + "speedup": 47.77823371818758 + }, + { + "group": "scalar", + "case": "chi_square_cdf", + "n": 20000, + "fastdist_s": 0.008192099980078638, + "baseline_s": 0.4177081000234466, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.412531956447161, + "baseline_noise_pct": 3.7904220560601867, + "max_abs_diff": 1.177946629127291e-13, + "speedup": 50.98913600166254 + }, + { + "group": "scalar", + "case": "beta_cdf", + "n": 20000, + "fastdist_s": 0.010213500005193055, + "baseline_s": 0.4712754999927711, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.845253837890077, + "baseline_noise_pct": 1.428888203239127, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 46.14240953181094 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 100000, + "fastdist_s": 0.0271978999953717, + "baseline_s": 0.0008704999927431345, + "baseline_name": "numpy", + "fastdist_noise_pct": 2.5704190301836105, + "baseline_noise_pct": 1.7001752448466325, + "max_abs_diff": null, + "speedup": 0.03200614727207867 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 100000, + "fastdist_s": 0.02320339999278076, + "baseline_s": 0.0002272999845445156, + "baseline_name": "numpy", + "fastdist_noise_pct": 5.532809833899264, + "baseline_noise_pct": 0.08800292139481238, + "max_abs_diff": null, + "speedup": 0.009795977512572947 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 1000000, + "fastdist_s": 0.28518350000376813, + "baseline_s": 0.010094299999764189, + "baseline_name": "numpy", + "fastdist_noise_pct": 1.5593819367165542, + "baseline_noise_pct": 4.571887030660729, + "max_abs_diff": null, + "speedup": 0.035395806558341604 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 1000000, + "fastdist_s": 0.24961210001492873, + "baseline_s": 0.003158200008329004, + "baseline_name": "numpy", + "fastdist_noise_pct": 1.0643314058091122, + "baseline_noise_pct": 2.3716036999370944, + "max_abs_diff": null, + "speedup": 0.012652431545346237 + } + ] +} diff --git a/benchmarks/results/0.1.0_20260905T070807+0000_8f0f46fdd0.json b/benchmarks/results/0.1.0_20260905T070807+0000_8f0f46fdd0.json new file mode 100644 index 0000000..ee8ffed --- /dev/null +++ b/benchmarks/results/0.1.0_20260905T070807+0000_8f0f46fdd0.json @@ -0,0 +1,486 @@ +{ + "environment": { + "timestamp_utc": "2026-09-05T07:08:07+00:00", + "fastdist_version": "0.1.0", + "git_commit": "8f0f46fdd0221ac2a772c26e1c508b217d50fdc5", + "git_branch": "fix/flaky-rng-tolerances", + "git_dirty": true, + "cuda_available": false, + "python": "3.14.2", + "numpy": "2.5.2", + "scipy": "1.18.1", + "platform": "Windows-11-10.0.26200-SP0", + "processor": "AMD Ryzen 7 7700 8-Core Processor", + "machine": "AMD64" + }, + "results": [ + { + "group": "batch", + "case": "normal_pdf", + "n": 1000, + "fastdist_s": 5.215999553911388e-06, + "baseline_s": 3.129800024908036e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.07669886650167344, + "baseline_noise_pct": 1.7445188080875789, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 6.000383996507541 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000, + "fastdist_s": 6.008000345900655e-06, + "baseline_s": 3.0255999881774187e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.19972550942841416, + "baseline_noise_pct": 2.135114571900112, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 5.035951754300129 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000, + "fastdist_s": 1.8759997328743339e-06, + "baseline_s": 3.1969999545253815e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.31986293272949645, + "baseline_noise_pct": 3.528309215266958, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 17.04158000932688 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000, + "fastdist_s": 4.166000289842486e-06, + "baseline_s": 2.955000032670796e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.608730254316881, + "baseline_noise_pct": 5.482231850703793, + "max_abs_diff": 0.0, + "speedup": 7.093134486513736 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000, + "fastdist_s": 4.199999966658652e-06, + "baseline_s": 3.039000032003969e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.28571646355196106, + "baseline_noise_pct": 19.921024903780676, + "max_abs_diff": 8.326672684688674e-17, + "speedup": 7.235714419354324 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000, + "fastdist_s": 2.0399998174980283e-06, + "baseline_s": 3.366800025105476e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.19607994195212294, + "baseline_noise_pct": 1.847450079860052, + "max_abs_diff": 0.0, + "speedup": 16.50392316816337 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000, + "fastdist_s": 2.102000289596617e-06, + "baseline_s": 3.246000036597252e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.28541687839171587, + "baseline_noise_pct": 3.9679602515201995, + "max_abs_diff": 0.0, + "speedup": 15.442433822024704 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000, + "fastdist_s": 4.0668000001460316e-05, + "baseline_s": 4.1037999908439813e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.3622507717939523, + "baseline_noise_pct": 64.01384048198858, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 1.00909806007097 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000, + "fastdist_s": 9.935999987646937e-06, + "baseline_s": 7.203200017102063e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.2077328080896443, + "baseline_noise_pct": 4.484118363959397, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 7.249597449735846 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000, + "fastdist_s": 2.0000000949949025e-06, + "baseline_s": 4.943799984175712e-05, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.2000015133991153, + "baseline_noise_pct": 10.878271047615565, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 24.718998746789122 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 100000, + "fastdist_s": 0.00041470001451671124, + "baseline_s": 0.0014407999988179654, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.1687909818693038, + "baseline_noise_pct": 22.008606274158044, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 3.474318660193597 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 100000, + "fastdist_s": 0.0005292000132612884, + "baseline_s": 0.002089199988404289, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.6991628634252743, + "baseline_noise_pct": 6.121961237659482, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.947845684147333 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 100000, + "fastdist_s": 7.759997970424592e-05, + "baseline_s": 0.0016646000149194151, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.3196058980278025, + "baseline_noise_pct": 6.8785292393572535, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 21.451036730468832 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 100000, + "fastdist_s": 0.0003087999939452857, + "baseline_s": 0.0015406000020448118, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.09715100072222363, + "baseline_noise_pct": 8.256522060527768, + "max_abs_diff": 0.0, + "speedup": 4.988989741747797 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 100000, + "fastdist_s": 0.000312299991492182, + "baseline_s": 0.0016438000020571053, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.913868502124682, + "baseline_noise_pct": 5.274363710912244, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 5.263528808319726 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 100000, + "fastdist_s": 9.369998588226736e-05, + "baseline_s": 0.001505599997472018, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.21347988993345876, + "baseline_noise_pct": 10.082358659787163, + "max_abs_diff": 0.0, + "speedup": 16.0683054890081 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 100000, + "fastdist_s": 9.929999941959977e-05, + "baseline_s": 0.0015924000181257725, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.503528509737327, + "baseline_noise_pct": 3.7427755827857387, + "max_abs_diff": 0.0, + "speedup": 16.03625405270109 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 100000, + "fastdist_s": 0.0040179000061471015, + "baseline_s": 0.0032318999874405563, + "baseline_name": "scipy", + "fastdist_noise_pct": 4.726349349417224, + "baseline_noise_pct": 4.857824559310516, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.8043754156390102 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 100000, + "fastdist_s": 0.0013130000152159482, + "baseline_s": 0.006176499999128282, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.1751704785110773, + "baseline_noise_pct": 4.617501641811325, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.704112663785794 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 100000, + "fastdist_s": 0.0002992000081576407, + "baseline_s": 0.00353660000837408, + "baseline_name": "scipy", + "fastdist_noise_pct": 7.219247747711818, + "baseline_noise_pct": 3.1018489803508493, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 11.82018687148811 + }, + { + "group": "batch", + "case": "normal_pdf", + "n": 1000000, + "fastdist_s": 0.004816299973754212, + "baseline_s": 0.020377599983476102, + "baseline_name": "scipy", + "fastdist_noise_pct": 1.337126529971661, + "baseline_noise_pct": 2.947353979125874, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 4.23096569867349 + }, + { + "group": "batch", + "case": "normal_cdf", + "n": 1000000, + "fastdist_s": 0.0059634000062942505, + "baseline_s": 0.022553299990249798, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.030720573242678, + "baseline_noise_pct": 1.881764536672156, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 3.7819532425202462 + }, + { + "group": "batch", + "case": "normal_logpdf", + "n": 1000000, + "fastdist_s": 0.0013123000098858029, + "baseline_s": 0.02150829997844994, + "baseline_name": "scipy", + "fastdist_noise_pct": 3.040462041878866, + "baseline_noise_pct": 0.9884557047448108, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 16.38977353990998 + }, + { + "group": "batch", + "case": "exponential_pdf", + "n": 1000000, + "fastdist_s": 0.0036939000128768384, + "baseline_s": 0.017347600020002574, + "baseline_name": "scipy", + "fastdist_noise_pct": 12.133517585894355, + "baseline_noise_pct": 3.5987685903535342, + "max_abs_diff": 0.0, + "speedup": 4.6962830503070725 + }, + { + "group": "batch", + "case": "exponential_cdf", + "n": 1000000, + "fastdist_s": 0.0038022999942768365, + "baseline_s": 0.019786299992119893, + "baseline_name": "scipy", + "fastdist_noise_pct": 6.724877746574977, + "baseline_noise_pct": 2.2798603115570915, + "max_abs_diff": 1.6653345369377348e-16, + "speedup": 5.203771407280311 + }, + { + "group": "batch", + "case": "uniform_pdf", + "n": 1000000, + "fastdist_s": 0.0013889999827370048, + "baseline_s": 0.018563900026492774, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.11159234287221, + "baseline_noise_pct": 2.0405193598675533, + "max_abs_diff": 0.0, + "speedup": 13.364938990073185 + }, + { + "group": "batch", + "case": "uniform_cdf", + "n": 1000000, + "fastdist_s": 0.001486199995269999, + "baseline_s": 0.019019999977899715, + "baseline_name": "scipy", + "fastdist_noise_pct": 7.630197196174329, + "baseline_noise_pct": 4.63932716179416, + "max_abs_diff": 0.0, + "speedup": 12.797739226505877 + }, + { + "group": "batch", + "case": "poisson_pmf", + "n": 1000000, + "fastdist_s": 0.042776000016601756, + "baseline_s": 0.03877519999514334, + "baseline_name": "scipy", + "fastdist_noise_pct": 0.5238918767416122, + "baseline_noise_pct": 0.864470051174673, + "max_abs_diff": 1.9651164376299768e-19, + "speedup": 0.9064709178065815 + }, + { + "group": "batch", + "case": "poisson_cdf", + "n": 1000000, + "fastdist_s": 0.01414640000439249, + "baseline_s": 0.06390479998663068, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.0648363697598042, + "baseline_noise_pct": 3.065653936484622, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 4.517389580867789 + }, + { + "group": "batch", + "case": "bernoulli_pmf", + "n": 1000000, + "fastdist_s": 0.0035083999973721802, + "baseline_s": 0.038803799980087206, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.3800021544997296, + "baseline_noise_pct": 0.7883248576790491, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 11.0602553896795 + }, + { + "group": "scalar", + "case": "normal_pdf", + "n": 20000, + "fastdist_s": 0.007224400003906339, + "baseline_s": 0.4412822999875061, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.2866951769042276, + "baseline_noise_pct": 1.417097407235027, + "max_abs_diff": 1.1102230246251565e-16, + "speedup": 61.082207484206066 + }, + { + "group": "scalar", + "case": "normal_cdf", + "n": 20000, + "fastdist_s": 0.0072438000061083585, + "baseline_s": 0.4257381999923382, + "baseline_name": "scipy", + "fastdist_noise_pct": 2.645020618789027, + "baseline_noise_pct": 3.7515543640784395, + "max_abs_diff": 2.220446049250313e-16, + "speedup": 58.77277114681977 + }, + { + "group": "scalar", + "case": "gamma_cdf", + "n": 20000, + "fastdist_s": 0.008775800000876188, + "baseline_s": 0.4267418999806978, + "baseline_name": "scipy", + "fastdist_noise_pct": 3.625880012970469, + "baseline_noise_pct": 2.405950771697355, + "max_abs_diff": 1.177946629127291e-13, + "speedup": 48.62712230658074 + }, + { + "group": "scalar", + "case": "chi_square_cdf", + "n": 20000, + "fastdist_s": 0.007912200002465397, + "baseline_s": 0.44508679999853484, + "baseline_name": "scipy", + "fastdist_noise_pct": 5.624225666220858, + "baseline_noise_pct": 1.007151862223517, + "max_abs_diff": 1.177946629127291e-13, + "speedup": 56.25322917265088 + }, + { + "group": "scalar", + "case": "beta_cdf", + "n": 20000, + "fastdist_s": 0.009476500010350719, + "baseline_s": 0.467804799991427, + "baseline_name": "scipy", + "fastdist_noise_pct": 6.900226894401768, + "baseline_noise_pct": 2.2513663850786214, + "max_abs_diff": 8.881784197001252e-16, + "speedup": 49.36472320798466 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 100000, + "fastdist_s": 0.028224000008776784, + "baseline_s": 0.0008936000231187791, + "baseline_name": "numpy", + "fastdist_noise_pct": 9.085175721429552, + "baseline_noise_pct": 0.358101519668714, + "max_abs_diff": null, + "speedup": 0.03166099854169846 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 100000, + "fastdist_s": 0.023780899995472282, + "baseline_s": 0.00022689998149871826, + "baseline_name": "numpy", + "fastdist_noise_pct": 3.8652869088133985, + "baseline_noise_pct": 0.0881452354370526, + "max_abs_diff": null, + "speedup": 0.00954126973924109 + }, + { + "group": "sample", + "case": "normal_sample", + "n": 1000000, + "fastdist_s": 0.2951441000041086, + "baseline_s": 0.009936599992215633, + "baseline_name": "numpy", + "fastdist_noise_pct": 0.3300421711928733, + "baseline_noise_pct": 2.979892365443683, + "max_abs_diff": null, + "speedup": 0.03366694435727263 + }, + { + "group": "sample", + "case": "uniform_sample", + "n": 1000000, + "fastdist_s": 0.2582290999998804, + "baseline_s": 0.003140300017548725, + "baseline_name": "numpy", + "fastdist_noise_pct": 0.869344307332093, + "baseline_noise_pct": 2.5124984607259133, + "max_abs_diff": null, + "speedup": 0.012160906797685386 + } + ] +} diff --git a/benchmarks/run.py b/benchmarks/run.py index c193647..fb6d758 100644 --- a/benchmarks/run.py +++ b/benchmarks/run.py @@ -227,7 +227,12 @@ def run(sizes, sample_sizes) -> list[Result]: print(f" batch {case:<18} n={n:<9,} {_fmt(results[-1])}") for case, n, fd, sp in scalar_cases(): - results.append(measure("scalar", case, n, fd, sp, "scipy")) + # More rounds than the batch cases get. These 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 + # near the truth. At 7 rounds two runs of an unchanged normal_cdf + # differed by 8%, which is enough to look like a regression. + results.append(measure("scalar", case, n, fd, sp, "scipy", repeat=21)) print(f" scalar {case:<18} n={n:<9,} {_fmt(results[-1])}") for case, n, gpu, cpu in cuda_cases(sizes):