#!/usr/bin/env python3
"""
DARK MATTER GEOMETRIC CORRELATION TEST
=======================================
Project: Prometheus / Time Ledger Theory / Cipher Cross-Scale Validation
Date: 2026-03-26

HYPOTHESIS: Galaxy rotation curve "excess" (dark matter fraction) correlates
with geometric parameters (pitch angle, spiral type, morphology) because it
is geometric amplification, not missing mass.

DATA:
  - SPARC: 175 galaxies with rotation curves (Lelli, McGaugh, Schombert 2016)
  - Diaz-Garcia 2019: 391 galaxies with pitch angles from S4G
  - Cross-match on galaxy name

OUTPUT-AGNOSTIC. DATA SHOWS WHAT IT SHOWS.
"""

import json
import os
import sys
import warnings
from datetime import datetime

import numpy as np

warnings.filterwarnings('ignore')

# ============================================================
# DATA PARSING
# ============================================================

def parse_sparc_galaxies(filepath):
    """Parse SPARC_Lelli2016c.mrt — galaxy properties table."""
    galaxies = {}
    with open(filepath, 'r') as f:
        lines = f.readlines()

    # Find where data starts (after the LAST separator line)
    data_start = 0
    for i, line in enumerate(lines):
        if line.startswith('---'):
            data_start = i + 1

    for line in lines[data_start:]:
        if len(line.strip()) < 20:
            continue
        try:
            # Use split-based parsing (MRT fixed-width has alignment issues with 2-digit T)
            parts = line.split()
            if len(parts) < 15:
                continue
            name = parts[0]
            T = int(parts[1])
            D = float(parts[2])
            # e_D = parts[3]
            # f_D = parts[4]
            Inc = float(parts[5])
            # e_Inc = parts[6]
            L36 = float(parts[7])
            # e_L36 = parts[8]
            Reff = float(parts[9])
            SBeff = float(parts[10])
            Rdisk = float(parts[11])
            # SBdisk = parts[12]
            MHI = float(parts[13])
            # RHI = parts[14]
            Vflat = float(parts[15]) if len(parts) > 15 else 0
            e_Vflat = float(parts[16]) if len(parts) > 16 else 0
            Q = int(parts[17]) if len(parts) > 17 else 3

            galaxies[name] = {
                'name': name,
                'hubble_type': T,
                'distance_mpc': D,
                'inclination_deg': Inc,
                'luminosity_1e9Lsun': L36,
                'eff_radius_kpc': Reff,
                'eff_SB': SBeff,
                'disk_scale_kpc': Rdisk,
                'HI_mass_1e9Msun': MHI,
                'Vflat_kms': Vflat,
                'e_Vflat': e_Vflat,
                'quality': Q,
            }
        except (ValueError, IndexError):
            continue

    return galaxies


def parse_sparc_rotation_curves(filepath):
    """Parse MassModels_Lelli2016c.mrt — rotation curve data."""
    curves = {}
    with open(filepath, 'r') as f:
        lines = f.readlines()

    # Find data start (after LAST separator)
    data_start = 0
    for i, line in enumerate(lines):
        if line.startswith('---'):
            data_start = i + 1

    for line in lines[data_start:]:
        if len(line.strip()) < 20:
            continue
        try:
            parts = line.split()
            if len(parts) < 7:
                continue
            name = parts[0]
            # parts[1] = D (distance)
            R = float(parts[2])
            Vobs = float(parts[3])
            e_Vobs = float(parts[4])
            Vgas = float(parts[5])
            Vdisk = float(parts[6])
            Vbul = float(parts[7]) if len(parts) > 7 else 0

            if name not in curves:
                curves[name] = {'R': [], 'Vobs': [], 'e_Vobs': [],
                                'Vgas': [], 'Vdisk': [], 'Vbul': []}

            curves[name]['R'].append(R)
            curves[name]['Vobs'].append(Vobs)
            curves[name]['e_Vobs'].append(e_Vobs)
            curves[name]['Vgas'].append(Vgas)
            curves[name]['Vdisk'].append(Vdisk)
            curves[name]['Vbul'].append(Vbul)
        except (ValueError, IndexError):
            continue

    # Convert to numpy
    for name in curves:
        for key in curves[name]:
            curves[name][key] = np.array(curves[name][key])

    return curves


