#!/usr/bin/env python3
"""Registered engineering pipeline for GAL-FINGERPRINT-I protocol v0.1.

The script does not alter STRUC-I. It prepares deterministic pilot ladders, creates
six-slot chamber batches, audits chamber exports, consolidates complete profiles, and
computes the registered source-retrieval distances against the frozen R1 gallery.
"""

from __future__ import annotations

import argparse
import csv
import hashlib
import json
import math
import shutil
import sys
import zipfile
from dataclasses import dataclass
from pathlib import Path
from typing import Iterable, Sequence

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

PROTOCOL_VERSION = "0.1"
INSTRUMENT = "STRUC-I"
INSTRUMENT_VERSION = "1.0.4"
MC_RUNS = 2000
KAPPA_MIN = 0.01
KAPPA_MAX = 1.0
KAPPA_STEPS = 40
KPC_M = 3.0856775814913673e19
ACC_FACTOR = 1.0e6 / KPC_M
BASE_DISK = 0.5
BASE_BULGE = 0.7

PILOT_PANEL = [
    "NGC4010", "UGC00731", "NGC4214", "NGC3769",
    "UGC08286", "NGC4183", "UGC05721", "UGC02885",
    "UGC07524", "NGC5985", "NGC4013", "UGC03580",
]

RUN_COLUMNS = [
    "protocol_version", "branch", "realization_id", "repeat_id", "galaxy_id",
    "source_galaxy_id", "generation_rule", "generation_seed", "n_source",
    "n_realization", "radial_coverage_min", "radial_coverage_max", "upsilon_disk",
    "upsilon_bulge", "input_path", "input_sha256", "json_export_path",
    "json_sha256", "csv_export_path", "csv_sha256", "chamber_timestamp",
    "instrument_version", "mc_runs", "kappa_min", "kappa_max", "kappa_steps",
    "execution_status", "exclusion_reason", "ladder_name", "batch_id",
    "batch_slot", "variant_details",
]

PROFILE_COLUMNS = [
    "protocol_version", "branch", "realization_id", "galaxy_id", "source_galaxy_id",
    "ladder_name", "source_json", "source_csv", "chamber_timestamp", "kappa_index",
    "kappa", "epsilon", "nu_V", "vuln_count", "inv", "rho", "Ak",
]


def project_root() -> Path:
    return Path(__file__).resolve().parents[1]


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


def stable_rng(seed: int, key: str) -> np.random.Generator:
    key_int = int(hashlib.sha256(key.encode("utf-8")).hexdigest()[:16], 16)
    sequence = np.random.SeedSequence([int(seed), key_int & 0xFFFFFFFF, key_int >> 32])
    return np.random.Generator(np.random.PCG64(sequence))


def canonical_kappa_grid() -> np.ndarray:
    return np.geomspace(KAPPA_MIN, KAPPA_MAX, KAPPA_STEPS)


def read_primary_ids(root: Path) -> list[str]:
    path = root / "inputs/sample_lists/primary_121.txt"
    values = [line.strip() for line in path.read_text(encoding="utf-8").splitlines() if line.strip()]
    if len(values) != 121 or len(set(values)) != 121:
        raise ValueError("Primary sample is not the frozen 121 unique galaxy list")
    return values


def load_profiles(root: Path) -> pd.DataFrame:
    frame = pd.read_csv(root / "data/prepared/baryonic_profiles.csv")
    required = {
        "galaxy_id", "radius_kpc", "radius_normalized", "v_gas_kms", "v_disk_kms",
        "v_bulge_kms", "valid_point", "log10_baryonic_acceleration",
    }
    missing = sorted(required.difference(frame.columns))
    if missing:
        raise ValueError(f"Prepared profile table lacks columns: {missing}")
    if frame["valid_point"].dtype != bool:
        frame["valid_point"] = frame["valid_point"].astype(str).str.lower().eq("true")
    return frame


def valid_radial_rows(frame: pd.DataFrame, galaxy_id: str) -> pd.DataFrame:
    rows = frame.loc[(frame["galaxy_id"] == galaxy_id) & frame["valid_point"]].copy()
    rows = rows.sort_values(["radius_kpc", "log10_baryonic_acceleration"], kind="mergesort")
    if len(rows) < 2:
        raise ValueError(f"{galaxy_id}: fewer than two valid radial rows")
    return rows.reset_index(drop=True)


def values_from_rows(rows: pd.DataFrame, disk: float = BASE_DISK, bulge: float = BASE_BULGE) -> np.ndarray:
    g2 = (
        rows["v_gas_kms"].to_numpy(float) * np.abs(rows["v_gas_kms"].to_numpy(float))
        + float(disk) * rows["v_disk_kms"].to_numpy(float) ** 2
        + float(bulge) * rows["v_bulge_kms"].to_numpy(float) ** 2
    )
    radius = rows["radius_kpc"].to_numpy(float)
    acc = (g2 / radius) * ACC_FACTOR
    keep = np.isfinite(acc) & np.isfinite(radius) & (radius > 0) & (g2 > 0) & (acc > 0)
    values = np.unique(np.sort(np.log10(acc[keep])))
    if values.size < 2:
        raise ValueError("Generated ladder has fewer than two finite positive unique values")
    return values


