SimurgArray/src/simurg/pipeline/cn.py

"""``simurg cn build-model`` / ``simurg cn call`` — genome-wide CNV/LOH/mosaic calling (spec CYTO-001 … CYTO-040).

I/O and orchestration only; the numerics live in ``simurg.cn``.
"""

from __future__ import annotations

import hashlib
import io
import json
import math
import zipfile
from dataclasses import dataclass, field
from pathlib import Path
from typing import Any

import numpy as np
import numpy.typing as npt
import pyarrow as pa
import pyarrow.parquet as pq

from simurg import __version__
from simurg.core.errors import ExitCode, SimurgError
from simurg.core.log import Reporter
from simurg.formats.bpm import Manifest, read_bpm
from simurg.formats.fasta import ContigMap, Genome

MODEL_SCHEMA = 1
F64 = npt.NDArray[np.float64]


def probe_set_hash(manifest: Manifest) -> str:
    h = hashlib.sha256()
    for name, chrom, pos in zip(manifest.names, manifest.chroms, manifest.map_infos.tolist(), strict=True):
        h.update(f"{name}\t{chrom}\t{pos}\n".encode())
    return h.hexdigest()


@dataclass
class CytoModel:
    header: dict[str, Any]
    names: list[str]
    gc: npt.NDArray[np.float32]  # (n, 3): 1 kb, 100 kb, 1 Mb
    pfb: npt.NDArray[np.float32]
    ref_lrr: npt.NDArray[np.float32] | None = None  # batch reference: per-probe median LRR
    ref_lrr_f: npt.NDArray[np.float32] | None = None  # same, X/Y from female reference samples
    ref_lrr_m: npt.NDArray[np.float32] | None = None  # same, X/Y from male reference samples
    ref_het: npt.NDArray[np.float32] | None = None  # per-probe median heterozygous BAF


def run_build_model(
    bpm: Path,
    fasta: Path,
    out: Path,
    rep: Reporter,
    egt: Path | None = None,
    contig_alias: str = "auto",
    reference_gtcs: Path | None = None,
    max_reference: int = 100,
) -> int:
    """CYTO-038: GC covariates (and PFB when a cluster file is given) per manifest probe → ``.simurg-cyto``."""
    from simurg.cn.model import GC_WINDOWS, gc_fractions, population_bfreq

    man = read_bpm(bpm)
    genome = Genome(fasta)
    cmap = ContigMap([e.name for e in genome.entries], contig_alias)
    chroms = np.array(man.chroms)
    pos = man.map_infos.astype(np.int64)
    gc = np.full((man.num_loci, len(GC_WINDOWS)), np.nan, dtype=np.float32)
    for chrom in sorted(set(man.chroms)):
        contig = cmap.get(chrom)
        if contig is None:
            continue
        sel = np.flatnonzero((chroms == chrom) & (pos > 0))
        if len(sel):
            gc[sel] = gc_fractions(genome._contig_array(contig), pos[sel])
            rep.info(f"GC windows for chromosome {chrom}: {len(sel)} probes")
    pfb = population_bfreq(None, man.num_loci)
    egt_name = ""
    if egt is not None:
        from simurg.formats.egt import read_egt
        from simurg.genotype.clusters import align_clusters

        clusters = align_clusters(read_egt(egt), man, allow_missing=True)
        pfb = population_bfreq(np.where(clusters.present[:, None], clusters.n, 0), man.num_loci)
        egt_name = Path(egt).name
    ref_cols: dict[str, Any] = {}
    n_ref = 0
    if reference_gtcs is not None:
        from simurg.cn.model import reference_stats
        from simurg.formats.discovery import discover_gtcs
        from simurg.formats.gtc import read_gtc

        paths = sorted(discover_gtcs(reference_gtcs).values())  # E_NO_INPUT when there is none
        if len(paths) > max_reference:
            idx = np.random.default_rng(42).choice(len(paths), max_reference, replace=False)
            paths = [paths[i] for i in sorted(idx)]
        if len(paths) < 24:
            rep.warn("W_QC_TOO_FEW", f"batch reference from only {len(paths)} GTCs (plan: >= 24)")
        lrr = np.empty((len(paths), man.num_loci), dtype=np.float32)
        baf = np.empty((len(paths), man.num_loci), dtype=np.float32)
        sexes = []
        for i, p in enumerate(paths):
            g = read_gtc(p)
            if len(g.logr_ratios) != man.num_loci:
                raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{p} has {len(g.logr_ratios)} loci, the manifest "
                                  f"{man.num_loci}")  # fmt: skip
            lrr[i], baf[i] = g.logr_ratios, g.b_allele_freqs
            sexes.append(g.gender.decode(errors="replace")[:1].upper() if g.gender else "U")
        ref_lrr, ref_het, n_het = reference_stats(lrr, baf)
        # X and Y from same-sex reference samples only (CYTO-005 baselines differ by sex)
        sx = np.isin(chroms, ["X", "XY", "Y"])
        ref_by_sex = {}
        for sex in ("F", "M"):
            rows = np.array([s_ == sex for s_ in sexes])
            col = ref_lrr.copy()
            if rows.sum() >= 3:
                col[sx] = reference_stats(lrr[rows][:, sx], baf[rows][:, sx])[0]
            ref_by_sex[sex] = col
        n_f, n_m = sexes.count("F"), sexes.count("M")
        rep.info(f"batch reference sexes: {n_f} F, {n_m} M, {len(sexes) - n_f - n_m} U")
        ref_cols = {"ref_lrr": ref_lrr, "ref_lrr_f": ref_by_sex["F"], "ref_lrr_m": ref_by_sex["M"],
                    "ref_het_baf": ref_het, "ref_n_het": n_het}  # fmt: skip
        n_ref = len(paths)
        rep.info(f"batch reference from {n_ref} GTC(s)")
    header = {
        "schema": "simurg.cyto_model",
        "schema_version": MODEL_SCHEMA,
        "software": f"SimurgArray {__version__}",
        "manifest": man.manifest_name,
        "probe_set_sha256": probe_set_hash(man),
        "n_probes": man.num_loci,
        "genome_build": genome.build,
        "fasta_sha256": genome.sha256,
        "gc_windows_bp": list(GC_WINDOWS),
        "pfb_source": f"cluster file {egt_name}" if egt_name else "none (0.5)",
        "reference_samples": n_ref,
    }
    table = pa.table({
        "name": man.names,
        "gc_1kb": gc[:, 0], "gc_100kb": gc[:, 1], "gc_1mb": gc[:, 2],
        "pfb": pfb,
        **ref_cols,
    })  # fmt: skip
    buf = io.BytesIO()
    pq.write_table(table, buf, compression="zstd")
    out.parent.mkdir(parents=True, exist_ok=True)
    tmp = out.with_suffix(out.suffix + ".part")
    with zipfile.ZipFile(tmp, "w", compression=zipfile.ZIP_STORED) as z:
        z.writestr("header.json", json.dumps(header, indent=1))
        z.writestr("probes.parquet", buf.getvalue())
    tmp.replace(out)
    rep.info(f"cyto model for {man.num_loci} probes written to {out}")
    return ExitCode.OK


