Add WC stability metrics, live drift monitor, sandbox configs, and locked geometry generator

openwave/xperiments/m4_ewt/_launcher.py

@@ -549,10 +549,15 @@ def compute_wave_oscillation(state):
                 state.wave_field,
                 state.trackers,
             )
-        if state.frame % 60 == 0 or state.frame == 10:
+            mean_drift, active_wc = instrument.log_stability_metrics(
+                state.frame, state.wave_center
+            )
+            sampling.sample_stability_metrics(state.frame, mean_drift, active_wc)
+        if state.frame % 60 == 0 or state.frame == 1:
             instrument.log_timestep_data(
                 state.frame, state.wave_field, state.trackers, state.wave_center
             )
+            sampling.save_monitor_data()
 
     state.amp_global_rms = state.trackers.amp_global_emarms_am[None] * constants.ATTOMETER
     state.freq_global_avg = state.trackers.freq_global_avg_rHz[None] / constants.RONTOSECOND

openwave/xperiments/m4_ewt/utils/instrumentation.py

@@ -133,3 +133,65 @@ def to_float(val):
         log_entry["wc_positions"] = positions
 
     json_logger.log_timestep(log_entry)
+
+    import numpy as np
+
+
+_pairwise_ref = None
+_pairwise_ref_set = False
+
+
+def log_stability_metrics(timestep: int, wave_center) -> tuple:
+    """
+    Log WC stability metrics: mean pairwise distance drift and active WC count.
+    Returns (mean_drift, n_active) so that the caller can feed the live monitor.
+    """
+    import numpy as np
+
+    global _pairwise_ref, _pairwise_ref_set
+
+    positions = []
+    for i in range(wave_center.num_sources):
+        if wave_center.active[i] == 0:
+            continue
+        pos = wave_center.position_float[i]
+        x = float(pos[0])
+        y = float(pos[1])
+        z = float(pos[2])
+        if any(np.isnan([x, y, z])):
+            continue
+        positions.append([x, y, z])
+
+    n_active = len(positions)
+
+    if n_active < 2:
+        json_logger.log_timestep(
+            {
+                "timestep": timestep,
+                "mean_drift": None,
+                "active_wc": n_active,
+            }
+        )
+        return None, n_active
+
+    positions_np = np.array(positions, dtype=np.float64)
+    diff = positions_np[:, np.newaxis, :] - positions_np[np.newaxis, :, :]
+    dist = np.sqrt(np.sum(diff**2, axis=-1))
+
+    if not _pairwise_ref_set:
+        _pairwise_ref = dist.copy()
+        _pairwise_ref_set = True
+        mean_drift = 0.0
+    else:
+        drift_matrix = np.abs(dist - _pairwise_ref)
+        i_upper = np.triu_indices_from(drift_matrix, k=1)
+        mean_drift = float(np.mean(drift_matrix[i_upper]))
+
+    json_logger.log_timestep(
+        {
+            "timestep": timestep,
+            "mean_drift": mean_drift,
+            "active_wc": n_active,
+        }
+    )
+    return mean_drift, n_active

openwave/xperiments/m4_ewt/utils/live_monitor_viewer.py

@@ -30,9 +30,10 @@ def main(data_path_str):
 
     plt.ion()
     plt.style.use("dark_background")
-    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), facecolor=colormap.DARK_GRAY[1])
+    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 10), facecolor=colormap.DARK_GRAY[1])
     fig.suptitle("OPENWAVE Live Monitor (external)", fontsize=20, family="Monospace")
 
+    # Panel 1: Displacement & Amplitude
     (line_disp,) = ax1.plot(
         [], [], color=colormap.viridis_palette[2][1], linewidth=2, label="DISPLACEMENT (am)"
     )
@@ -54,6 +55,7 @@ def main(data_path_str):
     ax1.legend(loc="upper right")
     ax1.set_ylim(auto=True)
 
+    # Panel 2: Frequency
     (line_freq,) = ax2.plot(
         [], [], color=colormap.blueprint_palette[2][1], linewidth=2, label="FREQUENCY (rHz)"
     )
@@ -72,11 +74,19 @@ def main(data_path_str):
     ax2.legend(loc="upper right")
     ax2.set_ylim(auto=True)
 
+    # Panel 3: Pairwise Distance Drift
+    (line_drift,) = ax3.plot(
+        [], [], color=colormap.ironbow_palette[2][1], linewidth=2, label="MEAN DRIFT (vox)"
+    )
+    ax3.set_xlabel("Timestep", family="Monospace")
+    ax3.set_ylabel("Mean Pairwise Drift [vox]", family="Monospace")
+    ax3.set_title("PAIRWISE DISTANCE DRIFT", family="Monospace")
+    ax3.grid(True, alpha=0.3)
+    ax3.legend(loc="upper right")
+
     plt.tight_layout()
     fig.show()
 