def evenly_spaced_indices(n: int, k: int) -> np.ndarray:
    if n < k:
        raise ValueError(f"Cannot select {k} ranks from n={n}")
    indices = np.rint(np.linspace(0, n - 1, k)).astype(int)
    if len(np.unique(indices)) != k:
        raise ValueError(f"Rounded rank rule produced duplicates for n={n}, k={k}")
    return indices


def deletion_indices(n: int) -> list[int]:
    if n <= 12:
        return list(range(n))
    return sorted({int(round(q * (n - 1))) for q in (0.0, 0.25, 0.50, 0.75, 1.0)})


def representative_deletion_index(n: int, panel_index: int) -> int:
    if n <= 12:
        return int(round(0.5 * (n - 1)))
    return deletion_indices(n)[panel_index % 5]


def random_stratified_indices(n: int, seed: int, galaxy_id: str) -> np.ndarray:
    if n < 10:
        raise ValueError(f"Random n=10 thinning requires n>=10, got n={n}")
    if n == 10:
        return np.arange(10, dtype=int)
    interior = np.arange(1, n - 1, dtype=int)
    strata = [part for part in np.array_split(interior, 8) if len(part)]
    if len(strata) != 8:
        raise ValueError(f"Expected eight non-empty interior strata for n={n}")
    rng = stable_rng(seed, galaxy_id)
    selected = [0] + [int(rng.choice(part)) for part in strata] + [n - 1]
    selected = np.array(sorted(selected), dtype=int)
    if len(np.unique(selected)) != 10:
        raise ValueError(f"Random stratified selection is not ten unique ranks for {galaxy_id}")
    return selected


def gap_order_null(values: np.ndarray, seed: int, galaxy_id: str) -> np.ndarray:
    gaps = np.diff(values)
    if np.any(gaps <= 0):
        raise ValueError(f"{galaxy_id}: canonical ladder must have positive gaps")
    rng = stable_rng(seed, galaxy_id)
    permuted = rng.permutation(gaps)
    result = np.concatenate(([values[0]], values[0] + np.cumsum(permuted)))
    return result


def uniform_range_null(values: np.ndarray, seed: int, galaxy_id: str) -> np.ndarray:
    if len(values) == 2:
        return values.copy()
    rng = stable_rng(seed, galaxy_id)
    interior = np.sort(rng.uniform(values[0], values[-1], len(values) - 2))
    return np.concatenate(([values[0]], interior, [values[-1]]))


def write_ladder(path: Path, values: np.ndarray) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    pd.DataFrame({"value": np.asarray(values, dtype=float)}).to_csv(path, index=False, float_format="%.17g")


def relative(root: Path, path: Path) -> str:
    return path.resolve().relative_to(root.resolve()).as_posix()


@dataclass
class Variant:
    branch: str
    realization_id: str
    galaxy_id: str
    generation_rule: str
    seed: int | None
    values: np.ndarray
    n_source: int
    radial_min: float | None
    radial_max: float | None
    disk: float
    bulge: float
    details: str

    @property
    def ladder_name(self) -> str:
        safe = self.realization_id.replace("-", "_").replace(".", "p")
        return f"GF1_{self.galaxy_id}_{safe}"