def parse_pitch_angles(filepath):
    """Parse Diaz-Garcia 2019 pitch angle catalog."""
    pitches = {}
    with open(filepath, 'r') as f:
        for line in f:
            if len(line.strip()) < 20:
                continue
            try:
                name = line[0:10].strip()
                spiral_type = line[13].strip() if len(line) > 13 else ''
                flag = int(line[15]) if len(line) > 15 and line[15].strip() else 2
                n_segments = int(line[17]) if len(line) > 17 and line[17].strip() else 0
                Pmean = float(line[22:27].strip()) if line[22:27].strip() else 0

                # Parse error — may be -9999.99 (missing)
                e_Pmean_str = line[28:36].strip() if len(line) > 36 else ''
                e_Pmean = float(e_Pmean_str) if e_Pmean_str and float(e_Pmean_str) > -9000 else None

                Pmedian = float(line[49:54].strip()) if len(line) > 54 and line[49:54].strip() else 0
                Pinner = float(line[58:63].strip()) if len(line) > 63 and line[58:63].strip() else 0

                # Weighted pitch angle
                Pweighted_str = line[76:81].strip() if len(line) > 81 else ''
                Pweighted = float(Pweighted_str) if Pweighted_str else 0

                # Spiral strength (QT)
                QT_str = line[91:99].strip() if len(line) > 99 else ''
                QTspiral = float(QT_str) if QT_str and float(QT_str) > -9000 else None

                # A2 amplitude
                A2_str = line[136:144].strip() if len(line) > 144 else ''
                A2spiral = float(A2_str) if A2_str and float(A2_str) > -9000 else None

                pitches[name] = {
                    'name': name,
                    'spiral_type': spiral_type,  # G=grand-design, M=multi-armed, F=flocculent
                    'quality_flag': flag,
                    'n_segments': n_segments,
                    'pitch_mean_deg': Pmean,
                    'pitch_err_deg': e_Pmean,
                    'pitch_median_deg': Pmedian,
                    'pitch_inner_deg': Pinner,
                    'pitch_weighted_deg': Pweighted,
                    'QT_spiral': QTspiral,
                    'A2_spiral': A2spiral,
                }
            except (ValueError, IndexError):
                continue

    return pitches


# ============================================================
# DARK MATTER FRACTION CALCULATION
# ============================================================

def compute_dm_fraction(curves, ML_disk=0.5, ML_bulge=0.7):
    """
    Compute dark matter fraction for each galaxy.

    V_bar^2 = Vgas^2 + ML_disk * Vdisk^2 + ML_bulge * Vbul^2
    V_DM^2 = Vobs^2 - V_bar^2  (where positive)
    DM_fraction = sum(V_DM^2) / sum(Vobs^2) over all radii

    ML_disk = 0.5 is the standard SPARC assumption for 3.6um.
    """
    dm_fractions = {}

    for name, data in curves.items():
        Vobs = data['Vobs']
        Vgas = data['Vgas']
        Vdisk = data['Vdisk']
        Vbul = data['Vbul']

        # Baryonic velocity (quadrature sum)
        Vbar2 = Vgas**2 + ML_disk * Vdisk**2 + ML_bulge * Vbul**2

        # Observed velocity squared
        Vobs2 = Vobs**2

        # Dark matter contribution (must be non-negative)
        VDM2 = np.maximum(Vobs2 - Vbar2, 0)

        # DM fraction: fraction of total kinetic support from "dark matter"
        total_Vobs2 = np.sum(Vobs2)
        total_VDM2 = np.sum(VDM2)

        if total_Vobs2 > 0:
            dm_frac = total_VDM2 / total_Vobs2
        else:
            dm_frac = 0

        # Also compute at outermost radius (where DM dominates most)
        if len(Vobs) > 3:
            outer_idx = slice(-3, None)  # last 3 points
            outer_Vobs2 = np.mean(Vobs[outer_idx]**2)
            outer_Vbar2 = np.mean(Vbar2[outer_idx])
            outer_VDM2 = max(outer_Vobs2 - outer_Vbar2, 0)
            outer_dm_frac = outer_VDM2 / outer_Vobs2 if outer_Vobs2 > 0 else 0
        else:
            outer_dm_frac = dm_frac

        dm_fractions[name] = {
            'dm_fraction_total': float(dm_frac),
            'dm_fraction_outer': float(outer_dm_frac),
            'n_points': len(Vobs),
            'max_radius_kpc': float(data['R'][-1]) if len(data['R']) > 0 else 0,
        }

    return dm_fractions


