From 50d7c8329f92bc3db446cd965f7c324a34231dc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ahmet=20Bar=C4=B1=C5=9F=20G=C3=BCnayd=C4=B1n?= Date: Fri, 31 Jul 2026 20:22:49 +0700 Subject: [PATCH] feat(vnc): resolve motor output per muscle, and test the rhythm question MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three things: an instrumentation change, a builder fix that is inert until the blob is rebuilt, and the analysis that answers whether this connectome can walk. MUSCLE-RESOLVED READOUT. main.ts averaged each leg's ~60 motor neurons into one scalar, which cancels every flexor against its extensor by construction. vnc.meta.json already partitions the 369 leg motor neurons into 6 legs x up to 20 named muscle targets, and motor_by_target was loaded and never read. The readout now resolves per muscle and computes antagonist ratios for the three joints with a true pair in all six legs — tibia, femur/trochanter, coxa_twist — pooled per side and normalised, since the pools are lopsided (14-15 flexor MNs against 2 extensor at the tibia). tergotrochanter is excluded: it is the escape-jump muscle and does not belong in a walking drive. This is instrumentation only; driveLegs is untouched and nothing consumes the numbers yet. Measured live under a DNa01 stim: per-muscle values vary over time and differ between legs, so the signal is alive. GLUTAMATE SIGN. build_vnc.py signed glutamate -1 for every neuron. That is right centrally, where fly glutamate is inhibitory via GluCl, and wrong at the neuromuscular junction, where it is the excitatory transmitter. Motor neurons are now gated to +1. build_csr.py is deliberately unchanged — the brain has no NMJ. Scope, measured rather than assumed: of 709 motor neurons in the shipped bin, 580 (81.8%) carry sign -1. But that sign is read in exactly one place in src/ — a boot log line — and the weights motor neurons RECEIVE, which is what drives the readout, come from their presynaptic partners and are untouched. Motor neurons emit 0.126% of total synaptic mass. The maximal counterfactual moves the structural result 0.67% and the rhythm sweep not at all. This is a labelling and biophysics correction, not a dynamics fix. INERT UNTIL REBUILT: data/manc is not in the checkout, so vnc.bin was not regenerated. Re-run download_manc.sh, build_vnc.py, write_manifest.py and upload_to_r2.sh to land it, and re-verify walking afterwards — the drive constants in main.ts were tuned against the current blob. tools/vnc_rhythm.py answers the question the readout was built to ask, and the answer is no. Findings, all reproducible from that script against the shipped bin: - The network DOES oscillate — 5-12 Hz, spectral prominence 58x to 1.8e5x over background. An earlier claim that it never oscillates was wrong. - But it is global synchrony, not gait: all six legs peak at one frequency in phase, Kuramoto R_sync 0.475-0.910 against R_tripod 0.049-0.331, and across every probed setting the best-fitting 3/3 split of the legs was never the true tripod. - Tripod alternation is not reachable by tuning, and this is a proof rather than a failed search: the excitatory fraction is 0.517-0.632 in all 36 leg-pair couplings, so nulling across-tripod coupling needs an inhibition multiplier above 1.716 while keeping within-tripod positive needs it below 1.326. The interval is empty. - 600-setting sweep (100x weight, 20x drive, 16x E/I): 0/600 whole-leg rhythm, 0/600 muscle-resolved under the strict criterion. Resolving antagonists does not reveal a hidden rhythm — that hypothesis is dead. - The tripod signal in the wiring is fully explained by segment distance and ipsi/contra; residual after removing anatomy is identically zero. The connectome is organised by proximity, not by gait. Pool-level averaging is the wrong lens, and any future attempt needs per-neuron structure. Co-Authored-By: Claude Opus 5 --- src/main.ts | 75 ++++++ tools/build_vnc.py | 17 +- tools/vnc_rhythm.py | 557 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 647 insertions(+), 2 deletions(-) create mode 100644 tools/vnc_rhythm.py diff --git a/src/main.ts b/src/main.ts index 2137824..783a97d 100644 --- a/src/main.ts +++ b/src/main.ts @@ -367,11 +367,28 @@ async function main() { */ interface MancMotor { legs: Record; // T1_left … T3_right + muscles: Record>; // leg → muscle → mean rate + antagonists: Record>; // leg → joint → ext-vs-flex in [-1,1] wing: number; // mean wing motor rate neck: number; // mean neck motor rate abdomen: number; // mean abdomen motor rate } + const LEG_READOUT_ORDER = ["T1_left", "T1_right", "T2_left", "T2_right", "T3_left", "T3_right"]; + // Antagonist muscle pairs, verified present in all six legs of + // vnc.meta.json's motor groups. `legs` above averages a leg's ~60 MNs + // into one number, which cancels every flexor against its extensor; + // these pools keep the opposition. + // + // tergotr (tergotrochanter) is deliberately absent from every pool: + // it is the escape-jump muscle driven by the giant fiber, not a + // walking muscle, so it must never enter a walking-related drive. + const ANTAGONISTS: Record = { + tibia: { ext: ["ti_extensor"], flex: ["ti_flexor", "acc_ti_flexor"] }, + femur: { ext: ["tr_extensor", "sternotrochanter"], flex: ["tr_flexor", "acc_tr_flexor"] }, + coxa_twist: { ext: ["sternal_anterior_rotator"], flex: ["sternal_posterior_rotator"] }, + }; + async function runRealSpine(brainRate: Float32Array): Promise { if (!vncSim || !vncMeta) return null; const N = vncSim.brain.header.numNeurons; @@ -393,10 +410,39 @@ async function main() { const vncRate = await vncSim.captureRollingRate(20); // Aggregate motor neuron rates per leg+side. const legs: Record = {}; + const muscles: Record> = {}; + const antagonists: Record> = {}; for (const [legKey, group] of Object.entries(vncMeta.motor)) { let sum = 0; for (const i of group.all) sum += vncRate[i]; legs[legKey] = group.all.length ? sum / group.all.length : 0; + // Pool mean, not mean-of-subclass-means: the pools are very + // unbalanced (T1 has 14–15 tibia flexor MNs against 2 extensors), + // so each side is divided by its own neuron count before the + // ratio, otherwise the difference just reports pool size. + const poolMean = (subs: string[]): number => { + let s = 0, n = 0; + for (const sub of subs) { + const idxs = group[sub]; + if (!idxs) continue; + for (const i of idxs) s += vncRate[i]; + n += idxs.length; + } + return n ? s / n : 0; + }; + const perMuscle: Record = {}; + for (const sub of Object.keys(group)) { + if (sub === "all") continue; + perMuscle[sub] = poolMean([sub]); + } + muscles[legKey] = perMuscle; + const joints: Record = {}; + for (const [joint, pair] of Object.entries(ANTAGONISTS)) { + const aExt = poolMean(pair.ext); + const aFlex = poolMean(pair.flex); + joints[joint] = (aExt - aFlex) / (aExt + aFlex + 1e-9); + } + antagonists[legKey] = joints; } // Subclass aggregates (wm = wing, nm = neck, ad = abdomen). const aggSub = (key: string): number => { @@ -408,6 +454,8 @@ async function main() { }; return { legs, + muscles, + antagonists, wing: aggSub("wm"), neck: aggSub("nm"), abdomen: aggSub("ad"), @@ -417,6 +465,13 @@ async function main() { let driveFwd = 0, driveTurn = 0; // smoothed let lastSpineMode = "synthetic"; // tracked for the diagnostic readout let lastMancActivity = { legs: 0, wing: 0, neck: 0, abdomen: 0 }; + // Muscle-resolved readout. Instrumentation only — nothing downstream + // consumes it yet; driveLegs still runs off fwd/turn. + const mancMuscleProbe: { + muscles: Record>; + antagonists: Record>; + } = { muscles: {}, antagonists: {} }; + (window as unknown as { __mancMuscles: typeof mancMuscleProbe }).__mancMuscles = mancMuscleProbe; // Build VNC stand-in context once; reused by every drive update. const vncCtx: MotorContext = { famousDns, @@ -452,6 +507,8 @@ async function main() { mancTotal = meanL + meanR; mancAsym = meanR - meanL; lastMancActivity = { legs: mancTotal, wing: realMotor.wing, neck: realMotor.neck, abdomen: realMotor.abdomen }; + mancMuscleProbe.muscles = realMotor.muscles; + mancMuscleProbe.antagonists = realMotor.antagonists; } const visualActive = visual && Number.isFinite(visual.angle) && visual.area > 0; if (visualActive) { @@ -742,6 +799,7 @@ async function main() { const n = Math.max(0, Math.min(w, Math.round(r * w * 3))); return "█".repeat(n) + "·".repeat(w - n); }; + const signed = (v: number) => (v >= 0 ? "+" : "") + v.toFixed(2); const lines = [ `fwd ${driveFwd.toFixed(2)} turn ${driveTurn.toFixed(2)} speed ${sp.toFixed(2)} cm/s [${lastSpineMode}]`, `synth: fwd ${bar(vnc.fwdRate)} bwd ${bar(vnc.bwdRate)} esc ${bar(vnc.escape)}`, @@ -749,6 +807,23 @@ async function main() { if (vncMeta) { const m = lastMancActivity; lines.push(`MANC : leg ${bar(m.legs)} wing ${bar(m.wing)} abd ${bar(m.abdomen)}`); + const ant = mancMuscleProbe.antagonists; + // A ratio of 0 is ambiguous — balanced pools and two silent pools + // both give 0 — so silence prints as a dash instead of "+0.00". + const silent = (leg: string) => { + const mus = mancMuscleProbe.muscles[leg] ?? {}; + return [...ANTAGONISTS.tibia.ext, ...ANTAGONISTS.tibia.flex] + .every((sub) => (mus[sub] ?? 0) <= 0); + }; + const cells = LEG_READOUT_ORDER + .filter((k) => k in ant) + .map((k) => { + const lbl = k.replace("_left", "L").replace("_right", "R"); + return `${lbl} ${silent(k) ? " — " : signed(ant[k].tibia)}`; + }); + if (cells.length) { + lines.push(`MANC ti ext−flex: ${cells.join(" ")}`); + } } driveReadout.innerHTML = lines.join("
"); }, 50); diff --git a/tools/build_vnc.py b/tools/build_vnc.py index 8cc443b..aca3fc3 100644 --- a/tools/build_vnc.py +++ b/tools/build_vnc.py @@ -80,6 +80,13 @@ "unclear": 0, "unknown": 0, } + +# Motor neurons are the exception to the central glutamate rule: at the +# neuromuscular junction fly muscle expresses ionotropic GluRs, not GluCl, +# so a glutamatergic motor neuron is EXCITATORY onto its muscle. The sign +# here also governs the motor neuron's (few) intra-VNC output edges. +MOTOR_NT_SIGN = {**NT_SIGN, "glutamate": +1} + NT_CONF_MIN = 0.5 @@ -248,10 +255,16 @@ def main() -> None: print("signing weights by presynaptic neurotransmitter …") pre_nt_by_idx = df["predictedNt"].fillna("unknown").astype(str).str.lower().to_numpy() pre_conf = pd.to_numeric(df["predictedNtProb"], errors="coerce").fillna(0.0).to_numpy() + class_arr = df.apply(classify, axis=1).to_numpy() sign_arr = np.zeros(N, dtype=np.float32) for i, nt in enumerate(pre_nt_by_idx): if pre_conf[i] >= NT_CONF_MIN: - sign_arr[i] = NT_SIGN.get(nt, 0) + table = MOTOR_NT_SIGN if class_arr[i] == CLASS_MOTOR else NT_SIGN + sign_arr[i] = table.get(nt, 0) + + n_nmj = int(((class_arr == CLASS_MOTOR) & (pre_nt_by_idx == "glutamate") + & (pre_conf >= NT_CONF_MIN)).sum()) + print(f" {n_nmj:,} glutamatergic motor neurons signed +1 (NMJ)") pre_idx = edges["pre_idx"].to_numpy(dtype=np.uint32) post_idx = edges["post_idx"].to_numpy(dtype=np.uint32) @@ -283,7 +296,7 @@ def main() -> None: soma = row.get("position") x, y, z = parse_position(soma) pos[i] = (x, y, z) - cell_class[i] = classify(row) + cell_class[i] = class_arr[i] leg_seg[i] = leg_segment_packed(row) # predictedNtProb gives float confidence; default 0 nt_conf[i] = float(pre_conf[i]) diff --git a/tools/vnc_rhythm.py b/tools/vnc_rhythm.py new file mode 100644 index 0000000..33c9274 --- /dev/null +++ b/tools/vnc_rhythm.py @@ -0,0 +1,557 @@ +#!/usr/bin/env python3 +""" +vnc_rhythm.py — offline test of whether a walking rhythm exists in the +shipped MANC connectome (public/vnc.bin + public/vnc.meta.json). + +Read-only. Touches no repo file. Run with: + uv run --with numpy --with scipy python3 vnc_rhythm.py struct + uv run --with numpy --with scipy python3 vnc_rhythm.py sweep + +Dynamics replicate src/shaders/lif.wgsl exactly, with the constants the +VNC sim is actually created with (src/main.ts:594 -> DEFAULT_PARAMS in +src/sim.ts:31-40): + dt = 1 ms, tau_m = 20 ms -> alpha = exp(-1/20) + v_thresh = -45, v_reset = v_rest = -52 + refractory 2.2 ms -> 2 steps + ext_gain = 2.0, w_syn = 0.005 + alpha synapse A_SYN = 0.81873 (lif.wgsl:52) +""" +from __future__ import annotations + +import json +import struct +import sys +import time +from pathlib import Path + +import numpy as np + +ROOT = Path("/Users/ahmetbarisgunaydin/dev/webgpu-fly") +BIN = ROOT / "public" / "vnc.bin" +META = ROOT / "public" / "vnc.meta.json" + +CLASS_UNKNOWN, CLASS_DN_INPUT, CLASS_MOTOR, CLASS_SENSORY, CLASS_INTRINSIC = 0, 1, 2, 3, 4 + +# --- LIF constants (see docstring) --- +ALPHA = float(np.exp(-1.0 / 20.0)) +V_TH, V_RESET, V_REST = -45.0, -52.0, -52.0 +REFRAC_STEPS = 2 +EXT_GAIN = 2.0 +W_SYN = 0.005 +A_SYN = 0.81873 + +# src/vnc.ts:51-54. LEG_KEYS order [T1L,T1R,T2L,T2R,T3L,T3R] = idx 0..5; +# tripod A = idx 0,3,4 ; tripod B = idx 1,2,5. +LEG_KEYS = ["T1_left", "T1_right", "T2_left", "T2_right", "T3_left", "T3_right"] +TRIPOD_A = ["T1_left", "T2_right", "T3_left"] +TRIPOD_B = ["T1_right", "T2_left", "T3_right"] + +# src/main.ts:386-390 +ANTAGONISTS = { + "tibia": (["ti_extensor"], ["ti_flexor", "acc_ti_flexor"]), + "femur": (["tr_extensor", "sternotrochanter"], ["tr_flexor", "acc_tr_flexor"]), + "coxa_twist": (["sternal_anterior_rotator"], ["sternal_posterior_rotator"]), +} + + +# ---------------------------------------------------------------- loading +def load_vnc(): + raw = np.memmap(BIN, dtype=np.uint8, mode="r") + magic = bytes(raw[:8]) + assert magic == b"WGFLYVNC", magic + version, N, E, flags = struct.unpack(" post in + targets[a] (default targets = pools). Returns (raw, per-pair-density).""" + row_ptr, col_idx = d["row_ptr"], d["col_idx"] + if targets is None: + targets = pools + n = len(LEG_KEYS) + raw = np.zeros((n, n)) + dens = np.zeros((n, n)) + memb = {} + for j, leg in enumerate(LEG_KEYS): + m = np.zeros(d["N"], dtype=bool) + m[pools[leg]] = True + memb[j] = m + for a, la in enumerate(LEG_KEYS): + tgt = targets[la] + # concatenate all incoming edges of the target set + segs_c, segs_w = [], [] + for t in tgt: + s, e = row_ptr[t], row_ptr[t + 1] + segs_c.append(col_idx[s:e]) + segs_w.append(weight[s:e]) + c = np.concatenate(segs_c) if segs_c else np.zeros(0, np.int64) + w = np.concatenate(segs_w) if segs_w else np.zeros(0, np.float32) + for b in range(len(LEG_KEYS)): + sel = memb[b][c] + raw[a, b] = w[sel].sum() + dens[a, b] = raw[a, b] / (len(tgt) * len(pools[LEG_KEYS[b]])) + return raw, dens + + +def tripod_split(M, groupA=None): + """Return (within-tripod off-diagonal mean, across-tripod mean).""" + if groupA is None: + groupA = TRIPOD_A + tri = {k: ("A" if k in groupA else "B") for k in LEG_KEYS} + within, across = [], [] + for a, la in enumerate(LEG_KEYS): + for b, lb in enumerate(LEG_KEYS): + if a == b: + continue + (within if tri[la] == tri[lb] else across).append(M[a, b]) + return float(np.mean(within)), float(np.mean(across)), within, across + + +def tripod_permutation_control(M): + """All 10 distinct 3/3 splits of the six legs. Where does the true + tripod rank on (across - within)?""" + from itertools import combinations + seen, out = set(), [] + for c in combinations(LEG_KEYS, 3): + key = frozenset(c) + comp = frozenset(set(LEG_KEYS) - key) + if comp in seen: + continue + seen.add(key) + wi, ac, _, _ = tripod_split(M, groupA=list(c)) + out.append((ac - wi, tuple(sorted(c)))) + out.sort(reverse=True) + return out + + +def fmt_matrix(M, scale=1.0): + hdr = " " + "".join(f"{k:>12}" for k in LEG_KEYS) + lines = [hdr] + for a, la in enumerate(LEG_KEYS): + lines.append(f"{la:>12}" + "".join(f"{M[a,b]*scale:12.4g}" for b in range(6))) + return "\n".join(lines) + + +def stage_struct(d): + cc = d["cell_class"] + print(f"vnc.bin N={d['N']:,} E={d['E']:,} flags={d['flags']} version={d['version']}") + print("cell_class counts:", {int(k): int(v) for k, v in + zip(*np.unique(cc, return_counts=True))}) + w0 = d["weight"] + wf, nflip = glut_fix_weights(w0, d["col_idx"], cc) + print(f"glutamate-fix counterfactual flips {nflip:,} edges " + f"({100*nflip/d['E']:.3f}% of E), |w| moved = {2*np.abs(w0[w0<0][:0]).sum() if False else np.abs(w0[(cc[d['col_idx']]==CLASS_MOTOR)&(w0<0)]).sum():,.0f} " + f"of {np.abs(w0).sum():,.0f} total |w| ({100*np.abs(w0[(cc[d['col_idx']]==CLASS_MOTOR)&(w0<0)]).sum()/np.abs(w0).sum():.3f}%)") + + for include_motor in (False, True): + tag = "premotor INCLUDES motor neurons" if include_motor else "premotor EXCLUDES motor neurons" + pools = premotor_pools(d, include_motor) + sizes = {k: len(v) for k, v in pools.items()} + allp = np.concatenate([pools[k] for k in LEG_KEYS]) + uniq = len(np.unique(allp)) + print(f"\n===== {tag} =====") + print("pool sizes:", sizes, f" union={uniq:,} (sum={len(allp):,}, overlap factor {len(allp)/uniq:.2f}x)") + + for label, wts in (("SHIPPED", w0), ("GLUT-FIXED", wf)): + raw, dens = signed_matrix(d, pools, wts) + npos = int((raw > 0).sum()) + wi, ac, wl, al = tripod_split(dens) + wi_r, ac_r, _, _ = tripod_split(raw) + print(f"\n--- {label} premotor<-premotor signed net weight (raw sums) ---") + print(fmt_matrix(raw)) + print(f" positive entries: {npos}/36 negative: {36-npos}/36 " + f"min={raw.min():.4g} max={raw.max():.4g}") + print(f" per-pair density (raw / |Pa|*|Pb|), within-tripod mean={wi:.6g} " + f"across-tripod mean={ac:.6g} across-minus-within={ac-wi:+.6g}") + print(f" raw-sum within={wi_r:.6g} across={ac_r:.6g} across-minus-within={ac_r-wi_r:+.6g}") + + # direct: premotor -> that leg's motor pool + targets = {k: np.array(d["meta"]["motor"][k]["all"], dtype=np.int64) for k in LEG_KEYS} + rawD, densD = signed_matrix(d, pools, wts, targets=targets) + nposD = int((rawD > 0).sum()) + wiD, acD, _, _ = tripod_split(densD) + print(f"--- {label} motor<-premotor (leg a's MOTOR pool <- leg b's premotor) ---") + print(fmt_matrix(rawD)) + print(f" positive entries: {nposD}/36 min={rawD.min():.4g} max={rawD.max():.4g}") + print(f" density within-tripod={wiD:.6g} across-tripod={acD:.6g} " + f"across-minus-within={acD-wiD:+.6g}") + + # motor <- motor directly + rawM, densM = signed_matrix(d, targets, wts, targets=targets) + nposM = int((rawM > 0).sum()) + wiM, acM, _, _ = tripod_split(densM) + print(f"--- {label} motor<-motor ---") + print(fmt_matrix(rawM)) + print(f" positive entries: {nposM}/36 density within={wiM:.6g} across={acM:.6g} " + f"across-minus-within={acM-wiM:+.6g}") + + perm = tripod_permutation_control(dens) + rank = [i for i, (_, c) in enumerate(perm) + if c == tuple(sorted(TRIPOD_A))][0] + print(f" tripod permutation control on premotor<-premotor density " + f"(across-minus-within, all 10 splits):") + for i, (v, c) in enumerate(perm): + mark = " <== TRUE TRIPOD" if c == tuple(sorted(TRIPOD_A)) else "" + print(f" {i+1:2d}. {v:+.6g} {'/'.join(c)}{mark}") + print(f" true tripod rank {rank+1}/10") + + +# ----------------------------------------------------------- LIF sweep (3,4) +def build_outgoing(d): + """Outgoing (CSC-like) view: for each presynaptic neuron, its edges.""" + N, E = d["N"], d["E"] + pre = d["col_idx"] + post = post_of_edges(d["row_ptr"], E) + order = np.argsort(pre, kind="stable") + pre_s = pre[order] + rows = post[order].astype(np.int32) + eidx = order # maps outgoing slot -> original edge index + indptr = np.zeros(N + 1, dtype=np.int64) + np.add.at(indptr, pre_s + 1, 1) + np.cumsum(indptr, out=indptr) + return indptr, rows, eidx + + +def simulate(indptr, rows, w_out, ext, N, n_steps, rec_idx, csr=None, + noise=0.0, seed=0): + """csr: optional scipy CSR of the same signed weights in INCOMING form. + Used instead of the sparse gather when the network is densely active — + identical arithmetic, just a faster route for the saturated regime.""" + vm = np.full(N, V_REST, np.float32) + gx = np.zeros(N, np.float32) + gy = np.zeros(N, np.float32) + refrac = np.zeros(N, np.int32) + spikes = np.zeros(N, bool) + ext_in = (ext * EXT_GAIN).astype(np.float32) + rec = np.zeros((n_steps, len(rec_idx)), np.uint8) + pop = np.zeros(n_steps, np.int32) + arange_N = np.arange(N) + sv = np.zeros(N, np.float32) + rng = np.random.default_rng(seed) + for t in range(n_steps): + act = np.flatnonzero(spikes) + if csr is not None and act.size > 0.12 * N: + sv[:] = spikes + delta = csr.dot(sv).astype(np.float32) + elif act.size: + s = indptr[act] + c = indptr[act + 1] - s + tot = int(c.sum()) + if tot: + csum = np.cumsum(c) + base = np.repeat(s - (csum - c), c) + slots = base + np.arange(tot, dtype=np.int64) + delta = np.bincount(rows[slots], weights=w_out[slots], + minlength=N).astype(np.float32) + else: + delta = np.zeros(N, np.float32) + else: + delta = np.zeros(N, np.float32) + gy = gy * A_SYN + gx + gx = gx * A_SYN + delta + i_in = gy * W_SYN + ext_in + if noise > 0: + i_in = i_in + rng.normal(0.0, noise, N).astype(np.float32) + r = refrac > 0 + vm[r] = V_RESET + refrac[r] -= 1 + nr = ~r + idx_nr = arange_N[nr] + v = np.float32(V_REST) + np.float32(ALPHA) * (vm[idx_nr] - np.float32(V_REST)) + i_in[idx_nr] + sp = v >= V_TH + v[sp] = V_RESET + vm[idx_nr] = v + fired = idx_nr[sp] + spikes = np.zeros(N, bool) + spikes[fired] = True + refrac[fired] = REFRAC_STEPS + rec[t] = spikes[rec_idx] + pop[t] = len(fired) + return rec, pop + + +def boxcar(x, k=20): + """Sliding-window mean over axis 0, matching captureRollingRate(20).""" + cs = np.cumsum(np.vstack([np.zeros((1,) + x.shape[1:]), x]), axis=0) + return (cs[k:] - cs[:-k]) / k + + +def spectrum_stats(sig, fs=1000.0): + """Return (band_fraction 3-25 Hz, peak_freq_Hz, peak_bin).""" + x = np.asarray(sig, float) + x = x - x.mean() + if not np.any(np.abs(x) > 0): + return 0.0, 0.0, 0 + T = len(x) + win = np.hanning(T) + X = np.fft.rfft(x * win) + P = np.abs(X) ** 2 + f = np.fft.rfftfreq(T, d=1.0 / fs) + tot = P[1:].sum() + if tot <= 0: + return 0.0, 0.0, 0 + band = P[(f >= 3.0) & (f <= 25.0)].sum() + pk = int(np.argmax(P[1:]) + 1) + return float(band / tot), float(f[pk]), pk + + +def pearson(a, b): + a = np.asarray(a, float); b = np.asarray(b, float) + a = a - a.mean(); b = b - b.mean() + da, db = np.sqrt((a * a).sum()), np.sqrt((b * b).sum()) + if da < 1e-12 or db < 1e-12: + return np.nan + return float((a * b).sum() / (da * db)) + + +def stage_sweep(d, use_glut_fix=False, dn_set="walking", quick=False): + meta = d["meta"] + N = d["N"] + indptr, rows, eidx = build_outgoing(d) + w_base = (glut_fix_weights(d["weight"], d["col_idx"], d["cell_class"])[0] + if use_glut_fix else d["weight"])[eidx].astype(np.float32) + neg = w_base < 0 + + # ---- DN drive targets + brain_meta = json.loads((ROOT / "public" / "brain.meta.json").read_text()) + if dn_set == "walk": + # forward/backward walking command DNs only (no escape DNs) + names = ["DNa01", "DNa02", "DNb01", "DNg13", "MDN"] + elif dn_set == "famous": + names = list(brain_meta["famous_dns"].keys()) + list(brain_meta["walking_circuit_dns"].keys()) + else: + names = None + if names is None: + dn_idx = np.flatnonzero(d["cell_class"] == CLASS_DN_INPUT) + dn_label = f"ALL {len(dn_idx)} DN-input neurons" + else: + got = [] + used = [] + for n in names: + if n in meta["dn_inputs"]: + got.extend(meta["dn_inputs"][n]); used.append(n) + dn_idx = np.unique(np.array(got, dtype=np.int64)) + dn_label = f"{len(dn_idx)} neurons from {len(used)} DN types: {','.join(sorted(used))}" + + # ---- recorded neurons: the 369 leg motor neurons + leg_all = {k: np.array(meta["motor"][k]["all"], dtype=np.int64) for k in LEG_KEYS} + rec_idx = np.unique(np.concatenate([leg_all[k] for k in LEG_KEYS])) + pos_of = {int(g): i for i, g in enumerate(rec_idx)} + leg_local = {k: np.array([pos_of[int(i)] for i in v]) for k, v in leg_all.items()} + musc_local = {} + for k in LEG_KEYS: + musc_local[k] = {} + for sub, idxs in meta["motor"][k].items(): + if sub == "all": + continue + musc_local[k][sub] = np.array([pos_of[int(i)] for i in idxs]) + + import os + W_SCALES = [0.1, 0.3, 1.0, 3.0, 10.0] if not quick else [1.0] + if os.environ.get("VNC_WS"): + W_SCALES = [float(x) for x in os.environ["VNC_WS"].split(",")] + DRIVES = [0.1, 0.2, 0.35, 0.5, 1.0, 2.0] if not quick else [0.5] + EI = [0.25, 0.5, 1.0, 2.0, 4.0] if not quick else [1.0] + WARM, T = 500, 3000 + if quick: + WARM, T = 200, 1000 + + print(f"# sweep glut_fix={use_glut_fix} dn={dn_label}") + print(f"# w_scale={W_SCALES} drive={DRIVES} ei={EI} warm={WARM} T={T} (dt=1ms)") + print("# LOOSE rhythm := band(3-25Hz)/total > 0.5 AND peak bin != 0 AND tripod corr < -0.3") + print("# STRICT rhythm := same but peak freq must itself lie in [3,25] Hz") + hdr = ("wS drive ei | popHz legHz popBand popPk | leg: corr band pkHz raw_band | " + "ratio: joint corr band pkHz | extflexcorr | L_leg L_mus S_leg S_mus") + print(hdr) + results = [] + t0 = time.time() + from scipy.sparse import csr_matrix + w_in_base = (glut_fix_weights(d["weight"], d["col_idx"], d["cell_class"])[0] + if use_glut_fix else d["weight"]).astype(np.float32) + neg_in = w_in_base < 0 + for ws in W_SCALES: + w_run = (w_base * ws).astype(np.float32) + w_in_ws = (w_in_base * ws).astype(np.float32) + for ei in EI: + w_ei = w_run.copy() + w_ei[neg] *= ei + w_in = w_in_ws.copy() + w_in[neg_in] *= ei + csr = csr_matrix((w_in, d["col_idx"], d["row_ptr"]), shape=(N, N)) + for dr in DRIVES: + ext = np.zeros(N, np.float32) + ext[dn_idx] = dr + rec, pop = simulate(indptr, rows, w_ei, ext, N, WARM + T, rec_idx, csr=csr) + rec = rec[WARM:] + pop = pop[WARM:] + rf = rec.astype(np.float32) + r = boxcar(rf, 20) # (T-19, n_mn) + pop_hz = pop.mean() / N * 1000.0 + pop_band, pop_pk, _ = spectrum_stats(boxcar(pop.astype(float)[:, None], 20)[:, 0]) + leg_tr = {k: r[:, leg_local[k]].mean(axis=1) for k in LEG_KEYS} + leg_raw = {k: rf[:, leg_local[k]].mean(axis=1) for k in LEG_KEYS} + leg_hz = np.mean([leg_tr[k].mean() for k in LEG_KEYS]) * 1000.0 + A = np.mean([leg_tr[k] for k in TRIPOD_A], axis=0) + B = np.mean([leg_tr[k] for k in TRIPOD_B], axis=0) + corr_leg = pearson(A, B) + bands, pks, rawbands = [], [], [] + for k in LEG_KEYS: + b, f, pb = spectrum_stats(leg_tr[k]) + bands.append(b); pks.append((f, pb)) + rawbands.append(spectrum_stats(leg_raw[k])[0]) + band_leg = float(np.mean(bands)) + rawband_leg = float(np.mean(rawbands)) + pk_leg = float(np.median([p[0] for p in pks])) + pkbin_ok_leg = all(p[1] > 0 for p in pks) + pkfreq_ok_leg = all(3.0 <= p[0] <= 25.0 for p in pks) + + # ---- muscle-resolved + def pool(k, subs): + ii = [musc_local[k][s] for s in subs if s in musc_local[k]] + if not ii: + return np.zeros(r.shape[0]), 0 + ii = np.concatenate(ii) + return r[:, ii].mean(axis=1), len(ii) + + ratio_tr = {} + extflex = [] + for k in LEG_KEYS: + for joint, (ex, fl) in ANTAGONISTS.items(): + ae, ne = pool(k, ex) + af, nf = pool(k, fl) + ratio_tr[(k, joint)] = (ae - af) / (ae + af + 1e-9) + c = pearson(ae, af) + if not np.isnan(c): + extflex.append(c) + best = None + for joint in ANTAGONISTS: + RA = np.mean([ratio_tr[(k, joint)] for k in TRIPOD_A], axis=0) + RB = np.mean([ratio_tr[(k, joint)] for k in TRIPOD_B], axis=0) + c = pearson(RA, RB) + bb, ff, pb = [], [], [] + for k in LEG_KEYS: + x, y, z = spectrum_stats(ratio_tr[(k, joint)]) + bb.append(x); ff.append(y); pb.append(z) + entry = dict(joint=joint, corr=c, band=float(np.mean(bb)), + pk=float(np.median(ff)), + pkbin_ok=all(p > 0 for p in pb), + pkfreq_ok=all(3.0 <= q <= 25.0 for q in ff)) + loose = (not np.isnan(c)) and c < -0.3 and entry["band"] > 0.5 and entry["pkbin_ok"] + strict = loose and entry["pkfreq_ok"] + score = (strict, loose, -(c if not np.isnan(c) else 0)) + if best is None or score > best[0]: + best = (score, entry) + e = best[1] + joint, corr_m, band_m, pk_m = e["joint"], e["corr"], e["band"], e["pk"] + ef = float(np.mean(extflex)) if extflex else float("nan") + + L_leg = (not np.isnan(corr_leg)) and corr_leg < -0.3 and band_leg > 0.5 and pkbin_ok_leg + S_leg = L_leg and pkfreq_ok_leg + L_mus = best[0][1] + S_mus = best[0][0] + print(f"{ws:4g} {dr:5g} {ei:4g} | {pop_hz:6.2f} {leg_hz:6.2f} " + f"{pop_band:.3f} {pop_pk:6.2f} | " + f"{corr_leg:+.3f} {band_leg:.3f} {pk_leg:6.2f} {rawband_leg:.3f} | " + f"{joint:>10s} {corr_m:+.3f} {band_m:.3f} {pk_m:6.2f} | " + f"{ef:+.3f} | {int(L_leg)} {int(L_mus)} {int(S_leg)} {int(S_mus)}", + flush=True) + results.append(dict(ws=ws, dr=dr, ei=ei, pop_hz=pop_hz, leg_hz=leg_hz, + pop_band=pop_band, pop_pk=pop_pk, + corr_leg=corr_leg, band_leg=band_leg, pk_leg=pk_leg, + rawband_leg=rawband_leg, + joint=joint, corr_m=corr_m, band_m=band_m, pk_m=pk_m, + extflex=ef, L_leg=bool(L_leg), L_mus=bool(L_mus), + S_leg=bool(S_leg), S_mus=bool(S_mus))) + n = len(results) + alive = [r for r in results if r["leg_hz"] > 0.05] + cl = [r["corr_leg"] for r in alive if not np.isnan(r["corr_leg"])] + cm = [r["corr_m"] for r in alive if not np.isnan(r["corr_m"])] + ce = [r["extflex"] for r in alive if not np.isnan(r["extflex"])] + print(f"\n#### {n} settings, {len(alive)} with non-silent leg motor pools " + f"({100*len(alive)/n:.1f}%)") + for key, lbl in (("L_leg", "LOOSE whole-leg"), ("L_mus", "LOOSE muscle-resolved"), + ("S_leg", "STRICT whole-leg"), ("S_mus", "STRICT muscle-resolved")): + c = sum(r[key] for r in results) + print(f"#### rhythm {lbl:24s}: {c}/{n} = {100*c/n:.2f}%") + for key, lbl in (("corr_leg", "tripod corr whole-leg"), ("corr_m", "tripod corr muscle"), + ("extflex", "within-leg ext-vs-flex corr")): + v = [r[key] for r in alive if not np.isnan(r[key])] + if v: + print(f"#### {lbl:28s} range [{min(v):+.3f}, {max(v):+.3f}] median {np.median(v):+.3f}") + # component-wise pass rates: which criterion is the binding constraint? + for lbl, f in (("corr_leg < -0.3", lambda r: r["corr_leg"] < -0.3), + ("corr_musc < -0.3", lambda r: r["corr_m"] < -0.3), + ("band_leg > 0.5", lambda r: r["band_leg"] > 0.5), + ("band_musc > 0.5", lambda r: r["band_m"] > 0.5), + ("rawband_leg > 0.5", lambda r: r["rawband_leg"] > 0.5), + ("pk_leg in [3,25]", lambda r: 3 <= r["pk_leg"] <= 25), + ("pk_musc in [3,25]", lambda r: 3 <= r["pk_m"] <= 25)): + c = sum(1 for r in results if not np.isnan(r["corr_leg"]) and f(r)) + print(f"#### criterion {lbl:20s} passes {c}/{n}") + print(f"#### wall clock {time.time()-t0:.1f}s") + return results + + +if __name__ == "__main__": + stage = sys.argv[1] if len(sys.argv) > 1 else "struct" + d = load_vnc() + if stage == "struct": + stage_struct(d) + elif stage == "sweep": + gf = "--glutfix" in sys.argv + dn = "all" if "--all-dn" in sys.argv else ("walk" if "--walk" in sys.argv else "famous") + stage_sweep(d, use_glut_fix=gf, dn_set=dn, quick="--quick" in sys.argv)