def pilot_variants(root: Path) -> list[Variant]:
    frame = load_profiles(root)
    variants: list[Variant] = []
    for panel_index, galaxy_id in enumerate(PILOT_PANEL):
        rows = valid_radial_rows(frame, galaxy_id)
        n = len(rows)
        baseline = values_from_rows(rows)

        delete_index = representative_deletion_index(n, panel_index)
        kept = rows.drop(index=delete_index).reset_index(drop=True)
        variants.append(Variant(
            "L", f"L_DEL_R{delete_index + 1:03d}", galaxy_id,
            "pilot representative of registered single-radius deletion rule", None,
            values_from_rows(kept), n, float(kept.radius_normalized.min()),
            float(kept.radius_normalized.max()), BASE_DISK, BASE_BULGE,
            f"deleted 1-based radial rank {delete_index + 1} of {n}",
        ))

        det_idx = evenly_spaced_indices(n, 10)
        det_rows = rows.iloc[det_idx].reset_index(drop=True)
        variants.append(Variant(
            "T", "T_DET10", galaxy_id,
            "registered deterministic n=10 radial-rank thinning", None,
            values_from_rows(det_rows), n, float(det_rows.radius_normalized.min()),
            float(det_rows.radius_normalized.max()), BASE_DISK, BASE_BULGE,
            "1+round(j(n-1)/9), j=0..9",
        ))

        rand_idx = random_stratified_indices(n, 91001, galaxy_id)
        rand_rows = rows.iloc[rand_idx].reset_index(drop=True)
        variants.append(Variant(
            "T", "T_RAND10_S91001", galaxy_id,
            "registered random stratified n=10 thinning pilot seed", 91001,
            values_from_rows(rand_rows), n, float(rand_rows.radius_normalized.min()),
            float(rand_rows.radius_normalized.max()), BASE_DISK, BASE_BULGE,
            "endpoints plus one draw from each of eight contiguous interior strata",
        ))

        full_idx = evenly_spaced_indices(n, 8)
        full_rows = rows.iloc[full_idx].reset_index(drop=True)
        variants.append(Variant(
            "C", "C_FULL8", galaxy_id,
            "registered full-coverage n=8 radial-rank selection", None,
            values_from_rows(full_rows), n, float(full_rows.radius_normalized.min()),
            float(full_rows.radius_normalized.max()), BASE_DISK, BASE_BULGE,
            "eight evenly spaced ranks over full observed range",
        ))

        inner_rows_all = rows.loc[rows["radius_normalized"] <= 0.85].reset_index(drop=True)
        if len(inner_rows_all) < 8:
            raise ValueError(f"{galaxy_id}: fewer than eight rows at R/Rmax<=0.85")
        inner_idx = evenly_spaced_indices(len(inner_rows_all), 8)
        inner_rows = inner_rows_all.iloc[inner_idx].reset_index(drop=True)
        variants.append(Variant(
            "C", "C_INNER85_8", galaxy_id,
            "registered inner-coverage n=8 radial-rank selection", None,
            values_from_rows(inner_rows), n, float(inner_rows.radius_normalized.min()),
            float(inner_rows.radius_normalized.max()), BASE_DISK, BASE_BULGE,
            "eight evenly spaced ranks among R/Rmax<=0.85 rows",
        ))

        if panel_index % 2 == 0:
            disk, rid = 0.4, "M_D_LOW"
        else:
            disk, rid = 0.6, "M_D_HIGH"
        variants.append(Variant(
            "M", rid, galaxy_id,
            "pilot alternation across registered disk M/L sensitivity levels", None,
            values_from_rows(rows, disk=disk, bulge=BASE_BULGE), n,
            float(rows.radius_normalized.min()), float(rows.radius_normalized.max()),
            disk, BASE_BULGE, "signed gas preserved; baseline-valid radial rows retained",
        ))

        has_bulge = bool(np.any(np.abs(rows["v_bulge_kms"].to_numpy(float)) > 0))
        if has_bulge:
            bulge, rid = ((0.6, "M_B_LOW") if panel_index % 2 == 0 else (0.8, "M_B_HIGH"))
            variants.append(Variant(
                "M", rid, galaxy_id,
                "pilot alternation across registered bulge M/L sensitivity levels", None,
                values_from_rows(rows, disk=BASE_DISK, bulge=bulge), n,
                float(rows.radius_normalized.min()), float(rows.radius_normalized.max()),
                BASE_DISK, bulge, "signed gas preserved; baseline-valid radial rows retained",
            ))

        variants.append(Variant(
            "N_GAP", "N_GAP_S92001", galaxy_id,
            "registered gap-order null pilot seed", 92001,
            gap_order_null(baseline, 92001, galaxy_id), n,
            float(rows.radius_normalized.min()), float(rows.radius_normalized.max()),
            BASE_DISK, BASE_BULGE, "preserves n, endpoints, total range, and complete gap multiset",
        ))
        variants.append(Variant(
            "N_UNIFORM", "N_UNIFORM_S93001", galaxy_id,
            "registered uniform-range null pilot seed", 93001,
            uniform_range_null(baseline, 93001, galaxy_id), n,
            float(rows.radius_normalized.min()), float(rows.radius_normalized.max()),
            BASE_DISK, BASE_BULGE, "preserves n, minimum, and maximum",
        ))

    return variants


def clear_generated_pilot(root: Path) -> tuple[Path, Path]:
    input_root = root / "inputs/ladders/struc_i/gal_fingerprint_i/pilot"
    manifest_root = root / "runs/gal_fingerprint_i/manifests"
    if input_root.exists():
        shutil.rmtree(input_root)
    input_root.mkdir(parents=True, exist_ok=True)
    manifest_root.mkdir(parents=True, exist_ok=True)
    return input_root, manifest_root