def read_model(path: Path, manifest: Manifest) -> CytoModel:
    """IN-051: check schema and probe set."""
    with zipfile.ZipFile(path) as z:
        header = json.loads(z.read("header.json"))
        table = pq.read_table(io.BytesIO(z.read("probes.parquet")))
    if int(header.get("schema_version", 0)) < MODEL_SCHEMA:
        raise SimurgError("E_MODEL_VERSION", f"{path}: cyto model schema {header.get('schema_version')} is too old",
                          hint="rebuild it with simurg cn build-model")  # fmt: skip
    if header.get("probe_set_sha256") != probe_set_hash(manifest):
        raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{path} was built for manifest {header.get('manifest')!r}, "
                          f"not {manifest.manifest_name!r} (probe set differs)")  # fmt: skip
    gc = np.column_stack([table.column(c).to_numpy(zero_copy_only=False) for c in ("gc_1kb", "gc_100kb", "gc_1mb")])

    def col(c: str) -> npt.NDArray[np.float32] | None:
        return table.column(c).to_numpy(zero_copy_only=False).astype(np.float32) if c in table.column_names else None

    pfb = col("pfb")
    assert pfb is not None
    return CytoModel(header, table.column("name").to_pylist(), gc.astype(np.float32), pfb, ref_lrr=col("ref_lrr"),
                     ref_lrr_f=col("ref_lrr_f"), ref_lrr_m=col("ref_lrr_m"), ref_het=col("ref_het_baf"))  # fmt: skip


# ================================================================================================ cn call
PAR = {  # 1-based inclusive pseudo-autosomal regions on X
    "GRCh38": ((10_001, 2_781_479), (155_701_383, 156_030_895)),
    "GRCh37": ((60_001, 2_699_520), (154_931_044, 155_260_560)),
}
CALL_CHROMS = [str(i) for i in range(1, 23)] + ["X", "Y"]


@dataclass
class CnCallOptions:
    out: Path
    bpm: Path
    model: Path
    gtcs: Path | None = None
    sheet: Path | None = None
    summary: Path | None = None
    fasta: Path | None = None
    thresholds: Any = None  # simurg.cn.events.Thresholds
    smoothing: int = 5
    cn_cap: int = 4
    gc_correct: bool = True
    unknown_sex_policy: str = "skip-sex-chromosomes"
    gt_mode: str = "vcf"
    svlen_mode: str = "signed"
    mask_bed: Path | None = None
    bgzip: bool = True
    threads: int = 0
    force: bool = False
    argv: list[str] = field(default_factory=list)


def _builtin_file(name: str) -> bytes:
    from importlib import resources

    return resources.files("simurg.cn").joinpath("data", name).read_bytes()


def read_mask(path: Path | None, build: str) -> dict[str, list[tuple[int, int]]]:
    """CYTO-017: BED (0-based half-open) → per-chromosome intervals as 1-based inclusive."""
    text = path.read_text() if path is not None else _builtin_file(f"mask_{build}.bed").decode()
    out: dict[str, list[tuple[int, int]]] = {}
    for line in text.splitlines():
        if not line.strip() or line.startswith(("#", "track", "browser")):
            continue
        c, s, e = line.split("\t")[:3]
        c = c[3:] if c.lower().startswith("chr") else c
        out.setdefault(c, []).append((int(s) + 1, int(e)))
    return out


def moving_median(x: F64, k: int) -> F64:
    """CYTO-004: median over 2k+1 consecutive values, truncated (never wrapped) at the ends; NaN ignored."""
    if k <= 0 or len(x) == 0:
        return x.copy()
    from numpy.lib.stride_tricks import sliding_window_view

    pad = np.concatenate([np.full(k, np.nan), x, np.full(k, np.nan)])
    with np.errstate(all="ignore"):
        import warnings

        with warnings.catch_warnings():
            warnings.simplefilter("ignore", RuntimeWarning)
            med = np.nanmedian(sliding_window_view(pad, 2 * k + 1), axis=1)
    return np.where(np.isfinite(x), med, np.nan)


@dataclass
class SampleInput:
    sample_id: str
    sample_name: str
    gtc: str
    sex: str  # M | F | U
    logr_dev: float
    sex_source: str


def _robust_sd(x: F64, floor: float) -> float:
    x = x[np.isfinite(x)]
    if not len(x):
        return floor
    return max(floor, float(1.4826 * np.median(np.abs(x - np.median(x)))))


def _resolve_sex(si: SampleInput, chroms: npt.NDArray[np.str_], lrr: F64, geno: npt.NDArray[np.uint8],
                 policy: str) -> tuple[str, str, list[str]]:  # fmt: skip
    """CYTO-024."""
    if si.sex in ("M", "F"):
        return si.sex, si.sex_source, []
    x = chroms == "X"
    called = geno[x] > 0
    het = float((geno[x] == 2).sum() / called.sum()) if called.any() else math.nan
    med = float(np.nanmedian(lrr[x])) if x.any() else math.nan
    if math.isfinite(het) and math.isfinite(med):
        if het < 0.03 and med < -0.15:
            return "M", "resolved", []
        if het > 0.08 and med > -0.1:
            return "F", "resolved", []
    if policy == "female":
        return "F", "assumed", ["sex_assumed"]
    if policy == "male":
        return "M", "assumed", ["sex_assumed"]
    if policy == "error":
        raise SimurgError("E_CN_SEX_UNRESOLVED", f"{si.sample_id}: sex could not be resolved")
    return "U", "unresolved", ["sex_unresolved"]


