wcm_substrate.py

"""WCM Substrate Mixin — FHN wave physics, domain lifecycle, IPC, sparse frontier.

Provides SubstrateMixin, the chipset physics layer composed into WaveKernel.
Covers the FitzHugh-Nagumo (FHN) excitable-medium simulation, domain/port
topology, phase-windowed inter-process communication, process scheduling,
sparse-frontier optimisation, and attractor persistence.

Architectural role (see docs/DEV_MODE.md §2.1):
  This is pure chipset infrastructure — it steps the equations, manages
  diffusion topology, and delivers wave messages.  All policy decisions
  (what to display, how to route input) are made by WaveIL, not here.
"""

import numpy as np
from typing import Any, Dict, List, Optional, Set, Tuple

from wcm_constants import (
    A_FHN, B_FHN, DIFFUSE, DT, EPS,
    ACTIVE_BLOCK_FLOOR,
    BISTABLE_THRESHOLD,
    BLOCK_TO_RECLAIM_TICKS,
    DECAY_ENERGY_MULT,
    ENERGY_DECAY,
    GLOBAL_SPIKE_DAMP, GLOBAL_SPIKE_THRESH,
    HOMEO_GAIN, HOMEO_TARGET_ENERGY,
    INTEGRATOR_DECAY, INTEGRATOR_THRESH,
    MEMBRANE_EDGE_SCALE,
    PHASE_GATE_WIDTH,
    PORT_CAPACITY_MAX, PORT_CAPACITY_MIN,
    PORT_DEFER_DROP_AGE, PORT_EDGE_SCALE, PORT_ENERGY_CAP,
    PORT_MAX_DELIVER_PER_TICK, PORT_MAX_QUEUE, PORT_MSG_TTL,
    PORT_PHASE_WINDOW, PORT_PHASE_WINDOW_MAX, PORT_PHASE_WINDOW_MIN,
    PROC_BLOCKED, PROC_DECAYING, PROC_INIT, PROC_RECLAIMED, PROC_RUNNING,
    PUMP_ACTIVITY_HI, PUMP_ACTIVITY_LO,
    PUMP_GAIN_BASELINE, PUMP_GAIN_DOWN, PUMP_GAIN_MAX, PUMP_GAIN_MIN,
    PUMP_GAIN_RELAX, PUMP_GAIN_UP,
    PUMP_INJECT_RADIUS, PUMP_PHASE_COUPLING,
    REFRACT_DECAY,
    RES_RING_BOOST, RES_RING_DAMP,
    SUBSTEPS,
    V_GATE, V_HIGH, V_LOW, V_REST,
    W_REST,
    disk_indices,
    AttractorObject,
    Domain,
)