def prepare_pilot(root: Path) -> None:
    input_root, manifest_root = clear_generated_pilot(root)
    variants = pilot_variants(root)
    if len(variants) != 100:
        raise ValueError(f"Frozen pilot must contain 100 realizations, generated {len(variants)}")

    ordered = sorted(variants, key=lambda v: (v.branch, v.galaxy_id, v.realization_id))
    rows: list[dict[str, object]] = []
    batch_root = input_root / "chamber_batches"
    batch_root.mkdir(parents=True, exist_ok=True)

    for index, variant in enumerate(ordered):
        batch_number = index // 6 + 1
        slot = index % 6 + 1
        batch_id = f"PILOT_BATCH_{batch_number:03d}"
        batch_dir = batch_root / batch_id
        batch_dir.mkdir(parents=True, exist_ok=True)
        path = batch_dir / f"{variant.ladder_name}.csv"
        write_ladder(path, variant.values)
        rows.append({
            "protocol_version": PROTOCOL_VERSION,
            "branch": variant.branch,
            "realization_id": variant.realization_id,
            "repeat_id": "PILOT",
            "galaxy_id": variant.galaxy_id,
            "source_galaxy_id": variant.galaxy_id,
            "generation_rule": variant.generation_rule,
            "generation_seed": "" if variant.seed is None else variant.seed,
            "n_source": variant.n_source,
            "n_realization": len(variant.values),
            "radial_coverage_min": variant.radial_min,
            "radial_coverage_max": variant.radial_max,
            "upsilon_disk": variant.disk,
            "upsilon_bulge": variant.bulge,
            "input_path": relative(root, path),
            "input_sha256": sha256_path(path),
            "json_export_path": "",
            "json_sha256": "",
            "csv_export_path": "",
            "csv_sha256": "",
            "chamber_timestamp": "",
            "instrument_version": INSTRUMENT_VERSION,
            "mc_runs": MC_RUNS,
            "kappa_min": KAPPA_MIN,
            "kappa_max": KAPPA_MAX,
            "kappa_steps": KAPPA_STEPS,
            "execution_status": "INPUT_READY",
            "exclusion_reason": "",
            "ladder_name": variant.ladder_name,
            "batch_id": batch_id,
            "batch_slot": slot,
            "variant_details": variant.details,
        })

    manifest = pd.DataFrame(rows, columns=RUN_COLUMNS)
    manifest.to_csv(manifest_root / "run_manifest.csv", index=False)
    manifest.to_csv(manifest_root / "input_ladder_manifest.csv", index=False)

    batch_summary_rows = []
    for batch_dir in sorted(p for p in batch_root.iterdir() if p.is_dir()):
        files = sorted(batch_dir.glob("*.csv"))
        batch_zip = batch_root / f"{batch_dir.name}.zip"
        with zipfile.ZipFile(batch_zip, "w", compression=zipfile.ZIP_DEFLATED) as archive:
            for file in files:
                archive.write(file, arcname=file.name)
        batch_summary_rows.append({
            "batch_id": batch_dir.name,
            "slot_count": len(files),
            "batch_zip": relative(root, batch_zip),
            "batch_zip_sha256": sha256_path(batch_zip),
            "ladder_names": "|".join(file.stem for file in files),
        })
    pd.DataFrame(batch_summary_rows).to_csv(input_root / "pilot_batch_manifest.csv", index=False)

    master_zip = input_root / "GAL_FINGERPRINT_I_PILOT_BATCHES.zip"
    with zipfile.ZipFile(master_zip, "w", compression=zipfile.ZIP_DEFLATED) as archive:
        for batch_dir in sorted(p for p in batch_root.iterdir() if p.is_dir()):
            for file in sorted(batch_dir.glob("*.csv")):
                archive.write(file, arcname=f"{batch_dir.name}/{file.name}")
        archive.write(input_root / "pilot_batch_manifest.csv", arcname="pilot_batch_manifest.csv")
        archive.write(manifest_root / "run_manifest.csv", arcname="run_manifest.csv")

    print(f"PILOT_INPUTS_READY realizations={len(manifest)} batches={len(batch_summary_rows)}")
    print(f"MASTER_ZIP {relative(root, master_zip)} SHA256={sha256_path(master_zip)}")


def unique_files(paths: Iterable[Path]) -> tuple[list[Path], list[dict[str, str]]]:
    by_hash: dict[str, Path] = {}
    audit: list[dict[str, str]] = []
    for path in sorted(paths):
        digest = sha256_path(path)
        duplicate_of = ""
        if digest in by_hash:
            duplicate_of = by_hash[digest].name
        else:
            by_hash[digest] = path
        audit.append({"path": str(path), "sha256": digest, "duplicate_of": duplicate_of})
    return list(by_hash.values()), audit