def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> dict[str, Any]:
    """One sample: preprocessing, per-region HMM, events, VCF + segments. Returns the summary row + events."""
    from simurg.cn import hmm
    from simurg.cn.events import (
        annotate_stats,
        apply_thresholds,
        effective_size,
        extract,
        median_spacing,
        post_filter,
        set_filters,
    )
    from simurg.cn.gcwave import gc_correct
    from simurg.formats.gtc import read_gtc

    man_chroms: npt.NDArray[np.str_] = ctx["chroms"]
    pos_all: npt.NDArray[np.int64] = ctx["pos"]
    model: CytoModel = ctx["model"]
    th = opts.thresholds
    g = read_gtc(si.gtc)
    lrr = g.logr_ratios.astype(np.float64)
    baf = g.b_allele_freqs.astype(np.float64)
    geno = g.genotypes
    lrr[~np.isfinite(lrr)] = np.nan
    baf[~np.isfinite(baf)] = np.nan
    sex, sex_source, flags = _resolve_sex(si, man_chroms, lrr - float(np.nanmedian(lrr[ctx["auto"]])), geno,
                                          opts.unknown_sex_policy)  # fmt: skip
    ref = {"F": model.ref_lrr_f, "M": model.ref_lrr_m}.get(sex) if model.ref_lrr_f is not None else None
    if ref is None:
        ref = model.ref_lrr
    if ref is not None:
        lrr = lrr - ref
    if model.ref_het is not None:
        het = (baf > 0.1) & (baf < 0.9)
        baf = np.where(het, np.clip(baf - (model.ref_het - 0.5), 0.0, 1.0), baf)
    auto = ctx["auto"]
    gcfit = None
    if opts.gc_correct:
        lrr, gcfit = gc_correct(lrr, model.gc, auto)
    lrr = lrr - float(np.nanmedian(lrr[auto]))
    logr_dev = si.logr_dev if math.isfinite(si.logr_dev) else float(np.nanstd(lrr[auto]))
    # CYTO-014 gates on the noise of the LRR the caller actually uses (after the batch reference and GC correction):
    # the summary's logr_dev also measures cluster-file mismatch, which the reference removes (K-027)
    gate_sd = float(np.nanstd(lrr[auto])) if model.ref_lrr is not None else logr_dev
    noisy = gate_sd > 0.2
    whole_chrom_only = gate_sd > 0.3
    if noisy:
        flags.append("high_noise")
    if whole_chrom_only:
        flags.append("whole_chromosome_only")
    pfb = model.pfb.astype(np.float64)
    lrr_s = np.full_like(lrr, np.nan)
    for sel in ctx["by_chrom"].values():
        lrr_s[sel] = moving_median(lrr[sel], opts.smoothing)
    # emission SD from the unsmoothed LRR: smoothed values are autocorrelated, their own SD would overstate evidence
    seg_sd = _robust_sd(lrr[auto], 0.05)
    b_het = baf[auto & (baf > 0.3) & (baf < 0.7)]
    baf_sd = _robust_sd(b_het, 0.02)
    dip, hap = hmm.diploid_states(), hmm.haploid_states()
    counts: dict[str, int] = {}
    events: list[Any] = []
    segments: list[dict[str, Any]] = []
    ploidy_num = 0.0
    ploidy_den = 0
    par = PAR.get(ctx["build"], PAR["GRCh38"])
    for chrom in CALL_CHROMS:
        if chrom in ("X", "Y") and sex == "U":
            continue
        if chrom == "Y" and sex == "F":
            continue
        sel = ctx["by_chrom"].get(chrom)
        if sel is None or len(sel) == 0:
            continue
        pos = pos_all[sel]
        in_par = np.zeros(len(sel), dtype=np.bool_)
        if chrom == "X":
            for a, b in par:
                in_par |= (pos >= a) & (pos <= b)
        regions = [(np.ones(len(sel), dtype=np.bool_), False)]
        if chrom == "X" and sex == "M":
            regions = [(in_par, False), (~in_par, True)]
        elif chrom == "Y":
            regions = [(np.ones(len(sel), dtype=np.bool_), True)]
        m_sp = median_spacing(pos)
        for rmask, haploid in regions:
            if rmask.sum() < 2:
                continue
            idx = sel[rmask]
            p, x_raw, b, gt = pos_all[idx], lrr[idx], baf[idx], geno[idx]
            x_seg = lrr_s[idx]
            states = hap if haploid else dip
            prior = np.array([s.prior for s in states])
            prior /= prior.sum()
            offset = float(np.nanmedian(x_raw)) if haploid else 0.0
            if not math.isfinite(offset):
                offset = 0.0
            prm = hmm.Params(seg_sd, baf_sd, lrr_offset=offset)
            e = hmm.emissions(states, x_seg, np.full(len(b), np.nan) if haploid else b, pfb[idx], prm)
            sw = hmm.switch_probs(p)
            path = hmm.viterbi(e, sw, prior)
            post = hmm.forward_backward(e, sw, prior)
            evs = extract(chrom, path, states, p)
            for ev in evs:
                annotate_stats(ev, p, x_raw, b, gt, post, states, m_sp, ctx["gamma"], haploid)
            evs = post_filter(evs, [x_raw[ev.i0 : ev.i1] for ev in evs], th, counts)
            for s_, e_ in _runs(path):
                st = states[int(path[s_])]
                segments.append({"chrom": chrom, "start": int(p[s_]), "end": int(p[e_ - 1]), "probes": e_ - s_,
                                 "state": st.name, "call": st.call, "cn": st.cn, "haploid": haploid,
                                 "posterior": float(post[s_:e_, int(path[s_])].mean())})  # fmt: skip
            base_cn = 1 if haploid else 2
            eff = effective_size(p, m_sp)
            if whole_chrom_only:
                evs = [ev for ev in evs if ev.effective_size >= 0.9 * eff]
            ploidy_num += base_cn * eff + sum((ev.cn - base_cn) * ev.effective_size for ev in evs)
            ploidy_den += eff
            events.extend((ev, haploid) for ev in evs)
    kept = apply_thresholds([ev for ev, _ in events], th, counts)
    kept = _mask_events(kept, ctx["mask"], counts)
    set_filters(kept, th, noisy)
    hap_of = {id(ev): h for ev, h in events}
    row = {
        "sample_id": si.sample_id, "sample_name": si.sample_name, "status": "OK", "logr_dev": logr_dev,
        "overall_ploidy": ploidy_num / max(ploidy_den, 1), "gc_correct": bool(gcfit and gcfit.applied),
        "gcwf": gcfit.gcwf if gcfit else math.nan, "sex_used": sex, "sex_source": sex_source,
        "n_del": sum(e.call == "DEL" for e in kept), "n_dup": sum(e.call == "DUP" for e in kept),
        "n_loh": sum(e.call == "LOH" for e in kept), "n_gainloh": sum(e.call == "GAINLOH" for e in kept),
        "n_mosaic": sum(e.mosaic for e in kept), "n_filtered": sum(counts.values()),
        "qc_flags": ";".join(flags),
    }  # fmt: skip
    per_chr: dict[str, int] = {}
    for ev in kept:
        if ev.call in ("DEL", "DUP", "LOH"):
            per_chr[ev.chrom] = per_chr.get(ev.chrom, 0) + 1
    excess = [c for c, n in per_chr.items() if n > 10]
    if excess:
        row["qc_flags"] = ";".join([*flags, "excess_events:" + ",".join(excess)])
    row["lrr_sd_segmentation"] = seg_sd
    row["lrr_sd_gate"] = gate_sd
    row["baf_sd"] = baf_sd
    # CYTO-033 sample and chromosome statistics (unsmoothed LRR), for cn annotate
    la = lrr[auto & np.isfinite(lrr)]
    row["logr_dev_robust"] = _robust_sd(la, 0.0)
    row["call_rate"] = float(g.call_rate)
    outside = np.ones(len(lrr), dtype=np.bool_)
    for ev in kept:
        sel = ctx["by_chrom"].get(ev.chrom, np.zeros(0, np.int64))
        p_ = pos_all[sel]
        outside[sel[(p_ >= ev.start) & (p_ <= ev.end)]] = False
    bdev = {}
    for key, code in (("aa", 1), ("ab", 2), ("bb", 3)):
        v = baf[auto & outside & (geno == code) & np.isfinite(baf)]
        bdev[key] = float(np.std(v)) if len(v) > 1 else None
    row["baf_dev"] = bdev
    chrom_stats = []
    for chrom in CALL_CHROMS:
        sel = ctx["by_chrom"].get(chrom)
        if sel is None or len(sel) == 0:
            continue
        x = lrr[sel][np.isfinite(lrr[sel])]
        gt = geno[sel]
        called = gt > 0
        chrom_stats.append({
            "id": chrom, "size": ctx["contig_len"].get(chrom, int(pos_all[sel].max())),
            "effective_size": effective_size(pos_all[sel], median_spacing(pos_all[sel])),
            "percent_het": float(100 * (gt == 2).sum() / called.sum()) if called.any() else None,
            "lrr_median": float(np.median(x)) if len(x) else None, "lrr_mean": float(np.mean(x)) if len(x) else None,
            "lrr_dev": float(np.std(x)) if len(x) > 1 else None, "min_lrr": float(x.min()) if len(x) else None,
            "max_lrr": float(x.max()) if len(x) else None,
        })  # fmt: skip
    row["chromosomes"] = chrom_stats
    return {"row": row, "events": kept, "haploid": hap_of, "segments": segments, "counts": counts}