-    # If the launcher dies without tearing us down (a crash, a force-quit), we
-    # get reparented and would otherwise loop forever holding a window open.
     parent_pid = os.getppid()
 
     while data_path.exists():
@@ -87,6 +97,8 @@ def main(data_path_str):
         try:
             with open(data_path, "r") as f:
                 data = json.load(f)
+
+            # Update panels 1-2
             ts = data.get("timesteps", [])
             if ts:
                 line_disp.set_data(ts, data.get("displacements", []))
@@ -98,7 +110,15 @@ def main(data_path_str):
                 ax1.autoscale_view(scaley=True)
                 ax2.relim()
                 ax2.autoscale_view(scaley=True)
-                fig.canvas.draw()
+
+            # Update panels 3-4
+            st_ts = data.get("stability_timesteps", [])
+            if st_ts:
+                line_drift.set_data(st_ts, data.get("mean_drifts", []))
+                ax3.set_xlim(st_ts[0], max(st_ts[-1], st_ts[0] + 1))
+                ax3.relim()
+                ax3.autoscale_view(scaley=True)
+            fig.canvas.draw()
         except Exception:
             pass
         plt.pause(0.5)

openwave/xperiments/m4_ewt/utils/plotting.py

@@ -34,6 +34,8 @@ def _read_timestep_data():
         "frequencies": [],
     }
     for rec in doc["data"]:
+        if "displacement_am" not in rec:
+            continue
         data["timesteps"].append(rec["timestep"])
         data["displacements"].append(rec["displacement_am"])
         data["amplitudes"].append(rec["amp_local_emarms_am"])
@@ -219,9 +221,74 @@ def plot_live_values():
     print("\nLive monitor plot saved to:\n", save_path, "\n")
 
 
+def plot_wc_drift():
+    """Plot mean pairwise distance drift over time."""
+    import json
+
+    log_path = json_logger._data_dir / json_logger._filename
+    if not log_path.exists():
+        return
+    with open(log_path) as f:
+        doc = json.load(f)
+
+    drift_data = [d for d in doc["data"] if d.get("mean_drift") is not None]
+    if not drift_data:
+        print("[plot_wc_drift] No drift data in log.")
+        return
+
+    ts = [d["timestep"] for d in drift_data]
+    drift = [d["mean_drift"] for d in drift_data]
+
+    plt.style.use("dark_background")
+    fig, ax = plt.subplots(figsize=(10, 4), facecolor=colormap.DARK_GRAY[1])
+    ax.plot(ts, drift, color=colormap.viridis_palette[2][1], linewidth=2)
+    ax.set_xlabel("Timestep")
+    ax.set_ylabel("Mean Pairwise Drift [vox]")
+    ax.set_title("Pairwise Distance Drift")
+    ax.grid(True, alpha=0.3)
+    plt.tight_layout()
+    PLOT_DIR.mkdir(parents=True, exist_ok=True)
+    plt.savefig(PLOT_DIR / "wc_drift.png", dpi=150, bbox_inches="tight")
+    print(f"[plot_wc_drift] Saved to {PLOT_DIR / 'wc_drift.png'}")
+
+
+def plot_wc_active():
+    """Plot number of active wave centers over time."""
+    import json
+
+    log_path = json_logger._data_dir / json_logger._filename
+    if not log_path.exists():
+        return
+    with open(log_path) as f:
+        doc = json.load(f)
+
+    active_data = [d for d in doc["data"] if "active_wc" in d]
+    if not active_data:
+        print("[plot_wc_active] No active WC data in log.")
+        return
+
+    ts = [d["timestep"] for d in active_data]
+    active = [d["active_wc"] for d in active_data]
+
+    plt.style.use("dark_background")
+    fig, ax = plt.subplots(figsize=(10, 4), facecolor=colormap.DARK_GRAY[1])
+    ax.plot(ts, active, color=colormap.ironbow_palette[2][1], linewidth=2)
+    ax.set_xlabel("Timestep")
+    ax.set_ylabel("Active WC count")
+    ax.set_title("Active Wave Centers")
+    ax.grid(True, alpha=0.3)
+    ax.set_ylim(bottom=0)
+    plt.tight_layout()
+    PLOT_DIR.mkdir(parents=True, exist_ok=True)
+    plt.savefig(PLOT_DIR / "wc_active.png", dpi=150, bbox_inches="tight")
+    print(f"[plot_wc_active] Saved to {PLOT_DIR / 'wc_active.png'}")
+
+
 def generate_plots():
     """Generate all instrumentation plots. Called after simulation ends."""
     json_logger.finalize()
     plot_probe_values()
     plot_live_values()
+    plot_wc_drift()
+    # plot_wc_active()
     plt.show()

openwave/xperiments/m4_ewt/utils/sampling.py

