#!/usr/bin/env python3
"""Rebuild every structural curve-atlas table, figure, and report from frozen profiles."""

from __future__ import annotations

import argparse
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.metrics import pairwise_distances
from sklearn.preprocessing import StandardScaler

CHANNELS = ("rho", "Ak", "nu_V")


def first_kappa(kappa: np.ndarray, mask: np.ndarray) -> float:
    indices = np.flatnonzero(mask)
    return float(kappa[indices[0]]) if len(indices) else float("nan")


def build_descriptors(profiles: pd.DataFrame, summary: pd.DataFrame) -> pd.DataFrame:
    n_map = summary.set_index("galaxy_id")["n"].astype(int)
    rows = []
    for galaxy_id, group in profiles.groupby("galaxy_id", sort=True):
        group = group.sort_values("kappa")
        kappa = group["kappa"].to_numpy(float)
        logk = np.log10(kappa)
        rho = group["rho"].to_numpy(float)
        ak = group["Ak"].to_numpy(float)
        nu = group["nu_V"].to_numpy(float)
        vuln = group["vuln_count"].to_numpy(float)
        n = int(n_map.loc[galaxy_id])
        peak_index = int(np.argmax(rho))
        rows.append({
            "galaxy_id": galaxy_id,
            "ladder_size": n,
            "rho_onset_kappa": first_kappa(kappa, nu > 0),
            "rho_peak": float(rho[peak_index]),
            "rho_peak_kappa": float(kappa[peak_index]),
            "rho_final": float(rho[-1]),
            "rho_auc_logk": float(np.trapezoid(rho, logk)),
            "rho_early_mean": float(rho[kappa < 0.1].mean()),
            "rho_mid_mean": float(rho[(kappa >= 0.1) & (kappa < 0.5)].mean()),
            "rho_late_mean": float(rho[kappa >= 0.5].mean()),
            "rho_roughness": float(np.mean(np.abs(np.diff(rho)))),
            "Ak_min": float(np.min(ak)),
            "Ak_deficit_auc_logk": float(np.trapezoid(1.0 - ak, logk)),
            "Ak_below_one_steps": int(np.sum(ak < 1.0)),
            "nu_onset_kappa": first_kappa(kappa, nu > 0),
            "nu_final": float(nu[-1]),
            "nu_final_fraction": float(nu[-1] / n),
            "vulnerability_final": float(vuln[-1]),
            "vulnerability_final_fraction": float(vuln[-1] / n),
        })
    return pd.DataFrame(rows)


