#!/usr/bin/env python3
"""Build and live-verify an inventory of public complete WES BAMs.

Only official 1000 Genomes Phase 3 and NIST/NCBI GIAB indexes are used.
Each BAM and BAI is probed with an HTTP byte-range GET; no complete BAM is
downloaded.  A file is counted as downloadable only when both probes return
HTTP 206, expose a total byte size, and begin with the expected BAM/BGZF and
BAI magic bytes.
"""

from __future__ import annotations

import argparse
import csv
import gzip
import hashlib
import json
import re
import ssl
import statistics
import sys
import threading
import time
import urllib.error
import urllib.request
from concurrent.futures import ThreadPoolExecutor, as_completed
from datetime import datetime, timezone
from pathlib import Path
from typing import Any


ONE_KGP_INDEX_URL = (
    "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/"
    "20130502.phase3.exome.alignment.index"
)
ONE_KGP_EXOME_ONLY_INDEX_URL = (
    "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/"
    "phase3_EX_or_LC_only_alignment/20130502.exome_only.alignment.index"
)
ONE_KGP_RELEASE_URL = (
    "https://www.internationalgenome.org/announcements/"
    "official-release-phase3-alignment-data-available-2013-05-25/"
)
ONE_KGP_HSMETRICS_URL = (
    "https://1000genomes.s3.amazonaws.com/alignment_indices/"
    "20130502.exome.alignment.index.HsMetrics.gz"
)
ONE_KGP_S3_BASE = "https://1000genomes.s3.amazonaws.com/phase3/"
ONE_KGP_S3_ROOT = "https://1000genomes.s3.amazonaws.com"
ONE_KGP_EBI_PHASE3_BASE = "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/"
ONE_KGP_EBI_ROOT = "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp"
GIAB_REPO_URL = "https://github.com/genome-in-a-bottle/giab_data_indexes"
GIAB_FTP_ROOT = "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/"
GIAB_NA12878_ION_DIRECTORY_URL = (
    "https://ftp.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/ion_exome/"
)

GIAB_MASTER_INDEX_NAMES = {
    "alignment.index.AJtrio_IonTorrent_exome_TMAP_GRCh37_07292015",
    "alignment.index.AJtrio_OsloUniversityHospital_IlluminaExome_bwamem_GRCh37_11252015",
    "alignment.index.ChineseTrio_HG005_IonTorrent_exome_TMAP_GRCh37_09232015",
    "alignment.index.Chinesetrio_HG005_OsloUniversityHospital_IlluminaExome_bwamem_GRCh37_11252015",
    "alignment.index.NA12878_HiSeq_Exome_Garvan_GRCh37_09252015",
    "alignment.index.NA12878_IonTorrent_exome_TMAP_GRCh37_09252015",
    "alignment.index.NA12878_TruSeq_Exome_Nebraska_GRCh37_09252015",
}

CSV_FIELDS = [
    "repository",
    "collection",
    "dataset_id",
    "official_index_membership",
    "logical_alignment_key",
    "logical_duplicate_or_reheader_note",
    "sample_id",
    "sample_alias",
    "individual_key",
    "technology",
    "reference_build",
    "mean_target_coverage",
    "pct_target_bases_20x",
    "pct_target_bases_30x",
    "official_depth_metrics_available",
    "official_depth_metrics_url",
    "alignment_format",
    "file_name",
    "alignment_url",
    "alternate_alignment_mirror_url",
    "alignment_bytes",
    "alignment_md5_from_official_index",
    "alignment_http_status",
    "alignment_content_range",
    "alignment_accept_ranges",
    "alignment_first_bytes_hex",
    "alignment_magic_ok",
    "index_url",
    "alternate_index_mirror_url",
    "index_bytes",
    "index_md5_from_official_index",
    "index_http_status",
    "index_content_range",
    "index_accept_ranges",
    "index_first_bytes_hex",
    "index_magic_ok",
    "downloadable",
    "verification_attempts_alignment",
    "verification_attempts_index",
    "verification_error_alignment",
    "verification_error_index",
    "official_source_index_url",
    "official_source_index_file",
    "deduplication_and_mirror_note",
    "verified_at_utc",
]


def now_utc() -> str:
    return datetime.now(timezone.utc).replace(microsecond=0).isoformat()