@@ -5,6 +5,9 @@
 _plot_displacements = []
 _plot_amplitudes = []
 _plot_frequencies = []
+_stability_timesteps = []
+_stability_drifts = []
+_stability_active = []
 
 # File shared with the external live-monitor process
 # The model root, one level above utils/. This must match the path the launcher
@@ -13,6 +16,15 @@
 SAVE_EVERY = 10  # write the shared file every 10 samples
 
 
+def sample_stability_metrics(timestep, mean_drift, active_wc):
+    _stability_timesteps.append(timestep)
+    _stability_drifts.append(mean_drift if mean_drift is not None else 0.0)
+    _stability_active.append(active_wc)
+
+    if len(_stability_timesteps) % SAVE_EVERY == 0:
+        save_monitor_data()
+
+
 def sample_for_plots(timestep, wave_field, trackers):
     px = wave_field.nx // 2 + 1
     py = wave_field.ny // 2
@@ -49,6 +61,9 @@ def save_monitor_data():
         "displacements": _plot_displacements,
         "amplitudes": _plot_amplitudes,
         "frequencies": _plot_frequencies,
+        "stability_timesteps": _stability_timesteps,
+        "mean_drifts": _stability_drifts,
+        "active_wcs": _stability_active,
     }
     MONITOR_DATA_PATH.parent.mkdir(parents=True, exist_ok=True)
     with open(MONITOR_DATA_PATH, "w") as f:

openwave/xperiments/m4_ewt/xparameters/sandbox_k_selectivity/electron_k09_ewt_vmode10_neumann_boost0015.py

@@ -0,0 +1,70 @@
+# electron_k09_ewt_vmode10_neumann_boost0015.py
+"""
+K=9, tricapped trigonal prism (EWT geometry), neumann interaction, V_MODE=10, WC_BOOST=0.015.
+K-selectivity test: K=9 EWT geometry with V_MODE=10.
+"""
+
+from openwave.xperiments.m4_ewt.xparameters.geometry import generate_positions_by_EWT_geometry
+from openwave.common import constants
+
+UNIVERSE_EDGE = 2e-15
+TARGET_VOXELS = 55_000_000
+K = 9
+PERTURBATION = 0.02
+
+POSITIONS = generate_positions_by_EWT_geometry(
+    UNIVERSE_EDGE, K, center=(0.5, 0.5, 0.5), perturbation=PERTURBATION
+)
+PHASES = [180] * K
+
+XPARAMETERS = {
+    "meta": {
+        "X_NAME": "K=9 EWT tricapped prism neumann boost=0.015 V_MODE=10",
+        "DESCRIPTION": "K-selectivity test: K=9 EWT geometry with V_MODE=10 and neumann interaction",
+    },
+    "camera": {"INITIAL_POSITION": [1.20, 1.80, 1.20]},
+    "universe": {
+        "SIZE": [UNIVERSE_EDGE, UNIVERSE_EDGE, UNIVERSE_EDGE],
+        "TARGET_VOXELS": TARGET_VOXELS,
+    },
+    "wave_centers": {
+        "COUNT": K,
+        "POSITION": POSITIONS,
+        "PHASE_OFFSETS_DEG": PHASES,
+        "APPLY_MOTION": True,
+    },
+    "engine": {
+        "SEED_MODE": 2,
+        "SEED_BOOST": 0.01,
+        "V_MODE": 10,
+        "V_C1": -0.3,
+        "V_C2": 0.05,
+        "WC_INTERACT_MODE": 2,
+        "WC_BOOST": 0.015,
+        "WC_RADIUS": 2,
+        "WC_SIGMA": 1.5,
+        "VELOCITY_DAMPING": 0.990,
+        "R_WALL": 100.0,
+        "WALL_HEIGHT": 1.2,
+        "DEFICIT_DEPTH": 0.9,
+        "R_SOLITON": 35.0,
+        "SIGMA": 3.0,
+        "PRESSURE_STRENGTH": 0.0,
+        "CFL_SAFETY": 0.1,
+    },
+    "ui_defaults": {
+        "SHOW_AXIS": False,
+        "TICK_SPACING": 0.25,
+        "SHOW_GRID": False,
+        "SHOW_EDGES": False,
+        "FLUX_MESH_PLANES": [0.5, 0.5, 0.5],
+        "SHOW_FLUX_MESH": 1,
+        "WARP_MESH": 100,
+        "SHOW_GRANULES": False,
+        "PARTICLE_SHELL": False,
+        "SIM_SPEED": 1.0,
+        "PAUSED": False,
+    },
+    "color_defaults": {"COLOR_THEME": "OCEAN", "WAVE_MENU": 1},
+    "analytics": {"INSTRUMENTATION": True, "EXPORT_VIDEO": False, "VIDEO_FRAMES": 24},
+}