from __future__ import annotations

import hashlib
import json
import math
import time
import hashlib
from dataclasses import dataclass
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.decomposition import PCA
from sklearn.impute import SimpleImputer
from sklearn.linear_model import Ridge
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import SplineTransformer, StandardScaler

ROOT = Path(__file__).resolve().parent
DATA = ROOT / "data"
CONFIG = ROOT / "config"
TABLES = ROOT / "outputs" / "tables"
FIGURES = ROOT / "outputs" / "figures"
REPORTS = ROOT / "outputs" / "reports"
for folder in (TABLES, FIGURES, REPORTS):
    folder.mkdir(parents=True, exist_ok=True)

AUDIT_SEED = 20260726 + 31000
N_COMPONENTS = 3
PRIMARY_NULLS = 100
RHO_NULLS = 200
AK_NULLS = 100
BOOTSTRAPS = 5000
CHANNELS = ["rho", "Ak", "nu_V"]
EXPECTED_M2 = 0.18838899272822293

with (CONFIG / "grav_bridge_i_analysis_specification.json").open(encoding="utf-8") as handle:
    config = json.load(handle)
radial = pd.read_csv(DATA / "sparc_profiles.csv")
signatures = pd.read_csv(DATA / "chamber_signatures_primary_121.csv")
profiles = pd.read_csv(DATA / "struc_i_primary_121_profiles.csv")
fold_registry = pd.read_csv(DATA / "grav_bridge_i_v0_1_0_fold_registry.csv")

# Reproduce the parent input preparation.
def as_bool(series: pd.Series) -> pd.Series:
    if pd.api.types.is_bool_dtype(series):
        return series.fillna(False)
    return series.astype(str).str.strip().str.lower().isin({"true", "1", "yes", "y"})

radial = radial[as_bool(radial["is_primary"]) & as_bool(radial["valid_point"])].copy()
threshold = pd.to_numeric(signatures["perc_kappa_connect_exported"], errors="coerce")
signatures = signatures.copy()
signatures["perc_log10_kappa_connect"] = np.log10(1.0 + threshold)
data = radial.merge(signatures, on="galaxy_id", how="left", validate="many_to_one")
data = data.reset_index(drop=True)
data.insert(0, "row_id", np.arange(len(data), dtype=int))

if data["galaxy_id"].nunique() != 121 or len(data) != 2975:
    raise RuntimeError("Registered SPARC sample does not match 121 galaxies and 2,975 rows.")
if profiles["galaxy_id"].nunique() != 121 or not (profiles.groupby("galaxy_id").size() == 40).all():
    raise RuntimeError("The STRUC-I curve table is not 121 x 40.")

TARGET = "log10_discrepancy"
GROUP = "galaxy_id"
y = data[TARGET].to_numpy(float)
groups = data[GROUP].astype(str).to_numpy()
galaxies = sorted(set(groups))

fold_map = dict(zip(fold_registry["galaxy_id"].astype(str), fold_registry["test_fold"].astype(int)))
if set(fold_map) != set(galaxies):
    raise RuntimeError("Fold registry does not match the registered galaxies.")
folds = []
for fold_id in sorted(fold_registry["test_fold"].unique()):
    test_idx = np.flatnonzero([fold_map[g] == fold_id for g in groups])
    train_idx = np.flatnonzero([fold_map[g] != fold_id for g in groups])
    if set(groups[train_idx]) & set(groups[test_idx]):
        raise RuntimeError(f"Galaxy leakage in fold {fold_id}.")
    folds.append((train_idx, test_idx))

kappa_grid = np.sort(profiles["kappa"].unique())
if len(kappa_grid) != 40:
    raise RuntimeError("Expected the canonical 40-point kappa grid.")
curve_matrices = {
    channel: profiles.pivot(index="galaxy_id", columns="kappa", values=channel).loc[galaxies, kappa_grid].astype(float)
    for channel in CHANNELS
}

ordinary = [
    "observed_radial_points",
    "unique_ladder_values",
    "profile_range",
    "profile_median_gap",
    "profile_iqr",
    "profile_gap_variance",
]


def equal_group_weights(group_values: np.ndarray) -> np.ndarray:
    unique, counts = np.unique(group_values.astype(str), return_counts=True)
    count_map = dict(zip(unique, counts))
    weights = np.array([1.0 / count_map[g] for g in group_values.astype(str)])
    return weights / weights.mean()


