run_user_tests.py

"""
Mock user-testing session for lfm-physics v0.x → v1.0 assessment.

5 users, 3 experiments each = 15 experiments total.
Users chosen to represent diverse use cases and skill levels.
Simulations represent what users WANT to do, not what the library was
designed for — this surfaces real-world API gaps.

Run:
    python run_user_tests.py
"""

from __future__ import annotations

import sys

import numpy as np

import lfm

PASS = "PASS"
FAIL = "FAIL"
GAP = "GAP"  # feature works but API is awkward / missing convenience

results: list[dict] = []


def record(user, exp, label, status, notes=""):
    tag = f"[{status}]"
    print(f"  {tag:6s} {label}")
    if notes:
        print(f"         {notes}")
    results.append({"user": user, "exp": exp, "label": label, "status": status, "notes": notes})


def section(title):
    print()
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)


# =========================================================================
# USER 1 — Alice, high-school student, curious about space
# Goal: visual intuition, "what does a black hole look like?"
# =========================================================================
section("USER 1: Alice (high-school) — black holes & space visuals")

try:
    # Exp 1a: Create a single soliton and inspect its chi-well
    cfg1 = lfm.SimulationConfig(grid_size=32, field_level=lfm.FieldLevel.REAL)
    sim1 = lfm.Simulation(cfg1)
    sim1.place_soliton((16, 16, 16), amplitude=8.0, sigma=2.5)
    sim1.equilibrate()
    sim1.run(steps=300)
    stats = lfm.chi_statistics(sim1.chi)
    hz = lfm.find_apparent_horizon(sim1.chi)
    assert stats["chi_min"] < 19.0
    assert hz["r_horizon"] >= 0
    record(
        "Alice",
        "1a",
        "Single BH chi-well + horizon",
        PASS,
        f"chi_min={stats['chi_min']:.2f}, r_horizon={hz['r_horizon']:.1f} cells",
    )
except Exception as e:
    record("Alice", "1a", "Single BH chi-well + horizon", FAIL, str(e))

try:
    # Exp 1b: Create two solitons on a collision course; use collider_event_display
    # with a manually constructed result dict (the collider display is a reporting
    # utility — detect_collision_events is the data source when trajectories are
    # available; here we verify the display function itself works)
    from lfm import collider_event_display

    fake_events = [
        {
            "step": 120,
            "type": "approach",
            "separation": 4.2,
            "energy_before": 12.5,
            "energy_after": 11.8,
        },
        {
            "step": 230,
            "type": "merge",
            "separation": 1.1,
            "energy_before": 11.8,
            "energy_after": 10.2,
        },
    ]
    result_dict = {
        "events": fake_events,
        "n_solitons": 2,
        "total_steps": 400,
        "score": 0.82,
    }
    display = collider_event_display(result_dict)
    assert isinstance(display, str)
    assert len(display.splitlines()) >= 3
    record(
        "Alice",
        "1b",
        "collider_event_display ASCII box",
        PASS,
        f"display has {len(display.splitlines())} lines",
    )
except Exception as e:
    record("Alice", "1b", "collider_event_display ASCII box", FAIL, str(e))

try:
    # Exp 1c: galaxy_summary_plot (visual — only runs if matplotlib present)
    try:
        import matplotlib  # noqa
        from lfm.viz import galaxy_summary_plot
        import matplotlib.pyplot as plt

        row = lfm.sparc_load("NGC6503")["NGC6503"]
        sim3 = lfm.Simulation(lfm.SimulationConfig(grid_size=32, field_level=lfm.FieldLevel.REAL))
        sim3.place_soliton((16, 16, 16), amplitude=5.0)
        sim3.equilibrate()
        fig, axes = galaxy_summary_plot(sim3, row)
        assert len(axes) == 2
        plt.close(fig)
        record(
            "Alice",
            "1c",
            "galaxy_summary_plot returns 2-panel figure",
            PASS,
            "matplotlib present; figure created and closed",
        )
    except ImportError:
        record(
            "Alice",
            "1c",
            "galaxy_summary_plot graceful ImportError",
            PASS,
            "matplotlib absent — expected ImportError raised",
        )
except Exception as e:
    record("Alice", "1c", "galaxy_summary_plot", FAIL, str(e))


# =========================================================================
# USER 2 — Bob, grad student in astrophysics
# Goal: reproduce galaxy rotation curves and compare to SPARC
# =========================================================================
section("USER 2: Bob (astrophysics grad student) — rotation curves")

try:
    # Exp 2a: Load SPARC data and inspect
    all_sparc = lfm.sparc_load()
    assert len(all_sparc) == 5
    for name, row in all_sparc.items():
        assert "r_kpc" in row and "v_obs_kms" in row
    record(
        "Bob",
        "2a",
        "sparc_load() returns 5 galaxies with correct keys",
        PASS,
        f"galaxies: {list(all_sparc.keys())}",
    )
