#!/usr/bin/env python3
"""Build the public, aggregate age-by-sex HRV reference used by Feelmo docs.

The script reads the PhysioNet Autonomic Aging v1.0.0 WFDB records, analyzes
the first five minutes of resting ECG, and publishes aggregate statistics only.
Subject-level features stay in the research artifact directory because HRV can
carry identity information.

This is a research reference, not an Apple Watch normal range. Apple Watch
HealthKit SDNN samples and a standardized five-minute clinical ECG are not
interchangeable measurements.
"""

from __future__ import annotations

import argparse
import hashlib
import importlib.metadata
import json
import os
import shutil
import tempfile
from dataclasses import dataclass
from pathlib import Path

# Scientific plotting is imported transitively by NeuroKit. Keep its caches in
# writable temporary storage both in the main process and any worker.
_MPL_CACHE = Path("/private/tmp/feelmo-matplotlib-cache")
_XDG_CACHE = Path("/private/tmp/feelmo-xdg-cache")
_MPL_CACHE.mkdir(parents=True, exist_ok=True)
_XDG_CACHE.mkdir(parents=True, exist_ok=True)
os.environ.setdefault("MPLCONFIGDIR", str(_MPL_CACHE))
os.environ.setdefault("XDG_CACHE_HOME", str(_XDG_CACHE))

import numpy as np
import neurokit2 as nk
import pandas as pd
import wfdb
from joblib import Parallel, delayed


SCRIPT_PATH = Path(__file__).resolve()


def locate_docs_repo(script_path: Path) -> Path:
    """Locate the repository from either the source or published script copy."""
    for candidate in (script_path.parent, *script_path.parents):
        if (
            (candidate / "package.json").is_file()
            and (candidate / "docs" / ".vuepress").is_dir()
        ):
            return candidate
    return script_path.parent.parent


DOCS_REPO = locate_docs_repo(SCRIPT_PATH)
DEFAULT_PUBLIC_DIR = (
    DOCS_REPO / "docs" / ".vuepress" / "public" / "data" / "autonomic-aging"
)
DEFAULT_SUBJECT_OUTPUT = (
    DOCS_REPO.parent
    / "female-hrv-atlas"
    / "artifacts"
    / "autonomic-aging"
    / "feelmo_reference_features_5min.parquet"
)


DATASET_DOI = "https://doi.org/10.13026/2hsy-t491"
PAPER_DOI = "https://doi.org/10.1038/s41597-022-01202-y"
DATASET_VERSION = "1.0.0"
INDEPENDENT_QC_DOI = "https://doi.org/10.1038/s41598-023-40385-1"
SEX_LABELS = {0: "male", 1: "female"}
AGE_BANDS = (
    ("18–29", (1, 2, 3)),
    ("30–39", (4, 5)),
    ("40–49", (6, 7)),
    ("50–59", (8, 9)),
    ("60–69", (10, 11)),
    ("70–92", (12, 13, 14, 15)),
)
METRICS = ("SDNN", "RMSSD")
MIN_PUBLIC_CELL_N = 10
PUBLIC_SCRIPT_NAME = "build_autonomic_aging_reference.py"
COUNT_DISAGREEMENT_TOLERANCE = 0.10
QUALITY_CLASS_RANK = {
    "": -1,
    "Unacceptable": 0,
    "Barely acceptable": 1,
    "Excellent": 2,
}


@dataclass(frozen=True)
class Paths:
    data_dir: Path
    public_dir: Path
    subject_output: Path


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 package_version(name: str) -> str:
    try:
        return importlib.metadata.version(name)
    except importlib.metadata.PackageNotFoundError:
        return "not-installed"


def verify_source_checksums(data_dir: Path) -> dict:
    """Verify the release-provided SHA-256 manifest before analysis."""
    checksum_path = data_dir / "SHA256SUMS.txt"
    checked = 0
    failures: list[str] = []
    for raw_line in checksum_path.read_text(encoding="utf-8").splitlines():
        if not raw_line.strip():
            continue
        expected, relative_name = raw_line.split(maxsplit=1)
        source_path = data_dir / relative_name
        if not source_path.is_file() or sha256(source_path) != expected:
            failures.append(relative_name)
        checked += 1
    if failures:
        sample = ", ".join(failures[:5])
        raise RuntimeError(
            f"Source checksum verification failed for {len(failures)} file(s): {sample}"
        )
    return {
        "performed": True,
        "listed_files": checked,
        "failed_files": 0,
        "checksum_manifest_sha256": sha256(checksum_path),
    }