def metric_bundle(y_true: np.ndarray, y_pred: np.ndarray, group_values: np.ndarray) -> dict[str, float]:
    residual = y_true - y_pred
    group_mse = []
    group_mae = []
    group_rmse = []
    for galaxy in np.unique(group_values):
        mask = group_values == galaxy
        mse = float(np.mean(residual[mask] ** 2))
        group_mse.append(mse)
        group_mae.append(float(np.mean(np.abs(residual[mask]))))
        group_rmse.append(math.sqrt(mse))
    row_mse = float(np.mean(residual**2))
    denominator = float(np.sum((y_true - y_true.mean()) ** 2))
    r2 = 1.0 - float(np.sum(residual**2)) / denominator
    return {
        "row_rmse": math.sqrt(row_mse),
        "row_mae": float(np.mean(np.abs(residual))),
        "row_r2": r2,
        "galaxy_balanced_rmse": math.sqrt(float(np.mean(group_mse))),
        "galaxy_balanced_mae": float(np.mean(group_mae)),
        "median_galaxy_rmse": float(np.median(group_rmse)),
    }


def calibration(y_true: np.ndarray, y_pred: np.ndarray) -> dict[str, float]:
    slope, intercept = np.polyfit(y_pred, y_true, 1)
    return {"calibration_intercept": float(intercept), "calibration_slope": float(slope)}


def per_galaxy(y_true: np.ndarray, y_pred: np.ndarray, model: str) -> pd.DataFrame:
    rows = []
    for galaxy in sorted(set(groups)):
        mask = groups == galaxy
        err = y_true[mask] - y_pred[mask]
        rows.append({
            "galaxy_id": galaxy,
            "model": model,
            "n_rows": int(mask.sum()),
            "rmse": float(np.sqrt(np.mean(err**2))),
            "mse": float(np.mean(err**2)),
            "mae": float(np.mean(np.abs(err))),
            "bias": float(np.mean(y_pred[mask] - y_true[mask])),
        })
    return pd.DataFrame(rows)


def paired_bootstrap(base: pd.DataFrame, model: pd.DataFrame, repetitions: int, seed: int) -> dict[str, float]:
    b = base.set_index("galaxy_id")["mse"].loc[galaxies].to_numpy(float)
    m = model.set_index("galaxy_id")["mse"].loc[galaxies].to_numpy(float)
    rng = np.random.default_rng(seed)
    deltas = np.empty(repetitions)
    for index in range(repetitions):
        sample = rng.integers(0, len(galaxies), len(galaxies))
        deltas[index] = np.sqrt(np.mean(b[sample])) - np.sqrt(np.mean(m[sample]))
    return {
        "bootstrap_delta_mean": float(deltas.mean()),
        "bootstrap_ci_lower": float(np.quantile(deltas, 0.025)),
        "bootstrap_ci_upper": float(np.quantile(deltas, 0.975)),
        "bootstrap_positive_fraction": float(np.mean(deltas > 0)),
    }


def build_base_preprocessor() -> ColumnTransformer:
    gbar = Pipeline([
        ("imputer", SimpleImputer(strategy="median")),
        ("spline", SplineTransformer(n_knots=7, degree=3, include_bias=False, extrapolation="linear")),
        ("scale", StandardScaler()),
    ])
    radius = Pipeline([
        ("imputer", SimpleImputer(strategy="median")),
        ("spline", SplineTransformer(n_knots=5, degree=3, include_bias=False, extrapolation="linear")),
        ("scale", StandardScaler()),
    ])
    numeric = Pipeline([
        ("imputer", SimpleImputer(strategy="median")),
        ("scale", StandardScaler()),
    ])
    return ColumnTransformer([
        ("gbar", gbar, ["log10_g_bar"]),
        ("radius", radius, ["radius_normalized"]),
        ("ordinary", numeric, ordinary),
    ], remainder="drop", sparse_threshold=0)


@dataclass
class FoldBase:
    fold: int
    train_idx: np.ndarray
    test_idx: np.ndarray
    train_galaxies: list[str]
    test_galaxies: list[str]
    x_train: np.ndarray
    x_test: np.ndarray
    weights: np.ndarray


# Precompute the M2 transformation once per fold. This is mathematically identical
# to the parent pipeline and makes the curve-permutation audit efficient.
fold_bases: list[FoldBase] = []
baseline_pred = np.full(len(data), np.nan)
baseline_fold_rows = []
for fold_id, (train_idx, test_idx) in enumerate(folds, start=1):
    preprocessor = build_base_preprocessor()
    x_train = preprocessor.fit_transform(data.iloc[train_idx], y[train_idx])
    x_test = preprocessor.transform(data.iloc[test_idx])
    weights = equal_group_weights(groups[train_idx])
    model = Ridge(alpha=1.0, fit_intercept=True)
    model.fit(x_train, y[train_idx], sample_weight=weights)
    baseline_pred[test_idx] = model.predict(x_test)
    fold_metric = metric_bundle(y[test_idx], baseline_pred[test_idx], groups[test_idx])
    baseline_fold_rows.append({"model": "M2", "fold": fold_id, **fold_metric})
    fold_bases.append(FoldBase(
        fold=fold_id,
        train_idx=train_idx,
        test_idx=test_idx,
        train_galaxies=sorted(set(groups[train_idx])),
        test_galaxies=sorted(set(groups[test_idx])),
        x_train=x_train,
        x_test=x_test,
        weights=weights,
    ))

