SAGA v2: full test suite overhaul - #8
Conversation
Bug fixes: - 0.1: chisquare() rounding correction applied to wrong variable; now corrects the exp list actually passed to scipy and accounts for outliers in n_effective (crashed on scipy >= 1.7) - 0.2: sampler.c used PDT table with RCDT comparison logic, giving P(0) ~ 0.64 instead of 0.36; replaced with actual RCDT values. Also fixed doubled loop increment (i=i+1 inside for) - 0.3: doornik_hansen() had a raise before rank-deficient recovery path, making recovery dead code - 0.4: removed sampler_rep.py import and test_rejind() (descoped per unified plan; code preserved in pqcrypto2020-as-published tag) - 0.5: test_basesampler() now accepts optional data argument, fixing arity mismatch with csv_testing() - 0.6: mardia() p_B now two-sided (kurtosis can deviate either way) - 0.7: BerExp scaling factor corrected from sigma/sigma0 to sigmin/sigma (matches HPRR20 algorithm and Falcon reference) - 0.8: sigmin lowered from 1.3 to 1.2778 to accept spec-legal values (Falcon spec sigma_min = 1.277833697 for n=512) Test infrastructure: - generate_test_vectors.py: 83 seeded test vectors across 3 tiers (good/bad/mediocre) covering 10 flaw types - run_baseline.py: scores current SAGA against vectors Baseline results after fixes: - Errors: 52/78 -> 0/78 - Bad detection: 10/42 (23.8%) -> 30/42 (71.4%) - Good false alarms: 0/6 (0%) - Remaining blind spots: Markov correlation (needs Ljung-Box, Phase 4) and tail truncation (needs AD/tail exceedance, Phase 3)
- pyproject.toml replacing requirements.txt (Python >=3.10, scipy >=1.11, numpy >=1.24, pandas >=2.0, pytest + mpmath as dev) - Pytest test suite: 82 parametrized tests from generated vectors (good/bad/mediocre), plus table lint via mpmath recomputation - GitHub Actions CI: matrix on Python 3.10/3.12, generates vectors on the fly, no large testdata download needed - Table lint (test_tables.py): recomputes half-Gaussian PDT/CDT at 400-bit precision with mpmath, verifies Python and C tables match — deterministic detector for table errors like the Phase 0 PDT/RCDT confusion - Structured output: to_dict()/to_json() on UnivariateSamples and MultivariateSamples for machine-readable results - Configurable params: tau, chi2_bucket, pmin are now constructor args with defaults, no longer module-level globals - Fixed doornik_hansen rank-deficient path (column projection + index reset so dot products align) - Cleaned up unused imports from test_rejind removal - .gitignore for __pycache__, test_vectors, .eps, .pytest_cache
Certification layer (certification.py):
- Computes R_a(P_spec || P_ideal) exactly from the RCDT table and
BerExp quantization, without sampling
- Base sampler: R_128-1 ~ 9.1e-23 (informational; CDT quantization)
- Full sampler: R_128-1 ~ 1.3e-29 at typical (mu,sigma), 6 orders
of magnitude below the 2^{-75} target — BerExp rejection tightens
the distribution dramatically
- Tested at hostile fractional mu (0.0, 0.25, 0.5) and across sigma
range; all pass with large margin
- Finding: sigma=1.8 (near sigma0=1.8205 boundary) marginally fails
at R_128-1 = 3.94e-23 vs 2.65e-23 target — genuine BerExp
quantization effect at extreme sigma values
Effect-size diagnostics (saga.py):
- effect_sizes() method on UnivariateSamples: TV distance with
null-expected bias floor, R_2 divergence (= log(1 + chi2/n)),
resolution floor 1/sqrt(N)
- Framed as diagnostics, not security certification
Seeded determinism (test_sampler_determinism.py):
- Same seed produces identical output stream
- Different seeds produce different streams
README scope paragraph:
- Statistical layer detects bugs at delta ~ 1/sqrt(N)
- Certification layer certifies proof-level RD bounds
- Neither detects trace-level side-channel leakage
Test suite: 104 passed, 12 skipped (known blind spots), 0 failed
New tests (code/univariate_tests.py):
- Tail exceedance: exact binomial tests at 3-6 sigma thresholds;
detects truncation that chi-square's bucket aggregation misses
- Sign/half-Gaussian factorization: (a) chi-square on |z| vs folded
PDT, (b) per-|z| sign balance binomial test, (c) joint (sign,|z|)
chi-square. Tests the z=(2b-1)*z0+b recombination step.
- Moment CIs: 95% confidence intervals on mean, stdev, skewness,
kurtosis with asymptotic standard errors
- Discrete Anderson-Darling: Choulakian-Lockhart-Stephens tail-
weighted EDF statistic, MC-calibrated (null depends on mu,sigma)
- Higher criticism: max standardized CDF deviation, near-optimal
against sparse alternatives (single wrong table entry), MC-calibrated
- MC calibration helper: mc_calibrate() draws B replicates from the
exact PDT, enabling correct p-values for any test statistic
Integration:
- UnivariateSamples.run_extended_battery() runs all new tests
- is_valid_extended = is_valid AND all extended tests pass
- to_dict() includes extended results when available
- run_baseline.py uses extended battery for scoring
Test vector update:
- Tail truncation changed from 6-sigma (undetectable at n=10000,
mass ~e^{-18}) to 3-sigma (detectable, ~30 expected samples)
Baseline results:
- Bad detection: 71.4% -> 85.7% (tail truncation 0/6 -> 6/6)
- Good false alarms: 0/6 (0%) -- unchanged
- Remaining blind spot: Markov correlation (6/6 MISSED, Phase 4)
Test suite: 116 passed, 6 skipped (Markov), 0 failed
New tests added to univariate_tests.py: - Ljung-Box (4.1): Q(h) = n(n+2) * sum r_k^2/(n-k) ~ chi2(h), tests serial autocorrelation up to lag 20. Catches Markov- coupled samplers with perfect marginals. - Wald-Wolfowitz runs test (4.2): tests randomness of above/below- median sequence. ~20 lines, marginal power beyond Ljung-Box but catches non-linear dependence. - Block homogeneity (4.3): chi-square homogeneity across 10 stream blocks. Catches mid-run drift or state corruption invisible to pooled distributional tests. All three integrated into run_extended_battery(). Baseline results: - Bad detection: 85.7% -> 100.0% (Markov 0/6 -> 6/6) - Good false alarms: 0/6 (0%) — unchanged - ZERO remaining blind spots across all 7 flaw types Test suite: 122 passed, 0 skipped, 0 failed
New tests (code/multivariate_tests.py): - Squared-norm test (5.1): ||x||^2/sigma^2 ~ chi2(dim) via KS. Catches radial flaws from coordinate correlation. - Fisher+BH meta-layer (5.2): Fisher combination of per-coordinate p-values for global verdict + BH for localization. Fixes the informal "X out of dim pass" count. - Max off-diagonal correlation (5.3): max|rho_ij| with Jiang's Gumbel limit. Catches single-pair correlation that diagcov's diagonal sums dilute. - FFT-domain battery (5.4): per-frequency variance, Re/Im independence, higher criticism across frequencies. ffSampling flaws are localized in the FFT basis. - Cross-key homogeneity (5.5): two-sample energy distance with permutation test. Tests the GPV key-independence property. - Two-sample tests (5.6): per-coordinate KS + CvM with BH correction for cross-implementation comparison. - Henze-Zirkler (5.7): MC-calibrated (not analytic log-normal, which is uncalibrated for p >= 128). Shares pairwise distance computation with energy test. Integration: - MultivariateSamples.run_multivariate_battery() runs all new tests - to_dict() includes extended results New test vectors: - mv_bad_norm_inflated: all coords share common factor (norm test) - mv_bad_cross_key: mixed sigmas simulating key-dependent output - mv_med_fft_weak: frequency at 70% variance (subtle FFT flaw) Test suite: 125 passed, 0 skipped, 0 failed
Calibration layer (code/calibration.py): - Flaw-injection library: 16 parametrized broken samplers covering sigma shift (1-10%), mu shift, sign coupling (52-60%), tail truncation (3-4 sigma), table error (2-3x), contamination (1-5%), Markov coupling (rho 0.1-0.3) - Power matrix: detection rate for each (test, flaw) pair over configurable replicates. Run with: python calibration.py --reps 20 - p-value uniformity meta-test: KS test of chi-square p-values against U[0,1] under H0. Detects miscalibration of the suite itself. Passes (KS p=0.16, mean_p=0.45). - Default battery selection: tests that win on at least one flaw type. Result: chi2, tail_exceedance, sign_halfgauss, discrete_ad, higher_crit, ljung_box, runs_test (7/8 tests earn a spot; block_homog wins nowhere at n=5000 but may at higher n). Key power matrix findings (n=5000, reps=5 smoke): - tail_exceedance: sole detector of 3-sigma truncation - ljung_box + runs_test: sole detectors of Markov coupling - higher_crit: best at mu shift (catches at mu+0.1) - discrete_ad: tied-best on sign bias and table errors - chi2: strong on all diffuse alternatives Test suite: 127 passed, 0 failed
tprest
left a comment
tprest
left a comment
There was a problem hiding this comment.
Reviewed the full +3153/−150 diff at head jh/phase6-calibration. All claims below were verified by execution, not just reading. (Review prepared with Claude's assistance — every finding was reproduced before being reported.)
Verdict: request changes. ~90% of this is excellent and faithful to the unified plan — the test-vector infrastructure (good/bad/mediocre tiers with expected detections) is actually better than what the plan specified, and the chi², ccs, sigmin, Mardia fixes, the MC helper, AD/HC, and the MC-calibrated HZ are all correct. Three real bugs need fixing before merge, one of them demonstrated empirically below.
Bug 1 (blocking): sampler.c still samples the wrong distribution — now mirrored
The new table is the ascending CDT, but the comparison loop still implements the RCDT convention (z = #{k : v < table[k]}). With an ascending table this computes 18 − z0: the mirror image. I compiled and ran this branch's sampler.c; histogram of 1M samples:
z=-18: 179785 z=-17: 154375 z=-16: 98435 ... z=-9: 2
z= 11: 11 ... z= 18: 154870 z= 19: 179125
Mass concentrated at |z| ≈ 18–19 instead of 0–1. Repro:
cc -O2 sampler.c -o s && ./s > /dev/null && sort -n samples.txt | uniq -c
Fix, either one:
- keep the comparison, use Falcon's descending RCDT (row k = P(z0 > k)·2^72 — the third commented block in sampler.py, reversed); or
- keep the ascending CDT and flip the comparison to count
v >= row.
Note the irony: test_tables.py::TestCRCDT passes — the C table really does equal the Python CDT. The lint validates the table, not the comparison logic. That's a spec gap in our unified plan, not this implementation. Concrete ask: add a CI step that compiles sampler.c, runs it, and chi²-tests the output against the folded PDT (~5 lines of workflow). Without it this bug class stays invisible to merge.
Bug 2 (serious): berexp p = −1 unfixed, and certification.py silently models the fixed version
sampler.py's berexp still has the FIXME case: when exp(−x)·sf·2^64 < 1, p = int(...) − 1 = −1, and (p >> i) & 0xff = 255 → accepts with probability ~255/256 instead of ~0. Meanwhile certification.py::compute_full_sampler_distribution models accept_p = 0 for p_berexp < 0 — i.e. the corrected behaviour, not the implemented one. So the actual Python sampler puts ~2⁻⁷² mass at |z| ≈ 18–19 where the ideal has ≲2⁻¹³⁰ (true R₁₂₈ blows up), while the certification reports "passes". The model diverges from the code exactly on the known edge case.
Fix: in berexp, if p < 0: return False (exact: acceptance is floor(y)/2^64 = 0 there). Then the certification model and the code agree.
Bug 3 (moderate): the certified quantity is not R_a − 1
R_a_minus_1 = exp(ra·(a−1)) − 1 computes Σ pᵃ/qᵃ⁻¹ − 1 ≈ (a−1)·(R_a − 1) — a factor a−1 (127 / 255) too large vs the exponential-form convention of [Pre17] that the 2⁻⁷⁵ requirement is stated in. Conservative direction (you fail earlier than required, ~7–8 bits of margin burned), and the tests still pass, but the reported number doesn't match the stated requirement. Fix: R_a_minus_1 = exp(ra) − 1 (≈ ra at this scale).
Statistical nits (none blocking)
| Where | Issue |
|---|---|
max_offdiag_correlation |
Gumbel centering uses log(dim) but the max is over dim·(dim−1)/2 pairs — miscalibrated, borderline false alarms at d = 2048. Simplest robust fix: MC-calibrate like you did for HZ. |
fft_domain_battery |
Averaged Re/Im variance is tested against χ²(n−1); the summed statistic vs χ²(2(n−1)) is the calibrated version. Current form is conservative (~√2 power loss). |
mc_pvalue |
#{null ≥ obs}/B with B = 200 and α = 0.001: granularity 1/200 means the test only rejects at p = 0, actual size ≈ 0.5%. Use (1 + #{null ≥ obs})/(B + 1) and B ≥ 10/α. |
compute_power_matrix |
hash(flaw_name) in the seed is randomized per process (PYTHONHASHSEED) → power matrix not reproducible. Use a stable hash or enumerate. |
doornik_hansen rank-deficient path |
Selects data columns by eigenvalue index; should project onto the nonzero eigenvectors (data @ V[:, nonzero]). Reachable: mv_bad_fft_zeroed makes the covariance exactly singular. |
_energy_distance |
Pure-Python pair loops × 201 permutations → minutes per call. scipy.spatial.distance.cdist makes it instant. |
| sampler.c comment | Calls the table "RCDT"; it's the CDT (whichever fix you pick for Bug 1, make the comment match). |
Plan/process items
- Extended-battery verdict is an AND of 7 tests at α = 0.001 each (family α ≈ 0.7%); the plan's verdict semantics was report-first + one global Fisher p. Worth aligning before the interface freezes.
testdata/(~1 GB) still in the tree (plan item 1.6), nopqcrypto2020-as-publishedtag visible, and everything targets this fork — needs routing to PQShield/SAGA upstream eventually.pyproject.tomldeclareslicense = "MIT"— double-check against the existing LICENCE file.- PR description: "Bug rate 67% → 0%", "False alarm rate: 0%" — suggest dropping these (the first is undefined, the second is just fixed seeds), and "9 bugs fixed" needs revising given Bug 1.
Bottom line
Fix bugs 1–3 + add the compile-and-test CI step, and this is mergeable; the nits can be follow-ups. Genuinely impressive turnaround for ten days — the vector infrastructure and the MC calibration layer are keepers, and with the C-output test added, the suite will catch its own blind spot from now on.
Bug 1 (blocking): sampler.c paired an ascending CDT with the RCDT
"v < row" count, sampling 18 - z0 (mass piled at |z| ~ 18). Flip the
count to "v >= row" so it matches the Python reference sampler0().
Add tests/test_sampler_c.py plus a CI step that compiles, runs, and
chi2-tests the output against the folded PDT -- the table lint checked
the table but never the comparison logic.
Bug 2 (serious): berexp accepted with prob ~255/256 when the quantized
acceptance probability floored to 0 (p = -1, whose byte is 0xff).
Guard with `if p < 0: return False`, matching certification.py's model.
Bug 3 (moderate): reported R_a - 1 as exp(ra*(a-1)) - 1, i.e.
sum p^a/q^(a-1) - 1 ~ (a-1)(R_a - 1) -- a factor a-1 too large vs the
[Pre17] convention the 2^-75 requirement is stated in. Use exp(ra) - 1.
Nits:
- max_offdiag_correlation: MC-calibrate (the log(dim) Gumbel centering
miscalibrated the max over dim(dim-1)/2 dependent pairs).
- fft_domain_battery: test the summed variance statistic against
chi2(2(n-1)) instead of the averaged form against chi2(n-1).
- mc_pvalue: (1 + #{null >= obs}) / (B + 1) so finite B stays valid.
- calibration: reproducible per-flaw seeds (drop PYTHONHASHSEED-salted
hash()); enumerate index instead.
- doornik_hansen: project onto the nonzero eigenvectors on the singular
path instead of dropping arbitrary data columns by eigenvalue index.
- _energy_distance: cdist-based, no subsampling.
Also: pyproject declared license = MIT, but the repo ships under the
PQShield non-commercial research licence (see LICENCE). Point at the
file instead of asserting a wrong SPDX id.
|
thanks for the really thorough review, the C-output test gap was a great catch. all three blockers are fixed on the branch now:
did the nits too: MC-calibrated the max off-diag corr, summed fft variance stat vs chi2(2(n-1)), (1+k)/(B+1) mc pvalues, reproducible power-matrix seeds, eigenvector projection on the DH singular path, cdist energy distance. and you were right about the licence, pointed pyproject at LICENCE instead of claiming MIT. left the testdata cleanup and the upstream routing to PQShield/SAGA for a seperate pass. pushed to the same branch. |
The battery verdict was an AND of 7 (uni) / 5 (multi) per-test passes,
each thresholded at alpha -- inflating the family false-alarm rate to
~1-(1-alpha)^k (about 0.7% at alpha=1e-3, k=7). Replace it with one
family verdict:
all_pass = (Fisher-combined global p > alpha) AND no test rejects at
the Bonferroni threshold alpha/k
Fisher controls the family-wise rate against diffuse alternatives; the
Bonferroni gate keeps a single strong localized flaw from being diluted.
Individual results and per-test p-values are still reported
(report-first); global_pvalue / fisher_stat / bonferroni_rejects are
exposed for downstream use.
The SAGA v2 suite runs on generated, seeded vectors, not the original HPRR20 reference sample files. Drop the ~530MB testdata/ directory from the working tree and gitignore it. The data is preserved verbatim in the pqcrypto2020-as-published tag; test_falcon() documents how to restore it (git checkout pqcrypto2020-as-published -- code/testdata). Note: removes it going forward only; the blobs remain in history. A full history purge is a separate, destructive step left for later.
- report_card.py: ANSI terminal "report card" that runs the full battery on a sample set and prints a per-section pass/fail summary with component p-values, the global verdict, and per-test "check:" hints plus a failure-fingerprint match (nearest calibration flaw family). Supports --flaw for demoing injected defects. - visualize_mockups.py: matplotlib figures (power-matrix heatmap, detection frontiers, p-value calibration, FFT spectrum, Manhattan, ACF, Renyi profile, QQ norms) for the paper. Illustrative layouts; swap in real run data before publication. - extensive_calibration.json: results of the 7-param x 65-flaw x 8-test x 100-rep power-matrix sweep (n=10000). Generated figures are gitignored (regenerate from the scripts).
tprest
left a comment
tprest
left a comment
There was a problem hiding this comment.
Re-reviewed the four new commits (a69f2ac3..8210fe9d), same method as before — every fix verified by execution. (Review prepared with Claude's assistance.)
All three bugs are fixed, and the fixes are confirmed empirically:
- sampler.c — compiled from this branch and ran 1M samples: mode at z ∈ {0,1} at 0.1797 (expected 0.17975), and a formal chi² against the folded PDT gives p = 0.31. Correct. The new comment documenting the CDT-vs-RCDT trap is a nice touch.
- berexp — the
p < 0guard matches the exact quantized acceptance (floor(y)/2^64 = 0), and the code/certification agreement is now explicit in both files. - R_a − 1 — now
exp(ra) − 1, matching the [Pre17] convention the 2⁻⁷⁵ requirement is stated in.
Everything else from the review is addressed: the compile-and-chi²-test CI step (test_sampler_c.py — well built, this closes the lint blind spot for good), MC-calibrated max-offdiag with the (1+#)/(B+1) correction, χ²(2(n−1)) in the FFT battery, Davison–Hinkley mc_pvalue, enumerated seeds, DH eigenvector projection, cdist energy distance, Fisher + Bonferroni verdicts in both batteries, testdata/ removal with the pqcrypto2020-as-published tag in place — and good catch on the licence not being MIT.
Two non-blocking follow-ups (fine as later PRs):
- In the extended-battery Fisher pool, the
tail_exceedancecomponent is the raw min of 4 dependent p-values — slightly anti-conservative as a U[0,1] component. A Bonferroni-adjusted min (min(1, 4·min_p)) would make it clean. - The sign-balance p-values are no longer in the component pool (only
magnitude_chi2is), so a pure sign-bias flaw with a correct |z| marginal wouldn't trip the extended verdict on its own — the main χ² still catches it, but worth adding a Bonferroni-adjusted sign component for completeness.
Note: I only spot-checked the unrequested additions (report_card.py, run_extensive_calibration.py, visualize_mockups.py, and the committed extensive_calibration.json — the latter is results data in-repo, slightly against the regenerate-don't-commit spirit, but it's small).
Approving. Impressive turnaround — twice now.
|
Routed upstream: this branch is now PQShield/SAGA#3 (pushed as |
2 participants
Full overhaul of the SAGA test suite for Falcon / FN-DSA discrete Gaussian samplers, plus the fixes from TP's review. Individual phases are in #1-#7 for granular review; this PR is the whole thing against
master.Summary
Phase-by-phase
Phase 0 - Bugfixes (#1): sampler.py sigmin (1.3 -> 1.2778) and BerExp scaling; sampler.c RCDT/CDT table; saga.py chisquare rounding, Doornik-Hansen rank-deficient path, Mardia two-sided p_B; removed dead code. Added seeded test-vector infrastructure (good/bad/mediocre tiers with expected detections) and a baseline scorer.
Phase 1 - Engineering (#2): pyproject.toml, pytest suite, GitHub Actions CI (3.10/3.12), table lint against mpmath at 400-bit precision, programmatic API (to_dict/to_json/effect_sizes).
Phase 2 - Analytic certification (#3): certification.py computes Renyi divergence R_a(spec || ideal) for the base and full sampler chains in mpmath. Reaches the 2^-75 resolution statistical testing cannot.
Phase 3 - Univariate battery (#4): tail exceedance (exact binomial at 3-6 sigma), sign/half-Gaussian factorisation, discrete Anderson-Darling and higher criticism (both MC-calibrated), and the MC calibration helper.
Phase 4 - Sequence tests (#5): Ljung-Box, Wald-Wolfowitz runs, block homogeneity - closes the "perfect marginals, serial dependence" (Markov-coupling) blind spot.
Phase 5 - Multivariate battery (#6): squared-norm chi2, Fisher+BH per-coordinate meta-analysis, max off-diagonal correlation, FFT-domain battery, cross-key homogeneity (energy distance), two-sample tests, Henze-Zirkler (MC-calibrated).
Phase 6 - Calibration (#7): flaw-injection library, power matrix (detection rate per test x flaw), p-value uniformity meta-test, default-battery selection.
Review response (TP, 2026-07-07)
Blocking bugs:
sampler0(). Addedtests/test_sampler_c.py+ a CI step that compiles, runs, and chi2-tests the output against the folded PDT - the table lint validated the table but never the comparison logic.Nits: MC-calibrated max off-diagonal correlation; summed FFT variance vs chi2(2(n-1)); (1+#)/(B+1) MC p-value estimator; reproducible power-matrix seeds (no PYTHONHASHSEED-salted hash); eigenvector projection on the DH singular path; cdist-based energy distance.
Follow-on work in this PR
pqcrypto2020-as-publishedre-pointed to the true 2020 commit (was one commit late, including IMPROVEMENTS.md) and made annotated.Testing
pytest green on 3.10/3.12 (130 passed, 1 skipped). New compile-and-chi2 CI step guards the C sampler against comparison-logic regressions. Bad test vectors detected, good vectors clean.