diff --git a/CHANGELOG.md b/CHANGELOG.md index 11116d3..354cf85 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- **Accuracy characterization vs mpmath (#46)**: new + `tools/accuracy_sweep.cpp` (deterministic, bit-exact CSV of scalar and + FORCE_VECTORIZED batch `pdf`/`logpdf`/`cdf`/`quantile` over + support-aware grids, 19 distributions × 3 instances, tails to + `p = 1e-300`) and `tools/accuracy_vs_mpmath.py` (mp.dps = 50 oracle, + 42 unconditional self-checks, law-of-F tail normalization), generating + `docs/ACCURACY_CHARACTERIZATION.md`. Replaces the issue's pylibstats + approach, which pins to a released libstats and would characterize the + wrong code. The oracle carries its own large-parameter machinery where + mpmath's `betainc`/lower `gammainc` hang or raise (clean-room + continued-fraction incomplete beta, upper-gamma complement, far-tail + asymptotic guards, safeguarded log-log false-position quantile + inversion — each anchored by self-checks). Confirms the four pinned + accuracy gates and surfaces real findings: inconsistent batch NaN + propagation across 8 distributions, NaN/clamped returns at `+-inf` + inputs, quantile NaN/saturation at extreme `p`, and large-parameter + CDF accuracy limits (binomial `n = 1e6` off 1.3e-2 at the mean) — see + the doc's Findings section. + ### Changed - **LogNormal CDF lower tail no longer collapses (#49)**: every path diff --git a/PLAN.md b/PLAN.md index c7bb3d5..2407df2 100644 --- a/PLAN.md +++ b/PLAN.md @@ -152,7 +152,31 @@ history. covers the isStandardNormal_ path): fail-first max_rel 1.0 pre-fix; post-fix max 0.287 of the same law budget, batch-vs- scalar abs ≤ 1.11e-16. Correctness suite 53/53 on Zen 4 MSVC. - - #46 — Benchmark: SIMD accuracy characterization vs mpmath. Last. + - #46 — Accuracy characterization vs mpmath [DERIVED, in review]: + replaced the issue's pylibstats route (pins to released v2.2.0, + would characterize the wrong code) with tools/accuracy_sweep.cpp + (bit-exact deterministic CSV, 19 dists × 3 instances, scalar + + FORCE_VECTORIZED batch, tails to p=1e-300; two-process determinism + verified) + tools/accuracy_vs_mpmath.py (dps-50 oracle, 42 + self-checks) + docs/ACCURACY_CHARACTERIZATION.md. Oracle needed + substantial large-parameter hardening beyond the agents' build — + mpmath betainc/gammainc hang or raise for min(a,b) ≳ 5e3 even in + the central region: clean-room Lentz CF incomplete beta, upper- + gamma complement, far-tail lead guards (exp(lead) below every + double), safeguarded log-log false-position quantile solver + (secant plateaus and pure bisection both failed), asymptotic + normal seeds. CF/quadrature cross-validation at 1e-44. Gate + cross-check: gaussian law_frac 1.319 ↔ pinned 0.287·(1e-15 budget) + ≈ 1.29 bare-law equivalent; lognormal 1.784 ≲ 2.2 equivalent; + von Mises max_abs 1.7e-16/5.1e-16 vs gates 2.2e-16/8.9e-16. + Findings (doc Findings section + generated appendix): 86 contract + violations — batch NaN propagation inconsistent (8 dists), NaN at + ±inf where limits exist, batch logpdf(+inf) clamp −4605.0, + quantile NaN/saturation at extreme p, large-param CDF limits + (binomial n=1e6: 1.3e-2 at the mean → corvus #47/#52 remedy + class). [OPEN] follow-up issue candidates, on user's word: + (a) batch NaN propagation, (b) ±inf limit returns, (c) quantile + extreme-p contract. - **Parked on corvus adoption** (open, unmilestoned since 2026-08-20): #47 (bessel.h Tier 2 A&S fallback, ~10⁻⁷ on AppleClang) and #52 (Binomial CDF via PMF summation). Not "won't fix" — each tracks the diff --git a/docs/ACCURACY_CHARACTERIZATION.md b/docs/ACCURACY_CHARACTERIZATION.md new file mode 100644 index 0000000..8763958 --- /dev/null +++ b/docs/ACCURACY_CHARACTERIZATION.md @@ -0,0 +1,479 @@ +# Accuracy Characterization vs mpmath + +Issue #46. Full sweep of `pdf`/`logpdf`/`cdf`/`quantile` — scalar and batch +— across all 19 libstats distributions, measured against `mpmath` +(`mp.dps = 50`) reference implementations, generated by +`tools/accuracy_vs_mpmath.py`. + +> **PROVISIONAL.** Not the same standing as an audited per-tier accuracy +> claim. Currently one machine only: Zen 4 (AVX-512), MSVC Release; scalar +> libm baseline is MSVC's UCRT. Per-ISA results differ — the Kaby Lake +> (AVX2) and M1 (NEON) legs are pending and are not represented below. Do +> not read a bound in this document as holding on any other target until +> those legs land. + +## Regenerating + +```sh +cmake --preset release -G Ninja +cmake --build build --target accuracy_sweep +build/tools/accuracy_sweep > sweep.csv + tools/accuracy_vs_mpmath.py sweep.csv +``` + +`accuracy_sweep` (owned separately, `tools/accuracy_sweep.cpp` + +`tools/CMakeLists.txt`) emits the `dist,method,p1_bits,p2_bits,x_bits, +scalar_bits,batch_bits` CSV described in that tool's own header comment. +`accuracy_vs_mpmath.py` runs its oracle self-checks first (exits non-zero +without writing anything if any fails), then rewrites the generated block +below in place. Do not hand-edit anything between the `BEGIN GENERATED` / +`END GENERATED` markers — the next regeneration overwrites it. + +## Reading the tables + +- `max_abs` / `max_rel` / `p99_rel`: worst absolute error, worst relative + error, and 99th-percentile relative error across the swept rows for that + (distribution, method, source) group, against the mpmath reference. +- `law_frac(cdf)`: for `cdf` rows whose reference `F < 1e-3` only, the + worst observed relative error expressed as a **fraction of the + achievable-accuracy law** `rel(F) ~ |ln F| * 2^-52`, established in + `tests/test_lognormal_cdf_accuracy.cpp`'s `law_budget` comment (its + pinned gate budget adds headroom on top of this same bare law; this + column reports the bare law). A large `max_rel` in the deep tail + alongside a `law_frac` at or below roughly 1 is **the law, not a + defect** — a flat relative-error budget is mathematically unachievable + at double precision in that regime, for any implementation of these + closed forms. Treat a `law_frac` well above 1 as the signal worth + investigating, not the raw `max_rel` next to it. +- `batch_vs_scalar`: worst absolute/relative disagreement between the + batch (span-overload) and scalar paths on the same rows, shown on the + `batch` row. +- `worst_x`: the input at which `max_rel` was observed (for `quantile` + rows, this is the probability `p`, not a domain point — see the input + contract in `tools/accuracy_vs_mpmath.py`'s header comment). +- Any contract violation (NaN-in/NaN-out, `+-inf` boundary behavior) is + reported as its own row, separate from the error statistics above — it + is never folded into `max_abs`/`max_rel`/`p99_rel`. The full list is + the appendix at the end of the generated block. + +Four metric artifacts to read past — each is a property of the metric at +that row, not a defect: + +- **Von Mises `cdf` relative columns.** The Bessel-series CDF is + constructed and gated (#51) in **absolute** terms; its deep-tail + `max_rel`/`law_frac` values are astronomically large and meaningless + (an absolute-1e-16 error against `F ~ 1e-89` reads as `rel ~ 1e+73`). + Read `max_abs`, which must sit within the pinned gate budgets. +- **`logpdf` relative error near its zero crossing.** Where + `logpdf(x) ~ 0` (e.g. exponential at `lambda*x ~ ln(lambda)`), any + fixed absolute error reads as a large relative one. Read `max_abs`. +- **`max_abs` on singular-density rows.** Beta/Weibull with shape < 1 + have integrable singularities; near them the density reaches `1e+298` + and a `max_abs` of `1e+93` alongside `max_rel ~ 1e-10` is the value's + magnitude, not the error's. +- **Quantile rows at extreme `p`.** The oracle's reference quantiles run + to `1e-30000` and beyond, far outside double range; where the library + saturates (support bound, grid edge, `+-inf`), `max_rel ~ 1` records + the saturation point, not solver accuracy. The discrete-distribution + `rel = 1.0` rows at exact lattice probabilities (`F(k) == p`) are a + right-continuous-inverse convention difference, one integer off, not + an accuracy defect. + +## Findings + +What the sweep actually surfaced, beyond confirming the four pinned +gates (details row-by-row in the generated appendix): + +- **Batch NaN propagation is inconsistent** (35 rows, 8 distributions: + uniform, gamma, chi-squared, laplace, pareto, weibull, beta, + von Mises). Scalar paths propagate NaN; the batch paths return finite + values instead — uniform batch `pdf(NaN)` returns the in-support + density (the NaN comparison silently takes the in-range branch), + pareto batch `cdf(NaN)` returns `0.999...`. +- **`+-inf` inputs**: gamma/chi-squared/rayleigh/weibull scalar + `pdf`/`logpdf` return NaN at `+inf` where the limits (`0`, `-inf`) + exist; chi-squared/gamma **batch** `logpdf(+inf)` instead returns the + clamp constant `-4605.0` (`ln(1e-2000)`). Von Mises applies a + consistent saturation convention (`pdf -> 0`, `logpdf -> -inf`, + `cdf -> 0/1`) in both paths — recorded here as the convention. +- **Quantile extreme-`p` behavior**: gamma and chi-squared quantiles + return NaN for some deep-tail `p`; pareto returns `+inf`, student-t + `-inf`, where finite doubles exist; cauchy loses all accuracy at + `p = 1e-300` (`rel = 1.0`). +- **Large-parameter CDF accuracy**, oracle independently validated by + direct quadrature to `1e-44`: binomial `n = 1e6` is off by `1.3e-2` + relative **at the mean**; chi-squared `k = 1e5` by `8.3e-6`; beta + `(1e4, 1e4)` by `4.1e-9`. The corvus adoption line items (#47/#52, + incomplete gamma/beta cores) are the planned remedy for this class. + +## Consistency with the existing accuracy gates + +This document's numbers must not contradict the four ctest accuracy gates +already pinned in this repo; a mismatch here is a reason to re-check the +oracle or the gate, not to silently prefer one over the other: + +- `tests/test_trig_ulp_gates.cpp` (issue #95) — per-tier ULP budget for + `vector_cos`/`vector_sin` (max 1 ULP, mean 0.022-0.028 ULP). +- `tests/test_vonmises_cdf_accuracy.cpp` (issue #51) — Von Mises CDF, + flat budgets by kappa bucket (`2e-15` for kappa <= 100, `4e-15` for + kappa in (100, 1000]), against a 40-digit direct-quadrature oracle + independent of the library's Bessel-series implementation. +- `tests/test_lognormal_cdf_accuracy.cpp` (issue #49) — LogNormal CDF, + the `law_budget(F)` relative-error law this document's `law_frac(cdf)` + column also reports against. +- `tests/test_gaussian_cdf_accuracy.cpp` (issue #49, Gaussian instance) — + same law, own reference vectors, covers the `isStandardNormal_` fast + path at `(mean, sigma) = (0, 1)`. + + + +### beta + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 1.628e-9 | 4.113e-9 | 4.113e-9 | 4.468e+5 | - | 0.498166 | +| cdf | batch | 1.628e-9 | 4.113e-9 | 4.113e-9 | 4.468e+5 | abs=0, rel=0 | 0.498166 | +| logpdf | scalar | 4.811e-11 | 5.896e-10 | 5.896e-10 | - | - | 0.489076 | +| logpdf | batch | 2.993e-11 | 6.066e-10 | 6.066e-10 | - | abs=5.821e-11, rel=1.876e-11 | 0.510924 | +| logpdf | *(contract)* | 1 violation(s) -- see appendix | | | | | | +| pdf | scalar | 0.006387 | 3.051e-11 | 3.051e-11 | - | - | 0.497018 | +| pdf | batch | 0.006387 | 2.993e-11 | 2.993e-11 | - | abs=0.0004883, rel=1.819e-12 | 0.483195 | +| pdf | *(contract)* | 1 violation(s) -- see appendix | | | | | | +| quantile | scalar | 0.3712 | 8.018e+259 | 8.018e+259 | - | - | 0.001 | + +### binomial + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 0.0064 | 0.01279 | 0.01279 | 6.388e+5 | - | 300000 | +| cdf | batch | 0.0064 | 0.01279 | 0.01279 | 6.388e+5 | abs=0, rel=0 | 300000 | +| logpdf | scalar | 1.933e-9 | 1.716e-10 | 1.716e-10 | - | - | 300687 | +| logpdf | batch | 1.933e-9 | 1.716e-10 | 1.716e-10 | - | abs=0, rel=0 | 300687 | +| pdf | scalar | 4.006e-13 | 1.933e-9 | 1.933e-9 | - | - | 303666 | +| pdf | batch | 4.006e-13 | 1.933e-9 | 1.933e-9 | - | abs=0, rel=0 | 303666 | +| quantile | scalar | 14.0 | 4.667e-5 | 4.667e-5 | - | - | 0.5 | + +### cauchy + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 9.774e-17 | 2.908e-16 | 2.908e-16 | 0.06686 | - | 6.81795e+07 | +| cdf | batch | 9.774e-17 | 2.908e-16 | 2.908e-16 | 0.06686 | abs=0, rel=0 | 6.81795e+07 | +| logpdf | scalar | 7.235e-15 | 1.627e-15 | 1.627e-15 | - | - | -0.000318309 | +| logpdf | batch | 7.907e-15 | 6.426e-16 | 6.426e-16 | - | abs=1.421e-14, rel=1.552e-15 | -0.0031831 | +| pdf | scalar | 6.421e-11 | 6.839e-15 | 6.839e-15 | - | - | -3.15472e+20 | +| pdf | batch | 4.075e-11 | 7.566e-15 | 7.566e-15 | - | abs=5.821e-11, rel=1.424e-14 | -3.15472e+08 | +| quantile | scalar | 1.979e+57 | 1.0 | 1.0 | - | - | 1e-300 | + +### chi_squared + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 4.149e-6 | 8.298e-6 | 8.298e-6 | 1.614e+8 | - | 99999.3 | +| cdf | batch | 4.149e-6 | 8.298e-6 | 8.298e-6 | 1.614e+8 | abs=0, rel=0 | 99999.3 | +| cdf | *(contract)* | 6 violation(s) -- see appendix | | | | | | +| logpdf | scalar | 9.381e-11 | 1.197e-11 | 1.197e-11 | - | - | 99427.3 | +| logpdf | batch | 1.212e-10 | 1.541e-11 | 1.541e-11 | - | abs=1.164e-10, rel=1.627e-11 | 99622.7 | +| logpdf | *(contract)* | 7 violation(s) -- see appendix | | | | | | +| pdf | scalar | 1.923e+93 | 9.382e-11 | 9.382e-11 | - | - | 99427.3 | +| pdf | batch | 2.086e+93 | 1.212e-10 | 1.212e-10 | - | abs=4.009e+93, rel=1.164e-10 | 101672 | +| pdf | *(contract)* | 2 violation(s) -- see appendix | | | | | | +| quantile | scalar | 3540.0 | 4.136e+189 | 4.136e+189 | - | - | 1e-300 | +| quantile | *(contract)* | 1 violation(s) -- see appendix | | | | | | + +### discrete + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 5.306e-17 | 8.234e-17 | 8.234e-17 | 0.02553 | - | 288675 | +| cdf | batch | 9.508e-17 | 1.745e-16 | 1.745e-16 | 0.02631 | abs=1.11e-16, rel=2.22e-16 | 0 | +| logpdf | scalar | 7.829e-16 | 9.427e-17 | 9.427e-17 | - | - | 0 | +| logpdf | batch | 7.829e-16 | 9.427e-17 | 9.427e-17 | - | abs=0, rel=0 | 0 | +| pdf | scalar | 5.551e-18 | 7.833e-17 | 7.833e-17 | - | - | -1e+06 | +| pdf | batch | 5.551e-18 | 7.833e-17 | 7.833e-17 | - | abs=0, rel=0 | -1e+06 | +| quantile | scalar | 1.0 | 1.0 | 1.0 | - | - | 0.1 | + +### exponential + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 8.835e-17 | 1.741e-16 | 1.741e-16 | 0.03675 | - | 1.0005e-09 | +| cdf | batch | 6.79e-17 | 1.741e-16 | 1.741e-16 | 0.03675 | abs=1.11e-16, rel=1.708e-16 | 1.0005e-09 | +| logpdf | scalar | 3.941e-15 | 2.966e-5 | 2.966e-5 | - | - | 1.38155e-05 | +| logpdf | batch | 3.941e-15 | 2.966e-5 | 2.966e-5 | - | abs=0, rel=0 | 1.38155e-05 | +| pdf | scalar | 5.492e-11 | 3.438e-15 | 3.438e-15 | - | - | 3.45396e-05 | +| pdf | batch | 5.492e-11 | 3.438e-15 | 3.438e-15 | - | abs=5.821e-11, rel=1.735e-16 | 3.45396e-05 | +| quantile | scalar | 4.46e-12 | 0.0007993 | 0.0007993 | - | - | 1e-15 | + +### gamma + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 8.718e-8 | 1.768e-7 | 1.768e-7 | 6.558e+7 | - | 9.94789e+06 | +| cdf | batch | 8.718e-8 | 1.768e-7 | 1.768e-7 | 6.558e+7 | abs=0, rel=0 | 9.94789e+06 | +| cdf | *(contract)* | 6 violation(s) -- see appendix | | | | | | +| logpdf | scalar | 2.694e-11 | 2.144e-12 | 2.144e-12 | - | - | 9.94789e+06 | +| logpdf | batch | 3.42e-11 | 2.674e-12 | 2.674e-12 | - | abs=5.821e-11, rel=4.55e-12 | 1.00842e+07 | +| logpdf | *(contract)* | 7 violation(s) -- see appendix | | | | | | +| pdf | scalar | 1.54e+93 | 2.694e-11 | 2.694e-11 | - | - | 9.94789e+06 | +| pdf | batch | 6.271e+92 | 3.42e-11 | 3.42e-11 | - | abs=2.167e+93, rel=5.821e-11 | 1.00842e+07 | +| pdf | *(contract)* | 2 violation(s) -- see appendix | | | | | | +| quantile | scalar | 8.256e+5 | 7.071e+139 | 7.071e+139 | - | - | 1e-300 | +| quantile | *(contract)* | 1 violation(s) -- see appendix | | | | | | + +### gaussian + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 1.417e-16 | 1.877e-14 | 1.877e-14 | 1.319 | - | 889773 | +| cdf | batch | 6.733e-17 | 1.877e-14 | 1.877e-14 | 1.319 | abs=1.11e-16, rel=1.388e-16 | 889773 | +| logpdf | scalar | 1.696e-14 | 1.755e-15 | 1.755e-15 | - | - | 0.00471902 | +| logpdf | batch | 9.85e-15 | 9.733e-16 | 9.733e-16 | - | abs=7.105e-15, rel=9.584e-16 | -0.00271902 | +| pdf | scalar | 6.482e-14 | 1.555e-14 | 1.555e-14 | - | - | -0.0100227 | +| pdf | batch | 6.482e-14 | 8.43e-15 | 8.43e-15 | - | abs=2.842e-14, rel=7.123e-15 | -0.0100227 | +| quantile | scalar | 3.081e+4 | 0.4439 | 0.4439 | - | - | 1e-15 | + +### geometric + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 1.0 | 1.0 | 1.0 | 0.1689 | - | 4.001e+09 | +| cdf | batch | 1.0 | 1.0 | 1.0 | 0.1689 | abs=0, rel=0 | 4.001e+09 | +| logpdf | scalar | 2.878e-8 | 2.836e-11 | 2.836e-11 | - | - | 1.001e+09 | +| logpdf | batch | 2.878e-8 | 2.836e-11 | 2.836e-11 | - | abs=0, rel=0 | 1.001e+09 | +| logpdf | *(contract)* | 2 violation(s) -- see appendix | | | | | | +| pdf | scalar | 3.519e-17 | 1.872e-8 | 1.872e-8 | - | - | 6.51e+08 | +| pdf | batch | 3.519e-17 | 1.872e-8 | 1.872e-8 | - | abs=0, rel=0 | 6.51e+08 | +| quantile | scalar | 2.354e+7 | 0.6815 | 0.6815 | - | - | 1 | + +### laplace + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 6.021e-17 | 4.306e-14 | 2.716e-14 | 0.442 | - | -4.51927e+08 | +| cdf | batch | 5.79e-17 | 4.306e-14 | 3.504e-14 | 0.6266 | abs=1.11e-16, rel=5.686e-14 | -4.51927e+08 | +| cdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| logpdf | scalar | 7.249e-14 | 0.9175 | 1.17e-5 | - | - | -1.31224e-05 | +| logpdf | batch | 7.249e-14 | 0.9175 | 1.17e-5 | - | abs=0, rel=0 | -1.31224e-05 | +| logpdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| pdf | scalar | 5.159e-10 | 5.504e-14 | 5.501e-14 | - | - | -551.927 | +| pdf | batch | 5.159e-10 | 5.501e-14 | 5.499e-14 | - | abs=5.821e-11, rel=2.061e-16 | -597.979 | +| pdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| quantile | scalar | 8.498e-9 | 1.601e-16 | 1.601e-16 | - | - | 1e-15 | + +### lognormal + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 1.108e-16 | 1.88e-14 | 1.88e-14 | 1.784 | - | 0.895631 | +| cdf | batch | 1.108e-16 | 1.88e-14 | 1.88e-14 | 1.784 | abs=1.11e-16, rel=2.082e-15 | 0.895631 | +| logpdf | scalar | 1.2e-14 | 9.822e-16 | 9.822e-16 | - | - | 0.96957 | +| logpdf | batch | 1.2e-14 | 1.822e-15 | 1.822e-15 | - | abs=7.105e-15, rel=8.85e-16 | 0.96957 | +| pdf | scalar | 1.681e-14 | 1.194e-14 | 1.194e-14 | - | - | 3.75136e+09 | +| pdf | batch | 1.943e-14 | 1.194e-14 | 1.194e-14 | - | abs=1.421e-14, rel=7.145e-15 | 3.75136e+09 | +| quantile | scalar | 3.743e+9 | 473.6 | 473.6 | - | - | 1 | + +### negative_binomial + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 2.143e-9 | 4.477e-9 | 4.477e-9 | 5.311e+5 | - | 970100 | +| cdf | batch | 2.143e-9 | 4.477e-9 | 4.477e-9 | 5.311e+5 | abs=0, rel=0 | 970100 | +| logpdf | scalar | 8.992e-8 | 2.294e-10 | 2.294e-10 | - | - | 1.00492e+06 | +| logpdf | batch | 8.992e-8 | 2.294e-10 | 2.294e-10 | - | abs=0, rel=0 | 1.00492e+06 | +| pdf | scalar | 5.182e-14 | 3.423e-9 | 3.423e-9 | - | - | 890501 | +| pdf | batch | 5.182e-14 | 3.423e-9 | 3.423e-9 | - | abs=0, rel=0 | 890501 | +| quantile | scalar | 70.0 | 0.25 | 0.25 | - | - | 0.5 | + +### pareto + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 5.418e-15 | 0.0002983 | 0.0002983 | 3.89e+10 | - | 1e-06 | +| cdf | batch | 6.563e-14 | 0.0002983 | 0.0002983 | 3.89e+10 | abs=6.528e-14, rel=6.051e-8 | 1e-06 | +| cdf | *(contract)* | 1 violation(s) -- see appendix | | | | | | +| logpdf | scalar | 2.325e-13 | 2.351e-14 | 2.351e-14 | - | - | 1.00223e+06 | +| logpdf | batch | 1.563e-13 | 1.566e-14 | 1.566e-14 | - | abs=2.274e-13, rel=2.469e-14 | 1e+06 | +| logpdf | *(contract)* | 1 violation(s) -- see appendix | | | | | | +| pdf | scalar | 1.697e-11 | 2.325e-13 | 2.325e-13 | - | - | 1.00921e+06 | +| pdf | batch | 9.756e-12 | 1.562e-13 | 1.562e-13 | - | abs=1.819e-11, rel=2.275e-13 | 1.09648e+06 | +| pdf | *(contract)* | 1 violation(s) -- see appendix | | | | | | +| quantile | scalar | 9.731e+179 | 9.731e-15 | 9.731e-15 | - | - | 0.99 | +| quantile | *(contract)* | 3 violation(s) -- see appendix | | | | | | + +### poisson + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 0.0007912 | 0.00158 | 0.00158 | 4.912e+5 | - | 100000 | +| cdf | batch | 0.0007912 | 0.00158 | 0.00158 | 4.912e+5 | abs=0, rel=0 | 100000 | +| logpdf | scalar | 3.281e-9 | 3.328e-11 | 3.328e-11 | - | - | 100000 | +| logpdf | batch | 3.281e-9 | 3.328e-11 | 3.328e-11 | - | abs=0, rel=0 | 100000 | +| pdf | scalar | 7.285e-7 | 0.004219 | 0.004219 | - | - | 99368 | +| pdf | batch | 2.802e-13 | 2.221e-10 | 2.221e-10 | - | abs=7.285e-7, rel=0.004213 | 100000 | +| quantile | scalar | 3.0 | 1.0 | 1.0 | - | - | 1e-300 | + +### rayleigh + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 9.255e-17 | 3.223e-16 | 3.223e-16 | 0.1166 | - | 1414.21 | +| cdf | batch | 1.032e-16 | 3.223e-16 | 3.223e-16 | 0.1166 | abs=1.11e-16, rel=3.676e-16 | 1414.21 | +| logpdf | scalar | 4.621e-15 | 5.933e-16 | 5.933e-16 | - | - | 4.47035e-14 | +| logpdf | batch | 4.621e-15 | 8.357e-16 | 8.357e-16 | - | abs=4.441e-15, rel=1.429e-15 | 4.47035e-14 | +| logpdf | *(contract)* | 6 violation(s) -- see appendix | | | | | | +| pdf | scalar | 1.272e-9 | 4.603e-15 | 4.603e-15 | - | - | 6.78614e+06 | +| pdf | batch | 7.065e-10 | 4.603e-15 | 4.603e-15 | - | abs=1.048e-9, rel=4.501e-15 | 6.78614e+06 | +| pdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| quantile | scalar | 1.788e-5 | 0.0003997 | 0.0003997 | - | - | 1e-15 | + +### student_t + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 1.486e-7 | 0.003744 | 0.003744 | 2.631e+11 | - | -11.0227 | +| cdf | batch | 1.486e-7 | 0.003744 | 0.003744 | 2.631e+11 | abs=0, rel=0 | -11.0227 | +| logpdf | scalar | 2.073e-10 | 2.255e-10 | 2.255e-10 | - | - | 0 | +| logpdf | batch | 2.516e-10 | 2.623e-10 | 2.623e-10 | - | abs=5.543e-11, rel=4.214e-11 | 0.125661 | +| pdf | scalar | 8.268e-11 | 2.073e-10 | 2.073e-10 | - | - | -11.0227 | +| pdf | batch | 9.622e-11 | 2.516e-10 | 2.516e-10 | - | abs=1.548e-11, rel=5.543e-11 | -0.518671 | +| quantile | scalar | 3.082e+14 | 1.0 | 1.0 | - | - | 1 | +| quantile | *(contract)* | 3 violation(s) -- see appendix | | | | | | + +### uniform + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 6.985e-17 | 1.049e-16 | 1.001e-16 | 0.02343 | - | 0.0003 | +| cdf | batch | 6.985e-17 | 1.049e-16 | 1.001e-16 | 0.02343 | abs=0, rel=0 | 0.0003 | +| logpdf | scalar | 3.04e-16 | 4.891e-17 | 4.891e-17 | - | - | -0.001 | +| logpdf | batch | 3.04e-16 | 4.891e-17 | 4.891e-17 | - | abs=0, rel=0 | -0.001 | +| logpdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| pdf | scalar | 1.041e-14 | 2.092e-17 | 2.092e-17 | - | - | -1e+08 | +| pdf | batch | 1.041e-14 | 2.092e-17 | 2.092e-17 | - | abs=0, rel=0 | -1e+08 | +| pdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| quantile | scalar | 9.074e-9 | 1.908e-16 | 1.908e-16 | - | - | 0.999999 | + +### von_mises + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 1.702e-16 | 5.605e+23 | 5.605e+23 | 2.74e+37 | - | -1.45952 | +| cdf | batch | 5.101e-16 | 5.561e+63 | 5.561e+63 | 1.36e+77 | abs=5.551e-16, rel=1.0 | -2.52455 | +| logpdf | scalar | 1.426e-14 | 5.191e-15 | 5.191e-15 | - | - | -0.128412 | +| logpdf | batch | 1.426e-14 | 5.191e-15 | 5.191e-15 | - | abs=1.421e-14, rel=2.56e-15 | -0.128412 | +| logpdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| pdf | scalar | 2.148e-14 | 1.422e-14 | 1.422e-14 | - | - | -0.374577 | +| pdf | batch | 2.148e-14 | 1.422e-14 | 1.422e-14 | - | abs=5.551e-17, rel=1.43e-14 | -0.374577 | +| pdf | *(contract)* | 3 violation(s) -- see appendix | | | | | | +| quantile | scalar | 6.283 | 1.407e+15 | 1.407e+15 | - | - | 0.5 | + +### weibull + +| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x | +|---|---|---|---|---|---|---|---| +| cdf | scalar | 6.31e-14 | 1.147e-13 | 1.147e-13 | 71.4 | - | 9977.52 | +| cdf | batch | 6.31e-14 | 1.532e-13 | 1.532e-13 | 71.4 | abs=5.44e-14, rel=1.474e-13 | 9898.23 | +| cdf | *(contract)* | 1 violation(s) -- see appendix | | | | | | +| logpdf | scalar | 2.319e-12 | 1.03e-13 | 1.03e-13 | - | - | 10195.1 | +| logpdf | batch | 2.333e-12 | 1.056e-13 | 1.056e-13 | - | abs=8.527e-14, rel=1.453e-14 | 10195.1 | +| logpdf | *(contract)* | 6 violation(s) -- see appendix | | | | | | +| pdf | scalar | 1.563e+284 | 2.319e-12 | 2.319e-12 | - | - | 10360.6 | +| pdf | batch | 1.563e+284 | 2.333e-12 | 2.333e-12 | - | abs=9.12e+51, rel=8.524e-14 | 10360.6 | +| pdf | *(contract)* | 4 violation(s) -- see appendix | | | | | | +| quantile | scalar | 1.114e+136 | 0.0005329 | 0.0005329 | - | - | 1e-15 | + +### Contract findings (appendix) + +86 contract violations across the sweep. `csv_line` indexes the sweep CSV this report was generated from (see the commit/isa banner in the regeneration log). + +| dist | method | source | csv_line | finding | +|---|---|---|---|---| +| beta | logpdf | batch | 2831 | NaN input did not produce NaN batch output | +| beta | pdf | batch | 2814 | NaN input did not produce NaN batch output | +| chi_squared | cdf | batch | 3085 | NaN input did not produce NaN batch output | +| chi_squared | cdf | batch | 3086 | reference is finite (1.0), batch_bits decoded to nan | +| chi_squared | cdf | batch | 3190 | NaN input did not produce NaN batch output | +| chi_squared | cdf | batch | 3191 | reference is finite (1.0), batch_bits decoded to nan | +| chi_squared | cdf | batch | 3295 | NaN input did not produce NaN batch output | +| chi_squared | cdf | batch | 3296 | reference is finite (1.0), batch_bits decoded to nan | +| chi_squared | logpdf | scalar | 3055 | reference is -inf, scalar_bits decoded to nan | +| chi_squared | logpdf | batch | 3055 | reference is -inf, batch_bits decoded to nan | +| chi_squared | logpdf | batch | 3056 | reference is -inf, batch_bits decoded to -4605.0 | +| chi_squared | logpdf | batch | 3162 | reference is -inf, batch_bits decoded to -4605.0 | +| chi_squared | logpdf | scalar | 3266 | reference is -inf, scalar_bits decoded to nan | +| chi_squared | logpdf | batch | 3266 | reference is -inf, batch_bits decoded to nan | +| chi_squared | logpdf | batch | 3267 | reference is -inf, batch_bits decoded to -4605.0 | +| chi_squared | pdf | scalar | 3024 | reference is finite (0.0), scalar_bits decoded to nan | +| chi_squared | pdf | scalar | 3236 | reference is finite (0.0), scalar_bits decoded to nan | +| chi_squared | quantile | scalar | 3298 | reference is finite (84333.5), scalar_bits decoded to nan | +| gamma | cdf | batch | 1465 | NaN input did not produce NaN batch output | +| gamma | cdf | batch | 1466 | reference is finite (1.0), batch_bits decoded to nan | +| gamma | cdf | batch | 1570 | NaN input did not produce NaN batch output | +| gamma | cdf | batch | 1571 | reference is finite (1.0), batch_bits decoded to nan | +| gamma | cdf | batch | 1675 | NaN input did not produce NaN batch output | +| gamma | cdf | batch | 1676 | reference is finite (1.0), batch_bits decoded to nan | +| gamma | logpdf | scalar | 1434 | reference is -inf, scalar_bits decoded to nan | +| gamma | logpdf | batch | 1434 | reference is -inf, batch_bits decoded to nan | +| gamma | logpdf | batch | 1435 | reference is -inf, batch_bits decoded to -4605.0 | +| gamma | logpdf | batch | 1542 | reference is -inf, batch_bits decoded to -4605.0 | +| gamma | logpdf | scalar | 1646 | reference is -inf, scalar_bits decoded to nan | +| gamma | logpdf | batch | 1646 | reference is -inf, batch_bits decoded to nan | +| gamma | logpdf | batch | 1647 | reference is -inf, batch_bits decoded to -4605.0 | +| gamma | pdf | scalar | 1402 | reference is finite (0.0), scalar_bits decoded to nan | +| gamma | pdf | scalar | 1616 | reference is finite (0.0), scalar_bits decoded to nan | +| gamma | quantile | scalar | 1678 | reference is finite (6.73769e+6), scalar_bits decoded to nan | +| geometric | logpdf | scalar | 5887 | reference is finite (-4014.82), scalar_bits decoded to -inf | +| geometric | logpdf | scalar | 5888 | reference is finite (-16014.8), scalar_bits decoded to -inf | +| laplace | cdf | batch | 3445 | NaN input did not produce NaN batch output | +| laplace | cdf | batch | 3595 | NaN input did not produce NaN batch output | +| laplace | cdf | batch | 3745 | NaN input did not produce NaN batch output | +| laplace | logpdf | batch | 3400 | NaN input did not produce NaN batch output | +| laplace | logpdf | batch | 3550 | NaN input did not produce NaN batch output | +| laplace | logpdf | batch | 3700 | NaN input did not produce NaN batch output | +| laplace | pdf | batch | 3355 | NaN input did not produce NaN batch output | +| laplace | pdf | batch | 3505 | NaN input did not produce NaN batch output | +| laplace | pdf | batch | 3655 | NaN input did not produce NaN batch output | +| pareto | cdf | batch | 3946 | NaN input did not produce NaN batch output | +| pareto | logpdf | batch | 3920 | NaN input did not produce NaN batch output | +| pareto | pdf | batch | 3894 | NaN input did not produce NaN batch output | +| pareto | quantile | scalar | 3961 | reference is finite (1.0e+594), scalar_bits decoded to inf | +| pareto | quantile | scalar | 3962 | reference is finite (9.99992e+993), scalar_bits decoded to inf | +| pareto | quantile | scalar | 3963 | reference is finite (1.08324e+1494), scalar_bits decoded to inf | +| rayleigh | logpdf | scalar | 4129 | reference is -inf, scalar_bits decoded to nan | +| rayleigh | logpdf | batch | 4129 | reference is -inf, batch_bits decoded to nan | +| rayleigh | logpdf | scalar | 4237 | reference is -inf, scalar_bits decoded to nan | +| rayleigh | logpdf | batch | 4237 | reference is -inf, batch_bits decoded to nan | +| rayleigh | logpdf | scalar | 4345 | reference is -inf, scalar_bits decoded to nan | +| rayleigh | logpdf | batch | 4345 | reference is -inf, batch_bits decoded to nan | +| rayleigh | pdf | scalar | 4098 | reference is finite (0.0), scalar_bits decoded to nan | +| rayleigh | pdf | scalar | 4206 | reference is finite (0.0), scalar_bits decoded to nan | +| rayleigh | pdf | scalar | 4314 | reference is finite (0.0), scalar_bits decoded to nan | +| student_t | quantile | scalar | 1783 | reference is finite (-1.56839e+60), scalar_bits decoded to -inf | +| student_t | quantile | scalar | 1888 | reference is finite (-1.59857e+299), scalar_bits decoded to -inf | +| student_t | quantile | scalar | 1993 | reference is finite (-37.0598), scalar_bits decoded to -inf | +| uniform | logpdf | batch | 1060 | NaN input did not produce NaN batch output | +| uniform | logpdf | batch | 1199 | NaN input did not produce NaN batch output | +| uniform | logpdf | batch | 1320 | NaN input did not produce NaN batch output | +| uniform | pdf | batch | 1012 | NaN input did not produce NaN batch output | +| uniform | pdf | batch | 1161 | NaN input did not produce NaN batch output | +| uniform | pdf | batch | 1286 | NaN input did not produce NaN batch output | +| von_mises | logpdf | batch | 2391 | NaN input did not produce NaN batch output | +| von_mises | logpdf | batch | 2499 | NaN input did not produce NaN batch output | +| von_mises | logpdf | batch | 2617 | NaN input did not produce NaN batch output | +| von_mises | pdf | batch | 2360 | NaN input did not produce NaN batch output | +| von_mises | pdf | batch | 2468 | NaN input did not produce NaN batch output | +| von_mises | pdf | batch | 2581 | NaN input did not produce NaN batch output | +| weibull | cdf | batch | 4576 | NaN input did not produce NaN batch output | +| weibull | logpdf | scalar | 4453 | reference is -inf, scalar_bits decoded to nan | +| weibull | logpdf | batch | 4453 | reference is -inf, batch_bits decoded to nan | +| weibull | logpdf | batch | 4527 | reference is inf, batch_bits decoded to -inf | +| weibull | logpdf | batch | 4550 | NaN input did not produce NaN batch output | +| weibull | logpdf | scalar | 4654 | reference is -inf, scalar_bits decoded to nan | +| weibull | logpdf | batch | 4654 | reference is -inf, batch_bits decoded to nan | +| weibull | pdf | scalar | 4422 | reference is finite (0.0), scalar_bits decoded to nan | +| weibull | pdf | batch | 4501 | reference is inf, batch_bits decoded to 0.0 | +| weibull | pdf | batch | 4524 | NaN input did not produce NaN batch output | +| weibull | pdf | scalar | 4623 | reference is finite (0.0), scalar_bits decoded to nan | + + diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 22b9d87..764bf9a 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -48,6 +48,7 @@ add_libstats_tool(parameter_recovery_benchmark parameter_recovery_benchmark.cpp) add_libstats_tool(threshold_validator threshold_validator.cpp) add_libstats_tool(toctou_validator toctou_validator.cpp) add_libstats_tool(copy_move_stress copy_move_stress.cpp) +add_libstats_tool(accuracy_sweep accuracy_sweep.cpp) # Issue #33 opt-in dev tools -- gated separately since these are research/experiment tools, not part # of the standard tool set built by every developer. @@ -127,3 +128,7 @@ message( STATUS " - copy_move_stress: Concurrent copy/move throughput stress test across all 19 distributions (replaces test_copy_move_stress)" ) +message( + STATUS + " - accuracy_sweep: Issue #46 deterministic bit-exact scalar/batch pdf/logpdf/cdf/quantile sweep across all 19 distributions for the mpmath oracle comparison" +) diff --git a/tools/accuracy_sweep.cpp b/tools/accuracy_sweep.cpp new file mode 100644 index 0000000..731a760 --- /dev/null +++ b/tools/accuracy_sweep.cpp @@ -0,0 +1,642 @@ +/** + * @file accuracy_sweep.cpp + * @brief Issue #46: deterministic bit-exact accuracy sweep across all 19 distributions. + * + * Evaluates the scalar AND batch (span/SIMD) code paths of every distribution's + * pdf/logpdf/cdf/quantile methods over fixed, support-aware characterization grids + * and dumps the results as bit-exact hex to a CSV consumed by a sibling mpmath + * oracle (tools/accuracy_vs_mpmath.py, owned by a different agent). This tool + * replaces the issue's original pylibstats approach -- pylibstats pins to a + * *released* libstats and would characterize the wrong code; this binary always + * characterizes the code actually checked out. + * + * Determinism is the load-bearing property: no randomness anywhere. The x/k + * grids are built from each distribution's OWN getQuantile()/getMean()/ + * getVariance(), which makes grid construction support-aware for free (a + * quantile near p=0 or p=1 lands deep in whatever tail the support allows) + * without hand-deriving per-family tail formulas. + * + * Output contract (do not change without updating the sibling tool): + * # libstats accuracy_sweep v1 + * # commit= isa= date=YYYY-MM-DD + * dist,method,p1_bits,p2_bits,x_bits,scalar_bits,batch_bits + * ... + * # skipped_quantile (only emitted when n > 0) + * + * Usage: accuracy_sweep + */ + +// Use tool_utils.h for the consolidated (LIBSTATS_FULL_INTERFACE) libstats.h include, +// which pulls in all 19 distribution headers. +#include "tool_utils.h" + +#include "libstats/platform/simd_policy.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +//============================================================================== +// Bit-pattern formatting +//============================================================================== + +std::uint64_t bitsOf(double v) { + std::uint64_t u = 0; + std::memcpy(&u, &v, sizeof(u)); + return u; +} + +std::string hex64(std::uint64_t bits) { + std::ostringstream oss; + oss << "0x" << std::hex << std::setw(16) << std::setfill('0') << bits; + return oss.str(); +} + +std::string fmtBits(double v) { + return hex64(bitsOf(v)); +} + +//============================================================================== +// Row / CSV plumbing +//============================================================================== + +struct Row { + std::string dist; + std::string method; + std::string p1; + std::string p2; + std::string x; + std::string scalar; + std::string batch; + + std::string toCsvLine() const { + std::string line; + line.reserve(dist.size() + method.size() + p1.size() + p2.size() + x.size() + + scalar.size() + batch.size() + 8); + line += dist; + line += ','; + line += method; + line += ','; + line += p1; + line += ','; + line += p2; + line += ','; + line += x; + line += ','; + line += scalar; + line += ','; + line += batch; + return line; + } +}; + +Row makeRow(const std::string& dist, const std::string& method, double p1, double p2, bool hasP2, + double x, double scalar, double batch, bool hasBatch) { + Row r; + r.dist = dist; + r.method = method; + r.p1 = fmtBits(p1); + r.p2 = hasP2 ? fmtBits(p2) : std::string("0x0"); + r.x = fmtBits(x); + r.scalar = fmtBits(scalar); + r.batch = hasBatch ? fmtBits(batch) : std::string("-"); + return r; +} + +// Fixed, dependency-ordered dist name list (matches the issue #46 contract and +// drives both the summary table and the skip-count trailer, so both keep a +// deterministic dist order regardless of sweep call order). +constexpr std::array kDistNames = { + "gaussian", "lognormal", "exponential", "uniform", "poisson", + "gamma", "discrete", "student_t", "cauchy", "von_mises", + "binomial", "negative_binomial", "geometric", "beta", "chi_squared", + "laplace", "pareto", "rayleigh", "weibull"}; + +class Sink { + public: + Sink() { + for (const char* name : kDistNames) { + counts_.emplace_back(name, 0); + skips_.emplace_back(name, 0); + } + } + + void addRow(Row row) { + bump(counts_, row.dist); + rows_.push_back(std::move(row)); + } + + void addSkip(const std::string& dist) { + bump(skips_, dist); + } + + const std::vector& rows() const { + return rows_; + } + const std::vector>& counts() const { + return counts_; + } + const std::vector>& skips() const { + return skips_; + } + + private: + static void bump(std::vector>& v, const std::string& dist) { + for (auto& [name, n] : v) { + if (name == dist) { + ++n; + return; + } + } + v.emplace_back(dist, 1); + } + + std::vector rows_; + std::vector> counts_; + std::vector> skips_; +}; + +//============================================================================== +// Deterministic grids (no randomness anywhere) +//============================================================================== + +// Probability grid used to *construct* the continuous x-grid via getQuantile(p). +// Log-spaced deep into both tails (down to p=1e-300, i.e. F~1e-300 in the left +// tail / 1-F~1e-300 in the right tail where the support allows it reaching that +// far) plus a denser cluster around the center. Symmetric about p=0.5 by +// construction. ~57 points. +const std::vector& xConstructionPGrid() { + static const std::vector grid = [] { + std::vector lower; // all values strictly < 0.5 + for (int e : {300, 280, 260, 240, 220, 200, 180, 160, 140, 120, 100, 80, 60, 40, 20}) { + lower.push_back(std::pow(10.0, -e)); + } + for (double e : {15.0, 10.0, 8.0, 6.0, 4.0, 3.0, 2.0, 1.0, 0.7, 0.52}) { + lower.push_back(std::pow(10.0, -e)); + } + for (double p : {0.35, 0.4, 0.45}) { + lower.push_back(p); + } + // Strictly ascending by construction: both exponent lists descend + // (so the values ascend), and 10^-0.52 ~ 0.302 < 0.35. Assert + // rather than sort+unique -- GCC 13's -Werror=strict-overflow + // fires a pointer-wraparound false positive when it constexpr- + // expands std::sort inside this immediately-invoked initializer. + for (std::size_t i = 1; i < lower.size(); ++i) { + assert(lower[i - 1] < lower[i]); + } + + std::vector full; + full.reserve(lower.size() * 2 + 1); + for (double p : lower) { + full.push_back(p); + } + full.push_back(0.5); + for (auto it = lower.rbegin(); it != lower.rend(); ++it) { + full.push_back(1.0 - *it); + } + return full; + }(); + return grid; +} + +// Probability grid for method=quantile rows (the contract's fixed 15-point list). +const std::vector& quantilePGrid() { + static const std::vector grid = {1e-300, 1e-15, 1e-10, 1e-6, 1e-3, 0.01, + 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, + 1.0 - 1e-6, 1.0 - 1e-10, 1.0 - 1e-15}; + return grid; +} + +// Boundary epsilon offsets applied to both-sides-bounded supports (uniform, beta), +// per the issue's grid spec. +const std::array& boundaryEpsilons() { + static const std::array eps = {1e-15, 1e-9, 1e-4}; + return eps; +} + +// k-landmark multipliers (in units of standard deviation, both directions around +// the mean) used to build the discrete integer grid. Reaches far enough into the +// tail (mean +/- 16000*sd) that cdf reaches 1-1e-12 for every discrete instance +// used below; landmarks are clamped to the distribution's own support bounds. +const std::array& discreteTMultipliers() { + static const std::array t = {0.5, 1, 1.5, 2, 3, 4, 5, 6, + 8, 10, 15, 20, 30, 50, 75, 100, + 150, 250, 400, 650, 1000, 4000, 16000}; + return t; +} + +//============================================================================== +// Generic grid builders (templated on the concrete distribution type; every +// distribution exposes the same DistributionInterface surface used here) +//============================================================================== + +// Builds the continuous x-grid for one instance: quantile-derived points spanning +// both tails as deep as the support allows, plus explicit near-boundary points +// for both-sides-bounded supports (uniform, beta). Sorted, deduplicated. +template +std::vector buildContinuousXGrid(Dist& dist, bool boundedBothSides) { + std::vector xs; + for (double p : xConstructionPGrid()) { + try { + double x = dist.getQuantile(p); + if (std::isfinite(x)) { + xs.push_back(x); + } + } catch (...) { + // Support doesn't reach this deep at this p for this instance; skip. + } + } + if (boundedBothSides) { + double lo = dist.getSupportLowerBound(); + double hi = dist.getSupportUpperBound(); + for (double eps : boundaryEpsilons()) { + xs.push_back(lo + eps); + xs.push_back(hi - eps); + } + } + std::sort(xs.begin(), xs.end()); + xs.erase(std::unique(xs.begin(), xs.end()), xs.end()); + return xs; +} + +// Builds the discrete k-grid for one instance from the distribution's own +// mean/variance/support -- landmark points (0, 1, 2, mean, mean +/- t*sd for a +// deterministic multiplier ladder, and the support edges), not sampling. +template +std::vector buildDiscreteKGrid(Dist& dist) { + double lo = dist.getSupportLowerBound(); + double hi = dist.getSupportUpperBound(); + double mean = dist.getMean(); + double var = dist.getVariance(); + double sd = (std::isfinite(var) && var > 0.0) ? std::sqrt(var) : 1.0; + if (!std::isfinite(mean)) { + mean = 0.0; + } + + auto clamp = [&](double k) { + if (std::isfinite(lo) && k < lo) { + k = lo; + } + if (std::isfinite(hi) && k > hi) { + k = hi; + } + return k; + }; + + std::vector ks; + auto add = [&](double k) { ks.push_back(clamp(std::floor(k + 0.5))); }; + + add(0.0); + add(1.0); + add(2.0); + add(mean); + for (double t : discreteTMultipliers()) { + add(mean - t * sd); + add(mean + t * sd); + } + if (std::isfinite(lo)) { + add(lo); + add(lo + 1.0); + add(lo + 2.0); + } + if (std::isfinite(hi)) { + add(hi); + add(hi - 1.0); + add(hi - 2.0); + } + + std::sort(ks.begin(), ks.end()); + ks.erase(std::unique(ks.begin(), ks.end()), ks.end()); + return ks; +} + +//============================================================================== +// Generic row emission (identical method surface across all 19 distributions) +//============================================================================== + +template +void emitPdfLogpdfCdfRows(Dist& dist, const std::string& distName, double p1, double p2, + bool hasP2, std::vector xs, bool includeSpecials, Sink& sink) { + if (includeSpecials) { + xs.push_back(std::numeric_limits::quiet_NaN()); + xs.push_back(std::numeric_limits::infinity()); + xs.push_back(-std::numeric_limits::infinity()); + } + + using ScalarFn = double (Dist::*)(double) const; + using BatchFn = void (Dist::*)(std::span, std::span, + const stats::detail::PerformanceHint&) const; + + struct Method { + const char* name; + ScalarFn scalarFn; + BatchFn batchFn; + }; + const std::array methods = { + Method{"pdf", static_cast(&Dist::getProbability), + static_cast(&Dist::getProbability)}, + Method{"logpdf", static_cast(&Dist::getLogProbability), + static_cast(&Dist::getLogProbability)}, + Method{"cdf", static_cast(&Dist::getCumulativeProbability), + static_cast(&Dist::getCumulativeProbability)}, + }; + + stats::detail::PerformanceHint hint; + hint.strategy = stats::detail::PerformanceHint::PreferredStrategy::FORCE_VECTORIZED; + + std::vector batchOut(xs.size()); + for (const auto& m : methods) { + std::fill(batchOut.begin(), batchOut.end(), 0.0); + (dist.*(m.batchFn))(std::span(xs), std::span(batchOut), hint); + for (std::size_t i = 0; i < xs.size(); ++i) { + double scalar = (dist.*(m.scalarFn))(xs[i]); + sink.addRow(makeRow(distName, m.name, p1, p2, hasP2, xs[i], scalar, batchOut[i], + /*hasBatch=*/true)); + } + } +} + +template +void emitQuantileRows(Dist& dist, const std::string& distName, double p1, double p2, bool hasP2, + Sink& sink) { + for (double p : quantilePGrid()) { + try { + double q = dist.getQuantile(p); + sink.addRow( + makeRow(distName, "quantile", p1, p2, hasP2, p, q, 0.0, /*hasBatch=*/false)); + } catch (...) { + sink.addSkip(distName); + } + } +} + +// Drives one continuous distribution across its 3 parameter instances. +template +void sweepContinuous(const std::string& distName, Factory factory, + const std::array, 3>& instances, bool hasP2, + bool boundedBothSides, Sink& sink) { + for (auto [p1, p2] : instances) { + Dist dist = factory(p1, p2); + std::vector xs = buildContinuousXGrid(dist, boundedBothSides); + emitPdfLogpdfCdfRows(dist, distName, p1, p2, hasP2, xs, /*includeSpecials=*/true, sink); + emitQuantileRows(dist, distName, p1, p2, hasP2, sink); + } +} + +// Drives one discrete distribution across its 3 parameter instances. +template +void sweepDiscrete(const std::string& distName, Factory factory, + const std::array, 3>& instances, bool hasP2, + Sink& sink) { + for (auto [p1, p2] : instances) { + Dist dist = factory(p1, p2); + std::vector ks = buildDiscreteKGrid(dist); + emitPdfLogpdfCdfRows(dist, distName, p1, p2, hasP2, ks, /*includeSpecials=*/false, sink); + emitQuantileRows(dist, distName, p1, p2, hasP2, sink); + } +} + +//============================================================================== +// Per-distribution instance tables (typical / small-boundary / stressed) and +// factories. Parameter values are chosen within each header's documented valid +// range (see the @param Doxygen comments checked against directly). +//============================================================================== + +void sweepAll(Sink& sink) { + using stats::BetaDistribution; + using stats::BinomialDistribution; + using stats::CauchyDistribution; + using stats::ChiSquaredDistribution; + using stats::DiscreteDistribution; + using stats::ExponentialDistribution; + using stats::GammaDistribution; + using stats::GaussianDistribution; + using stats::GeometricDistribution; + using stats::LaplaceDistribution; + using stats::LogNormalDistribution; + using stats::NegativeBinomialDistribution; + using stats::ParetoDistribution; + using stats::PoissonDistribution; + using stats::RayleighDistribution; + using stats::StudentTDistribution; + using stats::UniformDistribution; + using stats::VonMisesDistribution; + using stats::WeibullDistribution; + + // gaussian(mean, sigma): sigma > 0 + sweepContinuous( + "gaussian", [](double a, double b) { return GaussianDistribution(a, b); }, + {{{0.0, 1.0}, {1e-3, 1e-3}, {1e6, 1e4}}}, true, false, sink); + + // lognormal(mu, sigma): sigma > 0 + sweepContinuous( + "lognormal", [](double a, double b) { return LogNormalDistribution(a, b); }, + {{{0.0, 1.0}, {0.0, 0.01}, {0.0, 2.0}}}, true, false, sink); + + // exponential(lambda): lambda > 0 + sweepContinuous( + "exponential", [](double a, double) { return ExponentialDistribution(a); }, + {{{1.0, 0.0}, {1e-3, 0.0}, {1e6, 0.0}}}, false, false, sink); + + // uniform(a, b): b > a + sweepContinuous( + "uniform", [](double a, double b) { return UniformDistribution(a, b); }, + {{{0.0, 1.0}, {-1e-3, 1e-3}, {-1e8, 1e8}}}, true, true, sink); + + // gamma(alpha, beta): both > 0 (beta is a rate parameter) + sweepContinuous( + "gamma", [](double a, double b) { return GammaDistribution(a, b); }, + {{{2.0, 1.0}, {0.01, 0.01}, {1e4, 1e-3}}}, true, false, sink); + + // student_t(nu): nu > 0 + sweepContinuous( + "student_t", [](double a, double) { return StudentTDistribution(a); }, + {{{5.0, 0.0}, {1.001, 0.0}, {1e6, 0.0}}}, false, false, sink); + + // cauchy(x0, gamma): gamma > 0 + sweepContinuous( + "cauchy", [](double a, double b) { return CauchyDistribution(a, b); }, + {{{0.0, 1.0}, {0.0, 1e-6}, {1e8, 1e6}}}, true, false, sink); + + // von_mises(mu, kappa): kappa >= 0 + sweepContinuous( + "von_mises", [](double a, double b) { return VonMisesDistribution(a, b); }, + {{{0.0, 1.0}, {0.0, 1e-6}, {0.0, 100.0}}}, true, false, sink); + + // beta(alpha, beta): both > 0; support fixed at [0, 1] + sweepContinuous( + "beta", [](double a, double b) { return BetaDistribution(a, b); }, + {{{2.0, 3.0}, {0.01, 0.01}, {1e4, 1e4}}}, true, true, sink); + + // chi_squared(k): k > 0 + sweepContinuous( + "chi_squared", [](double a, double) { return ChiSquaredDistribution(a); }, + {{{3.0, 0.0}, {0.01, 0.0}, {1e5, 0.0}}}, false, false, sink); + + // laplace(mu, b): b > 0 + sweepContinuous( + "laplace", [](double a, double b) { return LaplaceDistribution(a, b); }, + {{{0.0, 1.0}, {0.0, 1e-6}, {1e8, 1e6}}}, true, false, sink); + + // pareto(scale=x_m, alpha): both > 0 + sweepContinuous( + "pareto", [](double a, double b) { return ParetoDistribution(a, b); }, + {{{1.0, 2.0}, {1e-6, 0.01}, {1e6, 100.0}}}, true, false, sink); + + // rayleigh(sigma): sigma > 0 + sweepContinuous( + "rayleigh", [](double a, double) { return RayleighDistribution(a); }, + {{{1.0, 0.0}, {1e-6, 0.0}, {1e6, 0.0}}}, false, false, sink); + + // weibull(shape, scale): both > 0 + sweepContinuous( + "weibull", [](double a, double b) { return WeibullDistribution(a, b); }, + {{{1.5, 1.0}, {0.01, 1e-3}, {100.0, 1e4}}}, true, false, sink); + + // --- Discrete distributions --- + + // poisson(lambda): lambda > 0 + sweepDiscrete( + "poisson", [](double a, double) { return PoissonDistribution(a); }, + {{{4.0, 0.0}, {1e-3, 0.0}, {1e5, 0.0}}}, false, sink); + + // discrete(a, b): int bounds, b > a; p1/p2 store the double bits of the ints used + sweepDiscrete( + "discrete", + [](double a, double b) { + return DiscreteDistribution(static_cast(a), static_cast(b)); + }, + {{{0.0, 9.0}, {0.0, 1.0}, {-1000000.0, 1000000.0}}}, true, sink); + + // binomial(n, p): n positive int, p in [0,1]; p1 stores double bits of the int n used + sweepDiscrete( + "binomial", + [](double a, double b) { return BinomialDistribution(static_cast(a), b); }, + {{{20.0, 0.5}, {1.0, 0.5}, {1000000.0, 0.3}}}, true, sink); + + // negative_binomial(r, p): r > 0 (real-valued), p in (0,1] + sweepDiscrete( + "negative_binomial", + [](double a, double b) { return NegativeBinomialDistribution(a, b); }, + {{{5.0, 0.5}, {0.01, 0.9}, {10000.0, 0.01}}}, true, sink); + + // geometric(p): p in (0,1] + sweepDiscrete( + "geometric", [](double a, double) { return GeometricDistribution(a); }, + {{{0.3, 0.0}, {1.0, 0.0}, {1e-6, 0.0}}}, false, sink); +} + +//============================================================================== +// Environment metadata (commit / isa / date) +//============================================================================== + +std::string getGitCommitShort() { +#if defined(_WIN32) + FILE* pipe = _popen("git rev-parse --short HEAD 2>NUL", "r"); +#else + FILE* pipe = popen("git rev-parse --short HEAD 2>/dev/null", "r"); +#endif + if (!pipe) { + return "unknown"; + } + std::string result; + char buffer[128]; + while (fgets(buffer, sizeof(buffer), pipe) != nullptr) { + result += buffer; + } +#if defined(_WIN32) + _pclose(pipe); +#else + pclose(pipe); +#endif + while (!result.empty() && (result.back() == '\n' || result.back() == '\r')) { + result.pop_back(); + } + return result.empty() ? "unknown" : result; +} + +std::string getDateStamp() { + std::time_t now = std::time(nullptr); + std::tm tmv{}; +#if defined(_WIN32) + localtime_s(&tmv, &now); +#else + localtime_r(&now, &tmv); +#endif + std::ostringstream oss; + oss << std::put_time(&tmv, "%Y-%m-%d"); + return oss.str(); +} + +} // namespace + +int main(int argc, char** argv) { + if (argc != 2) { + std::cerr << "Usage: accuracy_sweep \n"; + return 1; + } + const std::string outPath = argv[1]; + + Sink sink; + sweepAll(sink); + + std::ofstream out(outPath, std::ios::binary | std::ios::trunc); + if (!out) { + std::cerr << "Failed to open output file: " << outPath << "\n"; + return 1; + } + + out << "# libstats accuracy_sweep v1\n"; + out << "# commit=" << getGitCommitShort() + << " isa=" << stats::arch::simd::SIMDPolicy::getLevelString() + << " date=" << getDateStamp() << "\n"; + out << "dist,method,p1_bits,p2_bits,x_bits,scalar_bits,batch_bits\n"; + for (const auto& row : sink.rows()) { + out << row.toCsvLine() << "\n"; + } + for (const auto& [distName, n] : sink.skips()) { + if (n > 0) { + out << "# skipped_quantile " << distName << " " << n << "\n"; + } + } + out.close(); + + std::size_t total = sink.rows().size(); + std::cout << "accuracy_sweep: wrote " << total << " rows to " << outPath << "\n"; + std::cout << "Per-distribution row counts:\n"; + for (const auto& [distName, n] : sink.counts()) { + std::cout << " " << distName << ": " << n << "\n"; + } + bool anySkips = false; + for (const auto& [distName, n] : sink.skips()) { + if (n > 0) { + if (!anySkips) { + std::cout << "Skipped quantile rows (getQuantile threw):\n"; + anySkips = true; + } + std::cout << " " << distName << ": " << n << "\n"; + } + } + if (!anySkips) { + std::cout << "No quantile rows skipped.\n"; + } + + return 0; +} diff --git a/tools/accuracy_vs_mpmath.py b/tools/accuracy_vs_mpmath.py new file mode 100644 index 0000000..d51453f --- /dev/null +++ b/tools/accuracy_vs_mpmath.py @@ -0,0 +1,1795 @@ +#!/usr/bin/env python3 +"""Compare libstats' accuracy_sweep output against mpmath (dps=50) reference +implementations, for issue #46 (accuracy characterization across all 19 +libstats distributions). + +This is the oracle/report half of #46: the accuracy_sweep C++ tool (owned by +a sibling change, tools/accuracy_sweep.cpp) enumerates (dist, method, +params, x, scalar_result, batch_result) rows to stdout/CSV; this script is +the independent reference and report generator that consumes that CSV. + +INPUT CONTRACT -- one row per (dist, method, params, x) sample, comma +separated, no header required beyond the informational `#`-comment banner +the sweep tool emits (commit/isa/date). Columns: + + dist,method,p1_bits,p2_bits,x_bits,scalar_bits,batch_bits + + - All *_bits fields are 0x-hex uint64 IEEE-754 double bit patterns. + - method in {pdf, logpdf, cdf, quantile}. For quantile, x_bits carries the + probability p (not a domain point), and batch_bits may be the literal + string "-" (no batch quantile path). + - dist is the header basename (see DIST_PARAMS below for the parameter + mapping used here, read from each include/libstats/distributions/.h + Doxygen constructor comment). + - p1_bits/p2_bits are the two constructor parameters in DECLARATION order; + single-parameter distributions still carry a p2_bits field (its value is + unused by the reference functions below). + - Lines starting with `#` (the banner, and any trailing + `# skipped_quantile ...` notes) are comments and are skipped. + +ORACLE DOCTRINE (matches this repo's established generators -- see +scripts/gen_vonmises_cdf_vectors.py, scripts/gen_lognormal_cdf_vectors.py, +scripts/gen_gaussian_cdf_vectors.py, and tests/test_lognormal_cdf_accuracy.cpp's +law_budget comment, all read before writing this file): + + - mpmath at mp.dps = 50 (wider than the 40-digit CDF generators above, + since this script also differentiates for pdf/logpdf and root-finds for + quantile, both of which cost a few extra guard digits). + - Gaussian/LogNormal CDF: erfc closed form, exactly as those two + generators compute it (F = erfc(-z/sqrt(2))/2). + - Von Mises CDF: direct quadrature of the density (independent of the + library's Bessel-series/trapezoidal implementation), mu-centered wrap + replicated in DOUBLE arithmetic exactly as gen_vonmises_cdf_vectors.py's + wrap_to_pi -- reusing that file's PI_D anchoring: the wrap threshold + must be the double-rounded value of pi (math.pi), NOT mpmath's wider + mp.pi, or exact-endpoint rows spuriously disagree by a full branch-cut + fold. See wrap_to_pi's docstring below for the copied rationale. + - Discrete distributions (poisson, binomial, negative_binomial, + geometric, discrete): pmf/cdf via mpmath's regularized incomplete + gamma/beta identities where available (closed form), direct summation + only where no closed form is simpler; quantile by monotone search + against the closed-form CDF. + - Deep-tail CDF rows (reference F < 1e-3) additionally report law_frac: + the observed relative error as a fraction of the achievable-accuracy + law rel(F) ~ |ln F| * 2^-52 established in + tests/test_lognormal_cdf_accuracy.cpp's law_budget comment -- a flat + relative-error budget is unachievable in the deep tail for ANY + double-precision implementation of these closed forms, so a large + max_rel alongside a small law_frac (<= ~1) is the law, not a defect. + - NaN/Inf rows are contract-checked (NaN in -> NaN out, etc.) separately + from the error statistics, never folded into max_abs/max_rel/p99_rel. + +ORACLE SELF-CHECKS run unconditionally before any row is compared (see +`run_self_checks`) and this script exits non-zero if any fails -- a broken +oracle must never silently produce a report. + +Usage: + tools/accuracy_vs_mpmath.py [--out PATH] + +Writes the stdout summary table, then rewrites the generated block of +docs/ACCURACY_CHARACTERIZATION.md (between the `` +and `` markers) with per-distribution tables built +from the same statistics. +""" + +from __future__ import annotations + +import argparse +import math +import os +import struct +import sys +import time +from dataclasses import dataclass, field + +import mpmath as mp + +mp.mp.dps = 50 + +D = struct.Struct(" float helpers (matches the convention used by every gen_*.py +# generator in scripts/). +# --------------------------------------------------------------------------- + + +def bits_to_f64(b: int) -> float: + return D.unpack(Q.pack(b))[0] + + +def f64_to_bits(x: float) -> int: + return Q.unpack(D.pack(x))[0] + + +def parse_bits_field(s: str) -> "float | None": + """Parse a 0x-hex uint64 bit pattern into a double, or None for the + literal '-' (batch_bits sentinel for quantile rows with no batch path).""" + s = s.strip() + if s == "-": + return None + return bits_to_f64(int(s, 16)) + + +# --------------------------------------------------------------------------- +# Parameter mapping: constructor parameter order per include/libstats/ +# distributions/.h (read directly from each header's Doxygen comment +# before writing this table -- see the final report for any judgment calls). +# NPARAMS = 1 means p2_bits is present in the row but unused here. +# --------------------------------------------------------------------------- + +DIST_PARAMS = { + "gaussian": ("mean", "standardDeviation", 2), + "lognormal": ("mu", "sigma", 2), + "exponential": ("lambda", None, 1), + "uniform": ("a", "b", 2), + "poisson": ("lambda", None, 1), + "gamma": ("alpha_shape", "beta_rate", 2), + "discrete": ("a", "b", 2), # integer bounds, discrete uniform on {a,...,b} + "student_t": ("nu", None, 1), + "cauchy": ("x0", "gamma_scale", 2), + "von_mises": ("mu", "kappa", 2), + "binomial": ("n", "p", 2), # n is an integer trial count + "negative_binomial": ("r", "p", 2), + "geometric": ("p", None, 1), + "beta": ("alpha", "beta", 2), + "chi_squared": ("k", None, 1), + "laplace": ("mu", "b_scale", 2), + "pareto": ("scale_xm", "alpha", 2), + "rayleigh": ("sigma", None, 1), + "weibull": ("shape_k", "scale_lambda", 2), +} + +DISCRETE_DISTS = {"poisson", "binomial", "negative_binomial", "geometric", "discrete"} + +TWO_PI = 2 * mp.pi +SQRT2 = mp.sqrt(2) +LOG_2PI = mp.log(2 * mp.pi) + + +# --------------------------------------------------------------------------- +# Von Mises wrap: bit-exact DOUBLE replica of VonMisesDistribution::wrapAngle, +# copied (rationale included) from scripts/gen_vonmises_cdf_vectors.py. That +# generator's docstring explains at length why the wrap threshold must be +# the DOUBLE-rounded pi (math.pi, bit-identical to detail::PI) rather than +# mpmath's wider mp.pi: the library's own branch compare is a double compare +# against that same rounded constant, so anchoring the oracle's wrap on +# mp.pi instead would leave the exact double x = -PI_D just inside the +# interior of (-mp.pi, mp.pi] and disagree with the library's own fold-to- +# +PI_D convention at every exact-endpoint row -- a spurious ~1.0 CDF +# disagreement with nothing to do with kernel accuracy. This repo already +# litigated this seam (issue #51); reuse rather than rediscover it. +# --------------------------------------------------------------------------- + +PI_D = math.pi +TWO_PI_D = 2.0 * math.pi + + +def wrap_to_pi(x: float, mu: float) -> "mp.mpf": + t = math.fmod(x - mu, TWO_PI_D) + if t <= -PI_D: + t += TWO_PI_D + if t > PI_D: + t -= TWO_PI_D + return mp.mpf(t) + + +_vm_denom_cache: dict = {} + + +def _vm_circle_integral(kappa: "mp.mpf", lo: "mp.mpf", hi: "mp.mpf") -> "mp.mpf": + f = lambda p: mp.e ** (kappa * mp.cos(p)) + if lo <= 0 <= hi: + return mp.quad(f, [lo, mp.mpf(0), hi]) + return mp.quad(f, [lo, hi]) + + +def _vm_denom(kappa: "mp.mpf") -> "mp.mpf": + key = float(kappa) + if key not in _vm_denom_cache: + _vm_denom_cache[key] = _vm_circle_integral(kappa, -mp.pi, mp.pi) + return _vm_denom_cache[key] + + +def vonmises_cdf(mu: float, kappa: float, x: float) -> "mp.mpf": + kap = mp.mpf(kappa) + t = wrap_to_pi(x, mu) + num = _vm_circle_integral(kap, -mp.pi, t) + F = num / _vm_denom(kap) + if F < 0: + F = mp.mpf(0) + if F > 1: + F = mp.mpf(1) + return F + + +def vonmises_quantile(mu: float, kappa: float, p: float) -> "mp.mpf": + """Bisect t in [-pi, pi] on the mu=0-centered CDF for F(t) = p, then + shift by mu and wrap. This mirrors the SHAPE of the library's own + approach (build a mu=0 CDF, then shift/wrap for the actual mu) but not + its 2049-point grid + linear interpolation, so scalar_bits should be + expected to differ from this oracle by up to the grid's own + interpolation error, not just double-precision rounding -- see the + final report's judgment-call notes.""" + kap = mp.mpf(kappa) + Z = _vm_denom(kap) + pm = mp.mpf(p) + + def f(t): + return _vm_circle_integral(kap, -mp.pi, t) / Z - pm + + if pm <= 0: + t = -mp.pi + elif pm >= 1: + t = mp.pi + else: + # mpmath's bisect solver takes the bracket AS the start value; + # passing a scalar start plus an x0 kwarg is a TypeError. + t = mp.findroot(f, (-mp.pi, mp.pi), solver="bisect") + result = mp.mpf(mu) + t + # Wrap into (-pi, pi] at mpf precision (quantile has no library-side + # bit-exact double wrap to replicate; this is a reference convention). + while result > mp.pi: + result -= TWO_PI + while result <= -mp.pi: + result += TWO_PI + return result + + +# --------------------------------------------------------------------------- +# Reference implementations. Each entry: pdf(p1,p2,x), logpdf(p1,p2,x), +# cdf(p1,p2,x), quantile(p1,p2,p) -> mp.mpf. p1/p2 are passed as plain +# Python floats (the exact double values decoded from the row); each +# function lifts them to mpf itself. +# --------------------------------------------------------------------------- + + +def _tail_logspace_bisect(fn, target, seed, increasing=True, hi_cap=None): + """Root of fn(x) = target for fn strictly monotone on (0, hi_cap or inf), + by geometric-midpoint (log-space) bisection at mpf precision. + + Replaces the secant mp.findroot calls originally used for the gamma / + student-t / beta quantile oracles. A secant step that lands in a + deep-tail plateau -- where fn is constant to 50 digits -- stalls there + permanently (observed on the real #46 sweep: gamma quantile near + p = 1 - 1e-16 stuck at |fn - target| ~ 1e-15 against a 2.6e-54 + tolerance). Why not pure log-space bisection (the first replacement): + correct, but a fixed ~200 fn evaluations per row is hours of wall time + when fn is an incomplete beta with a huge first parameter -- student-t + nu = 1e6 costs tens of ms PER EVALUATION at dps 50, and minutes at the + extreme arguments the expansion probes. The hybrid below keeps + bisection's bracket invariant (cannot stall; anti-stagnation forcing + bounds it at ~2x bisection worst case) while false position in + (log x, log fn) coordinates converges superlinearly -- exactly, in one + step, wherever the tail is a power law, since F ~ C*x^a is a straight + line in those coordinates. Typical rows finish in 10-25 evaluations. + """ + seed = mp.mpf(seed) + if (not mp.isfinite(seed)) or seed <= 0: + seed = mp.mpf(1) + if hi_cap is not None and seed >= hi_cap: + seed = mp.mpf(hi_cap) / 2 + sgn = 1 if increasing else -1 + # Bracket expansion doubles the step each miss (with an exponent cap so + # mpf exponents stay small ints): a FIXED step cannot reach the bracket + # for small shape parameters -- gamma alpha=0.01 at p ~ 1e-300 has its + # quantile near 1e-30000, ~2800 fixed /4096 steps away but only ~13 + # doubling steps. Log-space bisection is indifferent to the overshoot: + # it converges in log-width, which the doubling only grows linearly. + step_cap = mp.mpf(2) ** 1000000 + lo = seed + step = mp.mpf(4096) + for _ in range(200): + if sgn * (fn(lo) - target) < 0: + break + lo /= step + step = min(step * step, step_cap) + else: + raise RuntimeError("quantile bracket expansion failed (low side)") + if hi_cap is not None: + hi = mp.mpf(hi_cap) + else: + hi = seed + step = mp.mpf(4096) + for _ in range(200): + if sgn * (fn(hi) - target) >= 0: + break + hi *= step + step = min(step * step, step_cap) + else: + raise RuntimeError("quantile bracket expansion failed (high side)") + # Safeguarded false position in u = log(x), h = log(fn) - log(target) + # coordinates. Bracket invariant: side(ulo) < 0 <= side(uhi) after sign + # folding. h is used only to PLACE the interpolated step (it is the + # near-linear coordinate in the tails); the plain sign decides which + # bracket end moves, so an underflowed fn (exact 0, no log) degrades + # a step to bisection instead of breaking anything. + ltarget = mp.log(target) + + def _eval(u): + v = fn(mp.exp(u)) + side = sgn * (v - target) + smooth = sgn * (mp.log(v) - ltarget) if v > 0 else None + return side, smooth + + ulo, uhi = mp.log(lo), mp.log(hi) + _, flo = _eval(ulo) + _, fhi = _eval(uhi) + tol = mp.mpf("1e-45") # in log-space = relative precision of x + last_moved = 0 # +1 lo moved, -1 hi moved; two in a row forces bisection + force_bisect = False + for _ in range(200): + width = uhi - ulo + if width < tol: + break + u_new = None + if not force_bisect and flo is not None and fhi is not None and fhi != flo: + u_try = ulo - flo * width / (fhi - flo) + margin = width / 64 + if ulo + margin < u_try < uhi - margin: + u_new = u_try + if u_new is None: + u_new = ulo + width / 2 + force_bisect = False + side_new, f_new = _eval(u_new) + if side_new < 0: + ulo, flo = u_new, f_new + force_bisect = last_moved == 1 + last_moved = 1 + else: + uhi, fhi = u_new, f_new + force_bisect = last_moved == -1 + last_moved = -1 + return mp.exp((ulo + uhi) / 2) + + +try: + from mpmath.libmp.libhyper import NoConvergence as _MPNoConvergence +except ImportError: # pragma: no cover - future mpmath relocation guard + + class _MPNoConvergence(Exception): + pass + + +def _gamma_cdf(alpha: "mp.mpf", beta: "mp.mpf", x: "mp.mpf") -> "mp.mpf": + """Regularized lower incomplete gamma P(alpha, beta*x), hardened for + large alpha the same way _betainc_reg is for large parameters: + + - far-lower-tail guard: leading term P ~ y^alpha e^-y / (alpha + Gamma(alpha)) below every positive double short-circuits the + special-function call entirely (and is the exact log-coordinate the + quantile solver interpolates on); + - right of the mean (y > alpha) the LOWER-gamma series converges too + slowly at large alpha and mpmath raises NoConvergence (observed at + chi-squared k = 1e5, quantile solve); the UPPER-gamma continued + fraction is the convergent representation there, so compute the + complement. NoConvergence is mpmath's own exception class, NOT a + ValueError -- catching it explicitly matters. + """ + if x <= 0: + return mp.mpf(0) + y = beta * x + if y > alpha: + # Right of the mean: complement via the upper-gamma continued + # fraction. P >= ~0.5 here, so the tiny-value guard never applies. + return 1 - mp.gammainc(alpha, y, mp.inf, regularized=True) + # Left of the mean only: the y -> 0 leading term is a valid guard. + lead = alpha * mp.log(y) - y - mp.log(alpha) - mp.loggamma(alpha) + if lead < -750: + return mp.exp(lead) + try: + return mp.gammainc(alpha, 0, y, regularized=True) + except (ValueError, _MPNoConvergence): + if lead < -100: + return mp.exp(lead) + raise + + +def _gamma_quantile(alpha: float, beta: float, p: float) -> "mp.mpf": + a = mp.mpf(alpha) + b = mp.mpf(beta) + pm = mp.mpf(p) + if pm <= 0: + return mp.mpf(0) + if pm >= 1: + return mp.inf + # Wilson-Hilferty approximation for the STANDARD (beta=1) gamma, used + # only to seed the bracket; deep-tail p rounds 2p-1 to +-1 at dps 50 + # and the erfinv blows up, so fall back to the mean and let the + # bracket expansion walk out. + z = mp.sqrt(2) * mp.erfinv(2 * pm - 1) + guess = a * (1 - 1 / (9 * a) + z / (3 * mp.sqrt(a))) ** 3 + if guess <= 0 or not mp.isfinite(guess): + guess = a + return _tail_logspace_bisect( + lambda xx: _gamma_cdf(a, b, xx), pm, guess / b) + + +def _discrete_search_quantile(cdf_int, lo: int, hi: int, p: "mp.mpf") -> int: + """Smallest integer k in [lo, hi] with cdf_int(k) >= p, by binary search + (cdf_int assumed monotone non-decreasing on integers).""" + if p <= 0: + return lo + while lo < hi: + mid = (lo + hi) // 2 + if cdf_int(mid) >= p: + hi = mid + else: + lo = mid + 1 + return lo + + +def _poisson_cdf(lam: "mp.mpf", k: int) -> "mp.mpf": + if k < 0: + return mp.mpf(0) + return mp.gammainc(k + 1, lam, mp.inf, regularized=True) + + +def _betainc_reg(a: "mp.mpf", b: "mp.mpf", x: "mp.mpf") -> "mp.mpf": + """Regularized incomplete beta I_x(a, b) with a far-tail asymptotic + guard for huge parameters. + + mp.betainc routes through hyp2f1, and for a huge first parameter at an + argument far from the transition region hypercomb FAILS TO CONVERGE -- + it escalates working precision (observed 189 -> 4577+ bits, with + million-bit exponents in play) for minutes and then raises ValueError. + That was the #46 oracle stall: student-t nu=1e6 probing t=100 dies + inside betainc(5e5, 0.5, 0, 0.99). + + Guard: the leading term of I_x(a,b) as x -> 0 is x^a / (a*B(a,b)), so + lead = a*ln x - ln a - ln B(a,b). When lead < -750 (value < ~1e-326, + below every positive double, including denormals -- no sweep target or + library-representable reference can live there) return exp(lead) + instead of calling betainc at all. The omitted 2F1 factor is O(1)-ish + there; irrelevant at 300+ orders of magnitude below any comparison, + while exp(lead) stays monotone for the quantile solvers' bracketing + and IS the correct leading log-coordinate their false-position phase + interpolates on. If betainc still raises inside the guard boundary + with lead < -100, fall back to the same asymptotic; a failure with a + non-tiny lead is a genuine oracle problem and re-raises. + """ + if x <= 0: + return mp.mpf(0) + if x >= 1: + return mp.mpf(1) + lead = a * mp.log(x) - mp.log(a) - ( + mp.loggamma(a) + mp.loggamma(b) - mp.loggamma(a + b)) + if lead < -750: + return mp.exp(lead) + if min(a, b) >= 5000: + # BOTH parameters large: mp.betainc dies here too, even in the + # central region where the value is ~0.5 and the lead guard + # rightly does not fire (observed: betainc(7e5, 3e5, 0, 0.7) + # hangs hyp2f1 for 300+ s). The continued fraction converges in + # ~O(sqrt(ab/(a+b))) cheap mpf iterations in exactly that regime. + return _betainc_cf(a, b, x) + try: + return mp.betainc(a, b, 0, x, regularized=True) + except (ValueError, _MPNoConvergence): + if lead < -100: + return mp.exp(lead) + return _betainc_cf(a, b, x) + + +def _betainc_cf(a: "mp.mpf", b: "mp.mpf", x: "mp.mpf") -> "mp.mpf": + """I_x(a, b) by the standard continued-fraction expansion (modified + Lentz iteration), written directly from the textbook formula: + + I_x(a,b) = x^a (1-x)^b / (a B(a,b)) * 1 / (1 + d1/(1 + d2/(1 + ...))) + d_{2m} = m (b - m) x / ((a + 2m - 1)(a + 2m)) + d_{2m+1} = -(a + m)(a + b + m) x / ((a + 2m)(a + 2m + 1)) + + valid (rapidly convergent) for x < (a+1)/(a+b+2); the complement + I_x(a,b) = 1 - I_{1-x}(b,a) covers the other side. Runs at elevated + working precision so the ambient dps-50 result keeps full accuracy. + Self-checked against mp.betainc on moderate parameters and against + the exact symmetry point I_{1/2}(a,a) = 1/2 in run_self_checks. + """ + if x > (a + 1) / (a + b + 2): + return 1 - _betainc_cf(b, a, 1 - x) + with mp.workprec(mp.mp.prec + 40): + tiny = mp.mpf(2) ** (-2 * mp.mp.prec) + eps = mp.mpf(2) ** (-mp.mp.prec + 5) + # Modified Lentz for the CF part. + c = mp.mpf(1) + d = mp.mpf(1) - (a + b) * x / (a + 1) + if abs(d) < tiny: + d = tiny + d = 1 / d + h = d + converged = False + for m in range(1, 100000): + m2 = 2 * m + num = m * (b - m) * x / ((a + m2 - 1) * (a + m2)) + d = 1 + num * d + if abs(d) < tiny: + d = tiny + c = 1 + num / c + if abs(c) < tiny: + c = tiny + d = 1 / d + h *= d * c + num = -(a + m) * (a + b + m) * x / ((a + m2) * (a + m2 + 1)) + d = 1 + num * d + if abs(d) < tiny: + d = tiny + c = 1 + num / c + if abs(c) < tiny: + c = tiny + d = 1 / d + delta = d * c + h *= delta + if abs(delta - 1) < eps: + converged = True + break + if not converged: + raise RuntimeError( + f"incomplete beta CF failed to converge (a={a}, b={b}, x={x})") + log_pre = (a * mp.log(x) + b * mp.log(1 - x) - mp.log(a) + - (mp.loggamma(a) + mp.loggamma(b) - mp.loggamma(a + b))) + result = mp.exp(log_pre) * h + return +result + + +def _binomial_cdf(n: int, p: "mp.mpf", k: int) -> "mp.mpf": + if k < 0: + return mp.mpf(0) + if k >= n: + return mp.mpf(1) + return _betainc_reg(mp.mpf(n - k), mp.mpf(k + 1), 1 - p) + + +def _negbinom_cdf(r: "mp.mpf", p: "mp.mpf", k: int) -> "mp.mpf": + if k < 0: + return mp.mpf(0) + return _betainc_reg(r, mp.mpf(k + 1), p) + + +def _student_t_cdf(nu: "mp.mpf", t: "mp.mpf") -> "mp.mpf": + if t == 0: + return mp.mpf("0.5") + xt = nu / (nu + t * t) + ib = _betainc_reg(nu / 2, mp.mpf("0.5"), xt) + return mp.mpf("0.5") * ib if t < 0 else 1 - mp.mpf("0.5") * ib + + +def _student_t_quantile(nu: "mp.mpf", p: "mp.mpf") -> "mp.mpf": + if p == mp.mpf("0.5"): + return mp.mpf(0) + guess = SQRT2 * mp.erfinv(2 * p - 1) # normal approx, seeds the bracket + if not mp.isfinite(guess): + # Deep tail: 2p-1 rounds to +-1 at dps 50 (p ~ 1e-300) and erfinv + # blows up. Use the asymptotic normal quantile z^2 ~ 2L - ln(2*pi* + # z^2), L = -ln(min(p, 1-p)), instead: without a finite seed the + # bracket expansion starts at 1 and its doubling steps probe + # extreme |t| where the large-nu incomplete beta costs MINUTES per + # evaluation (the observed #46 stall at nu = 1e6). A ~1% seed is + # plenty: expansion only needs to straddle the root, and the + # false-position phase converges from there in a handful of steps. + # (For small nu the true t-quantile is far beyond this normal-scale + # seed -- t ~ 1e300 at nu=1 -- and the doubling expansion walks out + # to it cheaply, small-nu evaluations being fast.) + tail = p if p < mp.mpf("0.5") else 1 - p + big_l = -mp.log(tail) + z2 = 2 * big_l + z2 = 2 * big_l - mp.log(2 * mp.pi * z2) + guess = mp.sqrt(z2) + if p < mp.mpf("0.5"): + guess = -guess + if p > mp.mpf("0.5"): + # t* > 0; F is increasing from 0.5 to 1 on (0, inf). + return _tail_logspace_bisect( + lambda t: _student_t_cdf(nu, t), p, guess) + # p < 0.5: solve on the negative axis directly through the + # cancellation-free incomplete-beta branch (F(-u) = 0.5*I(...) is + # DECREASING in u on (0, inf), range (0, 0.5)). Never via the + # symmetry Q(p) = -Q(1-p): for p ~ 1e-300, 1-p rounds to exactly 1 + # at dps 50 and the reflected problem degenerates. + return -_tail_logspace_bisect( + lambda u: _student_t_cdf(nu, -u), p, -guess, increasing=False) + + +REFS: dict = {} + + +def _reg(name): + def deco(cls): + REFS[name] = cls + return cls + + return deco + + +class Ref: + """Base: subclasses set pdf/logpdf/cdf/quantile as staticmethods taking + (p1, p2, x) -> mp.mpf (x is p for quantile).""" + + +@_reg("gaussian") +class GaussianRef(Ref): + @staticmethod + def pdf(mean, sigma, x): + m, s = mp.mpf(mean), mp.mpf(sigma) + z = (mp.mpf(x) - m) / s + return mp.e ** (-z * z / 2) / (s * mp.sqrt(2 * mp.pi)) + + @staticmethod + def logpdf(mean, sigma, x): + m, s = mp.mpf(mean), mp.mpf(sigma) + z = (mp.mpf(x) - m) / s + return -z * z / 2 - mp.log(s) - LOG_2PI / 2 + + @staticmethod + def cdf(mean, sigma, x): + z = (mp.mpf(x) - mp.mpf(mean)) / mp.mpf(sigma) + return mp.erfc(-z / SQRT2) / 2 + + @staticmethod + def quantile(mean, sigma, p): + return mp.mpf(mean) + mp.mpf(sigma) * SQRT2 * mp.erfinv(2 * mp.mpf(p) - 1) + + +@_reg("lognormal") +class LogNormalRef(Ref): + @staticmethod + def pdf(mu, sigma, x): + xm, mum, sm = mp.mpf(x), mp.mpf(mu), mp.mpf(sigma) + if xm <= 0: + return mp.mpf(0) + z = (mp.log(xm) - mum) / sm + return mp.e ** (-z * z / 2) / (xm * sm * mp.sqrt(2 * mp.pi)) + + @staticmethod + def logpdf(mu, sigma, x): + xm, mum, sm = mp.mpf(x), mp.mpf(mu), mp.mpf(sigma) + if xm <= 0: + return -mp.inf + z = (mp.log(xm) - mum) / sm + return -z * z / 2 - mp.log(xm) - mp.log(sm) - LOG_2PI / 2 + + @staticmethod + def cdf(mu, sigma, x): + xm, mum, sm = mp.mpf(x), mp.mpf(mu), mp.mpf(sigma) + if xm <= 0: + return mp.mpf(0) + z = (mp.log(xm) - mum) / sm + return mp.erfc(-z / SQRT2) / 2 + + @staticmethod + def quantile(mu, sigma, p): + z = SQRT2 * mp.erfinv(2 * mp.mpf(p) - 1) + return mp.e ** (mp.mpf(mu) + mp.mpf(sigma) * z) + + +@_reg("exponential") +class ExponentialRef(Ref): + @staticmethod + def pdf(lam, _p2, x): + l, xm = mp.mpf(lam), mp.mpf(x) + return l * mp.e ** (-l * xm) if xm >= 0 else mp.mpf(0) + + @staticmethod + def logpdf(lam, _p2, x): + l, xm = mp.mpf(lam), mp.mpf(x) + return mp.log(l) - l * xm if xm >= 0 else -mp.inf + + @staticmethod + def cdf(lam, _p2, x): + l, xm = mp.mpf(lam), mp.mpf(x) + return 1 - mp.e ** (-l * xm) if xm >= 0 else mp.mpf(0) + + @staticmethod + def quantile(lam, _p2, p): + return -mp.log(1 - mp.mpf(p)) / mp.mpf(lam) + + +@_reg("uniform") +class UniformRef(Ref): + @staticmethod + def pdf(a, b, x): + a, b, xm = mp.mpf(a), mp.mpf(b), mp.mpf(x) + return 1 / (b - a) if a <= xm <= b else mp.mpf(0) + + @staticmethod + def logpdf(a, b, x): + a, b, xm = mp.mpf(a), mp.mpf(b), mp.mpf(x) + return -mp.log(b - a) if a <= xm <= b else -mp.inf + + @staticmethod + def cdf(a, b, x): + a, b, xm = mp.mpf(a), mp.mpf(b), mp.mpf(x) + if xm < a: + return mp.mpf(0) + if xm > b: + return mp.mpf(1) + return (xm - a) / (b - a) + + @staticmethod + def quantile(a, b, p): + return mp.mpf(a) + mp.mpf(p) * (mp.mpf(b) - mp.mpf(a)) + + +@_reg("gamma") +class GammaRef(Ref): + @staticmethod + def pdf(alpha, beta, x): + a, b, xm = mp.mpf(alpha), mp.mpf(beta), mp.mpf(x) + if xm <= 0: + return mp.mpf(0) + return b**a / mp.gamma(a) * xm ** (a - 1) * mp.e ** (-b * xm) + + @staticmethod + def logpdf(alpha, beta, x): + a, b, xm = mp.mpf(alpha), mp.mpf(beta), mp.mpf(x) + if xm <= 0: + return -mp.inf + return a * mp.log(b) - mp.loggamma(a) + (a - 1) * mp.log(xm) - b * xm + + @staticmethod + def cdf(alpha, beta, x): + return _gamma_cdf(mp.mpf(alpha), mp.mpf(beta), mp.mpf(x)) + + @staticmethod + def quantile(alpha, beta, p): + return _gamma_quantile(alpha, beta, p) + + +@_reg("discrete") +class DiscreteRef(Ref): + @staticmethod + def pdf(a, b, x): + a, b = int(round(a)), int(round(b)) + xm = mp.mpf(x) + n = b - a + 1 + return 1 / mp.mpf(n) if xm == mp.floor(xm) and a <= xm <= b else mp.mpf(0) + + @staticmethod + def logpdf(a, b, x): + p = DiscreteRef.pdf(a, b, x) + return mp.log(p) if p > 0 else -mp.inf + + @staticmethod + def cdf(a, b, x): + a, b = int(round(a)), int(round(b)) + xm = mp.mpf(x) + n = b - a + 1 + if xm < a: + return mp.mpf(0) + if xm >= b: + return mp.mpf(1) + return (mp.floor(xm) - a + 1) / mp.mpf(n) + + @staticmethod + def quantile(a, b, p): + a, b = int(round(a)), int(round(b)) + n = b - a + 1 + pm = mp.mpf(p) + if pm <= 0: + return mp.mpf(a) + if pm >= 1: + return mp.mpf(b) + k = a - 1 + int(mp.ceil(pm * n)) + return mp.mpf(min(max(k, a), b)) + + +@_reg("student_t") +class StudentTRef(Ref): + @staticmethod + def pdf(nu, _p2, x): + n, xm = mp.mpf(nu), mp.mpf(x) + return ( + mp.gamma((n + 1) / 2) + / (mp.sqrt(n * mp.pi) * mp.gamma(n / 2)) + * (1 + xm * xm / n) ** (-(n + 1) / 2) + ) + + @staticmethod + def logpdf(nu, _p2, x): + n, xm = mp.mpf(nu), mp.mpf(x) + return ( + mp.loggamma((n + 1) / 2) + - mp.loggamma(n / 2) + - mp.log(n * mp.pi) / 2 + - (n + 1) / 2 * mp.log(1 + xm * xm / n) + ) + + @staticmethod + def cdf(nu, _p2, x): + return _student_t_cdf(mp.mpf(nu), mp.mpf(x)) + + @staticmethod + def quantile(nu, _p2, p): + return _student_t_quantile(mp.mpf(nu), mp.mpf(p)) + + +@_reg("cauchy") +class CauchyRef(Ref): + @staticmethod + def pdf(x0, gamma, x): + x0, g, xm = mp.mpf(x0), mp.mpf(gamma), mp.mpf(x) + z = (xm - x0) / g + return 1 / (mp.pi * g * (1 + z * z)) + + @staticmethod + def logpdf(x0, gamma, x): + x0, g, xm = mp.mpf(x0), mp.mpf(gamma), mp.mpf(x) + z = (xm - x0) / g + return -mp.log(mp.pi * g) - mp.log1p(z * z) + + @staticmethod + def cdf(x0, gamma, x): + x0, g, xm = mp.mpf(x0), mp.mpf(gamma), mp.mpf(x) + return mp.mpf("0.5") + mp.atan((xm - x0) / g) / mp.pi + + @staticmethod + def quantile(x0, gamma, p): + return mp.mpf(x0) + mp.mpf(gamma) * mp.tan(mp.pi * (mp.mpf(p) - mp.mpf("0.5"))) + + +@_reg("von_mises") +class VonMisesRef(Ref): + @staticmethod + def pdf(mu, kappa, x): + m, k, xm = mp.mpf(mu), mp.mpf(kappa), mp.mpf(x) + return mp.e ** (k * mp.cos(xm - m)) / (2 * mp.pi * mp.besseli(0, k)) + + @staticmethod + def logpdf(mu, kappa, x): + m, k, xm = mp.mpf(mu), mp.mpf(kappa), mp.mpf(x) + return k * mp.cos(xm - m) - mp.log(2 * mp.pi * mp.besseli(0, k)) + + @staticmethod + def cdf(mu, kappa, x): + return vonmises_cdf(mu, kappa, x) + + @staticmethod + def quantile(mu, kappa, p): + return vonmises_quantile(mu, kappa, p) + + +@_reg("binomial") +class BinomialRef(Ref): + @staticmethod + def pdf(n, p, x): + n = int(round(n)) + k = int(round(x)) + pm = mp.mpf(p) + if k < 0 or k > n: + return mp.mpf(0) + return mp.binomial(n, k) * pm**k * (1 - pm) ** (n - k) + + @staticmethod + def logpdf(n, p, x): + v = BinomialRef.pdf(n, p, x) + return mp.log(v) if v > 0 else -mp.inf + + @staticmethod + def cdf(n, p, x): + n = int(round(n)) + k = int(math.floor(x)) + return _binomial_cdf(n, mp.mpf(p), k) + + @staticmethod + def quantile(n, p, prob): + n = int(round(n)) + pm = mp.mpf(p) + pr = mp.mpf(prob) + return mp.mpf(_discrete_search_quantile(lambda k: _binomial_cdf(n, pm, k), 0, n, pr)) + + +@_reg("negative_binomial") +class NegativeBinomialRef(Ref): + @staticmethod + def pdf(r, p, x): + r, pm = mp.mpf(r), mp.mpf(p) + k = int(round(x)) + if k < 0: + return mp.mpf(0) + return mp.gamma(k + r) / (mp.gamma(r) * mp.factorial(k)) * pm**r * (1 - pm) ** k + + @staticmethod + def logpdf(r, p, x): + v = NegativeBinomialRef.pdf(r, p, x) + return mp.log(v) if v > 0 else -mp.inf + + @staticmethod + def cdf(r, p, x): + k = int(math.floor(x)) + return _negbinom_cdf(mp.mpf(r), mp.mpf(p), k) + + @staticmethod + def quantile(r, p, prob): + rm, pm, pr = mp.mpf(r), mp.mpf(p), mp.mpf(prob) + # Doubling search for an upper bound, then binary search (no fixed + # support ceiling for negative binomial's failure count). + hi = 1 + while _negbinom_cdf(rm, pm, hi) < pr: + hi *= 2 + if hi > 1 << 40: + break + return mp.mpf(_discrete_search_quantile(lambda k: _negbinom_cdf(rm, pm, k), 0, hi, pr)) + + +@_reg("geometric") +class GeometricRef(Ref): + @staticmethod + def pdf(p, _p2, x): + pm = mp.mpf(p) + k = int(round(x)) + return pm * (1 - pm) ** k if k >= 0 else mp.mpf(0) + + @staticmethod + def logpdf(p, _p2, x): + pm = mp.mpf(p) + k = int(round(x)) + if k < 0: + return -mp.inf + if k == 0: + # k * log1p(-p) is exactly 0 here; spelling it out avoids the + # 0 * (-inf) = NaN artifact at p == 1 (pmf(0) = p exactly). + return mp.log(pm) + return mp.log(pm) + k * mp.log1p(-pm) + + @staticmethod + def cdf(p, _p2, x): + pm = mp.mpf(p) + k = int(math.floor(x)) + if k < 0: + return mp.mpf(0) + return 1 - (1 - pm) ** (k + 1) + + @staticmethod + def quantile(p, _p2, prob): + pm, pr = mp.mpf(p), mp.mpf(prob) + if pr <= 0: + return mp.mpf(0) + if pm >= 1: + return mp.mpf(0) + k = mp.ceil(mp.log1p(-pr) / mp.log1p(-pm) - 1) + k = max(k, 0) + # Closed form can be off by one at exact boundaries; nudge via the + # exact CDF (monotone, so a local search of +/-2 always corrects it). + while GeometricRef.cdf(p, _p2, k - 1) >= pr and k > 0: + k -= 1 + while GeometricRef.cdf(p, _p2, k) < pr: + k += 1 + return mp.mpf(k) + + +@_reg("beta") +class BetaRef(Ref): + @staticmethod + def _edge_pdf(edge_shape, other_shape): + # Density AT an endpoint: the endpoint's shape exponent governs. + # shape < 1 is an integrable singularity (+inf is the limit, and + # what the library returns); shape == 1 leaves the finite value + # 1/B(1, other) = other; shape > 1 pins the density to 0. + if edge_shape < 1: + return mp.inf + if edge_shape == 1: + return other_shape + return mp.mpf(0) + + @staticmethod + def pdf(alpha, beta, x): + a, b, xm = mp.mpf(alpha), mp.mpf(beta), mp.mpf(x) + if xm < 0 or xm > 1: + return mp.mpf(0) + if xm == 0: + return BetaRef._edge_pdf(a, b) + if xm == 1: + return BetaRef._edge_pdf(b, a) + return xm ** (a - 1) * (1 - xm) ** (b - 1) / mp.beta(a, b) + + @staticmethod + def logpdf(alpha, beta, x): + a, b, xm = mp.mpf(alpha), mp.mpf(beta), mp.mpf(x) + if xm < 0 or xm > 1: + return -mp.inf + if xm == 0 or xm == 1: + edge = BetaRef._edge_pdf(a if xm == 0 else b, b if xm == 0 else a) + return mp.log(edge) if edge > 0 else -mp.inf + return (a - 1) * mp.log(xm) + (b - 1) * mp.log1p(-xm) - mp.log(mp.beta(a, b)) + + @staticmethod + def cdf(alpha, beta, x): + a, b, xm = mp.mpf(alpha), mp.mpf(beta), mp.mpf(x) + if xm <= 0: + return mp.mpf(0) + if xm >= 1: + return mp.mpf(1) + return _betainc_reg(a, b, xm) + + @staticmethod + def quantile(alpha, beta, p): + a, b, pm = mp.mpf(alpha), mp.mpf(beta), mp.mpf(p) + if pm <= 0: + return mp.mpf(0) + if pm >= 1: + return mp.mpf(1) + + def q_lower(aa, bb, ppm): + # Lower-half solve on (0, 1]; hi_cap=1 is always a valid upper + # bracket (I_1 = 1 >= ppm), the low side log-expands toward 0. + return _tail_logspace_bisect( + lambda xx: _betainc_reg(aa, bb, xx), + ppm, aa / (aa + bb), hi_cap=1) + + if pm > mp.mpf("0.5"): + # Reflect: Q(a,b,p) = 1 - Q(b,a,1-p). Here 1-pm IS exact at + # dps 50 (pm is a lifted double in (0.5, 1), so 1-pm needs + # < 53 bits), unlike the p < 0.5 deep-tail direction -- and it + # moves the solve to the log-space-friendly lower tail. + return 1 - q_lower(b, a, 1 - pm) + return q_lower(a, b, pm) + + +@_reg("chi_squared") +class ChiSquaredRef(Ref): + @staticmethod + def pdf(k, _p2, x): + kk, xm = mp.mpf(k), mp.mpf(x) + if xm <= 0: + return mp.mpf(0) + return xm ** (kk / 2 - 1) * mp.e ** (-xm / 2) / (2 ** (kk / 2) * mp.gamma(kk / 2)) + + @staticmethod + def logpdf(k, _p2, x): + kk, xm = mp.mpf(k), mp.mpf(x) + if xm <= 0: + return -mp.inf + return (kk / 2 - 1) * mp.log(xm) - xm / 2 - (kk / 2) * mp.log(2) - mp.loggamma(kk / 2) + + @staticmethod + def cdf(k, _p2, x): + return _gamma_cdf(mp.mpf(k) / 2, mp.mpf("0.5"), mp.mpf(x)) + + @staticmethod + def quantile(k, _p2, p): + return _gamma_quantile(float(k) / 2.0, 0.5, p) + + +@_reg("laplace") +class LaplaceRef(Ref): + @staticmethod + def pdf(mu, b, x): + m, bb, xm = mp.mpf(mu), mp.mpf(b), mp.mpf(x) + return mp.e ** (-abs(xm - m) / bb) / (2 * bb) + + @staticmethod + def logpdf(mu, b, x): + m, bb, xm = mp.mpf(mu), mp.mpf(b), mp.mpf(x) + return -abs(xm - m) / bb - mp.log(2 * bb) + + @staticmethod + def cdf(mu, b, x): + m, bb, xm = mp.mpf(mu), mp.mpf(b), mp.mpf(x) + if xm < m: + return mp.mpf("0.5") * mp.e ** ((xm - m) / bb) + return 1 - mp.mpf("0.5") * mp.e ** (-(xm - m) / bb) + + @staticmethod + def quantile(mu, b, p): + m, bb, pm = mp.mpf(mu), mp.mpf(b), mp.mpf(p) + if pm < mp.mpf("0.5"): + return m + bb * mp.log(2 * pm) + return m - bb * mp.log(2 * (1 - pm)) + + +@_reg("pareto") +class ParetoRef(Ref): + @staticmethod + def pdf(scale, alpha, x): + xm0, a, xm = mp.mpf(scale), mp.mpf(alpha), mp.mpf(x) + if xm < xm0: + return mp.mpf(0) + return a * xm0**a / xm ** (a + 1) + + @staticmethod + def logpdf(scale, alpha, x): + xm0, a, xm = mp.mpf(scale), mp.mpf(alpha), mp.mpf(x) + if xm < xm0: + return -mp.inf + return mp.log(a) + a * mp.log(xm0) - (a + 1) * mp.log(xm) + + @staticmethod + def cdf(scale, alpha, x): + xm0, a, xm = mp.mpf(scale), mp.mpf(alpha), mp.mpf(x) + if xm < xm0: + return mp.mpf(0) + return 1 - (xm0 / xm) ** a + + @staticmethod + def quantile(scale, alpha, p): + xm0, a, pm = mp.mpf(scale), mp.mpf(alpha), mp.mpf(p) + return xm0 / (1 - pm) ** (1 / a) + + +@_reg("rayleigh") +class RayleighRef(Ref): + @staticmethod + def pdf(sigma, _p2, x): + s, xm = mp.mpf(sigma), mp.mpf(x) + if xm < 0: + return mp.mpf(0) + return xm / s**2 * mp.e ** (-xm * xm / (2 * s * s)) + + @staticmethod + def logpdf(sigma, _p2, x): + s, xm = mp.mpf(sigma), mp.mpf(x) + if xm <= 0: + return -mp.inf + return mp.log(xm) - 2 * mp.log(s) - xm * xm / (2 * s * s) + + @staticmethod + def cdf(sigma, _p2, x): + s, xm = mp.mpf(sigma), mp.mpf(x) + if xm < 0: + return mp.mpf(0) + return 1 - mp.e ** (-xm * xm / (2 * s * s)) + + @staticmethod + def quantile(sigma, _p2, p): + s, pm = mp.mpf(sigma), mp.mpf(p) + return s * mp.sqrt(-2 * mp.log(1 - pm)) + + +@_reg("weibull") +class WeibullRef(Ref): + @staticmethod + def pdf(shape, scale, x): + k, lam, xm = mp.mpf(shape), mp.mpf(scale), mp.mpf(x) + if xm < 0: + return mp.mpf(0) + return (k / lam) * (xm / lam) ** (k - 1) * mp.e ** (-((xm / lam) ** k)) + + @staticmethod + def logpdf(shape, scale, x): + k, lam, xm = mp.mpf(shape), mp.mpf(scale), mp.mpf(x) + if xm < 0: + return -mp.inf + if xm == 0: + # Density at 0: +inf for k < 1 (integrable singularity, what + # the library returns), k/lam = 1/lam for k == 1, 0 for k > 1. + if k < 1: + return mp.inf + if k == 1: + return -mp.log(lam) + return -mp.inf + return mp.log(k) - mp.log(lam) + (k - 1) * (mp.log(xm) - mp.log(lam)) - (xm / lam) ** k + + @staticmethod + def cdf(shape, scale, x): + k, lam, xm = mp.mpf(shape), mp.mpf(scale), mp.mpf(x) + if xm < 0: + return mp.mpf(0) + return 1 - mp.e ** (-((xm / lam) ** k)) + + @staticmethod + def quantile(shape, scale, p): + k, lam, pm = mp.mpf(shape), mp.mpf(scale), mp.mpf(p) + return lam * (-mp.log(1 - pm)) ** (1 / k) + + +@_reg("poisson") +class PoissonRef(Ref): + @staticmethod + def pdf(lam, _p2, x): + l = mp.mpf(lam) + k = int(round(x)) + if k < 0: + return mp.mpf(0) + return mp.e ** (-l) * l**k / mp.factorial(k) + + @staticmethod + def logpdf(lam, _p2, x): + l = mp.mpf(lam) + k = int(round(x)) + if k < 0: + return -mp.inf + return -l + k * mp.log(l) - mp.loggamma(k + 1) + + @staticmethod + def cdf(lam, _p2, x): + k = int(math.floor(x)) + return _poisson_cdf(mp.mpf(lam), k) + + @staticmethod + def quantile(lam, _p2, p): + l, pr = mp.mpf(lam), mp.mpf(p) + hi = max(1, int(float(l) * 2) + 10) + while _poisson_cdf(l, hi) < pr: + hi *= 2 + if hi > 1 << 40: + break + return mp.mpf(_discrete_search_quantile(lambda k: _poisson_cdf(l, k), 0, hi, pr)) + + +# --------------------------------------------------------------------------- +# Oracle self-checks. Run before any comparison; assertion failure -> the +# process exits non-zero (see main()). At least two independent known-value +# checks per distribution: a median/symmetry identity and a closed-form +# spot value, matching this repo's generator doctrine (house rule: "generated +# references are trusted over comments" -- so the CHECKS, not comments, are +# what's trusted here). +# --------------------------------------------------------------------------- + + +def run_self_checks() -> list: + """Returns the list of (name, ok, detail) tuples; raises AssertionError + immediately on the first failing check (fail-fast, matches the + generators' doctrine).""" + results = [] + + def check(name, cond, detail=""): + assert cond, f"SELF-CHECK FAILED: {name} {detail}" + results.append((name, True, detail)) + + tol = mp.mpf("1e-30") + + # --- incomplete-beta continued fraction (large-parameter path) --- + # The CF replaces mp.betainc wherever min(a,b) >= 5000 (see + # _betainc_reg); anchor it against mp.betainc where mpmath is healthy, + # and against the exact symmetry point I_{1/2}(a,a) = 1/2 in the + # large-parameter regime mpmath cannot reach. + for _a, _b, _x in ((2.5, 7.0, 0.2), (30.0, 4.0, 0.9), (100.0, 250.0, 0.31)): + _ref = mp.betainc(mp.mpf(_a), mp.mpf(_b), 0, mp.mpf(_x), regularized=True) + _cf = _betainc_cf(mp.mpf(_a), mp.mpf(_b), mp.mpf(_x)) + check( + f"betainc_cf.vs_mpmath({_a},{_b},{_x})", + abs(_cf - _ref) <= mp.mpf("1e-45") * _ref, + ) + check( + "betainc_cf.symmetry_point_large", + abs(_betainc_cf(mp.mpf(10000), mp.mpf(10000), mp.mpf("0.5")) - mp.mpf("0.5")) + <= mp.mpf("1e-45"), + ) + + # --- gaussian --- + check("gaussian.cdf(mean)==0.5", abs(GaussianRef.cdf(0.3, 2.0, 0.3) - mp.mpf("0.5")) <= tol) + check( + "gaussian.symmetry", + abs(GaussianRef.cdf(0.0, 1.0, 1.5) - (1 - GaussianRef.cdf(0.0, 1.0, -1.5))) <= tol, + ) + + # --- lognormal --- + check( + "lognormal.cdf(exp(mu))==0.5", + abs(LogNormalRef.cdf(0.5, 1.0, float(mp.e ** mp.mpf("0.5"))) - mp.mpf("0.5")) <= mp.mpf("1e-15"), + ) + check( + "lognormal.pdf_positive_support", + LogNormalRef.pdf(0.0, 1.0, -1.0) == 0 and LogNormalRef.pdf(0.0, 1.0, 1.0) > 0, + ) + + # --- exponential --- + check( + "exponential.cdf(1/lambda)==1-e^-1", + abs(ExponentialRef.cdf(2.0, 0, 0.5) - (1 - mp.e ** -1)) <= tol, + ) + check("exponential.quantile_inverts_cdf", abs(ExponentialRef.quantile(3.0, 0, ExponentialRef.cdf(3.0, 0, 0.7)) + - mp.mpf(0.7)) <= mp.mpf("1e-15")) + + # --- uniform --- + check("uniform.cdf(mid)==0.5", UniformRef.cdf(2.0, 6.0, 4.0) == mp.mpf("0.5")) + check("uniform.pdf_outside_zero", UniformRef.pdf(0.0, 1.0, 5.0) == 0) + + # --- poisson --- + check("poisson.pmf(0;lambda)==e^-lambda", abs(PoissonRef.pdf(2.5, 0, 0) - mp.e ** mp.mpf("-2.5")) <= tol) + check( + "poisson.cdf_sums_to_pmf", + abs(PoissonRef.cdf(3.0, 0, 2) - sum(PoissonRef.pdf(3.0, 0, k) for k in range(3))) <= tol, + ) + + # --- gamma --- + check( + "gamma.cdf(0)==0_and_cdf(inf)->1", + GammaRef.cdf(2.0, 1.0, 0.0) == 0 and GammaRef.cdf(2.0, 1.0, 1e6) > mp.mpf("0.999999"), + ) + check( + "gamma.quantile_inverts_cdf", + abs(GammaRef.quantile(3.0, 2.0, GammaRef.cdf(3.0, 2.0, 1.5)) - mp.mpf(1.5)) <= mp.mpf("1e-15"), + ) + + # --- discrete --- + check("discrete.pdf_uniform", DiscreteRef.pdf(2, 5, 3) == mp.mpf(1) / 4) + check("discrete.cdf(b)==1", DiscreteRef.cdf(2, 5, 5) == 1) + + # --- student_t --- + check("student_t.cdf(0)==0.5", StudentTRef.cdf(5.0, 0, 0.0) == mp.mpf("0.5")) + check( + "student_t.symmetry", + abs(StudentTRef.cdf(5.0, 0, 1.2) - (1 - StudentTRef.cdf(5.0, 0, -1.2))) <= tol, + ) + + # --- cauchy --- + check("cauchy.cdf(x0)==0.5", CauchyRef.cdf(1.0, 2.0, 1.0) == mp.mpf("0.5")) + check( + "cauchy.quantile_inverts_cdf", + abs(CauchyRef.quantile(0.0, 1.0, CauchyRef.cdf(0.0, 1.0, 3.0)) - mp.mpf(3.0)) <= mp.mpf("1e-15"), + ) + + # --- von_mises --- + check( + "vonmises.cdf(mu)==0.5", + abs(VonMisesRef.cdf(0.7, 2.0, 0.7) - mp.mpf("0.5")) <= mp.mpf("1e-25"), + ) + check( + "vonmises.kappa0_uniform_closed_form", + abs(VonMisesRef.cdf(0.0, 1e-12, 1.0) - (mp.mpf(1.0) + mp.pi) / TWO_PI) <= mp.mpf("1e-6"), + ) + + # --- binomial --- + check("binomial.pmf_sums_to_one", abs(sum(BinomialRef.pdf(5, 0.3, k) for k in range(6)) - 1) <= tol) + check("binomial.cdf(n)==1", BinomialRef.cdf(5, 0.3, 5) == 1) + + # --- negative_binomial --- + check( + "negbinom.pmf(0;r,p)==p^r", + abs(NegativeBinomialRef.pdf(3.0, 0.4, 0) - mp.mpf(0.4) ** 3) <= tol, + ) + check( + "negbinom.cdf_monotone", + NegativeBinomialRef.cdf(3.0, 0.4, 5) >= NegativeBinomialRef.cdf(3.0, 0.4, 2), + ) + + # --- geometric --- + check("geometric.pmf(0)==p", GeometricRef.pdf(0.3, 0, 0) == mp.mpf(0.3)) + check( + "geometric.cdf_closed_form", + abs(GeometricRef.cdf(0.3, 0, 4) - (1 - (1 - mp.mpf(0.3)) ** 5)) <= tol, + ) + + # --- beta --- + check("beta.cdf(0)==0_cdf(1)==1", BetaRef.cdf(2.0, 3.0, 0.0) == 0 and BetaRef.cdf(2.0, 3.0, 1.0) == 1) + check( + "beta.symmetric_alpha_eq_beta", + abs(BetaRef.cdf(2.0, 2.0, 0.5) - mp.mpf("0.5")) <= tol, + ) + + # --- chi_squared --- + check( + "chisq.equals_gamma_k_over_2_half", + abs(ChiSquaredRef.cdf(4.0, 0, 3.0) - GammaRef.cdf(2.0, 0.5, 3.0)) <= tol, + ) + check("chisq.cdf(0)==0", ChiSquaredRef.cdf(4.0, 0, 0.0) == 0) + + # --- laplace --- + check("laplace.cdf(mu)==0.5", LaplaceRef.cdf(1.0, 2.0, 1.0) == mp.mpf("0.5")) + check( + "laplace.symmetry", + abs(LaplaceRef.cdf(0.0, 1.0, 2.0) - (1 - LaplaceRef.cdf(0.0, 1.0, -2.0))) <= tol, + ) + + # --- pareto --- + check("pareto.cdf(scale)==0", ParetoRef.cdf(1.0, 3.0, 1.0) == 0) + check( + "pareto.quantile_inverts_cdf", + abs(ParetoRef.quantile(1.0, 3.0, ParetoRef.cdf(1.0, 3.0, 2.5)) - mp.mpf(2.5)) <= mp.mpf("1e-15"), + ) + + # --- rayleigh --- + check( + "rayleigh.cdf_closed_form", + abs(RayleighRef.cdf(2.0, 0, 2.0) - (1 - mp.e ** mp.mpf("-0.5"))) <= tol, + ) + check("rayleigh.cdf(0)==0", RayleighRef.cdf(2.0, 0, 0.0) == 0) + + # --- weibull --- + check( + "weibull.reduces_to_exponential_at_shape1", + abs(WeibullRef.cdf(1.0, 2.0, 3.0) - ExponentialRef.cdf(0.5, 0, 3.0)) <= tol, + ) + check("weibull.cdf(0)==0", WeibullRef.cdf(2.0, 1.0, 0.0) == 0) + + return results + + +# --------------------------------------------------------------------------- +# CSV parsing +# --------------------------------------------------------------------------- + + +@dataclass +class Row: + dist: str + method: str + p1: float + p2: float + x: float # for quantile rows, this is p + scalar: float + batch: "float | None" + lineno: int + + +def parse_csv(path: str) -> "tuple[list[Row], dict]": + rows: list[Row] = [] + meta: dict = {} + with open(path, "r", encoding="utf-8") as f: + for lineno, line in enumerate(f, start=1): + line = line.strip() + if not line: + continue + if line.startswith("#"): + if "commit=" in line and "isa=" in line: + meta["banner"] = line.lstrip("#").strip() + continue + if line == "dist,method,p1_bits,p2_bits,x_bits,scalar_bits,batch_bits": + # Literal column-header line the sweep tool emits (see the + # input contract in this file's module docstring) -- not a + # `#`-comment, but not data either. + continue + parts = line.split(",") + if len(parts) != 7: + raise ValueError(f"{path}:{lineno}: expected 7 fields, got {len(parts)}: {line!r}") + dist, method, p1s, p2s, xs, scs, bts = parts + p1 = parse_bits_field(p1s) + p2 = parse_bits_field(p2s) + x = parse_bits_field(xs) + scalar = parse_bits_field(scs) + batch = parse_bits_field(bts) + if p1 is None or x is None or scalar is None: + raise ValueError(f"{path}:{lineno}: p1/x/scalar must not be '-': {line!r}") + rows.append( + Row(dist.strip(), method.strip(), p1, p2 if p2 is not None else 0.0, x, scalar, batch, lineno) + ) + return rows, meta + + +# --------------------------------------------------------------------------- +# Comparison / statistics +# --------------------------------------------------------------------------- + + +def law_of_f(f_ref: "mp.mpf") -> "mp.mpf": + """Achievable double-precision relative-error law in the CDF's deep + tail: rel(F) ~ |ln F| * 2^-52. See tests/test_lognormal_cdf_accuracy.cpp's + law_budget comment -- this is the UNPADDED law (that test's pinned + budget adds >2x headroom on top of this for gating; this report uses + the bare law as the yardstick for "law-limited, not a defect").""" + if f_ref <= 0: + return TWO_POW_NEG52 * mp.mpf(20) + return abs(mp.log(f_ref)) * TWO_POW_NEG52 + + +def relerr(got: "mp.mpf", ref: "mp.mpf") -> "mp.mpf": + if abs(ref) < mp.mpf("1e-290"): + return abs(got - ref) + return abs(got - ref) / abs(ref) + + +@dataclass +class SourceStats: + n: int = 0 + max_abs: "mp.mpf" = field(default_factory=lambda: mp.mpf(0)) + max_rel: "mp.mpf" = field(default_factory=lambda: mp.mpf(0)) + rels: list = field(default_factory=list) + worst_x: float = 0.0 + max_law_frac: "mp.mpf" = field(default_factory=lambda: mp.mpf(0)) + law_worst_x: float = 0.0 + has_law_rows: bool = False + + def observe(self, got: "mp.mpf", ref: "mp.mpf", x: float, is_cdf_tail: bool): + a = abs(got - ref) + r = relerr(got, ref) + self.n += 1 + self.rels.append(r) + if a > self.max_abs: + self.max_abs = a + if r > self.max_rel: + self.max_rel = r + self.worst_x = x + if is_cdf_tail: + self.has_law_rows = True + lf = r / law_of_f(ref) + if lf > self.max_law_frac: + self.max_law_frac = lf + self.law_worst_x = x + + def p99_rel(self) -> "mp.mpf": + if not self.rels: + return mp.mpf(0) + s = sorted(self.rels) + idx = min(len(s) - 1, int(math.ceil(0.99 * len(s))) - 1) + return s[max(idx, 0)] + + +@dataclass +class GroupResult: + dist: str + method: str + scalar: SourceStats = field(default_factory=SourceStats) + batch: SourceStats = field(default_factory=SourceStats) + batch_vs_scalar_max_abs: "mp.mpf" = field(default_factory=lambda: mp.mpf(0)) + batch_vs_scalar_max_rel: "mp.mpf" = field(default_factory=lambda: mp.mpf(0)) + n_batch_rows: int = 0 + violations: list = field(default_factory=list) + + +def is_nan(v: float) -> bool: + return isinstance(v, float) and math.isnan(v) + + +def compare(rows: list) -> "tuple[dict, list]": + groups: dict = {} + skipped = [] + + n_done = 0 + for row in rows: + n_done += 1 + if n_done % 500 == 0: + print(f" [progress] {n_done}/{len(rows)} rows", file=sys.stderr, flush=True) + ref_fns = REFS.get(row.dist) + if ref_fns is None: + skipped.append((row.lineno, f"unknown dist {row.dist!r}")) + continue + fn = getattr(ref_fns, row.method, None) + if fn is None: + skipped.append((row.lineno, f"unknown method {row.method!r} for {row.dist}")) + continue + + key = (row.dist, row.method) + g = groups.setdefault(key, GroupResult(row.dist, row.method)) + + nan_input = is_nan(row.p1) or is_nan(row.p2) or is_nan(row.x) + inf_input = math.isinf(row.x) + + if nan_input: + # Contract check only: NaN in -> NaN out, both paths. + if not is_nan(row.scalar): + g.violations.append((row.lineno, "scalar", "NaN input did not produce NaN scalar output")) + if row.batch is not None and not is_nan(row.batch): + g.violations.append((row.lineno, "batch", "NaN input did not produce NaN batch output")) + continue + + if inf_input and row.method in ("pdf", "logpdf", "cdf"): + # +-inf inputs have universal limits for every distribution on + # (a subset of) the real line: pdf -> 0, logpdf -> -inf, + # cdf -> 1 at +inf / 0 at -inf. Evaluating the raw reference + # formulas AT the limit instead produces inf-inf / 0*inf NaN + # artifacts, which is an oracle artifact, not a finding. Von + # Mises is periodic -- no limit exists, the library's NaN is + # the right answer, so check NaN agreement like NaN-input rows. + if row.dist == "von_mises": + # Periodic distribution: no mathematical limit exists at + # +-inf, so no reference value can be asserted. The library + # uses a deliberate saturation convention (pdf -> 0, + # logpdf -> -inf, cdf -> 0/1), which the characterization + # doc records as a convention; the only checkable contract + # here is that scalar and batch agree with each other. + if row.batch is not None and not ( + (is_nan(row.scalar) and is_nan(row.batch)) + or row.scalar == row.batch + ): + g.violations.append((row.lineno, "batch", + f"scalar/batch disagree at +-inf input " + f"({row.scalar} vs {row.batch})")) + continue + if row.method == "pdf": + ref = mp.mpf(0) + elif row.method == "logpdf": + ref = mp.mpf("-inf") + else: + ref = mp.mpf(1) if row.x > 0 else mp.mpf(0) + else: + try: + _t0 = time.perf_counter() + ref = fn(row.p1, row.p2, row.x) + _dt = time.perf_counter() - _t0 + if _dt > 1.0: + print(f" [slow-row] line {row.lineno} {row.dist}/{row.method} " + f"p1={row.p1!r} p2={row.p2!r} x={row.x!r} took {_dt:.1f}s", + file=sys.stderr, flush=True) + except (ValueError, ZeroDivisionError, OverflowError) as exc: + if inf_input: + skipped.append((row.lineno, f"oracle skipped inf-input row ({exc})")) + continue + raise + + if not mp.isfinite(ref): + # Reference itself is +-inf (e.g. quantile at p=1 for unbounded + # support): contract-check finiteness/sign agreement rather than + # a relative error, which is not meaningful against an infinite + # reference. + if mp.isnan(ref): + # A NaN reference on a non-inf input row is an oracle bug, + # not a library finding; surface it as such. + g.violations.append( + (row.lineno, "oracle", f"oracle produced NaN reference (x={row.x!r})") + ) + continue + scalar_ok = math.isinf(row.scalar) and (row.scalar > 0) == (ref > 0) + if not scalar_ok: + g.violations.append( + (row.lineno, "scalar", f"reference is {ref}, scalar_bits decoded to {row.scalar}") + ) + if row.batch is not None: + batch_ok = math.isinf(row.batch) and (row.batch > 0) == (ref > 0) + if not batch_ok: + g.violations.append( + (row.lineno, "batch", f"reference is {ref}, batch_bits decoded to {row.batch}") + ) + continue + + is_cdf_tail = row.method == "cdf" and ref < mp.mpf("1e-3") + if not math.isfinite(row.scalar): + g.violations.append( + (row.lineno, "scalar", f"reference is finite ({mp.nstr(ref, 6)}), scalar_bits decoded to {row.scalar}") + ) + continue + g.scalar.observe(mp.mpf(row.scalar), ref, row.x, is_cdf_tail) + + if row.batch is not None and not math.isfinite(row.batch): + g.violations.append( + (row.lineno, "batch", f"reference is finite ({mp.nstr(ref, 6)}), batch_bits decoded to {row.batch}") + ) + elif row.batch is not None: + g.n_batch_rows += 1 + g.batch.observe(mp.mpf(row.batch), ref, row.x, is_cdf_tail) + bs_abs = abs(mp.mpf(row.batch) - mp.mpf(row.scalar)) + bs_rel = relerr(mp.mpf(row.batch), mp.mpf(row.scalar)) + g.batch_vs_scalar_max_abs = max(g.batch_vs_scalar_max_abs, bs_abs) + g.batch_vs_scalar_max_rel = max(g.batch_vs_scalar_max_rel, bs_rel) + + return groups, skipped + + +# --------------------------------------------------------------------------- +# Reporting +# --------------------------------------------------------------------------- + + +def fmt(x) -> str: + if isinstance(x, mp.mpf): + if x == 0: + return "0" + return mp.nstr(x, 4, strip_zeros=True) + return str(x) + + +def print_stdout_summary(groups: dict, skipped: list, violations_total: int) -> None: + print(f"{'dist':<18} {'method':<9} {'src':<6} {'n':>5} {'max_abs':>12} {'max_rel':>12} " + f"{'p99_rel':>12} {'law_frac':>10} {'worst_x':>14}") + for (dist, method), g in sorted(groups.items()): + for src_name, s in (("scalar", g.scalar), ("batch", g.batch)): + if s.n == 0: + continue + law_col = fmt(s.max_law_frac) if s.has_law_rows else "-" + print( + f"{dist:<18} {method:<9} {src_name:<6} {s.n:>5} {fmt(s.max_abs):>12} " + f"{fmt(s.max_rel):>12} {fmt(s.p99_rel()):>12} {law_col:>10} {s.worst_x:>14.6g}" + ) + if g.n_batch_rows: + print( + f"{'':<18} {'':<9} {'b_vs_s':<6} {g.n_batch_rows:>5} " + f"{fmt(g.batch_vs_scalar_max_abs):>12} {fmt(g.batch_vs_scalar_max_rel):>12}" + ) + if g.violations: + for lineno, src, msg in g.violations: + print(f" VIOLATION line {lineno} ({src}): {msg}") + + if skipped: + print(f"\n{len(skipped)} row(s) skipped:") + for lineno, msg in skipped[:20]: + print(f" line {lineno}: {msg}") + if len(skipped) > 20: + print(f" ... and {len(skipped) - 20} more") + + print(f"\ntotal contract violations: {violations_total}") + + +def render_markdown(groups: dict) -> str: + lines = [] + dists = sorted({d for d, _ in groups.keys()}) + for dist in dists: + lines.append(f"### {dist}\n") + lines.append("| method | source | max_abs | max_rel | p99_rel | law_frac(cdf) | batch_vs_scalar | worst_x |") + lines.append("|---|---|---|---|---|---|---|---|") + methods = sorted(m for d, m in groups.keys() if d == dist) + for method in methods: + g = groups[(dist, method)] + for src_name, s in (("scalar", g.scalar), ("batch", g.batch)): + if s.n == 0: + continue + law_col = fmt(s.max_law_frac) if s.has_law_rows else "-" + bvs = ( + f"abs={fmt(g.batch_vs_scalar_max_abs)}, rel={fmt(g.batch_vs_scalar_max_rel)}" + if src_name == "batch" and g.n_batch_rows + else "-" + ) + lines.append( + f"| {method} | {src_name} | {fmt(s.max_abs)} | {fmt(s.max_rel)} | " + f"{fmt(s.p99_rel())} | {law_col} | {bvs} | {s.worst_x:.6g} |" + ) + if g.violations: + lines.append( + f"| {method} | *(contract)* | {len(g.violations)} violation(s) -- see appendix | | | | | |" + ) + lines.append("") + # Durable appendix: every contract violation, so the checked-in doc + # stands alone without the stdout log. + all_v = [ + (dist, method, lineno, who, msg) + for (dist, method), g in sorted(groups.items()) + for (lineno, who, msg) in g.violations + ] + if all_v: + lines.append("### Contract findings (appendix)" + "\n") + lines.append( + f"{len(all_v)} contract violations across the sweep. `csv_line` " + "indexes the sweep CSV this report was generated from (see the " + "commit/isa banner in the regeneration log)." + ) + lines.append("") + lines.append("| dist | method | source | csv_line | finding |") + lines.append("|---|---|---|---|---|") + for dist, method, lineno, who, msg in all_v: + lines.append(f"| {dist} | {method} | {who} | {lineno} | {msg} |") + lines.append("") + return "\n".join(lines) + + +GEN_BEGIN = "" +GEN_END = "" + + +def rewrite_doc(doc_path: str, generated_body: str) -> None: + with open(doc_path, "r", encoding="utf-8") as f: + text = f.read() + if GEN_BEGIN not in text or GEN_END not in text: + raise ValueError(f"{doc_path}: missing {GEN_BEGIN}/{GEN_END} markers") + pre, rest = text.split(GEN_BEGIN, 1) + _, post = rest.split(GEN_END, 1) + new_text = pre + GEN_BEGIN + "\n\n" + generated_body + "\n" + GEN_END + post + with open(doc_path, "w", encoding="utf-8") as f: + f.write(new_text) + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("csv", help="path to the accuracy_sweep CSV output") + parser.add_argument( + "--out", + default=None, + help="path to docs/ACCURACY_CHARACTERIZATION.md (default: alongside this script's repo root)", + ) + args = parser.parse_args(argv) + + repo_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + doc_path = args.out or os.path.join(repo_root, "docs", "ACCURACY_CHARACTERIZATION.md") + + print("Running oracle self-checks (mp.dps=50)...") + try: + checks = run_self_checks() + except AssertionError as exc: + print(f"ORACLE SELF-CHECK FAILED: {exc}", file=sys.stderr) + return 1 + print(f" {len(checks)} self-checks passed.\n") + + rows, meta = parse_csv(args.csv) + if "banner" in meta: + print(f"Input banner: {meta['banner']}") + print(f"Parsed {len(rows)} rows from {args.csv}\n") + + groups, skipped = compare(rows) + violations_total = sum(len(g.violations) for g in groups.values()) + + print_stdout_summary(groups, skipped, violations_total) + + body = render_markdown(groups) + rewrite_doc(doc_path, body) + print(f"\nWrote generated block to {doc_path}") + + return 0 + + +if __name__ == "__main__": + sys.exit(main())