def _runs(path: npt.NDArray[np.int32]) -> list[tuple[int, int]]:
    if len(path) == 0:
        return []
    br = np.flatnonzero(np.diff(path) != 0) + 1
    return list(zip(np.r_[0, br].tolist(), np.r_[br, len(path)].tolist(), strict=True))


def _mask_events(events: list[Any], mask: dict[str, list[tuple[int, int]]], counts: dict[str, int]) -> list[Any]:
    """CYTO-017: trim to unmasked sequence; drop events ≥ 50 % masked."""
    out = []
    for ev in events:
        covered = 0
        for a, b in mask.get(ev.chrom, []):
            lo, hi = max(a, ev.start), min(b, ev.end)
            if lo <= hi:
                covered += hi - lo + 1
                if a <= ev.start <= b:
                    ev.start = b + 1
                    ev.flags.append("centromere_trimmed")
                if a <= ev.end <= b:
                    ev.end = a - 1
                    ev.flags.append("centromere_trimmed")
        if covered >= 0.5 * (ev.end - ev.start + 1 + covered) or ev.end < ev.start:
            counts["masked"] = counts.get("masked", 0) + 1
            continue
        out.append(ev)
    return out


# ------------------------------------------------------------------------------------------------ VCF (CYTO-028..031)
def _gt(ev: Any, haploid: bool, mode: str) -> str:
    if mode == "compat":
        return "1/1"
    if haploid:
        return "1"
    if ev.mosaic:
        return "0/1"
    if ev.call in ("LOH", "GAINLOH"):
        return "1/1"
    if ev.call == "DEL":
        return "1/1" if ev.cn == 0 else "0/1"
    return "./1" if ev.cn >= 4 else "0/1"


def _fmt(v: float | None, nd: int) -> str:
    return "." if v is None or not math.isfinite(v) else f"{v:.{nd}f}"