def ftp_to_https(url: str) -> str:
    return re.sub(
        r"^ftp://(?:ftp-trace|ftp)\.ncbi\.(?:nlm\.)?nih\.gov/",
        "https://ftp.ncbi.nlm.nih.gov/",
        url,
    )


def giab_sample(url: str) -> tuple[str, str, str]:
    match = re.search(r"/(HG00[2-7])_[^/]+/", url)
    if match:
        sample = match.group(1)
        return sample, "", sample
    if "/NA12878/" in url:
        return "NA12878", "HG001", "HG001/NA12878"
    raise ValueError(f"Cannot determine GIAB sample from {url}")


def giab_technology(index_name: str) -> str:
    if "IonTorrent" in index_name:
        return "Ion Proton exome"
    if "OsloUniversityHospital" in index_name:
        return "Illumina exome (Oslo University Hospital)"
    if "HiSeq_Exome_Garvan" in index_name:
        return "Illumina HiSeq exome (Garvan)"
    if "TruSeq_Exome_Nebraska" in index_name:
        return "Illumina TruSeq exome (Nebraska)"
    return "Exome sequencing"


def parse_1000g(
    index_path: Path,
    *,
    exome_only: bool = False,
    hsmetrics: dict[str, dict[str, str]] | None = None,
) -> list[dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    with index_path.open("r", encoding="utf-8", newline="") as handle:
        reader = csv.DictReader(handle, delimiter="\t")
        for raw in reader:
            bam_path = raw["BAM FILE"]
            name = Path(bam_path).name
            if ".mapped." not in name or ".chrom" in name or ".unmapped." in name:
                continue
            if exome_only:
                sample_match = re.search(r"/data/([^/]+)/exome_alignment/", bam_path)
                if not sample_match:
                    raise RuntimeError(f"Unexpected exome-only BAM path: {bam_path}")
                sample = sample_match.group(1)
                alignment_url = ONE_KGP_S3_ROOT + bam_path
                index_url = ONE_KGP_S3_ROOT + raw["BAI FILE"]
                alternate_alignment_url = ONE_KGP_EBI_ROOT + bam_path
                alternate_index_url = ONE_KGP_EBI_ROOT + raw["BAI FILE"]
                collection = "Phase 3 QA-passed exome-only alignments"
                dataset_id = "1KGP_Phase3_20130502_exome_only"
                source_index_url = ONE_KGP_EXOME_ONLY_INDEX_URL
            else:
                sample = bam_path.split("/")[1]
                alignment_url = ONE_KGP_S3_BASE + bam_path
                index_url = ONE_KGP_S3_BASE + raw["BAI FILE"]
                alternate_alignment_url = ONE_KGP_EBI_PHASE3_BASE + bam_path
                alternate_index_url = ONE_KGP_EBI_PHASE3_BASE + raw["BAI FILE"]
                collection = "Phase 3 official exome alignments"
                dataset_id = "1KGP_Phase3_20130502_exome"
                source_index_url = ONE_KGP_INDEX_URL
            rows.append(
                {
                    "repository": "1000 Genomes Project",
                    "collection": collection,
                    "dataset_id": dataset_id,
                    "official_index_membership": True,
                    "logical_alignment_key": f"1KGP:{sample}:{raw['BAM MD5']}",
                    "logical_duplicate_or_reheader_note": "",
                    "sample_id": sample,
                    "sample_alias": "HG001" if sample == "NA12878" else "",
                    "individual_key": "HG001/NA12878" if sample == "NA12878" else sample,
                    "technology": "Illumina exome",
                    "reference_build": "GRCh37",
                    "mean_target_coverage": (
                        hsmetrics.get(name, {}).get("MEAN_TARGET_COVERAGE", "")
                        if hsmetrics
                        else ""
                    ),
                    "pct_target_bases_20x": (
                        hsmetrics.get(name, {}).get("PCT_TARGET_BASES_20X", "")
                        if hsmetrics
                        else ""
                    ),
                    "pct_target_bases_30x": (
                        hsmetrics.get(name, {}).get("PCT_TARGET_BASES_30X", "")
                        if hsmetrics
                        else ""
                    ),
                    "official_depth_metrics_available": bool(hsmetrics and name in hsmetrics),
                    "official_depth_metrics_url": ONE_KGP_HSMETRICS_URL if hsmetrics else "",
                    "alignment_format": "BAM",
                    "file_name": name,
                    "alignment_url": alignment_url,
                    "alternate_alignment_mirror_url": alternate_alignment_url,
                    "alignment_md5_from_official_index": raw["BAM MD5"],
                    "index_url": index_url,
                    "alternate_index_mirror_url": alternate_index_url,
                    "index_md5_from_official_index": raw["BAI MD5"],
                    "official_source_index_url": source_index_url,
                    "official_source_index_file": index_path.name,
                    "deduplication_and_mirror_note": (
                        "One alignment object counted once. The S3 URL was live-probed; "
                        "the EBI URL is an alternate official mirror and is not an "
                        "additional alignment."
                    ),
                }
            )
    return rows


def parse_giab(giab_root: Path) -> list[dict[str, Any]]:
    rows: list[dict[str, Any]] = []
    discovered = {
        path.name: path
        for path in giab_root.rglob("alignment.index.*")
        if path.name in GIAB_MASTER_INDEX_NAMES
    }
    missing = sorted(GIAB_MASTER_INDEX_NAMES - discovered.keys())
    if missing:
        raise RuntimeError(f"Missing expected GIAB indexes: {missing}")

    for index_name in sorted(discovered):
        index_path = discovered[index_name]
        source_url = (
            f"{GIAB_REPO_URL}/blob/master/"
            f"{index_path.parent.name}/{index_name}"
        )
        with index_path.open("r", encoding="utf-8") as handle:
            lines = [line.strip() for line in handle if line.strip()]
        for line in lines[1:]:
            parts = line.split()
            if len(parts) != 4:
                raise RuntimeError(f"Unexpected GIAB index row in {index_name}: {line}")
            bam_url, bam_md5, bai_url, bai_md5 = parts
            bam_url = ftp_to_https(bam_url)
            bai_url = ftp_to_https(bai_url)
            sample, alias, individual_key = giab_sample(bam_url)
            file_name = bam_url.rsplit("/", 1)[-1]
            is_na12878_ion = file_name == "IonXpress_020_rawlib.b37.bam"
            rows.append(
                {
                    "repository": "NIST Genome in a Bottle / NCBI",
                    "collection": "GIAB public exome alignments",
                    "dataset_id": index_name.removeprefix("alignment.index."),
                    "official_index_membership": True,
                    "logical_alignment_key": (
                        "GIAB:NA12878:IonXpress_020_rawlib"
                        if is_na12878_ion
                        else f"GIAB:{file_name}"
                    ),
                    "logical_duplicate_or_reheader_note": (
                        "The indexed b37 object and legacy hg19 directory object are "
                        "two distinct downloadable BAM objects from the same "
                        "NA12878 Ion exome run; count once for logical-run totals."
                        if is_na12878_ion
                        else ""
                    ),
                    "sample_id": sample,
                    "sample_alias": alias,
                    "individual_key": individual_key,
                    "technology": giab_technology(index_name),
                    "reference_build": "GRCh37",
                    "mean_target_coverage": "",
                    "pct_target_bases_20x": "",
                    "pct_target_bases_30x": "",
                    "official_depth_metrics_available": False,
                    "official_depth_metrics_url": "",
                    "alignment_format": "BAM",
                    "file_name": file_name,
                    "alignment_url": bam_url,
                    "alternate_alignment_mirror_url": "",
                    "alignment_md5_from_official_index": bam_md5,
                    "index_url": bai_url,
                    "alternate_index_mirror_url": "",
                    "index_md5_from_official_index": bai_md5,
                    "official_source_index_url": source_url,
                    "official_source_index_file": index_name,
                    "deduplication_and_mirror_note": (
                        "Distinct complete exome run retained as one alignment file. "
                        "Multiple GIAB runs may belong to the same individual; person "
                        "counts are deduplicated using individual_key. The official "
                        "FTP URL was normalized to HTTPS, not counted as a mirror."
                    ),
                }
            )
    legacy_name = "IonXpress_020_rawlib.hg19.bam"
    legacy_url = GIAB_NA12878_ION_DIRECTORY_URL + legacy_name
    rows.append(
        {
            "repository": "NIST Genome in a Bottle / NCBI",
            "collection": "GIAB public exome alignments",
            "dataset_id": "NA12878_IonTorrent_exome_TMAP_hg19_direct_directory",
            "official_index_membership": False,
            "logical_alignment_key": "GIAB:NA12878:IonXpress_020_rawlib",
            "logical_duplicate_or_reheader_note": (
                "Legacy hg19 object paired with the indexed b37 object. They are two "
                "distinct downloadable BAM objects from the same NA12878 Ion exome "
                "run (sizes differ by 39 bytes); count once for logical-run totals."
            ),
            "sample_id": "NA12878",
            "sample_alias": "HG001",
            "individual_key": "HG001/NA12878",
            "technology": "Ion Proton exome",
            "reference_build": "hg19 (GRCh37)",
            "mean_target_coverage": "",
            "pct_target_bases_20x": "",
            "pct_target_bases_30x": "",
            "official_depth_metrics_available": False,
            "official_depth_metrics_url": "",
            "alignment_format": "BAM",
            "file_name": legacy_name,
            "alignment_url": legacy_url,
            "alternate_alignment_mirror_url": "",
            "alignment_md5_from_official_index": "",
            "index_url": legacy_url + ".bai",
            "alternate_index_mirror_url": "",
            "index_md5_from_official_index": "",
            "official_source_index_url": GIAB_NA12878_ION_DIRECTORY_URL,
            "official_source_index_file": "NCBI directory listing (not GIAB index)",
            "deduplication_and_mirror_note": (
                "Counted as a distinct downloadable server object, but grouped with "
                "the b37 object under one logical_alignment_key and one individual."
            ),
        }
    )
    return rows


def total_from_headers(headers: Any, status: int) -> int | None:
    content_range = headers.get("Content-Range", "")
    match = re.search(r"/(\d+)$", content_range)
    if match:
        return int(match.group(1))
    length = headers.get("Content-Length")
    if status == 200 and length and length.isdigit():
        return int(length)
    return None


def range_probe(url: str, kind: str, retries: int = 4) -> dict[str, Any]:
    headers = {
        "Range": "bytes=0-15",
        "User-Agent": "LGASieve-public-WES-audit/2026-08-01",
        "Accept-Encoding": "identity",
    }
    last_error = ""
    for attempt in range(1, retries + 1):
        try:
            request = urllib.request.Request(url, headers=headers, method="GET")
            with urllib.request.urlopen(request, timeout=45, context=SSL_CONTEXT) as response:
                status = int(response.status)
                first = response.read(16)
                returned_headers = response.headers
            magic_ok = (
                first.startswith(b"\x1f\x8b\x08")
                if kind == "bam"
                else first.startswith(b"BAI\x01")
            )
            return {
                "status": status,
                "content_range": returned_headers.get("Content-Range", ""),
                "accept_ranges": returned_headers.get("Accept-Ranges", ""),
                "bytes": total_from_headers(returned_headers, status),
                "first_bytes_hex": first.hex(),
                "magic_ok": magic_ok,
                "attempts": attempt,
                "error": "",
            }
        except Exception as exc:  # network errors are captured in the audit
            last_error = f"{type(exc).__name__}: {exc}"
            if attempt < retries:
                time.sleep(min(2 ** (attempt - 1), 4))
    return {
        "status": 0,
        "content_range": "",
        "accept_ranges": "",
        "bytes": None,
        "first_bytes_hex": "",
        "magic_ok": False,
        "attempts": retries,
        "error": last_error,
    }


def verify_row(row: dict[str, Any]) -> dict[str, Any]:
    bam = range_probe(row["alignment_url"], "bam")
    bai = range_probe(row["index_url"], "bai")
    output = dict(row)
    output.update(
        {
            "alignment_bytes": bam["bytes"] or "",
            "alignment_http_status": bam["status"],
            "alignment_content_range": bam["content_range"],
            "alignment_accept_ranges": bam["accept_ranges"],
            "alignment_first_bytes_hex": bam["first_bytes_hex"],
            "alignment_magic_ok": bam["magic_ok"],
            "index_bytes": bai["bytes"] or "",
            "index_http_status": bai["status"],
            "index_content_range": bai["content_range"],
            "index_accept_ranges": bai["accept_ranges"],
            "index_first_bytes_hex": bai["first_bytes_hex"],
            "index_magic_ok": bai["magic_ok"],
            "downloadable": (
                bam["status"] == 206
                and bai["status"] == 206
                and bam["bytes"] is not None
                and bai["bytes"] is not None
                and bam["magic_ok"]
                and bai["magic_ok"]
            ),
            "verification_attempts_alignment": bam["attempts"],
            "verification_attempts_index": bai["attempts"],
            "verification_error_alignment": bam["error"],
            "verification_error_index": bai["error"],
            "verified_at_utc": now_utc(),
        }
    )
    return output


def sha256_file(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as handle:
        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def load_hsmetrics(path: Path) -> dict[str, dict[str, str]]:
    with gzip.open(path, "rt", encoding="utf-8", newline="") as handle:
        return {row["File_name"]: row for row in csv.DictReader(handle, delimiter="\t")}


def size_bins_decimal_gb(rows: list[dict[str, Any]]) -> dict[str, int]:
    bins = {
        "lt_1_GB": 0,
        "1_to_lt_5_GB": 0,
        "5_to_lt_10_GB": 0,
        "10_to_lt_20_GB": 0,
        "20_to_lt_30_GB": 0,
        "ge_30_GB": 0,
    }
    for row in rows:
        size = int(row["alignment_bytes"])
        if size < 1_000_000_000:
            bins["lt_1_GB"] += 1
        elif size < 5_000_000_000:
            bins["1_to_lt_5_GB"] += 1
        elif size < 10_000_000_000:
            bins["5_to_lt_10_GB"] += 1
        elif size < 20_000_000_000:
            bins["10_to_lt_20_GB"] += 1
        elif size < 30_000_000_000:
            bins["20_to_lt_30_GB"] += 1
        else:
            bins["ge_30_GB"] += 1
    return bins


def write_outputs(rows: list[dict[str, Any]], out_csv: Path, out_json: Path) -> None:
    rows.sort(key=lambda row: (row["repository"], row["sample_id"], row["file_name"]))
    with out_csv.open("w", encoding="utf-8", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=CSV_FIELDS, extrasaction="ignore")
        writer.writeheader()
        writer.writerows(rows)

    downloadable = [row for row in rows if row["downloadable"]]
    repository_summary: dict[str, Any] = {}
    for repository in sorted({row["repository"] for row in rows}):
        subset = [row for row in rows if row["repository"] == repository]
        good = [row for row in subset if row["downloadable"]]
        sizes = [int(row["alignment_bytes"]) for row in good]
        repository_summary[repository] = {
            "listed_alignment_files": len(subset),
            "verified_downloadable_alignment_files": len(good),
            "unique_individual_keys": len({row["individual_key"] for row in good}),
            "alignment_bytes_total": sum(sizes),
            "alignment_bytes_min": min(sizes) if sizes else None,
            "alignment_bytes_median": statistics.median(sizes) if sizes else None,
            "alignment_bytes_max": max(sizes) if sizes else None,
            "alignment_size_bins_decimal_GB": size_bins_decimal_gb(good),
        }

    dataset_summary: dict[str, Any] = {}
    for dataset_id in sorted({row["dataset_id"] for row in rows}):
        subset = [row for row in rows if row["dataset_id"] == dataset_id]
        good = [row for row in subset if row["downloadable"]]
        dataset_summary[dataset_id] = {
            "repository": subset[0]["repository"],
            "collection": subset[0]["collection"],
            "listed_alignment_files": len(subset),
            "verified_downloadable_alignment_files": len(good),
            "unique_individual_keys": len({row["individual_key"] for row in good}),
            "alignment_bytes_total": sum(int(row["alignment_bytes"]) for row in good),
            "alignment_size_bins_decimal_GB": size_bins_decimal_gb(good),
        }

    depth_rows = [
        row
        for row in downloadable
        if row.get("official_depth_metrics_available")
        and row.get("mean_target_coverage") not in ("", None)
    ]
    target_depth = [float(row["mean_target_coverage"]) for row in depth_rows]
    target_20x = [float(row["pct_target_bases_20x"]) for row in depth_rows]
    target_30x = [float(row["pct_target_bases_30x"]) for row in depth_rows]
    depth_summary = {
        "files_with_official_hsmetrics": len(depth_rows),
        "mean_target_coverage_median": statistics.median(target_depth),
        "mean_target_coverage_mean": statistics.mean(target_depth),
        "mean_target_coverage_min": min(target_depth),
        "mean_target_coverage_max": max(target_depth),
        "pct_target_bases_20x_median": statistics.median(target_20x),
        "pct_target_bases_30x_median": statistics.median(target_30x),
        "scope_note": (
            "Official HsMetrics matched all 2,535 main Phase 3 exomes. The "
            "16 exome-only technical samples and GIAB runs do not have entries "
            "in this 1000 Genomes HsMetrics table."
        ),
    }

    payload = {
        "schema_version": "1.0",
        "title": "Verified public complete human whole-exome alignment inventory",
        "generated_at_utc": now_utc(),
        "definition": (
            "A complete exome alignment is an official repository BAM representing "
            "an entire exome capture/run or merged exome sample, not a regional slice, "
            "targeted gene panel, raw FASTQ-only record, or VCF-only record."
        ),
        "counting_rule": (
            "Counted downloadable only when both BAM and BAI returned HTTP 206 to "
            "a bytes=0-15 request, reported total sizes, and had BAM/BGZF and BAI magic."
        ),
        "important_interpretation": {
            "alignment_file_count": len(downloadable),
            "verified_primary_alignment_url_count": len(downloadable),
            "verified_primary_index_url_count": len(downloadable),
            "verified_downloadable_object_count_bam_plus_index": 2 * len(downloadable),
            "logical_alignment_count_after_legacy_alias_deduplication": len(
                {row["logical_alignment_key"] for row in downloadable}
            ),
            "alternate_mirror_urls_not_counted_as_files": True,
            "unique_individual_count": len({row["individual_key"] for row in downloadable}),
            "technical_replicates": (
                "GIAB contains multiple WES runs for five benchmark individuals; file "
                "count is therefore larger than the number of people."
            ),
            "no_full_files_downloaded": True,
            "probe_bytes_per_endpoint": 16,
        },
        "summary": {
            "listed_alignment_files": len(rows),
            "verified_downloadable_alignment_files": len(downloadable),
            "failed_or_unverified_alignment_files": len(rows) - len(downloadable),
            "unique_individual_keys": len({row["individual_key"] for row in downloadable}),
            "alignment_bytes_total": sum(int(row["alignment_bytes"]) for row in downloadable),
            "index_bytes_total": sum(int(row["index_bytes"]) for row in downloadable),
            "alignment_size_bins_decimal_GB": size_bins_decimal_gb(downloadable),
            "repositories": repository_summary,
            "datasets": dataset_summary,
            "logical_alignment_count_after_legacy_alias_deduplication": len(
                {row["logical_alignment_key"] for row in downloadable}
            ),
            "official_depth_metrics": depth_summary,
        },
        "official_sources": [
            {
                "name": "1000 Genomes Phase 3 exome alignment index",
                "url": ONE_KGP_INDEX_URL,
            },
            {
                "name": "1000 Genomes Phase 3 exome-only alignment index",
                "url": ONE_KGP_EXOME_ONLY_INDEX_URL,
            },
            {
                "name": "1000 Genomes official Phase 3 alignment release notice",
                "url": ONE_KGP_RELEASE_URL,
            },
            {
                "name": "1000 Genomes official exome HsMetrics",
                "url": ONE_KGP_HSMETRICS_URL,
            },
            {
                "name": "NIST GIAB data indexes",
                "url": GIAB_REPO_URL,
            },
            {
                "name": "NCBI GIAB repository",
                "url": GIAB_FTP_ROOT,
            },
            {
                "name": "NCBI GIAB NA12878 Ion-exome directory listing",
                "url": GIAB_NA12878_ION_DIRECTORY_URL,
            },
        ],
        "exclusions": [
            {
                "category": "1000 Genomes chromosome-split BAMs",
                "count_in_official_indexes": 5102,
                "reason": "Fragments of samples, not complete exome alignments.",
            },
            {
                "category": "1000 Genomes unmapped BAMs",
                "count_in_official_indexes": 2551,
                "reason": "Unmapped-read supplements, not complete mapped exomes.",
            },
            {
                "category": "Regional BAM slices and targeted gene panels",
                "count_in_inventory": 0,
                "reason": "Explicitly outside the complete-WES definition.",
            },
            {
                "category": "Raw FASTQ/SRA and VCF-only public records",
                "count_in_inventory": 0,
                "reason": "No directly downloadable indexed complete BAM/CRAM pair.",
            },
            {
                "category": "Controlled-access human WES repositories",
                "count_in_inventory": 0,
                "reason": "Cannot verify anonymous direct downloadability.",
            },
        ],
        "inventory_csv": out_csv.name,
        "inventory_csv_sha256": sha256_file(out_csv),
        "rows": rows,
    }
    out_json.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8")


def main() -> int:
    parser = argparse.ArgumentParser()
    parser.add_argument("--one-kgp-index", required=True, type=Path)
    parser.add_argument("--one-kgp-exome-only-index", required=True, type=Path)
    parser.add_argument("--one-kgp-hsmetrics", required=True, type=Path)
    parser.add_argument("--giab-root", required=True, type=Path)
    parser.add_argument("--out-csv", required=True, type=Path)
    parser.add_argument("--out-json", required=True, type=Path)
    parser.add_argument("--workers", type=int, default=32)
    parser.add_argument("--limit", type=int, default=0)
    parser.add_argument(
        "--reuse-verified-csv",
        type=Path,
        help="Reuse rows already marked downloadable and retry only failures.",
    )
    args = parser.parse_args()

    hsmetrics = load_hsmetrics(args.one_kgp_hsmetrics)
    rows = (
        parse_1000g(args.one_kgp_index, hsmetrics=hsmetrics)
        + parse_1000g(args.one_kgp_exome_only_index, exome_only=True)
        + parse_giab(args.giab_root)
    )
    if args.limit:
        rows = rows[: args.limit]

    reused: list[dict[str, Any]] = []
    if args.reuse_verified_csv and args.reuse_verified_csv.exists():
        with args.reuse_verified_csv.open("r", encoding="utf-8", newline="") as handle:
            prior_rows = {
                row["alignment_url"]: row for row in csv.DictReader(handle)
            }
        pending: list[dict[str, Any]] = []
        for row in rows:
            prior = prior_rows.get(row["alignment_url"])
            if prior and prior.get("downloadable", "").lower() == "true":
                for key in (
                    "alignment_magic_ok",
                    "index_magic_ok",
                    "downloadable",
                    "official_depth_metrics_available",
                ):
                    prior[key] = prior.get(key, "").lower() == "true"
                for key in (
                    "alignment_bytes",
                    "alignment_http_status",
                    "index_bytes",
                    "index_http_status",
                    "verification_attempts_alignment",
                    "verification_attempts_index",
                ):
                    prior[key] = int(prior[key])
                reused.append({**prior, **row})
            else:
                pending.append(row)
        rows = pending
    print(
        f"Parsed {len(rows) + len(reused)} complete WES BAM/index pairs; "
        f"reused={len(reused)}, probing={len(rows)}",
        flush=True,
    )

    verified: list[dict[str, Any]] = list(reused)
    progress_lock = threading.Lock()
    completed = 0
    with ThreadPoolExecutor(max_workers=args.workers) as executor:
        futures = {executor.submit(verify_row, row): row for row in rows}
        for future in as_completed(futures):
            result = future.result()
            verified.append(result)
            with progress_lock:
                completed += 1
                if completed % 100 == 0 or completed == len(rows):
                    good = sum(1 for row in verified if row["downloadable"])
                    print(
                        f"Probed {completed}/{len(rows)} pending; "
                        f"downloadable total={good}",
                        flush=True,
                    )

    args.out_csv.parent.mkdir(parents=True, exist_ok=True)
    write_outputs(verified, args.out_csv, args.out_json)
    good = sum(1 for row in verified if row["downloadable"])
    print(f"Finished: {good}/{len(verified)} downloadable", flush=True)
    print(args.out_csv)
    print(args.out_json)
    return 0 if good == len(verified) else 2


SSL_CONTEXT = ssl.create_default_context()


if __name__ == "__main__":
    sys.exit(main())