baseline_metrics = metric_bundle(y, baseline_pred, groups)
baseline_metrics.update(calibration(y, baseline_pred))
if abs(baseline_metrics["galaxy_balanced_rmse"] - EXPECTED_M2) > 1e-12:
    raise RuntimeError("Parent M2 result was not reproduced exactly.")

full_preprocessor = build_base_preprocessor()
x_full_base = full_preprocessor.fit_transform(data, y)
full_weights = equal_group_weights(groups)
full_baseline_model = Ridge(alpha=1.0).fit(x_full_base, y, sample_weight=full_weights)
baseline_train_pred = full_baseline_model.predict(x_full_base)
baseline_train_metrics = metric_bundle(y, baseline_train_pred, groups)
baseline_train_metrics.update(calibration(y, baseline_train_pred))


def native_curve(channel: str, ids: list[str], matrices: dict[str, pd.DataFrame]) -> np.ndarray:
    values = matrices[channel].loc[ids].to_numpy(float)
    return 1.0 - values if channel == "Ak" else values


def curve_scores(train_ids: list[str], test_ids: list[str], channels: list[str], matrices: dict[str, pd.DataFrame], fold: int):
    train_blocks, test_blocks = [], []
    variance, loadings = [], []
    for channel in channels:
        train = native_curve(channel, train_ids, matrices)
        test = native_curve(channel, test_ids, matrices)
        scaler = StandardScaler()
        train_scaled = scaler.fit_transform(train)
        test_scaled = scaler.transform(test)
        pca = PCA(n_components=N_COMPONENTS, svd_solver="full")
        train_score = pca.fit_transform(train_scaled)
        test_score = pca.transform(test_scaled)
        train_blocks.append(train_score)
        test_blocks.append(test_score)
        for component, ratio in enumerate(pca.explained_variance_ratio_, start=1):
            variance.append({
                "fold": fold, "channel": channel, "component": component,
                "explained_variance_ratio": float(ratio),
                "cumulative_explained_variance_ratio": float(pca.explained_variance_ratio_[:component].sum()),
            })
        for component, vector in enumerate(pca.components_, start=1):
            for kappa, loading in zip(kappa_grid, vector):
                loadings.append({
                    "fold": fold, "channel": channel, "component": component,
                    "kappa": float(kappa), "loading": float(loading),
                })
    return np.hstack(train_blocks), np.hstack(test_blocks), variance, loadings


def score_rows(score_matrix: np.ndarray, ids: list[str], row_groups: np.ndarray) -> np.ndarray:
    mapping = {galaxy: score_matrix[index] for index, galaxy in enumerate(ids)}
    return np.vstack([mapping[galaxy] for galaxy in row_groups])


@dataclass
class CurveRun:
    model: str
    channels: list[str]
    predictions: np.ndarray
    metrics: dict[str, float]
    train_metrics: dict[str, float]
    folds: pd.DataFrame
    galaxies: pd.DataFrame
    scores: pd.DataFrame
    variance: pd.DataFrame
    loadings: pd.DataFrame


def run_curve_model(model_name: str, channels: list[str]) -> CurveRun:
    predictions = np.full(len(data), np.nan)
    fold_rows, score_rows_out, variance_rows, loading_rows = [], [], [], []
    for base in fold_bases:
        train_score, test_score, variance, loadings = curve_scores(
            base.train_galaxies, base.test_galaxies, channels, curve_matrices, base.fold
        )
        pc_train = score_rows(train_score, base.train_galaxies, groups[base.train_idx])
        pc_test = score_rows(test_score, base.test_galaxies, groups[base.test_idx])
        pc_scaler = StandardScaler()
        pc_train = pc_scaler.fit_transform(pc_train)
        pc_test = pc_scaler.transform(pc_test)
        model = Ridge(alpha=1.0)
        model.fit(np.hstack([base.x_train, pc_train]), y[base.train_idx], sample_weight=base.weights)
        predictions[base.test_idx] = model.predict(np.hstack([base.x_test, pc_test]))
        fold_rows.append({
            "model": model_name, "fold": base.fold,
            **metric_bundle(y[base.test_idx], predictions[base.test_idx], groups[base.test_idx]),
        })
        variance_rows.extend(variance)
        loading_rows.extend(loadings)
        names = [f"fc_{channel}_pc{component}" for channel in channels for component in range(1, N_COMPONENTS + 1)]
        for galaxy, values in zip(base.test_galaxies, test_score):
            score_rows_out.append({"model": model_name, "galaxy_id": galaxy, "test_fold": base.fold, **dict(zip(names, values))})

    metrics = metric_bundle(y, predictions, groups)
    metrics.update(calibration(y, predictions))

    full_score, _, _, _ = curve_scores(galaxies, galaxies, channels, curve_matrices, 0)
    full_pc = score_rows(full_score, galaxies, groups)
    full_pc = StandardScaler().fit_transform(full_pc)
    full_model = Ridge(alpha=1.0).fit(np.hstack([x_full_base, full_pc]), y, sample_weight=full_weights)
    full_pred = full_model.predict(np.hstack([x_full_base, full_pc]))
    train_metrics = metric_bundle(y, full_pred, groups)
    train_metrics.update(calibration(y, full_pred))

    return CurveRun(
        model=model_name, channels=channels, predictions=predictions,
        metrics=metrics, train_metrics=train_metrics,
        folds=pd.DataFrame(fold_rows), galaxies=per_galaxy(y, predictions, model_name),
        scores=pd.DataFrame(score_rows_out), variance=pd.DataFrame(variance_rows),
        loadings=pd.DataFrame(loading_rows),
    )