def detect_rpeaks(
    ecg: np.ndarray,
    fs: int,
    detector_method: str = "neurokit",
) -> tuple[np.ndarray, np.ndarray]:
    """Detect R peaks with one pre-declared NeuroKit detector.

    The previous positive/negative polarity fallback selected whichever side
    produced more peaks. A large T wave can therefore win by only one or two
    detections while producing a clinically implausible RR series. NeuroKit's
    default detector uses the cleaned ECG morphology and avoids that arbitrary
    choice. Elgendi and Pan-Tompkins are retained as technical comparators.
    """
    supported = {"neurokit", "elgendi2010", "pantompkins1985"}
    if detector_method not in supported:
        raise ValueError(f"unsupported detector method: {detector_method}")
    cleaned = nk.ecg_clean(ecg, sampling_rate=fs, method=detector_method)
    _, info = nk.ecg_peaks(
        cleaned,
        sampling_rate=fs,
        method=detector_method,
        correct_artifacts=False,
    )
    return cleaned, np.asarray(info["ECG_R_Peaks"], dtype=np.int64)


def correct_rpeaks(rpeaks: np.ndarray, fs: int) -> tuple[dict, np.ndarray]:
    """Apply Lipponen–Tarvainen/Kubios correction at the native sample rate."""
    artifacts, corrected = nk.signal_fixpeaks(
        rpeaks,
        sampling_rate=fs,
        iterative=True,
        method="kubios",
    )
    return artifacts, np.asarray(corrected, dtype=np.int64)


def normalize_processing_error(error: Exception) -> str:
    """Map implementation-specific exceptions to stable public QC labels."""
    if isinstance(error, ValueError) and "reshape" in str(error).lower():
        return "peak_correction_failed"
    return "processing_error"


def extract_subject(
    row: dict,
    data_dir: Path,
    segment_seconds: int,
    ecg_channel_offset: int = 0,
    detector_method: str = "neurokit",
) -> dict:
    subject_id = str(row["ID"]).zfill(4)
    result = {
        "subject_id": subject_id,
        "age_group": row.get("Age_group"),
        "sex": row.get("Sex"),
        "device": row.get("Device"),
        "segment_seconds": segment_seconds,
        "detector_method": detector_method,
        "ecg_channel_index": np.nan,
        "ecg_channel_name": "",
        "ecg_channel_count": np.nan,
        "n_rpeaks": np.nan,
        "n_rr_clean": np.nan,
        "ecg_quality_class": "",
        "n_artifact_flags": np.nan,
        "artifact_flag_fraction": np.nan,
        "MEAN_RR": np.nan,
        "SDNN": np.nan,
        "RMSSD": np.nan,
        "qc_status": "error",
        "qc_reason": "unknown",
        "processing_error_type": "",
    }

    try:
        header = wfdb.rdheader(str(data_dir / subject_id))
        ecg_channels = [
            index for index, name in enumerate(header.sig_name)
            if "ecg" in name.lower()
        ]
        result["ecg_channel_count"] = len(ecg_channels)
        if not ecg_channels:
            result["qc_reason"] = "no_ecg_channel"
            return result
        if ecg_channel_offset >= len(ecg_channels):
            result["qc_reason"] = "requested_ecg_channel_unavailable"
            return result

        ecg_channel_index = ecg_channels[ecg_channel_offset]
        result["ecg_channel_index"] = ecg_channel_index
        result["ecg_channel_name"] = str(header.sig_name[ecg_channel_index])

        fs = int(header.fs)
        segment_samples = min(int(header.sig_len), fs * segment_seconds)
        record = wfdb.rdrecord(
            str(data_dir / subject_id),
            sampfrom=0,
            sampto=segment_samples,
            channels=[ecg_channel_index],
        )
        ecg = np.asarray(record.p_signal[:, 0], dtype=np.float64)
        if not np.isfinite(ecg).all():
            result["qc_reason"] = "nonfinite_ecg"
            return result

        cleaned, rpeaks = detect_rpeaks(ecg, fs, detector_method)
        quality_class = nk.ecg_quality(
            cleaned,
            rpeaks=rpeaks,
            sampling_rate=fs,
            method="zhao2018",
            approach="fuzzy",
        )
        artifacts, corrected_rpeaks = correct_rpeaks(rpeaks, fs)
        rr_clean = np.diff(corrected_rpeaks) / fs * 1_000.0
        rr_clean = rr_clean[rr_clean > 0]
        artifact_indices = {
            int(index)
            for category in ("ectopic", "missed", "extra", "longshort")
            for index in artifacts.get(category, [])
        }
        result["n_rpeaks"] = int(rpeaks.size)
        result["n_rr_clean"] = int(rr_clean.size)
        result["ecg_quality_class"] = str(quality_class)
        result["n_artifact_flags"] = len(artifact_indices)
        result["artifact_flag_fraction"] = (
            len(artifact_indices) / max(int(rpeaks.size), 1)
        )

        if rr_clean.size < 150:
            result["qc_reason"] = "too_few_clean_intervals"
            return result

        mean_rr = float(np.mean(rr_clean))
        sdnn = float(np.std(rr_clean, ddof=1))
        rmssd = float(np.sqrt(np.mean(np.diff(rr_clean) ** 2)))
        result.update({"MEAN_RR": mean_rr, "SDNN": sdnn, "RMSSD": rmssd})

        if not 400 <= mean_rr <= 1_500:
            result["qc_reason"] = "implausible_mean_rr"
            return result
        if not 5 <= sdnn <= 300:
            result["qc_reason"] = "implausible_sdnn"
            return result
        if not 5 <= rmssd <= 300:
            result["qc_reason"] = "implausible_rmssd"
            return result

        result.update({"qc_status": "pass", "qc_reason": ""})
        return result
    except Exception as error:  # retain the reason while allowing the cohort to finish
        result["qc_reason"] = normalize_processing_error(error)
        result["processing_error_type"] = type(error).__name__
        return result


