#!/usr/bin/env python3 """Phase 4E: auditable engineering visual-input calibration for MaleCNS-LIF.""" from __future__ import annotations import argparse import gc import hashlib import json import math import os import platform import shutil import sys import time from collections import Counter, deque from pathlib import Path from typing import Any import numpy as np import pyarrow as pa import pyarrow.parquet as pq import torch import yaml from scipy import sparse from scipy.sparse import csgraph from brain_runtime import BrainRuntime, InputBatch, InputEvent def sha256(path: Path, chunk: int = 8 << 20) -> str: h = hashlib.sha256() with path.open("rb") as f: while data := f.read(chunk): h.update(data) return h.hexdigest() def write_json(path: Path, value: Any) -> None: path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(value, indent=2, sort_keys=True, default=np_json) + "\n", encoding="utf-8") def np_json(value: Any) -> Any: if isinstance(value, np.ndarray): return value.tolist() if isinstance(value, np.generic): return value.item() if isinstance(value, Path): return str(value) raise TypeError(type(value).__name__) def stable_seed(base: int, text: str) -> int: return base + int.from_bytes(hashlib.sha256(text.encode()).digest()[:4], "little") % 1_000_000 def verify_checksum_file(root: Path, checksum_file: Path) -> dict[str, Any]: rows, failures = 0, [] for line in checksum_file.read_text(encoding="utf-8").splitlines(): if not line.strip(): continue expected, rel = line.split(None, 1) path = root / rel.strip().lstrip("*") rows += 1 if not path.is_file() or sha256(path) != expected: failures.append(str(path)) return {"checksum_file": str(checksum_file), "checked": rows, "failures": failures, "passed": not failures} def immutable_audit(cfg: dict[str, Any]) -> tuple[list[Path], dict[str, str], dict[str, Any]]: p = cfg["paths"] controls = [ Path(p["reference_root"]) / "checksums.sha256", Path(p["runtime_root"]) / "checksums.sha256", Path(p["io_root"]) / "checksums.sha256", Path(p["phase4c_root"]) / "manifests/checksums.sha256", Path(p["phase4d_root"]) / "manifests/checksums.sha256", ] missing = [str(x) for x in controls if not x.is_file()] if missing: raise FileNotFoundError(f"immutable checksum controls missing: {missing}") snapshots = {str(x): sha256(x) for x in controls} checks = [] for x in controls: # Model seals are relative to their model directory; experiment seals # live in a manifests/ subdirectory but list paths relative to the # experiment root. root = x.parent.parent if x.parent.name == "manifests" else x.parent checks.append(verify_checksum_file(root, x)) if not all(x["passed"] for x in checks): raise RuntimeError(f"immutable baseline checksum failure: {checks}") return controls, snapshots, {"controls": checks, "all_passed": True} def load_populations(cfg: dict[str, Any], runtime: BrainRuntime, exp: Path) -> tuple[dict[str, np.ndarray], dict[str, Any]]: root = Path(cfg["paths"]["phase4c_root"]) / "visual-audit" populations: dict[str, np.ndarray] = {} summary: dict[str, Any] = {} for name, spec in cfg["visual_populations"].items(): if spec["source"] == "hierarchy": path = root / f"hierarchy-{spec['value']}.parquet" nodes = pq.read_table(path, columns=["neuron_idx"])["neuron_idx"].to_numpy().astype(np.int32) expression = f"Phase4C audited hierarchy={spec['value']}" else: field = spec["field"] if field not in runtime.table.column_names: raise RuntimeError(f"metadata field absent: {field}") values = np.asarray(runtime.table[field].to_pylist(), dtype=object) nodes = np.flatnonzero(np.isin(values, spec["values"])).astype(np.int32) expression = f"{field} in {spec['values']}" if not len(nodes) or len(np.unique(nodes)) != len(nodes): raise RuntimeError(f"invalid population {name}") populations[name] = nodes path = exp / "populations" / f"{name}.parquet" pq.write_table(pa.table({ "neuron_idx": pa.array(nodes, type=pa.int32()), "body_id": pa.array(runtime.body_id[nodes], type=pa.uint64()), "superclass": pa.array([runtime.superclass[i] for i in nodes], type=pa.string()), "class": pa.array([runtime.neuron_class[i] for i in nodes], type=pa.string()), "type": pa.array([runtime.neuron_type[i] for i in nodes], type=pa.string()), }), path, compression="zstd") summary[name] = {"count": int(len(nodes)), "selection": expression, **spec, "membership": str(path)} write_json(exp / "populations/summary.json", summary) return populations, summary def reverse_bfs(indptr: np.ndarray, indices: np.ndarray, sources: np.ndarray, n: int) -> np.ndarray: dist = np.full(n, -1, dtype=np.int16) dist[sources] = 0 frontier = sources.astype(np.int64, copy=False) depth = 0 while len(frontier): depth += 1 mask = np.zeros(n, dtype=bool) for node in frontier: mask[indices[int(indptr[node]):int(indptr[node + 1])]] = True mask[dist >= 0] = False frontier = np.flatnonzero(mask).astype(np.int64) if len(frontier): dist[frontier] = depth return dist def graph_reachability(cfg: dict[str, Any], runtime: BrainRuntime, populations: dict[str, np.ndarray], exp: Path) -> dict[str, Any]: cache = exp / "reachability/visual-output-reachability.parquet" summary_path = exp / "reachability/summary.json" if cache.is_file() and summary_path.is_file(): return json.loads(summary_path.read_text(encoding="utf-8")) n = runtime.n data = np.ones(len(runtime.indices_np), dtype=np.uint8) graph = sparse.csr_matrix((data, runtime.indices_np.astype(np.int32, copy=False), runtime.indptr_np), shape=(n, n)) reverse = graph.transpose().tocsr() desc = runtime.output_populations["descending"].astype(np.int64) motor = runtime.output_populations["motor_related"].astype(np.int64) desc_dist = reverse_bfs(reverse.indptr, reverse.indices, desc, n) motor_dist = reverse_bfs(reverse.indptr, reverse.indices, motor, n) component_count, labels = csgraph.connected_components(graph, directed=True, connection="strong", return_labels=True) src = np.repeat(np.arange(n, dtype=np.int32), np.diff(runtime.indptr_np)) csrc, cdst = labels[src], labels[runtime.indices_np] keep = csrc != cdst condensation = sparse.csr_matrix((np.ones(np.count_nonzero(keep), np.uint8), (csrc[keep], cdst[keep])), shape=(component_count, component_count)) condensation.sum_duplicates(); condensation.data[:] = 1 indegree = np.bincount(condensation.indices, minlength=component_count).astype(np.int64) queue = deque(np.flatnonzero(indegree == 0).tolist()); topo: list[int] = [] while queue: c = queue.popleft(); topo.append(c) for nxt in condensation.indices[condensation.indptr[c]:condensation.indptr[c + 1]]: indegree[nxt] -= 1 if indegree[nxt] == 0: queue.append(int(nxt)) if len(topo) != component_count: raise RuntimeError("SCC condensation is unexpectedly cyclic") output = np.concatenate([desc, motor]) bits = [0] * component_count for bit, node in enumerate(output): bits[int(labels[node])] |= 1 << bit for c in reversed(topo): value = bits[c] for nxt in condensation.indices[condensation.indptr[c]:condensation.indptr[c + 1]]: value |= bits[int(nxt)] bits[c] = value desc_mask = (1 << len(desc)) - 1 reach_desc_comp = np.fromiter(((x & desc_mask).bit_count() for x in bits), dtype=np.int32, count=component_count) reach_motor_comp = np.fromiter(((x >> len(desc)).bit_count() for x in bits), dtype=np.int32, count=component_count) rows: list[pa.Table] = [] by_layer: dict[str, Any] = {} raw = runtime.raw_np for layer, nodes in populations.items(): starts, stops = runtime.indptr_np[nodes], runtime.indptr_np[nodes + 1] edge_count = int(np.sum(stops - starts, dtype=np.int64)) weight = int(sum(int(raw[int(a):int(b)].sum(dtype=np.int64)) for a, b in zip(starts, stops))) rd, rm = reach_desc_comp[labels[nodes]], reach_motor_comp[labels[nodes]] dd, dm = desc_dist[nodes], motor_dist[nodes] rows.append(pa.table({ "neuron_idx": pa.array(nodes, type=pa.int32()), "body_id": pa.array(runtime.body_id[nodes], type=pa.uint64()), "visual_layer": pa.array([layer] * len(nodes), type=pa.string()), "min_hops_descending": pa.array(dd, type=pa.int16()), "min_hops_motor": pa.array(dm, type=pa.int16()), "reachable_descending_count": pa.array(rd, type=pa.int32()), "reachable_motor_count": pa.array(rm, type=pa.int32()), })) by_layer[layer] = { "count": int(len(nodes)), "outgoing_edges": edge_count, "outgoing_synapse_weight": weight, "can_reach_descending": int(np.count_nonzero(rd)), "can_reach_motor": int(np.count_nonzero(rm)), "min_hops_descending_distribution": dict(sorted(Counter(map(int, dd)).items())), "min_hops_motor_distribution": dict(sorted(Counter(map(int, dm)).items())), "reachable_descending_count_min_max": [int(rd.min(initial=0)), int(rd.max(initial=0))], "reachable_motor_count_min_max": [int(rm.min(initial=0)), int(rm.max(initial=0))], } pq.write_table(pa.concat_tables(rows), cache, compression="zstd") summary = { "method": "exact directed reachability via SCC-condensation DAG bitsets; min hops via reverse multi-source BFS", "thresholding": False, "strongly_connected_components": int(component_count), "descending_count": int(len(desc)), "motor_count": int(len(motor)), "layers": by_layer, "table": str(cache), } write_json(summary_path, summary) return summary def forward_distances(runtime: BrainRuntime, populations: dict[str, np.ndarray], exp: Path) -> dict[str, np.ndarray]: result = {} for name, nodes in populations.items(): path = exp / "reachability" / f"forward-distance-{name}.int16.npy" if path.is_file(): result[name] = np.load(path) else: result[name] = runtime.runtime_module.p4b.directed_distances(runtime.indptr_np, runtime.indices_np, nodes.astype(np.int64), runtime.n) np.save(path, result[name]) return result def subset(nodes: np.ndarray, percent: float, seed: int, offset: int = 0) -> np.ndarray: count = max(1, int(round(len(nodes) * percent / 100.0))) order = np.random.default_rng(seed).permutation(len(nodes)) if offset + count <= len(nodes): chosen = order[offset:offset + count] else: chosen = np.roll(order, -offset)[:count] return np.sort(nodes[chosen].astype(np.int64)) def make_batch(layer: str, nodes: np.ndarray, rate: float, amplitude: float, duration: float, source: str, stimulation_mode: str = "external_synaptic_drive", events: tuple[InputEvent, ...] | None = None) -> InputBatch: if events is None: event = InputEvent( target_neuron_idx=nodes.astype(np.int64), stimulation_mode=stimulation_mode, value_kind="poisson_rate_hz", value=float(rate), start_ms=0.0, duration_ms=float(duration), source_id=source, channel_id=layer, amplitude_mv=float(amplitude), metadata={"visual_layer": layer, "engineering_interface": stimulation_mode == "external_synaptic_drive", "connectome_modified": False, "external_drive_is_connectome_transmission": False}, ) events = (event,) return InputBatch(events=events, batch_id=f"{source}:{layer}", metadata={"visual_layer": layer}) def peak_step_fraction(runtime: BrainRuntime) -> float: events = runtime.read_spikes() if not len(events["step"]): return 0.0 return float(np.bincount(events["step"].astype(np.int64)).max(initial=0) / runtime.n) def run_trial(runtime: BrainRuntime, batch: InputBatch, targets: np.ndarray, duration: float, seed: int, distance: np.ndarray, cfg: dict[str, Any], capture: bool = False) -> tuple[dict[str, Any], dict[str, np.ndarray] | None]: runtime.reset(seed) chunk = min(float(cfg["calibration"]["chunk_ms"]), duration) elapsed, wall, aborted = 0.0, 0.0, False while elapsed < duration - 1e-9: step_ms = min(chunk, duration - elapsed) t0 = time.perf_counter(); runtime.step(batch, step_ms); wall += time.perf_counter() - t0 active = int(torch.count_nonzero(runtime.spike_counts).item()) total = int(runtime.spike_counts.sum().item()) actual_s = runtime.simulation_time_ms / 1000.0 runaway = cfg["calibration"]["runaway"] if (active / runtime.n >= float(runaway["active_fraction_gte"]) or total / (runtime.n * actual_s) >= float(runaway["mean_firing_rate_hz_gte"]) or peak_step_fraction(runtime) >= float(runaway["peak_step_fraction_gte"])): aborted = True break elapsed += step_ms events = runtime.read_spikes(); nodes = events["neuron_idx"] unique = np.unique(nodes) target_mask = np.zeros(runtime.n, bool); target_mask[targets] = True downstream = unique[~target_mask[unique]] if len(unique) else unique d = runtime.read_population("descending", runtime.simulation_time_ms) m = runtime.read_population("motor_related", runtime.simulation_time_ms) valid_depth = distance[unique] if len(unique) else np.empty(0, np.int16) valid_depth = valid_depth[valid_depth >= 0] descending_mask = np.isin(nodes, runtime.output_populations["descending"]) motor_mask = np.isin(nodes, runtime.output_populations["motor_related"]) out_mask = descending_mask | motor_mask out_latency = float(events["timestamp_ms"][out_mask].min()) if np.any(out_mask) else None total = int(len(nodes)); actual_s = runtime.simulation_time_ms / 1000.0 record = { "seed": seed, "requested_biological_ms": duration, "simulated_biological_ms": runtime.simulation_time_ms, "wall_seconds": wall, "aborted_runaway": aborted, "total_spikes": total, "active_neurons": int(len(unique)), "active_fraction": float(len(unique) / runtime.n), "mean_firing_rate_hz": float(total / (runtime.n * actual_s)), "peak_step_fraction": peak_step_fraction(runtime), "target_count": int(len(targets)), "target_active": int(np.count_nonzero(target_mask[unique])) if len(unique) else 0, "non_target_active": int(len(downstream)), "propagation_depth_from_full_layer": int(valid_depth.max(initial=0)), "descending_active": d["active_neurons"], "descending_spikes": d["total_spikes"], "motor_active": m["active_neurons"], "motor_spikes": m["total_spikes"], "output_first_spike_latency_ms": out_latency, "descending_first_spike_latency_ms": float(events["timestamp_ms"][descending_mask].min()) if np.any(descending_mask) else None, "motor_first_spike_latency_ms": float(events["timestamp_ms"][motor_mask].min()) if np.any(motor_mask) else None, "external_drive_events": int(sum(len(x) for x in runtime.drive_nodes)), "numerically_finite": bool(torch.isfinite(runtime.v).all().item() and torch.isfinite(runtime.g).all().item()), } success = (record["descending_active"] + record["motor_active"] > 0 and record["non_target_active"] > 0 and record["propagation_depth_from_full_layer"] >= int(cfg["calibration"]["propagation_success"]["minimum_graph_hops"])) record["multi_hop_output_propagation"] = bool(success) captured = None if capture: captured = { "step": events["step"].copy(), "neuron_idx": nodes.copy(), "counts": runtime.spike_counts.cpu().numpy().copy(), "v": runtime.v.cpu().numpy().copy(), "g": runtime.g.cpu().numpy().copy(), "feature": np.concatenate([d["spike_counts"], m["spike_counts"]]).astype(np.float32), } return record, captured def save_trial(exp: Path, category: str, key: str, record: dict[str, Any]) -> None: write_json(exp / "calibration" / category / f"{key}.json", record) def calibration_sweeps(cfg: dict[str, Any], runtime: BrainRuntime, populations: dict[str, np.ndarray], distances: dict[str, np.ndarray], exp: Path) -> tuple[dict[str, Any], list[dict[str, Any]]]: cal = cfg["calibration"]; base = float(cfg["model"]["w_syn_mv_per_anatomical_synapse"]) layers = cfg["candidate_selection"]["hierarchy_order"] records: list[dict[str, Any]] = [] amplitude_choice: dict[str, float | None] = {} comparison = [] for i, layer in enumerate(populations): nodes = subset(populations[layer], float(cal["fixed_sparsity_percent"]), stable_seed(int(cal["seed"]), layer)) mode = cfg["visual_populations"][layer]["stimulation_mode"] amplitude = 0.0 if mode == "spike_injection" else 8.0 batch = make_batch(layer, nodes, float(cal["fixed_rate_hz"]), amplitude, float(cal["trial_duration_ms"]), "layer-comparison", mode) rec, _ = run_trial(runtime, batch, nodes, float(cal["trial_duration_ms"]), int(cal["seed"]) + i, distances[layer], cfg) rec.update({"sweep": "five-layer-comparison", "layer": layer, "rate_hz": cal["fixed_rate_hz"], "sparsity_percent": cal["fixed_sparsity_percent"], "amplitude_mv": amplitude, "amplitude_multiple_w_syn": amplitude / base if amplitude else None}) comparison.append(rec); records.append(rec); save_trial(exp, "five-layer-comparison", layer, rec) for li, layer in enumerate(layers): seed0 = stable_seed(int(cal["seed"]), layer) nodes = subset(populations[layer], float(cal["fixed_sparsity_percent"]), seed0) layer_records = [] multipliers = list(map(float, cal["amplitude_multiples"])) found = False for multiplier in multipliers: amp = base * multiplier batch = make_batch(layer, nodes, float(cal["fixed_rate_hz"]), amp, float(cal["trial_duration_ms"]), "amplitude-sweep") rec, _ = run_trial(runtime, batch, nodes, float(cal["trial_duration_ms"]), seed0 + int(multiplier * 10), distances[layer], cfg) rec.update({"sweep": "amplitude", "layer": layer, "rate_hz": cal["fixed_rate_hz"], "sparsity_percent": cal["fixed_sparsity_percent"], "amplitude_mv": amp, "amplitude_multiple_w_syn": multiplier}) layer_records.append(rec); records.append(rec); save_trial(exp, "amplitude", f"{layer}-{multiplier:g}x", rec) found |= rec["multi_hop_output_propagation"] and not rec["aborted_runaway"] if not found: for multiplier in map(float, cal["extended_amplitude_multiples"]): amp = base * multiplier batch = make_batch(layer, nodes, float(cal["fixed_rate_hz"]), amp, float(cal["trial_duration_ms"]), "amplitude-sweep-extended") rec, _ = run_trial(runtime, batch, nodes, float(cal["trial_duration_ms"]), seed0 + int(multiplier * 10), distances[layer], cfg) rec.update({"sweep": "amplitude-extended", "layer": layer, "rate_hz": cal["fixed_rate_hz"], "sparsity_percent": cal["fixed_sparsity_percent"], "amplitude_mv": amp, "amplitude_multiple_w_syn": multiplier}) layer_records.append(rec); records.append(rec); save_trial(exp, "amplitude", f"{layer}-{multiplier:g}x", rec) if rec["multi_hop_output_propagation"] and not rec["aborted_runaway"]: break viable = [x for x in layer_records if x["multi_hop_output_propagation"] and not x["aborted_runaway"] and x["numerically_finite"]] amplitude_choice[layer] = min((x["amplitude_multiple_w_syn"] for x in viable), default=None) for layer in layers: multiplier = amplitude_choice[layer] if multiplier is None: continue seed0 = stable_seed(int(cal["seed"]), layer) fixed_nodes = subset(populations[layer], float(cal["fixed_sparsity_percent"]), seed0) for rate in map(float, cal["rate_hz"]): batch = make_batch(layer, fixed_nodes, rate, base * multiplier, float(cal["trial_duration_ms"]), "rate-sweep") rec, _ = run_trial(runtime, batch, fixed_nodes, float(cal["trial_duration_ms"]), seed0 + int(rate), distances[layer], cfg) rec.update({"sweep": "rate", "layer": layer, "rate_hz": rate, "sparsity_percent": cal["fixed_sparsity_percent"], "amplitude_mv": base * multiplier, "amplitude_multiple_w_syn": multiplier}) records.append(rec); save_trial(exp, "rate", f"{layer}-{rate:g}hz", rec) for sp in map(float, cal["sparsity_percent"]): nodes = subset(populations[layer], sp, seed0) batch = make_batch(layer, nodes, float(cal["fixed_rate_hz"]), base * multiplier, float(cal["trial_duration_ms"]), "sparsity-sweep") rec, _ = run_trial(runtime, batch, nodes, float(cal["trial_duration_ms"]), seed0 + int(sp * 100), distances[layer], cfg) rec.update({"sweep": "sparsity", "layer": layer, "rate_hz": cal["fixed_rate_hz"], "sparsity_percent": sp, "amplitude_mv": base * multiplier, "amplitude_multiple_w_syn": multiplier}) records.append(rec); save_trial(exp, "sparsity", f"{layer}-{sp:g}pct", rec) viable_by_layer: dict[str, dict[str, Any]] = {} for layer in layers: candidates = [x for x in records if x.get("layer") == layer and x["multi_hop_output_propagation"] and not x["aborted_runaway"]] if candidates: candidates.sort(key=lambda x: (float(x.get("amplitude_multiple_w_syn") or 1e9), float(x.get("rate_hz", 1e9)), float(x.get("sparsity_percent", 1e9)))) viable_by_layer[layer] = candidates[0] for layer, chosen in viable_by_layer.items(): seed0 = stable_seed(int(cal["seed"]), layer) for window in map(float, cal["latency_windows_ms"]): nodes = subset(populations[layer], float(chosen["sparsity_percent"]), seed0) batch = make_batch(layer, nodes, float(chosen["rate_hz"]), float(chosen["amplitude_mv"]), window, "latency-sweep") rec, _ = run_trial(runtime, batch, nodes, window, seed0 + int(window * 10), distances[layer], cfg) rec.update({"sweep": "latency", "layer": layer, "rate_hz": chosen["rate_hz"], "sparsity_percent": chosen["sparsity_percent"], "amplitude_mv": chosen["amplitude_mv"], "amplitude_multiple_w_syn": chosen["amplitude_multiple_w_syn"], "window_ms": window}) records.append(rec); save_trial(exp, "latency", f"{layer}-{window:g}ms", rec) table = pa.Table.from_pylist(records) pq.write_table(table, exp / "calibration/all-trials.parquet", compression="zstd") summary = {"base_w_syn_mv": base, "phase4d_external_drive_mv": 8.0, "phase4d_external_drive_multiple_w_syn": 8.0 / base, "amplitude_choice_by_layer": amplitude_choice, "viable_condition_by_layer": viable_by_layer, "trial_count": len(records), "five_layer_comparison": comparison} write_json(exp / "calibration/summary.json", summary) return summary, records def pair_metric(a: np.ndarray, b: np.ndarray) -> dict[str, float]: an, bn = float(np.linalg.norm(a)), float(np.linalg.norm(b)) cosine = float(np.dot(a, b) / (an * bn)) if an and bn else (1.0 if an == bn else 0.0) active_a, active_b = a > 0, b > 0; union = int(np.count_nonzero(active_a | active_b)) return {"euclidean": float(np.linalg.norm(a - b)), "cosine_distance": 1.0 - cosine, "active_jaccard": float(np.count_nonzero(active_a & active_b) / union) if union else 1.0} def pattern_batch(layer: str, populations: dict[str, np.ndarray], kind: str, rate: float, amp: float, duration: float, seed: int, cfg: dict[str, Any], sparse_percent: float | None = None) -> tuple[InputBatch, np.ndarray]: all_nodes = populations[layer] sparse_pct = float(cfg["patterns"]["sparse_percent"] if sparse_percent is None else sparse_percent) count = max(1, int(round(len(all_nodes) * sparse_pct / 100))) order = np.random.default_rng(seed).permutation(len(all_nodes)); sparse_a = np.sort(all_nodes[order[:count]].astype(np.int64)) sparse_b = np.sort(all_nodes[order[count:2 * count]].astype(np.int64)) if kind == "sparse-A": nodes = sparse_a elif kind == "sparse-B": nodes = sparse_b elif kind == "dense": nodes = subset(all_nodes, float(cfg["patterns"]["dense_percent"]), seed) elif kind == "zero": nodes = np.empty(0, np.int64) elif kind == "moving": segs = int(cfg["patterns"]["moving_segments"]); segment_ms = duration / segs; events = [] union = [] for i in range(segs): nodes_i = np.sort(all_nodes[order[(2 + i) * count:(3 + i) * count]].astype(np.int64)); union.append(nodes_i) events.append(InputEvent(target_neuron_idx=nodes_i, stimulation_mode="external_synaptic_drive", value_kind="poisson_rate_hz", value=rate, start_ms=i * segment_ms, duration_ms=segment_ms, source_id="pattern", channel_id=f"moving-{i}", amplitude_mv=amp, metadata={"pattern": "moving", "segment": i})) nodes = np.unique(np.concatenate(union)); return make_batch(layer, nodes, rate, amp, duration, "pattern", events=tuple(events)), nodes else: raise KeyError(kind) return make_batch(layer, nodes, rate, amp, duration, "pattern"), nodes def separability_and_determinism(cfg: dict[str, Any], runtime: BrainRuntime, populations: dict[str, np.ndarray], distances: dict[str, np.ndarray], viable: dict[str, Any], exp: Path) -> dict[str, Any]: result: dict[str, Any] = {} duration = float(cfg["patterns"]["trial_duration_ms"]) for layer, condition in viable.items(): features: dict[str, list[np.ndarray]] = {k: [] for k in ("sparse-A", "sparse-B", "dense", "moving")} trials, drive_counts = [], {} for kind in features: for seed in map(int, cfg["patterns"]["seeds"]): batch, nodes = pattern_batch(layer, populations, kind, float(condition["rate_hz"]), float(condition["amplitude_mv"]), duration, seed, cfg, float(condition.get("sparsity_percent", cfg["patterns"]["sparse_percent"]))) rec, cap = run_trial(runtime, batch, nodes, duration, seed, distances[layer], cfg, capture=True) assert cap is not None; features[kind].append(cap["feature"]); trials.append({"pattern": kind, **rec}) drive_counts[f"{kind}-{seed}"] = rec["external_drive_events"] means = {k: np.mean(v, axis=0) for k, v in features.items()} pairs = {}; names = sorted(means) for i, a in enumerate(names): for b in names[i + 1:]: pairs[f"{a}__{b}"] = pair_metric(means[a], means[b]) det_seed = int(cfg["patterns"]["seeds"][0]) batch, nodes = pattern_batch(layer, populations, "sparse-A", float(condition["rate_hz"]), float(condition["amplitude_mv"]), duration, det_seed, cfg, float(condition.get("sparsity_percent", cfg["patterns"]["sparse_percent"]))) _, c1 = run_trial(runtime, batch, nodes, duration, det_seed, distances[layer], cfg, capture=True) _, c2 = run_trial(runtime, batch, nodes, duration, det_seed, distances[layer], cfg, capture=True) assert c1 is not None and c2 is not None deterministic = bool(np.array_equal(c1["step"], c2["step"]) and np.array_equal(c1["neuron_idx"], c2["neuron_idx"]) and np.array_equal(c1["v"], c2["v"]) and np.array_equal(c1["g"], c2["g"])) matched = all(drive_counts[f"sparse-A-{s}"] == drive_counts[f"sparse-B-{s}"] for s in map(int, cfg["patterns"]["seeds"])) separable = all(x["euclidean"] > 0 for x in pairs.values()) result[layer] = {"condition": condition, "pairwise_output_features": pairs, "trials": trials, "matched_input_count_sparse_A_B": matched, "deterministic_exact": deterministic, "all_mean_output_features_distinct": separable, "candidate_eligible": deterministic and separable and matched} write_json(exp / "separability" / f"{layer}.json", result[layer]) write_json(exp / "separability/summary.json", result) return result def select_gateway(cfg: dict[str, Any], reach: dict[str, Any], sweep: dict[str, Any], sep: dict[str, Any]) -> dict[str, Any]: eligible = [] for rank, layer in enumerate(cfg["candidate_selection"]["hierarchy_order"]): condition = sweep["viable_condition_by_layer"].get(layer) if condition: hop_hist_d = reach["layers"][layer]["min_hops_descending_distribution"] hop_hist_m = reach["layers"][layer]["min_hops_motor_distribution"] positive = [int(k) for k in list(hop_hist_d) + list(hop_hist_m) if int(k) >= 0] min_hops = min(positive) if positive else 32767 eligible.append((float(condition["amplitude_multiple_w_syn"]), min_hops, rank, layer, condition)) eligible.sort(key=lambda x: x[:3]) if not eligible: return {"selected": False, "reason": "no condition produced multi-hop, non-runaway output propagation", "rule": cfg["candidate_selection"]["rule"]} _, hops, _, layer, condition = eligible[0] return {"selected": True, "name": "visual_gateway_candidate", "layer": layer, "condition": condition, "minimum_output_hops": hops, "rule": cfg["candidate_selection"]["rule"], "handpicked": False, "adapter_v2_eligible": bool(sep.get(layer, {}).get("candidate_eligible")), "adapter_v2_gate_note": "gateway selection uses anatomy/path/propagation; formal adapter eligibility additionally requires deterministic separability"} def continuous_test(cfg: dict[str, Any], runtime: BrainRuntime, populations: dict[str, np.ndarray], distances: dict[str, np.ndarray], gateway: dict[str, Any], exp: Path) -> dict[str, Any]: if not gateway["selected"]: result = {"executed": False, "reason": gateway["reason"]}; write_json(exp / "continuous/summary.json", result); return result layer, cond = gateway["layer"], gateway["condition"]; result: dict[str, Any] = {} duration = float(cfg["continuous"]["biological_duration_ms"]) for period in map(float, cfg["continuous"]["periods_ms"]): runtime.reset(int(cfg["continuous"]["seed"])); intervals = []; previous = None sequence = cfg["continuous"]["sequence"] for i, start in enumerate(np.arange(0.0, duration, period)): kind = sequence[i % len(sequence)] batch0, nodes = pattern_batch(layer, populations, kind, float(cond["rate_hz"]), float(cond["amplitude_mv"]), period, int(cfg["continuous"]["seed"]), cfg, float(cond.get("sparsity_percent", cfg["patterns"]["sparse_percent"]))) shifted = tuple(InputEvent(target_neuron_idx=e.target_neuron_idx, stimulation_mode=e.stimulation_mode, value_kind=e.value_kind, value=e.value, start_ms=e.start_ms + float(start), duration_ms=e.duration_ms, source_id="continuous", channel_id=e.channel_id, amplitude_mv=e.amplitude_mv, metadata=e.metadata) for e in batch0.events) batch = InputBatch(events=shifted, batch_id=f"continuous:{i}") t0 = time.perf_counter(); step = runtime.step(batch, period); wall = time.perf_counter() - t0 d = runtime.read_population("descending", period); m = runtime.read_population("motor_related", period) feature = np.concatenate([d["spike_counts"], m["spike_counts"]]).astype(np.float32) recent = runtime.read_spikes(period); out = np.isin(recent["neuron_idx"], np.concatenate([d["neuron_idx"], m["neuron_idx"]])) latency = float(recent["timestamp_ms"][out].min() - start) if np.any(out) else None intervals.append({"index": i, "start_ms": start, "pattern": kind, "target_count": len(nodes), "wall_seconds": wall, "drive_events": step["new_external_drive_events"], "descending_spikes": d["total_spikes"], "motor_spikes": m["total_spikes"], "output_latency_ms": latency, "output_feature_change_l2": None if previous is None else float(np.linalg.norm(feature - previous))}) previous = feature zero = [x for x in intervals if x["pattern"] == "zero"] result[str(int(period))] = {"period_ms": period, "interval_count": len(intervals), "continuous_state_without_reset": True, "wall_seconds": float(sum(x["wall_seconds"] for x in intervals)), "zero_period_output_spikes": int(sum(x["descending_spikes"] + x["motor_spikes"] for x in zero)), "intervals": intervals} write_json(exp / "continuous" / f"period-{int(period)}ms.json", result[str(int(period))]) write_json(exp / "continuous/summary.json", result) return result def regression_suite(cfg: dict[str, Any], runtime: BrainRuntime, exp: Path) -> dict[str, Any]: io_root = Path(cfg["paths"]["io_root"]); sys.path.insert(0, str(io_root)) import phase4d as p4d reg = p4d.regression(exp, cfg, runtime) api = p4d.api_tests(exp, runtime) stream = p4d.stream_contract(exp, runtime) nodes = np.flatnonzero(runtime.population_masks["auditory"]).astype(np.int64) batch = p4d.direct_batch(nodes, "auditory", 10, 0.0, 20.0, "phase4e-continuous-contract") runtime.reset(45999); runtime.step(batch, 10.0); checkpoint = runtime.checkpoint_state(); runtime.step(batch, 10.0) first = runtime.read_spikes(); first_v = runtime.v.cpu().numpy().copy(); runtime.restore_state(checkpoint); runtime.step(batch, 10.0) replay = runtime.read_spikes(); replay_v = runtime.v.cpu().numpy().copy() continuous = {"checkpoint_replay_events_exact": bool(np.array_equal(first["step"], replay["step"]) and np.array_equal(first["neuron_idx"], replay["neuron_idx"])), "checkpoint_replay_voltage_exact": bool(np.array_equal(first_v, replay_v)), "clock_ms": runtime.simulation_time_ms} continuous["passed"] = continuous["checkpoint_replay_events_exact"] and continuous["checkpoint_replay_voltage_exact"] write_json(exp / "audit/phase4d-continuous-contract.json", continuous) result = {"phase4c_five_trials_passed": all(x["passed"] for x in reg.values()), "phase4d_api_passed": api["all_passed"], "phase4d_continuous_checkpoint_passed": continuous["passed"], "activity_replay_passed": stream["replay_exact"], "trials": reg, "api": api, "continuous": continuous, "activity_stream": stream} result["all_passed"] = all(result[k] for k in ("phase4c_five_trials_passed", "phase4d_api_passed", "phase4d_continuous_checkpoint_passed", "activity_replay_passed")) write_json(exp / "audit/regression-summary.json", result) return result def freeze_candidate(cfg: dict[str, Any], gateway: dict[str, Any], sep: dict[str, Any], continuous: dict[str, Any], regression: dict[str, Any], code: Path, config: Path, exp: Path) -> dict[str, Any]: root = Path(cfg["paths"]["candidate_root"]); root.mkdir(parents=True, exist_ok=True) for src in (code, config): dst = root / src.name if src.resolve() != dst.resolve(): shutil.copy2(src, dst) passed = bool(gateway["selected"] and gateway.get("adapter_v2_eligible") and regression["all_passed"] and continuous.get("executed", True)) version = {"name": cfg["names"]["candidate"], "status": "CANDIDATE-PASS" if passed else "CANDIDATE-REJECTED", "formal_v2_frozen": False, "gateway": gateway, "connectome_changed": False, "lif_semantics_changed": False, "transmitter_sign_changed": False, "phase4f_started": False} write_json(root / "VERSION.json", version) files = sorted(p for p in root.iterdir() if p.is_file() and p.name not in {"MANIFEST.json", "checksums.sha256"}) write_json(root / "MANIFEST.json", {"artifact": cfg["names"]["candidate"], "files": [{"path": p.name, "bytes": p.stat().st_size, "sha256": sha256(p)} for p in files]}) targets = files + [root / "MANIFEST.json"] (root / "checksums.sha256").write_text("".join(f"{sha256(p)} {p.name}\n" for p in targets), encoding="utf-8") return version def materialize_requested_layout(exp: Path) -> None: """Expose the requested top-level sweep/gateway layout as real files.""" aliases = { "calibration/amplitude": "amplitude-sweep", "calibration/rate": "rate-sweep", "calibration/sparsity": "sparsity-sweep", "calibration/latency": "latency", } for source, destination in aliases.items(): dst = exp / destination; dst.mkdir(parents=True, exist_ok=True) for path in (exp / source).glob("*.json"): shutil.copy2(path, dst / path.name) gateway_dir = exp / "gateway"; gateway_dir.mkdir(parents=True, exist_ok=True) shutil.copy2(exp / "audit/visual-gateway-candidate.json", gateway_dir / "visual-gateway-candidate.json") def seal(cfg: dict[str, Any], exp: Path, code: Path, config: Path, controls: list[Path], summary: dict[str, Any]) -> None: excluded = {exp / "manifests/manifest.json", exp / "manifests/checksums.sha256"} outputs = sorted(p for p in exp.rglob("*") if p.is_file() and p not in excluded) manifest = {"phase": "4E", "status": summary["status"], "created_unix_time": time.time(), "software": {"python": platform.python_version(), "torch": torch.__version__, "numpy": np.__version__, "scipy": __import__("scipy").__version__, "pyarrow": pa.__version__}, "code": [{"path": str(x), "sha256": sha256(x)} for x in (code, config)], "immutable_inputs": [{"path": str(x), "sha256": sha256(x)} for x in controls], "outputs": [{"path": str(x), "bytes": x.stat().st_size, "sha256": sha256(x)} for x in outputs], "prohibitions": cfg["prohibitions"]} write_json(exp / "manifests/manifest.json", manifest) targets = outputs + [exp / "manifests/manifest.json"] (exp / "manifests/checksums.sha256").write_text("".join(f"{sha256(x)} {x.relative_to(exp)}\n" for x in targets), encoding="utf-8") def main() -> None: ap = argparse.ArgumentParser(); ap.add_argument("--config", type=Path, required=True); args = ap.parse_args() cfg = yaml.safe_load(args.config.read_text(encoding="utf-8")); exp = Path(cfg["paths"]["experiment_root"]) for sub in ("configs", "populations", "reachability", "calibration", "amplitude-sweep", "rate-sweep", "sparsity-sweep", "latency", "gateway", "separability", "continuous", "regression", "activity-stream", "audit", "manifests", "compiler-reports/extension-build"): (exp / sub).mkdir(parents=True, exist_ok=True) shutil.copy2(args.config, exp / "configs/phase4e.yaml") controls, before, immutable = immutable_audit(cfg); write_json(exp / "audit/immutable-inputs-before.json", immutable) # Frozen runtime-v1 imports the shared, unchanged Phase 4B graph helpers # from the model root. sys.path.insert(0, cfg["paths"]["model_root"]) runtime = BrainRuntime(args.config, exp / "compiler-reports/extension-build") populations, population_summary = load_populations(cfg, runtime, exp) semantics = {"phase4d_external_drive_application": "g[target] += amplitude_mv per selected external Poisson event, after g decay/connectome arrival and before membrane integration", "units": "mV-equivalent conductance state increment", "phase4d_amplitude_mv": 8.0, "w_syn_mv_per_anatomical_synapse": cfg["model"]["w_syn_mv_per_anatomical_synapse"], "ratio_phase4d_drive_to_one_anatomical_synapse": 8.0 / float(cfg["model"]["w_syn_mv_per_anatomical_synapse"]), "connectome_transmission": False, "raw_csr_weight_modified": False} write_json(exp / "audit/external-drive-semantics.json", semantics) reach = graph_reachability(cfg, runtime, populations, exp) distances = forward_distances(runtime, populations, exp) sweep, records = calibration_sweeps(cfg, runtime, populations, distances, exp) viable = sweep["viable_condition_by_layer"] sep = separability_and_determinism(cfg, runtime, populations, distances, viable, exp) gateway = select_gateway(cfg, reach, sweep, sep); write_json(exp / "audit/visual-gateway-candidate.json", gateway) continuous = continuous_test(cfg, runtime, populations, distances, gateway, exp) # Run frozen regressions in a fresh runtime, isolated from the long visual # stress sweeps and their allocator history. del runtime gc.collect(); torch.cuda.empty_cache() regression_runtime = BrainRuntime(args.config, exp / "compiler-reports/extension-build") regression = regression_suite(cfg, regression_runtime, exp) del regression_runtime gc.collect(); torch.cuda.empty_cache() materialize_requested_layout(exp) after = {str(x): sha256(x) for x in controls}; unchanged = before == after candidate = freeze_candidate(cfg, gateway, sep, continuous, regression, Path(__file__).resolve(), args.config.resolve(), exp) stable = all(not x["aborted_runaway"] and x["numerically_finite"] for x in records if x.get("multi_hop_output_propagation")) invariant_pass = regression["all_passed"] and unchanged and stable status = ("PASS" if invariant_pass and gateway.get("adapter_v2_eligible") else "FAIL-NO-V2-CANDIDATE" if invariant_pass else "FAILED") summary = {"status": status, "reference_phase4d_unchanged": unchanged, "immutable_checks": immutable, "visual_population_summary": population_summary, "external_drive_semantics": semantics, "reachability": reach, "calibration": sweep, "gateway": gateway, "separability": sep, "continuous": continuous, "regression": regression, "candidate_artifact": candidate, "propagating_trials_non_runaway_and_finite": stable, "phase4f_started": False, "game_rl_web_3d_started": False} write_json(exp / "audit/summary.json", summary) seal(cfg, exp, Path(__file__).resolve(), args.config.resolve(), controls, summary) if status == "FAILED": raise RuntimeError("Phase 4E invariant/regression failure") if __name__ == "__main__": main()