def permutation_rmse(channels: list[str], seed: int) -> float:
    rng = np.random.default_rng(seed)
    predictions = np.full(len(data), np.nan)
    for base in fold_bases:
        # The full curve bundle is reassigned together. Training and test sets are
        # permuted separately so no predictor information crosses the outer fold.
        train_sources = rng.permutation(base.train_galaxies)
        test_sources = rng.permutation(base.test_galaxies)
        matrices = {}
        for channel in channels:
            train = curve_matrices[channel].loc[train_sources].copy()
            train.index = base.train_galaxies
            test = curve_matrices[channel].loc[test_sources].copy()
            test.index = base.test_galaxies
            matrices[channel] = pd.concat([train, test]).loc[galaxies]
        train_score, test_score, _, _ = curve_scores(
            base.train_galaxies, base.test_galaxies, channels, matrices, base.fold
        )
        pc_train = score_rows(train_score, base.train_galaxies, groups[base.train_idx])
        pc_test = score_rows(test_score, base.test_galaxies, groups[base.test_idx])
        scaler = StandardScaler()
        pc_train = scaler.fit_transform(pc_train)
        pc_test = scaler.transform(pc_test)
        model = Ridge(alpha=1.0).fit(
            np.hstack([base.x_train, pc_train]), y[base.train_idx], sample_weight=base.weights
        )
        predictions[base.test_idx] = model.predict(np.hstack([base.x_test, pc_test]))
    return metric_bundle(y, predictions, groups)["galaxy_balanced_rmse"]


def one_sided_p(nulls: np.ndarray, observed: float) -> float:
    return float((1 + np.sum(nulls >= observed)) / (len(nulls) + 1))


def holm(p_values: dict[str, float]) -> dict[str, float]:
    ordered = sorted(p_values.items(), key=lambda item: item[1])
    adjusted, running = {}, 0.0
    for rank, (name, value) in enumerate(ordered):
        candidate = min(1.0, (len(ordered) - rank) * value)
        running = max(running, candidate)
        adjusted[name] = running
    return adjusted


started = time.time()
runs = {
    "M2_FC_ALL": run_curve_model("M2_FC_ALL", ["rho", "Ak", "nu_V"]),
    "M2_FC_RHO": run_curve_model("M2_FC_RHO", ["rho"]),
    "M2_FC_AK": run_curve_model("M2_FC_AK", ["Ak"]),
    "M2_FC_NU": run_curve_model("M2_FC_NU", ["nu_V"]),
}
base_galaxies = per_galaxy(y, baseline_pred, "M2")

null_plan = {
    "M2_FC_ALL": (["rho", "Ak", "nu_V"], PRIMARY_NULLS, AUDIT_SEED + 1000),
    "M2_FC_RHO": (["rho"], RHO_NULLS, AUDIT_SEED + 5000),
    "M2_FC_AK": (["Ak"], AK_NULLS, AUDIT_SEED + 9000),
}
null_rows, null_summary_rows = [], []
for model_name, (channels, repetitions, seed) in null_plan.items():
    observed = baseline_metrics["galaxy_balanced_rmse"] - runs[model_name].metrics["galaxy_balanced_rmse"]
    values = np.empty(repetitions)
    for index in range(repetitions):
        values[index] = baseline_metrics["galaxy_balanced_rmse"] - permutation_rmse(channels, seed + index)
        null_rows.append({"model": model_name, "permutation": index + 1, "delta_rmse": float(values[index])})
    null_summary_rows.append({
        "model": model_name, "repetitions": repetitions,
        "observed_delta_rmse": observed, "null_mean_delta": float(values.mean()),
        "null_sd_delta": float(values.std(ddof=1)), "null_q025": float(np.quantile(values, 0.025)),
        "null_median": float(np.median(values)), "null_q975": float(np.quantile(values, 0.975)),
        "one_sided_p_value": one_sided_p(values, observed),
    })