def eligible_records(features: pd.DataFrame) -> pd.DataFrame:
    """Return records eligible for age-by-sex aggregation."""
    eligible = features.loc[
        (features["qc_status"] == "pass")
        & features["sex"].isin([0, 1])
        & features["age_group"].between(1, 15)
    ].copy()
    eligible["sex"] = eligible["sex"].astype(int)
    eligible["age_group"] = eligible["age_group"].astype(int)
    return eligible


def aggregate(features: pd.DataFrame) -> pd.DataFrame:
    rows: list[dict] = []
    passed = eligible_records(features)

    for age_label, source_groups in AGE_BANDS:
        age_rows = passed.loc[passed["age_group"].isin(source_groups)]
        for sex_code, sex_label in SEX_LABELS.items():
            sex_rows = age_rows.loc[age_rows["sex"] == sex_code]
            for metric in METRICS:
                values = sex_rows[metric].dropna().to_numpy(dtype=float)
                publish = values.size >= MIN_PUBLIC_CELL_N
                rows.append({
                    "age_band": age_label,
                    "source_age_groups": ",".join(map(str, source_groups)),
                    "sex": sex_label,
                    "metric": metric,
                    "n": int(values.size),
                    "suppressed": not publish,
                    "median_ms": float(np.median(values)) if publish else np.nan,
                    "q25_ms": float(np.quantile(values, 0.25)) if publish else np.nan,
                    "q75_ms": float(np.quantile(values, 0.75)) if publish else np.nan,
                    "mean_ms": float(np.mean(values)) if publish else np.nan,
                    "sd_ms": float(np.std(values, ddof=1)) if publish else np.nan,
                })
    return pd.DataFrame(rows)


def aggregate_sensitivity_summary(
    primary_features: pd.DataFrame,
    sensitivity_features: pd.DataFrame,
) -> dict:
    """Compare every publishable cell without exposing subject-level values."""
    primary_eligible = eligible_records(primary_features)
    sensitivity_eligible = eligible_records(sensitivity_features)
    primary_table = aggregate(primary_eligible)
    sensitivity_table = aggregate(sensitivity_eligible)
    merged = primary_table.merge(
        sensitivity_table,
        on=["age_band", "source_age_groups", "sex", "metric"],
        suffixes=("_primary", "_sensitivity"),
    )
    comparable = merged.loc[
        ~merged["suppressed_primary"] & ~merged["suppressed_sensitivity"]
    ].copy()
    comparable["median_shift_ms"] = (
        comparable["median_ms_sensitivity"] - comparable["median_ms_primary"]
    )
    max_row = None
    if not comparable.empty:
        max_index = comparable["median_shift_ms"].abs().idxmax()
        max_row = comparable.loc[max_index]
    cell_shifts = [
        {
            "age_band": str(row["age_band"]),
            "sex": str(row["sex"]),
            "metric": str(row["metric"]),
            "n_primary": int(row["n_primary"]),
            "n_sensitivity": int(row["n_sensitivity"]),
            "primary_median_ms": round(float(row["median_ms_primary"]), 3),
            "sensitivity_median_ms": round(
                float(row["median_ms_sensitivity"]), 3
            ),
            "median_shift_ms": round(float(row["median_shift_ms"]), 3),
        }
        for _, row in comparable.iterrows()
    ]
    return {
        "primary_eligible_records": int(len(primary_eligible)),
        "sensitivity_eligible_records": int(len(sensitivity_eligible)),
        "eligible_record_change": int(
            len(sensitivity_eligible) - len(primary_eligible)
        ),
        "maximum_absolute_age_sex_median_shift_ms": (
            round(abs(float(max_row["median_shift_ms"])), 3)
            if max_row is not None
            else None
        ),
        "maximum_shift_cell": (
            {
                "age_band": str(max_row["age_band"]),
                "sex": str(max_row["sex"]),
                "metric": str(max_row["metric"]),
                "median_shift_ms": round(float(max_row["median_shift_ms"]), 3),
            }
            if max_row is not None
            else None
        ),
        "published_cell_median_shifts": cell_shifts,
    }