except Exception as e:
    record("Bob", "2a", "sparc_load() returns 5 galaxies", FAIL, str(e))

try:
    # Exp 2b: Set up disk simulation and compute rotation curve
    sim4 = lfm.Simulation(
        lfm.SimulationConfig(grid_size=48, field_level=lfm.FieldLevel.REAL, chi0=19.0, kappa=1 / 63)
    )
    lfm.initialize_disk(sim4, n_solitons=30, r_inner=4.0, r_outer=18.0, amplitude=4.0)
    sim4.equilibrate()
    sim4.run(steps=200)
    rc = lfm.rotation_curve(sim4.chi, sim4.energy_density)  # correct arg order
    assert "r" in rc and "v_circ" in rc
    assert len(rc["r"]) > 5
    v_max = float(np.max(np.abs(np.asarray(rc["v_circ"])[np.isfinite(np.asarray(rc["v_circ"]))])))
    record(
        "Bob",
        "2b",
        "Disk sim rotation_curve()",
        PASS,
        f"v_max={v_max:.4f} LFM units, {len(rc['r'])} radial bins",
    )
except Exception as e:
    record("Bob", "2b", "Disk sim rotation_curve()", FAIL, str(e))

try:
    # Exp 2c: Fit sim rotation curve to SPARC DDO154
    sim5 = lfm.Simulation(lfm.SimulationConfig(grid_size=48, field_level=lfm.FieldLevel.REAL))
    lfm.initialize_disk(sim5, n_solitons=20, r_inner=3.0, r_outer=15.0, amplitude=3.0)
    sim5.equilibrate()
    sim5.run(steps=200)
    rc5 = lfm.rotation_curve(sim5.chi, sim5.energy_density)
    row_ddo = lfm.sparc_load("DDO154")["DDO154"]
    fit = lfm.rotation_curve_fit(row_ddo, rc5["r"], rc5["v_circ"], n_tau=15)
    assert "tau_best" in fit and "chi2" in fit
    assert np.isfinite(fit["tau_best"])
    record(
        "Bob",
        "2c",
        "rotation_curve_fit() on DDO154",
        PASS,
        f"tau_best={fit['tau_best']:.1f}, chi2={fit['chi2']:.3f}",
    )
except Exception as e:
    record("Bob", "2c", "rotation_curve_fit() on DDO154", FAIL, str(e))


# =========================================================================
# USER 3 — Carol, ML researcher wanting to use LFM as a physics simulator
# Goal: run parameter sweeps, get CSV-like results
# =========================================================================
section("USER 3: Carol (ML researcher) — parameter sweeps")

try:
    # Exp 3a: Sweep soliton amplitude and measure chi_min
    def run_and_measure(amplitude):
        s = lfm.Simulation(lfm.SimulationConfig(grid_size=24, field_level=lfm.FieldLevel.REAL))
        s.place_soliton((12, 12, 12), amplitude=amplitude, sigma=2.0)
        s.equilibrate()
        return lfm.chi_statistics(s.chi)["chi_min"]

    amplitudes = [3.0, 5.0, 7.0]
    chi_mins = [run_and_measure(a) for a in amplitudes]
    assert chi_mins[0] > chi_mins[1] > chi_mins[2], f"Expected decreasing chi_min, got {chi_mins}"
    record(
        "Carol",
        "3a",
        "Manual amplitude sweep, chi_min decreases",
        PASS,
        f"chi_min at amps {amplitudes}: {[f'{v:.2f}' for v in chi_mins]}",
    )
except Exception as e:
    record("Carol", "3a", "Manual amplitude sweep", FAIL, str(e))

try:
    # Exp 3b: Use lfm.sweep() for a one-parameter grid search
    sweep_results = lfm.sweep(
        config=lfm.SimulationConfig(grid_size=16, field_level=lfm.FieldLevel.REAL),
        param="chi0",
        values=[19.0, 18.5, 18.0],
        steps=50,
        metric_names=["chi_min"],
        soliton={"position": (8, 8, 8), "amplitude": 4.0},
    )
    assert len(sweep_results) == 3
    for row in sweep_results:
        assert "chi_min" in row
    record(
        "Carol",
        "3b",
        "lfm.sweep() 3-point chi0 grid",
        PASS,
        f"results: {[(r['chi0'], round(r['chi_min'], 2)) for r in sweep_results]}",
    )
except Exception as e:
    record("Carol", "3b", "lfm.sweep()", FAIL, str(e))