# ============================================================
# CROSS-MATCHING
# ============================================================

def cross_match(sparc_galaxies, pitch_angles):
    """
    Match SPARC galaxies with pitch angle catalog.
    Try exact match first, then common name variants.
    """
    matches = {}
    unmatched_sparc = []

    # Build lookup with normalized names
    pitch_lookup = {}
    for name in pitch_angles:
        # Normalize: uppercase, remove spaces, hyphens
        norm = name.upper().replace(' ', '').replace('-', '')
        pitch_lookup[norm] = name
        # Also try without leading zeros in numbers
        import re
        compact = re.sub(r'0+(\d)', r'\1', norm)
        pitch_lookup[compact] = name

    for sparc_name in sparc_galaxies:
        norm_sparc = sparc_name.upper().replace(' ', '').replace('-', '')

        matched = False
        # Try direct match
        if norm_sparc in pitch_lookup:
            pitch_name = pitch_lookup[norm_sparc]
            matches[sparc_name] = pitch_name
            matched = True
        else:
            # Try partial matches
            for norm_pitch, pitch_name in pitch_lookup.items():
                if norm_sparc in norm_pitch or norm_pitch in norm_sparc:
                    matches[sparc_name] = pitch_name
                    matched = True
                    break

        if not matched:
            unmatched_sparc.append(sparc_name)

    return matches, unmatched_sparc


# ============================================================
# CORRELATION ANALYSIS
# ============================================================