class SubstrateMixin:

	# --- Process scheduler ---------------------------------------------------

	def _scheduler_pass(self, cached_activities=None):
		"""Tick the process state machine (INIT→RUNNING→BLOCKED/DECAYING→RECLAIMED)."""
		self._refresh_process_endpoints()
		for pid in sorted(self.process_table):
			p = self.process_table[pid]
			if p.state == PROC_RECLAIMED:
				continue
			if p.domain_name not in self.domains:
				self._transition_process_state(pid, PROC_RECLAIMED, "domain-missing")
				continue

			d = self.domains[p.domain_name]
			if cached_activities is not None:
				did = self.domain_id.get(p.domain_name)
				activity = float(cached_activities[did]) if did is not None and did < len(cached_activities) else 0.0
			else:
				m = d.mask
				activity = float(np.mean(np.abs(self.v[m] - np.float32(V_REST))))
			energy_ratio = activity / max(1e-6, p.energy_budget)
			has_work = len(d.inbox) > 0 or len(d.outbox) > 0
			has_work = has_work or len(d.user_mailbox) > 0

			p.total_ticks += 1
			p.ticks_in_state += 1

			if p.state == PROC_INIT:
				self._transition_process_state(pid, PROC_RUNNING, "init-complete")
				continue

			if p.state == PROC_RUNNING:
				if energy_ratio > DECAY_ENERGY_MULT:
					self._transition_process_state(pid, PROC_DECAYING, "over-budget")
				elif (activity < ACTIVE_BLOCK_FLOOR) and (not has_work):
					self._transition_process_state(pid, PROC_BLOCKED, "quiescent")
				continue

			if p.state == PROC_BLOCKED:
				if p.block_reason == "recv-wait":
					timed_out = p.recv_timeout_tick >= 0 and self.tick_count >= p.recv_timeout_tick
					if len(d.user_mailbox) > 0:
						self._transition_process_state(pid, PROC_RUNNING, "wakeup")
					elif timed_out:
						self._transition_process_state(pid, PROC_RUNNING, "recv-timeout")
				elif has_work or activity >= ACTIVE_BLOCK_FLOOR:
					self._transition_process_state(pid, PROC_RUNNING, "wakeup")
				elif p.ticks_in_state > BLOCK_TO_RECLAIM_TICKS:
					self._transition_process_state(pid, PROC_RECLAIMED, "idle-timeout")
				continue

			if p.state == PROC_DECAYING:
				self.v[m] = np.float32(V_REST) + (self.v[m] - np.float32(V_REST)) * np.float32(0.975)
				if energy_ratio <= 1.05:
					self._transition_process_state(pid, PROC_RUNNING, "budget-recovered")

	# --- Domain & port setup --------------------------------------------------

	def add_domain(self, name: str, mask: np.ndarray, phase_freq: float, phase_offset: float = 0.0):
		"""Register a named region on the grid.  Raises on overlap."""
		if name in self.domain_id:
			raise ValueError(f"Duplicate domain: {name}")
		did = len(self.domain_id)
		self.domain_id[name] = did

		overlap = mask & (self.owner >= 0)
		if np.any(overlap):
			raise ValueError(f"Domain {name} overlaps an existing domain")
		self.owner[mask] = did
		self.domain_union_mask |= mask
		self.edge_scale_dirty = True

		self.domains[name] = Domain(
			name=name,
			mask=mask,
			phase_freq=phase_freq,
			phase=0.0,
			phase_offset=phase_offset,
			energy_budget=1.0,
		)

	def add_port(self, domain_a: str, domain_b: str, mask: np.ndarray, strength: float = PORT_EDGE_SCALE):
		"""Create an IPC channel between two domains with its own diffusion corridor."""
		if domain_a not in self.domain_id or domain_b not in self.domain_id:
			raise ValueError("Both domains must exist before adding a port")
		self.ports.append({
			"id": self.next_port_id,
			"a": domain_a,
			"b": domain_b,
			"mask": mask.astype(bool),
			"strength": float(strength),
			"queue": [],
			"phase_window": float(PORT_PHASE_WINDOW),
			"base_phase_window": float(PORT_PHASE_WINDOW),
			"max_queue": int(PORT_MAX_QUEUE),
			"capacity": int(PORT_MAX_DELIVER_PER_TICK),
			"base_capacity": int(PORT_MAX_DELIVER_PER_TICK),
			"capacity_min": int(PORT_CAPACITY_MIN),
			"capacity_max": int(PORT_CAPACITY_MAX),
			"energy_cap": float(PORT_ENERGY_CAP),
			"stats": {
				"enqueued": 0,
				"delivered": 0,
				"deferred": 0,
				"dropped": 0,
			},
			"recent": {
				"enqueued": [],
				"delivered": [],
				"deferred": [],
				"dropped": [],
			},
		})
		self.next_port_id += 1
		self.edge_scale_dirty = True

	def _adapt_port_controls(self, port: dict):
		"""Adaptive flow control: tune delivery capacity and phase window."""
		q_len = len(port["queue"])
		q_frac = q_len / max(1, int(port["max_queue"]))
		pe = self._port_energy(port)
		cap = int(port["capacity"])
		base_cap = int(port["base_capacity"])
		cap_min = int(port["capacity_min"])
		cap_max = int(port["capacity_max"])
		pw = float(port["phase_window"])
		base_pw = float(port["base_phase_window"])
		ecap = float(port["energy_cap"])

		# Congested port: shrink openness and throughput.
		if pe > ecap:
			cap = max(cap_min, cap - 1)
			pw = max(PORT_PHASE_WINDOW_MIN, pw - 0.01)
		# Under energetic limits but queueing: open wider and deliver more.
		elif q_frac > 0.65 and pe < 0.75 * ecap:
			cap = min(cap_max, cap + 1)
			pw = min(PORT_PHASE_WINDOW_MAX, pw + 0.015)
		# Relax toward baseline when stable.
		else:
			if cap > base_cap:
				cap -= 1
			elif cap < base_cap:
				cap += 1
			if pw > base_pw:
				pw = max(base_pw, pw - 0.005)
			elif pw < base_pw:
				pw = min(base_pw, pw + 0.005)

		port["capacity"] = cap
		port["phase_window"] = pw

	def _phase_open(self, domain_name: str, width: float) -> bool:
		d = self.domains[domain_name]
		phase = (d.phase + d.phase_offset) % 1.0
		return phase < width or phase > (1.0 - width)

	def _port_energy(self, port: dict) -> float:
		m = port["mask"]
		if not np.any(m):
			return 0.0
		return float(np.mean(np.abs(self.v[m] - np.float32(V_REST))))

	def _port_center_for_domain(self, port: dict, domain_name: str) -> Tuple[int, int]:
		dm = self.domains[domain_name].mask & port["mask"]
		rr, cc = np.where(dm)
		if rr.size == 0:
			rr, cc = np.where(port["mask"])
			if rr.size == 0:
				return self.rows // 2, self.cols // 2
		return int(np.round(rr.mean())), int(np.round(cc.mean()))

	# --- IPC / channel transfer -----------------------------------------------

	def send_wave(self, src_domain: str, dst_domain: str, amp: float, radius: int = 3) -> bool:
		"""Enqueue a wave message on a matching explicit port."""
		for p in self.ports:
			a = p["a"]
			b = p["b"]
			if not ((src_domain == a and dst_domain == b) or (src_domain == b and dst_domain == a)):
				continue
			q = p["queue"]
			stats = p["stats"]
			if len(q) >= p["max_queue"]:
				stats["dropped"] += 1
				return False
			q.append({
				"src": src_domain,
				"dst": dst_domain,
				"amp": float(amp),
				"radius": int(radius),
				"age": 0,
				"ttl": int(PORT_MSG_TTL),
			})
			stats["enqueued"] += 1
			self._record_port_event(p, "enqueued")
			return True
		return False

	def process_channels(self):
		"""Phase-windowed IPC transfer with congestion-aware backpressure."""
		for p in self.ports:
			self._adapt_port_controls(p)
			q = p["queue"]
			stats = p["stats"]
			if not q:
				continue

			# Age and drop expired messages.
			alive = []
			for msg in q:
				msg["age"] += 1
				if msg["age"] <= msg["ttl"]:
					alive.append(msg)
				else:
					stats["dropped"] += 1
					self._record_port_event(p, "dropped")
			q[:] = alive
			if not q:
				continue

			phase_ok = self._phase_open(p["a"], p["phase_window"]) and self._phase_open(p["b"], p["phase_window"])
			congested = self._port_energy(p) > p["energy_cap"]

			if (not phase_ok) or congested:
				deferred_now = min(len(q), p["capacity"])
				stats["deferred"] += deferred_now
				for _ in range(deferred_now):
					self._record_port_event(p, "deferred")
				if q[0]["age"] > PORT_DEFER_DROP_AGE:
					q.pop(0)
					stats["dropped"] += 1
					self._record_port_event(p, "dropped")
				continue

			n = min(len(q), p["capacity"])
			for _ in range(n):
				msg = q.pop(0)
				dst = msg["dst"]
				if dst not in self.domains:
					stats["dropped"] += 1
					self._record_port_event(p, "dropped")
					continue
				r, c = self._port_center_for_domain(p, dst)
				payload = {
					"kind": "emit",
					"r": r,
					"c": c,
					"amp": msg["amp"],
					"radius": msg["radius"],
					"src": msg["src"],
					"dst": msg["dst"],
				}
				self.domains[dst].inbox.append(payload)
				# User-space mailbox keeps a copy so recv can observe channel traffic
				# even when physical inbox gets consumed by the kernel mail processor.
				self.domains[dst].user_mailbox.append(dict(payload))
				stats["delivered"] += 1
				self._record_port_event(p, "delivered")

	def port_stats_summary(self) -> List[str]:
		lines = []
		for p in self.ports:
			s = p["stats"]
			lines.append(
				f"{p['a']}<->{p['b']}: q={len(p['queue'])} "
				f"enq={s['enqueued']} del={s['delivered']} def={s['deferred']} drop={s['dropped']} "
				f"cap={p['capacity']} pw={p['phase_window']:.2f} pe={self._port_energy(p):.3f}"
			)
		return lines

	# --- Edge topology --------------------------------------------------------

	def rebuild_edge_scales(self):
		"""Rebuild diffusion membrane: low coupling at domain borders, open at port corridors."""
		self.edge_scale[:] = np.float32(1.0)
		self.transport_mask[:] = False
		o = self.owner

		down_cross = (o[:-1, :] != o[1:, :]) & ((o[:-1, :] >= 0) | (o[1:, :] >= 0))
		right_cross = (o[:, :-1] != o[:, 1:]) & ((o[:, :-1] >= 0) | (o[:, 1:] >= 0))

		self.edge_scale[0, :-1, :][down_cross] = np.float32(MEMBRANE_EDGE_SCALE)
		self.edge_scale[1, 1:, :][down_cross] = np.float32(MEMBRANE_EDGE_SCALE)
		self.edge_scale[2, :, :-1][right_cross] = np.float32(MEMBRANE_EDGE_SCALE)
		self.edge_scale[3, :, 1:][right_cross] = np.float32(MEMBRANE_EDGE_SCALE)

		for p in self.ports:
			ida = self.domain_id[p["a"]]
			idb = self.domain_id[p["b"]]
			mask = p["mask"]
			s = np.float32(p["strength"])

			# Transport corridor used for accounting (port plus a small halo).
			corr = np.array(mask, copy=True)
			for _ in range(2):
				expanded = np.array(corr, copy=True)
				expanded[1:] |= corr[:-1]
				expanded[:-1] |= corr[1:]
				expanded[:, 1:] |= corr[:, :-1]
				expanded[:, :-1] |= corr[:, 1:]
				corr = expanded
			self.transport_mask |= corr

			pair_down = (
				(((o[:-1, :] == ida) & (o[1:, :] == idb)) | ((o[:-1, :] == idb) & (o[1:, :] == ida)))
				& down_cross
				& (mask[:-1, :] | mask[1:, :])
			)
			pair_right = (
				(((o[:, :-1] == ida) & (o[:, 1:] == idb)) | ((o[:, :-1] == idb) & (o[:, 1:] == ida)))
				& right_cross
				& (mask[:, :-1] | mask[:, 1:])
			)

			self.edge_scale[0, :-1, :][pair_down] = s
			self.edge_scale[1, 1:, :][pair_down] = s
			self.edge_scale[2, :, :-1][pair_right] = s
			self.edge_scale[3, :, 1:][pair_right] = s
		self.edge_scale_dirty = False

	# --- Emit, route, resonator -----------------------------------------------

	def add_resonator(self, domain_name: str, r: int, c: int, radius: int = 5):
		"""Pin an anchor point that sustains oscillation in a domain."""
		self.domains[domain_name].anchors.append((r, c, radius))

	def emit_packet(self, domain_name: str, r: int, c: int, amp: float, radius: int = 3):
		"""Inject a Gaussian voltage pulse centred at (r, c) within a domain."""
		d = self.domains[domain_name]
		rr, cc = disk_indices(self.rows, self.cols, r, c, radius)
		domain_hit = d.mask[rr, cc]
		rr, cc = rr[domain_hit], cc[domain_hit]
		if rr.size == 0:
			return
		dr = rr - r
		dc = cc - c
		d2 = (dr * dr + dc * dc).astype(np.float32)
		pulse = np.float32(amp) * np.exp(-d2 / np.float32(2.0))
		self.v[rr, cc] += pulse
		new_cells = ~self.active_frontier_mask[rr, cc]
		if np.any(new_cells):
			self.active_frontier_count += int(np.count_nonzero(new_cells))
		self.active_frontier_mask[rr, cc] = True
		# Cache emitted coordinates for fast quiescent damping
		if not hasattr(self, "_emitted_frontier_rows"):
			self._emitted_frontier_rows = []
			self._emitted_frontier_cols = []
		self._emitted_frontier_rows.append(rr)
		self._emitted_frontier_cols.append(cc)
		if self.sparse_materialize_frontier_set:
			for r2, c2 in zip(rr.tolist(), cc.tolist()):
				self.active_frontier.add((int(r2), int(c2)))

	def route_operator(self, src_domain: str, dst_domain: str, strength: float = 0.20):
		"""Bias diffusion from source border toward destination border."""
		a = self.domains[src_domain].mask
		b = self.domains[dst_domain].mask

		a_idx = np.argwhere(a)
		b_idx = np.argwhere(b)
		if a_idx.size == 0 or b_idx.size == 0:
			return

		ar, ac = a_idx.mean(axis=0)
		br, bc = b_idx.mean(axis=0)
		dr = br - ar
		dc = bc - ac

		if abs(dr) >= abs(dc):
			if dr > 0:
				self.d[0, a] = np.float32(DIFFUSE + strength)
			else:
				self.d[1, a] = np.float32(DIFFUSE + strength)
		else:
			if dc > 0:
				self.d[2, a] = np.float32(DIFFUSE + strength)
			else:
				self.d[3, a] = np.float32(DIFFUSE + strength)

	# --- Fused whole-grid operators -------------------------------------------
	def _compute_domain_activities_bulk(self, compute_region_mask=None):
		"""Compute per-domain activity in ONE pass over the grid using bincount."""
		n_domains = len(self.domain_id)
		if n_domains == 0:
			return np.zeros(0, dtype=np.float64)
		valid = self.owner >= 0
		if isinstance(compute_region_mask, np.ndarray) and compute_region_mask.shape == valid.shape:
			valid = valid & compute_region_mask
		owner_flat = self.owner[valid].astype(np.intp)
		abs_dev = np.abs(self.v[valid] - np.float32(V_REST))
		sums = np.bincount(owner_flat, weights=abs_dev, minlength=n_domains)
		counts = np.bincount(owner_flat, minlength=n_domains)
		return sums / np.maximum(counts, 1)

	def _fused_domain_operators(self, runnable_names, compute_region_mask=None, cached_activities=None):
		"""Apply gatekeeper + homeostasis + integrator + metabolic_pump + resonator
		as whole-grid field operations instead of per-domain masked loops.

		Gatekeeper and homeostasis are fused into per-cell scale+bias arrays
		indexed by owner, applied in one vectorized pass.
		Integrator runs as one whole-grid pass.
		Resonator and pump injection stay as sparse local writes (already cheap).
		"""
		n_domains = len(self.domain_id)
		if n_domains == 0:
			return
		# Build runnable lookup
		runnable_lut = np.zeros(n_domains, dtype=bool)
		for name in runnable_names:
			did = self.domain_id.get(name)
			if did is not None:
				runnable_lut[did] = True
		if not np.any(runnable_lut):
			return

		# ── 1. Bulk activity (reuse cached if available) ──
		activities = cached_activities if cached_activities is not None else self._compute_domain_activities_bulk(compute_region_mask)

		# ── 2. Build per-domain scale and bias (scalar Python loop, ~8 iters) ──
		# Non-runnable domains keep scale=1, bias=0 (identity transform)
		scale = np.ones(n_domains, dtype=np.float32)
		bias = np.zeros(n_domains, dtype=np.float32)

		for name in runnable_names:
			did = self.domain_id.get(name)
			if did is None:
				continue
			d = self.domains[name]
			activity = float(activities[did])

			# — Gatekeeper: closed-phase domains get damped —
			phase = (d.phase + d.phase_offset) % 1.0
			is_open = phase < PHASE_GATE_WIDTH or phase > (1.0 - PHASE_GATE_WIDTH)
			if not is_open:
				scale[did] *= np.float32(0.985)

			# — Homeostasis: energy stabilization —
			budget_target = float(HOMEO_TARGET_ENERGY * max(0.1, d.energy_budget))
			e = activity + abs(float(V_REST))
			delta = e - budget_target
			if delta > 0:
				scale[did] *= np.float32(max(0.90, 1.0 - HOMEO_GAIN * delta))
			else:
				bias[did] += np.float32(HOMEO_GAIN * (-delta) * 0.02)

			# — Metabolic pump gain adaptation (scalar, no grid access) —
			if activity < PUMP_ACTIVITY_LO:
				d.pump_gain = min(PUMP_GAIN_MAX, d.pump_gain + PUMP_GAIN_UP)
			elif activity > PUMP_ACTIVITY_HI:
				d.pump_gain = max(PUMP_GAIN_MIN, d.pump_gain - PUMP_GAIN_DOWN)
			else:
				if d.pump_gain > PUMP_GAIN_BASELINE:
					d.pump_gain = max(PUMP_GAIN_BASELINE, d.pump_gain - PUMP_GAIN_RELAX)
				else:
					d.pump_gain = min(PUMP_GAIN_BASELINE, d.pump_gain + PUMP_GAIN_RELAX)

		# ── 3. Apply scale+bias per domain using cached masks (no temp arrays) ──
		has_crm = isinstance(compute_region_mask, np.ndarray) and compute_region_mask.shape == self.v.shape
		for name in runnable_names:
			did = self.domain_id.get(name)
			if did is None:
				continue
			s, b = scale[did], bias[did]
			if s == 1.0 and b == 0.0:
				continue  # identity — skip entirely
			d = self.domains[name]
			m = d.mask
			if has_crm:
				m = m & compute_region_mask
			if b == 0.0:
				self.v[m] *= s
			else:
				self.v[m] = self.v[m] * s + b

		# ── 4. Integrator: whole-grid pass over owned cells ──
		union = self.domain_union_mask
		imask = union & compute_region_mask if has_crm else union
		self.integrator[imask] = (
			self.integrator[imask] * np.float32(INTEGRATOR_DECAY)
			+ np.abs(self.v[imask]) * np.float32(0.01)
		)
		fired = imask & (self.integrator > INTEGRATOR_THRESH)
		if np.any(fired):
			if self.sparse_eventized_operators:
				id_to_name = {v: k for k, v in self.domain_id.items()}
				fired_owners = self.owner[fired]
				for did in np.unique(fired_owners):
					did_int = int(did)
					name = id_to_name.get(did_int)
					if name is None:
						continue
					dom_fired = fired & (self.owner == did_int)
					rr, cc = np.where(dom_fired)
					if rr.size:
						self._queue_sparse_event(name, int(np.mean(rr)), int(np.mean(cc)), amp=0.45, radius=2)
			else:
				self.v[fired] += np.float32(0.45)
			self.integrator[fired] *= np.float32(0.2)

		# ── 5. Sparse local operators (resonator + pump injection) ──
		for name in runnable_names:
			d = self.domains.get(name)
			if d is None or not d.anchors:
				continue
			eff_mask = d.mask
			if isinstance(compute_region_mask, np.ndarray) and compute_region_mask.shape == d.mask.shape:
				eff_mask = d.mask & compute_region_mask
				if not np.any(eff_mask):
					continue
			# Resonator
			if self.sparse_eventized_operators:
				for (r, c, rad) in d.anchors:
					self._queue_sparse_event(name, int(r), int(c),
						amp=float(RES_RING_BOOST) * 2.0,
						radius=max(2, int(rad) // 3), mask_override=eff_mask)
			else:
				for (r, c, rad) in d.anchors:
					rr, cc = disk_indices(self.rows, self.cols, r, c, rad)
					in_dom = eff_mask[rr, cc]
					rr, cc = rr[in_dom], cc[in_dom]
					if rr.size == 0:
						continue
					ring = np.abs((rr - r) ** 2 + (cc - c) ** 2 - rad * rad) <= (rad + 1)
					rr2, cc2 = rr[ring], cc[ring]
					if rr2.size:
						self.v[rr2, cc2] += np.float32(RES_RING_BOOST)
						self.v[rr2, cc2] *= np.float32(RES_RING_DAMP)
			# Metabolic pump injection
			if d.pump_gain < 1e-4:
				continue
			phase = (d.phase + d.phase_offset) % 1.0
			phase_mod = 0.5 + PUMP_PHASE_COUPLING * np.cos(2.0 * np.pi * phase)
			amp = d.pump_gain * phase_mod
			for (r, c, rad) in d.anchors:
				rr, cc = disk_indices(self.rows, self.cols, r, c, PUMP_INJECT_RADIUS)
				in_dom = eff_mask[rr, cc]
				rr, cc = rr[in_dom], cc[in_dom]
				if rr.size == 0:
					continue
				dr = rr - r
				dc = cc - c
				d2 = (dr * dr + dc * dc).astype(np.float32)
				pulse = np.float32(amp) * np.exp(-d2 / np.float32(3.0))
				self.v[rr, cc] += pulse

	def _fused_global_spike_damp(self, activities=None, compute_region_mask=None):
		"""Global spike damping using pre-computed activities instead of re-reading grid."""
		union = self.domain_union_mask
		if not np.any(union):
			return
		# Use pre-computed activities if available
		if activities is not None and len(activities) > 0:
			# Weighted mean of domain activities
			n_domains = len(self.domain_id)
			valid = self.owner >= 0
			if isinstance(compute_region_mask, np.ndarray):
				valid = valid & compute_region_mask
			counts = np.bincount(self.owner[valid].astype(np.intp), minlength=n_domains)
			total = float(np.sum(counts))
			g_act = float(np.sum(activities * counts)) / max(1.0, total) if total > 0 else 0.0
		else:
			g_act = float(np.mean(np.abs(self.v[union] - np.float32(V_REST))))
		if g_act > GLOBAL_SPIKE_THRESH:
			has_freeze = np.any(self.storage_freeze_mask)
			if has_freeze:
				damp_mask = union & ~self.storage_freeze_mask
			else:
				damp_mask = union
			self.v[damp_mask] *= np.float32(GLOBAL_SPIKE_DAMP)

	# --- Phase, mailbox, sparse events ----------------------------------------

	def update_phases(self):
		"""Advance each domain's oscillation phase by its frequency."""
		for d in self.domains.values():
			d.phase = (d.phase + d.phase_freq) % 1.0

	def process_mailboxes(self):
		"""Consume bounded inbox messages and emit their wave payloads."""
		for d in self.domains.values():
			# Consume bounded inbox amount each tick.
			for _ in range(min(2, len(d.inbox))):
				msg = d.inbox.pop(0)
				if msg["kind"] == "emit":
					self.emit_packet(d.name, msg["r"], msg["c"], msg["amp"], msg.get("radius", 3))

	def _queue_sparse_event(self, domain_name: str, r: int, c: int, amp: float, radius: int = 2, mask_override: Optional[np.ndarray] = None):
		if domain_name not in self.domains:
			return
		d = self.domains[domain_name]
		eff_mask = d.mask
		if isinstance(mask_override, np.ndarray) and mask_override.shape == d.mask.shape:
			eff_mask = d.mask & mask_override
		if not np.any(eff_mask):
			return
		r0 = int(max(0, min(self.rows - 1, r)))
		c0 = int(max(0, min(self.cols - 1, c)))
		if not bool(eff_mask[r0, c0]):
			rr, cc = np.where(eff_mask)
			if rr.size == 0:
				return
			r0 = int(rr[0])
			c0 = int(cc[0])
		self.sparse_event_queue.append((domain_name, r0, c0, float(amp), int(max(1, radius))))

	def _apply_sparse_events(self):
		if not self.sparse_event_queue:
			return
		for domain_name, r0, c0, amp, radius in self.sparse_event_queue:
			self.emit_packet(domain_name, r0, c0, amp=float(amp), radius=int(radius))
		self.sparse_event_queue = []

	def _frontier_mask_to_set(self, mask: np.ndarray) -> Set[Tuple[int, int]]:
		rr, cc = np.where(mask)
		return {(int(r), int(c)) for r, c in zip(rr.tolist(), cc.tolist())}

	def _shift4_or(self, mask: np.ndarray) -> np.ndarray:
		out = np.zeros_like(mask)
		out[1:, :] |= mask[:-1, :]
		out[:-1, :] |= mask[1:, :]
		out[:, 1:] |= mask[:, :-1]
		out[:, :-1] |= mask[:, 1:]
		return out

	def apply_sparse_control_packet(self, packet: Dict[str, Any]) -> bool:
		"""Control-plane entrypoint for sparse policy and event steering."""
		if not isinstance(packet, dict):
			return False

		policy = packet.get("policy", None)
		if isinstance(policy, dict):
			enabled = bool(policy.get("enabled", self.substrate_hal_name == "sparse-numpy"))
			self.set_sparse_substrate_mode(
				enabled=enabled,
				active_eps=float(policy.get("active_eps", self.sparse_active_eps)),
				delta_only_frontier=bool(policy.get("delta_only_frontier", self.sparse_delta_only_frontier)),
				eventized_operators=bool(policy.get("eventized_operators", self.sparse_eventized_operators)),
				frontier_budget_ratio=float(policy.get("frontier_budget_ratio", self.sparse_frontier_budget_ratio)),
				passive_relax_interval=int(policy.get("passive_relax_interval", self.sparse_passive_relax_interval)),
			)

		events = packet.get("events", [])
		if isinstance(events, list):
			for evt in events:
				if not isinstance(evt, dict):
					continue
				domain = str(evt.get("domain", ""))
				if not domain:
					continue
				self._queue_sparse_event(
					domain_name=domain,
					r=int(evt.get("r", 0)),
					c=int(evt.get("c", 0)),
					amp=float(evt.get("amp", 0.8)),
					radius=int(evt.get("radius", 2)),
				)

		if bool(packet.get("commit_events", False)):
			self._apply_sparse_events()
		return True

	def _rebuild_active_frontier_from_state(self):
		delta_v = np.abs(self.v - np.float32(V_REST))
		delta_w = np.abs(self.w - np.float32(W_REST))
		delta = delta_v + delta_w
		mask = delta > np.float32(max(1e-6, float(self.sparse_active_eps)))
		self.active_frontier_mask = mask
		self.active_frontier_count = int(np.count_nonzero(mask))
		if self.sparse_materialize_frontier_set:
			rr, cc = np.where(mask)
			self.active_frontier = {(int(r), int(c)) for r, c in zip(rr.tolist(), cc.tolist())}
		else:
			self.active_frontier = set()

	def _damp_emitted_frontier(self):
		"""Damp recently emitted frontier cells toward rest and clear frontier."""
		if self._emitted_frontier_rows:
			fr_r = np.concatenate(self._emitted_frontier_rows)
			fr_c = np.concatenate(self._emitted_frontier_cols)
			self.v[fr_r, fr_c] = np.float32(V_REST) + (self.v[fr_r, fr_c] - np.float32(V_REST)) * np.float32(0.5)
			self.w[fr_r, fr_c] = np.float32(W_REST) + (self.w[fr_r, fr_c] - np.float32(W_REST)) * np.float32(0.5)
			self.active_frontier_mask[fr_r, fr_c] = False
		else:
			fr = np.where(self.active_frontier_mask)
			if fr[0].size:
				self.v[fr] = np.float32(V_REST) + (self.v[fr] - np.float32(V_REST)) * np.float32(0.5)
				self.w[fr] = np.float32(W_REST) + (self.w[fr] - np.float32(W_REST)) * np.float32(0.5)
				self.active_frontier_mask[fr] = False
		self.active_frontier_count = 0

	# --- Sparse-hybrid FHN stepper --------------------------------------------

	def _kernel_step_numpy_sparse_hybrid(self, ticks: int = 1):
		"""Main physics tick loop.

		Three tiers: quiescent (skip grid), near-quiescent (damp only),
		active (frontier + halo FHN integration with budget cap).
		"""
		for _ in range(ticks):
			self.substrate_step_stats["total_ticks"] = float(self.substrate_step_stats.get("total_ticks", 0.0)) + 1.0
			self.tick_count += 1
			self._emitted_frontier_rows = []
			self._emitted_frontier_cols = []

			# Bulk activity computation — one bincount pass for all domains
			compute_region = str(self.region_runtime_policy.get("compute_region", "")).strip()
			compute_region_mask = self.chip_region_masks.get(compute_region, None) if compute_region else None

			# ── Quiescent fast-path: skip heavy grid ops when frontier is empty ──
			consecutive_idle = int(getattr(self, "_consecutive_idle_ticks", 0))
			is_quiescent = self.active_frontier_count <= 0 and consecutive_idle >= 1
			# Near-quiescent: tiny frontier from mailbox emissions — skip heavy
			# preamble but still run FHN physics on the small frontier.
			near_quiescent = (not is_quiescent
				and self.active_frontier_count > 0
				and self.active_frontier_count < 200
				and consecutive_idle >= 0)

			if is_quiescent or near_quiescent:
				# Re-use cached zero activities — no need to scan 2M cells
				n_domains = len(self.domain_id)
				cached_activities = getattr(self, "_cached_zero_activities", None)
				if cached_activities is None or len(cached_activities) != n_domains:
					cached_activities = np.zeros(n_domains, dtype=np.float64)
					self._cached_zero_activities = cached_activities
			else:
				cached_activities = self._compute_domain_activities_bulk(compute_region_mask)

			self._scheduler_pass(cached_activities=cached_activities)
			self.update_phases()
			self.process_channels()
			self.process_mailboxes()
			if bool(self.edge_scale_dirty):
				self.rebuild_edge_scales()

			# If we entered quiescent, stay on fast path — any cells emitted by
			# process_mailboxes will be picked up next tick (1-tick delay, acceptable).
			if is_quiescent:
				if self.active_frontier_count > 0 and self.active_frontier_count < 200:
					self._damp_emitted_frontier()
				self._emitted_frontier_rows = []
				self._emitted_frontier_cols = []
				self._consecutive_idle_ticks = consecutive_idle + 1
				continue

			if not is_quiescent and not near_quiescent:
				runnable = []
				for p in self.process_table.values():
					if p.state in (PROC_RUNNING, PROC_DECAYING):
						runnable.extend([p.domain_name] * max(1, p.priority))
				if not runnable:
					runnable = list(self.domains.keys())

				# Fused whole-grid operators instead of per-domain loops
				self._fused_domain_operators(runnable, compute_region_mask, cached_activities=cached_activities)
				if self.sparse_eventized_operators and self.substrate_hal_name == "sparse-numpy":
					self._apply_sparse_events()

				# Global spike damp using cached activities
				self._fused_global_spike_damp(activities=cached_activities, compute_region_mask=compute_region_mask)
			elif near_quiescent:
				if self.active_frontier_count > 0:
					self._damp_emitted_frontier()
				self._emitted_frontier_rows = []
				self._emitted_frontier_cols = []
				self._consecutive_idle_ticks = 1
				continue

			union = self.domain_union_mask
			frontier_ratio = float(self.active_frontier_count) / float(max(1, self.rows * self.cols))
			self.substrate_step_stats["frontier_ratio_sum"] = float(self.substrate_step_stats.get("frontier_ratio_sum", 0.0)) + frontier_ratio
			self.substrate_step_stats["frontier_ratio_peak"] = max(float(self.substrate_step_stats.get("frontier_ratio_peak", 0.0)), frontier_ratio)

			if self.active_frontier_count <= 0:
				self._consecutive_idle_ticks = consecutive_idle + 1
				# Skip passive relax once fully settled (v already at V_REST)
				if consecutive_idle < 10:
					outside = ~union
					if np.any(self.storage_freeze_mask):
						outside = outside & ~self.storage_freeze_mask
					self.v[outside] = np.float32(V_REST) + (self.v[outside] - np.float32(V_REST)) * np.float32(0.94)
				continue
			else:
				self._consecutive_idle_ticks = 0

			self.substrate_step_stats["sparse_ticks"] = float(self.substrate_step_stats.get("sparse_ticks", 0.0)) + 1.0

			max_r, max_c = self.rows - 1, self.cols - 1
			eps = np.float32(max(1e-6, float(self.sparse_active_eps)))

			frontier_mask = np.array(self.active_frontier_mask, copy=True)
			halo_mask = self._shift4_or(frontier_mask) & (~frontier_mask)
			compute_mask = frontier_mask | halo_mask
			if isinstance(compute_region_mask, np.ndarray) and compute_region_mask.shape == compute_mask.shape:
				compute_mask &= compute_region_mask
			if np.any(self.storage_freeze_mask):
				compute_mask &= ~self.storage_freeze_mask
			if not np.any(compute_mask):
				self.active_frontier = set()
				self.active_frontier_mask[:] = False
				self.active_frontier_count = 0
				continue
			rows, cols = np.where(compute_mask)
			rows = rows.astype(np.int32, copy=False)
			cols = cols.astype(np.int32, copy=False)
			active_sel = frontier_mask[rows, cols]
			n_compute = int(rows.size)
			n_active = int(np.count_nonzero(active_sel))

			for _ in range(SUBSTEPS):
				v_act = self.v[rows, cols]
				w_act = self.w[rows, cols]

				# Discrete Laplacian: sum of neighbour voltage differences
				# weighted by directional diffusion (d) and membrane scale.
				lap = np.zeros(n_compute, dtype=np.float32)
				nr = np.minimum(rows + 1, max_r)
				lap += self.d[0, rows, cols] * self.edge_scale[0, rows, cols] * (self.v[nr, cols] - v_act)
				nr = np.maximum(rows - 1, 0)
				lap += self.d[1, rows, cols] * self.edge_scale[1, rows, cols] * (self.v[nr, cols] - v_act)
				nc = np.minimum(cols + 1, max_c)
				lap += self.d[2, rows, cols] * self.edge_scale[2, rows, cols] * (self.v[rows, nc] - v_act)
				nc = np.maximum(cols - 1, 0)
				lap += self.d[3, rows, cols] * self.edge_scale[3, rows, cols] * (self.v[rows, nc] - v_act)

				v_nl = np.clip(np.nan_to_num(v_act, nan=np.float32(V_REST), posinf=np.float32(5.0), neginf=np.float32(-5.0)), -5.0, 5.0)
				w_nl = np.clip(np.nan_to_num(w_act, nan=np.float32(W_REST), posinf=np.float32(5.0), neginf=np.float32(-5.0)), -5.0, 5.0)
				lap = np.nan_to_num(lap, nan=np.float32(0.0), posinf=np.float32(0.0), neginf=np.float32(0.0))

				# FHN voltage: cubic nullcline + recovery + spatial coupling
				dv = v_nl - v_nl * v_nl * v_nl / np.float32(3.0) - w_nl + lap
				a_act = self.a_map[rows, cols]
				b_act = self.b_map[rows, cols]
				eps_act = self.eps_map[rows, cols]
				# FHN recovery: slow variable that quenches excitation
				dw = eps_act * (v_nl + a_act - b_act * w_nl)

				self.v[rows, cols] += np.float32(DT) * dv
				self.w[rows, cols] += np.float32(DT) * dw
				self.refractory_debt[rows, cols] *= np.float32(REFRACT_DECAY)
				self.v[rows, cols] *= np.float32(ENERGY_DECAY)
				np.clip(self.v[rows, cols], -5.0, 5.0, out=self.v[rows, cols])
				np.clip(self.w[rows, cols], -5.0, 5.0, out=self.w[rows, cols])

			fr_rows = rows[active_sel]
			fr_cols = cols[active_sel]
			next_mask = np.zeros((self.rows, self.cols), dtype=bool)
			if fr_rows.size:
				e_rest = np.abs(self.v[fr_rows, fr_cols] - np.float32(V_REST)) + np.abs(self.w[fr_rows, fr_cols] - np.float32(W_REST))
				keep = e_rest >= eps
				next_mask[fr_rows[keep], fr_cols[keep]] = True

			active_fired = self.v[fr_rows, fr_cols] > np.float32(V_GATE)
			if np.any(active_fired):
				f_rows = fr_rows[active_fired]
				f_cols = fr_cols[active_fired]
				next_mask[f_rows, f_cols] = True
				if not self.sparse_delta_only_frontier:
					fired_mask = np.zeros((self.rows, self.cols), dtype=bool)
					fired_mask[f_rows, f_cols] = True
					next_mask |= self._shift4_or(fired_mask)

			h_rows = rows[~active_sel]
			h_cols = cols[~active_sel]
			if h_rows.size:
				halo_fired = self.v[h_rows, h_cols] > np.float32(V_GATE)
				if np.any(halo_fired):
					next_mask[h_rows[halo_fired], h_cols[halo_fired]] = True

			# Exclude storage-frozen cells from next frontier
			if np.any(self.storage_freeze_mask):
				next_mask &= ~self.storage_freeze_mask

			# Keep non-compute domain cells softly anchored to rest, but not every
			# tick to avoid turning sparse into near-dense memory traffic.
			if self.sparse_passive_relax_interval > 0 and (self.tick_count % self.sparse_passive_relax_interval == 0):
				passive = np.array(union, copy=True)
				passive[rows, cols] = False
				if np.any(self.storage_freeze_mask):
					passive &= ~self.storage_freeze_mask
				if np.any(passive):
					self.v[passive] = np.float32(V_REST) + (self.v[passive] - np.float32(V_REST)) * np.float32(0.985)
					self.w[passive] = np.float32(W_REST) + (self.w[passive] - np.float32(W_REST)) * np.float32(0.985)

			candidate_count = int(np.count_nonzero(next_mask))
			if isinstance(compute_region_mask, np.ndarray) and compute_region_mask.shape == next_mask.shape:
				next_mask &= compute_region_mask
				candidate_count = int(np.count_nonzero(next_mask))
			admitted_count = candidate_count
			dropped_count = 0
			if self.sparse_frontier_budget_ratio > 0.0 and candidate_count > 0:
				max_cells = int(float(self.rows * self.cols) * float(self.sparse_frontier_budget_ratio))
				if max_cells > 0 and candidate_count > max_cells:
					cand_rows, cand_cols = np.where(next_mask)
					cand_rows = cand_rows.astype(np.int32, copy=False)
					cand_cols = cand_cols.astype(np.int32, copy=False)
					scores = np.abs(self.v[cand_rows, cand_cols] - np.float32(V_REST)) + np.abs(self.w[cand_rows, cand_cols] - np.float32(W_REST))
					keep_idx = np.argpartition(scores, -max_cells)[-max_cells:]
					keep_mask = np.zeros(candidate_count, dtype=bool)
					keep_mask[keep_idx] = True
					drop_mask = ~keep_mask
					if np.any(drop_mask):
						dr = cand_rows[drop_mask]
						dc = cand_cols[drop_mask]
						self.v[dr, dc] = np.float32(V_REST) + (self.v[dr, dc] - np.float32(V_REST)) * np.float32(0.90)
						self.w[dr, dc] = np.float32(W_REST) + (self.w[dr, dc] - np.float32(W_REST)) * np.float32(0.90)
					next_mask[:] = False
					next_mask[cand_rows[keep_mask], cand_cols[keep_mask]] = True
					admitted_count = int(np.count_nonzero(next_mask))
					dropped_count = int(max(0, candidate_count - admitted_count))

			self.substrate_step_stats["frontier_candidates_sum"] = float(self.substrate_step_stats.get("frontier_candidates_sum", 0.0)) + float(candidate_count)
			self.substrate_step_stats["frontier_admitted_sum"] = float(self.substrate_step_stats.get("frontier_admitted_sum", 0.0)) + float(admitted_count)
			self.substrate_step_stats["frontier_dropped_sum"] = float(self.substrate_step_stats.get("frontier_dropped_sum", 0.0)) + float(dropped_count)

			if np.any(next_mask):
				fr_rows, fr_cols = np.where(next_mask)
				fr_rows = fr_rows.astype(np.int32, copy=False)
				fr_cols = fr_cols.astype(np.int32, copy=False)
				outside_mask = ~union[fr_rows, fr_cols]
				if np.any(outside_mask):
					rr = fr_rows[outside_mask]
					cc = fr_cols[outside_mask]
					self.v[rr, cc] = np.float32(V_REST) + (self.v[rr, cc] - np.float32(V_REST)) * np.float32(0.94)

			self.active_frontier_mask = next_mask
			self.active_frontier_count = int(np.count_nonzero(next_mask))
			if self.sparse_materialize_frontier_set:
				self.active_frontier = self._frontier_mask_to_set(next_mask)
			else:
				self.active_frontier = set()

	def kernel_step(self, ticks: int = 1):
		self.substrate_hal.step_kernel(self, ticks=max(1, int(ticks)))

	# --- Attractor persistence ------------------------------------------------

	def save_attractor(self, key: str, domain_name: str):
		"""Persist a domain's voltage state as bistable latches (or dict fallback)."""
		d = self.domains[domain_name]
		# Wave-native: write as bistable latches
		if self.wave_native_storage:
			zone_name = self._ensure_domain_storage_zone(domain_name)
			if zone_name is not None:
				zone = self.param_zones[zone_name]
				n_zone = zone["zone_rows"] * zone["zone_cols"]
				v_domain = self.v[d.mask]
				bits = (v_domain > np.float32(BISTABLE_THRESHOLD)).astype(np.int32)
				if bits.size < n_zone:
					bits = np.pad(bits, (0, n_zone - bits.size), constant_values=0)
				else:
					bits = bits[:n_zone]
				self.save_to_grid(key, bits, zone_name)
		self.store.save(key, self.v[d.mask], d.mask)

	def damp_domain(self, domain_name: str, factor: float = 0.7):
		d = self.domains[domain_name]
		self.v[d.mask] *= np.float32(factor)
		self.w[d.mask] += np.float32(0.05)

	def load_attractor(self, key: str, domain_name: str, gain: float = 0.6):
		"""Restore saved state into a domain by blending latch values at gain."""
		d = self.domains[domain_name]
		# Wave-native: read from grid latches
		if self.wave_native_storage:
			zone_name = f"domain_store_{domain_name}"
			if zone_name in self.param_zones:
				bits = self.load_from_grid(key, zone_name)
				if bits is not None:
					n_domain = int(np.count_nonzero(d.mask))
					bits = bits[:n_domain]
					v_restored = np.where(bits, np.float32(V_HIGH), np.float32(V_LOW))
					self.v[d.mask] = self.v[d.mask] * np.float32(1.0 - gain) + v_restored * np.float32(gain)
					return True
		# Dict fallback
		if not self.store.has(key):
			return False
		tpl = self.store.load(key)
		vals = tpl["v"]
		self.v[d.mask] = self.v[d.mask] * np.float32(1.0 - gain) + vals * np.float32(gain)
		return True

	def domain_energy(self, domain_name: str) -> float:
		d = self.domains[domain_name]
		return float(np.mean(np.abs(self.v[d.mask])))

	def leakage_ratio(self) -> float:
		union = self.domain_union_mask
		# Leakage should measure propagated activity, not resting baseline potential.
		delta = np.abs(self.v - np.float32(V_REST))
		in_dom = float(np.sum(delta[union])) if np.any(union) else 1e-6
		outside_strict = ~(union | self.transport_mask)
		out_dom = float(np.sum(delta[outside_strict])) if np.any(outside_strict) else 0.0
		return out_dom / (in_dom + out_dom + 1e-6)

	# --- Utilities ------------------------------------------------------------

	def render_surface(self, stride: int = 4) -> str:
		"""Simple text renderer for quick state inspection."""
		sample = self.v[::stride, ::stride]
		out = []
		for row in sample:
			chars = []
			for x in row:
				if x > 0.8:
					chars.append("#")
				elif x > 0.1:
					chars.append("+")
				elif x > -0.5:
					chars.append(".")
				else:
					chars.append(" ")
			out.append("".join(chars))
		return "\n".join(out)