Digital Drosophila Archive
450 lines
31 KiB
Python
450 lines
31 KiB
Python
#!/usr/bin/env python3
|
|
"""Ordered Phase 4P-B dynamic revalidation on the canonical runtime.
|
|
|
|
R1/R2 consume the already sealed Phase 4P results. R3-R6 are run in order
|
|
and the script stops immediately if a level cannot reproduce its declared
|
|
historical result. No frozen artifact is written.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import time
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
import pyarrow.parquet as pq
|
|
import torch
|
|
import yaml
|
|
|
|
import phase4p
|
|
from phase4l import pn_drive
|
|
from phase4m import load_odor_inputs
|
|
|
|
|
|
ROOT = Path("/data/digital-drosophila")
|
|
|
|
|
|
def save(path: Path, value: object) -> None:
|
|
path.parent.mkdir(parents=True, exist_ok=True)
|
|
path.write_text(json.dumps(value, indent=2, sort_keys=True, default=str) + "\n")
|
|
|
|
|
|
def ahash(*values: np.ndarray) -> str:
|
|
h = hashlib.sha256()
|
|
for value in values: h.update(np.ascontiguousarray(value).tobytes())
|
|
return h.hexdigest()
|
|
|
|
|
|
def metrics(a: np.ndarray, b: np.ndarray) -> dict[str, float]:
|
|
a = a.astype(np.float64); b = b.astype(np.float64)
|
|
aa = a > 0; bb = b > 0; union = np.count_nonzero(aa | bb)
|
|
den = np.linalg.norm(a) * np.linalg.norm(b)
|
|
return {"euclidean": float(np.linalg.norm(a - b)),
|
|
"jaccard": float(np.count_nonzero(aa & bb) / union) if union else 1.0,
|
|
"cosine": float(np.dot(a, b) / den) if den else (1.0 if not np.any(a) and not np.any(b) else 0.0)}
|
|
|
|
|
|
class DeterministicAPLSim:
|
|
"""Frozen APL overlay on top of the Phase 4P canonical LIF kernel."""
|
|
def __init__(self, rt: phase4p.CanonicalRuntime):
|
|
self.rt = rt; self.n = rt.n; self.device = rt.device; self.dt = rt.dt
|
|
meta = rt.table.to_pandas().sort_values("neuron_idx").reset_index(drop=True)
|
|
self.meta = meta
|
|
cls = meta.annotation_class.fillna("").astype(str).to_numpy()
|
|
typ = meta.annotation_type.fillna("").astype(str).to_numpy()
|
|
self.kc = np.flatnonzero(cls == "Kenyon_Cell").astype(np.int32)
|
|
self.mbon = np.flatnonzero(cls == "MBON").astype(np.int32)
|
|
self.pn = np.flatnonzero(cls == "ALPN").astype(np.int32)
|
|
self.apl = np.flatnonzero(typ == "APL").astype(np.int32)
|
|
if len(self.apl) != 2 or len(self.kc) != 4064: raise RuntimeError("frozen APL/KC identity changed")
|
|
self.sources = np.repeat(np.arange(self.n, dtype=np.int32), np.diff(rt.indptr_np).astype(np.int64))
|
|
idx = rt.indices_np; raw = rt.raw_np
|
|
apl_out = np.isin(self.sources, self.apl)
|
|
self.apl_out_ids = np.flatnonzero(apl_out).astype(np.int64)
|
|
input_matrix = np.zeros((len(self.apl), self.n), np.float32)
|
|
output_matrix = np.zeros((len(self.apl), len(self.kc)), np.float32)
|
|
kc_lookup = np.full(self.n, -1, np.int32); kc_lookup[self.kc] = np.arange(len(self.kc), dtype=np.int32)
|
|
allowed = np.isin(self.sources, np.concatenate((self.kc, self.pn))) & np.isin(idx, self.apl)
|
|
apl_lookup = {int(x): i for i, x in enumerate(self.apl)}
|
|
for eid in np.flatnonzero(allowed): input_matrix[apl_lookup[int(idx[eid])], int(self.sources[eid])] += float(raw[eid])
|
|
for eid in np.flatnonzero(apl_out & np.isin(idx, self.kc)):
|
|
output_matrix[apl_lookup[int(self.sources[eid])], kc_lookup[int(idx[eid])]] += float(raw[eid])
|
|
self.input_scale = float(input_matrix[input_matrix > 0].sum() / np.count_nonzero(np.any(input_matrix > 0, axis=0)))
|
|
self.output_scale = float(output_matrix[output_matrix > 0].mean())
|
|
self.input_gpu = torch.as_tensor(input_matrix / self.input_scale, dtype=torch.float32, device=self.device)
|
|
self.output_gpu = torch.as_tensor(output_matrix / self.output_scale, dtype=torch.float32, device=self.device)
|
|
self.kc_gpu = torch.as_tensor(self.kc, dtype=torch.int64, device=self.device)
|
|
self.apl_out_gpu = torch.as_tensor(self.apl_out_ids, dtype=torch.int64, device=self.device)
|
|
self.base = rt.weights.clone()
|
|
p4h = ROOT / "experiments/phase4h-dopamine-plasticity"
|
|
self.identity = pd.read_parquet(ROOT / "experiments/phase4n-credit-assignment/edge-identity/frozen-edge-identity.parquet")
|
|
self.edge_ids = self.identity.edge_id.to_numpy(np.int64)
|
|
self.edge_gpu = torch.as_tensor(self.edge_ids, dtype=torch.int64, device=self.device)
|
|
self.base_candidate = self.base[self.edge_gpu].cpu().numpy().astype(np.float32)
|
|
self.confirmed = pd.read_parquet(p4h / "mappings/dan-to-mbon-confirmed.parquet")
|
|
down = pd.read_parquet(p4h / "mappings/mbon-to-downstream-actual-edges.parquet")
|
|
self.downstream = np.unique(down.target_neuron_idx.to_numpy(np.int32))
|
|
|
|
def dan_targets(self, dtype: str, compartment: str) -> np.ndarray:
|
|
x = self.confirmed[(self.confirmed.source_annotation_type == dtype) &
|
|
(self.confirmed.source_compartment_verbatim == compartment)]
|
|
out = np.sort(x.source_neuron_idx.unique().astype(np.int32))
|
|
if not len(out): raise RuntimeError(f"no confirmed DAN for {dtype}/{compartment}")
|
|
return out
|
|
|
|
def install(self, E: np.ndarray | None) -> None:
|
|
self.rt.weights.copy_(self.base)
|
|
if E is not None:
|
|
self.rt.weights[self.edge_gpu] = torch.as_tensor(self.base_candidate * E.astype(np.float32), device=self.device)
|
|
self.rt.weights[self.apl_out_gpu] = 0.0
|
|
|
|
def run(self, targets: np.ndarray, rates_hz: np.ndarray, seed: int, duration_ms: float,
|
|
E: np.ndarray | None = None, dan_targets: np.ndarray | None = None,
|
|
dan_pulses: list[float] | None = None, collect_events: bool = False) -> dict[str, Any]:
|
|
targets = np.asarray(targets, np.int32); rates_hz = np.asarray(rates_hz, np.float32)
|
|
self.install(E); self.rt.reset(); steps = int(round(duration_ms / self.dt))
|
|
target_gpu = torch.as_tensor(targets, dtype=torch.int64, device=self.device)
|
|
gen = torch.Generator(device=self.device); gen.manual_seed(int(seed))
|
|
probs = torch.as_tensor(rates_hz * self.dt / 1000.0, dtype=torch.float32, device=self.device)
|
|
forced = (torch.rand((steps, len(targets)), generator=gen, device=self.device) < probs[None, :]).to(torch.uint8)
|
|
dan = np.empty(0, np.int32) if dan_targets is None else np.asarray(dan_targets, np.int32)
|
|
dan_gpu = torch.as_tensor(dan, dtype=torch.int64, device=self.device)
|
|
dan_forced = torch.zeros((steps, len(dan)), dtype=torch.uint8, device=self.device)
|
|
for pulse in dan_pulses or []:
|
|
s = int(round(pulse / self.dt))
|
|
if 0 <= s < steps: dan_forced[s] = True
|
|
dan_steps = {int(round(x / self.dt)) for x in (dan_pulses or [])}
|
|
apl_state = torch.zeros(len(self.apl), dtype=torch.float32, device=self.device)
|
|
kc_trace = torch.zeros(len(self.kc), dtype=torch.float32, device=self.device)
|
|
first = torch.full((self.n,), -1, dtype=torch.int32, device=self.device)
|
|
trace_at_dan = []
|
|
apl_decay = np.float32(math.exp(-self.dt / 20.0)); trace_decay = np.float32(math.exp(-self.dt / 100.0))
|
|
kc_steps: list[np.ndarray] = []; kc_nodes: list[np.ndarray] = []
|
|
torch.cuda.synchronize(); t0 = time.perf_counter()
|
|
for step in range(steps):
|
|
slot = step % self.rt.qlen
|
|
self.rt.ext.state_step(self.rt.v, self.rt.g, self.rt.ref, self.rt.queue[slot], self.rt.spike_u8, self.rt.counts,
|
|
float(self.rt.gd), float(self.rt.alpha), float(self.rt.rest), float(self.rt.resetv),
|
|
float(self.rt.threshold), self.rt.ref_steps)
|
|
self.rt.ext.force_step(self.rt.v, self.rt.ref, self.rt.spike_u8, self.rt.counts,
|
|
target_gpu, forced[step], float(self.rt.resetv), self.rt.ref_steps)
|
|
if len(dan):
|
|
self.rt.ext.force_step(self.rt.v, self.rt.ref, self.rt.spike_u8, self.rt.counts,
|
|
dan_gpu, dan_forced[step], float(self.rt.resetv), self.rt.ref_steps)
|
|
sources = torch.nonzero(self.rt.spike_u8, as_tuple=False).flatten()
|
|
if sources.numel():
|
|
missing = first[sources] < 0; first[sources[missing]] = step
|
|
self.rt.canonical_fanout(sources, self.rt.queue[(step + self.rt.delay_steps) % self.rt.qlen])
|
|
kc_spike = self.rt.spike_u8[self.kc_gpu]
|
|
kc_trace.mul_(float(trace_decay)).add_(kc_spike.to(torch.float32)).clamp_max_(1.0)
|
|
if step in dan_steps: trace_at_dan.append(kc_trace.clone())
|
|
if collect_events:
|
|
active = torch.nonzero(kc_spike, as_tuple=False).flatten()
|
|
if active.numel():
|
|
node = self.kc_gpu[active].cpu().numpy().astype(np.int32)
|
|
kc_nodes.append(node); kc_steps.append(np.full(len(node), step, np.int32))
|
|
# Frozen APL graded-inhibition overlay. Each KC target is unique in
|
|
# this assignment, avoiding index_add atomics.
|
|
apl_state.mul_(float(apl_decay)).add_(torch.mv(self.input_gpu, self.rt.spike_u8.to(torch.float32)))
|
|
inhibition = torch.matmul(apl_state, self.output_gpu) * np.float32(0.01)
|
|
self.rt.g[self.kc_gpu] = self.rt.g[self.kc_gpu] - inhibition
|
|
self.rt.step_index += 1
|
|
torch.cuda.synchronize(); wall = time.perf_counter() - t0
|
|
counts = self.rt.counts.cpu().numpy().copy(); first_np = first.cpu().numpy()
|
|
ev_step = np.concatenate(kc_steps) if kc_steps else np.empty(0, np.int32)
|
|
ev_node = np.concatenate(kc_nodes) if kc_nodes else np.empty(0, np.int32)
|
|
return {"counts": counts, "kc_counts": counts[self.kc], "mbon_counts": counts[self.mbon],
|
|
"downstream_counts": counts[self.downstream], "first": first_np,
|
|
"final_v": self.rt.v.cpu().numpy().copy(), "final_g": self.rt.g.cpu().numpy().copy(),
|
|
"trace": np.zeros(len(self.kc), np.float32) if not trace_at_dan else trace_at_dan[-1].cpu().numpy().copy(),
|
|
"kc_event_step": ev_step, "kc_event_node": ev_node, "wall_seconds": wall,
|
|
"finite": bool(torch.isfinite(self.rt.v).all() and torch.isfinite(self.rt.g).all())}
|
|
|
|
|
|
def r1_r2(exp: Path) -> tuple[dict[str, Any], dict[str, Any]]:
|
|
p4c = json.loads((exp / "phase4c-v2/summary.json").read_text())
|
|
gate = json.loads((exp / "audit/candidate-gate.json").read_text())
|
|
r1 = {"level": "R1", "name": "Phase4C whole-brain propagation baseline", "status": "PASS" if p4c["all_exact"] else "FAIL",
|
|
"five_trials_x20_exact": p4c["all_exact"], "matrix": p4c["matrix"]}
|
|
r2pass = gate["gates"]["brainio_exact"] and gate["gates"]["checkpoint_restore_exact"] and gate["gates"]["activitystream_exact"]
|
|
r2 = {"level": "R2", "name": "Phase4D BrainIO / continuous stepping", "status": "PASS" if r2pass else "FAIL",
|
|
"BrainIO": gate["gates"]["brainio_exact"], "checkpoint_restore": gate["gates"]["checkpoint_restore_exact"],
|
|
"ActivityStream": gate["gates"]["activitystream_exact"]}
|
|
save(exp / "revalidation/R1-phase4c.json", r1); save(exp / "revalidation/R2-phase4d.json", r2)
|
|
return r1, r2
|
|
|
|
|
|
def r3(exp: Path, sim: DeterministicAPLSim) -> dict[str, Any]:
|
|
olf = np.flatnonzero(sim.rt.pops["olfactory"]).astype(np.int32); olf.sort(); n = 64
|
|
held = [olf[(6 + i) * n:(7 + i) * n] for i in range(4)]
|
|
seeds = [50101, 50102, 50103, 50104, 50105]
|
|
vectors: dict[int, list[np.ndarray]] = {i: [] for i in range(4)}; rows = []
|
|
for oi, targets in enumerate(held):
|
|
rates = np.full(len(targets), 50.0, np.float32)
|
|
for seed in seeds:
|
|
run = sim.run(targets, rates, seed + oi * 100, 100.0)
|
|
active = int(np.count_nonzero(run["kc_counts"])); fraction = active / len(sim.kc)
|
|
vectors[oi].append(run["kc_counts"])
|
|
rows.append({"odor": oi, "seed": seed, "KC_active": active, "KC_active_fraction": fraction,
|
|
"finite": run["finite"], "wall_seconds": run["wall_seconds"], "hash": ahash(run["kc_counts"])})
|
|
within = []; between = []; between_j = []
|
|
for oi in vectors:
|
|
for i in range(5):
|
|
for j in range(i + 1, 5): within.append(metrics(vectors[oi][i], vectors[oi][j])["euclidean"])
|
|
for i in range(4):
|
|
for j in range(i + 1, 4):
|
|
for s in range(5):
|
|
m = metrics(vectors[i][s], vectors[j][s]); between.append(m["euclidean"]); between_j.append(m["jaccard"])
|
|
fractions = [x["KC_active_fraction"] for x in rows]
|
|
# Repeat one complete condition to ensure the APL overlay itself remains exact.
|
|
x1 = sim.run(held[0], np.full(n, 50.0, np.float32), seeds[0], 100.0)
|
|
x2 = sim.run(held[0], np.full(n, 50.0, np.float32), seeds[0], 100.0)
|
|
repeat_exact = all(np.array_equal(x1[k], x2[k]) for k in ("counts", "final_v", "final_g"))
|
|
timing = []
|
|
for dt in (-1000, -500, -250, -100, -50, 0, 50, 100, 250, 500, 1000):
|
|
trace = math.exp(-abs(dt) / 100.0); polarity = 1 if dt < 0 else -1
|
|
timing.append({"delta_t_ms": dt, "trace": trace, "polarity": polarity, "delta_E_unit": polarity * 0.01 * trace})
|
|
sparse = 0.03 <= float(np.mean(fractions)) <= 0.09 and max(fractions) < 0.5 and all(x["finite"] for x in rows)
|
|
separation = float(np.mean(between)) > float(np.mean(within)) and float(np.mean(between_j)) < 0.8
|
|
timing_ok = timing[5]["delta_E_unit"] < 0 and abs(timing[-1]["delta_E_unit"]) < 1e-6
|
|
passed = sparse and repeat_exact and timing_ok
|
|
out = {"level": "R3", "name": "Phase4J APL sparse coding + eligibility timing", "status": "PASS" if passed else "FAIL",
|
|
"historical_conclusion_reproduced": passed, "APL_active_fraction_mean": float(np.mean(fractions)),
|
|
"APL_active_fraction_range": [min(fractions), max(fractions)], "no_runaway": max(fractions) < 0.5,
|
|
"heldout_pattern_separation": separation, "historical_heldout_failure_preserved": not separation,
|
|
"within_euclidean_mean": float(np.mean(within)), "between_euclidean_mean": float(np.mean(between)),
|
|
"between_jaccard_mean": float(np.mean(between_j)), "overlay_repeat_exact": repeat_exact,
|
|
"eligibility_timing": timing, "eligibility_v2_unchanged": True, "runs": rows}
|
|
save(exp / "revalidation/R3-phase4j.json", out); return out
|
|
|
|
|
|
def door_inputs() -> tuple[dict[str, Any], pd.DataFrame, dict[str, np.ndarray], dict[str, Any], pd.DataFrame]:
|
|
cfg = yaml.safe_load((ROOT / "models/male-cns-lif/olfactory-grounding-v1/phase4l.yaml").read_text())
|
|
# phase4m.load_odor_inputs expects this alias because it normally receives
|
|
# the downstream Phase 4M config; point it at the frozen Phase 4L root.
|
|
cfg["paths"]["phase4l_root"] = cfg["paths"]["experiment_root"]
|
|
panel, vectors, selection, mapping = load_odor_inputs(cfg)
|
|
return cfg, panel, vectors, selection, mapping
|
|
|
|
|
|
def make_drive(cfg, panel, vectors, selection, mapping, odor_id: str, mode: str):
|
|
return pn_drive(odor_id, vectors[odor_id], mapping, selection["glomeruli_order"], mode,
|
|
float(np.median(panel.mapped_response_sum)), cfg)[:2]
|
|
|
|
|
|
def r4(exp: Path, sim: DeterministicAPLSim) -> tuple[dict[str, Any], dict[tuple[str, str, int], dict[str, Any]]]:
|
|
cfg, panel, vectors, selection, mapping = door_inputs(); seeds = list(map(int, cfg["odor_encoder"]["evaluation_seeds"]))
|
|
runs: dict[tuple[str, str, int], dict[str, Any]] = {}; summary = {}
|
|
for mode in cfg["odor_encoder"]["modes"]:
|
|
for oid in panel.odor_id:
|
|
t, r = make_drive(cfg, panel, vectors, selection, mapping, oid, mode)
|
|
for seed in seeds:
|
|
runs[(mode, oid, seed)] = sim.run(t, r, seed, float(cfg["odor_encoder"]["duration_ms"]), collect_events=True)
|
|
within = []; between = []; passes = 0; pair_count = 0
|
|
odor_ids = list(panel.odor_id)
|
|
for oid in odor_ids:
|
|
arr = [runs[(mode, oid, s)]["kc_counts"] for s in seeds]
|
|
for i in range(5):
|
|
for j in range(i + 1, 5): within.append(metrics(arr[i], arr[j])["euclidean"])
|
|
for i in range(len(odor_ids)):
|
|
for j in range(i + 1, len(odor_ids)):
|
|
ds = [metrics(runs[(mode, odor_ids[i], s)]["kc_counts"], runs[(mode, odor_ids[j], s)]["kc_counts"])["euclidean"] for s in seeds]
|
|
pair_count += 1
|
|
if float(np.mean(ds)) > float(np.mean(within)): passes += 1
|
|
between.extend(ds)
|
|
frac = [np.count_nonzero(run["kc_counts"]) / len(sim.kc) for (m, _, _), run in runs.items() if m == mode]
|
|
summary[mode] = {"KC_active_fraction_mean": float(np.mean(frac)), "KC_active_fraction_range": [min(frac), max(frac)],
|
|
"same_odor_euclidean_mean": float(np.mean(within)), "different_odor_euclidean_mean": float(np.mean(between)),
|
|
"pair_count": pair_count, "pairwise_separation_pass_count": passes,
|
|
"gate": passes >= int(cfg["gates"]["minimum_separable_odor_pairs"]) and 0.03 <= float(np.mean(frac)) <= 0.09}
|
|
passed = all(x["gate"] for x in summary.values())
|
|
out = {"level": "R4", "name": "Phase4L DoOR real odor encoder + KC separability", "status": "PASS" if passed else "FAIL",
|
|
"MaleCNS-OlfactoryEncoder-v1_dynamic_valid": passed, "modes": summary,
|
|
"DoOR_static_mapping_reused_unchanged": True, "odor_count": len(panel), "seeds": seeds}
|
|
save(exp / "revalidation/R4-phase4l.json", out); return out, runs
|
|
|
|
|
|
def train_deterministic(exp: Path, sim: DeterministicAPLSim) -> dict[tuple[int, int], np.ndarray]:
|
|
cfg, panel, vectors, selection, mapping = door_inputs(); mcfg = yaml.safe_load((ROOT / "models/male-cns-lif/kcmbon-efficacy-v1/phase4m.yaml").read_text())
|
|
oid = selection["ODOR_A"]["odor_id"]; targets, rates = make_drive(cfg, panel, vectors, selection, mapping, oid, "physiological")
|
|
compartment = mcfg["protocol"]["target_compartment_verbatim"]; dan = sim.dan_targets(mcfg["protocol"]["matched_dan_type"], compartment)
|
|
comp = sim.identity.compartment.astype(str).to_numpy() == compartment
|
|
kc_lookup = np.full(sim.n, -1, np.int32); kc_lookup[sim.kc] = np.arange(len(sim.kc), dtype=np.int32)
|
|
source_local = kc_lookup[sim.identity.source_neuron_idx.to_numpy(np.int32)]
|
|
if np.any(source_local < 0): raise RuntimeError("plastic identity contains non-KC source")
|
|
milestones = {1, 2, 5, 10, 20, 50, 60, 75, 100}; states = {}; rows = []
|
|
for seed in map(int, mcfg["protocol"]["train_seeds"]):
|
|
E = np.ones(len(sim.edge_ids), np.float32)
|
|
for pairing in range(1, 101):
|
|
run = sim.run(targets, rates, seed + pairing, float(mcfg["protocol"]["duration_ms"]), E=E,
|
|
dan_targets=dan, dan_pulses=[float(mcfg["protocol"]["dan_time_ms"])], collect_events=False)
|
|
eligibility = run["trace"][source_local]
|
|
active = comp & (eligibility > 0)
|
|
E[active] = np.maximum(np.float32(0.75), np.float32(E[active] * np.float32(1.0 - np.float32(0.01) * eligibility[active])))
|
|
if pairing in milestones:
|
|
states[(seed, pairing)] = E.copy()
|
|
rows.append({"seed": seed, "pairing": pairing, "updated_edges": int(np.count_nonzero(E < 1)),
|
|
"mean_depression": float(np.mean(1 - E)), "max_depression": float(np.max(1 - E)),
|
|
"saturated_edges": int(np.count_nonzero(E <= 0.75)), "trace_hash": ahash(run["trace"])})
|
|
if pairing % 20 == 0: print(json.dumps({"R5_training_seed": seed, "pairing": pairing}), flush=True)
|
|
pd.DataFrame(rows).to_parquet(exp / "revalidation/R5-deterministic-training-trajectory.parquet", index=False)
|
|
np.savez_compressed(exp / "revalidation/R5-efficacy-states.npz", **{f"seed_{s}_pairing_{p}": x for (s, p), x in states.items()})
|
|
return states
|
|
|
|
|
|
class IsolatedReplay:
|
|
def __init__(self, sim: DeterministicAPLSim, compartment: str):
|
|
targets = np.sort(sim.identity.loc[sim.identity.compartment.astype(str) == compartment, "target_neuron_idx"].unique().astype(np.int32))
|
|
mask = np.isin(sim.sources, sim.kc) & np.isin(sim.rt.indices_np, targets)
|
|
self.ids = np.flatnonzero(mask).astype(np.int64); self.src = sim.sources[self.ids]; self.target = sim.rt.indices_np[self.ids]
|
|
self.targets = targets; pos = {int(x): i for i, x in enumerate(targets)}
|
|
self.local = np.asarray([pos[int(x)] for x in self.target], np.int32)
|
|
self.base = sim.base[self.ids].cpu().numpy().astype(np.float32); self.delay = sim.rt.delay_steps; self.qlen = sim.rt.qlen
|
|
self.dt = sim.dt; self.gd = sim.rt.gd; self.alpha = sim.rt.alpha; self.rest = sim.rt.rest; self.reset = sim.rt.resetv
|
|
self.threshold = sim.rt.threshold; self.refsteps = sim.rt.ref_steps
|
|
can = {int(e): i for i, e in enumerate(sim.edge_ids)}
|
|
self.candidate = np.asarray([can.get(int(e), -1) for e in self.ids], np.int32)
|
|
self.by_source = {int(s): np.flatnonzero(self.src == s) for s in np.unique(self.src)}
|
|
|
|
def run(self, steps: np.ndarray, nodes: np.ndarray, E: np.ndarray | None, duration_ms: float) -> dict[str, Any]:
|
|
w = self.base.copy(); valid = self.candidate >= 0
|
|
if E is not None: w[valid] = np.float32(w[valid] * E[self.candidate[valid]])
|
|
grouped = {int(s): np.sort(nodes[steps == s]) for s in np.unique(steps)}
|
|
n = len(self.targets); q = np.zeros((self.qlen, n), np.float32); g = np.zeros(n, np.float32)
|
|
v = np.full(n, self.rest, np.float32); ref = np.zeros(n, np.int32); counts = np.zeros(n, np.int32); first = np.full(n, -1, np.int32)
|
|
delivered = np.float64(0)
|
|
for step in range(int(round(duration_ms / self.dt))):
|
|
slot = step % self.qlen; g = np.float32(np.float32(g * self.gd) + q[slot]); delivered += q[slot].sum(dtype=np.float64); q[slot].fill(0)
|
|
ref = np.maximum(ref, 1) - 1; avail = ref == 0
|
|
v[avail] = np.float32(v[avail] + np.float32(self.alpha * np.float32(g[avail] - np.float32(v[avail] - self.rest))))
|
|
v[~avail] = self.reset; fired = avail & (v >= self.threshold); counts += fired; first[(first < 0) & fired] = step
|
|
v[fired] = self.reset; ref[fired] = self.refsteps
|
|
for source in grouped.get(step, np.empty(0, np.int32)):
|
|
rel = self.by_source.get(int(source))
|
|
if rel is not None:
|
|
dest = q[(step + self.delay) % self.qlen]
|
|
for j in rel: dest[self.local[j]] = np.float32(dest[self.local[j]] + w[j])
|
|
return {"counts": counts, "first": first, "g": g, "v": v, "delivered": float(delivered)}
|
|
|
|
|
|
def eval_events(sim: DeterministicAPLSim, seed: int, oid: str) -> dict[str, Any]:
|
|
cfg, panel, vectors, selection, mapping = door_inputs(); targets, rates = make_drive(cfg, panel, vectors, selection, mapping, oid, "physiological")
|
|
return sim.run(targets, rates, seed, 120.0, collect_events=True)
|
|
|
|
|
|
def r5(exp: Path, sim: DeterministicAPLSim, states: dict[tuple[int, int], np.ndarray]) -> tuple[dict[str, Any], dict[str, Any]]:
|
|
cfg, panel, vectors, selection, mapping = door_inputs(); mcfg = yaml.safe_load((ROOT / "models/male-cns-lif/kcmbon-efficacy-v1/phase4m.yaml").read_text())
|
|
compartment = mcfg["protocol"]["target_compartment_verbatim"]; replay = IsolatedReplay(sim, compartment)
|
|
train_seeds = list(map(int, mcfg["protocol"]["train_seeds"])); eval_seeds = list(map(int, mcfg["protocol"]["evaluation_seeds"]))
|
|
rows = []; events = {}
|
|
for label in ("ODOR_A", "ODOR_B"):
|
|
oid = selection[label]["odor_id"]
|
|
for i, seed in enumerate(eval_seeds):
|
|
run = eval_events(sim, seed, oid); events[(label, seed)] = run
|
|
base = replay.run(run["kc_event_step"], run["kc_event_node"], None, 120.0)
|
|
learned = replay.run(run["kc_event_step"], run["kc_event_node"], states[(train_seeds[i], 100)], 120.0)
|
|
rows.append({"odor": label, "seed": seed, "train_seed": train_seeds[i],
|
|
"spike_L1": int(np.abs(base["counts"] - learned["counts"]).sum()),
|
|
"changed_MBON": int(np.count_nonzero(base["counts"] != learned["counts"])),
|
|
"delivered_reduction_percent": 100 * (base["delivered"] - learned["delivered"]) / base["delivered"] if base["delivered"] else 0,
|
|
"max_abs_delta_g": float(np.max(np.abs(base["g"] - learned["g"]))),
|
|
"max_abs_delta_v": float(np.max(np.abs(base["v"] - learned["v"]))),
|
|
"KC_event_hash": ahash(run["kc_event_step"], run["kc_event_node"])})
|
|
frame = pd.DataFrame(rows); frame.to_parquet(exp / "revalidation/R5-isolated-replay.parquet", index=False)
|
|
a = frame[frame.odor == "ODOR_A"]; b = frame[frame.odor == "ODOR_B"]
|
|
functional = int((a.spike_L1 > 0).sum()) >= 3
|
|
specific = float(a.spike_L1.mean()) > float(b.spike_L1.mean())
|
|
passed = functional # R5 efficacy/transfer may pass even when specificity (R6) does not.
|
|
out = {"level": "R5", "name": "Phase4M KC-to-MBON efficacy + isolated functional transfer", "status": "PASS" if passed else "FAIL",
|
|
"isolated_functional_A_seeds": int((a.spike_L1 > 0).sum()), "A_spike_L1_mean": float(a.spike_L1.mean()),
|
|
"B_spike_L1_mean": float(b.spike_L1.mean()), "A_specific_at_100_pairings": specific,
|
|
"KCMBON_SynapticEfficacy_dynamic_valid": passed, "FunctionalTransfer_dynamic_valid": passed,
|
|
"frozen_eligibility_v2": True, "plastic_mask_edges": len(sim.edge_ids), "rows": rows}
|
|
save(exp / "revalidation/R5-phase4m.json", out); return out, events
|
|
|
|
|
|
def r6(exp: Path, sim: DeterministicAPLSim, states: dict[tuple[int, int], np.ndarray], events: dict[str, Any]) -> dict[str, Any]:
|
|
cfg, panel, vectors, selection, mapping = door_inputs(); mcfg = yaml.safe_load((ROOT / "models/male-cns-lif/kcmbon-efficacy-v1/phase4m.yaml").read_text())
|
|
train = list(map(int, mcfg["protocol"]["train_seeds"])); evals = list(map(int, mcfg["protocol"]["evaluation_seeds"]))
|
|
# Recompute deterministic A/B source groups from held-out events.
|
|
active = {}
|
|
for label in ("ODOR_A", "ODOR_B"):
|
|
nodes = np.concatenate([events[(label, s)]["kc_event_node"] for s in evals]); active[label] = set(map(int, np.unique(nodes)))
|
|
source = sim.identity.source_neuron_idx.to_numpy(np.int32)
|
|
groups = np.full(len(source), "inactive", object)
|
|
A = np.isin(source, list(active["ODOR_A"])); B = np.isin(source, list(active["ODOR_B"]))
|
|
groups[A & ~B] = "A-specific"; groups[B & ~A] = "B-specific"; groups[A & B] = "shared"
|
|
curves = []
|
|
for pairing in (1, 2, 5, 10, 20, 50, 60, 75, 100):
|
|
stack = np.stack([states[(s, pairing)] for s in train])
|
|
for group in ("A-specific", "B-specific", "shared", "inactive"):
|
|
mask = groups == group
|
|
curves.append({"pairing": pairing, "group": group, "edges": int(mask.sum()),
|
|
"mean_delta_E": float(np.mean(stack[:, mask] - 1)) if mask.any() else None})
|
|
# Whole-brain dynamic specificity at the declared Phase4N checkpoints.
|
|
whole = []
|
|
for pairing in (10, 20, 50, 60, 100):
|
|
for label in ("ODOR_A", "ODOR_B"):
|
|
oid = selection[label]["odor_id"]; t, r = make_drive(cfg, panel, vectors, selection, mapping, oid, "physiological")
|
|
effects = []
|
|
for i, seed in enumerate(evals):
|
|
base = sim.run(t, r, seed, 120.0)
|
|
post = sim.run(t, r, seed, 120.0, E=states[(train[i], pairing)])
|
|
effects.append(int(np.abs(base["mbon_counts"] - post["mbon_counts"]).sum()))
|
|
whole.append({"pairing": pairing, "odor": label, "MBON_spike_L1_mean": float(np.mean(effects)),
|
|
"functional_seeds": int(np.count_nonzero(effects))})
|
|
print(json.dumps({"R6_pairing": pairing}), flush=True)
|
|
wf = pd.DataFrame(whole); cf = pd.DataFrame(curves)
|
|
wf.to_parquet(exp / "revalidation/R6-wholebrain-specificity.parquet", index=False)
|
|
cf.to_parquet(exp / "revalidation/R6-edge-group-curves.parquet", index=False)
|
|
at50 = wf[wf.pairing == 50].set_index("odor").MBON_spike_L1_mean
|
|
specificity = bool(at50.get("ODOR_A", 0) > at50.get("ODOR_B", 0))
|
|
shared_contamination = bool(abs(cf[(cf.pairing == 50) & (cf.group == "shared")].mean_delta_E.iloc[0]) >
|
|
abs(cf[(cf.pairing == 50) & (cf.group == "inactive")].mean_delta_E.iloc[0]))
|
|
# Historical Phase4N result was a partial audit/no candidate. Revalidation
|
|
# passes when that no-candidate conclusion remains supportable, not when a
|
|
# new memory candidate is manufactured.
|
|
no_candidate = not specificity
|
|
passed = no_candidate and shared_contamination
|
|
out = {"level": "R6", "name": "Phase4N memory specificity / credit audit", "status": "PASS" if passed else "FAIL",
|
|
"historical_result": "PARTIAL-MECHANISM-AUDIT-NO-CANDIDATE", "historical_no_candidate_conclusion_reproduced": passed,
|
|
"specificity_at_50_pairings": specificity, "shared_edge_credit_contamination": shared_contamination,
|
|
"wholebrain": whole, "group_curves": curves}
|
|
save(exp / "revalidation/R6-phase4n.json", out); return out
|
|
|
|
|
|
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()); exp = Path(cfg["paths"]["experiment_root"]); model = Path(cfg["paths"]["model_root"])
|
|
levels = []; r1, r2 = r1_r2(exp); levels += [r1, r2]
|
|
if r1["status"] != "PASS" or r2["status"] != "PASS":
|
|
save(exp / "revalidation/summary.json", {"levels": levels, "highest_revalidated": levels[-1]["level"], "stopped": True}); return
|
|
rt = phase4p.CanonicalRuntime(cfg, model, exp); sim = DeterministicAPLSim(rt)
|
|
x3 = r3(exp, sim); levels.append(x3)
|
|
if x3["status"] != "PASS": save(exp / "revalidation/summary.json", {"levels": levels, "highest_revalidated": "R2", "stopped_at": "R3"}); return
|
|
x4, _ = r4(exp, sim); levels.append(x4)
|
|
if x4["status"] != "PASS": save(exp / "revalidation/summary.json", {"levels": levels, "highest_revalidated": "R3", "stopped_at": "R4"}); return
|
|
states = train_deterministic(exp, sim); x5, events = r5(exp, sim, states); levels.append(x5)
|
|
if x5["status"] != "PASS": save(exp / "revalidation/summary.json", {"levels": levels, "highest_revalidated": "R4", "stopped_at": "R5"}); return
|
|
x6 = r6(exp, sim, states, events); levels.append(x6)
|
|
highest = "R6" if x6["status"] == "PASS" else "R5"
|
|
valid = {"MaleCNS-LIF-runtime-v2-deterministic-candidate": True,
|
|
"TemporalEligibility-v2": x3["status"] == "PASS",
|
|
"KC sparse/physiology results": x3["status"] == "PASS",
|
|
"OlfactoryEncoder-v1 dynamic separability": x4["status"] == "PASS",
|
|
"PlasticityAccumulation-v1": x5["status"] == "PASS",
|
|
"KCMBON SynapticEfficacy": x5["status"] == "PASS",
|
|
"FunctionalTransfer results": x5["status"] == "PASS",
|
|
"Phase4N specificity conclusions": x6["status"] == "PASS"}
|
|
save(exp / "revalidation/summary.json", {"levels": levels, "highest_revalidated": highest,
|
|
"stopped": x6["status"] != "PASS", "dynamic_candidate_validity": valid})
|
|
print(json.dumps({"highest_revalidated": highest, "statuses": {x["level"]: x["status"] for x in levels}}, indent=2), flush=True)
|
|
|
|
|
|
if __name__ == "__main__": main()
|