def correlation_analysis(matched_data):
    """
    Run all four pre-registered correlation tests.
    Returns results dict with statistics.
    """
    from scipy import stats

    results = {}

    # Extract arrays
    dm_frac = np.array([d['dm_fraction_total'] for d in matched_data])
    dm_frac_outer = np.array([d['dm_fraction_outer'] for d in matched_data])
    pitch = np.array([d['pitch_mean_deg'] for d in matched_data])
    hubble_T = np.array([d['hubble_type'] for d in matched_data])
    luminosity = np.array([d['luminosity_1e9Lsun'] for d in matched_data])
    Vflat = np.array([d['Vflat_kms'] for d in matched_data])
    spiral_type = np.array([d['spiral_type'] for d in matched_data])

    # Filter valid data
    valid = (pitch > 0) & (dm_frac > 0) & (dm_frac < 1) & (Vflat > 0)
    n_valid = np.sum(valid)

    print(f"\n  Valid matched galaxies for analysis: {n_valid}")

    if n_valid < 10:
        print("  WARNING: Too few matches for reliable statistics")
        results['n_valid'] = int(n_valid)
        results['warning'] = 'Too few matches'
        return results

    dm_v = dm_frac[valid]
    dm_v_outer = dm_frac_outer[valid]
    pitch_v = pitch[valid]
    hubble_v = hubble_T[valid]
    lum_v = luminosity[valid]
    vflat_v = Vflat[valid]
    spiral_v = spiral_type[valid]

    # ── P1: PITCH ANGLE vs DM FRACTION ──
    print("\n  ── P1: PITCH ANGLE vs DM FRACTION ──")
    r_pearson, p_pearson = stats.pearsonr(pitch_v, dm_v)
    r_spearman, p_spearman = stats.spearmanr(pitch_v, dm_v)
    print(f"    Pearson:  r = {r_pearson:+.4f}, p = {p_pearson:.4e}")
    print(f"    Spearman: r = {r_spearman:+.4f}, p = {p_spearman:.4e}")

    # Same for outer DM fraction
    r_outer, p_outer = stats.spearmanr(pitch_v, dm_v_outer)
    print(f"    Spearman (outer): r = {r_outer:+.4f}, p = {p_outer:.4e}")

    sig_p1 = p_spearman < 0.05
    print(f"    Significant (p<0.05): {'YES' if sig_p1 else 'NO'}")
    if r_spearman < 0:
        print(f"    Direction: NEGATIVE (tighter spiral → more DM) — MATCHES PREDICTION P1")
    else:
        print(f"    Direction: POSITIVE (tighter spiral → less DM) — CONTRADICTS PREDICTION P1")

    results['P1_pitch_vs_dm'] = {
        'pearson_r': float(r_pearson), 'pearson_p': float(p_pearson),
        'spearman_r': float(r_spearman), 'spearman_p': float(p_spearman),
        'spearman_outer_r': float(r_outer), 'spearman_outer_p': float(p_outer),
        'significant': bool(sig_p1),
        'direction': 'negative' if r_spearman < 0 else 'positive',
        'matches_prediction': bool(r_spearman < 0),
    }

    # ── P2: SPIRAL TYPE vs DM FRACTION ──
    print("\n  ── P2: SPIRAL TYPE vs DM FRACTION ──")
    for stype in ['G', 'M', 'F']:
        mask = spiral_v == stype
        n_type = np.sum(mask)
        if n_type > 2:
            mean_dm = np.mean(dm_v[mask])
            std_dm = np.std(dm_v[mask])
            label = {'G': 'Grand-design', 'M': 'Multi-armed', 'F': 'Flocculent'}[stype]
            print(f"    {label} ({stype}): n={n_type}, DM_frac={mean_dm:.3f} ± {std_dm:.3f}")

    # ANOVA across spiral types
    groups = [dm_v[spiral_v == t] for t in ['G', 'M', 'F'] if np.sum(spiral_v == t) > 2]
    if len(groups) >= 2:
        F_stat, p_anova = stats.f_oneway(*groups)
        print(f"    ANOVA: F = {F_stat:.3f}, p = {p_anova:.4e}")
        sig_p2 = p_anova < 0.05
        print(f"    Significant (p<0.05): {'YES' if sig_p2 else 'NO'}")
        results['P2_spiral_type_vs_dm'] = {
            'F_stat': float(F_stat), 'p_value': float(p_anova),
            'significant': bool(sig_p2),
        }
    else:
        results['P2_spiral_type_vs_dm'] = {'note': 'Insufficient groups'}

    # ── P3: SURFACE BRIGHTNESS (symmetry proxy) vs DM FRACTION ──
    print("\n  ── P3: SURFACE BRIGHTNESS vs DM FRACTION ──")
    # Higher surface brightness = more concentrated = more symmetric
    SB_valid = np.array([d['eff_SB'] for d in matched_data])[valid]
    if np.sum(SB_valid > 0) > 10:
        sb_mask = SB_valid > 0
        r_sb, p_sb = stats.spearmanr(SB_valid[sb_mask], dm_v[sb_mask])
        print(f"    Spearman: r = {r_sb:+.4f}, p = {p_sb:.4e}")
        sig_p3 = p_sb < 0.05
        print(f"    Significant (p<0.05): {'YES' if sig_p3 else 'NO'}")
        results['P3_SB_vs_dm'] = {
            'spearman_r': float(r_sb), 'spearman_p': float(p_sb),
            'significant': bool(sig_p3),
        }

    # ── P4: HUBBLE TYPE vs DM FRACTION ──
    print("\n  ── P4: HUBBLE TYPE vs DM FRACTION ──")
    r_hubble, p_hubble = stats.spearmanr(hubble_v, dm_v)
    print(f"    Spearman: r = {r_hubble:+.4f}, p = {p_hubble:.4e}")
    sig_p4 = p_hubble < 0.05
    print(f"    Significant (p<0.05): {'YES' if sig_p4 else 'NO'}")

    # Group by Hubble type
    for T in sorted(set(hubble_v)):
        mask = hubble_v == T
        n_T = np.sum(mask)
        if n_T > 1:
            T_labels = {0:'S0', 1:'Sa', 2:'Sab', 3:'Sb', 4:'Sbc',
                        5:'Sc', 6:'Scd', 7:'Sd', 8:'Sdm', 9:'Sm', 10:'Im', 11:'BCD'}
            label = T_labels.get(int(T), f'T={int(T)}')
            print(f"    {label:>4}: n={n_T:2d}, DM_frac={np.mean(dm_v[mask]):.3f} ± {np.std(dm_v[mask]):.3f}")

    results['P4_hubble_vs_dm'] = {
        'spearman_r': float(r_hubble), 'spearman_p': float(p_hubble),
        'significant': bool(sig_p4),
    }

    # ── CONTROLS: Mass and Vflat vs DM fraction ──
    print("\n  ── CONTROLS ──")

    # Vflat (proxy for total mass)
    r_vflat, p_vflat = stats.spearmanr(vflat_v, dm_v)
    print(f"    Vflat vs DM: r = {r_vflat:+.4f}, p = {p_vflat:.4e}")

    # Luminosity (proxy for stellar mass)
    lum_mask = lum_v > 0
    if np.sum(lum_mask) > 10:
        r_lum, p_lum = stats.spearmanr(np.log10(lum_v[lum_mask]), dm_v[lum_mask])
        print(f"    log(L) vs DM: r = {r_lum:+.4f}, p = {p_lum:.4e}")
    else:
        r_lum, p_lum = 0, 1

    results['controls'] = {
        'Vflat_r': float(r_vflat), 'Vflat_p': float(p_vflat),
        'logL_r': float(r_lum), 'logL_p': float(p_lum),
    }

    # ── KEY COMPARISON ──
    print("\n  ── KEY COMPARISON: GEOMETRY vs MASS ──")
    print(f"    Pitch angle vs DM:  |r| = {abs(r_spearman):.4f}")
    print(f"    Vflat vs DM:        |r| = {abs(r_vflat):.4f}")
    print(f"    log(L) vs DM:       |r| = {abs(r_lum):.4f}")

    if abs(r_spearman) > abs(r_vflat) and abs(r_spearman) > abs(r_lum):
        print(f"    >>> GEOMETRY correlates STRONGER than mass <<<")
        results['geometry_beats_mass'] = True
    else:
        print(f"    Mass correlates stronger than geometry")
        results['geometry_beats_mass'] = False

    results['n_valid'] = int(n_valid)
    return results