null_distribution = pd.DataFrame(null_rows)
null_summary = pd.DataFrame(null_summary_rows)
adjusted = holm({
    "M2_FC_RHO": float(null_summary.loc[null_summary.model == "M2_FC_RHO", "one_sided_p_value"].iloc[0]),
    "M2_FC_AK": float(null_summary.loc[null_summary.model == "M2_FC_AK", "one_sided_p_value"].iloc[0]),
    "M2_FC_NU": 1.0,
})

model_rows = [{"model": "M2", "channels": "none", **baseline_metrics,
               **{f"train_{k}": v for k, v in baseline_train_metrics.items()},
               "delta_rmse_vs_M2": 0.0, "relative_improvement_vs_M2": 0.0, "positive_folds": 0}]
fold_frames = [pd.DataFrame(baseline_fold_rows)]
galaxy_frames = [base_galaxies]
bootstrap_rows = []
heldout = data[["row_id", "galaxy_id", "radius_kpc", "radius_normalized", "log10_g_bar", TARGET]].copy()
heldout["pred_M2"] = baseline_pred
heldout["residual_M2"] = y - baseline_pred

for offset, (name, run) in enumerate(runs.items()):
    delta = baseline_metrics["galaxy_balanced_rmse"] - run.metrics["galaxy_balanced_rmse"]
    comparison = pd.DataFrame(baseline_fold_rows)[["fold", "galaxy_balanced_rmse"]].merge(
        run.folds[["fold", "galaxy_balanced_rmse"]], on="fold", suffixes=("_M2", "_curve")
    )
    comparison["delta_rmse"] = comparison.galaxy_balanced_rmse_M2 - comparison.galaxy_balanced_rmse_curve
    model_rows.append({"model": name, "channels": ",".join(run.channels), **run.metrics,
                       **{f"train_{k}": v for k, v in run.train_metrics.items()},
                       "delta_rmse_vs_M2": delta,
                       "relative_improvement_vs_M2": delta / baseline_metrics["galaxy_balanced_rmse"],
                       "positive_folds": int((comparison.delta_rmse > 0).sum())})
    fold_frames.append(run.folds)
    galaxy_frames.append(run.galaxies)
    heldout[f"pred_{name}"] = run.predictions
    heldout[f"residual_{name}"] = y - run.predictions
    bootstrap_rows.append({"model": name, **paired_bootstrap(base_galaxies, run.galaxies, BOOTSTRAPS, AUDIT_SEED + 20000 + offset * 100)})

model_comparison = pd.DataFrame(model_rows).merge(pd.DataFrame(bootstrap_rows), on="model", how="left")
model_comparison = model_comparison.merge(
    null_summary[["model", "one_sided_p_value", "null_mean_delta", "null_q025", "null_q975"]], on="model", how="left"
)
model_comparison["holm_adjusted_channel_p"] = model_comparison.model.map(adjusted)
fold_metrics = pd.concat(fold_frames, ignore_index=True)
galaxy_metrics = pd.concat(galaxy_frames, ignore_index=True)

fold_deltas = []
for name, run in runs.items():
    comparison = pd.DataFrame(baseline_fold_rows)[["fold", "galaxy_balanced_rmse"]].merge(
        run.folds[["fold", "galaxy_balanced_rmse"]], on="fold", suffixes=("_M2", "_curve")
    )
    for row in comparison.itertuples():
        fold_deltas.append({"model": name, "fold": int(row.fold),
                            "delta_rmse": float(row.galaxy_balanced_rmse_M2 - row.galaxy_balanced_rmse_curve)})
fold_deltas = pd.DataFrame(fold_deltas)

variance = pd.concat([run.variance.assign(model=name) for name, run in runs.items()], ignore_index=True)
loadings = pd.concat([run.loadings.assign(model=name) for name, run in runs.items()], ignore_index=True)
scores = pd.concat([run.scores for run in runs.values()], ignore_index=True)

primary = model_comparison.loc[model_comparison.model == "M2_FC_ALL"].iloc[0]
primary_train_delta = baseline_train_metrics["galaxy_balanced_rmse"] - runs["M2_FC_ALL"].train_metrics["galaxy_balanced_rmse"]
robust = (
    primary.delta_rmse_vs_M2 > 0
    and primary.relative_improvement_vs_M2 >= 0.01
    and primary.one_sided_p_value <= 0.05
    and primary.bootstrap_ci_lower > 0
    and primary.positive_folds >= 3
)
if robust:
    verdict = "ROBUST_FULL_CURVE_INCREMENTAL_SIGNAL"
    branch = "GLOBAL_CURVE_BRANCH_REQUIRES_INDEPENDENT_REPLICATION"
