Digital Drosophila Archive
316 lines
23 KiB
Python
316 lines
23 KiB
Python
"""CPU-only observer for the unchanged Phase 4AH replay_core callback.
|
|
|
|
All input events are fixed observations, not interventions. This module writes
|
|
no files and never changes a callback array. Incoming records are antecedents,
|
|
not a selected cause, shortest path, necessity test, or whole-brain simulation.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import numpy as np
|
|
|
|
|
|
DELAY = 18
|
|
COEFF_KEYS = ("gd", "c0", "cg", "cv", "threshold", "reset", "seed_volt")
|
|
FLOAT_KEYS = ("old_v", "old_g", "v_integrated", "g_decayed", "v_threshold",
|
|
"v_after_reset", "g_after_reset", "v_pre_apl", "g_pre_apl", "delta_si", "v", "g")
|
|
INT_KEYS = ("ref_before", "ref_after_reset", "ref_pre_apl", "ref")
|
|
BOOL_KEYS = ("available", "request", "request_accepted", "request_rejected", "fired")
|
|
|
|
|
|
def require(condition, message):
|
|
if not condition:
|
|
raise ValueError(message)
|
|
|
|
|
|
def bit_equal(a, b):
|
|
return a.dtype == b.dtype and a.shape == b.shape and np.array_equal(
|
|
np.ascontiguousarray(a).view(np.uint8), np.ascontiguousarray(b).view(np.uint8))
|
|
|
|
|
|
class PNObserver:
|
|
def __init__(self, target_globals, design_target_globals, incoming, edge_groups,
|
|
group_names, coeff, steps, *, max_first_event_rows, event_sink=None):
|
|
targets = np.asarray(target_globals)
|
|
require(targets.dtype == np.int32 and targets.ndim == 1 and len(targets) > 0 and
|
|
np.all((targets >= 0) & (targets < 166700)) and np.all(np.diff(targets) > 0),
|
|
"Targets must be sorted unique global int32 IDs")
|
|
design = np.asarray(design_target_globals)
|
|
require(design.dtype == np.int32 and design.ndim == 1 and len(np.unique(design)) == len(design) and
|
|
np.isin(design, targets).all(), "Design target IDs must be unique int32 subset of observed PNs")
|
|
require(type(steps) is int and steps > 0, "steps must be a positive integer")
|
|
require(type(max_first_event_rows) is int and max_first_event_rows > 0, "Explicit positive row budget required")
|
|
require(event_sink is None or callable(event_sink), "event_sink must be callable or None")
|
|
names = tuple(group_names)
|
|
require(names and all(isinstance(name, str) and name for name in names) and len(set(names)) == len(names),
|
|
"Source group names must be nonempty and unique")
|
|
codes = np.asarray(edge_groups)
|
|
e = len(incoming.edge_id)
|
|
require(codes.ndim == 1 and codes.shape == (e,) and codes.dtype.kind in "iu" and
|
|
np.all((codes >= 0) & (codes < len(names))), "Invalid source-group code per incoming edge")
|
|
require(incoming.target_count == len(targets) and incoming.weights.dtype == np.float64 and
|
|
np.isfinite(incoming.weights).all(), "Incoming index/weights mismatch")
|
|
require(all(name in coeff and np.isfinite(coeff[name]) for name in COEFF_KEYS), "Missing/nonfinite coefficient")
|
|
self.targets = targets.copy()
|
|
self.design = np.isin(targets, design)
|
|
self.incoming, self.groups, self.names = incoming, codes.astype(np.int64, copy=True), names
|
|
self.coeff = {name: np.float64(coeff[name]) for name in COEFF_KEYS}
|
|
self.steps, self.n = steps, len(targets)
|
|
self.max_rows, self.sink = max_first_event_rows, event_sink
|
|
self.callback_count, self.record_count = 0, 0
|
|
self.first_arrival_records = []
|
|
self.first = np.full(self.n, -1, np.int32)
|
|
self.first_request = np.full(self.n, -1, np.int32)
|
|
self.first_accepted_request = np.full(self.n, -1, np.int32)
|
|
self.first_positive_arrival = np.full(self.n, -1, np.int32)
|
|
self.counts = np.zeros(self.n, np.int32)
|
|
self.requested = np.zeros(self.n, np.int64)
|
|
self.request_accepted = np.zeros(self.n, np.int64)
|
|
self.request_rejected = np.zeros(self.n, np.int64)
|
|
self.first_request_coincident = np.zeros(self.n, bool)
|
|
self.first_state = {name: np.full(self.n, np.nan, np.float64) for name in
|
|
("old_v", "old_g", "v_integrated", "v_threshold", "threshold_margin", "ref_before")}
|
|
shape = (len(names), self.n)
|
|
self.stats = {name: np.zeros(shape, np.int64) for name in
|
|
("due_edges", "accepted_edges", "rejected_edges", "accepted_zero_edges", "rejected_zero_edges",
|
|
"accepted_positive_edges", "accepted_negative_edges", "rejected_positive_edges", "rejected_negative_edges",
|
|
"accepted_nonzero_before_first_edges", "accepted_positive_before_first_edges", "accepted_negative_before_first_edges")}
|
|
self.sums = {name: np.zeros(shape, np.longdouble) for name in
|
|
("accepted_positive_V", "accepted_negative_abs_V", "rejected_positive_V", "rejected_negative_abs_V",
|
|
"accepted_positive_before_first_V", "accepted_negative_before_first_abs_V")}
|
|
self.edge_prefirst_count = np.zeros(e, np.int64)
|
|
self.edge_first_arrival = np.full(e, -1, np.int32)
|
|
self.edge_last_arrival = np.full(e, -1, np.int32)
|
|
# Non-feedback signed-source decomposition. Residual corrections are
|
|
# explicit lanes, never silently assigned to an anatomical population.
|
|
self.lane_names = [name + "/positive_weight" for name in names] + [name + "/negative_weight" for name in names]
|
|
self.lane_names += ["initial_and_affine_baseline", "request_or_reset", "roundoff_residual"]
|
|
self.baseline_lane, self.request_lane, self.round_lane = 2 * len(names), 2 * len(names) + 1, 2 * len(names) + 2
|
|
self.lv = np.zeros((len(self.lane_names), self.n), np.longdouble)
|
|
self.lg = np.zeros_like(self.lv)
|
|
self.first_integrated_lanes = np.full_like(self.lv, np.nan)
|
|
self.first_threshold_lanes = np.full_like(self.lv, np.nan)
|
|
self.first_ledger_residual = np.full(self.n, np.nan, np.longdouble)
|
|
self.integration_roundoff_v = np.zeros(self.n, np.longdouble)
|
|
self.integration_roundoff_g = np.zeros(self.n, np.longdouble)
|
|
self.arrival_roundoff_g = np.zeros(self.n, np.longdouble)
|
|
self.max_ledger_residual = np.longdouble(0)
|
|
self.previous = None
|
|
self.finished = False
|
|
|
|
def _validate_stage(self, tick, s):
|
|
require(not self.finished and type(tick) in (int, np.int32, np.int64) and
|
|
int(tick) == self.callback_count and tick < self.steps, "Missing, duplicate or out-of-order callback tick")
|
|
for names, dtype in ((FLOAT_KEYS, np.float64), (INT_KEYS, np.int32), (BOOL_KEYS, bool)):
|
|
for name in names:
|
|
require(name in s and isinstance(s[name], np.ndarray) and s[name].shape == (self.n,) and
|
|
s[name].dtype == dtype, "Missing/wrong callback array: " + name)
|
|
require(np.isfinite(s[name]).all(), "Nonfinite callback array: " + name)
|
|
require("h32" in s and s["h32"].dtype == np.float32 and s["h32"].shape == (self.n,) and
|
|
np.all(s["h32"] == 0) and np.all(s["delta_si"] == 0), "PN observer requires zero graded APL h/delta")
|
|
require(bit_equal(s["g"], s["g_pre_apl"]) and bit_equal(s["v"], s["v_pre_apl"]) and
|
|
bit_equal(s["ref"], s["ref_pre_apl"]), "Nonzero or unexpected post-APL PN state change")
|
|
require(np.all((s["ref_before"] >= 0) & (s["ref_before"] <= 22)) and
|
|
np.all((s["ref"] >= 0) & (s["ref"] <= 22)), "Refractory range mismatch")
|
|
require(np.array_equal(s["available"], s["ref_before"] <= 1), "Availability does not match pre-step refractory")
|
|
require(np.array_equal(s["request_accepted"], s["request"] & s["available"]) and
|
|
np.array_equal(s["request_rejected"], s["request"] & ~s["available"]), "Request accounting mismatch")
|
|
require(not np.any(s["request"] & ~self.design), "Observed request outside frozen odor design targets")
|
|
require(np.array_equal(s["fired"], s["available"] & (s["v_threshold"] > self.coeff["threshold"])),
|
|
"Callback firing does not match strict threshold")
|
|
if self.previous is not None:
|
|
require(all(bit_equal(s[key], self.previous[old]) for key, old in
|
|
(("old_v", "v"), ("old_g", "g"), ("ref_before", "ref"))), "Discontinuous callback state")
|
|
due = s.get("due_edge_indices")
|
|
require(isinstance(due, np.ndarray) and due.dtype == np.int64 and due.ndim == 1 and
|
|
np.all((due >= 0) & (due < len(self.incoming.edge_id))) and np.all(np.diff(due) > 0),
|
|
"Incoming queue must retain unique canonical edge rows")
|
|
require("due_edge_id" in s and bit_equal(s["due_edge_id"], self.incoming.edge_id[due]), "Incoming edge-ID mismatch")
|
|
accepted = s.get("arrival_accepted")
|
|
require(isinstance(accepted, np.ndarray) and accepted.dtype == bool and accepted.shape == due.shape and
|
|
np.array_equal(accepted, s["ref_after_reset"][self.incoming.dst_local[due]] == 0), "Arrival guard mismatch")
|
|
require(tick >= DELAY or not len(due), "Incoming event before the fixed 18-step delay")
|
|
return due, accepted
|
|
|
|
def __call__(self, tick, s):
|
|
due, accepted = self._validate_stage(tick, s)
|
|
tick = int(tick)
|
|
dst, weights, groups = self.incoming.dst_local[due], self.incoming.weights[due], self.groups[due]
|
|
before = self.first < 0
|
|
first_now = before & s["fired"]
|
|
# Arrivals happen AFTER this tick's threshold. Same-tick arrivals are
|
|
# never antecedents to a first spike, even if malformed data said accepted.
|
|
strictly_before = before[dst] & ~s["fired"][dst]
|
|
record_mask = accepted & (weights != 0) & strictly_before
|
|
new_records = int(np.count_nonzero(record_mask))
|
|
require(self.record_count + new_records <= self.max_rows,
|
|
"First-arrival record resource limit exceeded; no truncation or event selection permitted")
|
|
if self.previous is None:
|
|
self.lv[self.baseline_lane] = s["old_v"].astype(np.longdouble)
|
|
self.lg[self.baseline_lane] = s["old_g"].astype(np.longdouble)
|
|
available = s["available"]
|
|
# Use old lane-g for V, as in the author update; do NOT use decayed g.
|
|
self.lv[:, available] = (np.longdouble(self.coeff["cv"]) * self.lv[:, available] +
|
|
np.longdouble(self.coeff["cg"]) * self.lg[:, available])
|
|
self.lv[self.baseline_lane, available] += np.longdouble(self.coeff["c0"])
|
|
self.lg[:, available] *= np.longdouble(self.coeff["gd"])
|
|
rv = s["v_integrated"].astype(np.longdouble) - self.lv.sum(axis=0, dtype=np.longdouble)
|
|
rg = s["g_decayed"].astype(np.longdouble) - self.lg.sum(axis=0, dtype=np.longdouble)
|
|
self.integration_roundoff_v += rv
|
|
self.integration_roundoff_g += rg
|
|
self.lv[self.round_lane] += rv
|
|
self.lg[self.round_lane] += rg
|
|
self.first_integrated_lanes[:, first_now] = self.lv[:, first_now]
|
|
req = s["request_accepted"]
|
|
self.lv[:, req] = 0
|
|
self.lv[self.request_lane, req] = np.longdouble(self.coeff["seed_volt"])
|
|
self.first_threshold_lanes[:, first_now] = self.lv[:, first_now]
|
|
vres = s["v_threshold"].astype(np.longdouble) - self.lv.sum(axis=0, dtype=np.longdouble)
|
|
self.first_ledger_residual[first_now] = vres[first_now]
|
|
self.max_ledger_residual = max(self.max_ledger_residual, np.max(np.abs(vres), initial=np.longdouble(0)))
|
|
for name in self.first_state:
|
|
values = (self.coeff["threshold"] - s["v_threshold"]) if name == "threshold_margin" else s[name]
|
|
self.first_state[name][first_now] = values[first_now]
|
|
self.first_request_coincident[first_now] = req[first_now]
|
|
self.first[first_now] = tick
|
|
for array, selected in ((self.first_request, s["request"]), (self.first_accepted_request, req)):
|
|
array[(array < 0) & selected] = tick
|
|
self.requested += s["request"].astype(np.int64)
|
|
self.request_accepted += req.astype(np.int64)
|
|
self.request_rejected += s["request_rejected"].astype(np.int64)
|
|
self.counts += s["fired"].astype(np.int32)
|
|
fired = s["fired"]
|
|
self.lv[:, fired] = 0
|
|
self.lv[self.request_lane, fired] = np.longdouble(self.coeff["reset"])
|
|
self.lg[:, fired] = 0
|
|
positive, negative, zero = weights > 0, weights < 0, weights == 0
|
|
picks = {"due_edges": np.ones(len(due), bool), "accepted_edges": accepted,
|
|
"rejected_edges": ~accepted, "accepted_zero_edges": accepted & zero,
|
|
"rejected_zero_edges": ~accepted & zero,
|
|
"accepted_positive_edges": accepted & positive, "accepted_negative_edges": accepted & negative,
|
|
"rejected_positive_edges": ~accepted & positive, "rejected_negative_edges": ~accepted & negative,
|
|
"accepted_nonzero_before_first_edges": record_mask,
|
|
"accepted_positive_before_first_edges": record_mask & positive,
|
|
"accepted_negative_before_first_edges": record_mask & negative}
|
|
for name, mask in picks.items():
|
|
np.add.at(self.stats[name], (groups[mask], dst[mask]), 1)
|
|
for name, mask, sign in (
|
|
("accepted_positive_V", accepted & positive, 1), ("accepted_negative_abs_V", accepted & negative, -1),
|
|
("rejected_positive_V", ~accepted & positive, 1), ("rejected_negative_abs_V", ~accepted & negative, -1),
|
|
("accepted_positive_before_first_V", record_mask & positive, 1),
|
|
("accepted_negative_before_first_abs_V", record_mask & negative, -1)):
|
|
np.add.at(self.sums[name], (groups[mask], dst[mask]), sign * weights[mask].astype(np.longdouble))
|
|
for offset, mask in ((0, accepted & positive), (len(self.names), accepted & negative)):
|
|
np.add.at(self.lg, (groups[mask] + offset, dst[mask]), weights[mask].astype(np.longdouble))
|
|
ag = s["g_pre_apl"].astype(np.longdouble) - self.lg.sum(axis=0, dtype=np.longdouble)
|
|
self.arrival_roundoff_g += ag
|
|
self.lg[self.round_lane] += ag
|
|
gres = s["g"].astype(np.longdouble) - self.lg.sum(axis=0, dtype=np.longdouble)
|
|
self.max_ledger_residual = max(self.max_ledger_residual, np.max(np.abs(gres), initial=np.longdouble(0)))
|
|
require(np.isfinite(self.lv).all() and np.isfinite(self.lg).all(), "Nonfinite non-feedback lineage ledger")
|
|
positive_targets = np.unique(dst[accepted & positive])
|
|
self.first_positive_arrival[positive_targets[self.first_positive_arrival[positive_targets] < 0]] = tick
|
|
selected_rows = due[record_mask]
|
|
np.add.at(self.edge_prefirst_count, selected_rows, 1)
|
|
self.edge_first_arrival[selected_rows[self.edge_first_arrival[selected_rows] < 0]] = tick
|
|
self.edge_last_arrival[selected_rows] = tick
|
|
rows = []
|
|
for edge in selected_rows:
|
|
target = int(self.incoming.dst_local[edge])
|
|
weight = float(self.incoming.weights[edge])
|
|
rows.append({"edge_id": int(self.incoming.edge_id[edge]), "incoming_edge_row": int(edge),
|
|
"source_neuron_idx": int(self.incoming.src_global[edge]), "target_neuron_idx": int(self.targets[target]),
|
|
"source_group": self.names[int(self.groups[edge])], "source_spike_tick": tick - DELAY,
|
|
"arrival_tick": tick, "target_is_design": bool(self.design[target]), "accepted": True,
|
|
"weight_SI": weight, "weight_sign": 1 if weight > 0 else -1,
|
|
"ref_before": int(s["ref_before"][target]), "ref_at_arrival": int(s["ref_after_reset"][target]),
|
|
"available_at_integration": bool(s["available"][target]),
|
|
"old_v_SI": float(s["old_v"][target]), "old_g_SI": float(s["old_g"][target]),
|
|
"v_integrated_SI": float(s["v_integrated"][target]), "v_prethreshold_SI": float(s["v_threshold"][target]),
|
|
"threshold_margin_SI": float(self.coeff["threshold"] - s["v_threshold"][target]),
|
|
"g_before_arrivals_SI": float(s["g_after_reset"][target]),
|
|
"record_scope": "accepted nonzero input strictly before any observed first spike; join final first_spike_tick"})
|
|
if rows:
|
|
if self.sink is None:
|
|
self.first_arrival_records.extend(rows)
|
|
else:
|
|
self.sink(rows)
|
|
self.record_count += len(rows)
|
|
self.previous = {name: s[name].copy() for name in ("v", "g", "ref")}
|
|
self.callback_count += 1
|
|
|
|
def edge_summary(self):
|
|
"""All edges with stored pre-first records, including never-active PNs."""
|
|
selected = np.flatnonzero(self.edge_prefirst_count)
|
|
dst = self.incoming.dst_local[selected]
|
|
return {"edge_id": self.incoming.edge_id[selected].copy(),
|
|
"source_neuron_idx": self.incoming.src_global[selected].copy(), "target_neuron_idx": self.targets[dst].copy(),
|
|
"source_group": np.asarray(self.names, dtype=str)[self.groups[selected]],
|
|
"weight_SI": self.incoming.weights[selected].copy(),
|
|
"accepted_event_count": self.edge_prefirst_count[selected].copy(),
|
|
"nominal_weight_times_event_count_SI": (self.incoming.weights[selected].astype(np.longdouble) *
|
|
self.edge_prefirst_count[selected]).astype(np.float64),
|
|
"first_arrival_tick": self.edge_first_arrival[selected].copy(),
|
|
"last_arrival_tick": self.edge_last_arrival[selected].copy(),
|
|
"target_first_spike_tick": self.first[dst].copy(), "target_is_design": self.design[dst].copy()}
|
|
|
|
def finish(self, result):
|
|
"""Return (per-PN arrays, JSON report); caller verifies original events/state.
|
|
|
|
first_* floating fields are NaN when first_spike_observed is false.
|
|
They must not be imputed to zero or interpreted as observed voltages.
|
|
"""
|
|
require(not self.finished and self.callback_count == self.steps and self.previous is not None, "Observer did not cover every tick")
|
|
require(bit_equal(result["target_globals"], self.targets) and bit_equal(result["first_spike"], self.first) and
|
|
bit_equal(result["counts"], self.counts), "Replay result identity/count/first spike differs from observer")
|
|
require(all(bit_equal(result[name], self.previous[name]) for name in ("v", "g", "ref")), "Replay final state differs from last callback")
|
|
require(int(self.edge_prefirst_count.sum()) == self.record_count ==
|
|
int(self.stats["accepted_nonzero_before_first_edges"].sum()), "Pre-first event coverage mismatch")
|
|
require(np.array_equal(self.stats["due_edges"], self.stats["accepted_edges"] + self.stats["rejected_edges"]),
|
|
"Due/accepted/rejected coverage mismatch")
|
|
data = {"neuron_idx": self.targets.copy(), "is_design_target": self.design.copy(),
|
|
"spike_count": self.counts.copy(), "first_spike_tick": self.first.copy(), "first_spike_observed": self.first >= 0,
|
|
"first_request_tick": self.first_request.copy(), "first_accepted_request_tick": self.first_accepted_request.copy(),
|
|
"first_accepted_positive_arrival_tick": self.first_positive_arrival.copy(),
|
|
"requested": self.requested.copy(), "accepted_requests": self.request_accepted.copy(),
|
|
"rejected_requests": self.request_rejected.copy(), "first_spike_coincident_with_accepted_request": self.first_request_coincident.copy(),
|
|
"first_voltage_ledger_residual_SI": self.first_ledger_residual.astype(np.float64),
|
|
"integration_roundoff_v_SI": self.integration_roundoff_v.astype(np.float64),
|
|
"integration_roundoff_g_SI": self.integration_roundoff_g.astype(np.float64),
|
|
"arrival_roundoff_g_SI": self.arrival_roundoff_g.astype(np.float64)}
|
|
data.update({"first_" + name + ("_SI" if name != "ref_before" else ""): value.copy()
|
|
for name, value in self.first_state.items()})
|
|
for gi, name in enumerate(self.names):
|
|
for key, value in {**self.stats, **self.sums}.items():
|
|
data[name + "_" + key] = value[gi].copy() if key in self.stats else value[gi].astype(np.float64)
|
|
for index, name in enumerate(self.lane_names):
|
|
data["first_integrated_V_from_" + name] = self.first_integrated_lanes[index].astype(np.float64)
|
|
data["first_prethreshold_V_from_" + name] = self.first_threshold_lanes[index].astype(np.float64)
|
|
aggregates = {}
|
|
for name, mask in (("all_PN", np.ones(self.n, bool)), ("design_PN", self.design),
|
|
("non_design_PN", ~self.design), ("never_active_PN", self.first < 0)):
|
|
aggregates[name] = {"neurons": int(mask.sum()), "active": int(np.count_nonzero(mask & (self.first >= 0))),
|
|
"spikes": int(self.counts[mask].sum()), "source_groups": {
|
|
group: {**{key: int(value[gi, mask].sum()) for key, value in self.stats.items()},
|
|
**{key: float(value[gi, mask].sum(dtype=np.longdouble)) for key, value in self.sums.items()}}
|
|
for gi, group in enumerate(self.names)}}
|
|
report = {"pass_all": True, "scope": "Observer coverage/internal consistency only; caller must verify recorded PN events and final state exactly",
|
|
"callback_count": self.callback_count, "steps": self.steps, "PN_count": self.n,
|
|
"source_group_names": list(self.names), "first_arrival_record_count": self.record_count,
|
|
"first_arrival_record_limit": self.max_rows, "records_truncated": False,
|
|
"edge_summary_count": int(np.count_nonzero(self.edge_prefirst_count)), "aggregates": aggregates,
|
|
"time_contract": "source_spike_tick=arrival_tick-18; arrivals follow threshold; first-spike tick arrivals excluded",
|
|
"record_contract": "Every accepted nonzero edge event strictly before first spike for every PN; never-active PNs retain full-window records with final first=-1",
|
|
"sum_contract": "Observer nominal signed weight sums, not the sequentially rounded runtime g; weight-times-count edge summaries are descriptive",
|
|
"voltage_decomposition": {"nonfeedback": True, "precision": str(np.dtype(np.longdouble)),
|
|
"lanes": list(self.lane_names), "maximum_reconstruction_residual_SI": float(self.max_ledger_residual),
|
|
"roundoff_lane_explicit": True, "fitted_tolerance": False,
|
|
"reset_clears_prior_lanes": True, "request_override_has_separate_lane": True,
|
|
"first_integrated_is_before_request_override": True,
|
|
"never_spiked_first_fields_are_NaN": True,
|
|
"threshold_crossing_contributions_are_not_necessity_or_unique_cause": True},
|
|
"graded_APL_to_PN": "Required h32=0 and delta_si=0 every tick; fast APL weights supplied already frozen-masked by caller",
|
|
"new_brain_run": False, "parameter_changes": False, "shortest_path_claim": False,
|
|
"unique_cause_claim": False, "whole_brain_intervention_claim": False}
|
|
self.finished = True
|
|
return data, report
|