def artifact_sensitivity(features: pd.DataFrame) -> dict:
    """Report how much a stricter artifact screen changes published medians."""
    eligible = eligible_records(features)
    strict = eligible.loc[eligible["artifact_flag_fraction"] <= 0.05]
    return {
        "screen": "exclude records with >5% Kubios artifact flags",
        "additional_exclusions": int(len(eligible) - len(strict)),
        **aggregate_sensitivity_summary(eligible, strict),
    }


def quality_class_sensitivity(features: pd.DataFrame) -> dict:
    """Quantify the effect of retaining only Zhao2018 Excellent signals."""
    eligible = eligible_records(features)
    strict = eligible.loc[eligible["ecg_quality_class"] == "Excellent"]
    return {
        "screen": (
            "retain only NeuroKit Zhao2018 fuzzy quality class Excellent; "
            "exclude Barely acceptable"
        ),
        "additional_exclusions": int(len(eligible) - len(strict)),
        **aggregate_sensitivity_summary(eligible, strict),
    }


def channel_preference_key(row: pd.Series) -> tuple:
    """Rank a channel without using its SDNN or RMSSD value directly."""
    artifact_fraction = row.get("artifact_flag_fraction")
    artifact_score = (
        -float(artifact_fraction) if pd.notna(artifact_fraction) else -np.inf
    )
    return (
        row.get("qc_status") == "pass",
        QUALITY_CLASS_RANK.get(str(row.get("ecg_quality_class", "")), -1),
        artifact_score,
    )


def replace_with_alternate_channels(
    primary: pd.DataFrame,
    alternate: pd.DataFrame,
    *,
    quality_preferred: bool,
) -> tuple[pd.DataFrame, int]:
    """Replace multi-lead primary rows, optionally only when quality ranks higher."""
    selected = primary.set_index("subject_id", drop=False).copy()
    alternate_selected = 0
    for _, alternate_row in alternate.iterrows():
        subject_id = str(alternate_row["subject_id"])
        primary_row = selected.loc[subject_id]
        if quality_preferred and not (
            channel_preference_key(alternate_row) > channel_preference_key(primary_row)
        ):
            continue
        selected.loc[subject_id, selected.columns] = alternate_row[selected.columns]
        alternate_selected += 1
    return selected.reset_index(drop=True), alternate_selected


def channel_sensitivity(
    primary: pd.DataFrame,
    alternate: pd.DataFrame,
) -> dict:
    """Compare fixed lead II, lead I, and automated quality preference."""
    primary_multi = primary.loc[
        primary["subject_id"].isin(alternate["subject_id"])
    ].set_index("subject_id")
    alternate_multi = alternate.set_index("subject_id")
    paired = primary_multi.join(
        alternate_multi,
        how="inner",
        lsuffix="_primary",
        rsuffix="_alternate",
    )
    status_cross = (
        paired.groupby(["qc_status_primary", "qc_status_alternate"])
        .size()
        .to_dict()
    )
    both_pass = paired.loc[
        (paired["qc_status_primary"] == "pass")
        & (paired["qc_status_alternate"] == "pass")
    ]
    paired_metric_differences = {}
    for metric in ("MEAN_RR", *METRICS):
        absolute_difference = (
            both_pass[f"{metric}_alternate"] - both_pass[f"{metric}_primary"]
        ).abs()
        paired_metric_differences[metric] = {
            "median_absolute_difference_ms": round(
                float(absolute_difference.median()), 3
            ),
            "p95_absolute_difference_ms": round(
                float(absolute_difference.quantile(0.95)), 3
            ),
            "maximum_absolute_difference_ms": round(
                float(absolute_difference.max()), 3
            ),
        }

    quality_preferred, alternate_selected = replace_with_alternate_channels(
        primary,
        alternate,
        quality_preferred=True,
    )
    lead_i, lead_i_selected = replace_with_alternate_channels(
        primary,
        alternate,
        quality_preferred=False,
    )
    quality_summary = aggregate_sensitivity_summary(primary, quality_preferred)
    lead_i_summary = aggregate_sensitivity_summary(primary, lead_i)
    return {
        "multi_ecg_records": int(len(alternate)),
        "source_channel_mapping": {
            "ECG1": "lead I for Task Force Monitor records (device 0)",
            "ECG2": "lead II for Task Force Monitor records (device 0)",
            "ECG": "lead II for MP150/CNAP records (device 1)",
        },
        "primary_policy": (
            "fixed lead II: ECG2 for device 0 and ECG for device 1"
        ),
        "primary_qc_pass_multi": int(
            (primary_multi["qc_status"] == "pass").sum()
        ),
        "alternate_qc_pass_multi": int(
            (alternate_multi["qc_status"] == "pass").sum()
        ),
        "qc_status_cross_table": {
            f"primary_{primary_status}__alternate_{alternate_status}": int(count)
            for (primary_status, alternate_status), count in status_cross.items()
        },
        "both_channels_qc_pass": int(len(both_pass)),
        "paired_metric_differences_among_both_pass": paired_metric_differences,
        "quality_preferred_sensitivity": {
            "selection_rule": (
                "prefer QC pass, then Zhao2018 quality class, then lower Kubios "
                "artifact fraction; ties retain the primary channel; SDNN and "
                "RMSSD are not direct ranking inputs"
            ),
            "primary_channels_selected_among_multi": int(
                len(alternate) - alternate_selected
            ),
            "alternate_channels_selected_among_multi": int(alternate_selected),
            "qc_pass_records": int(
                (quality_preferred["qc_status"] == "pass").sum()
            ),
            **quality_summary,
        },
        "lead_i_sensitivity": {
            "selection_rule": (
                "use ECG1/lead I for all device-0 records; device-1 records "
                "remain ECG/lead II because only one ECG signal is available"
            ),
            "alternate_channels_selected_among_multi": int(lead_i_selected),
            "qc_pass_records": int((lead_i["qc_status"] == "pass").sum()),
            **lead_i_summary,
        },
        "interpretation": (
            "Aggregate medians are robust to channel policy, but individual-channel "
            "differences can be large. Automated channel preference is a sensitivity "
            "analysis, not independently validated ground truth."
        ),
    }


