analyze_acm_vs_mond.py

#!/usr/bin/env python
"""
Compare the retained ACM galaxy trunk against a fixed-a0 MOND baseline on the
same SPARC sample and geometry-correction chain.
"""

import os
import sys
from pathlib import Path

import numpy as np
import pandas as pd

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

from src.data_loader.load_sparc import load_sparc_galaxy_list, load_sparc_rotation_curve
from src.fitting.fit_sparc import (
    chi2_sparc_galaxy_direct_local,
    fit_sparc_gated_background_model,
)
from src.models.eta_path_integral import eta_local_gated_background_profile
from src.models.inclination_correction import get_inclination_correction_for_galaxy


MAIN7_GALAXIES = [
    "NGC3521",
    "NGC2841",
    "NGC3198",
    "NGC6946",
    "NGC2403",
    "NGC2903",
    "NGC5055",
]

G_CRIT = 3.0e-11
GATE_POWER = 1.0
BG_CONTROL_MODE = "shape_depth"
SHAPE_DEPTH_MODE = "c31_vflat"
SHAPE_C_CRIT = 2.6
SHAPE_DEPTH_CRIT = 2.1
MOND_A0 = 1.2e-10
RIGID_ETA_BASE = 1.0e-29
RIGID_BETA_DENSITY = 4.220349186054929e-28
RIGID_BETA_BG = 7.56705328360024e-29
RIGID_LAMBDA_SUP = 0.1689888140144467


def load_sparc_dict(galaxies=None):
    if galaxies is None:
        table = load_sparc_galaxy_list()
        galaxies = table["Galaxy"].astype(str).tolist()
    sparc_data = {}
    for galaxy in galaxies:
        rc = load_sparc_rotation_curve(galaxy)
        if rc is not None and len(rc) > 5:
            sparc_data[galaxy] = rc
    return sparc_data


def mond_chi2_per_point(rc_data, a0=MOND_A0):
    if "Vbar" not in rc_data.columns:
        return np.nan, 0

    r_kpc = rc_data["Rad"].to_numpy(dtype=float)
    v_obs = rc_data["Vobs"].to_numpy(dtype=float)
    errv = rc_data["errV"].to_numpy(dtype=float)
    vbar = rc_data["Vbar"].to_numpy(dtype=float)

    valid = (
        np.isfinite(r_kpc)
        & np.isfinite(v_obs)
        & np.isfinite(errv)
        & np.isfinite(vbar)
        & (r_kpc > 0)
        & (errv > 0)
        & (vbar > 0)
    )
    if np.count_nonzero(valid) < 3:
        return np.nan, 0

    r_valid = r_kpc[valid]
    r_m = r_valid * 3.08567758e19
    vbar_ms = vbar[valid] * 1000.0
    g_bar = vbar_ms**2 / r_m

    # Standard one-parameter MOND using the "simple" interpolation-equivalent
    # algebraic relation: g = 0.5*(g_N + sqrt(g_N^2 + 4 g_N a0))
    g_mond = 0.5 * (g_bar + np.sqrt(g_bar**2 + 4.0 * g_bar * float(a0)))
    v_mond = np.sqrt(r_m * g_mond) / 1000.0

    galaxy_name = rc_data["Galaxy"].iloc[0] if "Galaxy" in rc_data.columns and len(rc_data) else None
    inc_info = get_inclination_correction_for_galaxy(galaxy_name)
    v_mond = v_mond * float(inc_info.get("k_inc", 1.0))

    chi2 = np.sum(((v_obs[valid] - v_mond) / errv[valid]) ** 2)
    return float(chi2 / np.count_nonzero(valid)), int(np.count_nonzero(valid))


def main():
    output_dir = Path("analysis_outputs")
    output_dir.mkdir(exist_ok=True)

    sparc_full = load_sparc_dict()

    eta_base = RIGID_ETA_BASE
    beta_density = RIGID_BETA_DENSITY
    beta_bg = RIGID_BETA_BG
    lambda_sup = RIGID_LAMBDA_SUP

    rows = []
    for galaxy, rc in sparc_full.items():
        acm_profile = eta_local_gated_background_profile(
            rc,
            eta_base=eta_base,
            beta_density=beta_density,
            beta_bg=beta_bg,
            lambda_sup=lambda_sup,
            galaxy_name=galaxy,
            g_crit=G_CRIT,
            gate_power=GATE_POWER,
            bg_control_mode=BG_CONTROL_MODE,
            shape_depth_mode=SHAPE_DEPTH_MODE,
            shape_c_crit=SHAPE_C_CRIT,
            shape_depth_crit=SHAPE_DEPTH_CRIT,
        )
        if acm_profile is None:
            continue

        acm_cpp = chi2_sparc_galaxy_direct_local(rc, acm_profile) / max(1, len(rc))
        mond_cpp, n_valid = mond_chi2_per_point(rc, a0=MOND_A0)
        if not np.isfinite(mond_cpp):
            continue

        rows.append(
            {
                "Galaxy": galaxy,
                "n_points": len(rc),
                "n_valid_mond": n_valid,
                "acm_cpp": float(acm_cpp),
                "mond_cpp": float(mond_cpp),
                "delta_cpp_mond_minus_acm": float(mond_cpp - acm_cpp),
            }
        )

    per_galaxy = pd.DataFrame(rows).sort_values("delta_cpp_mond_minus_acm", ascending=False)
    per_galaxy.to_csv(output_dir / "acm_vs_mond_per_galaxy.csv", index=False)

    def summarize(series):
        arr = np.asarray(series, dtype=float)
        return float(np.nanmedian(arr)), float(np.nanmean(arr))

    acm_median, acm_mean = summarize(per_galaxy["acm_cpp"])
    mond_median, mond_mean = summarize(per_galaxy["mond_cpp"])

    summary = pd.DataFrame(
        [
            {
                "sample_n_galaxies": int(len(per_galaxy)),
                "mond_a0": MOND_A0,
                "acm_eta_base": eta_base,
                "acm_beta_density": beta_density,
                "acm_beta_bg": beta_bg,
                "acm_lambda_sup": lambda_sup,
                "acm_median_cpp": acm_median,
                "acm_mean_cpp": acm_mean,
                "mond_median_cpp": mond_median,
                "mond_mean_cpp": mond_mean,
                "median_cpp_gain_mond_minus_acm": mond_median - acm_median,
                "mean_cpp_gain_mond_minus_acm": mond_mean - acm_mean,
                "n_acm_better": int(np.sum(per_galaxy["delta_cpp_mond_minus_acm"] > 0)),
                "n_mond_better": int(np.sum(per_galaxy["delta_cpp_mond_minus_acm"] < 0)),
            }
        ]
    )
    summary.to_csv(output_dir / "acm_vs_mond_summary.csv", index=False)

    print("Saved:")
    print(output_dir / "acm_vs_mond_summary.csv")
    print(output_dir / "acm_vs_mond_per_galaxy.csv")
    print(f"ACM median_cpp  = {acm_median:.3f}")
    print(f"MOND median_cpp = {mond_median:.3f}")
    print(f"gain (MOND-ACM) = {mond_median - acm_median:.3f}")


if __name__ == "__main__":
    main()