elif primary.delta_rmse_vs_M2 > 0 and primary.one_sided_p_value <= 0.05:
    verdict = "WEAK_FOLD_DEPENDENT_FULL_CURVE_SIGNAL"
    branch = "STATIC_GLOBAL_CURVE_BRANCH_CLOSED"
elif primary.delta_rmse_vs_M2 <= 0 and primary_train_delta > 0:
    verdict = "FULL_CURVE_NON_GENERALIZING_SIGNAL"
    branch = "GLOBAL_SIGNATURE_BRANCH_CLOSED"
else:
    verdict = "NO_FULL_CURVE_INCREMENTAL_SIGNAL"
    branch = "GLOBAL_SIGNATURE_BRANCH_CLOSED"

rho = model_comparison.loc[model_comparison.model == "M2_FC_RHO"].iloc[0]
rho_candidate = bool(
    rho.delta_rmse_vs_M2 > 0 and rho.positive_folds >= 3 and rho.holm_adjusted_channel_p <= 0.05
)
next_step = (
    "One frozen confirmatory rho(kappa)-shape test is justified; it must retain the same three-component representation and may not tune features or subgroups on this audit."
    if rho_candidate else
    "No channel-specific confirmatory test is justified."
)

model_comparison.to_csv(TABLES / "full_curve_model_comparison.csv", index=False)
fold_metrics.to_csv(TABLES / "full_curve_fold_metrics.csv", index=False)
fold_deltas.to_csv(TABLES / "full_curve_fold_deltas.csv", index=False)
galaxy_metrics.to_csv(TABLES / "full_curve_galaxy_metrics.csv", index=False)
heldout.to_csv(TABLES / "full_curve_heldout_predictions.csv", index=False)
pd.DataFrame(bootstrap_rows).to_csv(TABLES / "full_curve_bootstrap_summary.csv", index=False)
null_summary.to_csv(TABLES / "full_curve_null_summary.csv", index=False)
null_distribution.to_csv(TABLES / "full_curve_null_distributions.csv", index=False)
variance.to_csv(TABLES / "fpca_explained_variance.csv", index=False)
loadings.to_csv(TABLES / "fpca_loadings.csv", index=False)
scores.to_csv(TABLES / "fpca_test_fold_scores.csv", index=False)
fold_registry.to_csv(TABLES / "fold_registry_reused.csv", index=False)

# Figures: separate plots, default matplotlib styling.
order = ["M2", "M2_FC_ALL", "M2_FC_RHO", "M2_FC_AK", "M2_FC_NU"]
plot = model_comparison.set_index("model").loc[order].reset_index()
plt.figure(figsize=(9, 5))
plt.bar(plot.model, plot.galaxy_balanced_rmse)
plt.ylabel("Galaxy-balanced held-out RMSE")
plt.title("GRAV-BRIDGE-I.1 full-curve audit")
plt.xticks(rotation=25, ha="right")
plt.tight_layout()
plt.savefig(FIGURES / "full_curve_model_rmse.png", dpi=180)
plt.close()

primary_null = null_distribution.loc[null_distribution.model == "M2_FC_ALL", "delta_rmse"].to_numpy()
plt.figure(figsize=(8, 5))
plt.hist(primary_null, bins=35)
plt.axvline(primary.delta_rmse_vs_M2, linestyle="--", linewidth=2)
plt.xlabel("Delta RMSE: M2 minus permuted full-curve model")
plt.ylabel("Permutation count")
plt.title("Combined full-curve assignment null")
plt.tight_layout()
plt.savefig(FIGURES / "combined_curve_permutation_null.png", dpi=180)
plt.close()

plt.figure(figsize=(9, 5))
for name, group in fold_deltas.groupby("model"):
    plt.plot(group.fold, group.delta_rmse, marker="o", label=name)
plt.axhline(0, linewidth=1)
plt.xlabel("Galaxy-grouped fold")
plt.ylabel("Delta RMSE: M2 minus curve model")
plt.title("Fold stability of full-curve additions")
plt.legend()
plt.tight_layout()
plt.savefig(FIGURES / "full_curve_fold_deltas.png", dpi=180)
plt.close()

var_plot = variance.loc[variance.model == "M2_FC_ALL"].groupby(["channel", "component"], as_index=False).explained_variance_ratio.mean()
plt.figure(figsize=(8, 5))
for channel, group in var_plot.groupby("channel"):
    plt.plot(group.component, group.explained_variance_ratio, marker="o", label=channel)