try:
    # Exp 3c: Extract scalar metrics dict for ML pipeline logging
    sim6 = lfm.Simulation(lfm.SimulationConfig(grid_size=32, field_level=lfm.FieldLevel.REAL))
    sim6.place_soliton((16, 16, 16), amplitude=5.0)
    sim6.equilibrate()
    sim6.run(steps=100)
    # chi_statistics returns a clean flat dict — directly usable in ML pipelines
    stats = lfm.chi_statistics(sim6.chi)
    energy_stats = lfm.chi_statistics(sim6.energy_density)  # reuse for energy
    wf = lfm.well_fraction(sim6.chi)
    vf = lfm.void_fraction(sim6.chi)
    row = {
        **{f"chi_{k}": v for k, v in stats.items()},
        "well_fraction": float(wf),
        "void_fraction": float(vf),
        "energy_max": energy_stats["max"],
    }
    assert all(np.isfinite(v) for v in row.values() if isinstance(v, float))
    record(
        "Carol",
        "3c",
        "Flat scalar dict for ML pipeline",
        PASS,
        f"{len(row)} fields: {list(row.keys())}",
    )
except Exception as e:
    record("Carol", "3c", "Flat scalar dict for ML pipeline", FAIL, str(e))


# =========================================================================
# USER 4 — Dave, physics hobbyist wanting to simulate nuclear scattering
# Goal: two-soliton scatter, compare head-on vs glancing impact params
# =========================================================================
section("USER 4: Dave (hobbyist) — nuclear scattering")

try:
    # Exp 4a: Two solitons at different impact parameters
    def scatter_chi_min(b_offset):
        s = lfm.Simulation(lfm.SimulationConfig(grid_size=32, field_level=lfm.FieldLevel.REAL))
        s.place_soliton((8, 16, 16), amplitude=8.0, velocity=(0.04, 0.0, 0.0))
        s.place_soliton((24, 16 + b_offset, 16), amplitude=8.0, velocity=(-0.04, 0.0, 0.0))
        s.equilibrate()
        s.run(steps=400)
        return lfm.chi_statistics(s.chi)["chi_min"]

    chi_min_ho = scatter_chi_min(0)
    chi_min_gl = scatter_chi_min(4)
    assert chi_min_ho < 18.5
    assert chi_min_gl < 18.5
    record(
        "Dave",
        "4a",
        "Two-soliton scatter (head-on vs glancing)",
        PASS,
        f"head-on chi_min={chi_min_ho:.2f}, glancing chi_min={chi_min_gl:.2f}",
    )
except Exception as e:
    record("Dave", "4a", "Two-soliton scatter", FAIL, str(e))

try:
    # Exp 4b: compute_impact_parameter needs two trajectory dicts — API is not
    # discoverable without docs (no convenience to extract from Simulation)
    import inspect

    sig = inspect.signature(lfm.compute_impact_parameter)
    params = list(sig.parameters.keys())
    assert params == ["traj_i", "traj_j"], f"got {params}"
    record(
        "Dave",
        "4b",
        "compute_impact_parameter signature",
        GAP,
        "Needs traj_i, traj_j dicts; no convenience method to extract from Simulation",
    )
except Exception as e:
    record("Dave", "4b", "compute_impact_parameter", FAIL, str(e))

try:
    # Exp 4c: Track soliton peaks over time using n_peaks kwarg
    sim_traj = lfm.Simulation(lfm.SimulationConfig(grid_size=32, field_level=lfm.FieldLevel.REAL))
    sim_traj.place_soliton((8, 16, 16), amplitude=5.0, velocity=(0.02, 0.0, 0.0))
    sim_traj.place_soliton((24, 16, 16), amplitude=5.0, velocity=(-0.02, 0.0, 0.0))
    sim_traj.equilibrate()
    peaks_over_time = []
    for _ in range(10):
        sim_traj.run(steps=30)
        pk = lfm.find_peaks(sim_traj.energy_density, n_peaks=2)  # uses new alias
        peaks_over_time.append(pk)
    assert len(peaks_over_time) == 10
    record(
        "Dave",
        "4c",
        "Trajectories via find_peaks(n_peaks=2) alias",
        PASS,
        f"Tracked {len(peaks_over_time)} snapshots, {len(peaks_over_time[-1])} peaks last step",
    )
except Exception as e:
    record("Dave", "4c", "Trajectory tracking with find_peaks n_peaks alias", FAIL, str(e))


# =========================================================================
# USER 5 — Eve, cosmologist wanting to test early-universe structure
# Goal: random ICs, watch structure form, quantify dark matter halo
# =========================================================================
section("USER 5: Eve (cosmologist) — structure formation")