def vcf_text(
    sample_name: str,
    events: list[Any],
    haploid: dict[int, bool],
    contigs: list[tuple[str, int]],
    header: dict[str, str],
    opts: CnCallOptions,
    ref_base: Any,
) -> tuple[str, list[tuple[str, int, int, str]]]:
    """Header and records; also (chrom, beg0, end, line) for indexing."""
    th = opts.thresholds
    extra = {
        "smoothing": opts.smoothing,
        "max_mosaic_fraction": th.max_mosaic_fraction,
        "min_mosaic_fraction": f"{th.min_mosaic_fraction:.2f}",
        "filter_qual": th.filter_qual,
        "size_basis": "effective",
        "cn_cap": opts.cn_cap,
        "gt_mode": opts.gt_mode,
        "svlen_mode": opts.svlen_mode,
    }
    h = [
        "##fileformat=VCFv4.3",
        f"##source=SimurgArray {__version__} cn",
        f"##reference=file://{header['reference']}",
        f"##simurg_genome_build={header['build']}",
        f"##simurg_manifest={header['manifest']}",
        f"##simurg_cyto_model={header['model']}",
        "##simurg_thresholds=" + th.as_header(extra),
        f"##simurg_sex_used={header['sex']};source={header['sex_source']}",
        f"##OverallPloidy={float(header['ploidy']):.3f}",
        f"##GCCorrect={'true' if header['gc'] == 'True' else 'false'}",
        *(f"##contig=<ID={c},length={n}>" for c, n in contigs),
        '##ALT=<ID=DEL,Description="Copy number loss">',
        '##ALT=<ID=DUP,Description="Copy number gain">',
        '##ALT=<ID=LOH,Description="Absence/loss of heterozygosity (AOH/LOH/ROH)">',
        '##FILTER=<ID=PASS,Description="All filters passed">',
        f'##FILTER=<ID=Q{th.filter_qual},Description="Quality below {th.filter_qual}">',
        '##FILTER=<ID=SampleQuality,Description="Sample flagged as potentially low quality: LRR SD above 0.2">',
        f'##FILTER=<ID=LowFraction,Description="Mosaic event with estimated mosaic fraction below '
        f'{th.min_mosaic_fraction:.2f}">',
        '##INFO=<ID=END,Number=1,Type=Integer,Description="Last base of the event (1-based, inclusive)">',
        '##INFO=<ID=SVTYPE,Number=1,Type=String,Description="DEL, DUP or LOH">',
        '##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Event length in bp; negative for deletions">',
        '##INFO=<ID=EFFLEN,Number=1,Type=Integer,Description="Effective length in bp (gaps above 150x median probe '
        'spacing replaced by the median spacing)">',
        '##INFO=<ID=PROBE,Number=1,Type=Integer,Description="Number of probes in the event">',
        '##INFO=<ID=MOSAIC,Number=0,Type=Flag,Description="Mosaic event">',
        '##INFO=<ID=HIGHFRACTION,Number=0,Type=Flag,Description="Mosaic event promoted to constitutional because its '
        'mosaic fraction reached the maximum">',
        '##INFO=<ID=LOHTYPE,Number=1,Type=String,Description="AOH (constitutional, copy-neutral or gain), CNLOH '
        '(mosaic copy-neutral), GAINLOH (mosaic copy-gain)">',
        '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype (see SimurgArray docs, CYTO-031)">',
        f'##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number; the maximum value {opts.cn_cap} means '
        f'{opts.cn_cap} or more">',
        '##FORMAT=<ID=NR,Number=1,Type=Float,Description="Aggregated normalized intensity (median 2^LRR; about CN/2)">',
        '##FORMAT=<ID=LRD,Number=1,Type=Float,Description="Standard deviation of LRR within the event">',
        '##FORMAT=<ID=MF,Number=1,Type=Float,Description="Estimated mosaic fraction; missing if not mosaic or not '
        'estimable">',
        f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}",
    ]
    order = {c: i for i, (c, _) in enumerate(contigs)}
    recs = []
    for ev in sorted(events, key=lambda e: (order.get(e.chrom, 999), e.start, e.end)):
        s_, e_ = ev.start, ev.end
        pos = s_ - 1 if s_ > 1 else 1
        ref = ref_base(ev.chrom, pos)
        length = e_ - s_ + 1
        svlen = -length if ev.svtype == "DEL" and opts.svlen_mode == "signed" else length
        info = [f"END={e_}", f"SVTYPE={ev.svtype}", f"SVLEN={svlen}", f"EFFLEN={ev.effective_size}",
                f"PROBE={ev.probes}"]  # fmt: skip
        if ev.mosaic:
            info.append("MOSAIC")
        if ev.high_fraction:
            info.append("HIGHFRACTION")
        lt = {"LOH": "AOH", "GAINLOH": "AOH", "MLOH": "CNLOH"}.get(ev.call)
        if ev.call == "GAINLOH" and ev.mosaic:
            lt = "GAINLOH"
        if lt:
            info.append(f"LOHTYPE={lt}")
        cn = min(ev.cn, opts.cn_cap)
        mf = ev.mosaic_fraction if ev.mosaic or ev.high_fraction else None
        line = (f"{ev.chrom}\t{pos}\t{ev.id_prefix}:{ev.chrom}:{s_}-{e_}\t{ref}\t<{ev.svtype}>\t{int(ev.quality)}\t"
                f"{';'.join(ev.filters)}\t{';'.join(info)}\tGT:CN:NR:LRD:MF\t"
                f"{_gt(ev, haploid.get(id(ev), False), opts.gt_mode)}:{cn}:{_fmt(ev.nr, 4)}:{_fmt(ev.lrd, 4)}:"
                f"{_fmt(mf, 3)}")  # fmt: skip
        recs.append((ev.chrom, pos - 1, e_, line))
    return "\n".join(h) + "\n", recs


def write_vcf(path: Path, header: str, recs: list[tuple[str, int, int, str]], bgzip: bool) -> None:
    from simurg.core.atomic import atomic_open
    from simurg.formats.bgzf import BgzfWriter, TabixBuilder, write_tbi

    tbx = TabixBuilder()
    with atomic_open(path, "wb") as raw:
        if bgzip:
            w = BgzfWriter(raw)  # type: ignore[arg-type]
            w.write_record(header.encode())
            for chrom, beg, end, line in recs:
                vs, ve = w.write_record((line + "\n").encode())
                tbx.add(chrom, beg, end, vs, ve)
            w.close()
        else:
            raw.write((header + "".join(r[3] + "\n" for r in recs)).encode())  # type: ignore[arg-type]
    if bgzip:
        write_tbi(Path(str(path) + ".tbi"), tbx)


# ------------------------------------------------------------------------------------------------ batch driver
_W: dict[str, Any] = {}


def _init_worker(path: str) -> None:
    import pickle

    with open(path, "rb") as fh:
        _setup(pickle.load(fh))


def _setup(state: dict[str, Any]) -> None:
    _W.update(state)
    fasta = state["opts"].fasta
    if fasta is not None:
        genome = Genome(fasta)
        _W["genome"] = genome
        _W["cmap"] = ContigMap([e.name for e in genome.entries], "auto")