def read_json_export(path: Path) -> tuple[dict, pd.DataFrame]:
    payload = json.loads(path.read_text(encoding="utf-8"))
    for key in ("chamber", "version", "timestamp", "kappa_range", "kappa_steps", "mc_runs", "profiles"):
        if key not in payload:
            raise ValueError(f"{path.name}: missing JSON field {key}")
    if payload["chamber"] != INSTRUMENT or str(payload["version"]) != INSTRUMENT_VERSION:
        raise ValueError(f"{path.name}: wrong instrument/version")
    if int(payload["kappa_steps"]) != KAPPA_STEPS or int(payload["mc_runs"]) != MC_RUNS:
        raise ValueError(f"{path.name}: wrong κ-step or Monte Carlo setting")
    krange = np.asarray(payload["kappa_range"], dtype=float)
    if not np.allclose(krange, [KAPPA_MIN, KAPPA_MAX], rtol=0, atol=1e-12):
        raise ValueError(f"{path.name}: wrong κ range {krange.tolist()}")
    rows: list[dict[str, object]] = []
    name_counts: dict[str, int] = {}
    for profile in payload["profiles"]:
        name = str(profile["name"])
        profile_instance = name_counts.get(name, 0)
        name_counts[name] = profile_instance + 1
        curve = profile["curve"]
        if len(curve) != KAPPA_STEPS:
            raise ValueError(f"{path.name}/{name}: expected {KAPPA_STEPS} curve rows")
        for index, point in enumerate(curve):
            rows.append({
                "ladder_name": name,
                "profile_instance": profile_instance,
                "kappa_index": index,
                "kappa": float(point["kappa"]),
                "epsilon": float(point["epsilon"]) if "epsilon" in point else np.nan,
                "nu_V": float(point["nu"]),
                "vuln_count": float(point["vulnCount"]),
                "inv": float(point["inv"]),
                "rho": float(point["rho"]),
                "Ak": float(point["Ak"]),
            })
    return payload, pd.DataFrame(rows)


def read_csv_export(path: Path) -> pd.DataFrame:
    frame = pd.read_csv(path)
    required = {"ladder", "kappa", "epsilon", "nu_V", "vuln_count", "inv", "rho", "Ak"}
    missing = sorted(required.difference(frame.columns))
    if missing:
        raise ValueError(f"{path.name}: missing CSV columns {missing}")
    frame = frame.rename(columns={"ladder": "ladder_name"})
    frame["profile_instance"] = frame.groupby("ladder_name", sort=False).cumcount() // KAPPA_STEPS
    frame["kappa_index"] = frame.groupby(["ladder_name", "profile_instance"], sort=False).cumcount()
    return frame


def validate_kappa(frame: pd.DataFrame, context: str, require_epsilon: bool = True) -> None:
    expected = canonical_kappa_grid()
    group_keys = ["ladder_name"] + (["profile_instance"] if "profile_instance" in frame.columns else [])
    for key, group in frame.groupby(group_keys, sort=False):
        name = key[0] if isinstance(key, tuple) else key
        group = group.sort_values("kappa")
        values = group["kappa"].to_numpy(float)
        if len(values) != KAPPA_STEPS:
            raise ValueError(f"{context}/{name}: expected 40 κ points")
        if not np.allclose(values, expected, rtol=1e-6, atol=5e-7):
            raise ValueError(f"{context}/{name}: κ grid mismatch")
        columns = ["kappa", "nu_V", "vuln_count", "inv", "rho", "Ak"]
        if require_epsilon:
            columns.insert(1, "epsilon")
        numeric = group[columns].to_numpy(float)
        if not np.isfinite(numeric).all():
            raise ValueError(f"{context}/{name}: non-finite curve coordinate")


def pair_exports(json_paths: Sequence[Path], csv_paths: Sequence[Path]) -> list[tuple[Path, Path, dict, pd.DataFrame]]:
    csv_by_names: dict[tuple[str, ...], list[tuple[Path, pd.DataFrame]]] = {}
    for path in csv_paths:
        frame = read_csv_export(path)
        names = tuple(sorted(frame["ladder_name"].astype(str).unique()))
        csv_by_names.setdefault(names, []).append((path, frame))
    pairs = []
    used_csv: set[Path] = set()
    for json_path in json_paths:
        payload, jframe = read_json_export(json_path)
        names = tuple(sorted(jframe["ladder_name"].astype(str).unique()))
        candidates = [item for item in csv_by_names.get(names, []) if item[0] not in used_csv]
        if not candidates:
            raise ValueError(f"{json_path.name}: no CSV export with the same ladder-name set")
        csv_path, cframe = sorted(candidates, key=lambda item: item[0].name)[0]
        used_csv.add(csv_path)
        validate_kappa(jframe, json_path.name, require_epsilon=False)
        validate_kappa(cframe, csv_path.name, require_epsilon=True)
        # JSON retains the full κ precision; compare other exported columns after ordered alignment.
        j = jframe.sort_values(["ladder_name", "profile_instance", "kappa_index"]).reset_index(drop=True)
        c = cframe.sort_values(["ladder_name", "profile_instance", "kappa_index"]).reset_index(drop=True)
        if len(j) != len(c):
            raise ValueError(f"{json_path.name}/{csv_path.name}: row-count mismatch")
        if not np.array_equal(j["ladder_name"].to_numpy(str), c["ladder_name"].to_numpy(str)):
            raise ValueError(f"{json_path.name}/{csv_path.name}: ladder ordering mismatch")
        for column in ("nu_V", "vuln_count", "inv", "rho", "Ak"):
            if not np.allclose(j[column].to_numpy(float), c[column].to_numpy(float), rtol=1e-5, atol=5e-6):
                raise ValueError(f"{json_path.name}/{csv_path.name}: {column} mismatch")
        # The chamber JSON omits epsilon; preserve the paired CSV epsilon while retaining JSON κ precision.
        j["epsilon"] = c["epsilon"].to_numpy(float)
        pairs.append((json_path, csv_path, payload, j))
    if len(used_csv) != len(csv_paths):
        unused = sorted(path.name for path in csv_paths if path not in used_csv)
        raise ValueError(f"Unpaired unique CSV exports: {unused}")
    return pairs