def relative_count_difference(
    first: pd.Series,
    second: pd.Series,
) -> pd.Series:
    """Return absolute beat-count difference divided by the larger count."""
    denominator = pd.concat([first, second], axis=1).max(axis=1)
    return (first - second).abs() / denominator.replace(0, np.nan)


def apply_count_agreement_qc(
    primary: pd.DataFrame,
    lead_i: pd.DataFrame,
    elgendi: pd.DataFrame,
) -> pd.DataFrame:
    """Apply pre-declared detector and, where available, lead count agreement QC."""
    checked = primary.set_index("subject_id", drop=False).copy()
    elgendi_by_id = elgendi.set_index("subject_id")
    lead_i_by_id = lead_i.set_index("subject_id")

    checked["elgendi_n_rpeaks"] = elgendi_by_id["n_rpeaks"].reindex(checked.index)
    checked["detector_count_relative_difference"] = relative_count_difference(
        checked["n_rpeaks"],
        checked["elgendi_n_rpeaks"],
    )
    checked["lead_i_n_rpeaks"] = lead_i_by_id["n_rpeaks"].reindex(checked.index)
    checked["lead_count_relative_difference"] = relative_count_difference(
        checked["n_rpeaks"],
        checked["lead_i_n_rpeaks"],
    )

    detector_disagreement = (
        checked["detector_count_relative_difference"].isna()
        | (
            checked["detector_count_relative_difference"]
            > COUNT_DISAGREEMENT_TOLERANCE
        )
    )
    has_alternate_lead = checked["ecg_channel_count"] > 1
    lead_disagreement = has_alternate_lead & (
        checked["lead_count_relative_difference"].isna()
        | (
            checked["lead_count_relative_difference"]
            > COUNT_DISAGREEMENT_TOLERANCE
        )
    )
    was_pass = checked["qc_status"] == "pass"
    detector_exclusion = was_pass & detector_disagreement
    lead_exclusion = was_pass & ~detector_disagreement & lead_disagreement
    checked.loc[detector_exclusion, "qc_status"] = "error"
    checked.loc[detector_exclusion, "qc_reason"] = "detector_count_disagreement"
    checked.loc[lead_exclusion, "qc_status"] = "error"
    checked.loc[lead_exclusion, "qc_reason"] = "lead_count_disagreement"
    return checked.reset_index(drop=True)


def count_comparison_summary(
    reference: pd.DataFrame,
    comparison: pd.DataFrame,
) -> dict:
    """Summarize count disagreement without publishing subject identifiers."""
    reference_by_id = reference.set_index("subject_id")
    comparison_by_id = comparison.set_index("subject_id")
    paired = reference_by_id[["n_rpeaks"]].join(
        comparison_by_id[["n_rpeaks"]],
        how="inner",
        lsuffix="_reference",
        rsuffix="_comparison",
    )
    difference = relative_count_difference(
        paired["n_rpeaks_reference"],
        paired["n_rpeaks_comparison"],
    )
    finite = difference.dropna()
    return {
        "paired_records": int(len(paired)),
        "records_with_computable_counts": int(len(finite)),
        "relative_difference_quantiles": {
            "median": round(float(finite.quantile(0.50)), 4),
            "p90": round(float(finite.quantile(0.90)), 4),
            "p95": round(float(finite.quantile(0.95)), 4),
            "p99": round(float(finite.quantile(0.99)), 4),
            "maximum": round(float(finite.max()), 4),
        },
        "missing_or_over_10_percent": int(
            (difference.isna() | (difference > COUNT_DISAGREEMENT_TOLERANCE)).sum()
        ),
    }