def _work(si: SampleInput) -> dict[str, Any]:
    ctx, opts = _W["ctx"], _W["opts"]
    try:
        res = call_sample(si, ctx, opts)
    except Exception as e:  # CYTO-040: one sample never stops the batch
        code = e.code if isinstance(e, SimurgError) else "E_CN_FAILED"
        return {"row": {"sample_id": si.sample_id, "sample_name": si.sample_name, "status": "FAILED",
                        "qc_flags": f"{code}: {e}"}, "events": [], "error": code}  # fmt: skip
    row = res["row"]
    genome = _W.get("genome")

    def ref_base(chrom: str, pos: int) -> str:
        if genome is None:
            return "N"
        contig = _W["cmap"].get(chrom)
        return genome.fetch(contig, pos - 1, pos) or "N" if contig else "N"

    header = {**_W["vcf_header"], "sex": row["sex_used"], "sex_source": row["sex_source"],
              "ploidy": str(row["overall_ploidy"]), "gc": str(row["gc_correct"])}  # fmt: skip
    htxt, recs = vcf_text(si.sample_name, res["events"], res["haploid"], ctx["contigs"], header, opts, ref_base)
    out = opts.out / f"{si.sample_id}.cyto.cnv.vcf{'.gz' if opts.bgzip else ''}"
    write_vcf(out, htxt, recs, opts.bgzip)
    seg = pa.Table.from_pylist(res["segments"]) if res["segments"] else pa.table({"chrom": pa.array([], pa.string())})
    pq.write_table(seg, opts.out / f"{si.sample_id}.cn_segments.parquet", compression="zstd")
    evs = [
        {"id": f"{ev.id_prefix}:{ev.chrom}:{ev.start}-{ev.end}", "chrom": ev.chrom, "start": ev.start, "end": ev.end,
         "call_type": ev.svtype, "call": ev.call, "loh_type": ev.loh_type, "mosaic": ev.mosaic,
         "mosaic_fraction": ev.mosaic_fraction, "high_fraction": ev.high_fraction, "copy_number": ev.cn,
         "quality": round(ev.quality, 1), "filters": ev.filters, "size": ev.size, "effective_size": ev.effective_size,
         "probe_count": ev.probes, "percent_het": ev.percent_het, "lrr_median": ev.lrr_median, "lrr_dev": ev.lrd,
         "baf_dev": ev.baf_dev, "nr": ev.nr, "flags": ev.flags, "base_cn": 1 if res["haploid"].get(id(ev)) else 2}
        for ev in res["events"]
    ]  # fmt: skip
    return {"row": row, "events": evs, "counts": res["counts"]}


def _read_summary(path: Path | None) -> dict[str, dict[str, str]]:
    if path is None or not path.exists():
        return {}
    lines = path.read_text(encoding="utf-8").splitlines()
    head = lines[0].split("\t")
    return {r[0]: dict(zip(head, r, strict=False)) for r in (ln.split("\t") for ln in lines[1:] if ln)}


def run_cn_call(opts: CnCallOptions, rep: Reporter) -> int:
    """``simurg cn call`` (CLI-028)."""
    import multiprocessing as mp
    import os
    import pickle
    import tempfile
    from concurrent.futures import ProcessPoolExecutor

    from simurg.core.atomic import atomic_write_text, check_outputs_absent
    from simurg.core.provenance import RunManifest
    from simurg.formats.discovery import discover_gtcs
    from simurg.pipeline.qc import write_tsv

    rm = RunManifest(command="cn call", argv=opts.argv)
    rep.run_id = rm.run_id
    if opts.gtcs is None:
        raise SimurgError("E_NO_INPUT", "give --gtcs (folder of GTC files)")
    gtcs = discover_gtcs(opts.gtcs)
    summary_path = opts.summary or (opts.gtcs / "sample_summary.tsv")
    summary = _read_summary(summary_path)
    ext = ".cyto.cnv.vcf.gz" if opts.bgzip else ".cyto.cnv.vcf"
    check_outputs_absent(
        [opts.out / f"{k}{ext}" for k in gtcs] + [opts.out / "cn_sample_summary.tsv"], force=opts.force
    )
    th = opts.thresholds
    for k in list(th.probes):
        if th.probes[k] < 10:
            rep.warn("W_CN_PROBE_FLOOR", f"--min-{k}-probes {th.probes[k]} raised to 10")
            th.probes[k] = 10
    man = read_bpm(opts.bpm)
    model = read_model(opts.model, man)
    build = str(model.header.get("genome_build") or "GRCh38")
    chroms = np.array(["X" if c == "XY" else (c[3:] if c.lower().startswith("chr") else c) for c in man.chroms])
    pos = man.map_infos.astype(np.int64)
    by_chrom = {}
    for c in CALL_CHROMS:
        sel = np.flatnonzero((chroms == c) & (pos > 0))
        by_chrom[c] = sel[np.argsort(pos[sel], kind="stable")]
    auto = np.isin(chroms, CALL_CHROMS[:22]) & (pos > 0)
    genome = None
    cmap = None
    lengths: dict[str, int] = {}
    if opts.fasta is not None:
        genome = Genome(opts.fasta)
        cmap = ContigMap([e.name for e in genome.entries], "auto")
        for c in CALL_CHROMS:
            contig = cmap.get(c)
            if contig:
                lengths[c] = genome.contig_length(contig)
    else:
        rep.warn("W_CN_NO_FASTA", "no --fasta: VCF REF padding bases are written as N")
    contigs = [(c, lengths.get(c, int(pos[by_chrom[c]].max()) if len(by_chrom[c]) else 0)) for c in CALL_CHROMS
               if len(by_chrom[c])]  # fmt: skip
    ctx = {"chroms": chroms, "pos": pos, "model": model, "auto": auto, "by_chrom": by_chrom, "build": build,
           "mask": read_mask(opts.mask_bed, build), "gamma": 0.5, "contigs": contigs,
           "contig_len": dict(contigs)}  # fmt: skip
    model_sha = hashlib.sha256(opts.model.read_bytes()).hexdigest()
    man_sha = hashlib.sha256(opts.bpm.read_bytes()).hexdigest()
    schema = model.header.get("schema_version")
    vcf_header = {
        "reference": str(opts.fasta.resolve()) if opts.fasta else "none",
        "build": build,
        "manifest": f"{man.manifest_name};sha256={man_sha}",
        "model": f"{opts.model.name};sha256={model_sha};schema={schema}",
    }
    samples = []
    for k, path in gtcs.items():
        r = summary.get(k, {})
        sex = r.get("sex_est", "U") or "U"
        src = "summary" if sex in ("M", "F") else "none"
        if sex not in ("M", "F") and r.get("sex_reported") in ("M", "F"):
            sex, src = r["sex_reported"], "sheet"
        try:
            ld = float(r.get("logr_dev", "nan"))
        except ValueError:
            ld = math.nan
        samples.append(SampleInput(k, r.get("sample_name", k) or k, str(path), sex, ld, src))
    opts.out.mkdir(parents=True, exist_ok=True)
    state = {"ctx": ctx, "opts": opts, "vcf_header": vcf_header}
    workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(samples)))
    rep.info(f"CN calling of {len(samples)} sample(s) with {workers} worker(s)")
    if workers == 1:
        _setup(state)
        try:
            results = [_work(s) for s in samples]
        finally:
            _W.clear()
    else:
        fd, sp = tempfile.mkstemp(prefix=".simurg-workers-", suffix=".pkl", dir=opts.out)
        try:
            with os.fdopen(fd, "wb") as fh:
                pickle.dump(state, fh, protocol=5)
            with ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn"), initializer=_init_worker,
                                     initargs=(sp,)) as ex:  # fmt: skip
                results = list(ex.map(_work, samples))
        finally:
            Path(sp).unlink(missing_ok=True)
    cols = ["sample_id", "sample_name", "status", "logr_dev", "overall_ploidy", "gc_correct", "gcwf", "sex_used",
            "sex_source", "n_del", "n_dup", "n_loh", "n_gainloh", "n_mosaic", "n_filtered", "qc_flags"]  # fmt: skip
    rows = []
    for s, res in zip(samples, results, strict=True):
        row = {**res["row"], "sample_folder": str(Path(s.gtc).parent)}
        if row["status"] != "OK":
            rep.error(res.get("error", "E_CN_FAILED"), f"{s.sample_id}: {row['qc_flags']}", sample_id=s.sample_id)
        elif row["sex_source"] == "resolved":
            rep.info(f"Sample {s.sample_id} sex updated from U to {row['sex_used']}", "I_CN_SEX_RESOLVED")
        rows.append(row)
    write_tsv(opts.out / "cn_sample_summary.tsv", cols, rows)
    th_dict = {f"min_{k}_probes": th.probes[k] for k in th.probes} | {
        f"min_{k}_size_bp": th.size_bp[k] for k in th.size_bp
    }
    sidecar = {
        "schema": "simurg.cn_call/1", "software": f"SimurgArray {__version__}", "run_id": rm.run_id,
        "genome_build": build, "model": vcf_header["model"], "manifest": vcf_header["manifest"],
        "reference": vcf_header["reference"],
        "summary": str(summary_path) if summary else None,
        "thresholds": th_dict | {"max_mosaic_fraction": th.max_mosaic_fraction,
                                 "min_mosaic_fraction": th.min_mosaic_fraction, "filter_qual": th.filter_qual,
                                 "smoothing": opts.smoothing, "cn_cap": opts.cn_cap, "size_basis": "effective"},
        "samples": [{**r["row"], "events": r["events"]} for r in results],
    }  # fmt: skip
    atomic_write_text(opts.out / "cn_call_manifest.json", json.dumps(sidecar, indent=1, default=_json_default))
    ok = sum(r["row"]["status"] == "OK" for r in results)
    rm.summary = {"samples": len(results), "ok": ok}
    rm.status = "completed"
    rm.write(opts.out / "run_manifest.json")
    rep.info(f"done: {ok}/{len(results)} samples; {sum(len(r['events']) for r in results)} events")
    if ok == len(results):
        return ExitCode.OK
    return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED


