|
| 1 | +#!/usr/bin/env python3 |
| 2 | +"""Reproducible-features demo: train a model on `rosalind features` output and |
| 3 | +prove the inputs (and therefore the trained model) are bit-reproducible. |
| 4 | +
|
| 5 | +The headline is REPRODUCIBILITY, not biological novelty: identical inputs -> |
| 6 | +byte-identical features -> bit-identical trained weights, verifiable by the |
| 7 | +BLAKE3 receipt. Runs with numpy only (no pandas/sklearn/pyarrow needed). |
| 8 | +
|
| 9 | + python examples/reproducible_features_demo.py [path-to-rosalind-binary] |
| 10 | +""" |
| 11 | + |
| 12 | +from __future__ import annotations |
| 13 | + |
| 14 | +import hashlib |
| 15 | +import json |
| 16 | +import os |
| 17 | +import shutil |
| 18 | +import subprocess |
| 19 | +import sys |
| 20 | +import tempfile |
| 21 | + |
| 22 | +import numpy as np |
| 23 | + |
| 24 | +REPO = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) |
| 25 | +sys.path.insert(0, os.path.join(REPO, "python")) |
| 26 | +from rosalind import features # noqa: E402 |
| 27 | + |
| 28 | + |
| 29 | +def find_binary(argv: list[str]) -> str: |
| 30 | + if len(argv) > 1: |
| 31 | + return argv[1] |
| 32 | + for cand in (os.path.join(REPO, "target", "release", "rosalind"), |
| 33 | + os.path.join(REPO, "target", "debug", "rosalind")): |
| 34 | + if os.path.exists(cand): |
| 35 | + return cand |
| 36 | + found = shutil.which("rosalind") |
| 37 | + if found: |
| 38 | + return found |
| 39 | + sys.exit("rosalind binary not found — pass it as argv[1] or `cargo build --release`") |
| 40 | + |
| 41 | + |
| 42 | +def sh(cmd: list[str]) -> None: |
| 43 | + subprocess.run(cmd, check=True) |
| 44 | + |
| 45 | + |
| 46 | +def output_blake3(manifest_path: str) -> str: |
| 47 | + """Read the recorded BLAKE3 of the feature TSV from the run receipt.""" |
| 48 | + with open(manifest_path) as fh: |
| 49 | + m = json.load(fh) |
| 50 | + return m["outputs"][0]["blake3"] |
| 51 | + |
| 52 | + |
| 53 | +def train_logreg(x: np.ndarray, y: np.ndarray, iters: int = 300, lr: float = 0.2): |
| 54 | + """Deterministic pure-numpy logistic regression (zero init, fixed schedule).""" |
| 55 | + mu = x.mean(0) |
| 56 | + sd = x.std(0) |
| 57 | + sd[sd == 0.0] = 1.0 |
| 58 | + xs = (x - mu) / sd |
| 59 | + xb = np.hstack([np.ones((xs.shape[0], 1)), xs]) |
| 60 | + w = np.zeros(xb.shape[1]) |
| 61 | + n = xb.shape[0] |
| 62 | + for _ in range(iters): |
| 63 | + p = 1.0 / (1.0 + np.exp(-(xb @ w))) |
| 64 | + w -= lr * (xb.T @ (p - y)) / n |
| 65 | + return w, mu, sd |
| 66 | + |
| 67 | + |
| 68 | +def predict(w, mu, sd, x: np.ndarray) -> np.ndarray: |
| 69 | + sd = sd.copy() |
| 70 | + sd[sd == 0.0] = 1.0 |
| 71 | + xs = (x - mu) / sd |
| 72 | + xb = np.hstack([np.ones((xs.shape[0], 1)), xs]) |
| 73 | + return (1.0 / (1.0 + np.exp(-(xb @ w)))) >= 0.5 |
| 74 | + |
| 75 | + |
| 76 | +def main() -> int: |
| 77 | + binary = find_binary(sys.argv) |
| 78 | + work = tempfile.mkdtemp(prefix="rosalind-demo-") |
| 79 | + data = os.path.join(work, "toy") |
| 80 | + |
| 81 | + print("== building a toy dataset + index ==") |
| 82 | + sh([sys.executable, os.path.join(REPO, "scripts", "generate_toy_data.py"), data]) |
| 83 | + ref = os.path.join(data, "reference.fa") |
| 84 | + reads = os.path.join(data, "reads_R1.fastq") |
| 85 | + raw = os.path.join(work, "raw.bam") |
| 86 | + bam = os.path.join(work, "sorted.bam") |
| 87 | + idx = os.path.join(work, "ref.idx") |
| 88 | + sh([binary, "align", "--reference", ref, "--reads", reads, "--format", "bam", "--output", raw]) |
| 89 | + sh([binary, "sort", "--input", raw, "--output", bam]) |
| 90 | + sh([binary, "index", "--reference", ref, "--output", idx]) |
| 91 | + |
| 92 | + print("== extracting features TWICE (independent runs) ==") |
| 93 | + a = features(idx, bam, binary=binary, workdir=os.path.join(work, "a")) |
| 94 | + b = features(idx, bam, binary=binary, workdir=os.path.join(work, "b")) |
| 95 | + |
| 96 | + # --- Reproducibility proof ------------------------------------------------- |
| 97 | + tsv_a = os.path.join(work, "a", "features.tsv") |
| 98 | + tsv_b = os.path.join(work, "b", "features.tsv") |
| 99 | + bytes_identical = open(tsv_a, "rb").read() == open(tsv_b, "rb").read() |
| 100 | + h_a, h_b = output_blake3(a.manifest_path), output_blake3(b.manifest_path) |
| 101 | + local_hash = hashlib.sha256(open(tsv_a, "rb").read()).hexdigest()[:16] |
| 102 | + |
| 103 | + print(f" rows: {len(a)}; numeric features: {a.data.shape[1]}") |
| 104 | + print(f" features byte-identical across runs: {bytes_identical}") |
| 105 | + print(f" receipt BLAKE3 (run A): {h_a}") |
| 106 | + print(f" receipt BLAKE3 (run B): {h_b}") |
| 107 | + print(f" receipts match: {h_a == h_b}") |
| 108 | + assert bytes_identical, "feature TSVs differ across runs" |
| 109 | + assert h_a == h_b, "receipt hashes differ across runs" |
| 110 | + |
| 111 | + # --- A genuine supervised task: is the reference base a purine (A/G)? ------- |
| 112 | + # Read counts peak at the ref base, so this is learnable from the features. |
| 113 | + def label(ft) -> np.ndarray: |
| 114 | + return np.isin(ft.ref, [b"A", b"G"]).astype(np.float64) |
| 115 | + |
| 116 | + # Deterministic train/test split by position parity. |
| 117 | + def split(ft): |
| 118 | + test = (ft.pos % 2) == 1 |
| 119 | + return ~test, test |
| 120 | + |
| 121 | + ya, yb = label(a), label(b) |
| 122 | + tr_a, te_a = split(a) |
| 123 | + w_a, mu_a, sd_a = train_logreg(a.data[tr_a], ya[tr_a]) |
| 124 | + acc_a = float((predict(w_a, mu_a, sd_a, a.data[te_a]) == (ya[te_a] >= 0.5)).mean()) |
| 125 | + |
| 126 | + w_b, mu_b, sd_b = train_logreg(b.data[split(b)[0]], yb[split(b)[0]]) |
| 127 | + pred_a = predict(w_a, mu_a, sd_a, a.data[te_a]) |
| 128 | + pred_b = predict(w_b, mu_b, sd_b, b.data[split(b)[1]]) |
| 129 | + |
| 130 | + weights_identical = np.array_equal(w_a, w_b) |
| 131 | + preds_identical = np.array_equal(pred_a, pred_b) |
| 132 | + |
| 133 | + print("== model: predict is-purine(ref) from pileup features (logistic regression) ==") |
| 134 | + print(f" held-out accuracy: {acc_a:.4f}") |
| 135 | + print(f" trained weights bit-identical across the two extractions: {weights_identical}") |
| 136 | + print(f" test predictions bit-identical across the two extractions: {preds_identical}") |
| 137 | + assert weights_identical, "trained weights differ across extractions (not reproducible)" |
| 138 | + assert preds_identical, "predictions differ across extractions" |
| 139 | + |
| 140 | + print("\nPROOF: byte-identical features (hash %s) -> bit-identical trained model." |
| 141 | + " Rosalind feature extraction yields bit-reproducible ML training inputs." % local_hash) |
| 142 | + shutil.rmtree(work, ignore_errors=True) |
| 143 | + return 0 |
| 144 | + |
| 145 | + |
| 146 | +if __name__ == "__main__": |
| 147 | + raise SystemExit(main()) |
0 commit comments