def detector_and_lead_audit(
    primary_raw: pd.DataFrame,
    primary_checked: pd.DataFrame,
    lead_i: pd.DataFrame,
    elgendi: pd.DataFrame,
    pantompkins: pd.DataFrame,
) -> dict:
    """Document detector/lead choices and their aggregate sensitivity."""
    raw_by_id = primary_raw.set_index("subject_id")
    checked_by_id = primary_checked.set_index("subject_id")
    newly_excluded = (
        (raw_by_id["qc_status"] == "pass")
        & (checked_by_id["qc_status"] != "pass")
    )
    return {
        "primary_channel": (
            "lead II fixed before signal inspection: ECG2 for Task Force Monitor "
            "records and ECG for MP150/CNAP records"
        ),
        "primary_detector": (
            "NeuroKit morphology-aware detector on NeuroKit-cleaned ECG"
        ),
        "count_disagreement_threshold": (
            "exclude a primary-QC-pass record when beat counts differ by more than "
            "10% of the larger count, or when a comparison count is unavailable"
        ),
        "secondary_detector_used_for_qc": {
            "method": "Elgendi 2010 via NeuroKit",
            **count_comparison_summary(primary_raw, elgendi),
            "aggregate_metric_sensitivity": aggregate_sensitivity_summary(
                primary_raw,
                elgendi,
            ),
        },
        "two_lead_count_audit": {
            "method": (
                "same NeuroKit detector on ECG2/lead II versus ECG1/lead I for "
                "the Task Force Monitor subset"
            ),
            **count_comparison_summary(primary_raw, lead_i),
        },
        "pantompkins_sensitivity_not_used_for_qc": {
            "method": "Pan-Tompkins 1985 via NeuroKit on fixed lead II",
            **count_comparison_summary(primary_raw, pantompkins),
            "aggregate_metric_sensitivity": aggregate_sensitivity_summary(
                primary_raw,
                pantompkins,
            ),
            "reason_not_adopted": (
                "This implementation frequently double-detected T waves and its "
                "non-apex timestamps inflated short-term HRV. It is reported as a "
                "sensitivity comparator rather than selected after seeing HRV values."
            ),
        },
        "concordance_qc_effect": {
            "new_exclusions_among_primary_qc_pass": int(newly_excluded.sum()),
            **aggregate_sensitivity_summary(primary_raw, primary_checked),
        },
    }


def render_chart(table: pd.DataFrame, public_dir: Path) -> None:
    matplotlib_cache = Path(
        os.environ.get("MPLCONFIGDIR", "/private/tmp/feelmo-matplotlib-cache")
    )
    matplotlib_cache.mkdir(parents=True, exist_ok=True)
    os.environ["MPLCONFIGDIR"] = str(matplotlib_cache)

    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from matplotlib import font_manager

    japanese_fonts = list(Path("/System/Library/Fonts").glob("*W8.ttc"))
    font_family = "sans-serif"
    if japanese_fonts:
        font_manager.fontManager.addfont(str(japanese_fonts[0]))
        font_family = font_manager.FontProperties(fname=str(japanese_fonts[0])).get_name()

    plt.rcParams.update({
        "font.family": font_family,
        "font.size": 10,
        "svg.hashsalt": "feelmo-autonomic-aging-v1",
    })
    colors = {"female": "#B55C7A", "male": "#3E78A8"}
    labels = {"female": "女性", "male": "男性"}
    x = np.arange(len(AGE_BANDS))
    fig, axes = plt.subplots(1, 2, figsize=(10.4, 4.1), sharex=True)

    for axis, metric in zip(axes, METRICS):
        for sex in ("female", "male"):
            rows = table.loc[
                (table["metric"] == metric) & (table["sex"] == sex)
            ].set_index("age_band").reindex([label for label, _ in AGE_BANDS])
            median = rows["median_ms"].to_numpy(dtype=float)
            q25 = rows["q25_ms"].to_numpy(dtype=float)
            q75 = rows["q75_ms"].to_numpy(dtype=float)
            axis.plot(x, median, marker="o", linewidth=2, color=colors[sex], label=labels[sex])
            axis.fill_between(x, q25, q75, color=colors[sex], alpha=0.13)

        axis.set_title(f"{metric}（中央値と四分位範囲）")
        axis.set_ylabel("ミリ秒（ms）")
        axis.set_xticks(x, [label for label, _ in AGE_BANDS], rotation=25)
        axis.grid(axis="y", color="#DCE3EA", linewidth=0.8)
        axis.spines[["top", "right"]].set_visible(False)

    axes[0].legend(frameon=False)
    fig.suptitle("安静時5分ECGにおける年齢別・性別HRV分布", fontweight="bold")
    fig.text(
        0.5,
        0.01,
        "出典: Schumann & Bär, Autonomic Aging v1.0.0（Feelmoによる再解析）",
        ha="center",
        fontsize=8,
        color="#5B6472",
    )
    fig.tight_layout(rect=(0, 0.05, 1, 0.94))
    fig.savefig(
        public_dir / "autonomic-aging-hrv-by-age-sex.svg",
        bbox_inches="tight",
        metadata={"Date": None},
    )
    fig.savefig(
        public_dir / "autonomic-aging-hrv-by-age-sex.png",
        dpi=180,
        bbox_inches="tight",
        metadata={"Software": "Feelmo Autonomic Aging aggregate builder"},
    )
    plt.close(fig)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument(
        "--data-dir",
        type=Path,
        default=Path("/Volumes/SSD/datasets/public-hrv/autonomic-aging"),
    )
    parser.add_argument(
        "--public-dir",
        type=Path,
        default=DEFAULT_PUBLIC_DIR,
    )
    parser.add_argument(
        "--subject-output",
        type=Path,
        default=DEFAULT_SUBJECT_OUTPUT,
    )
    parser.add_argument("--jobs", type=int, default=8)
    parser.add_argument("--segment-seconds", type=int, default=300)
    parser.add_argument("--max-subjects", type=int)
    parser.add_argument(
        "--verify-source-checksums",
        action="store_true",
        help="Verify every file listed in the source release SHA256SUMS.txt.",
    )
    return parser.parse_args()