def _json_default(v: Any) -> Any:
    if isinstance(v, np.generic):
        return v.item()
    if isinstance(v, float) and not math.isfinite(v):
        return None
    raise TypeError(type(v))


# ================================================================================================ cn annotate
@dataclass
class CnAnnotateOptions:
    vcfs: Path
    bundle: Path | None
    out: Path | None = None
    min_qual: int = 20
    size_basis: str = "effective"
    iscn_min_size: int = 0
    tighten: dict[str, int] = field(default_factory=dict)  # min_<type>_probes / min_<type>_size_bp overrides
    force: bool = False


def _load_bundle(path: Path | None, build: str) -> tuple[Any, Any, dict[str, Any] | None]:
    import gzip

    from simurg.cn.annotate import Bands, Genes

    if path is None:
        text = gzip.decompress(_builtin_file(f"cytoBand_{build}.txt.gz")).decode()
        return Bands.from_ucsc(text), None, None
    with zipfile.ZipFile(path) as z:
        man = json.loads(z.read("manifest.json"))
        bands = Bands.from_ucsc(gzip.decompress(z.read("cytoBand.txt.gz")).decode())
        rows = []
        for line in gzip.decompress(z.read("genes.tsv.gz")).decode().splitlines()[1:]:
            c, s_, e_, gene, tx = line.split("\t")
            rows.append((c, int(s_), int(e_), gene or tx))
    man["file"] = path.name
    man["sha256"] = hashlib.sha256(path.read_bytes()).hexdigest()
    return bands, Genes.from_rows(rows), man


