From 858d6d2b9512b02a635ff64324156cee0d556429 Mon Sep 17 00:00:00 2001 From: Jeff Jaureguy <67065808+Jaureguy760@users.noreply.github.com> Date: Wed, 14 Oct 2026 18:16:28 -0700 Subject: [PATCH] chore: remove bulk parity investigation artifacts --- .github/workflows/ci.yml | 4 - README.md | 10 +- docs/BULK_PARITY_VALIDATION.md | 118 --------------------- scripts/compare_archive_analysis.py | 157 +--------------------------- scripts/diagnose_bulk_numerics.py | 144 ------------------------- 5 files changed, 2 insertions(+), 431 deletions(-) delete mode 100644 docs/BULK_PARITY_VALIDATION.md delete mode 100644 scripts/diagnose_bulk_numerics.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1e54ae2..6277af5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -120,10 +120,6 @@ jobs: if: matrix.python-version == '3.11' run: uv pip install --python .venv/bin/python numpy==1.26.4 scipy==1.13.1 - - name: Diagnose platform numerical arithmetic - if: matrix.python-version == '3.11' - run: .venv/bin/python scripts/diagnose_bulk_numerics.py - - name: Run tests run: | source .venv/bin/activate diff --git a/README.md b/README.md index 326e04d..a42e994 100644 --- a/README.md +++ b/README.md @@ -143,12 +143,6 @@ has equal reference and alternate counts has a symmetric likelihood: `mu` and `1 - mu` are equally likely. Numerical libraries may select opposite optima; the reported direction is not identifiable in that case. -Archive validation reports strict numerical matching separately from statistical -equivalence. A mirrored unphased proportion qualifies only with a balanced first -retained SNP, complementary proportions, independently verified equal archived -likelihoods, and passing counts, dispersion and all other statistical checks. -Reports retain the raw proportion differences and strict failures. - Feature analysis preserves one SNV in every overlapping feature. `peak` remains an alias for `feature`: @@ -181,9 +175,7 @@ truly phased heterozygous SNPs; otherwise omit `--phased`. Valid interior dispersion values are retained without clipping. Exact rho endpoints and probability-product underflow use documented finite extensions -because archived Python is undefined there. The full public ATAC comparison -and its reference adaptations are described in -[the validation report](docs/BULK_PARITY_VALIDATION.md). +because archived Python is undefined there. ## Single-Cell Example diff --git a/docs/BULK_PARITY_VALIDATION.md b/docs/BULK_PARITY_VALIDATION.md deleted file mode 100644 index 4aecb8c6..0000000 --- a/docs/BULK_PARITY_VALIDATION.md +++ /dev/null @@ -1,118 +0,0 @@ -# Bulk Python/Rust compatibility validation - -Validation date: 2026-10-14. Changes are proposed in PR #133, which includes -the compatibility changes from PR #135. Neither PR is merged or released. - -## Reference and data - -The reference is archived **Python**, commit -`41354fc087404ba1fa99542803ee9cb9b7e0e6bf`, with NumPy 1.26.4 and SciPy 1.13.1. -Tests execute its functions; Rust outputs are not used as expected results. - -Public real ATAC inputs are GM12878 experiment ENCSR637XSC, replicate-2 -BAM ENCFF962FMH, peaks ENCFF452ZBX, and matched HG001 GIAB NISTv3.3.1 -GRCh38 genotypes. Source MD5 checks passed. Chromosome 21 contains 1,760,836 -alignment records, 27,182 eligible SNPs and 556 peak-overlapping SNP rows. -Counting matches the unmodified archived function at thread settings 1 and 4. -The full counting CLI independently reproduced the 556 peak rows and memberships. -Source links, compressed count tables and checksums are in -[the public-data fixture](../tests/analysis/fixtures/public_atac/README.md). - -## Results and reproduction - -- The 32-run matrix passes locally on macOS arm64: 16 legacy runs across - SNP/peak, single/linear, phased/unphased and pseudocount 0/1, plus 16 donor - feature runs across global/per-donor dispersion. Reported raw mu is compared. -- Counts require exact agreement. Statistical tolerances remain - `rtol=atol=5e-5`; no proportion flipping or region-specific exceptions occur. -- Every legacy matrix case also executes CLI inference and checks its output. - Mixed-genotype tests compare the actual archive's successful and error paths. -- Heterogeneous donor tests independently fit dispersion with Python and compare - donor effects, chi-square LRT p-values and donor-local BH correction. They cover - overlapping features, both dispersion models, phase modes and pseudocounts. -- Local suite: 474 Python tests passed, 82 skipped, 1 deselected; 132 Rust tests - passed, 7 ignored. Ruff, formatting and diff checks pass. Rust 1.88 Clippy - passes with the pre-existing `uninlined_format_args` lint disabled. - -In a full clone with the built extension and pinned numerical dependencies: - -```bash -python -m pytest tests/analysis/test_public_atac_matrix.py \ - tests/analysis/test_bulk_edge_parity.py \ - tests/analysis/test_donor_archive_parity.py -q -``` - -Python 3.11 Linux CI pins the numerical dependencies and runs the complete real -matrix. Other Python versions run the remaining compatibility regressions. -Final-head Linux and container outcomes must be reviewed separately in GitHub -Actions; local results do not establish those outcomes. - -The first full Linux run at `cb7b834` exposed six dynamic-reference failures: -three public unphased peak cases differed in raw mu, and three heterogeneous -per-donor linear cases differed in fitted statistics. Earlier frozen regressions -passed, so they did not establish full Linux compatibility. Both native amd64 -and arm64 container builds, ATAC smoke tests and label checks passed on that head. -These strict failures remain integration gates while platform arithmetic is -investigated. `scripts/diagnose_bulk_numerics.py` compares fused/unfused kernels -and NumPy scalar/vector arithmetic without modifying production source. -The probe confirms that the pinned Linux x86_64 SciPy build uses separate -multiply/add arithmetic at the tested Cephes components, while the macOS arm64 -build uses fused arithmetic. Rust now follows that measured compiler policy; -the full strict matrix must verify the resulting inference on each platform. -At `9d3006d`, all 32 executed-reference Linux matrix runs and heterogeneous donor -comparisons pass. At `625252d`, Linux Python 3.11 passes 464 tests, 76 skipped, -1 deselected, including all 16 native conditional checks; Python 3.10 and 3.12 -also pass. Independently executed archived Linux Python returns -`mu=0.6928487541708976` for the fixed symmetric linear case where the Mac -reference returns `0.3071512458291038`. The Rust trace test records each measured -platform's raw outcome, without accepting complementary outputs. -Conditional tests execute the native archived oracle at identical supplied -dispersion parameters; the original Mac snapshot is unchanged and also checked -in its pinned environment. Matching version numbers alone does not establish -identical compiler/library arithmetic across platforms. - -## Deliberate extensions and reference repairs - -Global/per-donor dispersion is the intended extension. Shared dispersion does -not pool effect tests or multiple-testing correction across donors. Archived -linear array shapes require explicit reference adapters; phased linear tests -explicitly forward the phase flag that the archived linear wrapper ignored. -The shape adapter preserves the original sequential probability recurrence; -exact scalar/vector recurrence tests guard against replacing it with a -mathematically equivalent but differently rounded product reduction. - -Interior rho values are preserved, including tested values from `1e-12` to -`1-1e-12`. Exact endpoints retain a finite regularization because archived -beta-binomial shapes are undefined there. Long probability products use -log-space fallback after losing normal range; this extends the archive's -underflowing arithmetic rather than reproducing its failure. - -Archived phasing checks use distinct original GT tokens and the presence of -`0|1`/`1|0`, then parse every retained first allele. Numeric homozygous tokens -follow that behavior. Malformed mixed tokens raise an error instead of silently -changing a region's likelihood. Numeric-GT duplicates are removed using the parsed first allele, while the -original tokens remain available for the backend's pre-dedup phase decision. -Sixteen extra regressions cover both models, both dispersion scopes, both -orientations and pseudocounts, including fitting weights, QC and phase provenance. -One public peak SNP has `0/1`: the unfiltered -phased command is rejected by both implementations. Phased comparisons explicitly -select truly phased heterozygous SNPs; unphased comparisons retain all SNPs. - -## Integration and release gates - -Review compatibility PR #135 first, then refresh #133 against the resulting dev -branch and rerun required checks. No merge or tag is included in this validation. -All active version metadata currently uses 1.5.0. CI/release image builds pass -the package version and source SHA to OCI labels; container smoke tests verify -those labels and Python/Rust package versions and exercise ATAC grouping/phasing. - -Review the exact final head's tests and containers before integration. A release -dry run remains necessary after the reviewed code reaches the permitted release -branch, including the separate Singularity preflight. Existing published images -and packages will remain old until a release is actually published. - -These comparisons cover one real donor/chromosome and controlled heterogeneous -donors. They do not establish biological calibration, a new WASP remapping pass, -or universal bitwise agreement across numerical environments. A symmetric -unphased likelihood has no identifiable biological direction even when its -reported mu agrees numerically with the pinned reference. diff --git a/scripts/compare_archive_analysis.py b/scripts/compare_archive_analysis.py index 7cfadec..0bafbbe 100644 --- a/scripts/compare_archive_analysis.py +++ b/scripts/compare_archive_analysis.py @@ -1,8 +1,7 @@ #!/usr/bin/env python3 """Compare analysis with Aaron Ho's pinned, pre-migration Python implementation. -Run from a full repository clone with an installed Rust extension: - python scripts/compare_archive_analysis.py --output comparison.json +Reference helpers for regression tests in a full repository clone. The archived code is read from git without altering it. AnnData is transposed at the boundary to compare the historical SNP-by-cell layout with today's cell-by-SNP layout. This is a compatibility check, not statistical calibration. @@ -10,15 +9,11 @@ from __future__ import annotations -import argparse import contextlib import io -import json import subprocess import sys -import tempfile import types -import warnings from collections.abc import Callable from pathlib import Path @@ -397,153 +392,3 @@ def compare_bulk_equivalence( "failures": failures, "regions": evidence, } - - -def run_comparison(archive: dict[str, types.ModuleType]) -> dict: - import wasp2_rust - from analysis.as_analysis_sc import get_imbalance_sc - from analysis.compare_ai import get_compared_imbalance, get_compared_imbalance_diff_snps - from analysis.run_analysis import run_ai_analysis - from analysis.run_analysis_sc import process_adata_inputs - - cases = {} - with tempfile.TemporaryDirectory(prefix="wasp2-parity-") as tmp: - temp = Path(tmp) - for model in ["single", "linear"]: - for unit in ["snv", "atac", "rna"]: - for phased in [False, True]: - for pc in [0, 1]: - name = f"bulk/{model}/{unit}/phased={phased}/pc={pc}" - region = {"snv": None, "atac": "region", "rna": "gene"}[unit] - counts = bulk_counts().drop( - columns=["region", "gene"] - if unit == "snv" - else ["gene"] - if unit == "atac" - else ["region"] - ) - counts.to_csv(temp / "counts.tsv", sep="\t", index=False) - try: - expected = archive["as_analysis"].get_imbalance( - counts.copy(), - method=model, - phased=phased, - region_col=region, - pseudocount=pc, - ) - except Exception as exc: - cases[name] = {"status": "REFERENCE_ERROR", "reason": repr(exc)} - continue - try: - run_ai_analysis( - temp / "counts.tsv", - model=model, - phased=phased, - region_col=region, - pseudocount=pc, - out_file=str(temp / "results.tsv"), - ) - actual = pd.read_csv(temp / "results.tsv", sep="\t") - if region != "region": - expected = expected.rename(columns={region or "variant": "region"}) - cases[name] = compare_frames(expected, actual) - except Exception as exc: - cases[name] = {"status": "ERROR", "reason": repr(exc)} - for phased in [False, True]: - for pc in [0, 1]: - old = single_cell_counts() - new = old.T.copy() - options = { - "sample": "sample1", - "phased": phased, - "pseudocount": pc, - "min_count": 10, - } - for label, reference, current in [ - ("single-cell", archive["as_analysis_sc"].get_imbalance_sc, get_imbalance_sc), - ( - "differential/shared", - archive["compare_ai"].get_compared_imbalance, - get_compared_imbalance, - ), - ( - "differential/different", - archive["compare_ai"].get_compared_imbalance_diff_snps, - get_compared_imbalance_diff_snps, - ), - ]: - name = f"{label}/phased={phased}/pc={pc}" - try: - expected = reference(old.copy(), **options) - except Exception as exc: - cases[name] = {"status": "REFERENCE_ERROR", "reason": repr(exc)} - continue - try: - prepared = process_adata_inputs(new.copy(), sample="sample1", phased=phased) - actual = current(prepared.adata, **options) - details = { - str(key): compare_frames(frame, actual[key]) - for key, frame in expected.items() - } - cases[name] = { - "status": "PASS" - if expected.keys() == actual.keys() - and all(d["status"] == "PASS" for d in details.values()) - else "FAIL", - "groups": details, - } - except Exception as exc: - cases[name] = {"status": "ERROR", "reason": repr(exc)} - return { - "archive_commit": subprocess.check_output( - ["git", "rev-parse", ARCHIVE_REF], cwd=ROOT, text=True - ).strip(), - "current_commit": subprocess.check_output( - ["git", "rev-parse", "HEAD"], cwd=ROOT, text=True - ).strip(), - "extension": wasp2_rust.__file__, - "cases": cases, - } - - -def main() -> None: - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument("--output", type=Path, required=True) - parser.add_argument( - "--adapt-linear-reference", - action="store_true", - help="Explicitly repair archived scalar/array shapes; report this adaptation in JSON.", - ) - args = parser.parse_args() - # Import the installed extension before src's stub package can shadow it. - import wasp2_rust # noqa: F401 - - sys.path.insert(0, str(ROOT / "src")) - with warnings.catch_warnings(), contextlib.redirect_stdout(io.StringIO()): - warnings.simplefilter("ignore", FutureWarning) - archive = load_archive() - if args.adapt_linear_reference: - adapt_linear_reference(archive) - report = run_comparison(archive) - report["reference_adaptations"] = ( - ["Scalarize one-element objectives and align per-SNP dispersion arrays"] - if args.adapt_linear_reference - else [] - ) - report["known_reference_limitations"] = [ - "Archived linear_model does not forward phased to parse_opt; phased multi-SNP " - "linear comparisons therefore test different declared behavior." - ] - args.output.write_text(json.dumps(report, indent=2) + "\n") - summary = { - status: sum(c["status"] == status for c in report["cases"].values()) - for status in ["PASS", "FAIL", "ERROR", "REFERENCE_ERROR"] - } - print(json.dumps(summary)) - for name, result in report["cases"].items(): - if result["status"] != "PASS": - print(name, result.get("failures", result.get("reason", result.get("groups")))) - - -if __name__ == "__main__": - main() diff --git a/scripts/diagnose_bulk_numerics.py b/scripts/diagnose_bulk_numerics.py deleted file mode 100644 index 9b3c7c6..0000000 --- a/scripts/diagnose_bulk_numerics.py +++ /dev/null @@ -1,144 +0,0 @@ -"""Compare fused/unfused Rust Cephes arithmetic with the installed SciPy. - -Builds two isolated probe modules from the production source in a temporary -crate. The repository's source and numerical policy are never modified. -""" - -from __future__ import annotations - -import csv -import hashlib -import io -import json -import math -import os -import platform -import subprocess -import sys -import tempfile -from pathlib import Path - -import numpy as np -import scipy -from scipy.special import gamma, gammaln -from scipy.stats import betabinom - -ROOT = Path(__file__).resolve().parents[1] - - -def main() -> None: - source = (ROOT / "rust/src/scipy_compat.rs").read_text() - start = source.index("#[inline]\nfn cephes_mul_add") - end = source.index("\n}\n", start) + 3 - implementation = source[start:end] - signature = "#[inline]\nfn cephes_mul_add(a: f64, b: f64, c: f64) -> f64" - fused = source.replace(implementation, signature + " { a.mul_add(b, c) }\n") - plain = source.replace(implementation, signature + " { a * b + c }\n") - wrapper = "\npub(crate) fn gamma_probe(x: f64) -> (f64, f64) { (gamma_positive(x), lgamma_positive(x)) }\n" - rust_main = r""" -mod fused; -mod plain; -fn main() { - for x in [0.1, 9.548148325980558, 16.97709256606237, 1000.0, 1e8, 5e11] { - let f = fused::gamma_probe(x); - let p = plain::gamma_probe(x); - println!("gamma,{x:?},{:?},{:?},{:?},{:?}", f.0, p.0, f.1, p.1); - } - for rho in [1e-12, 1e-9, 0.01, 0.1, 0.99] { - for k in [0, 10, 20] { - let a = 0.5 * (1.0-rho)/rho; - let f = fused::betabinom_logpmf(k, 20, a, a).unwrap(); - let p = plain::betabinom_logpmf(k, 20, a, a).unwrap(); - println!("pmf,{k},20,{a:?},{a:?},{f:?},{p:?}"); - } - } - let (a,b) = (8.977092566062371,5.5481483259805575); - println!("pmf,8,12,{a:?},{b:?},{:?},{:?}", - fused::betabinom_logpmf(8,12,a,b).unwrap(), - plain::betabinom_logpmf(8,12,a,b).unwrap()); - // Keep both reduction functions live when compiling with warning checks. - assert_eq!(fused::numpy_sum(&[1.0,2.0]),plain::numpy_sum(&[1.0,2.0])); -} -""" - with tempfile.TemporaryDirectory(prefix="wasp2-numerical-probe-") as directory: - temp = Path(directory) - (temp / "src").mkdir() - (temp / "Cargo.toml").write_text( - '[package]\nname="wasp2-numerical-probe"\nversion="0.0.0"\nedition="2021"\n' - '[dependencies]\nanyhow="1"\n' - ) - (temp / "src/fused.rs").write_text(fused + wrapper) - (temp / "src/plain.rs").write_text(plain + wrapper) - (temp / "src/main.rs").write_text(rust_main) - env = dict(os.environ) - env["CARGO_TARGET_DIR"] = str(ROOT / "rust/target/numerical-probe") - completed = subprocess.run( - ["cargo", "run", "--release", "--manifest-path", str(temp / "Cargo.toml")], - env=env, - capture_output=True, - text=True, - check=True, - ) - rows = [] - logs = [] - for row in csv.reader(io.StringIO(completed.stdout)): - if row[0] == "gamma": - x, fg, pg, fl, pl = map(float, row[1:]) - rows.append( - { - "kind": "gamma", - "x": x, - "python_gamma": float(gamma(x)), - "fused_gamma": fg, - "plain_gamma": pg, - "python_lgamma": float(gammaln(x)), - "fused_lgamma": fl, - "plain_lgamma": pl, - } - ) - else: - k, n, a, b, fused, plain_value = map(float, row[1:]) - reference = float(betabinom.logpmf(k, n, a, b)) - logs.append(reference) - rows.append( - { - "kind": "logpmf", - "k": k, - "n": n, - "alpha": a, - "beta": b, - "python": reference, - "fused": fused, - "plain": plain_value, - } - ) - values = np.array(logs * 32) - scalar_exp = np.array([math.exp(x) for x in values]) - vector_exp = np.exp(values) - numpy_so = Path(np.core._multiarray_umath.__file__) - report = { - "platform": platform.platform(), - "numpy": np.__version__, - "scipy": scipy.__version__, - "numpy_binary_sha256": hashlib.sha256(numpy_so.read_bytes()).hexdigest(), - "numpy_cpu_features": np.core._multiarray_umath.__cpu_features__, - "scalar_vector_exp_differences": int(np.count_nonzero(scalar_exp != vector_exp)), - "scalar_vector_exp_max_delta": float(np.max(np.abs(scalar_exp - vector_exp))), - "rows": rows, - } - print(json.dumps(report, indent=2)) - sys.path.insert(0, str(ROOT)) - from scripts.compare_archive_analysis import fixed_bulk_reference - - reference = json.loads( - (ROOT / "tests/analysis/fixtures/real_atac_optimizer_reference.json").read_text() - ) - # These are newly executed Python outcomes, not adjusted Rust expectations. - fixed_cases = { - case["name"]: fixed_bulk_reference(case).to_dict("records") for case in reference["cases"] - } - print("ARCHIVED_FIXED_NUISANCE_REFERENCES=" + json.dumps(fixed_cases)) - - -if __name__ == "__main__": - main()