From d7fa2fd50c9d39797ae3334997d66e401d37ddfb Mon Sep 17 00:00:00 2001 From: Ruge Lin Date: Mon, 7 Sep 2026 10:25:54 +0800 Subject: [PATCH] Harden binomial tails and validate every committed result Add stable probability/log-probability APIs, explicit positive-underflow handling, comprehensive result-field and checksum validation, deliberate corruption tests, and Python 3.10-3.13 CI. Preserve the August scientific correction, PDF, results, and six historical scripts. Version 1.1.1. --- .github/workflows/tests.yml | 22 ++- CHANGELOG.md | 9 + CHECKSUMS.sha256 | 29 ++-- CITATION.cff | 4 +- README.md | 9 + REPRODUCIBILITY.md | 78 ++++++++- VALIDATION_RECORD.md | 13 ++ pyproject.toml | 2 +- scripts/generate_results.py | 226 +++---------------------- scripts/validate_committed_results.py | 58 ++----- scripts/verify_checksums.py | 8 + src/dcp_challenge/__init__.py | 2 +- src/dcp_challenge/integrity.py | 63 +++++++ src/dcp_challenge/reference_results.py | 219 ++++++++++++++++++++++++ src/dcp_challenge/result_validation.py | 112 ++++++++++++ src/dcp_challenge/statistics.py | 115 ++++++++++--- tests/test_integrity.py | 35 ++++ tests/test_result_validation.py | 153 +++++++++++++++++ tests/test_stable_statistics.py | 75 ++++++++ 19 files changed, 933 insertions(+), 299 deletions(-) create mode 100644 scripts/verify_checksums.py create mode 100644 src/dcp_challenge/integrity.py create mode 100644 src/dcp_challenge/reference_results.py create mode 100644 src/dcp_challenge/result_validation.py create mode 100644 tests/test_integrity.py create mode 100644 tests/test_result_validation.py create mode 100644 tests/test_stable_statistics.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 9b8d7a3..2c2c600 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -3,23 +3,33 @@ name: corrected-analysis-tests on: push: pull_request: + workflow_dispatch: + +permissions: + contents: read jobs: test: runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.10", "3.11", "3.12", "3.13"] steps: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 with: - python-version: "3.11" + python-version: ${{ matrix.python-version }} cache: pip - name: Install corrected package - run: python -m pip install -e .[test] - - name: Run exact and numerical regression tests - run: pytest -q + run: python -m pip install -e '.[test]' + - name: Verify complete checksum manifest + run: python scripts/verify_checksums.py + - name: Run exact, numerical, and corruption regression tests + run: python -m pytest -q - name: Verify historical root scripts run: python scripts/verify_historical_files.py - - name: Validate committed scientific results + - name: Validate every committed scientific result field run: python scripts/validate_committed_results.py - name: Reproduce decisive exact counterexample run: python examples/reproduce_counterexample.py @@ -29,3 +39,5 @@ jobs: run: python examples/reproduce_figure5.py - name: Check minimal four-qubit circuit run: python examples/simulate_minimal_four_qubit_circuit.py + - name: Confirm validation did not modify tracked files + run: git diff --exit-code diff --git a/CHANGELOG.md b/CHANGELOG.md index 460bb48..80ecf4f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Changelog +## 1.1.1-author-correction - 2026-09-07 + +- Replaced overflow-prone binomial evaluation with scaled recurrence and log-tail APIs; positive underflow raises explicitly. +- Expanded result validation to all five result files and 458 fields, including schemas, archived Monte Carlo metadata, and Wilson arithmetic. +- Added tests corrupting every result field, checksum integrity tests, and exact rational statistical references. +- Added complete manifest verification and Python 3.10-3.13 CI coverage. +- Kept the August correction PDF, mathematics, results, historical scripts, and original release unchanged. +- Fast checks explicitly distinguish archived Monte Carlo records from fresh simulation. + ## 1.1.0-author-correction - 2026-08-17 ### Referee-level revision diff --git a/CHECKSUMS.sha256 b/CHECKSUMS.sha256 index 1e5bef6..abd73a4 100644 --- a/CHECKSUMS.sha256 +++ b/CHECKSUMS.sha256 @@ -1,14 +1,14 @@ -f6d1edcdabe82a8c9a6b8676470dca4755bd72abc667f9a542c0996794c8cef0 ./.github/workflows/tests.yml +3f8c91e32f752f38c51834eeb0a38abfd791921d0be7e247ac840bb48a93ee26 ./.github/workflows/tests.yml f0e43839ec4fe3f086a0997f1aa238429543c8784776b05e019a719629ee1d90 ./.gitignore 03353aa4d79773163b3c22bacc52f368149ef0b37abf8858a904169ebc7455c0 ./AI_ASSISTED_REANALYSIS.md -73e68d44953e332f88f65ac908a872fad3fe8e39473c5c6fbbcc470af25e31d7 ./CHANGELOG.md -7aaede28f7b53fdfa7f997bdb208fadb1c2339f3908812edbe9553003446d1d0 ./CITATION.cff +e9e71e870915420013b31dd5f23d5c1c41ba87cbdca3eb4f752a50a91c5cb747 ./CHANGELOG.md +ca826a2121e718472fff0f5d148f246467a519a7e4aa1d52bebe0ca7845a2cc6 ./CITATION.cff 0c4959711013692777b7e98a35a0c8aa9e6fb85dc2583d617b303ed0cf3b0bfb ./IBM.py aec03ca89a1ad1db2feba926dc30a6414e297a245987b3a77fdfa0af5efa24dc ./ORIGINAL_2022_CODE.md -8b2ce91ce9eb49e17cb281a9231aadbea1194c7cf6b71429631be4f156b6b430 ./README.md -0fa3c0c2a8785e07f02f6898180e58e7e3abe2b1a1bbffed2791c5c4f3a2c9da ./REPRODUCIBILITY.md +fcf40b8f974ba65a7ec932643886a951d0860c10656934fec757ffb03ded0261 ./README.md +7ad37e80a44290ad50c57c1aaa8479e52da51ccb952eea20cda8ee53ad1eb065 ./REPRODUCIBILITY.md 52a7415b170f3ce96863c792633132fdfc2cf69fbd78bddd7eff38c5553b0b7d ./SCIENTIFIC_STATUS.md -e949ae83ccc8224a136145832fed7d91a161a8e2cb13853e69225812498a5baf ./VALIDATION_RECORD.md +c432acd797b3f422af4a1101c90d4b6c8b5cce94a05492cc8e14fc142a1a7bd8 ./VALIDATION_RECORD.md 5a9b4e0129e44dbb52d368f9939a12cdedbd58e07d8e7bcd0945555377c8008d ./benchmarking.py 98b3c914732d1ff2b626e41de64619f058463956b997a7fad0eac6d41afb831d ./circuit.py b368fa72d84696531f234263a75b4f7fa4de528885f91cd205a679d7e61682d3 ./compare.py @@ -23,7 +23,7 @@ ccf0d83d58cf3645d17d8eef9b928471ee07958c1831a2a2fe47375dad7f4b2d ./examples/rep e9d3750e58d51dff274a505f82f24fa9f3b4540de3e696053443868f28c12f83 ./examples/reproduce_one_sample_optimality.py b0e02e239a277b63f7a4cd64182aea155052359e708d095a9be30e90bb872bdd ./examples/simulate_minimal_four_qubit_circuit.py 657a73bbc6425559834b35866ae103f8fe57f01417fd37ef7876fa0de6da01ed ./proba.py -0eea87da078b39bde7ebf921849fc4d0bd8d87f3a9062f1f0e13b730a03a4c7a ./pyproject.toml +a843ca2838df028e56d6eb19dfcdb9c3e8712f11a74fb3bf97e8998ea1bca729 ./pyproject.toml ece8e18e5ab0401e138785179855e47b4ac46560c4e578f6e1a5d44c45eb539f ./requirements-legacy.txt 06c182a1b9b0a31420d26fba2d45a224dc71a1c84bcc1fd55784830b89c9164f ./requirements.txt c4f52c8bf20bb65d0eba320aa1aef5d4b889e4c438156a42fc47307f251965f1 ./results/README.md @@ -32,21 +32,28 @@ fac5d6248ca774cb4d3fe506518bb36594bbdbaf31dc0c46aca0a0f3524a576a ./results/figu 6b366e9a1765ffa38e259c78729c80855e6162dcfa6fa5b258b10e18951b72e3 ./results/heralded_witness_results.json 02cb0b5ee8b241a556ce8c35a13829555aeeea14a0275f0bf869c06c2f58eeb6 ./results/one_sample_optimality.json da1935b2bd325cb523f313fa19a02270aa59e6acbcda00579a73a422b978294c ./results/validation_crosschecks.json -dec1ba0d29bc57db7507220340a34cf5ebe626212e8a670d80819604aee83e2f ./scripts/generate_results.py -3dfe8dc79a664ee9e8c11ae658981467fb54b8857080a8c66985e4d643fb1d2f ./scripts/validate_committed_results.py +1df70b11fc9a97b13aa9a91415508b0a0e40deaec8dc0f4fbfac0dd55125b458 ./scripts/generate_results.py +7f75e8208141a8f5e17418991900d55599dea78d123e5fe86687046d73e1c6f1 ./scripts/validate_committed_results.py +3b361b0be66dcda83398b443a4b39128548d5e9b3e2cb5623e2247cf1018ce70 ./scripts/verify_checksums.py 636c2fe2c246e6806e21982e1696675f4dadd295d46d00269bce26f7dafac211 ./scripts/verify_historical_files.py -134afb8cb64247983777b2698a5485686e88a2591936b55e3afb5945f05928d0 ./src/dcp_challenge/__init__.py +86b2a62c29e84eea3bf15f6476f81883a7ccf43952e65261d96535aadc414e3a ./src/dcp_challenge/__init__.py f520a4c550afa95ca11cf5a9b3cf656e4b0433c99b6f84a4088877d597062a35 ./src/dcp_challenge/exact_probabilities.py 0e5b18c70d8c71ebf1ee90b37e8970dba76658ec6d10507dc2cf1b0e5f573f16 ./src/dcp_challenge/heralded_witness.py +a6d6bfe51cc8d75a919487458e1ad622cd8ed5f2cb61559136d7318140e0ad13 ./src/dcp_challenge/integrity.py aa47265cc5fd3c4a2ca8f96b5b8c607ae35f74ec7db54d965c3dbfa7ea9e5f91 ./src/dcp_challenge/likelihood_decoder.py 4de797931f4767f99b2bce6f09d46ea115e337180465287dec66fec062f02293 ./src/dcp_challenge/minimal_circuit.py adc35aed3c8c0d8d905cb35e2afc4fb9200057bb7567c81cdcae74a8b223e1e9 ./src/dcp_challenge/original_protocol.py -d05eeb96dfb95747c49135ea9940ea31eb32e54db2d414723b91dbc9afd25f27 ./src/dcp_challenge/statistics.py +a8858a11540d6b009184a0511e3e5bbda6d04626373c8f900233bbd44e5408b7 ./src/dcp_challenge/reference_results.py +f92e394437cad537f44ef857f66b869bafa136fa66e0fb3f157d02e9c2731d13 ./src/dcp_challenge/result_validation.py +bd90bfbd09791bc87c218ff1ad2fa00b765c4d738198b5b66d1577b1a3161981 ./src/dcp_challenge/statistics.py 392a038abfb7ef5d14c445e11af42a942d97e75e322544ce16b045a9f372a5eb ./tests/test_exact_probabilities.py 61f2a3ef4c199710199f3dd378d85ba98e84ddaf5cca1017f7b9e0d86e1af36e ./tests/test_heralded_witness.py 862a3b357b0a2b9787274ee4186b060cbc7fe3d45c70ad510b46844138a98f12 ./tests/test_historical_files.py +1e02e495cfa86b53854279bd644382c1aa631a11dec92d691e1bf3e9b1dacaae ./tests/test_integrity.py d0202d343cb21e85e98b36241b2af4e261fd67c91603876a72f08404017b5f78 ./tests/test_likelihood_decoder.py ac3e495151c74a79eeb6a4cbcd135b6bddf7fd70e38ad6b1e3abd0bbb9a08075 ./tests/test_original_protocol.py +32d3c06a6756f27efc13f6cc25f386f75ba4b7ae2638a7e21dae66b93f0249a3 ./tests/test_result_validation.py +66b56beeb81973a82068e6c5f129e8a2f1f97d116ab90ef74167468174e9db54 ./tests/test_stable_statistics.py d482ef8a39d42dbfa8485abd79bb939699383928b544639abe5767697d83bd88 ./tests/test_statistics.py 1bccb732feb1a5403164d875e304e13446ebc7162cf86634141b59b11dc3c33c ./verification.py c4376bfd9dc74b1a1d781c062aab8d6a4c6150f89464e552138face53c1a654d ./witness/CIRCUIT.md diff --git a/CITATION.cff b/CITATION.cff index 5c26b14..a6ce1a7 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -5,8 +5,8 @@ message: >- corrected probability, soundness, or statistical analysis. title: "DCP Challenge: Author Correction and Follow-Up Analysis" type: software -version: 1.1.0-author-correction -date-released: 2026-08-17 +version: 1.1.1-author-correction +date-released: 2026-09-07 repository-code: "https://github.com/GoGoKo699/DCP-challenge" authors: - family-names: Lin diff --git a/README.md b/README.md index f4b2fa8..16a49c3 100644 --- a/README.md +++ b/README.md @@ -76,6 +76,7 @@ python examples/reproduce_counterexample.py python examples/reproduce_one_sample_optimality.py python examples/reproduce_figure5.py python scripts/verify_historical_files.py +python scripts/verify_checksums.py python scripts/validate_committed_results.py pytest -q ``` @@ -89,6 +90,14 @@ Same measurements, better decoder: 25/32 = 0.781250 See [REPRODUCIBILITY.md](REPRODUCIBILITY.md) for fixed seeds, exact-versus-numerical distinctions, and validation tolerances. The completed checks are summarized in [VALIDATION_RECORD.md](VALIDATION_RECORD.md). +### Numerical maintenance - 7 September 2026 + +Software version 1.1.1 adds stable binomial and log-probability calculations, +validation of every committed result field, deliberate-corruption tests, and +complete checksum verification. The August scientific correction, PDF, and +historical scripts are unchanged. See the [reproducibility guide](REPRODUCIBILITY.md) +for the distinction between recomputation and archived Monte Carlo checks. + ## Separate follow-up: phase-twirled heralded DCP witness The directory [witness/](witness/) contains a new 2026 follow-up. It is not part of the 2022 article and is not a general proof of quantum computation. diff --git a/REPRODUCIBILITY.md b/REPRODUCIBILITY.md index 621c345..0611f03 100644 --- a/REPRODUCIBILITY.md +++ b/REPRODUCIBILITY.md @@ -117,22 +117,84 @@ The tests cover: ```bash python scripts/validate_committed_results.py +python scripts/verify_checksums.py ``` -Validation policy: +The result validator checks all five JSON/CSV records, every field, duplicate +representations in the core JSON, schemas, and finite-number requirements. +It recomputes exact counterexamples, one-sample records, all Figure 5 analytical +probabilities, the exact Figure 5(a) likelihood result, witness summaries, +dense-matrix diagnostics, the 20,000-product-effect diagnostic, and tail examples. +Expected values are built without reading the files under validation. + +Integer and rational fields must agree exactly. Ordinary nonzero floating +quantities use relative tolerance `1e-12` and zero absolute tolerance. Only the +dense-matrix marginal residual permits absolute tolerance `1e-12`; this is not +used to round small probabilities to zero. Byte-identical floating-point +portability is not claimed. Validation uses explicit exceptions, so it also +runs under `python -O`. + +**Monte Carlo boundary:** fast validation does not repeat the million-trial +Figure 5(b) or the 60,000-trial Figure 5(c) experiment. It checks archived success +counts, seeds, trial counts, estimates, and recomputed Wilson intervals, including +the independent-seed cross-check records. These are archived evidence, not fresh +simulation results. Deliberate-corruption tests alter each of the 458 result +fields and require rejection, including changes made to duplicate records together. + +The checksum verifier additionally checks all tracked file bytes and complete +manifest coverage. For an extracted source ZIP without Git, it checks the file +tree excluding recognized local caches/build directories. Checksums detect +snapshot changes, not scientific correctness or authenticity. + +To regenerate outputs without overwriting the archived files: -- exact rational, integer, and combinatorial quantities are compared exactly; -- fixed Monte Carlo records include seeds, trial counts, estimates, and Wilson intervals; -- floating-point linear-algebra diagnostics are compared within the numerical tolerances encoded by the tests; -- byte-identical portability of floating-point JSON across all numerical libraries and processors is not claimed. +```bash +python scripts/generate_results.py --output /tmp/dcp-reference-results +``` -To regenerate the quick committed outputs: +To explicitly rerun the larger Figure 5 simulations, use: ```bash -python scripts/generate_results.py +python scripts/generate_results.py --monte-carlo --output /tmp/dcp-fresh-monte-carlo +``` + +The latter still retains the separate historical independent-seed cross-check +records. Fresh Monte Carlo output need not be byte-identical across numerical +environments; it must not silently replace the archived record or be confused +with the fast validator. Regeneration may change harmless last-bit diagnostics. + +## 8. Stable statistical evaluation + +The binomial routines use a scaled probability-mass recurrence and logarithms +instead of multiplying huge coefficients by tiny powers. They require no new +runtime dependency. For example, `binomial_lower_tail(2000, 400)` evaluates to +approximately `7.183777008e-8` instead of overflowing. + +```python +from dcp_challenge.statistics import ( + binomial_lower_tail, + binomial_log_lower_tail, + adaptive_sequential_log_lower_tail_bound, +) + +probability = binomial_lower_tail(2000, 400) +log_probability = binomial_log_lower_tail(10000, 0) +# log_probability is approximately -2876.8207245, not negative infinity. +log_bound = adaptive_sequential_log_lower_tail_bound(10000, 0) ``` -## 8. Historical software +All logarithms are natural. When a positive tail would underflow to zero, the +ordinary probability function raises `FloatingPointError` directing callers to +the log API. A zero return is reserved for an exact-zero endpoint. Counts must +be integers and probabilities must be finite and in `[0, 1]`. + +These are floating-point evaluations of the analytical bound, not certified +interval-arithmetic enclosures. They do not validate physical trust assumptions, +experimental stopping rules, or repeated selective restarts. Work and integer +storage increase with experiment size. Tests include exact rational references, +endpoint checks, tail monotonicity, and the previously overflowing cases. + +## 9. Historical software The original scripts require their 2022 Qibo environment and are not imported by the corrected test suite. `requirements-legacy.txt` records the historical dependency request. The proof-of-concept IBM experiment is preserved but is outside the present correction. diff --git a/VALIDATION_RECORD.md b/VALIDATION_RECORD.md index 6e2b613..20b2b7e 100644 --- a/VALIDATION_RECORD.md +++ b/VALIDATION_RECORD.md @@ -1,5 +1,18 @@ # Validation Record +## Software hardening - 7 September 2026 + +- Baseline tree: `93caf96908663f256de153feee54afccb358d342`, matching published `main` at `449548426b83fc22a9939fb2c2f6dede74dcd6f4`. +- Local suite: **59 tests passed**, including 458 individual result-field mutations. +- Previously overflowing binomial cases now agree with independent exact rational sums. +- Positive probability underflow is explicit; finite log tails remain available. +- All five committed result records and complete SHA-256 manifest coverage validate. +- Historical scripts and the August scientific documents/PDF are unchanged. +- The fast validator reuses archived large Monte Carlo counts and recomputes their arithmetic; it does not claim a fresh large simulation. +- The automated workflow repeats tests and integrity checks on Python 3.10, 3.11, 3.12, and 3.13. Its live execution status is recorded in GitHub Actions, not inferred from this document. + +## August scientific validation + **Validation date:** 17 August 2026 **Target:** author-correction repository state diff --git a/pyproject.toml b/pyproject.toml index 225d267..3ab0278 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "dcp-challenge-correction" -version = "1.1.0" +version = "1.1.1" description = "Corrected analytical tools and follow-up checks for the 2022 DCP challenge" readme = "README.md" requires-python = ">=3.10" diff --git a/scripts/generate_results.py b/scripts/generate_results.py index 1755095..4f7c0a4 100644 --- a/scripts/generate_results.py +++ b/scripts/generate_results.py @@ -1,216 +1,38 @@ -"""Regenerate the machine-readable results committed in ``results/``.""" +"""Regenerate reference outputs; large Monte Carlo runs are opt-in.""" from __future__ import annotations +import argparse import csv import json -import math -from fractions import Fraction from pathlib import Path -from typing import Any -from dcp_challenge.exact_probabilities import ( - exact_honest_success, - published_p_upper, - special_outcome_decoder_success, -) -from dcp_challenge.heralded_witness import ( - finite_sample_examples, - full_state_checks, - minimal_four_qubit_resources, - random_product_ratio_stress, - separable_saturation_probabilities, - witness_summary, -) -from dcp_challenge.likelihood_decoder import ( - exact_h_bayes_success, - explicit_two_sample_decoder_breakdown, - explicit_two_sample_decoder_success, - h_likelihood_monte_carlo, - one_sample_optimal_all_h_success, - one_sample_parity_difference, -) -from dcp_challenge.statistics import ( - adaptive_sequential_lower_tail_bound, - bernoulli_standard_error, - legacy_standard_error, -) +from dcp_challenge.reference_results import build_results ROOT = Path(__file__).resolve().parents[1] -RESULTS = ROOT / "results" -RESULTS.mkdir(exist_ok=True) -def fraction_record(value: Fraction) -> dict[str, float | int]: - return { - "numerator": value.numerator, - "denominator": value.denominator, - "decimal": float(value), - } - - -def canonicalize(value: Any) -> Any: - """Normalize floating diagnostics without pretending to exact bit portability.""" - if isinstance(value, dict): - return {key: canonicalize(item) for key, item in value.items()} - if isinstance(value, list): - return [canonicalize(item) for item in value] - if isinstance(value, tuple): - return [canonicalize(item) for item in value] - if isinstance(value, float): - if not math.isfinite(value): - return value - if abs(value) < 5e-15: - return 0.0 - return float(f"{value:.14g}") - return value - - -def build_figure5_rows(run_large_monte_carlo: bool) -> list[dict[str, object]]: - rows: list[dict[str, object]] = [] - for label, n, m, t in (("5a", 4, 6, 1), ("5b", 6, 9, 4), ("5c", 9, 21, 9)): - exact = exact_honest_success(n, m, t) - special = special_outcome_decoder_success(n, m, t) - row: dict[str, object] = { - "figure": label, - "n": n, - "N": 2**n, - "m": m, - "t": t, - "samples": m * t, - "p_B_special": float(special), - "p_upper_published": float(published_p_upper(n, m, t)), - "p_exact_honest": float(exact), - "exact_gap_honest_minus_p_B": float(exact - special), - } - if label == "5a": - result = exact_h_bayes_success(n, m * t) - row.update( - { - "all_H_full_likelihood": float(result), - "all_H_method": "exact", - "all_H_trials": "", - "all_H_seed": "", - "all_H_ci_low": "", - "all_H_ci_high": "", - } - ) - elif run_large_monte_carlo: - if label == "5b": - trials, seed, batch = 1_000_000, 20260817, 500 - else: - trials, seed, batch = 60_000, 20260818, 25 - result = h_likelihood_monte_carlo(n, m * t, trials, seed, batch) - row.update( - { - "all_H_full_likelihood": result.estimate, - "all_H_method": "Monte Carlo, Wilson 95% interval", - "all_H_trials": result.trials, - "all_H_seed": result.seed, - "all_H_ci_low": result.wilson_low, - "all_H_ci_high": result.wilson_high, - } - ) +def main(run_large_monte_carlo: bool = False, output: Path | None = None) -> None: + destination = output if output is not None else ROOT / "results" + destination.mkdir(parents=True, exist_ok=True) + records = build_results(run_large_monte_carlo) + for name, data in records.items(): + path = destination / name + if path.suffix == ".csv": + with path.open("w", newline="", encoding="utf-8") as handle: + writer = csv.DictWriter(handle, fieldnames=sorted(data[0])) + writer.writeheader() + writer.writerows(data) else: - reference = { - "5b": (0.937289, 1_000_000, 20260817, 0.9368121396967324, 0.9377625006608012), - "5c": (0.9571833333333334, 60_000, 20260818, 0.955533995616519, 0.9587741330997164), - }[label] - row.update( - { - "all_H_full_likelihood": reference[0], - "all_H_method": "frozen Monte Carlo, Wilson 95% interval", - "all_H_trials": reference[1], - "all_H_seed": reference[2], - "all_H_ci_low": reference[3], - "all_H_ci_high": reference[4], - } - ) - rows.append(row) - return rows - - -def one_sample_records() -> list[dict[str, object]]: - rows = [] - for n in range(1, 8): - N = 1 << n - difference = one_sample_parity_difference(n) - support = [index for index, value in enumerate(difference) if value] - rows.append( - { - "n": n, - "N": N, - "optimal_all_H_success": fraction_record(one_sample_optimal_all_h_success(n)), - "special_decoder_success": fraction_record(special_outcome_decoder_success(n, 1, 1)), - "difference_support_outcomes_rN_plus_y": support, - "difference_values": [fraction_record(difference[index]) for index in support], - } - ) - return rows - - -def main(run_large_monte_carlo: bool = False) -> None: - figure_rows = build_figure5_rows(run_large_monte_carlo) - fieldnames = sorted({key for row in figure_rows for key in row}) - with (RESULTS / "figure5_corrected.csv").open("w", newline="", encoding="utf-8") as handle: - writer = csv.DictWriter(handle, fieldnames=fieldnames) - writer.writeheader() - writer.writerows(figure_rows) - - p_figure5b = float(exact_honest_success(6, 9, 4)) - one_sample = one_sample_records() - witness = canonicalize( - { - "N2": witness_summary(1), - "N4": witness_summary(2), - "separable_saturating_product_effect_probabilities": separable_saturation_probabilities(), - "finite_sample_examples": finite_sample_examples(), - "adaptive_sequential_examples": { - "zero_errors_among_50": adaptive_sequential_lower_tail_bound(50, 0), - "at_most_20_errors_among_200": adaptive_sequential_lower_tail_bound(200, 20), - }, - "dense_full_state_N2": full_state_checks(1), - "dense_full_state_N4": full_state_checks(2), - "random_product_ratio_stress": random_product_ratio_stress(20_000, 20260817), - "minimal_four_qubit_resources": minimal_four_qubit_resources(), - } - ) - - core = { - "small_exact_counterexample": { - "published_special_decoder": fraction_record(special_outcome_decoder_success(2, 2, 1)), - "explicit_same_measurements_better_decoder": fraction_record(explicit_two_sample_decoder_success()), - "breakdown": { - key: fraction_record(value) - for key, value in explicit_two_sample_decoder_breakdown().items() - }, - }, - "one_sample_optimality": one_sample, - "figure5": figure_rows, - "standard_error_example_figure5b_r1000": { - "p": p_figure5b, - "legacy_formula": legacy_standard_error(p_figure5b, 1000), - "correct_bernoulli_formula": bernoulli_standard_error(p_figure5b, 1000), - }, - "heralded_witness": witness, - "validation_policy": { - "exact_quantities": "exact rational or integer comparison", - "monte_carlo": "committed seeds, trial counts, Wilson intervals, and independent-seed cross-checks", - "floating_linear_algebra": "numerical tolerances encoded in the tests; byte-identical portability is not claimed", - }, - } - - with (RESULTS / "correction_core_results.json").open("w", encoding="utf-8") as handle: - json.dump(core, handle, indent=2, sort_keys=True) - handle.write("\n") - with (RESULTS / "one_sample_optimality.json").open("w", encoding="utf-8") as handle: - json.dump(one_sample, handle, indent=2, sort_keys=True) - handle.write("\n") - with (RESULTS / "heralded_witness_results.json").open("w", encoding="utf-8") as handle: - json.dump(witness, handle, indent=2, sort_keys=True) - handle.write("\n") - - print(json.dumps(core, indent=2, sort_keys=True)) + path.write_text(json.dumps(data, indent=2, sort_keys=True, allow_nan=False) + "\n", + encoding="utf-8") + print(f"Wrote {len(records)} result files to {destination}.") + if not run_large_monte_carlo: + print("Monte Carlo counts are archived records, not newly simulated in this mode.") if __name__ == "__main__": - main() + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, help="separate directory for regenerated files") + parser.add_argument("--monte-carlo", action="store_true", help="rerun the large Figure 5 simulations") + args = parser.parse_args() + main(args.monte_carlo, args.output) diff --git a/scripts/validate_committed_results.py b/scripts/validate_committed_results.py index 5f879a2..ee6427e 100644 --- a/scripts/validate_committed_results.py +++ b/scripts/validate_committed_results.py @@ -1,58 +1,22 @@ -"""Fast validation of committed exact values and numerical tolerances.""" +"""Validate every committed scientific result field, without modifying files.""" from __future__ import annotations -import csv -import json -import math -from fractions import Fraction +import argparse from pathlib import Path -from dcp_challenge.exact_probabilities import exact_honest_success, special_outcome_decoder_success -from dcp_challenge.heralded_witness import full_state_checks, witness_summary -from dcp_challenge.likelihood_decoder import ( - explicit_two_sample_decoder_success, - one_sample_optimal_all_h_success, -) +from dcp_challenge.result_validation import validate_results ROOT = Path(__file__).resolve().parents[1] -RESULTS = ROOT / "results" -def assert_close(actual: float, expected: float, tolerance: float = 1e-12) -> None: - if not math.isclose(actual, expected, rel_tol=tolerance, abs_tol=tolerance): - raise AssertionError(f"{actual!r} != {expected!r} within {tolerance}") - - -def main() -> None: - core = json.loads((RESULTS / "correction_core_results.json").read_text()) - small = core["small_exact_counterexample"] - assert (small["published_special_decoder"]["numerator"], small["published_special_decoder"]["denominator"]) == (23, 32) - assert (small["explicit_same_measurements_better_decoder"]["numerator"], small["explicit_same_measurements_better_decoder"]["denominator"]) == (25, 32) - assert explicit_two_sample_decoder_success() == Fraction(25, 32) - - one_sample = json.loads((RESULTS / "one_sample_optimality.json").read_text()) - for row in one_sample: - n = int(row["n"]) - expected = one_sample_optimal_all_h_success(n) - stored = row["optimal_all_H_success"] - assert (stored["numerator"], stored["denominator"]) == (expected.numerator, expected.denominator) - assert expected == special_outcome_decoder_success(n, 1, 1) - - with (RESULTS / "figure5_corrected.csv").open(newline="", encoding="utf-8") as handle: - rows = list(csv.DictReader(handle)) - for row in rows: - n, m, t = int(row["n"]), int(row["m"]), int(row["t"]) - assert_close(float(row["p_exact_honest"]), float(exact_honest_success(n, m, t)), 1e-14) - - witness = json.loads((RESULTS / "heralded_witness_results.json").read_text()) - assert_close(witness["N2"]["honest_correct_conclusive"], witness_summary(1)["honest_correct_conclusive"]) - for n, key in ((1, "dense_full_state_N2"), (2, "dense_full_state_N4")): - fresh = full_state_checks(n) - assert_close(witness[key]["trace_norm_difference"], fresh["trace_norm_difference"], 1e-12) - assert witness[key]["local_marginal_max_error_from_maximally_mixed"] <= 1e-12 - - print("Committed result validation passed.") +def main(results_dir: Path = ROOT / "results") -> None: + leaves = validate_results(results_dir) + print(f"Committed result validation passed: 5 files, {leaves} fields.") + print("Exact values recomputed; floating diagnostics checked with declared tolerances.") + print("Monte Carlo records: archived counts/metadata and Wilson intervals checked; not re-simulated.") if __name__ == "__main__": - main() + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--results-dir", type=Path, default=ROOT / "results") + main(parser.parse_args().results_dir) diff --git a/scripts/verify_checksums.py b/scripts/verify_checksums.py new file mode 100644 index 0000000..0a56f29 --- /dev/null +++ b/scripts/verify_checksums.py @@ -0,0 +1,8 @@ +"""Verify all tracked file checksums and manifest coverage.""" +from pathlib import Path + +from dcp_challenge.integrity import verify_checksums + +if __name__ == "__main__": + count = verify_checksums(Path(__file__).resolve().parents[1]) + print(f"Checksum verification passed: {count} files; complete manifest coverage.") diff --git a/src/dcp_challenge/__init__.py b/src/dcp_challenge/__init__.py index b8e510f..8c294e9 100644 --- a/src/dcp_challenge/__init__.py +++ b/src/dcp_challenge/__init__.py @@ -26,4 +26,4 @@ "one_sample_parity_distributions", ] -__version__ = "1.1.0" +__version__ = "1.1.1" diff --git a/src/dcp_challenge/integrity.py b/src/dcp_challenge/integrity.py new file mode 100644 index 0000000..944bb83 --- /dev/null +++ b/src/dcp_challenge/integrity.py @@ -0,0 +1,63 @@ +"""Verify SHA-256 contents and manifest coverage, including exported ZIPs.""" +from __future__ import annotations + +import hashlib +import re +import subprocess +from pathlib import Path, PurePosixPath + + +class IntegrityError(ValueError): + """The manifest or repository content failed its integrity check.""" + + +def _repository_files(root: Path) -> set[str]: + try: + top = subprocess.run(["git", "-C", str(root), "rev-parse", "--show-toplevel"], + capture_output=True, text=True, check=True) + if Path(top.stdout.strip()).resolve() == root.resolve(): + files = subprocess.run(["git", "-C", str(root), "ls-files", "-z"], + capture_output=True, check=True).stdout.decode("utf-8") + return set(filter(None, files.split("\0"))) - {"CHECKSUMS.sha256"} + except (OSError, subprocess.CalledProcessError): + pass + # In source ZIPs, ignore only recognizable local build/cache directories. + ignored = {".git", "__pycache__", ".pytest_cache", ".venv", "venv", "build", "dist"} + return {p.relative_to(root).as_posix() for p in root.rglob("*") + if p.is_file() and p.relative_to(root).as_posix() != "CHECKSUMS.sha256" + and not any(part in ignored or part.endswith(".egg-info") + for part in p.relative_to(root).parts)} + + +def verify_checksums(root: Path) -> int: + """Check every manifest entry and require coverage of all tracked files. + + This checks snapshot integrity, not authenticity: an author can update both + a file and its checksum. Scientific consistency is checked separately. + """ + root = root.resolve() + manifest = root / "CHECKSUMS.sha256" + entries: dict[str, str] = {} + for number, line in enumerate(manifest.read_text(encoding="utf-8").splitlines(), 1): + match = re.fullmatch(r"([0-9a-f]{64}) (.+)", line) + if not match: + raise IntegrityError(f"invalid manifest line {number}") + digest, name = match.groups() + path = PurePosixPath(name) + if path.is_absolute() or ".." in path.parts or "\\" in name: + raise IntegrityError(f"unsafe manifest path: {name}") + name = path.as_posix() + if name in entries or name == "CHECKSUMS.sha256": + raise IntegrityError(f"duplicate or self-referential entry: {name}") + target = root / name + if (not target.is_file() or target.is_symlink() + or root not in target.resolve().parents): + raise IntegrityError(f"missing or unsafe file: {name}") + if hashlib.sha256(target.read_bytes()).hexdigest() != digest: + raise IntegrityError(f"checksum mismatch: {name}") + entries[name] = digest + expected = _repository_files(root) + if set(entries) != expected: + raise IntegrityError(f"manifest coverage mismatch; unlisted={sorted(expected-set(entries))}, " + f"extra={sorted(set(entries)-expected)}") + return len(entries) diff --git a/src/dcp_challenge/reference_results.py b/src/dcp_challenge/reference_results.py new file mode 100644 index 0000000..7d700d3 --- /dev/null +++ b/src/dcp_challenge/reference_results.py @@ -0,0 +1,219 @@ +"""Recompute reference results without reading the files being validated. + +Exact results and witness diagnostics are recomputed. The larger Monte Carlo +counts are archived reference data: default validation checks their provenance +and arithmetic, NOT a fresh replay of the million-trial simulations. +""" +from __future__ import annotations + +import math +from fractions import Fraction +from typing import Any + +from dcp_challenge.exact_probabilities import ( + exact_honest_success, + published_p_upper, + special_outcome_decoder_success, +) +from dcp_challenge.heralded_witness import ( + finite_sample_examples, + full_state_checks, + minimal_four_qubit_resources, + random_product_ratio_stress, + separable_saturation_probabilities, + witness_summary, +) +from dcp_challenge.likelihood_decoder import ( + exact_h_bayes_success, + explicit_two_sample_decoder_breakdown, + explicit_two_sample_decoder_success, + h_likelihood_monte_carlo, + one_sample_optimal_all_h_success, + one_sample_parity_difference, +) +from dcp_challenge.statistics import ( + adaptive_sequential_lower_tail_bound, + bernoulli_standard_error, + legacy_standard_error, + wilson_interval, +) + + + +def fraction_record(value: Fraction) -> dict[str, float | int]: + return { + "numerator": value.numerator, + "denominator": value.denominator, + "decimal": float(value), + } + + +def canonicalize(value: Any) -> Any: + """Normalize floating diagnostics without pretending to exact bit portability.""" + if isinstance(value, dict): + return {key: canonicalize(item) for key, item in value.items()} + if isinstance(value, list): + return [canonicalize(item) for item in value] + if isinstance(value, tuple): + return [canonicalize(item) for item in value] + if isinstance(value, float): + if not math.isfinite(value): + return value + return float(f"{value:.14g}") + return value + + +def build_figure5_rows(run_large_monte_carlo: bool) -> list[dict[str, object]]: + rows: list[dict[str, object]] = [] + for label, n, m, t in (("5a", 4, 6, 1), ("5b", 6, 9, 4), ("5c", 9, 21, 9)): + exact = exact_honest_success(n, m, t) + special = special_outcome_decoder_success(n, m, t) + row: dict[str, object] = { + "figure": label, + "n": n, + "N": 2**n, + "m": m, + "t": t, + "samples": m * t, + "p_B_special": float(special), + "p_upper_published": float(published_p_upper(n, m, t)), + "p_exact_honest": float(exact), + "exact_gap_honest_minus_p_B": float(exact - special), + } + if label == "5a": + result = exact_h_bayes_success(n, m * t) + row.update( + { + "all_H_full_likelihood": float(result), + "all_H_method": "exact", + "all_H_trials": "", + "all_H_seed": "", + "all_H_ci_low": "", + "all_H_ci_high": "", + } + ) + elif run_large_monte_carlo: + if label == "5b": + trials, seed, batch = 1_000_000, 20260817, 500 + else: + trials, seed, batch = 60_000, 20260818, 25 + result = h_likelihood_monte_carlo(n, m * t, trials, seed, batch) + row.update( + { + "all_H_full_likelihood": result.estimate, + "all_H_method": "Monte Carlo, Wilson 95% interval", + "all_H_trials": result.trials, + "all_H_seed": result.seed, + "all_H_ci_low": result.wilson_low, + "all_H_ci_high": result.wilson_high, + } + ) + else: + successes, trials, seed = { + "5b": (937289, 1_000_000, 20260817), + "5c": (57431, 60_000, 20260818), + }[label] + low, high = wilson_interval(successes, trials) + reference = (successes / trials, trials, seed, low, high) + row.update( + { + "all_H_full_likelihood": reference[0], + "all_H_method": "frozen Monte Carlo, Wilson 95% interval", + "all_H_trials": reference[1], + "all_H_seed": reference[2], + "all_H_ci_low": reference[3], + "all_H_ci_high": reference[4], + } + ) + rows.append(row) + return rows + + +def one_sample_records() -> list[dict[str, object]]: + rows = [] + for n in range(1, 8): + N = 1 << n + difference = one_sample_parity_difference(n) + support = [index for index, value in enumerate(difference) if value] + rows.append( + { + "n": n, + "N": N, + "optimal_all_H_success": fraction_record(one_sample_optimal_all_h_success(n)), + "special_decoder_success": fraction_record(special_outcome_decoder_success(n, 1, 1)), + "difference_support_outcomes_rN_plus_y": support, + "difference_values": [fraction_record(difference[index]) for index in support], + } + ) + return rows + + +def build_results(run_large_monte_carlo: bool = False) -> dict[str, Any]: + figure_rows = build_figure5_rows(run_large_monte_carlo) + p_figure5b = float(exact_honest_success(6, 9, 4)) + one_sample = one_sample_records() + witness = canonicalize( + { + "N2": witness_summary(1), + "N4": witness_summary(2), + "separable_saturating_product_effect_probabilities": separable_saturation_probabilities(), + "finite_sample_examples": finite_sample_examples(), + "adaptive_sequential_examples": { + "zero_errors_among_50": adaptive_sequential_lower_tail_bound(50, 0), + "at_most_20_errors_among_200": adaptive_sequential_lower_tail_bound(200, 20), + }, + "dense_full_state_N2": full_state_checks(1), + "dense_full_state_N4": full_state_checks(2), + "random_product_ratio_stress": random_product_ratio_stress(20_000, 20260817), + "minimal_four_qubit_resources": minimal_four_qubit_resources(), + } + ) + + core = { + "small_exact_counterexample": { + "published_special_decoder": fraction_record(special_outcome_decoder_success(2, 2, 1)), + "explicit_same_measurements_better_decoder": fraction_record(explicit_two_sample_decoder_success()), + "breakdown": { + key: fraction_record(value) + for key, value in explicit_two_sample_decoder_breakdown().items() + }, + }, + "one_sample_optimality": one_sample, + "figure5": figure_rows, + "standard_error_example_figure5b_r1000": { + "p": p_figure5b, + "legacy_formula": legacy_standard_error(p_figure5b, 1000), + "correct_bernoulli_formula": bernoulli_standard_error(p_figure5b, 1000), + }, + "heralded_witness": witness, + "validation_policy": { + "exact_quantities": "exact rational or integer comparison", + "monte_carlo": "committed seeds, trial counts, Wilson intervals, and independent-seed cross-checks", + "floating_linear_algebra": "numerical tolerances encoded in the tests; byte-identical portability is not claimed", + }, + } + + return { + "correction_core_results.json": core, + "one_sample_optimality.json": one_sample, + "heralded_witness_results.json": witness, + "figure5_corrected.csv": figure_rows, + "validation_crosschecks.json": archived_crosschecks(), + } + + +def archived_crosschecks() -> dict[str, Any]: + """Arithmetic reconstruction of archived independent-seed records.""" + record: dict[str, Any] = { + "purpose": "Independent-seed Monte Carlo cross-checks of the frozen Figure 5 likelihood-decoder estimates", + } + for label, n, samples, trials, seed, successes in ( + ("figure5b", 6, 36, 20000, 20260821, 18789), + ("figure5c", 9, 189, 5000, 20260822, 4789), + ): + low, high = wilson_interval(successes, trials) + record[label] = dict(n=n, samples=samples, trials=trials, seed=seed, + successes=successes, estimate=successes/trials, + wilson_low=low, wilson_high=high) + record["interpretation"] = "Both independent intervals contain or overlap the corresponding frozen estimates committed in figure5_corrected.csv." + return record diff --git a/src/dcp_challenge/result_validation.py b/src/dcp_challenge/result_validation.py new file mode 100644 index 0000000..dfa64cb --- /dev/null +++ b/src/dcp_challenge/result_validation.py @@ -0,0 +1,112 @@ +"""Strict, read-only validation of every field in the committed result records. + +The expected values are generated without reading the result files. Frozen +Monte Carlo counts are checked as archived evidence, not newly sampled data. +""" +from __future__ import annotations + +import csv +import json +import math +from functools import lru_cache +from pathlib import Path +from typing import Any + +from .reference_results import build_results + + +class ValidationError(ValueError): + """A result is missing, malformed, or inconsistent with its reference.""" + + +@lru_cache(maxsize=1) +def _expected_results() -> dict[str, Any]: + return build_results() + + +def compare_records(actual: Any, expected: Any, path: str = "result") -> int: + """Compare structure and every leaf; return the number of validated leaves. + + Integers and text are exact. Nonzero floats use relative tolerance 1e-12 + with no blanket absolute tolerance (important for small p-values). Only + dense-matrix marginal residuals use an absolute tolerance of 1e-12. + """ + if isinstance(expected, dict): + if not isinstance(actual, dict) or actual.keys() != expected.keys(): + raise ValidationError(f"{path}: missing or unexpected fields") + return sum(compare_records(actual[k], v, f"{path}.{k}") for k, v in expected.items()) + if isinstance(expected, (list, tuple)): + if not isinstance(actual, list) or len(actual) != len(expected): + raise ValidationError(f"{path}: wrong list length or type") + return sum(compare_records(a, e, f"{path}[{i}]") for i, (a, e) in enumerate(zip(actual, expected))) + if isinstance(expected, float): + if type(actual) not in (int, float) or not math.isfinite(actual): + raise ValidationError(f"{path}: expected a finite number") + residual = path.endswith(".local_marginal_max_error_from_maximally_mixed") + if residual and not 0 <= actual <= 1e-12: + raise ValidationError(f"{path}: invalid marginal residual") + if not math.isclose(actual, expected, rel_tol=1e-12, abs_tol=1e-12 if residual else 0.0): + raise ValidationError(f"{path}: {actual!r} != {expected!r}") + elif type(actual) is not type(expected) or actual != expected: + raise ValidationError(f"{path}: {actual!r} != {expected!r}") + return 1 + + +def _unique_object(pairs): + result = {} + for key, value in pairs: + if key in result: + raise ValidationError(f"duplicate JSON key: {key}") + result[key] = value + return result + + +def _reject_constant(value): + raise ValidationError(f"non-finite JSON constant: {value}") + + +def _load_csv(path: Path, expected: list[dict]) -> list[dict]: + with path.open(newline="", encoding="utf-8") as handle: + reader = csv.DictReader(handle) + headers = reader.fieldnames or [] + if len(headers) != len(set(headers)) or set(headers) != set(expected[0]): + raise ValidationError(f"{path.name}: wrong or duplicate CSV columns") + rows = list(reader) + if len(rows) != len(expected): + raise ValidationError(f"{path.name}: wrong row count") + for row, reference in zip(rows, expected): + if set(row) != set(reference): + raise ValidationError(f"{path.name}: malformed CSV row") + for key, target in reference.items(): + try: + if isinstance(target, int): + row[key] = int(row[key]) + elif isinstance(target, float): + row[key] = float(row[key]) + except (TypeError, ValueError) as exc: + raise ValidationError(f"{path.name}.{key}: invalid numeric cell") from exc + return rows + + +def validate_results(results_dir: Path) -> int: + """Validate all five result files without changing them or the repository.""" + expected = _expected_results() + actual_names = {p.name for p in results_dir.iterdir() if p.suffix in (".json", ".csv")} + if actual_names != set(expected): + raise ValidationError("missing or unexpected result files") + leaves = 0 + for name, reference in expected.items(): + path = results_dir / name + if path.is_symlink() or not path.is_file(): + raise ValidationError(f"{name}: expected a regular file") + try: + if path.suffix == ".csv": + actual = _load_csv(path, reference) + else: + actual = json.loads(path.read_text(encoding="utf-8"), + object_pairs_hook=_unique_object, + parse_constant=_reject_constant) + leaves += compare_records(actual, reference, name) + except (OSError, json.JSONDecodeError) as exc: + raise ValidationError(f"{name}: cannot read valid result data") from exc + return leaves diff --git a/src/dcp_challenge/statistics.py b/src/dcp_challenge/statistics.py index 0557050..49cfc42 100644 --- a/src/dcp_challenge/statistics.py +++ b/src/dcp_challenge/statistics.py @@ -1,7 +1,8 @@ """Statistical utilities used by the corrected analysis.""" from __future__ import annotations -from math import comb, sqrt +from math import comb, exp, expm1, fsum, isfinite, log, log1p, sqrt +from operator import index from typing import Tuple @@ -44,23 +45,89 @@ def wilson_interval( return center - radius, center + radius +def _tail_parameters(conclusive: int, errors: int, probability: float) -> tuple[int, int, float]: + """Reject fractional counts, booleans, NaNs, and out-of-domain inputs.""" + if isinstance(conclusive, bool) or isinstance(errors, bool): + raise TypeError("counts must be integers, not booleans") + try: + n, k = index(conclusive), index(errors) + except TypeError as exc: + raise TypeError("counts must be integers") from exc + if n < 0 or not 0 <= k <= n: + raise ValueError("require 0 <= errors <= conclusive") + p = float(probability) + if not isfinite(p) or not 0.0 <= p <= 1.0: + raise ValueError("null error probability must be finite and lie in [0, 1]") + return n, k, p + + +def _log_left_tail(n: int, k: int, log_p: float, log_q: float) -> float: + """Log of a lower tail below the mean, scaled by its largest term. + + ``comb`` remains a Python integer: only its logarithm is converted to a + float. The downward PMF ratios are at most one. Summation therefore cannot + overflow or lose the entire tail when the unscaled probability underflows. + """ + def scaled_terms(): + term = 1.0 + yield term + for j in range(k, 0, -1): + term *= exp(fsum((log(j), -log(n - j + 1), log_q, -log_p))) + if term == 0.0: + break + yield term + + log_mass = fsum((log(comb(n, k)), k * log_p, (n - k) * log_q)) + return min(0.0, fsum((log_mass, log(fsum(scaled_terms()))))) + + +def binomial_log_lower_tail( + conclusive: int, + errors: int, + null_error_probability: float = 0.25, +) -> float: + """Natural logarithm of ``P[Binomial(C,p) <= E]``. + + Uses a scaled PMF recurrence, or the reflected upper tail when above the + mean. No SciPy dependency is needed. Exact-zero endpoint probabilities + return ``-inf``; positive tails remain finite even below float range. + These are floating-point evaluations, not interval-arithmetic certificates. + Work and integer storage increase with the experiment size. + """ + n, k, p = _tail_parameters(conclusive, errors, null_error_probability) + if k == n or p == 0.0: + return 0.0 + if p == 1.0: + return -float("inf") + log_p, log_q = log(p), log1p(-p) + if k < n * p: + return _log_left_tail(n, k, log_p, log_q) + # P[X > k] = P[Binomial(n,1-p) <= n-k-1]. Swap logs rather than + # forming 1-p again, which could round to one for very small p. + log_upper = _log_left_tail(n, n - k - 1, log_q, log_p) + if log_upper < -log(2.0): + return log1p(-exp(log_upper)) + return log(-expm1(log_upper)) + + def binomial_lower_tail( conclusive: int, errors: int, null_error_probability: float = 0.25, ) -> float: - """Binomial lower tail used by the adaptive sequential witness bound.""" - if conclusive < 0 or not 0 <= errors <= conclusive: - raise ValueError("require 0 <= errors <= conclusive") - if not 0.0 <= null_error_probability <= 1.0: - raise ValueError("null error probability must lie in [0, 1]") - return sum( - comb(conclusive, j) - * null_error_probability**j - * (1.0 - null_error_probability) ** (conclusive - j) - for j in range(errors + 1) - ) + """Binomial lower tail, with explicit protection against false zeros. + If a positive probability would underflow to zero, raise + ``FloatingPointError`` and direct the caller to the log-probability API. + A returned zero is reserved for a mathematically impossible endpoint. + """ + log_probability = binomial_log_lower_tail(conclusive, errors, null_error_probability) + probability = exp(log_probability) + if probability == 0.0 and isfinite(log_probability): + raise FloatingPointError( + "positive tail is below floating-point range; use binomial_log_lower_tail" + ) + return probability def adaptive_sequential_lower_tail_bound( @@ -68,15 +135,19 @@ def adaptive_sequential_lower_tail_bound( errors: int, conditional_error_lower_bound: float = 0.25, ) -> float: - """Upper bound for an adaptive sequence with conditional error at least ``p``. + """Evaluate the adaptive bound under the sequential null in witness/PROOF.md. - Independence is not assumed. If every successive conclusive answer has - conditional error probability at least ``p`` given all previous public - history, then its error count stochastically dominates ``Binomial(C,p)``; - hence the probability of at most ``E`` errors is bounded by this lower tail. + This evaluates the analytical binomial bound; it does not test the physical + assumptions or an arbitrary stopping/restarting rule. See the protocol's + predeclared (C,E,R) rule. Positive underflow raises FloatingPointError. """ - return binomial_lower_tail( - conclusive, - errors, - null_error_probability=conditional_error_lower_bound, - ) + return binomial_lower_tail(conclusive, errors, conditional_error_lower_bound) + + +def adaptive_sequential_log_lower_tail_bound( + conclusive: int, + errors: int, + conditional_error_lower_bound: float = 0.25, +) -> float: + """Natural-log version of the same analytical adaptive bound.""" + return binomial_log_lower_tail(conclusive, errors, conditional_error_lower_bound) diff --git a/tests/test_integrity.py b/tests/test_integrity.py new file mode 100644 index 0000000..cbc0374 --- /dev/null +++ b/tests/test_integrity.py @@ -0,0 +1,35 @@ +import hashlib +from pathlib import Path + +import pytest + +from dcp_challenge.integrity import IntegrityError, verify_checksums + + +def _fixture(root: Path): + (root / "README.md").write_text("verified bytes\n") + digest = hashlib.sha256((root / "README.md").read_bytes()).hexdigest() + (root / "CHECKSUMS.sha256").write_text(f"{digest} ./README.md\n") + + +def test_zip_manifest_validates_without_git(tmp_path): + _fixture(tmp_path) + assert verify_checksums(tmp_path) == 1 + + +@pytest.mark.parametrize("damage", ["changed", "missing", "unlisted", "duplicate", "unsafe"]) +def test_checksum_mutations_fail(tmp_path, damage): + _fixture(tmp_path) + if damage == "changed": + (tmp_path / "README.md").write_text("different bytes") + elif damage == "missing": + (tmp_path / "README.md").unlink() + elif damage == "unlisted": + (tmp_path / "unreviewed.py").write_text("pass") + elif damage == "duplicate": + manifest = tmp_path / "CHECKSUMS.sha256" + manifest.write_text(manifest.read_text()*2) + elif damage == "unsafe": + (tmp_path / "CHECKSUMS.sha256").write_text("0"*64 + " ../outside\n") + with pytest.raises(IntegrityError): + verify_checksums(tmp_path) diff --git a/tests/test_result_validation.py b/tests/test_result_validation.py new file mode 100644 index 0000000..07f13df --- /dev/null +++ b/tests/test_result_validation.py @@ -0,0 +1,153 @@ +"""Mutation tests: correct data must pass and corrupted data must fail.""" +import copy +import csv +import json +import shutil +import subprocess +import sys +from pathlib import Path + +import pytest + +from dcp_challenge.result_validation import ValidationError, validate_results + +ROOT = Path(__file__).resolve().parents[1] + + +@pytest.fixture +def records(tmp_path): + destination = tmp_path / "results" + shutil.copytree(ROOT / "results", destination) + return destination + + +def _leaves(value, prefix=()): + if isinstance(value, dict): + for key, child in value.items(): + yield from _leaves(child, prefix+(key,)) + elif isinstance(value, list): + for i, child in enumerate(value): + yield from _leaves(child, prefix+(i,)) + else: + yield prefix, value + + +def _alter(value): + if isinstance(value, int): + return value + max(1, abs(value)//10) + if isinstance(value, float): + return value + max(0.125, abs(value)*0.1) + return value + " [deliberate corruption]" + + +def test_real_records_pass(records): + assert validate_results(records) == 458 + + +def test_every_json_leaf_is_checked(records): + mutations = 0 + for path in sorted(records.glob("*.json")): + original = path.read_text() + data = json.loads(original) + for keys, value in _leaves(data): + corrupted = copy.deepcopy(data) + target = corrupted + for key in keys[:-1]: + target = target[key] + target[keys[-1]] = _alter(value) + path.write_text(json.dumps(corrupted)) + with pytest.raises(ValidationError): + validate_results(records) + mutations += 1 + path.write_text(original) + assert mutations == 410 + + +def test_every_csv_cell_is_checked(records): + path = records / "figure5_corrected.csv" + original = path.read_text() + reader = csv.DictReader(original.splitlines()) + rows, headers = list(reader), reader.fieldnames + for i, row in enumerate(rows): + for key in row: + corrupted = copy.deepcopy(rows) + corrupted[i][key] = "999999" if row[key] != "999999" else "0" + with path.open("w", newline="") as handle: + writer = csv.DictWriter(handle, fieldnames=headers) + writer.writeheader() + writer.writerows(corrupted) + with pytest.raises(ValidationError): + validate_results(records) + path.write_text(original) + + +def test_corrupting_duplicate_records_together_still_fails(records): + for name in ("heralded_witness_results.json", "correction_core_results.json"): + path = records / name + data = json.loads(path.read_text()) + witness = data if name.startswith("heralded") else data["heralded_witness"] + witness["N2"]["separable_confidence_upper_bound"] = 0.99 + path.write_text(json.dumps(data)) + with pytest.raises(ValidationError, match="confidence"): + validate_results(records) + + +@pytest.mark.parametrize("damage", ["missing_file", "extra_file", "missing_key", "extra_key", + "empty_rows", "duplicate_key", "nan", "infinity", "bool_as_count", "extra_csv_row"]) +def test_schema_and_nonfinite_corruption(records, damage): + path = records / "one_sample_optimality.json" + data = json.loads(path.read_text()) + if damage == "missing_file": + path.unlink() + elif damage == "extra_file": + (records / "unreviewed.json").write_text("{}") + elif damage == "missing_key": + del data[0]["n"] + elif damage == "extra_key": + data[0]["unreviewed"] = 1 + elif damage == "empty_rows": + data = [] + elif damage == "duplicate_key": + path.write_text('[{"n":1,"n":2}]') + elif damage == "nan": + data[0]["optimal_all_H_success"]["decimal"] = float("nan") + elif damage == "infinity": + data[0]["optimal_all_H_success"]["decimal"] = float("inf") + elif damage == "bool_as_count": + data[0]["n"] = True + elif damage == "extra_csv_row": + csv_path = records / "figure5_corrected.csv" + csv_path.write_text(csv_path.read_text()+csv_path.read_text().splitlines()[1]+"\n") + if damage not in ("missing_file", "duplicate_key"): + path.write_text(json.dumps(data)) + with pytest.raises(ValidationError): + validate_results(records) + + +def test_small_p_value_cannot_be_rounded_to_zero(records): + path = records / "heralded_witness_results.json" + data = json.loads(path.read_text()) + data["finite_sample_examples"]["at_most_20_errors_among_200"] = 0 + path.write_text(json.dumps(data)) + with pytest.raises(ValidationError): + validate_results(records) + + +def test_matrix_residual_tolerance_is_not_a_probability_tolerance(records): + path = records / "heralded_witness_results.json" + data = json.loads(path.read_text()) + data["dense_full_state_N2"]["local_marginal_max_error_from_maximally_mixed"] = 1e-13 + path.write_text(json.dumps(data)) + validate_results(records) + data["N2"]["honest_wrong_conclusive"] = 1e-13 + path.write_text(json.dumps(data)) + with pytest.raises(ValidationError): + validate_results(records) + + +def test_validation_remains_enabled_under_python_optimization(records): + (records / "one_sample_optimality.json").write_text("[]") + result = subprocess.run([sys.executable, "-O", str(ROOT / "scripts/validate_committed_results.py"), + "--results-dir", str(records)], capture_output=True, text=True) + assert result.returncode != 0 + assert "ValidationError" in result.stderr diff --git a/tests/test_stable_statistics.py b/tests/test_stable_statistics.py new file mode 100644 index 0000000..468da5b --- /dev/null +++ b/tests/test_stable_statistics.py @@ -0,0 +1,75 @@ +"""Independent rational references for overflow, underflow, and tail boundaries.""" +import math +from fractions import Fraction + +import pytest + +from dcp_challenge.statistics import ( + adaptive_sequential_log_lower_tail_bound, + adaptive_sequential_lower_tail_bound, + binomial_log_lower_tail, + binomial_lower_tail, +) + + +def exact_tail(n, k, p): + a, b = p.numerator, p.denominator + return Fraction(sum(math.comb(n, j) * a**j * (b-a)**(n-j) + for j in range(k+1)), b**n) + + +def test_small_tails_against_exact_rationals(): + for n in range(26): + for k in range(n+1): + for p in (Fraction(0), Fraction(1, 8), Fraction(1, 4), + Fraction(1, 2), Fraction(3, 4), Fraction(7, 8), Fraction(1)): + expected = float(exact_tail(n, k, p)) + value = binomial_lower_tail(n, k, float(p)) + assert 0 <= value <= 1 + assert math.isclose(value, expected, rel_tol=2e-13, abs_tol=0) + log_value = binomial_log_lower_tail(n, k, float(p)) + assert log_value <= 0 + assert math.isclose(math.exp(log_value), expected, rel_tol=2e-13, abs_tol=0) + + +@pytest.mark.parametrize("n,k", [(1500, 300), (2000, 400), (2000, 1600)]) +def test_previously_overflowing_cases(n, k): + expected = float(exact_tail(n, k, Fraction(1, 4))) + assert math.isclose(binomial_lower_tail(n, k), expected, rel_tol=2e-12) + assert adaptive_sequential_lower_tail_bound(n, k) == binomial_lower_tail(n, k) + + +def test_positive_underflow_is_not_reported_as_zero(): + expected_log = 10000 * math.log(0.75) + assert math.isclose(binomial_log_lower_tail(10000, 0), expected_log, rel_tol=1e-14) + assert math.isfinite(adaptive_sequential_log_lower_tail_bound(10000, 0)) + with pytest.raises(FloatingPointError, match="log_lower_tail"): + binomial_lower_tail(10000, 0) + with pytest.raises(FloatingPointError): + adaptive_sequential_lower_tail_bound(10000, 0) + # Exact-zero endpoints remain distinguishable from numerical underflow. + assert binomial_lower_tail(10, 0, 1.0) == 0.0 + assert binomial_log_lower_tail(10, 0, 1.0) == -math.inf + + +def test_monotonicity_symmetry_and_extreme_probabilities(): + for p in (1e-300, 1e-12, 0.25, 0.5, 0.9, math.nextafter(1.0, 0.0)): + values = [binomial_log_lower_tail(200, k, p) for k in range(201)] + assert all(a <= b + 1e-12 for a, b in zip(values, values[1:])) + for n in (101, 501, 2001): + assert math.isclose(binomial_lower_tail(n, n//2, 0.5), 0.5, rel_tol=2e-12) + # Preserve a tiny but representable near-one log-CDF via the complement. + assert math.isclose(binomial_log_lower_tail(1, 0, 1e-300), -1e-300, rel_tol=1e-14) + + +@pytest.mark.parametrize("n,k,p", [(-1, 0, .25), (1, -1, .25), (2, 3, .25), + (10, 1, float("nan")), (10, 1, float("inf")), (10, 1, -.1), (10, 1, 1.1)]) +def test_tail_rejects_invalid_parameters(n, k, p): + with pytest.raises(ValueError): + binomial_log_lower_tail(n, k, p) + + +@pytest.mark.parametrize("n,k", [(True, 0), (3, False), (3.5, 1), (3, 1.0)]) +def test_tail_requires_integer_counts(n, k): + with pytest.raises(TypeError): + binomial_lower_tail(n, k)