def lead_ii_channel_offset(row: dict) -> int:
    """Map the source Device field to the fixed lead-II channel position."""
    try:
        return 1 if int(row.get("Device")) == 0 else 0
    except (TypeError, ValueError):
        return 0


def main() -> int:
    args = parse_args()
    paths = Paths(args.data_dir, args.public_dir, args.subject_output)
    if (args.segment_seconds != 300 or args.max_subjects is not None) and (
        paths.public_dir.resolve() == DEFAULT_PUBLIC_DIR.resolve()
        or paths.subject_output.resolve() == DEFAULT_SUBJECT_OUTPUT.resolve()
    ):
        raise SystemExit(
            "Refusing to overwrite the production five-minute outputs with a "
            "nonstandard or partial run. Pass explicit --public-dir and "
            "--subject-output paths for experiments."
        )

    subject_info_path = paths.data_dir / "subject-info.csv"
    checksum_verification = (
        verify_source_checksums(paths.data_dir)
        if args.verify_source_checksums
        else {
            "performed": False,
            "listed_files": 0,
            "failed_files": None,
            "checksum_manifest_sha256": sha256(paths.data_dir / "SHA256SUMS.txt"),
        }
    )
    info = pd.read_csv(subject_info_path, dtype={"ID": str})
    if args.max_subjects:
        info = info.head(args.max_subjects)

    info_rows = info.to_dict(orient="records")

    # Threads avoid repeatedly importing the scientific stack in macOS spawn
    # workers; scipy/neurokit's numerical kernels release the GIL.
    primary_rows = Parallel(n_jobs=args.jobs, backend="threading", verbose=10)(
        delayed(extract_subject)(
            row,
            paths.data_dir,
            args.segment_seconds,
            lead_ii_channel_offset(row),
            "neurokit",
        )
        for row in info_rows
    )
    primary_raw = pd.DataFrame(primary_rows)
    multi_subjects = set(
        primary_raw.loc[
            primary_raw["ecg_channel_count"] > 1,
            "subject_id",
        ]
    )
    multi_rows = [
        row for row in info_rows if str(row["ID"]).zfill(4) in multi_subjects
    ]
    lead_i_rows = Parallel(n_jobs=args.jobs, backend="threading", verbose=10)(
        delayed(extract_subject)(
            row,
            paths.data_dir,
            args.segment_seconds,
            0,
            "neurokit",
        )
        for row in multi_rows
    )
    elgendi_rows = Parallel(n_jobs=args.jobs, backend="threading", verbose=10)(
        delayed(extract_subject)(
            row,
            paths.data_dir,
            args.segment_seconds,
            lead_ii_channel_offset(row),
            "elgendi2010",
        )
        for row in info_rows
    )
    pantompkins_rows = Parallel(n_jobs=args.jobs, backend="threading", verbose=10)(
        delayed(extract_subject)(
            row,
            paths.data_dir,
            args.segment_seconds,
            lead_ii_channel_offset(row),
            "pantompkins1985",
        )
        for row in info_rows
    )
    lead_i = pd.DataFrame(lead_i_rows)
    elgendi = pd.DataFrame(elgendi_rows)
    pantompkins = pd.DataFrame(pantompkins_rows)
    features = apply_count_agreement_qc(primary_raw, lead_i, elgendi)
    paths.subject_output.parent.mkdir(parents=True, exist_ok=True)
    with tempfile.NamedTemporaryFile(
        dir=paths.subject_output.parent,
        prefix=f".{paths.subject_output.stem}-",
        suffix=".parquet",
        delete=False,
    ) as subject_temp:
        subject_temp_path = Path(subject_temp.name)
    try:
        features.to_parquet(subject_temp_path, index=False)
        os.replace(subject_temp_path, paths.subject_output)
    finally:
        subject_temp_path.unlink(missing_ok=True)

    table = aggregate(features)
    paths.public_dir.parent.mkdir(parents=True, exist_ok=True)

    passed = features.loc[features["qc_status"] == "pass"]
    aggregate_eligible = passed.loc[
        passed["sex"].isin([0, 1]) & passed["age_group"].between(1, 15)
    ]
    manifest = {
        "dataset": "Autonomic Aging",
        "dataset_version": DATASET_VERSION,
        "dataset_doi": DATASET_DOI,
        "paper_doi": PAPER_DOI,
        "independent_qc_paper_doi": INDEPENDENT_QC_DOI,
        "source_license": "Open Database License 1.0 (ODbL)",
        "source_subject_info_sha256": sha256(subject_info_path),
        "source_release_checksum_verification": checksum_verification,
        "analysis_script_sha256": sha256(SCRIPT_PATH),
        "software_versions": {
            name: package_version(name)
            for name in (
                "numpy",
                "pandas",
                "scipy",
                "wfdb",
                "neurokit2",
                "joblib",
                "matplotlib",
                "pyarrow",
            )
        },
        "input_records": int(len(features)),
        "qc_pass_records": int(len(passed)),
        "qc_excluded_records": int(len(features) - len(passed)),
        "aggregate_eligible_records": int(len(aggregate_eligible)),
        "aggregate_missing_age_or_sex": int(len(passed) - len(aggregate_eligible)),
        "minimum_public_cell_n": MIN_PUBLIC_CELL_N,
        "qc_exclusion_reasons": features.loc[
            features["qc_status"] != "pass", "qc_reason"
        ].value_counts().to_dict(),
        "signal_quality_classes_among_qc_pass": (
            passed["ecg_quality_class"].value_counts().to_dict()
        ),
        "artifact_flag_sensitivity": artifact_sensitivity(features),
        "signal_quality_class_sensitivity": quality_class_sensitivity(features),
        "channel_sensitivity_before_count_agreement_qc": channel_sensitivity(
            primary_raw,
            lead_i,
        ),
        "detector_and_lead_audit": detector_and_lead_audit(
            primary_raw,
            features,
            lead_i,
            elgendi,
            pantompkins,
        ),
        "sex_metadata": {
            "source_field": "Sex",
            "encoding": {"0": "male", "1": "female"},
            "semantics": (
                "Binary male/female sex metadata supplied by the source dataset; "
                "it is not a measure of gender identity."
            ),
        },
        "analysis_window": f"first {args.segment_seconds} seconds of resting ECG",
        "r_peak_pipeline": (
            "Fixed lead II + NeuroKit ECG cleaning and morphology-aware peak "
            "detection + Lipponen-Tarvainen/Kubios correction at native sampling "
            "rate; exclude >10% beat-count disagreement with Elgendi2010 and, "
            "where available, the simultaneously recorded lead I"
        ),
        "methodological_distinction": (
            "This is an independent automated reanalysis of the first five minutes, "
            "not a reproduction of the source paper or Calderon-Juarez et al. 2023, "
            "which used a visually supervised final five-minute segment."
        ),
        "published_outputs": [
            "autonomic-aging-age-sex-hrv.csv",
            "autonomic-aging-hrv-by-age-sex.svg",
            "autonomic-aging-hrv-by-age-sex.png",
            PUBLIC_SCRIPT_NAME,
        ],
        "privacy": "Only aggregate statistics are public; subject-level features are not published.",
        "measurement_warning": (
            "These five-minute clinical ECG distributions are not interchangeable with "
            "intermittent Apple Watch HealthKit SDNN samples."
        ),
    }

    generated_names = [
        "autonomic-aging-age-sex-hrv.csv",
        "autonomic-aging-hrv-by-age-sex.svg",
        "autonomic-aging-hrv-by-age-sex.png",
        PUBLIC_SCRIPT_NAME,
    ]
    with tempfile.TemporaryDirectory(
        dir=paths.public_dir.parent,
        prefix=".autonomic-aging-build-",
    ) as temp_name:
        temp_dir = Path(temp_name)
        table_path = temp_dir / generated_names[0]
        table.to_csv(table_path, index=False, float_format="%.3f")
        render_chart(table, temp_dir)
        shutil.copy2(SCRIPT_PATH, temp_dir / PUBLIC_SCRIPT_NAME)
        manifest["published_output_sha256"] = {
            name: sha256(temp_dir / name) for name in generated_names
        }
        (temp_dir / "analysis-manifest.json").write_text(
            json.dumps(manifest, ensure_ascii=False, indent=2) + "\n",
            encoding="utf-8",
        )
        paths.public_dir.mkdir(parents=True, exist_ok=True)
        for name in [*generated_names, "analysis-manifest.json"]:
            os.replace(temp_dir / name, paths.public_dir / name)

    print(json.dumps(manifest, ensure_ascii=False, indent=2))
    print(f"Aggregate table: {paths.public_dir / generated_names[0]}")
    print(f"Subject-level artifact (not public): {paths.subject_output}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())