def reference_gallery(root: Path) -> tuple[list[str], np.ndarray, dict[str, pd.DataFrame]]:
    profiles = pd.read_csv(root / "runs/consolidated/struc_i_primary_121_profiles.csv")
    validate_kappa(profiles.rename(columns={"galaxy_id": "ladder_name"}), "R1 gallery", require_epsilon=False)
    ids = sorted(read_primary_ids(root))
    grid = canonical_kappa_grid()
    galleries: dict[str, pd.DataFrame] = {}
    for channel in ("rho", "Ak", "nu_V"):
        matrix = profiles.pivot(index="galaxy_id", columns="kappa", values=channel)
        # κ values in CSV are rounded; use sorted positional alignment after grid validation.
        matrix = matrix.loc[ids, sorted(matrix.columns)].copy()
        matrix.columns = range(KAPPA_STEPS)
        galleries[channel] = matrix
    return ids, grid, galleries


def retrieval_results(root: Path, profiles: pd.DataFrame, manifest: pd.DataFrame) -> pd.DataFrame:
    ids, _, gallery = reference_gallery(root)
    manifest_map = manifest.set_index("ladder_name")
    reference_rho = gallery["rho"].to_numpy(float)
    scaler_rho = StandardScaler().fit(reference_rho)
    ref_rho_z = scaler_rho.transform(reference_rho)

    channel_ref = []
    channel_scalers = []
    for channel in ("rho", "Ak", "nu_V"):
        values = gallery[channel].to_numpy(float)
        if channel == "Ak":
            values = 1.0 - values
        if channel == "nu_V":
            n_map = pd.read_csv(root / "runs/consolidated/struc_i_primary_121_summary.csv").set_index("galaxy_id")["n"]
            denom = np.ceil((n_map.loc[ids].to_numpy(float) - 1.0) / 2.0)[:, None]
            values = values / np.maximum(denom, 1.0)
        scaler = StandardScaler().fit(values)
        channel_scalers.append(scaler)
        channel_ref.append(scaler.transform(values))

    pca = PCA(n_components=3, svd_solver="full").fit(ref_rho_z)
    ref_pc = pca.transform(ref_rho_z)

    out = []
    for ladder_name, group in profiles.groupby("ladder_name", sort=True):
        if ladder_name not in manifest_map.index:
            raise ValueError(f"Exported ladder {ladder_name} is absent from pilot manifest")
        row = manifest_map.loc[ladder_name]
        group = group.sort_values("kappa_index")
        rho = group["rho"].to_numpy(float)[None, :]
        rho_z = scaler_rho.transform(rho)
        d_primary = np.sqrt(np.mean((ref_rho_z - rho_z) ** 2, axis=1))

        query_channels = []
        for idx, channel in enumerate(("rho", "Ak", "nu_V")):
            values = group[channel].to_numpy(float)[None, :]
            if channel == "Ak":
                values = 1.0 - values
            if channel == "nu_V":
                denom = math.ceil((float(row["n_realization"]) - 1.0) / 2.0)
                values = values / max(denom, 1)
            query_channels.append(channel_scalers[idx].transform(values))
        per_channel_mse = [np.mean((channel_ref[i] - query_channels[i]) ** 2, axis=1) for i in range(3)]
        d_combined = np.sqrt(np.mean(np.vstack(per_channel_mse), axis=0))
        query_pc = pca.transform(rho_z)
        d_pca = np.sqrt(np.sum((ref_pc - query_pc) ** 2, axis=1))

        for metric, distances in (("rho_rms", d_primary), ("combined_rms", d_combined), ("rho_pca3", d_pca)):
            order = np.argsort(distances, kind="mergesort")
            source = str(row["source_galaxy_id"])
            source_index = ids.index(source)
            source_distance = float(distances[source_index])
            wrong = [i for i in order if ids[i] != source]
            nearest_wrong_index = wrong[0]
            rank = int(np.flatnonzero(order == source_index)[0]) + 1
            minimum = float(distances[order[0]])
            strict_top1 = bool(order[0] == source_index and np.sum(np.isclose(distances, minimum, rtol=0, atol=1e-12)) == 1)
            out.append({
                "ladder_name": ladder_name,
                "branch": row["branch"],
                "realization_id": row["realization_id"],
                "source_galaxy_id": source,
                "metric": metric,
                "own_source_distance": source_distance,
                "nearest_wrong_source": ids[nearest_wrong_index],
                "nearest_wrong_distance": float(distances[nearest_wrong_index]),
                "own_to_nearest_wrong_ratio": source_distance / max(float(distances[nearest_wrong_index]), 1e-12),
                "distance_margin": float(distances[nearest_wrong_index]) - source_distance,
                "source_rank": rank,
                "strict_top1": strict_top1,
                "top5": bool(rank <= 5),
            })
    return pd.DataFrame(out)