def run_cn_annotate(opts: CnAnnotateOptions, rep: Reporter) -> int:
    """``simurg cn annotate`` (CLI-029): per-sample ``<sample>.cyto.json`` with ISCN, bands and genes."""
    from simurg.cn.annotate import CHROM_ORDER, iscn_microarray_event, iscn_sample, iscn_traditional_event
    from simurg.core.atomic import atomic_write_text, check_outputs_absent

    side_p = opts.vcfs / "cn_call_manifest.json"
    if not side_p.exists():
        raise SimurgError("E_CN_SIDECAR_MISSING", f"{side_p} not found; run simurg cn call first")
    side = json.loads(side_p.read_text())
    out = opts.out or opts.vcfs
    samples = [s for s in side["samples"] if s.get("status") == "OK"]
    for s in samples:
        v = opts.vcfs / f"{s['sample_id']}.cyto.cnv.vcf.gz"
        if not (v.exists() and Path(str(v) + ".tbi").exists()):
            raise SimurgError("E_CN_VCF_UNINDEXED", f"{v} (+ .tbi) not found; cn call must run without --no-bgzip")
    check_outputs_absent([out / f"{s['sample_id']}.cyto.json" for s in samples], force=opts.force)
    th = dict(side["thresholds"])
    for k, val in opts.tighten.items():  # CLI-029: can only tighten
        if val is not None and val > th.get(k, 0):
            th[k] = val
    build = side.get("genome_build", "GRCh38")
    bands, genes, bundle = _load_bundle(opts.bundle, build)
    if bundle is None:
        rep.warn("W_CN_NO_FASTA", "no --annotation-bundle: built-in cytobands only, no gene lists")
    key = {"DEL": "del", "DUP": "dup", "LOH": "loh", "GAINLOH": "gainloh", "MDEL": "mdel", "MDUP": "mdup",
           "MLOH": "mloh"}  # fmt: skip
    out.mkdir(parents=True, exist_ok=True)
    for s in samples:
        events = []
        for ev in s["events"]:
            k = key[ev["call"]]
            size = ev["effective_size"] if opts.size_basis == "effective" else ev["size"]
            if ev["probe_count"] < th[f"min_{k}_probes"] or size < th[f"min_{k}_size_bp"]:
                continue
            if ev["quality"] < opts.min_qual:
                continue
            e = dict(ev)
            e["start_cytoband"], e["end_cytoband"] = (
                bands.band(ev["chrom"], ev["start"]),
                bands.band(ev["chrom"], ev["end"]),
            )
            e["iscn_microarray"] = f"arr[{build}] " + iscn_microarray_event(ev, bands)
            e["iscn_traditional"] = iscn_traditional_event(ev, bands)
            gl = genes.overlapping(ev["chrom"], ev["start"], ev["end"]) if genes else []
            e["genes"], e["gene_count"] = gl, len(gl)
            e["call_type"], e["loh_type"] = ev["call_type"], ev["loh_type"]
            events.append(e)
        reported = [e for e in events if e["filters"] == ["PASS"] and e["size"] >= opts.iscn_min_size]
        trad, micro = iscn_sample(reported, bands, build, s["sex_used"])
        eff_genome = sum(c["effective_size"] for c in s.get("chromosomes", []))
        loh = [e for e in events if e["call_type"] == "LOH"]

        def loh_stats(evs: list[dict[str, Any]], denom: int) -> dict[str, Any]:
            big = [e for e in evs if e["size"] > 1_000_000]
            tot = sum(e["effective_size"] for e in big)
            pct = round(100 * tot / denom, 3) if denom else None
            return {"num_loh_over_1mb": len(big), "num_loh_over_8mb": sum(e["size"] > 8_000_000 for e in evs),
                    "total_size_loh_over_1mb": tot, "percent_loh": pct}  # fmt: skip

        chroms = []
        for c in s.get("chromosomes", []):
            ce = [e for e in events if e["chrom"] == c["id"]]
            base = 1 if (c["id"] == "Y" or (c["id"] == "X" and s["sex_used"] == "M")) else 2
            cov = sum(e["effective_size"] for e in ce if e["call_type"] != "LOH")
            mean_cn = (base * max(c["effective_size"] - cov, 0) + sum(e["copy_number"] * e["effective_size"] for e in ce
                       if e["call_type"] != "LOH")) / max(c["effective_size"], 1)  # fmt: skip
            mfs = [e["mosaic_fraction"] for e in ce if e.get("mosaic") and e.get("mosaic_fraction") is not None]
            chroms.append({**c, "has_mosaicism": any(e.get("mosaic") for e in ce),
                           **loh_stats([e for e in loh if e["chrom"] == c["id"]], c["effective_size"]),
                           "copy_number_median": float(base if cov < c["effective_size"] / 2 else
                                                       max((e["copy_number"] for e in ce), default=base)),
                           "copy_number_mean": round(mean_cn, 4),
                           "median_mosaic_fraction": float(np.median(mfs)) if mfs else None,
                           "number_del": sum(e["call_type"] == "DEL" and not e["mosaic"] for e in ce),
                           "number_dup": sum(e["call_type"] == "DUP" and not e["mosaic"] for e in ce),
                           "number_loh": sum(e["call_type"] == "LOH" for e in ce),
                           "number_gainloh": sum(e["call"] == "GAINLOH" for e in ce),
                           "number_mosaic": sum(bool(e["mosaic"]) for e in ce)})  # fmt: skip
        flags = [f for f in str(s.get("qc_flags", "")).split(";") if f]
        doc = {
            "schema": "simurg.cyto/1", "software_version": __version__,
            "annotation_bundle": {k: bundle[k] for k in ("file", "sha256", "sources")} if bundle else None,
            "reference_genome": side.get("reference"), "genome_build": build,
            "annotation_type": "constitutional", "iscn_version": "ISCN 2020",
            "sample_id": s["sample_id"], "sample_name": s.get("sample_name"),
            "gc_correct": bool(s.get("gc_correct")), "gcwf": s.get("gcwf"),
            "parameters": {**th, "min_qual": opts.min_qual, "size_basis": opts.size_basis,
                           "iscn_min_size": opts.iscn_min_size},
            "overall_ploidy": round(float(s["overall_ploidy"]), 3), "call_rate": s.get("call_rate"),
            "logr_dev": s.get("logr_dev"), "logr_dev_robust": s.get("logr_dev_robust"), "baf_dev": s.get("baf_dev"),
            **loh_stats(loh, eff_genome),
            "copy_number_median": _weighted_median_cn(chroms, events, s["sex_used"]),
            "sex_estimate": s["sex_used"], "sex_source": s["sex_source"], "qc_flags": flags,
            "iscn_traditional": trad, "iscn_microarray": micro,
            "chromosomes": chroms,
            "events": sorted(events, key=lambda e: (CHROM_ORDER.get(e["chrom"], 99), e["start"])),
        }  # fmt: skip
        atomic_write_text(out / f"{s['sample_id']}.cyto.json", json.dumps(doc, indent=1, default=_json_default))
    rep.info(f"annotated {len(samples)} sample(s) -> {out}")
    return ExitCode.OK


def _weighted_median_cn(chroms: list[dict[str, Any]], events: list[dict[str, Any]], sex: str) -> float:
    """CYTO-026: length-weighted median CN over all segments (baseline stretches and CNV events)."""
    parts: list[tuple[float, int]] = []
    for c in chroms:
        base = 1 if (c["id"] == "Y" or (c["id"] == "X" and sex == "M")) else 2
        ce = [e for e in events if e["chrom"] == c["id"] and e["call_type"] != "LOH"]
        cov = sum(e["effective_size"] for e in ce)
        parts.append((float(base), max(c["effective_size"] - cov, 0)))
        parts += [(float(e["copy_number"]), e["effective_size"]) for e in ce]
    if not parts:
        return 2.0
    parts.sort()
    half = sum(w for _, w in parts) / 2
    acc = 0
    for cn, w in parts:
        acc += w
        if acc >= half:
            return cn
    return parts[-1][0]