def build_pca(profiles: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
    galaxy_ids = sorted(profiles["galaxy_id"].unique())
    kappa_grid = sorted(profiles["kappa"].unique())
    score_frames = []
    variance_rows = []
    loading_rows = []

    for channel in CHANNELS:
        matrix = profiles.pivot(index="galaxy_id", columns="kappa", values=channel)
        matrix = matrix.loc[galaxy_ids, kappa_grid].astype(float)
        values = matrix.to_numpy()
        if channel == "Ak":
            values = 1.0 - values
        standardized = StandardScaler().fit_transform(values)
        pca = PCA(n_components=3, svd_solver="full")
        scores = pca.fit_transform(standardized)
        score_frame = pd.DataFrame(
            scores,
            index=galaxy_ids,
            columns=[f"{channel}_PC1", f"{channel}_PC2", f"{channel}_PC3"],
        )
        score_frames.append(score_frame)
        cumulative = np.cumsum(pca.explained_variance_ratio_)
        for component, (ratio, cum) in enumerate(zip(pca.explained_variance_ratio_, cumulative), start=1):
            variance_rows.append({
                "channel": channel,
                "component": component,
                "explained_variance_ratio": float(ratio),
                "cumulative_explained_variance_ratio": float(cum),
            })
            for kappa, loading in zip(kappa_grid, pca.components_[component - 1]):
                loading_rows.append({
                    "channel": channel,
                    "component": component,
                    "kappa": float(kappa),
                    "loading": float(loading),
                })

    scores = pd.concat(score_frames, axis=1)
    combined_input = StandardScaler().fit_transform(scores.to_numpy(float))
    combined_pca = PCA(n_components=3, svd_solver="full")
    combined_scores = combined_pca.fit_transform(combined_input)
    for component in range(3):
        scores[f"structural_PC{component + 1}"] = combined_scores[:, component]
    cumulative = np.cumsum(combined_pca.explained_variance_ratio_)
    for component, (ratio, cum) in enumerate(zip(combined_pca.explained_variance_ratio_, cumulative), start=1):
        variance_rows.append({
            "channel": "combined",
            "component": component,
            "explained_variance_ratio": float(ratio),
            "cumulative_explained_variance_ratio": float(cum),
        })

    scores.index.name = "galaxy_id"
    return scores.reset_index(), pd.DataFrame(variance_rows), pd.DataFrame(loading_rows)


def build_neighbors(scores: pd.DataFrame) -> pd.DataFrame:
    ordered = scores.sort_values("galaxy_id").reset_index(drop=True)
    values = ordered[["structural_PC1", "structural_PC2", "structural_PC3"]].to_numpy(float)
    distances = pairwise_distances(values)
    np.fill_diagonal(distances, np.inf)
    nearest = np.argmin(distances, axis=1)
    return pd.DataFrame({
        "galaxy_id": ordered["galaxy_id"],
        "nearest_structural_neighbor": ordered.iloc[nearest]["galaxy_id"].to_numpy(),
        "distance": distances[np.arange(len(ordered)), nearest],
    })


def write_figures(profiles: pd.DataFrame, scores: pd.DataFrame, variance: pd.DataFrame, figures: Path) -> None:
    figures.mkdir(parents=True, exist_ok=True)

    fig, ax = plt.subplots(figsize=(10, 6))
    for _, group in profiles.groupby("galaxy_id", sort=True):
        group = group.sort_values("kappa")
        ax.plot(group["kappa"], group["rho"], linewidth=0.8, alpha=0.45)
    ax.set_xscale("log")
    ax.set_xlabel("κ")
    ax.set_ylabel("ρ(κ)")
    ax.set_title("STRUC-I ρ curve atlas — primary 121 galaxies")
    fig.tight_layout()
    fig.savefig(figures / "rho_curve_atlas.png", dpi=180)
    plt.close(fig)

    fig, ax = plt.subplots(figsize=(8, 6))
    ax.scatter(scores["structural_PC1"], scores["structural_PC2"], s=25)
    ax.set_xlabel("Structural PC1")
    ax.set_ylabel("Structural PC2")
    ax.set_title("Galaxy structural curve coordinates")
    fig.tight_layout()
    fig.savefig(figures / "structural_pca_map.png", dpi=180)
    plt.close(fig)

    rho = variance.loc[variance["channel"] == "rho"].sort_values("component")
    fig, ax = plt.subplots(figsize=(7, 5))
    ax.bar(rho["component"].astype(str), rho["explained_variance_ratio"] * 100.0)
    ax.plot(rho["component"].astype(str), rho["cumulative_explained_variance_ratio"] * 100.0, marker="o")
    ax.set_xlabel("ρ PCA component")
    ax.set_ylabel("Explained variance (%)")
    ax.set_title("Low-dimensionality of the ρ curves")
    fig.tight_layout()
    fig.savefig(figures / "rho_pca_variance.png", dpi=180)
    plt.close(fig)


def write_report(profiles: pd.DataFrame, variance: pd.DataFrame, path: Path) -> None:
    rho = variance.loc[variance["channel"] == "rho"].sort_values("component")
    ratios = rho["explained_variance_ratio"].to_numpy(float)
    cumulative = rho["cumulative_explained_variance_ratio"].to_numpy(float)
    report = f"""# Structural curve atlas — foundation report

## Corpus

- Primary galaxies: **{profiles['galaxy_id'].nunique()}**
- STRUC-I profile rows: **{len(profiles)}**
- κ points per galaxy: **{int(profiles.groupby('galaxy_id').size().iloc[0])}**
- Curve channels retained: **ρ(κ), Aκ(κ), ν(Vκ)(κ)**

## What is now observable

Each galaxy is represented by a multiscale structural-response fingerprint. The atlas records:

- activation onset;
- ρ peak height and location;
- integrated ρ response across the logarithmic κ range;
- early-, middle-, and late-scale response;
- Aκ deficit burden;
- normalized vulnerability capacity;
- low-dimensional curve coordinates;
- nearest structural neighbours.

## Low-dimensionality

The first three ρ-curve components explain:

- PC1: **{ratios[0] * 100:.3f}%**
- PC1–PC2 cumulative: **{cumulative[1] * 100:.3f}%**
- PC1–PC3 cumulative: **{cumulative[2] * 100:.3f}%**

This supports treating the complete curves as a compact structural coordinate system, while
not yet treating any coordinate as an intrinsic physical property.

## Immediate research gate

The registered next experiment is **GAL-FINGERPRINT-I — Galaxy Structural Identity and
Stability Audit**. It tests native repeat stability, single-radius deletion, matched point-count
thinning, fixed-count radial-coverage truncation, declared baryonic-preparation sensitivity,
and registered null models.

The structural-fingerprint interpretation is supported only when within-galaxy variation is
substantially smaller than between-galaxy variation and source identity survives the registered
sampling controls.
"""
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(report, encoding="utf-8")


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output-root", type=Path, default=None, help="alternate atlas outputs directory")
    args = parser.parse_args()

    root = Path(__file__).resolve().parents[3]
    profiles = pd.read_csv(root / "runs/consolidated/struc_i_primary_121_profiles.csv")
    summary = pd.read_csv(root / "runs/consolidated/struc_i_primary_121_summary.csv")
    if len(profiles) != 4840 or profiles["galaxy_id"].nunique() != 121:
        raise ValueError("Frozen primary profile table is incomplete")
    counts = profiles.groupby("galaxy_id").size()
    if not (counts == 40).all():
        raise ValueError("Every primary galaxy must have exactly 40 κ rows")

    output_root = args.output_root or (root / "analysis/curve_atlas/outputs")
    tables = output_root / "tables"
    figures = output_root / "figures"
    reports = output_root / "reports"
    tables.mkdir(parents=True, exist_ok=True)

    descriptors = build_descriptors(profiles, summary)
    scores, variance, loadings = build_pca(profiles)
    neighbors = build_neighbors(scores)

    descriptors.to_csv(tables / "curve_descriptors_primary_121.csv", index=False)
    scores.to_csv(tables / "pca_scores_primary_121.csv", index=False)
    variance.to_csv(tables / "pca_explained_variance.csv", index=False)
    loadings.to_csv(tables / "pca_loadings.csv", index=False)
    neighbors.to_csv(tables / "nearest_structural_neighbors.csv", index=False)
    write_figures(profiles, scores, variance, figures)
    write_report(profiles, variance, reports / "CURVE_ATLAS_FOUNDATION_REPORT.md")
    print(f"CURVE_ATLAS_REBUILT tables=5 figures=3 reports=1 output={output_root}")
    return 0


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