"""
TLT-003: Progressive Compaction Test
=====================================
Tests whether the N=3 triangular interference pattern naturally develops
site differentiation when subjected to multi-cycle decoherence accumulation.

Hypothesis (pre-registered):
  H0: All maxima remain equivalent (CV ~ 0, k=1) at all parameter values
  H1: Site differentiation emerges (CV > 0, k >= 2) at some parameter values

Four variants tested with identical grids and measurement:
  A: Single burst (control)
  B: Perfect re-phasing (M identical bursts)
  C: Decoherence phase drift (theory-faithful, sweep t/T and M)
  D: Random phase reset (maximal decoherence, sweep M)

Variables: f (frequency), t (decoherence), c (speed of light).
No tuned parameters. The data determines the outcome.
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from pathlib import Path
from scipy import ndimage
from scipy.signal import find_peaks
import json
import time as timer

# ============================================================================
# CONSTANTS (identical to TLT-002 for continuity)
# ============================================================================
C_LIGHT = 299_792_458.0          # m/s (exact, SI definition)
H_PLANCK = 6.62607015e-34        # J·s (exact, 2019 SI)

OUTPUT_DIR = Path("/root/rhetroiluso/project_prometheus/time_ledger_theory/"
                  "tlt results/unaudited/TLT-003_compaction")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)


# ============================================================================
# CORE: N=3 INTERFERENCE WITH PHASE OFFSETS
# ============================================================================

def generate_n3_with_phases(wavelength, grid_extent, resolution,
                            phase_offsets=(0.0, 0.0, 0.0)):
    """
    Generate N=3 interference pattern with per-wave phase offsets.

    psi(x,y) = sum_{n=0}^{2} exp(i * (k*(x*cos(theta_n) + y*sin(theta_n)) + phi_n))

    Parameters:
        wavelength: wavelength (normalized units)
        grid_extent: half-width of domain
        resolution: grid points per axis
        phase_offsets: tuple of 3 phase values (radians), one per wave

    Returns:
        X, Y: coordinate meshgrids
        intensity: |psi|^2
        psi: complex field
    """
    x = np.linspace(-grid_extent, grid_extent, resolution)
    y = np.linspace(-grid_extent, grid_extent, resolution)
    X, Y = np.meshgrid(x, y)

    k = 2.0 * np.pi / wavelength
    psi = np.zeros_like(X, dtype=complex)

    for n in range(3):
        theta = 2.0 * np.pi * n / 3.0
        kx = k * np.cos(theta)
        ky = k * np.sin(theta)
        psi += np.exp(1j * (kx * X + ky * Y + phase_offsets[n]))

    intensity = np.abs(psi) ** 2
    return X, Y, intensity, psi


# ============================================================================
# MEASUREMENT: SITE CLASSIFICATION
# ============================================================================

def find_maxima_intensities(intensity, threshold_fraction=0.5):
    """
    Find all local maxima in 2D intensity pattern and return their values.

    Uses morphological approach: label connected regions above threshold,
    find the peak value within each region.

    Returns:
        maxima_values: array of intensity values at each maximum
        maxima_positions: list of (row, col) positions
    """
    I_max = np.max(intensity)
    I_min = np.min(intensity)
    threshold = I_min + threshold_fraction * (I_max - I_min)

    # Binary mask of high-intensity regions
    binary = intensity > threshold

    # Label connected regions
    labeled, n_features = ndimage.label(binary)

    maxima_values = []
    maxima_positions = []

    for i in range(1, n_features + 1):
        region = (labeled == i)
        region_intensities = intensity[region]
        peak_val = np.max(region_intensities)

        # Find the position of the peak within this region
        peak_mask = region & (intensity == peak_val)
        positions = np.argwhere(peak_mask)
        if len(positions) > 0:
            pos = positions[0]  # take first if multiple pixels at same value
            maxima_values.append(peak_val)
            maxima_positions.append((int(pos[0]), int(pos[1])))

    return np.array(maxima_values), maxima_positions


def classify_sites(maxima_values):
    """
    Classify intensity maxima into clusters using gap-based splitting.

    Returns:
        n_classes: number of distinct intensity classes found
        cv: coefficient of variation of maxima intensities
        class_labels: array assigning each maximum to a class (0-indexed)
        class_means: mean intensity of each class
        silhouette: silhouette-like score (inter-class / intra-class spread)
    """
    if len(maxima_values) < 2:
        return 1, 0.0, np.zeros(len(maxima_values), dtype=int), [np.mean(maxima_values)], 0.0

    cv = np.std(maxima_values) / np.mean(maxima_values) if np.mean(maxima_values) > 0 else 0.0

    # Sort values and look for natural gaps
    sorted_vals = np.sort(maxima_values)
    diffs = np.diff(sorted_vals)

    if len(diffs) == 0:
        return 1, cv, np.zeros(len(maxima_values), dtype=int), [np.mean(maxima_values)], 0.0

    # Significant gap = any gap > 3x the median gap
    median_gap = np.median(diffs)
    if median_gap == 0:
        # All values identical or near-identical
        # Use absolute threshold instead
        gap_threshold = 0.01 * np.mean(maxima_values)
    else:
        gap_threshold = 3.0 * median_gap

    gap_indices = np.where(diffs > gap_threshold)[0]

    if len(gap_indices) == 0:
        # No significant gaps — all one class
        return 1, cv, np.zeros(len(maxima_values), dtype=int), [np.mean(maxima_values)], 0.0

    # Split at gaps to form classes
    boundaries = [0] + list(gap_indices + 1) + [len(sorted_vals)]
    n_classes = len(boundaries) - 1

    # Map back to original indices
    sort_order = np.argsort(maxima_values)
    class_labels = np.zeros(len(maxima_values), dtype=int)
    class_means = []

    for c in range(n_classes):
        start, end = boundaries[c], boundaries[c + 1]
        original_indices = sort_order[start:end]
        class_labels[original_indices] = c
        class_means.append(np.mean(sorted_vals[start:end]))

    # Compute silhouette-like score
    if n_classes >= 2:
        # Inter-class distance (between class means)
        inter = np.max(class_means) - np.min(class_means)
        # Intra-class spread (average within-class std)
        intra_stds = []
        for c in range(n_classes):
            vals = maxima_values[class_labels == c]
            if len(vals) > 1:
                intra_stds.append(np.std(vals))
            else:
                intra_stds.append(0.0)
        intra = np.mean(intra_stds) if intra_stds else 0.0
        silhouette = inter / (inter + intra) if (inter + intra) > 0 else 0.0
    else:
        silhouette = 0.0

    return n_classes, cv, class_labels, class_means, silhouette


# ============================================================================
# VARIANT A: CONTROL (single burst)
# ============================================================================

def run_variant_a(wavelength, grid_extent, resolution):
    """Single coherent N=3 burst. Baseline measurement."""
    print("\n" + "=" * 70)
    print("VARIANT A: CONTROL — Single N=3 burst")
    print("=" * 70)

    X, Y, intensity, _ = generate_n3_with_phases(
        wavelength, grid_extent, resolution, phase_offsets=(0, 0, 0)
    )

    maxima_vals, maxima_pos = find_maxima_intensities(intensity)
    n_classes, cv, labels, means, sil = classify_sites(maxima_vals)

    print(f"  Maxima found: {len(maxima_vals)}")
    print(f"  Intensity range of maxima: [{np.min(maxima_vals):.4f}, {np.max(maxima_vals):.4f}]")
    print(f"  CV = {cv:.6f}")
    print(f"  Classes found: {n_classes}")
    print(f"  Class means: {[f'{m:.4f}' for m in means]}")
    print(f"  Silhouette: {sil:.4f}")

    # Save plot
    fig, ax = plt.subplots(figsize=(10, 10))
    im = ax.pcolormesh(X, Y, intensity, cmap='inferno', shading='auto')
    fig.colorbar(im, ax=ax, shrink=0.8, label='Intensity |ψ|²')
    ax.set_title('Variant A: Single N=3 Burst (Control)', fontsize=16)
    ax.set_xlabel('x / λ', fontsize=14)
    ax.set_ylabel('y / λ', fontsize=14)
    ax.set_aspect('equal')
    plt.tight_layout()
    fig.savefig(OUTPUT_DIR / 'variant_A_control.png', dpi=150, bbox_inches='tight')
    plt.close(fig)
    print(f"  Saved: variant_A_control.png")

    return {
        'variant': 'A',
        'M': 1,
        'n_maxima': len(maxima_vals),
        'cv': float(cv),
        'n_classes': n_classes,
        'class_means': [float(m) for m in means],
        'silhouette': float(sil),
        'I_min_maxima': float(np.min(maxima_vals)),
        'I_max_maxima': float(np.max(maxima_vals))
    }


# ============================================================================
# VARIANT B: PERFECT RE-PHASING
# ============================================================================

def run_variant_b(wavelength, grid_extent, resolution, M_values):
    """M identical bursts accumulated. Tests pure accumulation effect."""
    print("\n" + "=" * 70)
    print("VARIANT B: PERFECT RE-PHASING — M identical bursts")
    print(f"Sweep: M = {M_values}")
    print("=" * 70)

    results = []

    for M in M_values:
        # M identical bursts = M * single burst intensity
        # Mathematically: accumulated = M * I_single
        # CV should be exactly 0 (all sites scale equally)
        # But we compute it explicitly to verify
        X, Y, intensity_single, _ = generate_n3_with_phases(
            wavelength, grid_extent, resolution, phase_offsets=(0, 0, 0)
        )

        # Accumulate M identical copies
        intensity_acc = M * intensity_single

        maxima_vals, maxima_pos = find_maxima_intensities(intensity_acc)
        n_classes, cv, labels, means, sil = classify_sites(maxima_vals)

        print(f"\n  M={M}: maxima={len(maxima_vals)}, CV={cv:.6f}, "
              f"classes={n_classes}, means={[f'{m:.2f}' for m in means]}")

        results.append({
            'variant': 'B',
            'M': M,
            'n_maxima': len(maxima_vals),
            'cv': float(cv),
            'n_classes': n_classes,
            'class_means': [float(m) for m in means],
            'silhouette': float(sil),
            'I_min_maxima': float(np.min(maxima_vals)),
            'I_max_maxima': float(np.max(maxima_vals))
        })

    return results


# ============================================================================
# VARIANT C: DECOHERENCE PHASE DRIFT (theory-faithful)
# ============================================================================

def run_variant_c(wavelength, grid_extent, resolution, M_values, tT_values):
    """
    M bursts with inter-burst phase drift governed by decoherence parameter.

    At each decoherence boundary, each wave component accumulates a phase
    offset. The drift per cycle for wave n is:

        delta_phi_n(cycle) = 2*pi * (n+1) * (t/T) * cycle_index / M

    where (t/T) is the decoherence-to-period ratio. This models partial
    phase memory loss: t/T=0 means perfect memory (= Variant B), t/T=1
    means one full period of drift per M cycles.

    The key: each wave component drifts at a DIFFERENT rate (proportional
    to its index n+1). This is the minimal assumption — the three directional
    components lose phase coherence at rates proportional to their directional
    index. If all three drifted equally, the pattern would simply translate
    (no site differentiation).
    """
    print("\n" + "=" * 70)
    print("VARIANT C: DECOHERENCE PHASE DRIFT (theory-faithful)")
    print(f"Sweep: M = {M_values}")
    print(f"Sweep: t/T = {tT_values}")
    print("=" * 70)

    results = []
    all_cv_data = np.zeros((len(M_values), len(tT_values)))

    for i_M, M in enumerate(M_values):
        for i_tT, tT in enumerate(tT_values):
            # Accumulate M bursts with progressive phase drift
            X, Y = None, None
            intensity_acc = None

            for cycle in range(M):
                # Phase drift for each wave component
                phase0 = 2.0 * np.pi * 1 * tT * cycle / M
                phase1 = 2.0 * np.pi * 2 * tT * cycle / M
                phase2 = 2.0 * np.pi * 3 * tT * cycle / M
                phases = (phase0, phase1, phase2)

                X_c, Y_c, intensity_c, _ = generate_n3_with_phases(
                    wavelength, grid_extent, resolution, phase_offsets=phases
                )

                if intensity_acc is None:
                    X, Y = X_c, Y_c
                    intensity_acc = intensity_c.copy()
                else:
                    intensity_acc += intensity_c

            # Normalize by M
            intensity_acc /= M

            maxima_vals, maxima_pos = find_maxima_intensities(intensity_acc)
            n_classes, cv, labels, means, sil = classify_sites(maxima_vals)

            all_cv_data[i_M, i_tT] = cv

            print(f"  M={M:4d}, t/T={tT:.2f}: maxima={len(maxima_vals):3d}, "
                  f"CV={cv:.6f}, classes={n_classes}, "
                  f"means={[f'{m:.2f}' for m in means]}")

            results.append({
                'variant': 'C',
                'M': M,
                'tT': float(tT),
                'n_maxima': len(maxima_vals),
                'cv': float(cv),
                'n_classes': n_classes,
                'class_means': [float(m) for m in means],
                'silhouette': float(sil),
                'I_min_maxima': float(np.min(maxima_vals)) if len(maxima_vals) > 0 else 0,
                'I_max_maxima': float(np.max(maxima_vals)) if len(maxima_vals) > 0 else 0
            })

            # Save pattern image for cases with site differentiation
            if n_classes >= 2:
                fig, ax = plt.subplots(figsize=(10, 10))
                im = ax.pcolormesh(X, Y, intensity_acc, cmap='inferno', shading='auto')
                fig.colorbar(im, ax=ax, shrink=0.8, label='Accumulated Intensity')
                ax.set_title(f'Variant C: M={M}, t/T={tT:.2f}\n'
                             f'{n_classes} classes, CV={cv:.4f}', fontsize=16)
                ax.set_xlabel('x / λ', fontsize=14)
                ax.set_ylabel('y / λ', fontsize=14)
                ax.set_aspect('equal')

                # Mark maxima colored by class
                for idx, (row, col) in enumerate(maxima_pos):
                    color = ['cyan', 'magenta', 'lime'][labels[idx] % 3]
                    ax.plot(X[row, col], Y[row, col], 'o', color=color,
                            markersize=4, markeredgecolor='white',
                            markeredgewidth=0.5)

                plt.tight_layout()
                fname = f'variant_C_M{M}_tT{tT:.2f}.png'
                fig.savefig(OUTPUT_DIR / fname, dpi=150, bbox_inches='tight')
                plt.close(fig)
                print(f"    ** Differentiation detected — saved: {fname}")

    # Save CV sweep heatmap
    fig, ax = plt.subplots(figsize=(12, 8))
    im = ax.imshow(all_cv_data, aspect='auto', cmap='viridis', origin='lower',
                   extent=[tT_values[0], tT_values[-1], 0, len(M_values)])
    ax.set_yticks(np.arange(len(M_values)) + 0.5)
    ax.set_yticklabels([str(m) for m in M_values])
    ax.set_xlabel('t/T (decoherence ratio)', fontsize=14)
    ax.set_ylabel('M (number of cycles)', fontsize=14)
    ax.set_title('Variant C: Coefficient of Variation vs Parameters\n'
                 'CV ~ 0 = all sites equivalent | CV > 0 = site differentiation',
                 fontsize=14)
    fig.colorbar(im, ax=ax, label='CV of maxima intensities')
    plt.tight_layout()
    fig.savefig(OUTPUT_DIR / 'variant_C_CV_sweep.png', dpi=150, bbox_inches='tight')
    plt.close(fig)
    print(f"\n  Saved: variant_C_CV_sweep.png")

    return results


# ============================================================================
# VARIANT D: RANDOM PHASE RESET
# ============================================================================

def run_variant_d(wavelength, grid_extent, resolution, M_values, n_seeds=10):
    """
    M bursts with independent random phase reset per wave at each boundary.
    Tests the maximal-decoherence limit.
    Multiple random seeds for statistical robustness.
    """
    print("\n" + "=" * 70)
    print("VARIANT D: RANDOM PHASE RESET (maximal decoherence)")
    print(f"Sweep: M = {M_values}")
    print(f"Random seeds per M: {n_seeds}")
    print("=" * 70)

    results = []

    for M in M_values:
        seed_cvs = []
        seed_classes = []

        for seed in range(n_seeds):
            rng = np.random.RandomState(seed * 1000 + M)

            X, Y = None, None
            intensity_acc = None

            for cycle in range(M):
                # Fully random phases for each wave component
                phases = tuple(rng.uniform(0, 2 * np.pi, size=3))

                X_c, Y_c, intensity_c, _ = generate_n3_with_phases(
                    wavelength, grid_extent, resolution, phase_offsets=phases
                )

                if intensity_acc is None:
                    X, Y = X_c, Y_c
                    intensity_acc = intensity_c.copy()
                else:
                    intensity_acc += intensity_c

            intensity_acc /= M

            maxima_vals, maxima_pos = find_maxima_intensities(intensity_acc)
            n_classes, cv, labels, means, sil = classify_sites(maxima_vals)

            seed_cvs.append(cv)
            seed_classes.append(n_classes)

        mean_cv = np.mean(seed_cvs)
        std_cv = np.std(seed_cvs)
        mean_classes = np.mean(seed_classes)

        print(f"  M={M:4d}: CV = {mean_cv:.6f} +/- {std_cv:.6f}, "
              f"mean classes = {mean_classes:.1f}, "
              f"class counts = {dict(zip(*np.unique(seed_classes, return_counts=True)))}")

        results.append({
            'variant': 'D',
            'M': M,
            'n_seeds': n_seeds,
            'cv_mean': float(mean_cv),
            'cv_std': float(std_cv),
            'class_counts': {int(k): int(v) for k, v in
                             zip(*np.unique(seed_classes, return_counts=True))},
            'all_cvs': [float(c) for c in seed_cvs]
        })

        # Save one representative pattern for each M
        # Use seed=0 for reproducibility
        rng = np.random.RandomState(M)
        intensity_acc = None
        for cycle in range(M):
            phases = tuple(rng.uniform(0, 2 * np.pi, size=3))
            X_c, Y_c, intensity_c, _ = generate_n3_with_phases(
                wavelength, grid_extent, resolution, phase_offsets=phases
            )
            if intensity_acc is None:
                X, Y = X_c, Y_c
                intensity_acc = intensity_c.copy()
            else:
                intensity_acc += intensity_c
        intensity_acc /= M

        fig, ax = plt.subplots(figsize=(10, 10))
        im = ax.pcolormesh(X, Y, intensity_acc, cmap='inferno', shading='auto')
        fig.colorbar(im, ax=ax, shrink=0.8, label='Accumulated Intensity')
        ax.set_title(f'Variant D: M={M}, Random Phases (seed=M)\n'
                     f'CV={mean_cv:.4f} +/- {std_cv:.4f}', fontsize=16)
        ax.set_xlabel('x / λ', fontsize=14)
        ax.set_ylabel('y / λ', fontsize=14)
        ax.set_aspect('equal')
        plt.tight_layout()
        fig.savefig(OUTPUT_DIR / f'variant_D_M{M}_representative.png',
                    dpi=150, bbox_inches='tight')
        plt.close(fig)
        print(f"  Saved: variant_D_M{M}_representative.png")

    return results


# ============================================================================
# OVERLAY COMPARISON (blind — shows both honeycomb and triangular)
# ============================================================================

def generate_lattice_overlay(X, Y, intensity, lattice_type, spacing, title,
                             filename):
    """
    Overlay known lattice positions onto the accumulated pattern.
    Shows BOTH honeycomb and triangular lattice simultaneously
    in different colors so neither is privileged.
    """
    fig, axes = plt.subplots(1, 3, figsize=(30, 10))

    for ax in axes:
        im = ax.pcolormesh(X, Y, intensity, cmap='inferno', shading='auto')
        ax.set_aspect('equal')
        ax.set_xlabel('x / λ', fontsize=12)
        ax.set_ylabel('y / λ', fontsize=12)

    # Panel 1: Pattern only (no overlay)
    axes[0].set_title(f'{title}\nNo overlay', fontsize=14)

    # Panel 2: Triangular lattice overlay
    a_tri = spacing
    # Triangular lattice: all vertices of equilateral triangles
    n_range = 20
    tri_x, tri_y = [], []
    for i in range(-n_range, n_range + 1):
        for j in range(-n_range, n_range + 1):
            px = i * a_tri + j * a_tri * 0.5
            py = j * a_tri * np.sqrt(3) / 2
            if abs(px) <= np.max(X) and abs(py) <= np.max(Y):
                tri_x.append(px)
                tri_y.append(py)
    axes[1].plot(tri_x, tri_y, 'o', color='cyan', markersize=5,
                 markeredgecolor='white', markeredgewidth=0.5, alpha=0.8)
    axes[1].set_title(f'{title}\nTriangular overlay (cyan)', fontsize=14)

    # Panel 3: Honeycomb lattice overlay
    # Honeycomb = two interpenetrating triangular sublattices offset by (a/sqrt(3), 0)
    a_hon = spacing
    hon_ax, hon_ay = [], []
    hon_bx, hon_by = [], []
    offset = a_hon / np.sqrt(3)
    for i in range(-n_range, n_range + 1):
        for j in range(-n_range, n_range + 1):
            # Sublattice A
            ax_pos = i * a_hon + j * a_hon * 0.5
            ay_pos = j * a_hon * np.sqrt(3) / 2
            if abs(ax_pos) <= np.max(X) and abs(ay_pos) <= np.max(Y):
                hon_ax.append(ax_pos)
                hon_ay.append(ay_pos)
            # Sublattice B (offset)
            bx_pos = ax_pos + offset * np.cos(np.pi / 6)
            by_pos = ay_pos + offset * np.sin(np.pi / 6)
            if abs(bx_pos) <= np.max(X) and abs(by_pos) <= np.max(Y):
                hon_bx.append(bx_pos)
                hon_by.append(by_pos)
    axes[2].plot(hon_ax, hon_ay, 'o', color='cyan', markersize=5,
                 markeredgecolor='white', markeredgewidth=0.5, alpha=0.8,
                 label='Sublattice A')
    axes[2].plot(hon_bx, hon_by, 's', color='magenta', markersize=5,
                 markeredgecolor='white', markeredgewidth=0.5, alpha=0.8,
                 label='Sublattice B')
    axes[2].legend(fontsize=12, loc='upper right')
    axes[2].set_title(f'{title}\nHoneycomb overlay (cyan A + magenta B)', fontsize=14)

    plt.tight_layout()
    fig.savefig(OUTPUT_DIR / filename, dpi=150, bbox_inches='tight')
    plt.close(fig)
    print(f"  Saved: {filename}")


# ============================================================================
# SUMMARY ANALYSIS
# ============================================================================

def generate_summary_plots(results_a, results_b, results_c, results_d):
    """Generate cross-variant comparison plots."""

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

    # Panel 1: Variant B — CV vs M
    M_b = [r['M'] for r in results_b]
    cv_b = [r['cv'] for r in results_b]
    axes[0, 0].plot(M_b, cv_b, 'o-', color='blue', markersize=8, linewidth=2)
    axes[0, 0].axhline(y=results_a['cv'], color='red', linestyle='--',
                        label=f'Variant A (control) CV={results_a["cv"]:.6f}')
    axes[0, 0].set_xlabel('M (number of cycles)', fontsize=12)
    axes[0, 0].set_ylabel('CV of maxima intensities', fontsize=12)
    axes[0, 0].set_title('Variant B: Perfect Re-phasing', fontsize=14)
    axes[0, 0].set_xscale('log')
    axes[0, 0].legend(fontsize=10)
    axes[0, 0].grid(True, alpha=0.3)

    # Panel 2: Variant C — CV vs t/T for different M
    M_vals_c = sorted(set(r['M'] for r in results_c))
    colors_c = plt.cm.viridis(np.linspace(0, 0.8, len(M_vals_c)))
    for idx, M in enumerate(M_vals_c):
        subset = [r for r in results_c if r['M'] == M]
        tT = [r['tT'] for r in subset]
        cv = [r['cv'] for r in subset]
        axes[0, 1].plot(tT, cv, 'o-', color=colors_c[idx], markersize=6,
                        linewidth=2, label=f'M={M}')
    axes[0, 1].axhline(y=results_a['cv'], color='red', linestyle='--',
                        label='Control', alpha=0.5)
    axes[0, 1].set_xlabel('t/T (decoherence ratio)', fontsize=12)
    axes[0, 1].set_ylabel('CV of maxima intensities', fontsize=12)
    axes[0, 1].set_title('Variant C: Decoherence Phase Drift', fontsize=14)
    axes[0, 1].legend(fontsize=9)
    axes[0, 1].grid(True, alpha=0.3)

    # Panel 3: Variant D — CV vs M (with error bars)
    M_d = [r['M'] for r in results_d]
    cv_d = [r['cv_mean'] for r in results_d]
    cv_err = [r['cv_std'] for r in results_d]
    axes[1, 0].errorbar(M_d, cv_d, yerr=cv_err, fmt='o-', color='green',
                        markersize=8, linewidth=2, capsize=5)
    axes[1, 0].axhline(y=results_a['cv'], color='red', linestyle='--',
                        label='Control', alpha=0.5)
    axes[1, 0].set_xlabel('M (number of cycles)', fontsize=12)
    axes[1, 0].set_ylabel('CV of maxima intensities', fontsize=12)
    axes[1, 0].set_title('Variant D: Random Phase Reset', fontsize=14)
    axes[1, 0].set_xscale('log')
    axes[1, 0].legend(fontsize=10)
    axes[1, 0].grid(True, alpha=0.3)

    # Panel 4: Number of classes across all variants
    # B
    classes_b = [r['n_classes'] for r in results_b]
    axes[1, 1].plot(M_b, classes_b, 's-', color='blue', markersize=8,
                    linewidth=2, label='B: Perfect')
    # D
    # Plot modal class count
    modal_d = []
    for r in results_d:
        counts = r['class_counts']
        modal_d.append(max(counts, key=counts.get))
    axes[1, 1].plot(M_d, modal_d, '^-', color='green', markersize=8,
                    linewidth=2, label='D: Random (modal)')
    # C at highest t/T
    max_tT = max(r['tT'] for r in results_c)
    c_at_max = [r for r in results_c if r['tT'] == max_tT]
    M_c_max = [r['M'] for r in c_at_max]
    classes_c_max = [r['n_classes'] for r in c_at_max]
    axes[1, 1].plot(M_c_max, classes_c_max, 'D-', color='orange', markersize=8,
                    linewidth=2, label=f'C: t/T={max_tT:.1f}')

    axes[1, 1].axhline(y=1, color='red', linestyle='--', alpha=0.3,
                        label='k=1 (no differentiation)')
    axes[1, 1].axhline(y=2, color='magenta', linestyle='--', alpha=0.3,
                        label='k=2 (honeycomb candidate)')
    axes[1, 1].set_xlabel('M (number of cycles)', fontsize=12)
    axes[1, 1].set_ylabel('Number of intensity classes', fontsize=12)
    axes[1, 1].set_title('Site Classification: All Variants', fontsize=14)
    axes[1, 1].set_xscale('log')
    axes[1, 1].legend(fontsize=9)
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].set_ylim(0.5, 4.5)
    axes[1, 1].set_yticks([1, 2, 3, 4])

    plt.suptitle('TLT-003: Progressive Compaction Test — Cross-Variant Summary',
                 fontsize=16, y=1.02)
    plt.tight_layout()
    fig.savefig(OUTPUT_DIR / 'TLT003_summary.png', dpi=150, bbox_inches='tight')
    plt.close(fig)
    print(f"\n  Saved: TLT003_summary.png")


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

def main():
    print("=" * 70)
    print("TLT-003: PROGRESSIVE COMPACTION TEST")
    print("=" * 70)
    print("Testing whether N=3 triangular pattern develops site")
    print("differentiation through multi-cycle decoherence accumulation.")
    print()
    print("H0: All maxima remain equivalent (CV ~ 0, k=1)")
    print("H1: Site differentiation emerges (CV > 0, k >= 2)")
    print()
    print("Variables: f, t, c — no tuned parameters.")
    print("=" * 70)

    start_time = timer.time()

    # Grid parameters (same scale as TLT-002 for continuity)
    wavelength = 1.0
    grid_extent = 5.0
    resolution = 800  # slightly lower than TLT-002's 1200 for speed across sweeps

    # N=3 peak spacing (from TLT-002): 2*lambda/3
    peak_spacing = 2.0 * wavelength / 3.0

    # ----------------------------
    # Run all variants
    # ----------------------------

    # Variant A: Control
    results_a = run_variant_a(wavelength, grid_extent, resolution)

    # Variant B: Perfect re-phasing
    M_values_b = [1, 5, 10, 50, 100, 500]
    results_b = run_variant_b(wavelength, grid_extent, resolution, M_values_b)

    # Variant C: Decoherence phase drift
    M_values_c = [10, 50, 100, 500]
    tT_values = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
    results_c = run_variant_c(wavelength, grid_extent, resolution,
                              M_values_c, tT_values)

    # Variant D: Random phase reset
    M_values_d = [10, 50, 100, 500]
    results_d = run_variant_d(wavelength, grid_extent, resolution,
                              M_values_d, n_seeds=10)

    # ----------------------------
    # Overlay comparisons for interesting cases
    # ----------------------------
    print("\n" + "=" * 70)
    print("OVERLAY COMPARISONS")
    print("=" * 70)

    # Always generate overlay for control (baseline)
    X, Y, intensity_a, _ = generate_n3_with_phases(
        wavelength, grid_extent, resolution, phase_offsets=(0, 0, 0)
    )
    generate_lattice_overlay(
        X, Y, intensity_a, 'both', peak_spacing,
        'Variant A: Control (single burst)',
        'overlay_variant_A.png'
    )

    # Generate overlay for Variant C at maximum differentiation
    c_differentiated = [r for r in results_c if r['n_classes'] >= 2]
    if c_differentiated:
        # Pick the case with highest CV
        best_c = max(c_differentiated, key=lambda r: r['cv'])
        M_best = best_c['M']
        tT_best = best_c['tT']
        print(f"\n  Best Variant C differentiation: M={M_best}, t/T={tT_best}, "
              f"CV={best_c['cv']:.6f}, classes={best_c['n_classes']}")

        # Regenerate that specific pattern
        intensity_best = None
        for cycle in range(M_best):
            phase0 = 2.0 * np.pi * 1 * tT_best * cycle / M_best
            phase1 = 2.0 * np.pi * 2 * tT_best * cycle / M_best
            phase2 = 2.0 * np.pi * 3 * tT_best * cycle / M_best
            _, _, ic, _ = generate_n3_with_phases(
                wavelength, grid_extent, resolution,
                phase_offsets=(phase0, phase1, phase2)
            )
            if intensity_best is None:
                intensity_best = ic.copy()
            else:
                intensity_best += ic
        intensity_best /= M_best

        generate_lattice_overlay(
            X, Y, intensity_best, 'both', peak_spacing,
            f'Variant C: M={M_best}, t/T={tT_best:.2f} (best differentiation)',
            f'overlay_variant_C_best.png'
        )
    else:
        print("\n  No Variant C cases showed differentiation (k >= 2)")

    # Generate overlay for Variant D representative
    rng = np.random.RandomState(100)
    M_d_best = M_values_d[-1]  # use highest M
    intensity_d = None
    for cycle in range(M_d_best):
        phases = tuple(rng.uniform(0, 2 * np.pi, size=3))
        _, _, ic, _ = generate_n3_with_phases(
            wavelength, grid_extent, resolution, phase_offsets=phases
        )
        if intensity_d is None:
            intensity_d = ic.copy()
        else:
            intensity_d += ic
    intensity_d /= M_d_best

    generate_lattice_overlay(
        X, Y, intensity_d, 'both', peak_spacing,
        f'Variant D: M={M_d_best}, Random Phases',
        f'overlay_variant_D_M{M_d_best}.png'
    )

    # ----------------------------
    # Summary plots
    # ----------------------------
    print("\n" + "=" * 70)
    print("GENERATING SUMMARY PLOTS")
    print("=" * 70)
    generate_summary_plots(results_a, results_b, results_c, results_d)

    # ----------------------------
    # Write data output
    # ----------------------------
    elapsed = timer.time() - start_time

    print("\n" + "=" * 70)
    print("WRITING DATA OUTPUT")
    print("=" * 70)

    data_file = OUTPUT_DIR / 'TLT003_compaction_data.txt'
    with open(data_file, 'w') as f:
        f.write("=" * 70 + "\n")
        f.write("TLT-003: PROGRESSIVE COMPACTION TEST — RAW DATA\n")
        f.write("=" * 70 + "\n")
        f.write(f"Date: 2026-03-13\n")
        f.write(f"Runtime: {elapsed:.1f} seconds\n")
        f.write(f"Grid: {resolution}x{resolution}, extent={grid_extent}, "
                f"wavelength={wavelength}\n")
        f.write(f"N=3 peak spacing: {peak_spacing:.4f} lambda\n\n")

        f.write("HYPOTHESIS:\n")
        f.write("  H0: All maxima equivalent (CV ~ 0, k=1) at all parameters\n")
        f.write("  H1: Site differentiation (CV > 0, k >= 2) at some parameters\n\n")

        # Variant A
        f.write("-" * 70 + "\n")
        f.write("VARIANT A: CONTROL\n")
        f.write("-" * 70 + "\n")
        f.write(json.dumps(results_a, indent=2) + "\n\n")

        # Variant B
        f.write("-" * 70 + "\n")
        f.write("VARIANT B: PERFECT RE-PHASING\n")
        f.write("-" * 70 + "\n")
        for r in results_b:
            f.write(f"  M={r['M']:4d}: CV={r['cv']:.6f}, classes={r['n_classes']}, "
                    f"means={r['class_means']}\n")
        f.write("\n")

        # Variant C
        f.write("-" * 70 + "\n")
        f.write("VARIANT C: DECOHERENCE PHASE DRIFT\n")
        f.write("-" * 70 + "\n")
        for r in results_c:
            f.write(f"  M={r['M']:4d}, t/T={r['tT']:.2f}: CV={r['cv']:.6f}, "
                    f"classes={r['n_classes']}, means={r['class_means']}\n")
        f.write("\n")

        # Variant D
        f.write("-" * 70 + "\n")
        f.write("VARIANT D: RANDOM PHASE RESET\n")
        f.write("-" * 70 + "\n")
        for r in results_d:
            f.write(f"  M={r['M']:4d}: CV={r['cv_mean']:.6f} +/- {r['cv_std']:.6f}, "
                    f"classes={r['class_counts']}\n")
        f.write("\n")

        # Verdict
        f.write("=" * 70 + "\n")
        f.write("VERDICT\n")
        f.write("=" * 70 + "\n")

        any_differentiation = False
        for r in results_c:
            if r['n_classes'] >= 2:
                any_differentiation = True
                break
        for r in results_d:
            if any(v >= 2 for v in r['class_counts'].keys()):
                any_differentiation = True
                break

        if any_differentiation:
            f.write("H1 SUPPORTED at some parameter values.\n")
            f.write("Site differentiation was observed. See images for spatial analysis.\n")
            f.write("Whether the differentiated pattern matches honeycomb topology\n")
            f.write("requires examination of the overlay comparison images.\n")
        else:
            f.write("H0 NOT REJECTED.\n")
            f.write("No site differentiation observed across all tested parameters.\n")
            f.write("The triangular pattern remains stable against decoherence cycling.\n")
            f.write("Progressive compaction does NOT occur from these mechanisms alone.\n")

        f.write("\nThis is an honest result, reported as-is.\n")
        f.write("=" * 70 + "\n")

    print(f"  Saved: {data_file}")

    # ----------------------------
    # Console verdict
    # ----------------------------
    print("\n" + "=" * 70)
    print("TLT-003 COMPLETE")
    print("=" * 70)
    print(f"Runtime: {elapsed:.1f} seconds")
    print(f"Output directory: {OUTPUT_DIR}")
    print(f"\nFiles generated:")
    for p in sorted(OUTPUT_DIR.iterdir()):
        print(f"  {p.name}")

    if any_differentiation:
        print("\n  >>> H1 SUPPORTED: Site differentiation detected <<<")
        print("  Examine overlay images for honeycomb topology comparison.")
    else:
        print("\n  >>> H0 NOT REJECTED: No differentiation detected <<<")
        print("  Triangular pattern persists across all parameter sweeps.")


if __name__ == '__main__':
    main()