def consolidate_pilot(root: Path, export_dir: Path | None = None) -> None:
    manifest_path = root / "runs/gal_fingerprint_i/manifests/run_manifest.csv"
    if not manifest_path.exists():
        raise FileNotFoundError("Pilot run manifest is missing; run prepare-pilot first")
    manifest = pd.read_csv(manifest_path, keep_default_na=False)
    export_dir = export_dir or (root / "outputs/exports/struc_i/gal_fingerprint_i/pilot")
    if not export_dir.exists():
        raise FileNotFoundError(f"Pilot export directory does not exist: {export_dir}")

    json_all = list(export_dir.rglob("*.json"))
    csv_all = [p for p in export_dir.rglob("*.csv") if "manifest" not in p.name.lower()]
    json_unique, json_audit = unique_files(json_all)
    csv_unique, csv_audit = unique_files(csv_all)
    if not json_unique or not csv_unique:
        raise ValueError("Both JSON and CSV chamber exports are required")
    pairs = pair_exports(json_unique, csv_unique)

    profile_rows: list[dict[str, object]] = []
    updated = manifest.copy().set_index("ladder_name", drop=False)
    seen_names: set[str] = set()
    file_rows: list[dict[str, object]] = []
    for json_path, csv_path, payload, frame in pairs:
        names = sorted(frame["ladder_name"].astype(str).unique())
        for name in names:
            if name in seen_names:
                raise ValueError(f"Pilot ladder {name} appears in more than one unique JSON export")
            seen_names.add(name)
            if name not in updated.index:
                raise ValueError(f"Unknown ladder in export: {name}")
            updated.loc[name, "json_export_path"] = relative(root, json_path)
            updated.loc[name, "json_sha256"] = sha256_path(json_path)
            updated.loc[name, "csv_export_path"] = relative(root, csv_path)
            updated.loc[name, "csv_sha256"] = sha256_path(csv_path)
            updated.loc[name, "chamber_timestamp"] = payload["timestamp"]
            updated.loc[name, "execution_status"] = "VALID_COMPLETE"
            sub = frame.loc[frame["ladder_name"] == name].sort_values("kappa_index")
            mrow = updated.loc[name]
            for _, point in sub.iterrows():
                profile_rows.append({
                    "protocol_version": PROTOCOL_VERSION,
                    "branch": mrow["branch"],
                    "realization_id": mrow["realization_id"],
                    "galaxy_id": mrow["galaxy_id"],
                    "source_galaxy_id": mrow["source_galaxy_id"],
                    "ladder_name": name,
                    "source_json": relative(root, json_path),
                    "source_csv": relative(root, csv_path),
                    "chamber_timestamp": payload["timestamp"],
                    "kappa_index": int(point["kappa_index"]),
                    "kappa": float(point["kappa"]),
                    "epsilon": float(point["epsilon"]) if "epsilon" in point else np.nan,
                    "nu_V": float(point["nu_V"]),
                    "vuln_count": float(point["vuln_count"]),
                    "inv": float(point["inv"]),
                    "rho": float(point["rho"]),
                    "Ak": float(point["Ak"]),
                })
        file_rows.append({
            "json_path": relative(root, json_path), "json_sha256": sha256_path(json_path),
            "csv_path": relative(root, csv_path), "csv_sha256": sha256_path(csv_path),
            "chamber_timestamp": payload["timestamp"], "ladder_count": len(names),
            "ladder_names": "|".join(names),
        })

    missing = sorted(set(updated.index) - seen_names)
    for name in missing:
        updated.loc[name, "execution_status"] = "MISSING_EXPORT"
    output = root / "runs/gal_fingerprint_i/consolidated"
    output.mkdir(parents=True, exist_ok=True)
    updated.reset_index(drop=True)[RUN_COLUMNS].to_csv(manifest_path, index=False)
    pd.DataFrame(file_rows).to_csv(root / "runs/gal_fingerprint_i/manifests/export_file_manifest.csv", index=False)
    pd.DataFrame(json_audit + csv_audit).to_csv(root / "runs/gal_fingerprint_i/manifests/export_audit.csv", index=False)
    profiles = pd.DataFrame(profile_rows, columns=PROFILE_COLUMNS)
    profiles.to_csv(output / "profiles_long.csv", index=False)
    if missing:
        print(f"PILOT_PARTIAL complete={len(seen_names)} missing={len(missing)}")
        return
    results = retrieval_results(root, profiles, updated.reset_index(drop=True))
    results.to_csv(output / "source_retrieval_results.csv", index=False)
    summary = (
        results.groupby(["branch", "metric"], as_index=False)
        .agg(
            realizations=("ladder_name", "nunique"),
            strict_top1_accuracy=("strict_top1", "mean"),
            top5_accuracy=("top5", "mean"),
            median_own_to_nearest_wrong_ratio=("own_to_nearest_wrong_ratio", "median"),
            median_distance_margin=("distance_margin", "median"),
        )
    )
    summary.to_csv(output / "branch_summary.csv", index=False)
    print(f"PILOT_CONSOLIDATED profiles={len(profiles)} realizations={profiles['ladder_name'].nunique()}")


