Digital Drosophila Archive
669 lines
42 KiB
Python
669 lines
42 KiB
Python
#!/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()
|