try:
    # Exp 5a: Random initial conditions, watch structure form
    rng = np.random.default_rng(42)
    sim7 = lfm.Simulation(
        lfm.SimulationConfig(grid_size=32, field_level=lfm.FieldLevel.REAL, chi0=19.0)
    )
    N = sim7.config.grid_size
    sim7.chi[:] = 19.0 + 0.1 * rng.standard_normal((N, N, N))
    for _ in range(4):
        pos = tuple(int(rng.integers(8, N - 8)) for _ in range(3))
        sim7.place_soliton(pos, amplitude=3.0 + rng.random() * 2.0, sigma=2.0)
    sim7.equilibrate()
    sim7.run(steps=300)
    wf = lfm.well_fraction(sim7.chi)
    vf = lfm.void_fraction(sim7.chi)
    clusters = lfm.count_clusters(sim7.chi)
    assert 0.0 <= wf <= 1.0
    assert 0.0 <= vf <= 1.0
    record(
        "Eve",
        "5a",
        "Structure formation: wells + voids + clusters",
        PASS,
        f"well_frac={wf:.3f}, void_frac={vf:.3f}, clusters={clusters}",
    )
except Exception as e:
    record("Eve", "5a", "Structure formation from random ICs", FAIL, str(e))

try:
    # Exp 5b: Dark matter halo — single central soliton, radial chi profile
    sim8 = lfm.Simulation(lfm.SimulationConfig(grid_size=48, field_level=lfm.FieldLevel.REAL))
    N8 = 48
    sim8.place_soliton((N8 // 2, N8 // 2, N8 // 2), amplitude=8.0, sigma=3.0)
    sim8.equilibrate()
    sim8.run(steps=400)
    prof = lfm.radial_profile(sim8.chi, center=(N8 // 2, N8 // 2, N8 // 2))
    assert "mean" in prof and "profile" in prof  # both aliases present
    chi_mean = np.asarray(prof["mean"])
    assert chi_mean[0] < 18.5, f"No DM well at centre? chi_mean[0]={chi_mean[0]:.2f}"
    record(
        "Eve",
        "5b",
        "DM halo radial_profile has 'mean' alias",
        PASS,
        f"chi dips to {np.min(chi_mean):.2f} near centre",
    )
except Exception as e:
    record("Eve", "5b", "Dark matter halo radial profile", FAIL, str(e))

try:
    # Exp 5c: Inclined disk via b_cells
    sim9 = lfm.Simulation(lfm.SimulationConfig(grid_size=48, field_level=lfm.FieldLevel.REAL))
    pos_flat = lfm.initialize_disk(
        sim9, n_solitons=20, r_inner=4.0, r_outer=16.0, amplitude=3.5, seed=7, b_cells=0.0
    )
    sim9b = lfm.Simulation(lfm.SimulationConfig(grid_size=48, field_level=lfm.FieldLevel.REAL))
    pos_incl = lfm.initialize_disk(
        sim9b, n_solitons=20, r_inner=4.0, r_outer=16.0, amplitude=3.5, seed=7, b_cells=5.0
    )
    center_flat = float(np.median(pos_flat[:, 2]))
    center_incl = float(np.median(pos_incl[:, 2]))
    assert abs(center_incl - center_flat - 5.0) < 0.5, (
        f"b_cells offset not applied: flat={center_flat:.1f}, incl={center_incl:.1f}"
    )
    record(
        "Eve",
        "5c",
        "Inclined disk via b_cells",
        PASS,
        f"z-plane flat={center_flat:.1f}, inclined={center_incl:.1f} (+5 cells)",
    )
except Exception as e:
    record("Eve", "5c", "Inclined disk via b_cells", FAIL, str(e))


# =========================================================================
# Summary
# =========================================================================
print()
print("=" * 60)
print("  MOCK USER TESTING SUMMARY")
print("=" * 60)
passed = [r for r in results if r["status"] == PASS]
failed = [r for r in results if r["status"] == FAIL]
gaps = [r for r in results if r["status"] == GAP]

print(f"  PASS : {len(passed):2d}/15")
print(f"  FAIL : {len(failed):2d}/15")
print(f"  GAP  : {len(gaps):2d}/15")

if failed:
    print()
    print("  Failures:")
    for r in failed:
        print(f"    [{r['user']:5s} {r['exp']}] {r['label']}")
        print(f"          {r['notes']}")

if gaps:
    print()
    print("  API gaps (feature works but needs improvement):")
    for r in gaps:
        print(f"    [{r['user']:5s} {r['exp']}] {r['label']}")
        print(f"          {r['notes']}")

print()
v1_criteria = len(failed) == 0
print(f"  v1.0 readiness: {'READY' if v1_criteria else 'NOT YET'}")
print("  (criterion: 0 failures across 15 user experiments)")
print()
sys.exit(0 if v1_criteria else 1)