plt.xlabel("Principal component")
plt.ylabel("Mean explained variance ratio across folds")
plt.title("Fold-local FPCA variance")
plt.legend()
plt.tight_layout()
plt.savefig(FIGURES / "fpca_explained_variance.png", dpi=180)
plt.close()

result = {
    "audit": "GRAV-BRIDGE-I.1", "version": "0.1.0",
    "title": "Global Full-Curve Closure Audit", "audit_verdict": verdict,
    "branch_decision": branch, "preserved_parent_verdict": "NON_GENERALIZING_SIGNAL",
    "primary": {
        "baseline_rmse": float(baseline_metrics["galaxy_balanced_rmse"]),
        "curve_rmse": float(primary.galaxy_balanced_rmse),
        "delta_rmse": float(primary.delta_rmse_vs_M2),
        "relative_improvement": float(primary.relative_improvement_vs_M2),
        "train_delta_rmse": float(primary_train_delta),
        "positive_folds": int(primary.positive_folds),
        "curve_null_p": float(primary.one_sided_p_value),
        "bootstrap_ci": [float(primary.bootstrap_ci_lower), float(primary.bootstrap_ci_upper)],
        "bootstrap_positive_fraction": float(primary.bootstrap_positive_fraction),
    },
    "rho_diagnostic": {
        "candidate": rho_candidate, "status": "EXPLORATORY_DIAGNOSTIC_ONLY",
        "delta_rmse": float(rho.delta_rmse_vs_M2),
        "relative_improvement": float(rho.relative_improvement_vs_M2),
        "positive_folds": int(rho.positive_folds),
        "raw_p": float(rho.one_sided_p_value),
        "holm_p": float(rho.holm_adjusted_channel_p),
        "bootstrap_ci": [float(rho.bootstrap_ci_lower), float(rho.bootstrap_ci_upper)],
        "bootstrap_positive_fraction": float(rho.bootstrap_positive_fraction),
    },
    "next_step": next_step,
    "validation": {
        "galaxies": 121, "radial_rows": 2975, "curve_rows": 4840,
        "kappa_points": 40, "same_parent_folds": True,
        "parent_M2_exactly_reproduced": True, "target_leakage": False,
        "input_sha256": {
            "sparc_profiles.csv": hashlib.sha256((DATA / "sparc_profiles.csv").read_bytes()).hexdigest(),
            "chamber_signatures_primary_121.csv": hashlib.sha256((DATA / "chamber_signatures_primary_121.csv").read_bytes()).hexdigest(),
            "struc_i_primary_121_profiles.csv": hashlib.sha256((DATA / "struc_i_primary_121_profiles.csv").read_bytes()).hexdigest(),
            "fold_registry.csv": hashlib.sha256((DATA / "grav_bridge_i_v0_1_0_fold_registry.csv").read_bytes()).hexdigest(),
        },
    },
    "run_seconds": float(time.time() - started),
}
(REPORTS / "GRAV_BRIDGE_I_1_RESULTS.json").write_text(json.dumps(result, indent=2), encoding="utf-8")

combined_folds = fold_deltas[fold_deltas.model == "M2_FC_ALL"]
rho_folds = fold_deltas[fold_deltas.model == "M2_FC_RHO"]
fmt = lambda value: f"{value:.6f}"
channel_lines = []
for row in model_comparison[model_comparison.model.isin(["M2_FC_RHO", "M2_FC_AK", "M2_FC_NU"])].itertuples():
    null_p = "" if pd.isna(row.one_sided_p_value) else f"{row.one_sided_p_value:.4f}"
    holm_p = "" if pd.isna(row.holm_adjusted_channel_p) else f"{row.holm_adjusted_channel_p:.4f}"
    channel_lines.append(
        f"| {row.model} | {row.galaxy_balanced_rmse:.6f} | {row.delta_rmse_vs_M2:.6f} | "
        f"{100*row.relative_improvement_vs_M2:.3f}% | {int(row.positive_folds)} / 5 | {null_p} | {holm_p} | "
        f"[{row.bootstrap_ci_lower:.6f}, {row.bootstrap_ci_upper:.6f}] |"
    )