def self_test(root: Path) -> None:
    primary = read_primary_ids(root)
    if sha256_path(root / "inputs/sample_lists/primary_121.txt") != "a2cff5356218021afe284f90926623fbcd6329e7b830b3cdea6c4bc41ec56b64":
        raise AssertionError("Frozen primary sample hash changed")
    frame = load_profiles(root)
    count = 0
    for galaxy_id in primary:
        count += len(deletion_indices(len(valid_radial_rows(frame, galaxy_id))))
    if count != 782:
        raise AssertionError(f"Registered deletion realization count changed: {count}")
    variants = pilot_variants(root)
    if len(variants) != 100:
        raise AssertionError("Pilot variant count must be 100")
    if len({v.ladder_name for v in variants}) != 100:
        raise AssertionError("Pilot ladder names are not unique")

    for variant in variants:
        values = variant.values
        if len(values) < 2 or np.any(np.diff(values) <= 0) or not np.isfinite(values).all():
            raise AssertionError(f"Invalid generated ladder {variant.ladder_name}")
        if variant.realization_id.startswith(("T_DET10", "T_RAND10")) and len(values) != 10:
            raise AssertionError(f"Thinning ladder does not contain 10 unique values: {variant.ladder_name}")
        if variant.realization_id.startswith("C_") and len(values) != 8:
            raise AssertionError(f"Coverage ladder does not contain 8 unique values: {variant.ladder_name}")

    # Validate the R1 export corpus and compare its primary curves with the frozen consolidated table.
    repeat = root / "outputs/exports/struc_i/repeat_01"
    json_unique, _ = unique_files(repeat.glob("*results*.json"))
    csv_unique, _ = unique_files(repeat.glob("*profiles*.csv"))
    pairs = pair_exports(json_unique, csv_unique)
    curve_by_name: dict[str, pd.DataFrame] = {}
    for _, _, _, profile_frame in pairs:
        for (name, _instance), group in profile_frame.groupby(["ladder_name", "profile_instance"]):
            group = group.sort_values("kappa_index").reset_index(drop=True)
            if name in curve_by_name:
                previous = curve_by_name[name]
                columns = ["kappa", "nu_V", "vuln_count", "inv", "rho", "Ak"]
                if not np.allclose(previous[columns], group[columns], rtol=1e-9, atol=1e-12):
                    if name in primary:
                        raise AssertionError(f"Conflicting duplicate primary R1 profile: {name}")
            else:
                curve_by_name[name] = group
    if set(primary) - set(curve_by_name):
        raise AssertionError(f"R1 misses primary galaxies: {sorted(set(primary)-set(curve_by_name))}")
    frozen = pd.read_csv(root / "runs/consolidated/struc_i_primary_121_profiles.csv")
    for galaxy_id in primary:
        a = curve_by_name[galaxy_id].sort_values("kappa_index")
        b = frozen.loc[frozen["galaxy_id"] == galaxy_id].sort_values("kappa")
        for column in ("nu_V", "vuln_count", "inv", "rho", "Ak"):
            if not np.allclose(a[column].to_numpy(float), b[column].to_numpy(float), rtol=1e-9, atol=1e-12):
                raise AssertionError(f"R1 mismatch for {galaxy_id}/{column}")
    print("GAL_FINGERPRINT_I_SELF_TEST_PASS primary=121 deletion=782 pilot=100 r1_profiles=4840")


def build_parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--root", type=Path, default=project_root(), help="project root")
    sub = parser.add_subparsers(dest="command", required=True)
    sub.add_parser("prepare-pilot", help="generate the frozen 100-realization engineering pilot")
    sub.add_parser("self-test", help="test registered generation rules and the preserved R1 corpus")
    consolidate = sub.add_parser("consolidate-pilot", help="audit and consolidate real pilot chamber exports")
    consolidate.add_argument("--export-dir", type=Path, default=None)
    return parser


def main(argv: Sequence[str] | None = None) -> int:
    args = build_parser().parse_args(argv)
    root = args.root.resolve()
    if args.command == "prepare-pilot":
        prepare_pilot(root)
    elif args.command == "self-test":
        self_test(root)
    elif args.command == "consolidate-pilot":
        consolidate_pilot(root, args.export_dir)
    else:
        raise AssertionError(args.command)
    return 0


if __name__ == "__main__":
    try:
        raise SystemExit(main())
    except Exception as exc:  # concise CLI failure with non-zero exit
        print(f"GAL_FINGERPRINT_I_ERROR {type(exc).__name__}: {exc}", file=sys.stderr)
        raise
