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/BENCHMARKS.md b/BENCHMARKS.md new file mode 100644 index 0000000..4eda950 --- /dev/null +++ b/BENCHMARKS.md @@ -0,0 +1,377 @@ +# 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. + +--- + +## 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. + +--- + +## 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 +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/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/README.md b/README.md index 7ca48c2..117051b 100644 --- a/README.md +++ b/README.md @@ -21,28 +21,106 @@ 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 +cmake -S . -B build -DCMAKE_BUILD_TYPE=Release +cmake --build build --target fastdist_tests --parallel +ctest --test-dir build --output-on-failure +``` + +On Windows the Visual Studio generator is multi-config, so `CMAKE_BUILD_TYPE` +is ignored -- pick the configuration at build and test time instead: + +```powershell +cmake -S . -B build +cmake --build build --target fastdist_tests --config Release --parallel +ctest --test-dir build -C Release --output-on-failure +``` + +**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 -python3 setup.py bdist_wheel +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. -**Important:** +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: -- This command must be run from the project root. +- 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`. -3. Install the generated wheel: +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 -pip install .\dist\fastdist--cpXXX-cpXXX-win_amd64.whl --force-reinstall +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 @@ -124,6 +202,106 @@ 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 +``` + +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 +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 +398,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/benchmarks/compare.py b/benchmarks/compare.py new file mode 100644 index 0000000..9edc68b --- /dev/null +++ b/benchmarks/compare.py @@ -0,0 +1,156 @@ +""" +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%, 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 + +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="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: + 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:,}" + + # 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 > limit: + flag, _ = " REGRESSED", regressions.append((label, change)) + 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}") + + 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/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 + } + ] +} 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 new file mode 100644 index 0000000..fb6d758 --- /dev/null +++ b/benchmarks/run.py @@ -0,0 +1,301 @@ +""" +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. + + 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 + 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] + + # 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), + 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)), + # 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)), + ] + + +# --------------------------------------------------------------------------- +# 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)) + + +# --------------------------------------------------------------------------- +# 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] = [] + + 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(): + # 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): + # 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 + # 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()) 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/include/fastdist/math/rng.h b/include/fastdist/math/rng.h new file mode 100644 index 0000000..da21e13 --- /dev/null +++ b/include/fastdist/math/rng.h @@ -0,0 +1,43 @@ +// 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. + // + // 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 + // 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::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 + // 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/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..7cc0d4d 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 @@ -28,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): @@ -45,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 df55b70..152d163 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 @@ -45,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/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..fcf2c6d 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 @@ -28,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/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..da06827 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 @@ -115,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 @@ -143,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/python/fastdist/distributions/utils.py b/python/fastdist/distributions/utils.py index c1d4f49..515b1e0 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 @@ -103,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/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/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/rng.cpp b/src/bindings/rng.cpp new file mode 100644 index 0000000..807e86b --- /dev/null +++ b/src/bindings/rng.cpp @@ -0,0 +1,32 @@ +// 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) { + // 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. + +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/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 1dc27aa..a63da08 100644 --- a/src/math/bernoulli.cpp +++ b/src/math/bernoulli.cpp @@ -1,6 +1,7 @@ -// Function declarations for Bernoulli distribution functions +// Function definitions 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..c3c6d60 100644 --- a/src/math/beta.cpp +++ b/src/math/beta.cpp @@ -1,7 +1,8 @@ -// Function declarations for beta distribution functions +// Function definitions for beta distribution functions #include #include #include +#include #include #include @@ -28,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) { @@ -41,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); } // ------------------------- @@ -63,26 +90,67 @@ namespace fastdist::math { // RNG // ------------------------- double beta_sample(const double alpha, const double beta) { - thread_local std::mt19937 rng{std::random_device{}()}; + // 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); - double b = gb(rng); + double a = ga(rng()); + double b = gb(rng()); return a / (a + b); } // ------------------------- // 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/src/math/binomial.cpp b/src/math/binomial.cpp index ac9c9e1..d2be9bb 100644 --- a/src/math/binomial.cpp +++ b/src/math/binomial.cpp @@ -1,7 +1,8 @@ -// Function declarations for binomial distribution functions +// Function definitions for binomial distribution functions #include #include #include +#include #include #include @@ -35,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); } @@ -89,10 +115,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/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 7425211..4fd9978 100644 --- a/src/math/discrete_uniform.cpp +++ b/src/math/discrete_uniform.cpp @@ -1,6 +1,7 @@ -// Function declarations for discrete uniform distribution functions +// Function definitions 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..8d99866 100644 --- a/src/math/exponential.cpp +++ b/src/math/exponential.cpp @@ -1,6 +1,8 @@ -// Function declarations for exponential distribution functions +// Function definitions for exponential distribution functions #include "fastdist/math/exponential.h" +#include #include +#include #include #include @@ -70,24 +72,43 @@ 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 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/gamma.cpp b/src/math/gamma.cpp index 87e0f83..d379075 100644 --- a/src/math/gamma.cpp +++ b/src/math/gamma.cpp @@ -1,7 +1,8 @@ -// Function declarations for gamma distribution functions +// Function definitions for gamma distribution functions #include #include #include +#include #include #include @@ -26,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 || @@ -84,9 +91,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()); } // ------------------------- @@ -115,7 +121,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/src/math/geometric.cpp b/src/math/geometric.cpp index 879b6ec..9435a06 100644 --- a/src/math/geometric.cpp +++ b/src/math/geometric.cpp @@ -1,6 +1,7 @@ -// Function declarations for geometric distribution functions +// Function definitions 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..4cccb47 100644 --- a/src/math/negative_binomial.cpp +++ b/src/math/negative_binomial.cpp @@ -1,7 +1,8 @@ -// Function declarations for negative binomial distribution functions +// Function definitions for negative binomial distribution functions #include #include #include +#include #include #include @@ -22,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); } // ------------------------- @@ -39,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); } @@ -108,9 +138,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..850bfbb 100644 --- a/src/math/normal.cpp +++ b/src/math/normal.cpp @@ -1,19 +1,46 @@ -// Function declarations for normal distribution functions +// Function definitions for normal distribution functions +#include #include #include +#include #include #include #include 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) { @@ -21,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) { @@ -31,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) { @@ -74,18 +98,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; } @@ -93,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 c9891a4..670b4d7 100644 --- a/src/math/poisson.cpp +++ b/src/math/poisson.cpp @@ -1,7 +1,8 @@ -// Function declarations for poisson distribution functions +// Function definitions for poisson distribution functions #include #include #include +#include #include #include @@ -11,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); @@ -39,13 +37,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); } @@ -86,21 +106,39 @@ 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 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/rng.cpp b/src/math/rng.cpp new file mode 100644 index 0000000..3233c25 --- /dev/null +++ b/src/math/rng.cpp @@ -0,0 +1,28 @@ +// 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::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{}()); } + +} // namespace fastdist::math diff --git a/src/math/uniform.cpp b/src/math/uniform.cpp index 28e2d53..e18af33 100644 --- a/src/math/uniform.cpp +++ b/src/math/uniform.cpp @@ -1,5 +1,7 @@ -// Function declarations for continuous uniform distribution functions +// Function definitions for continuous uniform distribution functions +#include #include +#include #include #include #include @@ -83,28 +85,62 @@ 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 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/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 9471ead..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,8 +106,8 @@ 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); + 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 8d76a0b..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,8 +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); - assert(std::abs(var - fastdist::math::geometric_variance(p)) < 0.15); + 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 9f22f5c..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,9 +107,9 @@ 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); + assert(std::abs(mean - fastdist::math::negative_binomial_mean(r, p)) < 0.075); - assert(std::abs(var - fastdist::math::negative_binomial_variance(r, p)) < 5e-1); + 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 ab65ffc..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,14 +123,16 @@ void test_normal() { const double mean = sum / N; const double var = sumsq / N - mean * mean; - assert(std::abs(mean - mu) < 5e-2); - assert(std::abs(var - sigma * sigma) < 5e-2); + 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; @@ -135,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 7d4b3c1..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,8 +103,8 @@ 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); + 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_beta.py b/tests/python/test_beta.py index f5e7e38..4b12932 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 @@ -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"): @@ -141,48 +138,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 +189,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 + ) # --------------------------------------------------------------------------- @@ -243,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_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_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_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)] 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) 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..2b25697 --- /dev/null +++ b/tests/python/test_rng.py @@ -0,0 +1,123 @@ +""" +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 pytest + +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 + + +@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__ + 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() 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: 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) 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." + )