#!/usr/bin/env python3 """ plasmoid_entrainment_sim.py (single-file, runnable, UPDATED) What changed in this version (to keep it simple but more “physics-like”): ✅ Default scan ratios widened: 0.6 .. 2.0 step 0.02 ✅ “More physics-like tongue” default: drive acts as THRESHOLD MODULATION (gating), not direct additive forcing ✅ Optional refractory time after bursts (default on) to prevent immediate re-triggering ✅ Q metric fixed: Q is now measured near the DRIVE frequency (Q_drive), not the global PSD peak (baseline f0 still estimated from global PSD peak) Dependencies: pip install numpy scipy matplotlib Run: python plasmoid_entrainment_sim.py """ import argparse import json import math import sys from dataclasses import dataclass from datetime import datetime from pathlib import Path from typing import Optional, Tuple, List from concurrent.futures import ProcessPoolExecutor, as_completed import numpy as np import matplotlib.pyplot as plt from scipy.signal import welch TWO_PI = 2.0 * np.pi # ============================================================ # Resonant math primitives (Candidate A/B + operational closure) # ============================================================ @dataclass class LatticeDeficitResult: C_lat: float m_star: np.ndarray residual: np.ndarray def closure_deficit_lattice(Delta: np.ndarray, n: int, search_radius: int = 1) -> LatticeDeficitResult: """ Candidate A: lattice-distance closure deficit in phase increments. Resonant math principle: C_lat(Δ) = min_{m∈Z^3} || nΔ - 2π m || """ Delta = np.asarray(Delta, dtype=float).reshape(-1) if Delta.shape[0] != 3: raise ValueError("Delta must be length 3: (ΔΘ, ΔΞ, ΔΦ).") if n <= 0: raise ValueError("n must be positive.") y = (n * Delta) / TWO_PI m0 = np.round(y).astype(int) best_m = None best_res = None best_norm = np.inf r = int(search_radius) offsets = np.array(np.meshgrid( np.arange(-r, r + 1), np.arange(-r, r + 1), np.arange(-r, r + 1), indexing="ij" )).reshape(3, -1).T for off in offsets: m = m0 + off res = n * Delta - TWO_PI * m norm = float(np.linalg.norm(res)) if norm < best_norm: best_norm = norm best_m = m best_res = res return LatticeDeficitResult(C_lat=best_norm, m_star=np.array(best_m, dtype=int), residual=np.array(best_res, dtype=float)) def closure_deficit_lattice_extended(Delta4: np.ndarray, n: int, search_radius: int = 1) -> LatticeDeficitResult: """ Extended lattice-distance closure deficit including a transverse increment. This is ONLY meaningful when a control protocol advances a transverse parameter in lockstep with intrinsic phases. C_lat^ext(Δ) = min_{m∈Z^4} || nΔ - 2π m ||, Δ=(ΔΘ,ΔΞ,ΔΦ,Δτ) """ Delta4 = np.asarray(Delta4, dtype=float).reshape(-1) if Delta4.shape[0] != 4: raise ValueError("Delta4 must be length 4: (ΔΘ, ΔΞ, ΔΦ, Δτ).") if n <= 0: raise ValueError("n must be positive.") y = (n * Delta4) / TWO_PI m0 = np.round(y).astype(int) best_m = None best_res = None best_norm = np.inf r = int(search_radius) offsets = np.array(np.meshgrid( np.arange(-r, r + 1), np.arange(-r, r + 1), np.arange(-r, r + 1), np.arange(-r, r + 1), indexing="ij" )).reshape(4, -1).T for off in offsets: m = m0 + off res = n * Delta4 - TWO_PI * m norm = float(np.linalg.norm(res)) if norm < best_norm: best_norm = norm best_m = m best_res = res return LatticeDeficitResult(C_lat=best_norm, m_star=np.array(best_m, dtype=int), residual=np.array(best_res, dtype=float)) def lattice_basin_boundary_eps(n: int, v_norm: float) -> float: """ Minimal perturbation model boundary estimate for Candidate A. Resonant math principle: |ε| ~ π / (n ||v||) """ if n <= 0 or v_norm <= 0: return float("nan") return float(np.pi / (n * v_norm)) def closure_deficit_geom(R: float, a: float, tau0: float, eps: float) -> float: """ Candidate B: geometric mismatch deficit for detuning k -> k + ε at fixed τ0. Resonant math principle: ρ(τ0) = R + a cos τ0 C_geom(τ0; ε) = 2 ρ(τ0) |sin(π ε)| """ rho = R + a * math.cos(tau0) return float(2.0 * rho * abs(math.sin(math.pi * eps))) def phase_locking_value(event_times: np.ndarray, fd: float) -> float: """ Operational closure on S^1: PLV measures phase-consistency of events. Resonant math principle: PLV = | (1/N) Σ exp(i 2π f_d t_k) | """ event_times = np.asarray(event_times, dtype=float) if fd <= 0 or event_times.size < 3: return float("nan") phi = (TWO_PI * fd * event_times) % TWO_PI return float(np.abs(np.mean(np.exp(1j * phi)))) def closure_deficit_phase(event_times: np.ndarray, fd: float) -> float: """ Resonant math principle: C_phase = 1 - PLV """ plv = phase_locking_value(event_times, fd) if not np.isfinite(plv): return float("nan") return float(1.0 - plv) # ========================= # Toy "plasmoid" simulator # ========================= @dataclass class SimOut: t: np.ndarray R: np.ndarray event_times: np.ndarray # threshold-crossing trigger times (TRUE events) def simulate_bursty_reconnection( T: float, dt: float, a: float, b: float, thresh: float, reset: float, burst_amp: float, burst_tau: float, noise_sigma: float, drive_type: str, drive_amp: float, fd: float, drive_phase: float, bb_band: Tuple[float, float], seed: int, drive_mode: str = "threshold", refractory: float = 0.25, drive_rms: Optional[float] = None, scramble_tau: float = 0.0 ) -> SimOut: """ Bursty reconnection-rate surrogate. drive_mode: - "additive": x += (... + drive) dt - "parametric": a_eff = a*(1+drive); x += (a_eff - b x) dt - "threshold": burst trigger threshold modulated: thresh_eff = thresh*(1-drive) refractory: minimum time (s) after a burst before another can trigger (physics-like gating). """ rng = np.random.default_rng(seed) n = int(T / dt) t = np.arange(n) * dt # Drive signal: # For sine/broadband we treat drive as a dimensionless modulation signal in [-A, A]. if drive_type == "none": drive = np.zeros(n) elif drive_type == "sine": drive = drive_amp * np.sin(TWO_PI * fd * t + drive_phase) elif drive_type == "broadband": wn = rng.normal(0.0, 1.0, size=n) W = np.fft.rfft(wn) f = np.fft.rfftfreq(n, dt) mask = (f >= bb_band[0]) & (f <= bb_band[1]) W[~mask] = 0.0 bb = np.fft.irfft(W, n=n) bb = bb / (np.std(bb) + 1e-12) drive = drive_amp * bb elif drive_type == "scrambled": # Narrowband carrier at fd whose phase is re-randomised every scramble_tau seconds. # Same centre frequency and same power as "sine"; long-range phase coherence removed. seg = max(1, int(round((scramble_tau if scramble_tau > 0 else 2.0 / max(fd, 1e-9)) / dt))) drive = np.zeros(n) for start in range(0, n, seg): stop = min(n, start + seg) drive[start:stop] = np.sin(TWO_PI * fd * t[start:stop] + rng.uniform(0.0, TWO_PI)) drive = drive_amp * drive elif drive_type == "dc": drive = drive_amp * np.ones(n) else: raise ValueError("drive_type must be one of: none, sine, broadband, scrambled, dc") # Exact power matching across drive shapes: fix mean-square, not peak amplitude. if drive_rms is not None and drive_type != "none": if drive_type == "dc": drive = drive_rms * np.ones(n) else: drive = drive - float(np.mean(drive)) drive = drive_rms * drive / (float(np.std(drive)) + 1e-12) # Output signal (proxy reconnection rate) R = np.zeros(n, dtype=float) # Internal "stress" variable that loads up and triggers bursts x = 0.0 event_times: List[float] = [] # Burst kernel (a short exponential decay added to R when an event triggers) klen = max(8, int(5.0 * burst_tau / dt)) kt = np.arange(klen) * dt kernel = burst_amp * np.exp(-kt / burst_tau) next_allowed = 0.0 if refractory > 0 else -np.inf for i in range(n): di = float(drive[i]) # --- Dynamics update --- if drive_mode == "additive": # (less physics-like when drive is strong; kept for reference) x += (a - b * x + di) * dt + noise_sigma * math.sqrt(dt) * rng.normal() elif drive_mode == "parametric": # drive modulates the loading rate a (dimensionless modulation) a_eff = a * (1.0 + di) x += (a_eff - b * x) * dt + noise_sigma * math.sqrt(dt) * rng.normal() elif drive_mode == "threshold": # no direct injection; instead, the drive "opens/closes the gate" x += (a - b * x) * dt + noise_sigma * math.sqrt(dt) * rng.normal() else: raise ValueError('drive_mode must be "additive", "parametric", or "threshold"') # --- Trigger condition --- if drive_mode == "threshold": # thresh_eff = thresh*(1 - drive) => when drive positive, threshold is lowered (more likely to burst). # clamp to avoid negative/too-small thresholds thresh_eff = thresh * (1.0 - di) thresh_eff = max(0.2 * thresh, min(2.0 * thresh, thresh_eff)) else: thresh_eff = thresh if (t[i] >= next_allowed) and (x >= thresh_eff): x = reset event_times.append(float(t[i])) next_allowed = float(t[i] + refractory) if refractory > 0 else next_allowed jmax = min(n, i + klen) R[i:jmax] += kernel[:(jmax - i)] # baseline offset R += 0.02 return SimOut(t=t, R=R, event_times=np.asarray(event_times, dtype=float)) # ========================= # Metrics / Analysis helpers # ========================= @dataclass class SpectralQ: f_peak: float q: float peak_power: float half_power_bw: float def psd_global_peak(t: np.ndarray, x: np.ndarray, dt: float) -> Tuple[np.ndarray, np.ndarray]: mu = float(np.mean(x)) fs = 1.0 / dt f, pxx = welch(x - mu, fs=fs, nperseg=min(4096, len(x))) return f, pxx def psd_q_at_target(t: np.ndarray, x: np.ndarray, dt: float, f_target: float, window_frac: float = 0.15) -> SpectralQ: """ Meaningful Q for driven systems: compute Q around a target frequency (typically fd), not global PSD peak. Steps: 1) compute PSD (Welch) 2) restrict to a window around f_target 3) find local peak inside that window 4) compute half-power bandwidth around that peak 5) Q = f_peak / bw """ if f_target <= 0: return SpectralQ(f_peak=np.nan, q=np.nan, peak_power=np.nan, half_power_bw=np.nan) f, pxx = psd_global_peak(t, x, dt) if len(f) < 8: return SpectralQ(f_peak=np.nan, q=np.nan, peak_power=np.nan, half_power_bw=np.nan) fmin = f_target * (1.0 - window_frac) fmax = f_target * (1.0 + window_frac) mask = (f >= fmin) & (f <= fmax) if not np.any(mask): return SpectralQ(f_peak=np.nan, q=np.nan, peak_power=np.nan, half_power_bw=np.nan) fi = f[mask] pi = pxx[mask] if fi.size < 3: return SpectralQ(f_peak=np.nan, q=np.nan, peak_power=np.nan, half_power_bw=np.nan) i_peak_local = int(np.argmax(pi)) f_peak = float(fi[i_peak_local]) p0 = float(pi[i_peak_local]) half = p0 / 2.0 iL = i_peak_local while iL > 0 and pi[iL] > half: iL -= 1 iR = i_peak_local while iR < len(pi) - 1 and pi[iR] > half: iR += 1 bw = float(fi[iR] - fi[iL]) if iR > iL else float("nan") q = float(f_peak / bw) if (np.isfinite(bw) and bw > 0) else float("nan") return SpectralQ(f_peak=f_peak, q=q, peak_power=p0, half_power_bw=bw) def psd_peak_and_q_global(t: np.ndarray, x: np.ndarray, dt: float) -> SpectralQ: """ Global-peak Q (kept only for baseline f0 estimation). """ f, pxx = psd_global_peak(t, x, dt) if len(f) < 4: return SpectralQ(f_peak=np.nan, q=np.nan, peak_power=np.nan, half_power_bw=np.nan) idx0 = 1 # skip DC i_peak = idx0 + int(np.argmax(pxx[idx0:])) f0 = float(f[i_peak]) p0 = float(pxx[i_peak]) half = p0 / 2.0 iL = i_peak while iL > idx0 and pxx[iL] > half: iL -= 1 iR = i_peak while iR < len(pxx) - 1 and pxx[iR] > half: iR += 1 bw = float(f[iR] - f[iL]) if iR > iL else float("nan") q = float(f0 / bw) if (np.isfinite(bw) and bw > 0) else float("nan") return SpectralQ(f_peak=f0, q=q, peak_power=p0, half_power_bw=bw) def intermittency_metrics(x: np.ndarray) -> Tuple[float, float, float]: mu = float(np.mean(x)) sd = float(np.std(x)) cv = float(sd / mu) if mu != 0 else float("nan") return mu, sd, cv def isi_cv(event_times: np.ndarray) -> float: """ Drive-independent organisation metric: coefficient of variation of the inter-event interval. Makes no reference to f_d, so it can be compared across drives of different spectral character. Lower = more regular. """ event_times = np.asarray(event_times, dtype=float) if event_times.size < 5: return float("nan") d = np.diff(event_times) m = float(np.mean(d)) return float(np.std(d) / m) if m > 0 else float("nan") def excess_kurtosis(x: np.ndarray) -> float: x = np.asarray(x, dtype=float) mu = np.mean(x) y = x - mu m2 = np.mean(y * y) m4 = np.mean((y * y) ** 2) if m2 <= 1e-30: return float("nan") return float(m4 / (m2 * m2) - 3.0) EFFICIENCY_ARMS = ("coherent", "scrambled", "off_resonant", "broadband", "dc") def _efficiency_cell(args): arm, rms, seed, f0, T, dt, noise_sigma, sim_params, drive_mode, refractory, offfactor = args drive_type = {"coherent": "sine", "scrambled": "scrambled", "off_resonant": "sine", "broadband": "broadband", "dc": "dc"}[arm] fd = f0 * (offfactor if arm == "off_resonant" else 1.0) out = simulate_bursty_reconnection( T=T, dt=dt, seed=seed, drive_type=drive_type, drive_amp=1.0, fd=fd, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=noise_sigma, drive_mode=drive_mode, refractory=refractory, drive_rms=rms, **sim_params ) return arm, rms, isi_cv(out.event_times), len(out.event_times) def run_efficiency(T, dt, f0, noise_sigma, sim_params, drive_mode, refractory, offfactor, rms_levels, nrep, seed0, targets): """ Equal-target energy comparison: minimum drive power needed to reach a drive-independent organisation target, per drive shape, at matched power. """ base = simulate_bursty_reconnection( T=T, dt=dt, seed=seed0, drive_type="none", drive_amp=0.0, fd=1.0, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=noise_sigma, drive_mode="additive", refractory=refractory, **sim_params ) baseline = isi_cv(base.event_times) jobs = [(arm, rms, seed0 + 1 + r, f0, T, dt, noise_sigma, sim_params, drive_mode, refractory, offfactor) for arm in EFFICIENCY_ARMS for rms in rms_levels for r in range(nrep)] results = {} with ProcessPoolExecutor() as ex: for arm, rms, cv, nev in ex.map(_efficiency_cell, jobs, chunksize=8): results.setdefault((arm, rms), []).append((cv, nev)) grid = {} for (arm, rms), vals in results.items(): cvs = [v[0] for v in vals if np.isfinite(v[0])] grid[(arm, rms)] = { "power": float(rms ** 2), "isi_cv": float(np.mean(cvs)) if cvs else float("nan"), "isi_cv_sd": float(np.std(cvs)) if cvs else float("nan"), "events": float(np.mean([v[1] for v in vals])), } table = {} for frac in targets: target = baseline * frac mins = {} for arm in EFFICIENCY_ARMS: hit = [grid[(arm, r)]["power"] for r in rms_levels if np.isfinite(grid[(arm, r)]["isi_cv"]) and grid[(arm, r)]["isi_cv"] <= target] mins[arm] = min(hit) if hit else None table[frac] = {"target_isi_cv": float(target), "min_power": mins} return {"baseline_isi_cv": float(baseline), "f0": float(f0), "rms_levels": list(rms_levels), "nrep": nrep, "grid": {f"{a}|{r}": v for (a, r), v in grid.items()}, "targets": {str(k): v for k, v in table.items()}} def _coherence_cell(args): tau, rms, seed, f0, T, dt, noise_sigma, sim_params, drive_mode, refractory = args drive_type = "sine" if tau is None else "scrambled" out = simulate_bursty_reconnection( T=T, dt=dt, seed=seed, drive_type=drive_type, drive_amp=1.0, fd=f0, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=noise_sigma, drive_mode=drive_mode, refractory=refractory, drive_rms=rms, scramble_tau=(tau or 0.0), **sim_params ) return tau, isi_cv(out.event_times) def run_coherence(T, dt, f0, noise_sigma, sim_params, drive_mode, refractory, rms, taus, nrep, seed0): """ Same power, same centre frequency, varying drive coherence time. Separates 'phase alignment matters' from 'my scramble period was too short'. """ base = simulate_bursty_reconnection( T=T, dt=dt, seed=seed0, drive_type="none", drive_amp=0.0, fd=1.0, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=noise_sigma, drive_mode="additive", refractory=refractory, **sim_params ) baseline = isi_cv(base.event_times) mean_isi = float(np.mean(np.diff(base.event_times))) jobs = [(tau, rms, seed0 + 1 + r, f0, T, dt, noise_sigma, sim_params, drive_mode, refractory) for tau in list(taus) + [None] for r in range(nrep)] acc = {} with ProcessPoolExecutor() as ex: for tau, cv in ex.map(_coherence_cell, jobs, chunksize=8): acc.setdefault(tau, []).append(cv) curve = [] for tau in list(taus) + [None]: vals = [v for v in acc[tau] if np.isfinite(v)] curve.append({"tau_c": tau, "tau_over_isi": (tau / mean_isi) if tau else None, "isi_cv": float(np.mean(vals)), "isi_cv_sd": float(np.std(vals))}) return {"baseline_isi_cv": float(baseline), "mean_isi": mean_isi, "power": float(rms ** 2), "f0": float(f0), "curve": curve} def print_coherence(res): print("\n[COHERENCE] fixed power, fixed centre frequency, varying drive coherence time") print(f" baseline ISI CV = {res['baseline_isi_cv']:.4f} power = {res['power']:.4f}" f" mean ISI = {res['mean_isi']:.3f}s") for row in res["curve"]: label = "inf (pure sine)" if row["tau_c"] is None else f"{row['tau_c']:g}s" ratio = "" if row["tau_over_isi"] is None else f" ({row['tau_over_isi']:.1f} x ISI)" print(f" tau_c = {label:<16} ISI CV = {row['isi_cv']:.4f} +/- {row['isi_cv_sd']:.4f}{ratio}") def print_efficiency(res): print("\n[EFFICIENCY] minimum drive power to reach a drive-independent target") print(f" baseline ISI CV = {res['baseline_isi_cv']:.4f} f0 = {res['f0']:.4f}") for frac, block in res["targets"].items(): mins = block["min_power"] ref = mins.get("coherent") print(f"\n target ISI CV <= {block['target_isi_cv']:.4f} ({float(frac):.0%} of baseline)") for arm in EFFICIENCY_ARMS: p = mins[arm] if p is None: cap = max(res["rms_levels"]) ** 2 rel = f">{cap / ref:.0f}x" if ref else "n/a" print(f" {arm:<14} not reached within swept range {rel}") else: rel = f"{p / ref:.1f}x" if ref else "1.0x" print(f" {arm:<14} power = {p:.5f} {rel}") def print_progress(msg: str) -> None: print(msg) sys.stdout.flush() # ========================= # Sweep runner # ========================= @dataclass class SweepResult: PLV: np.ndarray Cphase: np.ndarray CV: np.ndarray KurtEx: np.ndarray Qdrive: np.ndarray # Q near fd (meaningful) EventsMean: np.ndarray EventsStd: np.ndarray FRI: np.ndarray # False Resonance Index = Qdrive * (1 - PLV) def _sweep_cell_task(args): """Worker for one (A, ratio) grid cell. Returns (ia, ir, metrics...).""" (ia, ir, A, r, f0, T, dt, nrep, noise_sigma, sim_params, drive_mode, refractory, progress_prefix) = args fd = r * f0 plvs, cphs, cvs, ks, qds, evs = [], [], [], [], [], [] for rep in range(nrep): seed = 10000 + ia * 1000 + ir * 50 + rep + int(noise_sigma * 1000) out = simulate_bursty_reconnection( T=T, dt=dt, seed=seed, drive_type="sine", drive_amp=A, fd=fd, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=noise_sigma, drive_mode=drive_mode, refractory=refractory, **sim_params ) evs.append(len(out.event_times)) plv = phase_locking_value(out.event_times, fd) if np.isfinite(plv): plvs.append(plv) cphs.append(1.0 - plv) _, _, cv = intermittency_metrics(out.R) kx = excess_kurtosis(out.R) qd = psd_q_at_target(out.t, out.R, dt=dt, f_target=fd, window_frac=0.15).q if np.isfinite(cv): cvs.append(cv) if np.isfinite(kx): ks.append(kx) if np.isfinite(qd): qds.append(qd) plv_m = float(np.mean(plvs)) if plvs else np.nan cph_m = float(np.mean(cphs)) if cphs else np.nan cv_m = float(np.mean(cvs)) if cvs else np.nan k_m = float(np.mean(ks)) if ks else np.nan qd_m = float(np.mean(qds)) if qds else np.nan ev_m = float(np.mean(evs)) if evs else np.nan ev_s = float(np.std(evs, ddof=1)) if len(evs) > 1 else 0.0 fri = (qd_m) * (1.0 - plv_m) if np.isfinite(qd_m) and np.isfinite(plv_m) else np.nan return (ia, ir, plv_m, cph_m, cv_m, k_m, qd_m, ev_m, ev_s, float(fri) if np.isfinite(fri) else np.nan) def run_sweep( *, T: float, dt: float, f0: float, amps: List[float], ratios: List[float], nrep: int, noise_sigma: float, sim_params: dict, drive_mode: str, refractory: float, progress_prefix: str = "" ) -> SweepResult: PLV_map = np.full((len(amps), len(ratios)), np.nan, dtype=float) Cphase_map = np.full_like(PLV_map, np.nan) CV_map = np.full_like(PLV_map, np.nan) K_map = np.full_like(PLV_map, np.nan) Qd_map = np.full_like(PLV_map, np.nan) EventsMean_map = np.full_like(PLV_map, np.nan) EventsStd_map = np.full_like(PLV_map, np.nan) FRI_map = np.full_like(PLV_map, np.nan) # Parallelize over grid cells (A, fd/f0). Max workers fixed to 5. tasks = [] for ia, A in enumerate(amps): for ir, r in enumerate(ratios): tasks.append((ia, ir, A, r, f0, T, dt, nrep, noise_sigma, sim_params, drive_mode, refractory, progress_prefix)) total_cells = len(tasks) done = 0 step = max(1, total_cells // 10) # Print sweep header per amplitude (keeps readability similar to serial version) for ia, A in enumerate(amps): print_progress(f"{progress_prefix}Sweep: A={A:.3g} ({ia+1}/{len(amps)}) noise_sigma={noise_sigma:.3g}") with ProcessPoolExecutor(max_workers=16) as ex: futs = [ex.submit(_sweep_cell_task, t) for t in tasks] for fut in as_completed(futs): ia, ir, plv_m, cph_m, cv_m, k_m, qd_m, ev_m, ev_s, fri_m = fut.result() PLV_map[ia, ir] = plv_m Cphase_map[ia, ir] = cph_m CV_map[ia, ir] = cv_m K_map[ia, ir] = k_m Qd_map[ia, ir] = qd_m EventsMean_map[ia, ir] = ev_m EventsStd_map[ia, ir] = ev_s FRI_map[ia, ir] = fri_m done += 1 if done % step == 0 or done == total_cells: pct = 100.0 * done / total_cells print_progress(f"{progress_prefix} progress: {pct:.0f}% ({done}/{total_cells})") return SweepResult( PLV=PLV_map, Cphase=Cphase_map, CV=CV_map, KurtEx=K_map, Qdrive=Qd_map, EventsMean=EventsMean_map, EventsStd=EventsStd_map, FRI=FRI_map ) def summarize_event_counts(label: str, res: SweepResult, amps: List[float], ratios: List[float]) -> None: ev = res.EventsMean if ev is None: return finite = np.isfinite(ev) if not np.any(finite): print(f"[{label}] Event counts: no finite data") return ev_vals = ev[finite] print(f"[{label}] Event counts (mean per cell): min={np.min(ev_vals):.1f} median={np.median(ev_vals):.1f} max={np.max(ev_vals):.1f}") # Per-amplitude summary for i, A in enumerate(amps): row = ev[i, :] fin = np.isfinite(row) if np.any(fin): print(f"[{label}] A={A:.3g}: events_mean median={np.median(row[fin]):.1f} (range {np.min(row[fin]):.1f}–{np.max(row[fin]):.1f})") def tongue_widths(PLV_map: np.ndarray, amps: List[float], ratios: List[float], threshold: float = 0.8): ratios_arr = np.asarray(ratios, dtype=float) out = [] for i, A in enumerate(amps): row = PLV_map[i, :] mask = np.isfinite(row) & (row >= threshold) if not np.any(mask): out.append((A, 0.0, None, None)) else: rmin = float(np.min(ratios_arr[mask])) rmax = float(np.max(ratios_arr[mask])) out.append((A, rmax - rmin, rmin, rmax)) return out # ========================= # Plot helper: correct figure colorbars (no warnings) # ========================= def heat(fig, im_ax, Z, title, ratios, amps, vmin=None, vmax=None): im = im_ax.imshow( Z, aspect="auto", origin="lower", extent=[min(ratios), max(ratios), min(amps), max(amps)], vmin=vmin, vmax=vmax ) im_ax.set_title(title) im_ax.set_xlabel("fd / f0") im_ax.set_ylabel("drive amplitude A") fig.colorbar(im, ax=im_ax, fraction=0.046, pad=0.04) return im # ========================= # Main # ========================= def main(): ap = argparse.ArgumentParser() ap.add_argument("--T", type=float, default=200.0) ap.add_argument("--dt", type=float, default=1e-2) ap.add_argument("--no-plots", action="store_true") ap.add_argument("--save-validation", action="store_true", help="Also save validation-only plots (e.g., event-count heatmaps).") ap.add_argument("--seed0", type=int, default=1) ap.add_argument("--nrep", type=int, default=10) # Smaller amplitudes are more appropriate for threshold/parametric modulation. ap.add_argument("--amps", type=str, default="0.0,0.03,0.06,0.09,0.12") # Default ratios widened: 0.6 .. 2.0 step 0.02 (inclusive-ish by construction). ap.add_argument("--ratios", type=str, default=",".join(map(str, np.arange(0.6, 2.0, 0.02)))) ap.add_argument("--noise", type=float, default=0.12) # main sweep noise ap.add_argument("--noise-robust", type=float, default=0.18) # robustness sweep noise ap.add_argument("--demoA", type=float, default=0.09) # demo amplitude for one-off comparisons ap.add_argument("--offfactor", type=float, default=1.6) ap.add_argument("--outdir", type=str, default="") # if empty => timestamped # Simple controls for “physics-like” behavior ap.add_argument("--drive-mode", type=str, default="threshold", choices=["threshold", "parametric", "additive"]) ap.add_argument("--refractory", type=float, default=0.25) # Equal-target energy comparison (Paper II Sec. 10.4) ap.add_argument("--efficiency", action="store_true", help="Run the equal-target energy comparison and exit.") ap.add_argument("--eff-rms", type=str, default="0.01,0.02,0.03,0.045,0.06,0.08,0.10,0.13,0.17,0.22,0.28,0.35") ap.add_argument("--eff-targets", type=str, default="0.9,0.8,0.7,0.6") ap.add_argument("--eff-nrep", type=int, default=12) ap.add_argument("--coherence", action="store_true", help="Run the drive-coherence-time sweep and exit.") ap.add_argument("--coh-rms", type=float, default=0.10) ap.add_argument("--coh-taus", type=str, default="0.3,0.6,1.0,2.0,3.0,4.5,7.0,11.0,18.0,30.0,60.0") args = ap.parse_args() T = args.T dt = args.dt nrep = int(args.nrep) amps = [float(x) for x in args.amps.split(",")] ratios = [float(x) for x in args.ratios.split(",")] # Output directory if args.outdir.strip(): outdir = Path(args.outdir.strip()) else: run_id = datetime.now().strftime("%Y%m%d_%H%M%S") outdir = Path(f"run_{run_id}") outdir.mkdir(parents=True, exist_ok=True) print_progress(f"\nSaving outputs to: {outdir.resolve()}") print_progress(f"drive_mode={args.drive_mode} refractory={args.refractory:g}\n") # Toy parameters sim_params = dict( a=1.0, b=0.75, thresh=1.0, reset=0.0, burst_amp=2.0, burst_tau=0.18 ) # Baseline to estimate intrinsic cycle frequency f0 from GLOBAL PSD peak (fine) print_progress("Baseline run (no drive) to estimate intrinsic frequency f0...") base = simulate_bursty_reconnection( T=T, dt=dt, seed=args.seed0, drive_type="none", drive_amp=0.0, fd=1.0, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=float(args.noise), drive_mode="additive", # irrelevant when drive_type="none" refractory=float(args.refractory), **sim_params ) spec0 = psd_peak_and_q_global(base.t, base.R, dt=dt) f0 = spec0.f_peak if np.isfinite(spec0.f_peak) and spec0.f_peak > 0 else 1.0 mu0, sd0, cv0 = intermittency_metrics(base.R) kx0 = excess_kurtosis(base.R) print("[baseline]") print(f" mean={mu0:.6g} std={sd0:.6g} CV={cv0:.6g} kurtosis_excess={kx0:.6g}") print(f" PSD GLOBAL peak f0={spec0.f_peak:.6g} Q_global={spec0.q:.6g} true_events={len(base.event_times)}") print(f"\nEstimated intrinsic peak frequency f0 ~ {f0:.6g}\n") if args.coherence: coh = run_coherence( T=T, dt=dt, f0=f0, noise_sigma=float(args.noise), sim_params=sim_params, drive_mode=args.drive_mode, refractory=float(args.refractory), rms=float(args.coh_rms), taus=[float(x) for x in args.coh_taus.split(",")], nrep=int(args.eff_nrep), seed0=int(args.seed0), ) print_coherence(coh) with open(outdir / "coherence.json", "w", encoding="utf-8") as f: json.dump(coh, f, indent=2) print(f"\nSaved: {outdir / 'coherence.json'}") return if args.efficiency: eff = run_efficiency( T=T, dt=dt, f0=f0, noise_sigma=float(args.noise), sim_params=sim_params, drive_mode=args.drive_mode, refractory=float(args.refractory), offfactor=float(args.offfactor), rms_levels=[float(x) for x in args.eff_rms.split(",")], nrep=int(args.eff_nrep), seed0=int(args.seed0), targets=[float(x) for x in args.eff_targets.split(",")], ) print_efficiency(eff) with open(outdir / "efficiency.json", "w", encoding="utf-8") as f: json.dump(eff, f, indent=2) print(f"\nSaved: {outdir / 'efficiency.json'}") return # Main sweep print_progress("Running MAIN sweep (locking maps)...") main_sweep = run_sweep( T=T, dt=dt, f0=f0, amps=amps, ratios=ratios, nrep=nrep, noise_sigma=float(args.noise), sim_params=sim_params, drive_mode=args.drive_mode, refractory=float(args.refractory), progress_prefix="[MAIN] " ) # Robustness sweep print_progress("\nRunning ROBUSTNESS sweep (higher noise)...") robust_sweep = run_sweep( T=T, dt=dt, f0=f0, amps=amps, ratios=ratios, nrep=nrep, noise_sigma=float(args.noise_robust), sim_params=sim_params, drive_mode=args.drive_mode, refractory=float(args.refractory), progress_prefix="[ROBUST] " ) summarize_event_counts("MAIN", main_sweep, amps, ratios) summarize_event_counts("ROBUST", robust_sweep, amps, ratios) # Tongue widths table print("\n[Tongue widths: PLV ≥ 0.8] (using MAIN sweep)") for A, w, rmin, rmax in tongue_widths(main_sweep.PLV, amps, ratios, threshold=0.8): if rmin is None: print(f" A={A:g}: no locking above threshold") else: print(f" A={A:g}: width={w:.3g} band=[{rmin:.3g}, {rmax:.3g}]") # One-off demo comparisons A_demo = float(args.demoA) fd_res = 1.0 * f0 fd_off = float(args.offfactor) * f0 print_progress("\nRunning one-off comparisons (baseline vs resonant vs off-resonant vs broadband)...") res = simulate_bursty_reconnection( T=T, dt=dt, seed=2, drive_type="sine", drive_amp=A_demo, fd=fd_res, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=float(args.noise), drive_mode=args.drive_mode, refractory=float(args.refractory), **sim_params ) off = simulate_bursty_reconnection( T=T, dt=dt, seed=3, drive_type="sine", drive_amp=A_demo, fd=fd_off, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=float(args.noise), drive_mode=args.drive_mode, refractory=float(args.refractory), **sim_params ) bb = simulate_bursty_reconnection( T=T, dt=dt, seed=4, drive_type="broadband", drive_amp=A_demo, fd=f0, drive_phase=0.0, bb_band=(0.1, 5.0), noise_sigma=float(args.noise), drive_mode=args.drive_mode, refractory=float(args.refractory), **sim_params ) def summarize(label: str, out: SimOut, fd_for_phase: Optional[float]) -> None: mu, sd, cv = intermittency_metrics(out.R) kx = excess_kurtosis(out.R) specG = psd_peak_and_q_global(out.t, out.R, dt=dt) if fd_for_phase is not None: plv = phase_locking_value(out.event_times, fd_for_phase) cph = 1.0 - plv if np.isfinite(plv) else np.nan qd = psd_q_at_target(out.t, out.R, dt=dt, f_target=fd_for_phase, window_frac=0.15).q else: plv, cph, qd = np.nan, np.nan, np.nan print(f"[{label}]") print(f" mean={mu:.6g} std={sd:.6g} CV={cv:.6g} kurtosis_excess={kx:.6g}") print(f" PSD GLOBAL peak f={specG.f_peak:.6g} Q_global={specG.q:.6g} true_events={len(out.event_times)}") if np.isfinite(plv): print(f" PLV={plv:.6g} C_phase={cph:.6g} Q_drive={qd:.6g}") summarize("baseline (repeated)", base, None) summarize(f"resonant drive (A={A_demo}, fd≈f0)", res, fd_res) summarize(f"off-resonant drive (A={A_demo}, fd={args.offfactor} f0)", off, fd_off) summarize(f"broadband drive (A={A_demo})", bb, f0) # Candidate A/B illustrative deficits vs ε (still illustrative) eps_list = np.array([(r - 1.0) for r in ratios], dtype=float) n_lat = 64 k_int, xi_int, om_int = 7, 3, 5 Delta0 = np.array([ TWO_PI * k_int / n_lat, TWO_PI * xi_int / n_lat, TWO_PI * om_int / n_lat ], dtype=float) v = np.array([5.0 * (TWO_PI / n_lat), 0.0, 0.0], dtype=float) v_norm = float(np.linalg.norm(v)) eps_boundary = lattice_basin_boundary_eps(n_lat, v_norm) R_geom = 3.0 a_geom = 0.6 tau0 = 0.3 C_lat_list = [] C_geom_list = [] for eps in eps_list: Delta_eps = Delta0 + eps * v C_lat_list.append(closure_deficit_lattice(Delta_eps, n=n_lat, search_radius=1).C_lat) C_geom_list.append(closure_deficit_geom(R=R_geom, a=a_geom, tau0=tau0, eps=eps)) print("\n[Resonant math notes]") print(f" Candidate A illustrative basin boundary estimate: |ε| ~ {eps_boundary:.6g}") print(" Candidate A: C_lat(Δ)=min||nΔ-2πm|| (quantized closure basins in Δ-space)") print(" Candidate B: C_geom(τ0;ε)=2(R+a cosτ0)|sin(π ε)| (observable geometric mismatch)") print(" Operational entrainment closure: C_phase = 1 - PLV (closure on event-phase circle S^1)\n") # Save arrays np.savez( outdir / "sweeps.npz", amps=np.array(amps), ratios=np.array(ratios), f0=np.array([f0]), main_PLV=main_sweep.PLV, main_Cphase=main_sweep.Cphase, main_CV=main_sweep.CV, main_Qdrive=main_sweep.Qdrive, main_EventsMean=main_sweep.EventsMean, main_EventsStd=main_sweep.EventsStd, main_FRI=main_sweep.FRI, robust_PLV=robust_sweep.PLV, robust_Cphase=robust_sweep.Cphase, robust_CV=robust_sweep.CV, robust_Qdrive=robust_sweep.Qdrive, robust_EventsMean=robust_sweep.EventsMean, robust_EventsStd=robust_sweep.EventsStd, robust_FRI=robust_sweep.FRI, eps_list=eps_list, C_lat_list=np.array(C_lat_list), C_geom_list=np.array(C_geom_list), ) np.savez( outdir / "oneoff_timeseries.npz", t_base=base.t, R_base=base.R, ev_base=base.event_times, t_res=res.t, R_res=res.R, ev_res=res.event_times, t_off=off.t, R_off=off.R, ev_off=off.event_times, t_bb=bb.t, R_bb=bb.R, ev_bb=bb.event_times, f0=np.array([f0]), fd_res=np.array([fd_res]), fd_off=np.array([fd_off]), A_demo=np.array([A_demo]), ) if args.no_plots: print_progress(f"Done. Saved arrays to: {outdir.resolve()}") return def plot_heatmap_single(z: np.ndarray, amps: List[float], ratios: List[float], title: str, cbar_label: str, outpath: Path) -> None: fig, ax = plt.subplots(figsize=(7, 4)) im = ax.imshow(z, aspect="auto", origin="lower", extent=[min(ratios), max(ratios), min(amps), max(amps)]) ax.set_xlabel("fd / f0") ax.set_ylabel("drive amplitude A") ax.set_title(title) cbar = fig.colorbar(im, ax=ax) cbar.set_label(cbar_label) fig.tight_layout() fig.savefig(outpath, dpi=160) plt.close(fig) # --------------------- # Plot 1: Heatmaps MAIN # --------------------- fig, ax = plt.subplots(1, 3, figsize=(15, 4), constrained_layout=True) heat(fig, ax[0], main_sweep.PLV, "PLV (higher = stronger locking) [MAIN]", ratios, amps, vmin=0.0, vmax=1.0) heat(fig, ax[1], main_sweep.CV, "CV of R(t) (lower = less intermittent) [MAIN]", ratios, amps) heat(fig, ax[2], main_sweep.Qdrive, "Q_drive (near fd; higher = sharper line) [MAIN]", ratios, amps) fig.savefig(outdir / "heatmaps_main.png", dpi=200) plot_heatmap_single(main_sweep.FRI, amps, ratios, "False Resonance Index FRI = Q_drive*(1-PLV) [MAIN]", "FRI", outdir / "heatmap_fri_main.png") if args.save_validation: # Validation-only: event-count heatmap (disabled by default) plot_heatmap_single(main_sweep.EventsMean, amps, ratios, "Event count (mean per cell) [MAIN]", "events (mean)", outdir / "heatmap_events_main.png") # -------------------------- # Plot 2: Heatmaps ROBUSTNESS # -------------------------- figR, axR = plt.subplots(1, 3, figsize=(15, 4), constrained_layout=True) heat(figR, axR[0], robust_sweep.PLV, "PLV (higher = stronger locking) [ROBUST]", ratios, amps, vmin=0.0, vmax=1.0) heat(figR, axR[1], robust_sweep.CV, "CV of R(t) [ROBUST]", ratios, amps) heat(figR, axR[2], robust_sweep.Qdrive, "Q_drive (near fd) [ROBUST]", ratios, amps) figR.savefig(outdir / "heatmaps_robust.png", dpi=200) plot_heatmap_single(robust_sweep.FRI, amps, ratios, "False Resonance Index FRI = Q_drive*(1-PLV) [ROBUST]", "FRI", outdir / "heatmap_fri_robust.png") if args.save_validation: # Validation-only: event-count heatmap (disabled by default) plot_heatmap_single(robust_sweep.EventsMean, amps, ratios, "Event count (mean per cell) [ROBUST]", "events (mean)", outdir / "heatmap_events_robust.png") # --------------------------------------------------- # Plot 3: Locking curves (PLV vs fd/f0) # --------------------------------------------------- figC, axC = plt.subplots(1, 1, figsize=(10, 4), constrained_layout=True) for ia, A in enumerate(amps): axC.plot(ratios, main_sweep.PLV[ia, :], marker="o", label=f"A={A:g}") axC.set_xlabel("fd / f0") axC.set_ylabel("PLV") axC.set_title("Locking curves: PLV vs fd/f0 (higher = stronger locking)") axC.set_ylim(0.0, 1.02) axC.legend() figC.savefig(outdir / "locking_curves_plv.png", dpi=200) # --------------------------------------------------------- # Plot 4: Intermittency curves (CV vs fd/f0) # --------------------------------------------------------- figCV, axCV = plt.subplots(1, 1, figsize=(10, 4), constrained_layout=True) for ia, A in enumerate(amps): axCV.plot(ratios, main_sweep.CV[ia, :], marker="s", label=f"A={A:g}") axCV.set_xlabel("fd / f0") axCV.set_ylabel("CV of R(t)") axCV.set_title("Intermittency curves: CV vs fd/f0 (lower = less intermittent)") axCV.legend() figCV.savefig(outdir / "intermittency_curves_cv.png", dpi=200) # -------------------------- # Plot 5: Time series + PSD # -------------------------- fig2, ax2 = plt.subplots(2, 1, figsize=(12, 6), sharex=True, constrained_layout=True) ax2[0].plot(base.t, base.R, label="baseline") ax2[0].plot(res.t, res.R, label="resonant") ax2[0].plot(off.t, off.R, label="off-resonant") ax2[0].plot(bb.t, bb.R, label="broadband", alpha=0.85) ax2[0].set_ylabel("R(t) proxy") ax2[0].legend() def psd(tt, xx): mu = np.mean(xx) fs = 1.0 / dt f, pxx = welch(xx - mu, fs=fs, nperseg=min(4096, len(xx))) return f[1:], pxx[1:] for (tt, xx, lab) in [(base.t, base.R, "baseline"), (res.t, res.R, "resonant"), (off.t, off.R, "off-resonant"), (bb.t, bb.R, "broadband")]: f, pxx = psd(tt, xx) ax2[1].plot(f, pxx, label=lab) ax2[1].axvline(f0, linestyle="--", label="intrinsic f0") ax2[1].set_xlim(0, 5 * f0) ax2[1].set_xlabel("frequency") ax2[1].set_ylabel("PSD") ax2[1].legend() fig2.savefig(outdir / "timeseries_psd.png", dpi=200) # -------------------------- # Plot 6: Candidate A/B vs ε # -------------------------- fig3, ax3 = plt.subplots(1, 1, figsize=(10, 3.6), constrained_layout=True) ax3.plot(eps_list, C_lat_list, marker="o", label="Candidate A: C_lat(ε)") ax3.plot(eps_list, C_geom_list, marker="s", label="Candidate B: C_geom(ε)") ax3.axvline(+eps_boundary, linestyle="--") ax3.axvline(-eps_boundary, linestyle="--") ax3.set_xlabel("ε = fd/f0 - 1 (detuning coordinate)") ax3.set_ylabel("closure deficit (arb. units)") ax3.set_title("Illustrative closure deficits vs detuning ε (Candidate A/B)") ax3.legend() fig3.savefig(outdir / "candidateA_B_deficits.png", dpi=200) print_progress(f"\nSaved plots + arrays to: {outdir.resolve()}\n") plt.show() if __name__ == "__main__": main()