#!/usr/bin/env python3 """Build a public aggregate audit of Paper 3 same-record HRV linkability. The audit intentionally uses a simpler and stricter protocol than the earlier internal headline analysis: * Autonomic Aging subjects are fixed by the frozen held-out split. * Each person must have at least two non-overlapping 300-beat windows. * Window 0 is the enrollment template and window 1 is the query. * The query is never included in its own template (no self-contamination). * All 30 features use standardization parameters frozen on training data. * Only aggregate counts and metrics are written; subject identifiers and subject-level ranks are never published. This evaluates closed-set linkage between adjacent windows from the same recording session. It does not establish longitudinal identity, open-set authentication, Apple Watch performance, or real-world re-identification. """ from __future__ import annotations import argparse import hashlib import json from pathlib import Path import numpy as np import pandas as pd COHORT = "autonomic_aging" BOOTSTRAP_SEED = 20260825 PERMUTATION_SEED = 20260825 DEFAULT_REPLICATES = 20_000 TOP_K = (1, 5, 10) def sha256(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as source: for block in iter(lambda: source.read(1024 * 1024), b""): digest.update(block) return digest.hexdigest() def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--features", type=Path, required=True) parser.add_argument("--splits", type=Path, required=True) parser.add_argument("--standardization", type=Path, required=True) parser.add_argument("--feature-extractor", type=Path, required=True) parser.add_argument("--feature-definition", type=Path, required=True) parser.add_argument("--peak-correction", type=Path, required=True) parser.add_argument("--output", type=Path, required=True) parser.add_argument("--replicates", type=int, default=DEFAULT_REPLICATES) return parser.parse_args() def standardize(rows: pd.DataFrame, feature_names: list[str], params: dict) -> np.ndarray: matrix = rows[feature_names].to_numpy(dtype=np.float64) for index, name in enumerate(feature_names): feature_params = params[name] scale = float( feature_params.get( "scale", 1.4826 * float(feature_params.get("mad", 0.0)), ) ) if not np.isfinite(scale) or scale <= 0: raise ValueError(f"Invalid frozen scale for {name}: {scale}") matrix[:, index] = ( matrix[:, index] - float(feature_params["median"]) ) / scale matrix = np.where(np.isfinite(matrix), matrix, 0.0) norms = np.linalg.norm(matrix, axis=1, keepdims=True) if np.any(norms <= 0): raise ValueError("A standardized feature vector has zero norm") return matrix / norms def rank_hits(enrollment: np.ndarray, query: np.ndarray) -> dict[int, np.ndarray]: similarities = query @ enrollment.T ranks = np.argsort(-similarities, axis=1) truth = np.arange(len(enrollment)) return { k: (ranks[:, :k] == truth[:, None]).any(axis=1) for k in TOP_K } def fixed_candidate_bootstrap(hits: np.ndarray, replicates: int) -> dict: """Resample query subjects while keeping the candidate set unchanged. This avoids the invalid pseudo-ID construction in the superseded Paper 3 bootstrap, where repeated draws of the same person were made to compete as if they were different identities. """ rng = np.random.default_rng(BOOTSTRAP_SEED) n = len(hits) draws = rng.integers(0, n, size=(replicates, n)) estimates = hits[draws].mean(axis=1) low, median, high = np.quantile(estimates, [0.025, 0.5, 0.975]) return { "method": "fixed-candidate query-subject bootstrap", "replicates": replicates, "seed": BOOTSTRAP_SEED, "resampling_unit": "one query subject with one query window", "candidate_set_recomputed": False, "duplicate_draws_do_not_create_new_candidate_ids": True, "median": round(float(median), 6), "ci_95": [round(float(low), 6), round(float(high), 6)], "interpretation": ( "Conditional uncertainty across query subjects for this fixed " "202-person candidate set; not a population or deployment CI." ), } def wilson_interval(hit_count: int, n: int) -> list[float]: """Two-sided 95% Wilson score interval for a binomial proportion.""" z = 1.959963984540054 proportion = hit_count / n denominator = 1.0 + z * z / n centre = (proportion + z * z / (2.0 * n)) / denominator margin = ( z * np.sqrt(proportion * (1.0 - proportion) / n + z * z / (4.0 * n * n)) / denominator ) return [round(float(centre - margin), 6), round(float(centre + margin), 6)] def permutation_test(hits: dict[int, np.ndarray], enrollment: np.ndarray, query: np.ndarray, replicates: int) -> dict: """Test the no-linkage null by permuting query-to-identity assignments.""" similarities = query @ enrollment.T ranks = np.argsort(-similarities, axis=1) n = len(enrollment) observed = {k: float(hits[k].mean()) for k in TOP_K} exceedances = {k: 0 for k in TOP_K} rng = np.random.default_rng(PERMUTATION_SEED) truth = np.arange(n) for _ in range(replicates): permuted_truth = rng.permutation(truth) for k in TOP_K: null_accuracy = float( (ranks[:, :k] == permuted_truth[:, None]).any(axis=1).mean() ) exceedances[k] += int(null_accuracy >= observed[k]) return { "method": "bijection permutation of query identity labels", "replicates": replicates, "seed": PERMUTATION_SEED, "p_one_sided": { str(k): round((exceedances[k] + 1) / (replicates + 1), 8) for k in TOP_K }, } def main() -> int: args = parse_args() if args.replicates < 1_000: raise SystemExit("Use at least 1,000 bootstrap/permutation replicates") features = pd.read_parquet(args.features) split_document = json.loads(args.splits.read_text(encoding="utf-8")) standardization_document = json.loads( args.standardization.read_text(encoding="utf-8") ) params = standardization_document["params"] feature_names = list(params) if len(feature_names) != 30: raise ValueError(f"Expected 30 frozen features, found {len(feature_names)}") held_out_ids = split_document["splits"][COHORT]["test"] cohort_rows = features.loc[ (features["dataset"] == COHORT) & features["subject_id"].isin(held_out_ids) ].copy() if cohort_rows["subject_id"].nunique() != len(held_out_ids): raise ValueError( "The feature artifact does not contain every frozen held-out subject: " f"{cohort_rows['subject_id'].nunique()} of {len(held_out_ids)}" ) if "subject_low_sqi" in cohort_rows and cohort_rows["subject_low_sqi"].any(): raise ValueError("Held-out rows marked low-SQI must be resolved before audit") cohort_rows = cohort_rows.sort_values(["subject_id", "window_idx"]) window_counts = cohort_rows.groupby("subject_id").size() eligible_ids = sorted(window_counts.loc[window_counts >= 2].index) excluded_one_window = int((window_counts < 2).sum()) enrollment_rows = [] query_rows = [] for subject_id in eligible_ids: subject_rows = cohort_rows.loc[ cohort_rows["subject_id"] == subject_id ].sort_values("window_idx") enrollment_rows.append(subject_rows.iloc[0]) query_rows.append(subject_rows.iloc[1]) enrollment_frame = pd.DataFrame(enrollment_rows) query_frame = pd.DataFrame(query_rows) if not ( enrollment_frame["subject_id"].tolist() == query_frame["subject_id"].tolist() == eligible_ids ): raise AssertionError("Enrollment/query subject order is not identical") if not np.all(query_frame["window_idx"].to_numpy() > enrollment_frame["window_idx"].to_numpy()): raise AssertionError("Query windows must follow enrollment windows") selected_feature_values = pd.concat( [enrollment_frame[feature_names], query_frame[feature_names]], ignore_index=True, ).to_numpy(dtype=np.float64) nonfinite_feature_values = int((~np.isfinite(selected_feature_values)).sum()) enrollment = standardize(enrollment_frame, feature_names, params) query = standardize(query_frame, feature_names, params) forward_hits = rank_hits(enrollment, query) reverse_hits = rank_hits(query, enrollment) n = len(eligible_ids) metrics = {} for k in TOP_K: hit_count = int(forward_hits[k].sum()) metrics[f"top_{k}"] = { "hits": hit_count, "queries": n, "accuracy": round(hit_count / n, 6), "random_rank_chance": round(min(k / n, 1.0), 6), "wilson_95": wilson_interval(hit_count, n), "fixed_candidate_bootstrap": fixed_candidate_bootstrap( forward_hits[k], args.replicates ), } result = { "schema_version": "1.0.0", "analysis_date": "2026-08-25", "status": "unreviewed internal secondary analysis", "analysis_script_sha256": sha256(Path(__file__).resolve()), "input_sha256": { "features_pooled_parquet": sha256(args.features), "frozen_splits_json": sha256(args.splits), "training_frozen_standardization_json": sha256(args.standardization), "feature_extraction_script": sha256(args.feature_extractor), "feature_definition_script": sha256(args.feature_definition), "peak_correction_script": sha256(args.peak_correction), }, "cohort": COHORT, "frozen_held_out_subjects": len(held_out_ids), "eligible_subjects": n, "excluded_subjects_with_fewer_than_two_windows": excluded_one_window, "published_subject_identifiers": False, "features": { "count": len(feature_names), "names": feature_names, "standardization": "median and robust scale frozen on training data", "similarity": "cosine", "nonfinite_values_replaced_with_zero_after_standardization": nonfinite_feature_values, }, "upstream_feature_generation": { "ecg_channel": ( "ECG1; this is lead I in two-channel Task Force Monitor records " "and is a limitation of the existing Paper 3 feature artifact" ), "rr_physiological_gate_ms": [300, 2000], "peak_correction": "Lipponen-Tarvainen/Kubios before windowing", "windowing": "consecutive non-overlapping 300-beat windows", }, "protocol": { "candidate_set": f"fixed closed set of {n} held-out subjects", "enrollment": "first retained non-overlapping 300-beat window", "query": "immediately following retained non-overlapping 300-beat window", "windows_per_subject_used": 2, "query_in_enrollment_template": False, "model_fit_on_held_out_subjects": False, "one_query_per_subject": True, }, "metrics": metrics, "permutation_test": permutation_test( forward_hits, enrollment, query, args.replicates ), "direction_sensitivity": { "protocol": "swap the first and second windows", **{ f"top_{k}_accuracy": round(float(reverse_hits[k].mean()), 6) for k in TOP_K }, }, "supersedes_for_public_reporting": { "value": "top-5 81.10% on 211 subjects", "reason": ( "The earlier number is retained in the private audit trail but " "its current public-safe artifacts cannot reproduce the full " "evaluation chain. Code audit also found an incorrect " "leave-one-out centroid approximation in the point-estimate " "path and invalid pseudo-identity competition in its " "uncertainty procedure." ), }, "limits": [ "All enrollment and query windows come from the same recording session.", "This tests closed-set linkage, not open-set authentication or identity proof.", "It does not establish stability across days, devices, sleep/wake state, or Apple Watch.", "Shared session, sensor, preprocessing, or signal artifacts may contribute to linkage.", "The fixed-candidate bootstrap interval is not a population-performance interval.", ], } args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text( json.dumps(result, ensure_ascii=False, indent=2) + "\n", encoding="utf-8", ) print(json.dumps(result, ensure_ascii=False, indent=2)) return 0 if __name__ == "__main__": raise SystemExit(main())