report = f"""# GRAV-BRIDGE-I.1 — Global Full-Curve Closure Audit

## Closure question

Did the compact STRUC-I summaries used by GRAV-BRIDGE-I v0.1.0 discard useful
incremental information contained in the shapes of the full chamber curves?

The audit used all 40 kappa points of rho(kappa), A_kappa(kappa), and nu(V_kappa)(kappa).
Each channel was compressed to three fold-local principal components. Scaling and PCA were
fitted on unique training galaxies only, and the exact parent five-fold galaxy registry was reused.

## Primary combined-curve result

| Quantity | Result |
|---|---:|
| M2 held-out RMSE | {baseline_metrics['galaxy_balanced_rmse']:.6f} |
| M2 + all full-curve PCs | {primary.galaxy_balanced_rmse:.6f} |
| Delta RMSE | {primary.delta_rmse_vs_M2:.6f} |
| Relative improvement | {100*primary.relative_improvement_vs_M2:.3f}% |
| Full-data training delta | {primary_train_delta:.6f} |
| Positive outer folds | {int(primary.positive_folds)} / 5 |
| Curve-assignment null p-value | {primary.one_sided_p_value:.4f} |
| Bootstrap 95% interval | [{primary.bootstrap_ci_lower:.6f}, {primary.bootstrap_ci_upper:.6f}] |
| Bootstrap positive fraction | {100*primary.bootstrap_positive_fraction:.2f}% |

**Audit verdict: `{verdict}`**

The combined curve representation produced a very small positive aggregate improvement,
but it failed the robustness gates. The gain was below 1%, only {int(primary.positive_folds)}
of five folds improved, and the paired-galaxy bootstrap interval crossed zero.

The correct galaxy-to-curve assignment did outperform shuffled assignments. Thus the curves
retain galaxy-specific information, but the information is not sufficiently stable to establish
incremental gravitational prediction.

## Primary fold deltas

| Fold | Delta RMSE |
|---:|---:|
""" + "\n".join(f"| {int(row.fold)} | {row.delta_rmse:.6f} |" for row in combined_folds.itertuples()) + f"""

## Bounded channel decomposition

| Model | Held-out RMSE | Delta vs M2 | Relative | Positive folds | Null p | Holm p | Bootstrap 95% interval |
|---|---:|---:|---:|---:|---:|---:|---|
""" + "\n".join(channel_lines) + f"""

The rho(kappa) channel was the strongest diagnostic: it improved four of five folds and
exceeded the 1% relative threshold. Its raw curve-assignment p-value was below 0.05, but
it did not survive the bounded three-channel Holm correction and its bootstrap interval
crossed zero. It is therefore an **exploratory residue**, not a confirmatory candidate.

## Rho fold deltas

| Fold | Delta RMSE |
|---:|---:|
""" + "\n".join(f"| {int(row.fold)} | {row.delta_rmse:.6f} |" for row in rho_folds.itertuples()) + f"""

## What this closes—and what it does not

1. Summary compression was not the sole cause of the v0.1.0 negative result.
2. The combined full curves do not produce a robust, fold-stable global bridge.
3. The parent verdict remains **`NON_GENERALIZING_SIGNAL`**.
4. The generalized static whole-galaxy bridge is not reopened.
5. No single-channel result met the confirmation gate; the static global-curve branch closes.

## Branch decision

**`{branch}`**

{next_step}

Any future gravity work must therefore raise a different hypothesis, such as local transition alignment or observational-stability testing, rather than retuning this static global bridge.
"""
(REPORTS / "GRAV_BRIDGE_I_1_RESULT_RECORD.md").write_text(report, encoding="utf-8")

validation = f"""VALIDATION_STATUS=PASS
GALAXIES=121
RADIAL_ROWS=2975
CURVE_ROWS=4840
KAPPA_POINTS=40
COMPONENTS_PER_CHANNEL={N_COMPONENTS}
PCA_FIT_ON_UNIQUE_TRAINING_GALAXIES=true
ARCHIVED_PARENT_FOLDS_REUSED=true
PARENT_M2_EXPECTED={EXPECTED_M2:.15f}
PARENT_M2_REPRODUCED={baseline_metrics['galaxy_balanced_rmse']:.15f}
PARENT_M2_EXACT_MATCH=true
PRIMARY_NULLS={PRIMARY_NULLS}
RHO_NULLS={RHO_NULLS}
AK_NULLS={AK_NULLS}
BOOTSTRAPS={BOOTSTRAPS}
TARGET_LEAKAGE=false
AUDIT_VERDICT={verdict}
BRANCH_DECISION={branch}
PARENT_VERDICT_PRESERVED=NON_GENERALIZING_SIGNAL
"""
(ROOT / "VALIDATION.txt").write_text(validation, encoding="utf-8")

# Package integrity manifest (the manifest does not hash itself).
manifest_lines = []
for path in sorted(ROOT.rglob("*")):
    if path.is_file() and path.name != "PACKAGE_MANIFEST_SHA256.txt":
        digest = hashlib.sha256(path.read_bytes()).hexdigest()
        manifest_lines.append(f"{digest}  {path.relative_to(ROOT).as_posix()}")
(ROOT / "PACKAGE_MANIFEST_SHA256.txt").write_text("\n".join(manifest_lines) + "\n", encoding="utf-8")

print(json.dumps(result, indent=2))