# ============================================================
# MAIN
# ============================================================

def main():
    data_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')

    print("=" * 70)
    print("  DARK MATTER GEOMETRIC CORRELATION TEST")
    print("  Cipher cross-scale validation — {2,3} at galactic scale")
    print("=" * 70)

    # Parse data
    print("\n  Loading SPARC galaxy table...")
    galaxies = parse_sparc_galaxies(os.path.join(data_dir, 'SPARC_Lelli2016c.mrt'))
    print(f"    {len(galaxies)} galaxies loaded")

    print("  Loading SPARC rotation curves...")
    curves = parse_sparc_rotation_curves(os.path.join(data_dir, 'MassModels_Lelli2016c.mrt'))
    print(f"    {len(curves)} rotation curves loaded")

    print("  Loading Diaz-Garcia 2019 pitch angles...")
    pitches = parse_pitch_angles(os.path.join(data_dir, 'diaz_garcia_pitch_angles.dat'))
    print(f"    {len(pitches)} pitch angles loaded")

    # Cross-match
    print("\n  Cross-matching SPARC × pitch angle catalog...")
    matches, unmatched = cross_match(galaxies, pitches)
    print(f"    Matched: {len(matches)} galaxies")
    print(f"    Unmatched SPARC: {len(unmatched)} galaxies")

    if len(matches) < 5:
        print("\n  ERROR: Too few matches. Checking name formats...")
        print(f"  Sample SPARC names: {list(galaxies.keys())[:5]}")
        print(f"  Sample pitch names: {list(pitches.keys())[:5]}")
        return

    # Compute DM fractions
    print("\n  Computing dark matter fractions...")
    dm_fracs = compute_dm_fraction(curves)
    print(f"    {len(dm_fracs)} galaxies with DM fraction")

    # Build matched dataset
    matched_data = []
    for sparc_name, pitch_name in matches.items():
        if sparc_name in galaxies and sparc_name in dm_fracs and pitch_name in pitches:
            entry = {}
            entry.update(galaxies[sparc_name])
            entry.update(dm_fracs[sparc_name])
            entry.update(pitches[pitch_name])
            matched_data.append(entry)

    print(f"  Complete matched records: {len(matched_data)}")

    # Run analysis
    print("\n" + "=" * 70)
    print("  CORRELATION ANALYSIS")
    print("=" * 70)

    results = correlation_analysis(matched_data)

    # Summary
    print("\n" + "=" * 70)
    print("  SUMMARY")
    print("=" * 70)

    n_sig = sum(1 for k, v in results.items()
                if isinstance(v, dict) and v.get('significant', False))
    total_tests = sum(1 for k, v in results.items()
                      if isinstance(v, dict) and 'significant' in v)

    print(f"\n  Tests run: {total_tests}")
    print(f"  Significant (p<0.05): {n_sig}")
    print(f"  Geometry beats mass: {results.get('geometry_beats_mass', 'N/A')}")

    if n_sig == 0:
        print("\n  RESULT: No significant geometric correlations found.")
        print("  The dark matter fraction does not depend on galaxy geometry")
        print("  in this sample. The geometric amplification hypothesis is")
        print("  NOT SUPPORTED by this data.")
    elif n_sig >= 1 and not results.get('geometry_beats_mass'):
        print("\n  RESULT: Some geometric correlations exist, but mass")
        print("  correlates stronger. Geometry may be secondary to mass.")
        print("  The hypothesis is PARTIALLY SUPPORTED but not conclusive.")
    elif results.get('geometry_beats_mass'):
        print("\n  RESULT: Geometry correlates STRONGER than mass with")
        print("  the dark matter fraction. This is consistent with the")
        print("  geometric amplification hypothesis.")
        print("  The hypothesis is SUPPORTED by this data.")

    # Save results
    timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
    output = {
        'test': 'dark_matter_geometric_correlation',
        'timestamp': timestamp,
        'n_sparc': len(galaxies),
        'n_pitch': len(pitches),
        'n_matched': len(matches),
        'n_complete': len(matched_data),
        'ML_disk': 0.5,
        'ML_bulge': 0.7,
        'results': results,
        'matched_galaxies': [d['name'] for d in matched_data],
    }

    outfile = f'dm_correlation_results_{timestamp}.json'
    with open(outfile, 'w') as f:
        json.dump(output, f, indent=2, default=str)
    print(f"\n  Results saved: {outfile}")

    # Try to plot
    try:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt

        fig, axes = plt.subplots(2, 2, figsize=(14, 10))

        dm_arr = np.array([d['dm_fraction_total'] for d in matched_data])
        pitch_arr = np.array([d['pitch_mean_deg'] for d in matched_data])
        hubble_arr = np.array([d['hubble_type'] for d in matched_data])
        vflat_arr = np.array([d['Vflat_kms'] for d in matched_data])
        valid = (pitch_arr > 0) & (dm_arr > 0) & (dm_arr < 1) & (vflat_arr > 0)

        # P1: Pitch angle vs DM fraction
        ax = axes[0, 0]
        ax.scatter(pitch_arr[valid], dm_arr[valid], alpha=0.6, s=20)
        ax.set_xlabel('Pitch Angle (deg)')
        ax.set_ylabel('Dark Matter Fraction')
        ax.set_title('P1: Pitch Angle vs DM Fraction')
        ax.grid(True, alpha=0.3)

        # P2: Spiral type
        ax = axes[0, 1]
        spiral_arr = np.array([d['spiral_type'] for d in matched_data])
        for i, stype in enumerate(['G', 'M', 'F']):
            mask = valid & (spiral_arr == stype)
            if np.sum(mask) > 0:
                label = {'G': 'Grand-design', 'M': 'Multi-armed', 'F': 'Flocculent'}[stype]
                ax.scatter(pitch_arr[mask], dm_arr[mask], alpha=0.6, s=20, label=f'{label} (n={np.sum(mask)})')
        ax.legend()
        ax.set_xlabel('Pitch Angle (deg)')
        ax.set_ylabel('Dark Matter Fraction')
        ax.set_title('P2: Spiral Type Differentiation')
        ax.grid(True, alpha=0.3)

        # P4: Hubble type
        ax = axes[1, 0]
        ax.scatter(hubble_arr[valid], dm_arr[valid], alpha=0.6, s=20)
        ax.set_xlabel('Hubble Type T')
        ax.set_ylabel('Dark Matter Fraction')
        ax.set_title('P4: Hubble Type vs DM Fraction')
        ax.grid(True, alpha=0.3)

        # Control: Vflat
        ax = axes[1, 1]
        ax.scatter(vflat_arr[valid], dm_arr[valid], alpha=0.6, s=20, color='gray')
        ax.set_xlabel('Vflat (km/s)')
        ax.set_ylabel('Dark Matter Fraction')
        ax.set_title('Control: Mass (Vflat) vs DM Fraction')
        ax.grid(True, alpha=0.3)

        plt.suptitle('Dark Matter Geometric Correlation — Cipher Cross-Scale Test', fontsize=14)
        plt.tight_layout()
        plt.savefig(f'dm_correlation_plot_{timestamp}.png', dpi=150)
        print(f"  Plot saved: dm_correlation_plot_{timestamp}.png")

    except ImportError:
        print("  matplotlib not available — skipping plots")


if __name__ == "__main__":
    main()