Digital Drosophila Archive
503 lines
28 KiB
Python
503 lines
28 KiB
Python
"""Small, fixed-input Brian2 2.5.1 stress oracle for Phase 4AE.
|
|
|
|
The archived author source, equations, reset, parameters, and synapses are left
|
|
unchanged. Diagnostic initial states and the old common request interface are
|
|
explicit inputs. This file does not install packages or run whole-brain cases.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import hashlib
|
|
import importlib.metadata
|
|
import json
|
|
import platform
|
|
import sys
|
|
import types
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
|
|
AUTHOR_PATH = Path("/data/digital-drosophila/experiments/phase4ad-author-semantics/external-resources/model.py")
|
|
AUTHOR_SHA256 = "fc45837d7122c6ce2a7f3f2f23c515992e4b232aadb919efabb72337fac88e4e"
|
|
AUTHOR_COMMIT = "91bdd1e7dcf193f3e7ca5a8933497fcef63b7960"
|
|
LIO_NAMES = {"gd": "_lio_2", "c0": "_lio_8", "cg": "_lio_9", "cv": "_lio_10"}
|
|
|
|
|
|
def _save(path: Path, value) -> None:
|
|
with path.open("x", encoding="utf-8") as handle:
|
|
json.dump(value, handle, indent=2, sort_keys=True, allow_nan=False)
|
|
handle.write("\n")
|
|
|
|
|
|
def _save_arrays(path: Path, arrays: dict) -> None:
|
|
with path.open("xb") as handle:
|
|
np.savez_compressed(handle, **arrays)
|
|
|
|
|
|
def _digest(path: Path) -> str:
|
|
h = hashlib.sha256()
|
|
with path.open("rb") as handle:
|
|
for block in iter(lambda: handle.read(8 << 20), b""):
|
|
h.update(block)
|
|
return h.hexdigest()
|
|
|
|
|
|
def _check(condition, message: str) -> None:
|
|
if not condition:
|
|
raise ValueError(message)
|
|
|
|
|
|
def _load_case(case_dir: Path) -> tuple[dict, dict]:
|
|
spec = json.loads((case_dir / "spec.json").read_text(encoding="utf-8"))
|
|
with np.load(case_dir / "arrays.npz", allow_pickle=False) as archive:
|
|
arrays = {name: archive[name].copy() for name in archive.files}
|
|
expected_types = {"nodes": "int32", "edge_id": "int64", "src": "int32",
|
|
"dst": "int32", "signed_raw": "int64", "requests": "uint8",
|
|
"external_mv": "float64", "initial_v_mv": "float64",
|
|
"initial_g_mv": "float64"}
|
|
for name, dtype in expected_types.items():
|
|
_check(name in arrays, f"Missing input array: {name}")
|
|
_check(arrays[name].dtype == np.dtype(dtype), f"{name} must have dtype {dtype}")
|
|
nodes = arrays["nodes"]
|
|
_check(nodes.ndim == 1 and 1 <= len(nodes) <= 512, "Stress cases require 1..512 nodes")
|
|
_check(np.all(np.diff(nodes.astype(np.int64)) > 0), "nodes must be sorted and unique")
|
|
steps = int(spec["steps"])
|
|
_check(steps == spec["steps"] and 1 <= steps <= 1000, "steps must be an integer in 1..1000")
|
|
_check(float(spec["dt_ms"]) == .1, "The fixed stress protocol requires dt_ms=0.1")
|
|
if "duration_ms" in spec:
|
|
_check(np.isclose(float(spec["duration_ms"]), steps * .1, rtol=0, atol=1e-12),
|
|
"duration_ms disagrees with steps and dt_ms")
|
|
if "dt_ms" in spec["model"]:
|
|
_check(float(spec["model"]["dt_ms"]) == .1, "model dt_ms disagrees with case dt_ms")
|
|
edge_count = len(arrays["edge_id"])
|
|
# The unchanged author calls Synapses.connect(i=[], j=[]), which Brian2
|
|
# 2.5.1 does not support. Empty-network fixtures must explicitly provide
|
|
# a documented zero-weight synthetic edge in the input case.
|
|
_check(edge_count > 0, "Use an explicit zero-weight fixture edge instead of empty connectivity")
|
|
for name in ("edge_id", "src", "dst", "signed_raw"):
|
|
_check(arrays[name].shape == (edge_count,), f"{name} has an inconsistent shape")
|
|
_check(len(np.unique(arrays["edge_id"])) == edge_count, "edge_id must be unique")
|
|
for name in ("src", "dst"):
|
|
_check(np.all((arrays[name] >= 0) & (arrays[name] < len(nodes))), f"{name} is out of local range")
|
|
_check(np.array_equal(np.lexsort((arrays["edge_id"], arrays["src"])), np.arange(edge_count)),
|
|
"Edges must be in canonical source-then-edge order")
|
|
for name in ("requests", "external_mv"):
|
|
_check(arrays[name].shape == (steps, len(nodes)), f"{name} must have shape (steps, N)")
|
|
_check(np.all(arrays["requests"] <= 1), "requests must contain only 0 or 1")
|
|
for name in ("initial_v_mv", "initial_g_mv"):
|
|
_check(arrays[name].shape == (len(nodes),), f"{name} must have shape (N,)")
|
|
for name in ("external_mv", "initial_v_mv", "initial_g_mv"):
|
|
_check(np.isfinite(arrays[name]).all(), f"{name} contains a nonfinite value")
|
|
return spec, arrays
|
|
|
|
|
|
def _code_records(network) -> list:
|
|
records = []
|
|
for obj in network.sorted_objects:
|
|
row = {"name": obj.name, "type": type(obj).__name__,
|
|
"when": obj.when, "order": int(obj.order), "active": bool(obj.active)}
|
|
if hasattr(obj, "abstract_code"):
|
|
abstract = obj.abstract_code
|
|
row["abstract_code"] = abstract if isinstance(abstract, (str, dict)) else str(abstract)
|
|
row["code_objects"] = [
|
|
{"name": code.name, "generated_code": str(code.code)}
|
|
for code in getattr(obj, "code_objects", ())]
|
|
records.append(row)
|
|
return records
|
|
|
|
|
|
class _QueueReadProxy:
|
|
"""Forward the queue unchanged; copy only observations from peek()."""
|
|
|
|
def __init__(self, queue, get_tick, get_available, targets, observations, counters):
|
|
self._original_queue = queue
|
|
self._get_tick = get_tick
|
|
self._get_available = get_available
|
|
self._targets = targets
|
|
self._observations = observations
|
|
self._counters = counters
|
|
|
|
def peek(self, *args, **kwargs):
|
|
value = self._original_queue.peek(*args, **kwargs)
|
|
edge_indices = np.asarray(value, dtype=np.int64).copy()
|
|
if np.any(edge_indices < 0) or np.any(edge_indices >= len(self._targets)):
|
|
raise AssertionError("The anatomical queue returned an invalid local edge index")
|
|
tick = self._get_tick()
|
|
accepted = np.asarray(self._get_available(), dtype=bool)[self._targets[edge_indices]].astype(np.uint8)
|
|
self._observations.append((tick, edge_indices, accepted))
|
|
self._counters["peek_calls"] += 1
|
|
return value
|
|
|
|
def advance(self, *args, **kwargs):
|
|
self._counters["advance_calls"] += 1
|
|
return self._original_queue.advance(*args, **kwargs)
|
|
|
|
def __getattr__(self, name):
|
|
return getattr(self._original_queue, name)
|
|
|
|
|
|
def _scalar_namespace(code_runner, name: str):
|
|
for code in getattr(code_runner, "code_objects", ()):
|
|
if name in code.namespace:
|
|
value = np.asarray(code.namespace[name])
|
|
if value.ndim == 0 and np.isfinite(value):
|
|
return float(value), {"source": "CodeObject.namespace", "code_object": code.name,
|
|
"name": name, "namespace_dtype": str(value.dtype)}
|
|
return None, None
|
|
|
|
|
|
def _coefficients(neu, parameters, b2, dt_seconds: float) -> dict:
|
|
generated = "\n".join(str(code.code) for code in neu.state_updater.code_objects)
|
|
_check("_lio_2 * g" in generated and
|
|
"_lio_8 + ((_lio_9 * g) + (_lio_10 * v))" in generated,
|
|
"Generated state expression changed; the recorded _lio-to-GPU mapping is not valid")
|
|
dt = np.float64(dt_seconds)
|
|
tau = np.float64(parameters["tau"] / b2.second)
|
|
membrane = np.float64(parameters["t_mbr"] / b2.second)
|
|
rest = np.float64(parameters["v_0"] / b2.volt)
|
|
# This is the actual scalar expression order emitted by the pinned Brian
|
|
# state updater, not an algebraically rearranged analytic solution.
|
|
lio_1 = (-dt) / tau
|
|
lio_3 = tau / (membrane - tau)
|
|
lio_4 = dt / membrane
|
|
lio_5 = dt / tau
|
|
lio_6 = (-dt) / membrane
|
|
fallback = {
|
|
"gd": np.exp(lio_1), "c0": rest - (rest * np.exp(lio_6)),
|
|
"cg": ((lio_3 * ((-np.exp(lio_4)) + np.exp(lio_5))) * np.exp(lio_6)) * np.exp(lio_1),
|
|
"cv": np.exp(lio_6),
|
|
}
|
|
result, provenance = {}, {}
|
|
for key, name in LIO_NAMES.items():
|
|
value, record = _scalar_namespace(neu.state_updater, name)
|
|
if value is None:
|
|
value = float(fallback[key])
|
|
record = {"source": "same generated scalar expression evaluated with numpy.float64", "name": name}
|
|
result[key], provenance[key] = value, record
|
|
for key, runner, name, author_key in (
|
|
("threshold", neu.thresholder["spike"], "v_th", "v_th"),
|
|
("reset", neu.resetter["spike"], "v_rst", "v_rst"),
|
|
):
|
|
value, record = _scalar_namespace(runner, name)
|
|
if value is None:
|
|
value = float(parameters[author_key] / b2.volt)
|
|
record = {"source": "unchanged author parameter", "name": name}
|
|
result[key], provenance[key] = value, record
|
|
result["seed_volt"] = float((parameters["v_th"] + 1 * b2.mV) / b2.volt)
|
|
provenance["seed_volt"] = {"source": "common request interface: unchanged v_th + 1 mV"}
|
|
result["metadata"] = {
|
|
"units": "SI volts and dimensionless coefficients", "names": LIO_NAMES,
|
|
"provenance": provenance, "hex": {key: value.hex() for key, value in result.items()},
|
|
"voltage_expression": "c0 + ((cg * old_g) + (cv * old_v))",
|
|
"g_expression": "gd * old_g",
|
|
"constant_values_are_double": True,
|
|
"float32_execution_note": "Brian NumPy scalar/array promotion and operation order must be matched separately; double coefficients alone do not assert float64 intermediate array arithmetic.",
|
|
}
|
|
return result
|
|
|
|
|
|
def _parameters_match(parameters, model, b2) -> dict:
|
|
rows = {}
|
|
for config_key, author_key, unit in (
|
|
("membrane_tau_ms", "t_mbr", b2.ms), ("synaptic_tau_ms", "tau", b2.ms),
|
|
("refractory_ms", "t_rfc", b2.ms), ("synaptic_delay_ms", "t_dly", b2.ms),
|
|
("resting_potential_mv", "v_0", b2.mV), ("reset_potential_mv", "v_rst", b2.mV),
|
|
("threshold_mv", "v_th", b2.mV), ("w_syn_mv_per_anatomical_synapse", "w_syn", b2.mV),
|
|
):
|
|
requested, actual = float(model[config_key]), float(parameters[author_key] / unit)
|
|
_check(np.isfinite(requested) and np.isclose(requested, actual, rtol=0, atol=1e-12),
|
|
f"Model parameter {config_key} differs from the archived author")
|
|
rows[config_key] = {"case": requested, "author": actual}
|
|
return rows
|
|
|
|
|
|
def _run_once(spec, arrays, precision, output, source_bytes, trace_queue):
|
|
import brian2 as b2
|
|
from brian2.core.base import BrianObject
|
|
|
|
output.mkdir(parents=True, exist_ok=False)
|
|
saved = {key: b2.prefs[key] for key in (
|
|
"codegen.target", "core.default_float_dtype", "legacy.refractory_timing")}
|
|
old_dt = b2.defaultclock.dt
|
|
dtype = np.dtype(precision)
|
|
steps, nodes = int(spec["steps"]), arrays["nodes"]
|
|
dt_ms, dt_seconds = float(spec["dt_ms"]), float(spec["dt_ms"]) * .001
|
|
try:
|
|
b2.start_scope()
|
|
b2.prefs["codegen.target"] = "numpy"
|
|
b2.prefs["core.default_float_dtype"] = dtype.type
|
|
b2.prefs["legacy.refractory_timing"] = False
|
|
b2.defaultclock.dt = dt_ms * b2.ms
|
|
author = types.ModuleType("phase4ae_author_" + precision)
|
|
author.__file__ = str(AUTHOR_PATH)
|
|
exec(compile(source_bytes, str(AUTHOR_PATH), "exec"), author.__dict__)
|
|
parameters = dict(author.default_params)
|
|
parameter_checks = _parameters_match(parameters, spec["model"], b2)
|
|
completeness = pd.DataFrame({"local_index": np.arange(len(nodes), dtype=np.int32)},
|
|
index=pd.Index(nodes, name="global_neuron_idx"))
|
|
connectivity = pd.DataFrame({
|
|
"Presynaptic_Index": arrays["src"], "Postsynaptic_Index": arrays["dst"],
|
|
"Excitatory x Connectivity": arrays["signed_raw"], "edge_id": arrays["edge_id"],
|
|
"source_neuron_idx": nodes[arrays["src"]], "target_neuron_idx": nodes[arrays["dst"]],
|
|
})
|
|
comp_path, con_path = output / "completeness.csv", output / "connectivity.parquet"
|
|
completeness.to_csv(comp_path)
|
|
connectivity.to_parquet(con_path, engine="pyarrow", index=False)
|
|
neu, syn, spikes = author.create_model(comp_path, con_path, parameters)
|
|
_check(np.dtype(neu.variables["v"].dtype) == dtype and np.dtype(neu.variables["g"].dtype) == dtype,
|
|
"Brian state dtype disagrees with the precision arm")
|
|
_check(np.dtype(syn.variables["w"].dtype) == dtype, "Brian synaptic weight dtype disagrees")
|
|
np.testing.assert_array_equal(np.asarray(syn.i[:]), arrays["src"])
|
|
np.testing.assert_array_equal(np.asarray(syn.j[:]), arrays["dst"])
|
|
neu.v = arrays["initial_v_mv"] * b2.mV
|
|
neu.g = arrays["initial_g_mv"] * b2.mV
|
|
initial_v = np.asarray(neu.v_[:], dtype=dtype).copy()
|
|
initial_g = np.asarray(neu.g_[:], dtype=dtype).copy()
|
|
attempted = arrays["requests"].copy()
|
|
accepted = np.zeros_like(attempted)
|
|
rejected = np.zeros_like(attempted)
|
|
request_slots_seen = np.zeros(steps, dtype=np.uint8)
|
|
|
|
def current_tick():
|
|
return int(round(float(neu.clock.t / b2.second) / dt_seconds))
|
|
|
|
def requests_before_thresholds():
|
|
tick = current_tick()
|
|
if tick < 0 or tick >= steps or request_slots_seen[tick]:
|
|
raise AssertionError("Common request operation executed outside its unique tick")
|
|
request_slots_seen[tick] = 1
|
|
request = attempted[tick].astype(bool)
|
|
available = np.asarray(neu.not_refractory[:], dtype=bool)
|
|
permitted = request & available
|
|
accepted[tick] = permitted
|
|
rejected[tick] = request & ~available
|
|
if permitted.any():
|
|
neu.v[np.flatnonzero(permitted)] = parameters["v_th"] + 1 * b2.mV
|
|
|
|
request_op = b2.NetworkOperation(requests_before_thresholds, clock=neu.clock,
|
|
when="before_thresholds", order=0, name="oracle_common_requests")
|
|
objects = [neu, syn, spikes, request_op]
|
|
ext_ticks, ext_targets = np.nonzero(arrays["external_mv"])
|
|
external_weights = np.empty(0, dtype=dtype)
|
|
if len(ext_ticks):
|
|
ext_ids = np.arange(len(ext_ticks), dtype=np.int32)
|
|
ext_source = b2.SpikeGeneratorGroup(len(ext_ticks), ext_ids,
|
|
ext_ticks.astype(np.float64) * dt_ms * b2.ms, sorted=True,
|
|
clock=neu.clock, name="zz_external_sources")
|
|
ext_syn = b2.Synapses(ext_source, neu, "w : volt", on_pre="g_post += w",
|
|
delay=0 * b2.ms, name="zz_external_arrivals")
|
|
ext_syn.connect(i=ext_ids, j=ext_targets)
|
|
ext_syn.w = arrays["external_mv"][ext_ticks, ext_targets] * b2.mV
|
|
ext_syn.pre.order = 0
|
|
_check(syn.pre.when == ext_syn.pre.when == "synapses" and syn.pre.order < 0,
|
|
"External delivery must follow the original anatomical pathway")
|
|
external_weights = np.asarray(ext_syn.w_[:], dtype=dtype).copy()
|
|
objects.extend([ext_source, ext_syn])
|
|
|
|
observations = []
|
|
queue_counters = {"peek_calls": 0, "advance_calls": 0, "installations": []}
|
|
if trace_queue:
|
|
class QueueInstaller(BrianObject):
|
|
def __init__(self):
|
|
super().__init__(clock=neu.clock, when="end", order=1_000_000,
|
|
name="zz_oracle_queue_installer")
|
|
|
|
def before_run(self, run_namespace):
|
|
# Network.before_run traverses sorted objects. All original
|
|
# code runners have compiled before this last end-slot object.
|
|
installed = []
|
|
for code in syn.pre.code_objects:
|
|
if "_queue" not in code.namespace or "_queue.peek()" not in str(code.code):
|
|
continue
|
|
original = code.namespace["_queue"]
|
|
if isinstance(original, _QueueReadProxy):
|
|
original = original._original_queue
|
|
code.namespace["_queue"] = _QueueReadProxy(
|
|
original, current_tick, lambda: neu.not_refractory[:],
|
|
arrays["dst"], observations, queue_counters)
|
|
installed.append(code.name)
|
|
if len(installed) != 1:
|
|
raise AssertionError(f"Expected exactly one anatomical queue consumer, got {installed}")
|
|
queue_counters["installations"].append(installed[0])
|
|
|
|
objects.append(QueueInstaller())
|
|
|
|
trace = b2.StateMonitor(neu, ["v", "g", "lastspike", "not_refractory"],
|
|
record=True, when="end", order=0, name="oracle_end_state")
|
|
objects.append(trace)
|
|
network = b2.Network(*objects)
|
|
_save(output / "schedule.json", {"schedule": list(network.schedule),
|
|
"summary": str(network.scheduling_summary())})
|
|
_save(output / "input-audit.json", {
|
|
"source_sha256": AUTHOR_SHA256, "parameter_checks": parameter_checks,
|
|
"source_parameters_modified": False, "precision": precision,
|
|
"dt_ms": dt_ms, "steps": steps, "request_interface": "before_thresholds; v=v_th+1mV only if not_refractory",
|
|
"poisson_helper_used": False, "all_neurons_keep_original_rfc": True,
|
|
"queue_trace_enabled": trace_queue, "state_units": "SI volts, original dtype",
|
|
"completeness_sha256": _digest(comp_path), "connectivity_sha256": _digest(con_path),
|
|
})
|
|
network.run(0 * b2.ms, namespace={})
|
|
_save(output / "compiled-code.json", _code_records(network))
|
|
if request_slots_seen.any() or observations:
|
|
raise AssertionError("Zero-duration compilation unexpectedly advanced the case")
|
|
network.run(steps * dt_ms * b2.ms, namespace={})
|
|
_save(output / "generated-code.json", _code_records(network))
|
|
_check(np.all(request_slots_seen == 1), "Not every request tick was observed exactly once")
|
|
|
|
event_float = np.asarray(spikes.t[:] / b2.second, dtype=np.float64) / dt_seconds
|
|
event_ticks = np.rint(event_float).astype(np.int32)
|
|
_check(np.allclose(event_float, event_ticks, rtol=0, atol=1e-6), "Spike time is off the integer clock")
|
|
local_indices = np.asarray(spikes.i[:], dtype=np.int32).copy()
|
|
order = np.lexsort((nodes[local_indices], event_ticks))
|
|
event_ticks, local_indices = event_ticks[order], local_indices[order]
|
|
fired = np.zeros_like(attempted)
|
|
fired[event_ticks, local_indices] = 1
|
|
counts = np.asarray(spikes.count[:], dtype=np.int32).copy()
|
|
_check(np.array_equal(fired.sum(axis=0), counts), "Events and per-neuron counts disagree")
|
|
_check(not np.any(accepted & ~fired), "An applicable common request did not produce a real spike")
|
|
_check(np.array_equal(attempted, accepted + rejected), "Request accounting is incomplete")
|
|
state_ticks = np.asarray(trace.t[:] / b2.second, dtype=np.float64) / dt_seconds
|
|
_check(state_ticks.shape == (steps,) and np.allclose(state_ticks, np.arange(steps), rtol=0, atol=1e-6),
|
|
"StateMonitor does not cover the exact requested ticks")
|
|
voltage = np.asarray(trace.v_, dtype=dtype).T.copy()
|
|
conductance = np.asarray(trace.g_, dtype=dtype).T.copy()
|
|
last = np.asarray(trace.lastspike_, dtype=np.float64).T
|
|
seen = np.isfinite(last) & (last >= 0)
|
|
last_ticks = np.zeros(last.shape, dtype=np.int64)
|
|
last_ticks[seen] = np.rint(last[seen] / dt_seconds).astype(np.int64)
|
|
refractory_steps = int(round(float(parameters["t_rfc"] / b2.ms) / dt_ms))
|
|
remaining = np.where(seen, np.maximum(refractory_steps - (np.arange(steps)[:, None] - last_ticks), 0), 0).astype(np.int32)
|
|
available_end = np.asarray(trace.not_refractory, dtype=bool).T
|
|
_check(np.array_equal(remaining == 0, available_end), "Refractory flag and integer countdown disagree")
|
|
history = np.stack([voltage, conductance, remaining.astype(dtype)], axis=1)
|
|
_check(np.isfinite(history).all(), "Author state became nonfinite")
|
|
arrival_step = np.empty(0, dtype=np.int32)
|
|
arrival_edge = np.empty(0, dtype=np.int64)
|
|
arrival_accepted = np.empty(0, dtype=np.uint8)
|
|
if trace_queue:
|
|
_check(queue_counters["peek_calls"] == steps and queue_counters["advance_calls"] == steps,
|
|
"Queue proxy did not observe exactly one peek and advance per tick")
|
|
_check([entry[0] for entry in observations] == list(range(steps)), "Queue peek ticks are not complete and ordered")
|
|
arrival_step = np.concatenate([np.full(len(ids), tick, dtype=np.int32) for tick, ids, _ in observations])
|
|
arrival_edge = np.concatenate([ids for _, ids, _ in observations]).astype(np.int64, copy=False)
|
|
arrival_accepted = np.concatenate([mask for _, _, mask in observations]).astype(np.uint8, copy=False)
|
|
delay_steps = int(round(float(parameters["t_dly"] / b2.ms) / dt_ms))
|
|
if len(arrival_step):
|
|
_check(np.all(arrival_step >= delay_steps), "An anatomical event arrived before the fixed delay")
|
|
_check(np.all(fired[arrival_step - delay_steps, arrays["src"][arrival_edge]]),
|
|
"A queued arrival has no matching presynaptic spike at the author delay")
|
|
_check(np.array_equal(arrival_accepted,
|
|
available_end[arrival_step, arrays["dst"][arrival_edge]].astype(np.uint8)),
|
|
"Observed arrival gating disagrees with the post-threshold refractory state")
|
|
source_counts_due = fired[:max(0, steps - delay_steps)].sum(axis=0, dtype=np.int64)
|
|
expected_arrivals = int(source_counts_due[arrays["src"]].sum(dtype=np.int64))
|
|
_check(len(arrival_edge) == expected_arrivals, "Queue trace omitted or duplicated an anatomical arrival")
|
|
result = {
|
|
"nodes": nodes.copy(), "step": event_ticks, "local_idx": local_indices,
|
|
"neuron_idx": nodes[local_indices], "counts": counts,
|
|
"v": voltage[-1].copy(), "g": conductance[-1].copy(), "ref": remaining[-1].copy(),
|
|
"history": history, "ref_history": remaining,
|
|
"not_refractory": available_end.astype(np.uint8),
|
|
"weights": np.asarray(syn.w_[:], dtype=dtype).copy(),
|
|
"initial_v": initial_v, "initial_g": initial_g,
|
|
"request_attempted": attempted, "request_accepted": accepted, "request_rejected": rejected,
|
|
"request_count": attempted.sum(axis=1, dtype=np.int32),
|
|
"applicable_count": accepted.sum(axis=1, dtype=np.int32),
|
|
"rejected_count": rejected.sum(axis=1, dtype=np.int32), "fired": fired,
|
|
"arrival_step": arrival_step, "arrival_edge_index": arrival_edge,
|
|
"arrival_accepted": arrival_accepted,
|
|
"external_step": ext_ticks.astype(np.int32), "external_target": ext_targets.astype(np.int32),
|
|
"external_weights": external_weights,
|
|
"external_accepted": available_end[ext_ticks, ext_targets].astype(np.uint8),
|
|
}
|
|
coefficients = _coefficients(neu, parameters, b2, dt_seconds)
|
|
_save(output / "coefficients.json", coefficients)
|
|
_save(output / "queue-audit.json", {
|
|
**queue_counters, "trace_enabled": trace_queue, "due_events": len(arrival_edge),
|
|
"accepted_due_events": int(arrival_accepted.sum()),
|
|
"all_due_events_preserved": True if trace_queue else None,
|
|
"ordering": "raw anatomical queue.peek local edge-index sequence; no sorting or filtering",
|
|
"proxy_contract": "peek returns original object; advance and all other methods forward unchanged",
|
|
})
|
|
_save_arrays(output / "state.npz", result)
|
|
return result, coefficients
|
|
except Exception as error:
|
|
_save(output / "error.json", {"type": type(error).__name__, "message": str(error)})
|
|
raise
|
|
finally:
|
|
for key, value in saved.items():
|
|
b2.prefs[key] = value
|
|
b2.defaultclock.dt = old_dt
|
|
|
|
|
|
def run_case(case_dir: Path, out_dir: Path, precision: str) -> dict:
|
|
"""Run one <=512-node, <=1000-step fixed case with exclusive outputs.
|
|
|
|
Main outputs are ``state.npz`` and ``coefficients.json``. Queue tracing is
|
|
enabled by default. ``spec.trace_queue`` can disable it; when enabled,
|
|
``spec.validate_queue_proxy`` defaults to True and repeats the same case
|
|
without instrumentation. That repeat must match every non-arrival array
|
|
exactly. An opt-out is recorded as unverified, never as a successful check.
|
|
This function changes global Brian preferences temporarily; run it serially.
|
|
"""
|
|
_check(precision in ("float64", "float32"), "precision must be float64 or float32")
|
|
case_dir, out_dir = Path(case_dir).resolve(strict=True), Path(out_dir).resolve()
|
|
spec, arrays = _load_case(case_dir)
|
|
source_bytes = AUTHOR_PATH.read_bytes()
|
|
source_hash = hashlib.sha256(source_bytes).hexdigest()
|
|
_check(source_hash == AUTHOR_SHA256, "Archived author source hash mismatch")
|
|
trace_queue = spec.get("trace_queue", True)
|
|
validate_proxy = spec.get("validate_queue_proxy", trace_queue)
|
|
_check(isinstance(trace_queue, bool) and isinstance(validate_proxy, bool), "Queue switches must be JSON booleans")
|
|
_check(not validate_proxy or trace_queue, "Proxy validation requires queue tracing")
|
|
out_dir.mkdir(parents=True, exist_ok=False)
|
|
try:
|
|
versions = {}
|
|
for name in ("brian2", "numpy", "pandas", "scipy", "pyarrow", "sympy", "cython"):
|
|
versions[name] = importlib.metadata.version(name)
|
|
_save(out_dir / "environment.json", {"python": sys.version, "platform": platform.platform(),
|
|
"packages": versions, "backend": "Brian2 numpy runtime"})
|
|
_check(versions["brian2"] == "2.5.1", "The stress oracle requires Brian2 2.5.1")
|
|
_check(versions["numpy"] == "1.22.3", "The stress oracle requires the recorded NumPy 1.22.3 promotion rules")
|
|
_save(out_dir / "case.json", {"spec": spec, "precision": precision,
|
|
"author_path": str(AUTHOR_PATH), "author_commit": AUTHOR_COMMIT, "author_sha256": source_hash,
|
|
"case_spec_sha256": _digest(case_dir / "spec.json"), "case_arrays_sha256": _digest(case_dir / "arrays.npz"),
|
|
"trace_queue": trace_queue, "validate_queue_proxy": validate_proxy,
|
|
"common_interface_preserved": True, "original_poi_interface_adopted": False,
|
|
"whole_brain": False, "learning": False})
|
|
result, coefficients = _run_once(spec, arrays, precision, out_dir / "instrumented",
|
|
source_bytes, trace_queue)
|
|
# Preserve the primary result even if the independent proxy check fails.
|
|
_save_arrays(out_dir / "state.npz", result)
|
|
_save(out_dir / "coefficients.json", coefficients)
|
|
validation = {"performed": False, "pass": None,
|
|
"reason": "Disabled by explicit case switch" if not validate_proxy else None}
|
|
if validate_proxy:
|
|
baseline, baseline_coefficients = _run_once(spec, arrays, precision, out_dir / "unwrapped-audit",
|
|
source_bytes, False)
|
|
checks = {key: bool(np.array_equal(value, baseline[key]))
|
|
for key, value in result.items() if not key.startswith("arrival_")}
|
|
coefficient_keys = ["gd", "c0", "cg", "cv", "threshold", "reset", "seed_volt"]
|
|
checks["coefficients"] = all(coefficients[key] == baseline_coefficients[key] for key in coefficient_keys)
|
|
validation = {"performed": True, "pass": bool(all(checks.values())),
|
|
"exact_arrays": checks,
|
|
"arrival_arrays_excluded": "An unwrapped queue is deliberately not observed"}
|
|
_save(out_dir / "proxy-validation.json", validation)
|
|
_check(not validate_proxy or validation["pass"], "Queue proxy changed an author result")
|
|
_check(_digest(AUTHOR_PATH) == source_hash, "Author source changed during the stress case")
|
|
summary = {"name": spec["name"], "precision": precision, "steps": int(spec["steps"]),
|
|
"nodes": len(arrays["nodes"]), "edges": len(arrays["edge_id"]),
|
|
"spikes": len(result["step"]), "requests": int(result["request_count"].sum()),
|
|
"applicable_requests": int(result["applicable_count"].sum()),
|
|
"rejected_requests": int(result["rejected_count"].sum()),
|
|
"due_arrivals": len(result["arrival_step"]),
|
|
"queue_trace_enabled": trace_queue, "queue_proxy_validated": validation["pass"],
|
|
"source_unchanged": True, "state_file": str(out_dir / "state.npz"),
|
|
"coefficient_file": str(out_dir / "coefficients.json"), "complete": True}
|
|
_save(out_dir / "summary.json", summary)
|
|
return {"arrays": result, "coefficients": coefficients, "summary": summary}
|
|
except Exception as error:
|
|
_save(out_dir / "error.json", {"type": type(error).__name__, "message": str(error)})
|
|
raise
|