Decoding analyses of Watters, Gabel, Tenenbaum & Jazayeri (2026) on DANDI 000620¶
Paper: doi:10.64898/2026.01.27.702062 (bioRxiv preprint version 1); claims file claims/000620/10.64898-2026.01.27.702062.yaml; paper code https://github.com/jazlab/multi_object_memory_2025 at commit 6e4b333b78cd4461e6373a42378535e8fd392cdb (figures/supp_fig_decoding/). Dataset DANDI 000620 version 0.261019.2000 (doi:10.48324/dandi.000620/0.261019.2000), pinned as the most recent published version because the paper cites none and all published versions hold the same assets (see the main notebook, reproduce_watters_2026.ipynb, whose header gives the full reasoning and the table of all claims).
This notebook tests ext3ab-decoding and ext3c-decoder-orthogonality (Extended Data Fig. 3). The decoders are the paper's (models.py), written in numpy (wm.train_linear, wm.train_linear_2obj) with the configs' optimisers, learning rates, L1 weights, batch size and step counts, batched over the ten splits.
| claim | status |
|---|---|
| ext3ab-decoding | tested |
| ext3c-decoder-orthogonality | tested |
Variants that are not the code's (z-scored inputs, and for the Ring decoders without the output sigmoid) are in decoding_variants_watters_2026.ipynb; they do not decide an outcome.
import warnings
warnings.filterwarnings("ignore")
import time
import numpy as np, pandas as pd
import wm
t_start = time.time()
assets = wm.list_assets()
print(f"Dandiset {wm.DANDISET} version {wm.VERSION}: {len(assets)} assets,", assets.kind.value_counts().to_dict())
print("Raw ecephys files (not used):", f"{assets[assets.kind == 'ecephys']['size'].sum() / 1e12:.2f} TB")
keys = sorted(set(zip(assets.subject, assets.session)))
S = [wm.load_session(assets, subj, sess) for subj, sess in keys]
print(f"Loaded {len(S)} sessions in {time.time() - t_start:.0f} s")
Dandiset 000620 version 0.261019.2000: 142 assets, {'behavior+task': 48, 'ecephys': 47, 'spikesorting': 47}
Raw ecephys files (not used): 5.47 TB
Loaded 48 sessions in 268 s
import numpy as np, pandas as pd, scipy.stats as st
import matplotlib
import matplotlib.pyplot as plt
pd.set_option("display.width", 200); pd.set_option("display.max_columns", 30)
MONKEY = {"Perle": "P", "Elgar": "E"}
for x in S:
x["trials"]["subject"] = x["subject"]; x["trials"]["session"] = x["session"]
x["trials"]["task"] = x["task"]; x["trials"]["trial"] = np.arange(len(x["trials"]))
# Physiology sessions of each dataset: sessions with spike sorting. Ring sessions are those
# with 1- and 2-object trials (Methods, Behavioral Task); Perle 2022-06-06 has 2-object
# trials only and is kept apart (see the data check).
tri = [x for x in S if x["task"] == "triangle" and x["units"] is not None]
ring = [x for x in S if x["task"] == "ring" and x["units"] is not None
and (x["trials"].n_obj == 1).any()]
other = [x for x in S if x not in tri and x not in ring]
T = pd.concat([x["trials"] for x in S], ignore_index=True)
T["physiology"] = T.session.isin([x["session"] for x in tri + ring])
print("Triangle sessions:", {m: sum(x["subject"] == m for x in tri) for m in MONKEY})
print("Ring sessions:", {m: sum(x["subject"] == m for x in ring) for m in MONKEY})
print("Other sessions:", [(x["subject"], x["session"], x["task"], x["units"] is not None,
sorted(x["trials"].n_obj.unique())) for x in other])
Triangle sessions: {'Perle': 10, 'Elgar': 14}
Ring sessions: {'Perle': 10, 'Elgar': 12}
Other sessions: [('Perle', '2022-06-06', 'ring', True, [np.int64(2)]), ('Perle', '2022-06-08', 'ring', False, [np.int64(1), np.int64(2)])]
Data check¶
inv = pd.DataFrame([dict(monkey=MONKEY[x["subject"]], session=x["session"], task=x["task"], trials=len(x["trials"]),
completed=int(x["trials"].completed.sum()),
n1=int((x["trials"].n_obj == 1).sum()), n2=int((x["trials"].n_obj == 2).sum()), n3=int((x["trials"].n_obj == 3).sum()),
units=np.nan if x["units"] is None else len(x["units"])) for x in S])
print(inv.to_string(index=False))
r = np.concatenate([np.nanmean(x["rates"][w], 1) for x in tri + ring for w in ["stimulus", "delay"]])
print("Unit mean rates, stimulus and delay phases (Hz): min", r.min().round(3), "median", np.median(r).round(2), "max", r.max().round(1))
xy = T[["x0", "y0"]].dropna().values
print("Object positions (arena widths from centre): eccentricity", np.unique(np.hypot(*xy.T).round(4)))
print("Ring 2-object positions in decoder units (2x, 2y): range", (2 * T[T.task == "ring"][["x0", "y0", "x1", "y1"]]).min().min().round(3),
"to", (2 * T[T.task == "ring"][["x0", "y0", "x1", "y1"]]).max().max().round(3))
print("The main notebook (reproduce_watters_2026.ipynb) checks the clock and the metadata in full.")
monkey session task trials completed n1 n2 n3 units
E 2022-08-19 triangle 2847 1253 234 928 1685 323.0
E 2022-08-20 triangle 2619 1361 251 991 1377 414.0
E 2022-08-21 triangle 2928 1517 289 992 1647 229.0
E 2022-08-22 triangle 1964 1157 187 722 1055 234.0
E 2022-08-23 triangle 1816 1299 149 631 1036 217.0
E 2022-08-24 triangle 1401 849 130 501 770 283.0
E 2022-08-25 triangle 2083 1353 205 768 1110 164.0
E 2022-08-26 triangle 2132 1514 240 845 1047 141.0
E 2022-08-31 triangle 1578 1238 154 595 829 206.0
E 2022-09-01 triangle 1947 1373 189 720 1038 225.0
E 2022-09-02 triangle 1512 1193 148 574 790 151.0
E 2022-09-03 triangle 2267 1735 240 931 1096 203.0
E 2022-09-04 triangle 1954 1336 200 751 1003 290.0
E 2022-09-05 triangle 1598 1238 172 621 805 181.0
E 2022-09-06 ring 2080 1335 389 1691 0 197.0
E 2022-09-07 ring 2307 1297 437 1870 0 78.0
E 2022-09-09 ring 2186 1656 481 1705 0 127.0
E 2022-09-10 ring 1828 1448 376 1452 0 165.0
E 2022-09-11 ring 2435 1722 488 1947 0 148.0
E 2022-09-12 ring 2008 1402 428 1580 0 134.0
E 2022-09-13 ring 1922 1361 387 1535 0 71.0
E 2022-09-14 ring 2019 1350 414 1605 0 98.0
E 2022-09-15 ring 1849 1339 384 1465 0 93.0
E 2022-09-18 ring 1722 1323 315 1407 0 77.0
E 2022-09-19 ring 1598 1214 322 1276 0 114.0
E 2022-09-20 ring 1680 1242 306 1374 0 63.0
P 2022-05-26 triangle 1756 1002 241 587 928 101.0
P 2022-05-27 triangle 1629 1239 200 602 827 111.0
P 2022-05-28 triangle 1840 1412 198 711 931 137.0
P 2022-05-29 triangle 2040 1326 275 768 997 129.0
P 2022-05-30 triangle 1555 992 223 545 787 185.0
P 2022-05-31 triangle 1481 957 172 555 754 200.0
P 2022-06-01 triangle 1621 1097 209 648 764 226.0
P 2022-06-03 triangle 2091 1335 270 791 1030 198.0
P 2022-06-04 triangle 1862 1123 234 687 941 243.0
P 2022-06-05 triangle 1338 896 193 507 638 108.0
P 2022-06-06 ring 2378 1559 0 2378 0 77.0
P 2022-06-07 ring 2660 1481 538 2122 0 14.0
P 2022-06-08 ring 1517 1001 281 1236 0 NaN
P 2022-06-09 ring 2417 1416 488 1929 0 193.0
P 2022-06-10 ring 1476 1206 295 1181 0 216.0
P 2022-06-11 ring 1474 1170 280 1194 0 207.0
P 2022-06-12 ring 1724 1275 330 1394 0 231.0
P 2022-06-13 ring 1365 968 257 1108 0 232.0
P 2022-06-14 ring 1762 965 347 1415 0 182.0
P 2022-06-15 ring 2291 1516 465 1826 0 136.0
P 2022-06-16 ring 1372 968 283 1089 0 217.0
P 2022-06-19 ring 3678 1645 748 2930 0 193.0
Unit mean rates, stimulus and delay phases (Hz): min 0.009 median 6.0 max 133.9
Object positions (arena widths from centre): eccentricity [0.35]
Ring 2-object positions in decoder units (2x, 2y): range -0.7 to 0.7
The main notebook (reproduce_watters_2026.ipynb) checks the clock and the metadata in full.
Sessions and units of the decoders¶
# Units entering the models (T7): delay-phase AUC > 0.65 with p < 0.05 for position or identity,
# rate > 1 Hz. AUC on 1-object trials as in the main notebook (T6), Methods p-value (1,000
# shuffles); the sessions it ranks highest are the ones the code's decoding scripts name.
rng_sel = np.random.default_rng(620)
ID_NUM = {"a": 0, "b": 1, "c": 2}
sel = []
for x in tri + ring:
tr = x["trials"]; R = x["rates"]["delay"]
m = (tr.n_obj == 1).values & tr.completed.values
if x["task"] == "triangle":
m &= tr.on_triangle.values; pos = tr.loc0.values
else:
pos = wm.ring_sector(tr.theta0.values)
idn = tr.id0.map(ID_NUM).fillna(-1).astype(int).values
rate = np.nanmean(R, 1)
for k in np.where(rate > wm.MIN_RATE_HZ)[0]:
r = R[k]; ok = m & np.isfinite(r)
pa, pp = wm.unit_auc(r[ok], pos[ok], rng_sel, "methods")
ia, ip = wm.unit_auc(r[ok], idn[ok], rng_sel, "methods")
if (pa > 0.65 and pp < 0.05) or (ia > 0.65 and ip < 0.05):
sel.append(dict(subject=x["subject"], session=x["session"], task=x["task"], unit_index=k))
sel = pd.DataFrame(sel)
n_mod = sel.groupby(["task", "subject", "session"]).size().rename("modeled_units").reset_index()
top = {t: n_mod[n_mod.task == t].sort_values("modeled_units", ascending=False).head(3) for t in ["triangle", "ring"]}
for t in top:
print(f"{t}, three sessions with the most modeled units:"); print(top[t].to_string(index=False))
print("Code's sessions: Triangle Perle 2022-06-01, 2022-06-04, 2022-05-31 (run_1_triangle.py); Ring Perle 2022-06-09, 2022-06-10, 2022-06-11 (run_2_ring.py)")
SESS = {x["session"]: x for x in S}
UNITS = {s: g.unit_index.values for s, g in sel.groupby("session")}
triangle, three sessions with the most modeled units:
task subject session modeled_units
triangle Perle 2022-06-01 171
triangle Perle 2022-05-31 117
triangle Perle 2022-06-04 116
ring, three sessions with the most modeled units:
task subject session modeled_units
ring Perle 2022-06-10 113
ring Perle 2022-06-11 113
ring Perle 2022-06-09 89
Code's sessions: Triangle Perle 2022-06-01, 2022-06-04, 2022-05-31 (run_1_triangle.py); Ring Perle 2022-06-09, 2022-06-10, 2022-06-11 (run_2_ring.py)
ext3ab-decoding¶
Single-trial delay-phase activity carries information about object positions and identities, decoded above a shuffled-label control.
Source: Results, The Triangle Dataset, paragraph 1; Extended Data Fig. 3a and 3b. Paper: "We also found that neural activity contained information about object positions and identities ( Extended Data Fig. 3a-b )"
Criterion (claims file): As the Extended Data Fig. 3 legend and figures/supp_fig_decoding/ fix it, on the three Triangle sessions with the most units entering the models (T7; the code's analysis notebook uses Perle 2022-05-31, 2022-06-01 and 2022-06-04): a linear decoder (MSE loss to a one-hot target, L1 weight 0.003, SGD at learning rate 0.01, 10,000 steps, config_1_triangle.py) from single-trial rates of the stimulus and delay phases, ten 80/20 train/test splits, against the same with shuffled labels. Panel a decodes the location or identity of 1-object trials; panel b, on 2-object trials, the absent location or identity. The paper sets no rule for "above control"; readings (a) the neural accuracy exceeds the shuffled accuracy in all ten splits; (b) the mean neural accuracy exceeds the shuffled mean by more than two standard deviations of the shuffled splits. Reproduced when, under both readings, delay-phase position and identity decoding of panel a are above control in all three sessions; panel b is reported, with the legend's 66% mark.
# ext3ab-decoding. Panel a: 1-object Triangle trials, location or identity, SGD lr 0.01
# (config_1_triangle.py); panel b: 2-object trials, the absent location or identity, Adam lr
# 0.001 (config_2_triangle_complement.py). Ten 80/20 splits (seeds 0 to 9) on neural data
# and on data with shuffled labels; 10,000 steps, L1 0.003, batches of 512. Own generator.
# Inputs are rates in Hz, as the code's: dataset.py leaves its per-unit normalisation
# commented out. A z-scored variant is in decoding_variants_watters_2026.ipynb.
rng_dec = np.random.default_rng(3)
dec, filt = [], []
for kind, opt, lr in [("1", "sgd", 0.01), ("complement", "adam", 0.001)]:
for s in top["triangle"].session:
for phase in ["stimulus", "delay"]:
for feature in ["position", "identity"]:
X, Y, cnt = wm.decoding_data(SESS[s], UNITS[s], phase, kind, feature, info=True)
filt.append(dict(panel="a" if kind == "1" else "b", session=s, phase=phase, feature=feature, **cnt))
n = len(X); nt = int(0.2 * n)
Xs, Ys, Xt, Yt = [], [], [], []
for f in range(20):
r = np.random.default_rng(f % 10)
y = Y[np.random.default_rng(100 + f).permutation(n)] if f >= 10 else Y
perm = r.permutation(n)
Xs.append(X[perm[nt:]]); Ys.append(y[perm[nt:]]); Xt.append(X[perm[:nt]]); Yt.append(y[perm[:nt]])
W, b = wm.train_linear(np.array(Xs, np.float32), np.array(Ys, np.float32), rng_dec, opt=opt, lr=lr)
for f in range(20):
acc = float(np.mean(np.argmax(Xt[f] @ W[f] + b[f], 1) == np.argmax(Yt[f], 1)))
tr_acc = float(np.mean(np.argmax(Xs[f] @ W[f] + b[f], 1) == np.argmax(Ys[f], 1)))
dec.append(dict(panel="a" if kind == "1" else "b", session=s, phase=phase, feature=feature,
labels="shuffled" if f >= 10 else "neural", split=f % 10, accuracy=100 * acc,
train_accuracy=100 * tr_acc, n_trials=n, n_units=X.shape[1]))
dec = pd.DataFrame(dec)
filt = pd.DataFrame(filt)
print("Units and trials after each of the code's filters (units entering the models, then units with >= 50 trials,")
print("trials with >= 5 units, then the greedy simultaneity filter of dataset.py):")
print(filt.drop_duplicates(["panel", "session", "phase"]).drop(columns="feature").to_string(index=False))
summ = dec.groupby(["panel", "session", "phase", "feature", "labels"]).accuracy.agg(["mean", "std", "min", "max"]).unstack("labels")
print("\nTest accuracy (%):"); print(summ.round(1).to_string())
print("\nTraining accuracy (%), mean over splits:")
print(dec.groupby(["panel", "phase", "feature", "session", "labels"]).train_accuracy.mean().unstack("labels").round(1).to_string())
judge = []
for (p, s, ph, ft), g in dec.groupby(["panel", "session", "phase", "feature"]):
nn = g[g.labels == "neural"].sort_values("split").accuracy.values
sh = g[g.labels == "shuffled"].sort_values("split").accuracy.values
judge.append(dict(panel=p, session=s, phase=ph, feature=ft, neural_mean=nn.mean(), shuffled_mean=sh.mean(),
a_all_splits=bool((nn > sh).all()), b_two_sd=bool(nn.mean() > sh.mean() + 2 * sh.std(ddof=1))))
judge = pd.DataFrame(judge)
print(); print(judge.round(1).to_string(index=False))
ja = judge[(judge.panel == "a") & (judge.phase == "delay")]
print("\nPanel a, delay phase, position and identity above control in all three sessions: reading (a)",
bool(ja.a_all_splits.all()), "; reading (b)", bool(ja.b_two_sd.all()))
for ft in ["position", "identity"]:
jf = ja[ja.feature == ft]
print(f" {ft}: reading (a) in {int(jf.a_all_splits.sum())} of 3 sessions, reading (b) in {int(jf.b_two_sd.sum())} of 3")
jb = judge[judge.panel == "b"]
print("Panel b, mean neural accuracy above 66%:")
print(jb[["session", "phase", "feature", "neural_mean"]].assign(over_66=jb.neural_mean > 66).round(1).to_string(index=False))
fig, axes = plt.subplots(2, 2, figsize=(8, 5), sharey=True)
for i, p in enumerate(["a", "b"]):
for j, ft in enumerate(["position", "identity"]):
g = dec[(dec.panel == p) & (dec.feature == ft)].groupby(["phase", "labels"]).accuracy.mean().unstack()
g.plot.bar(ax=axes[i, j], rot=0); axes[i, j].set_title(f"panel {p}, {ft}"); axes[i, j].set_ylim(0, 100)
plt.tight_layout(); plt.show()
print("All accuracies lie in [0, 100]:", bool(dec[["accuracy", "train_accuracy"]].stack().between(0, 100).all()))
Units and trials after each of the code's filters (units entering the models, then units with >= 50 trials,
trials with >= 5 units, then the greedy simultaneity filter of dataset.py):
panel session phase units_in trials_in units_min_trials trials_min_units units_simultaneous trials_simultaneous
a 2022-06-01 stimulus 171 121 161 121 104 121
a 2022-06-01 delay 171 121 161 121 104 121
a 2022-05-31 stimulus 117 105 110 105 82 105
a 2022-05-31 delay 117 105 110 105 82 105
a 2022-06-04 stimulus 116 121 72 121 20 121
a 2022-06-04 delay 116 121 72 121 20 121
b 2022-06-01 stimulus 171 391 171 391 104 387
b 2022-06-01 delay 171 391 171 391 104 387
b 2022-05-31 stimulus 117 329 117 329 85 318
b 2022-05-31 delay 117 329 117 329 85 318
b 2022-06-04 stimulus 116 392 115 392 20 392
b 2022-06-04 delay 116 392 115 392 20 392
Test accuracy (%):
mean std min max
labels neural shuffled neural shuffled neural shuffled neural shuffled
panel session phase feature
a 2022-05-31 delay identity 46.2 32.4 3.2 9.7 42.9 14.3 52.4 47.6
position 95.7 33.3 2.7 11.2 90.5 14.3 100.0 47.6
stimulus identity 40.0 33.8 8.2 10.9 28.6 19.0 57.1 57.1
position 97.1 36.7 3.3 9.5 90.5 23.8 100.0 52.4
2022-06-01 delay identity 50.0 37.5 11.3 9.6 29.2 25.0 62.5 54.2
position 99.2 33.8 1.8 11.7 95.8 12.5 100.0 50.0
stimulus identity 42.5 32.5 8.1 6.7 29.2 20.8 54.2 41.7
position 99.6 32.9 1.3 16.9 95.8 4.2 100.0 62.5
2022-06-04 delay identity 30.0 30.8 7.6 8.4 12.5 16.7 37.5 45.8
position 82.5 32.1 7.0 10.0 70.8 16.7 95.8 50.0
stimulus identity 47.9 33.3 12.5 9.0 29.2 20.8 62.5 50.0
position 72.1 29.6 8.1 6.6 58.3 20.8 83.3 45.8
b 2022-05-31 delay identity 34.6 35.6 6.1 6.9 20.6 27.0 42.9 49.2
position 98.9 34.6 1.3 3.0 96.8 30.2 100.0 39.7
stimulus identity 38.6 34.9 3.7 5.7 33.3 25.4 44.4 42.9
position 99.7 32.4 0.7 4.5 98.4 27.0 100.0 39.7
2022-06-01 delay identity 38.4 31.7 5.7 5.2 29.9 22.1 45.5 39.0
position 99.5 35.7 0.9 5.2 97.4 31.2 100.0 46.8
stimulus identity 43.6 35.2 6.9 7.3 36.4 22.1 54.5 44.2
position 99.1 32.1 1.4 4.7 96.1 22.1 100.0 39.0
2022-06-04 delay identity 32.4 34.1 5.3 4.3 24.4 25.6 42.3 41.0
position 80.1 34.6 4.8 5.7 73.1 23.1 88.5 39.7
stimulus identity 40.9 30.8 4.9 4.1 35.9 24.4 52.6 38.5
position 78.5 35.0 4.8 5.3 69.2 28.2 83.3 42.3
Training accuracy (%), mean over splits:
labels neural shuffled
panel phase feature session
a delay identity 2022-05-31 92.5 89.5
2022-06-01 82.7 84.1
2022-06-04 63.6 66.3
position 2022-05-31 100.0 90.4
2022-06-01 100.0 83.0
2022-06-04 95.5 68.0
stimulus identity 2022-05-31 88.1 83.6
2022-06-01 93.2 88.8
2022-06-04 69.4 65.4
position 2022-05-31 100.0 88.8
2022-06-01 100.0 88.4
2022-06-04 88.7 63.2
b delay identity 2022-05-31 77.8 75.2
2022-06-01 77.1 74.1
2022-06-04 45.4 47.5
position 2022-05-31 100.0 74.7
2022-06-01 100.0 76.4
2022-06-04 82.1 45.4
stimulus identity 2022-05-31 79.5 74.0
2022-06-01 81.6 73.9
2022-06-04 48.5 47.2
position 2022-05-31 100.0 76.0
2022-06-01 100.0 75.3
2022-06-04 81.9 45.9
panel session phase feature neural_mean shuffled_mean a_all_splits b_two_sd
a 2022-05-31 delay identity 46.2 32.4 False False
a 2022-05-31 delay position 95.7 33.3 True True
a 2022-05-31 stimulus identity 40.0 33.8 False False
a 2022-05-31 stimulus position 97.1 36.7 True True
a 2022-06-01 delay identity 50.0 37.5 False False
a 2022-06-01 delay position 99.2 33.8 True True
a 2022-06-01 stimulus identity 42.5 32.5 False False
a 2022-06-01 stimulus position 99.6 32.9 True True
a 2022-06-04 delay identity 30.0 30.8 False False
a 2022-06-04 delay position 82.5 32.1 True True
a 2022-06-04 stimulus identity 47.9 33.3 False False
a 2022-06-04 stimulus position 72.1 29.6 True True
b 2022-05-31 delay identity 34.6 35.6 False False
b 2022-05-31 delay position 98.9 34.6 True True
b 2022-05-31 stimulus identity 38.6 34.9 False False
b 2022-05-31 stimulus position 99.7 32.4 True True
b 2022-06-01 delay identity 38.4 31.7 False False
b 2022-06-01 delay position 99.5 35.7 True True
b 2022-06-01 stimulus identity 43.6 35.2 True False
b 2022-06-01 stimulus position 99.1 32.1 True True
b 2022-06-04 delay identity 32.4 34.1 False False
b 2022-06-04 delay position 80.1 34.6 True True
b 2022-06-04 stimulus identity 40.9 30.8 False True
b 2022-06-04 stimulus position 78.5 35.0 True True
Panel a, delay phase, position and identity above control in all three sessions: reading (a) False ; reading (b) False
position: reading (a) in 3 of 3 sessions, reading (b) in 3 of 3
identity: reading (a) in 0 of 3 sessions, reading (b) in 0 of 3
Panel b, mean neural accuracy above 66%:
session phase feature neural_mean over_66
2022-05-31 delay identity 34.6 False
2022-05-31 delay position 98.9 True
2022-05-31 stimulus identity 38.6 False
2022-05-31 stimulus position 99.7 True
2022-06-01 delay identity 38.4 False
2022-06-01 delay position 99.5 True
2022-06-01 stimulus identity 43.6 False
2022-06-01 stimulus position 99.1 True
2022-06-04 delay identity 32.4 False
2022-06-04 delay position 80.1 True
2022-06-04 stimulus identity 40.9 False
2022-06-04 stimulus position 78.5 True
All accuracies lie in [0, 100]: True
Result: partially reproduced. Delay-phase position decoding of 1-object trials is above the shuffled control in all ten splits and by more than two standard deviations in all three sessions (82.5% to 99.2% against 32.1% to 33.8%), and so is stimulus-phase position decoding. Identity decoding is above control in none of the three sessions under either reading (0 of 3 under each): in the delay phase it is 46.2% (2022-05-31), 50.0% (2022-06-01) and 30.0% (2022-06-04) against 32.4%, 37.5% and 30.8% shuffled. Identity decoders fit their training trials hardly better than decoders trained on shuffled labels (delay, 92.5% against 89.5%, 82.7% against 84.1%, 63.6% against 66.3%), so they overfit. In panel b the absent location is decoded above 66% in every session and phase (78.5% to 99.7%), the absent identity in none (32.4% to 43.6%). The code's filters leave 105 to 121 1-object trials per session and, for 2022-06-04, 20 of 116 units (72 have 50 trials; the greedy simultaneity filter then keeps 20 units seen on all 121 trials). With z-scored inputs (variants notebook) position still passes in all three sessions and identity in one session under each reading, so identity still fails. Position holds and identity does not; the criterion defines only the reproduced case, and partially reproduced is this attempt's label for that.
ext3c-decoder-orthogonality¶
The coding axes of a two-object decoder come out nearly orthogonal even for synthetic data with no slot structure, so decoder orthogonality does not diagnose slots.
Source: Methods, Decoding Analyses, paragraphs 1 and 4; Extended Data Fig. 3c. Paper: "However, we found that this orthogonality resulted even on synthetic data that was specifically not Slot-like. In other words, orthogonality of these decoding subspaces is not diagnostic of object slots."
Criterion (claims file): The 2-object decoder of Methods, Decoding Analyses (two linear decoders and a per-trial assignment sigma pressed to be binary; figures/supp_fig_decoding/modeling/models.py, LinearModel2Objects) fitted on 2-object trials of the three Triangle sessions of ext3ab-decoding and of the three Ring sessions with the most units entering the models (the Extended Data Fig. 3 legend picks three sessions from each dataset; the code's run_2_ring.py uses Perle 2022-06-09, 2022-06-10 and 2022-06-11), and on the Receptive Field and Orthogonal synthetic datasets built from them as Methods describe; angles between corresponding columns of W0 and W1 for each location, over ten bootstraps, stimulus and delay phases. The paper gives no threshold for "orthogonal". Reproduced when the angles for the neural data and for the Receptive Field data are both centred near 90 degrees, read as (a) a median between 75 and 105 degrees and (b) a median within 10 degrees of the Orthogonal dataset's median; both cut-offs are the agent's. Not reproduced when neither reading holds; methods_underspecified when they differ.
# ext3c-decoder-orthogonality. The 2-object decoder (models.py LinearModel2Objects; configs
# config_2_triangle.py and config_2_ring.py) on 2-object trials of the three Triangle and three
# Ring sessions, neural data and the Receptive Field and Orthogonal synthetic data
# (dataset.py), ten seeds (0 to 9) with an 80/20 split each, stimulus and delay phase. Angles
# between corresponding columns of the two decoders. Rates in Hz and sigmoid outputs, as the
# code; variants are in decoding_variants_watters_2026.ipynb. Own generator.
rng_orth = np.random.default_rng(4)
ang, filt2 = [], []
for task in ["triangle", "ring"]:
for s in top[task].session:
for phase in ["stimulus", "delay"]:
X0, Y, cnt = wm.decoding_data(SESS[s], UNITS[s], phase, "2obj", "position", info=True)
filt2.append(dict(task=task, session=s, phase=phase, **cnt))
n = len(X0); nt = int(0.2 * n)
for mode in ["neural", "receptive_fields", "orthogonal"]:
Xs, Ys = [], []
for seed in range(10):
r = np.random.default_rng(seed)
X = X0 if mode == "neural" else wm.synthetic_2obj(X0, Y, task, mode, np.random.default_rng(1000 + seed))
perm = r.permutation(n)
Xs.append(X[perm[nt:]]); Ys.append(Y[perm[nt:]])
W0, W1 = wm.train_linear_2obj(np.array(Xs, np.float32), np.array(Ys, np.float32), rng_orth)
a = wm.angles_deg(W0, W1)
for seed in range(10):
for i, v in enumerate(a[seed]):
ang.append(dict(task=task, session=s, phase=phase, mode=mode, seed=seed, feature=i, angle=v,
n_trials=n, n_units=X0.shape[1]))
ang = pd.DataFrame(ang)
filt2 = pd.DataFrame(filt2)
print("Units and trials after each of the code's filters (2-object trials):")
print(filt2.to_string(index=False))
med = ang.groupby(["task", "phase", "mode"]).angle.median().unstack("mode")
print("\nMedian angle (degrees) over sessions, seeds and decoder columns:"); print(med.round(1).to_string())
print("\nRange of angles:", ang.groupby(["task", "mode"]).angle.agg(["min", "max"]).round(1).to_string())
def judge3c(med):
rows = []
for (t, ph), g in med.iterrows():
for mode in ["neural", "receptive_fields"]:
rows.append(dict(task=t, phase=ph, data=mode, median=g[mode],
a_75_105=bool(75 <= g[mode] <= 105), b_within_10_of_orth=bool(abs(g[mode] - g["orthogonal"]) <= 10)))
return pd.DataFrame(rows)
res3c = judge3c(med)
print(res3c.round(1).to_string(index=False))
print("Reading (a) holds everywhere:", bool(res3c.a_75_105.all()), "; reading (b):", bool(res3c.b_within_10_of_orth.all()))
for t in ["triangle", "ring"]:
for d in ["neural", "receptive_fields"]:
r = res3c[(res3c.task == t) & (res3c.data == d)]
print(f" {t}, {d}: reading (a) {bool(r.a_75_105.all())}, reading (b) {bool(r.b_within_10_of_orth.all())}")
print("\nMedian angle per session:")
print(ang.groupby(["task", "phase", "session", "mode"]).angle.median().unstack("mode").round(1).to_string())
print("\nTriangle without Perle 2022-06-04 (20 units after the simultaneity filter):")
r = judge3c(ang[(ang.task == "triangle") & (ang.session != "2022-06-04")].groupby(["task", "phase", "mode"]).angle.median().unstack("mode"))
print(r.round(1).to_string(index=False))
fig, axes = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
for ax, t in zip(axes, ["triangle", "ring"]):
g = ang[ang.task == t]
data = [g[(g["mode"] == m) & (g.phase == ph)].angle for m in ["neural", "orthogonal", "receptive_fields"] for ph in ["stimulus", "delay"]]
ax.boxplot(data, tick_labels=[f"{m[:5]}\n{ph[:4]}" for m in ["neural", "orthogonal", "receptive_fields"] for ph in ["stimulus", "delay"]])
ax.set_title(t); ax.axhline(90, c="k", lw=0.5)
axes[0].set_ylabel("angle (deg)"); plt.show()
Units and trials after each of the code's filters (2-object trials):
task session phase units_in trials_in units_min_trials trials_min_units units_simultaneous trials_simultaneous
triangle 2022-06-01 stimulus 171 391 171 391 104 387
triangle 2022-06-01 delay 171 391 171 391 104 387
triangle 2022-05-31 stimulus 117 329 117 329 85 318
triangle 2022-05-31 delay 117 329 117 329 85 318
triangle 2022-06-04 stimulus 116 392 115 392 20 392
triangle 2022-06-04 delay 116 392 115 392 20 392
ring 2022-06-10 stimulus 113 1010 113 1010 70 1010
ring 2022-06-10 delay 113 1010 113 1010 70 1010
ring 2022-06-11 stimulus 113 981 113 981 90 974
ring 2022-06-11 delay 113 981 113 981 90 974
ring 2022-06-09 stimulus 89 1152 89 1152 44 405
ring 2022-06-09 delay 89 1152 89 1152 44 405
Median angle (degrees) over sessions, seeds and decoder columns:
mode neural orthogonal receptive_fields
task phase
ring delay 88.800003 92.500000 96.400002
stimulus 93.199997 95.000000 96.699997
triangle delay 106.699997 91.000000 91.900002
stimulus 106.599998 90.800003 91.800003
Range of angles: min max
task mode
ring neural 77.800003 123.900002
orthogonal 68.300003 139.399994
receptive_fields 75.000000 138.000000
triangle neural 81.000000 144.100006
orthogonal 63.099998 144.300003
receptive_fields 59.500000 137.800003
task phase data median a_75_105 b_within_10_of_orth
ring delay neural 88.800003 True True
ring delay receptive_fields 96.400002 True True
ring stimulus neural 93.199997 True True
ring stimulus receptive_fields 96.699997 True True
triangle delay neural 106.699997 False False
triangle delay receptive_fields 91.900002 True True
triangle stimulus neural 106.599998 False False
triangle stimulus receptive_fields 91.800003 True True
Reading (a) holds everywhere: False ; reading (b): False
triangle, neural: reading (a) False, reading (b) False
triangle, receptive_fields: reading (a) True, reading (b) True
ring, neural: reading (a) True, reading (b) True
ring, receptive_fields: reading (a) True, reading (b) True
Median angle per session:
mode neural orthogonal receptive_fields
task phase session
ring delay 2022-06-09 88.800003 92.000000 92.099998
2022-06-10 87.699997 97.500000 96.000000
2022-06-11 89.599998 90.400002 96.599998
stimulus 2022-06-09 93.400002 91.900002 92.699997
2022-06-10 92.699997 96.099998 99.099998
2022-06-11 93.099998 97.199997 98.599998
triangle delay 2022-05-31 107.599998 90.099998 91.000000
2022-06-01 110.699997 93.500000 90.300003
2022-06-04 105.300003 94.199997 95.400002
stimulus 2022-05-31 96.800003 90.000000 94.400002
2022-06-01 108.300003 90.699997 90.699997
2022-06-04 114.599998 91.199997 92.099998
Triangle without Perle 2022-06-04 (20 units after the simultaneity filter):
task phase data median a_75_105 b_within_10_of_orth
triangle delay neural 109.099998 False False
triangle delay receptive_fields 90.300003 True True
triangle stimulus neural 99.300003 True True
triangle stimulus receptive_fields 91.699997 True True
Result: partially reproduced (assumed: the code's sigmoid output on the Ring decoders). With the code's decoder (rates in Hz, sigmoid outputs), the Receptive Field data give angles centred near 90 degrees in both datasets and phases (medians 91.8 to 96.7 degrees, within 10 degrees of the Orthogonal data's 90.8 to 95.0), so orthogonality does arise in data without slots, and the Ring neural data do too (88.8 and 93.2). The Triangle neural data do not meet either reading: their medians are 106.7 (delay) and 106.6 degrees (stimulus), outside 75 to 105 and 15.7 and 15.8 degrees from the Orthogonal data's 91.0 and 90.8. Per session the Triangle neural medians are 96.8 to 114.6 degrees, and all three delay-phase sessions are above 105 (105.3 to 110.7); without 2022-06-04 (20 units) the stimulus phase passes (99.3) and the delay phase does not (109.1). The criterion asks for both neural and Receptive Field data near 90 degrees and has no rule for a result that holds in one dataset and not the other; it holds for the Ring and fails for the Triangle neural data, so the claim holds in part, which is this attempt's label.
The Ring pass rests on the code's sigmoid output (models.py), whose range (0, 1) does not cover the Ring targets in [-0.7, 0.7]; the Methods describe linear decoders ("The variables W 0 , B 0 parameterize one linear decoder", Methods, Decoding Analyses). With a linear output the Ring neural medians are 140.2 and 142.1 degrees (variants notebook) and fail both readings, so the neural data would fail in both datasets, which by the criterion is not reproduced (the Triangle decoders were not run with a linear output). The outcome therefore rests on taking the code's decoder and is marked assumed. The variants (z-scored inputs; Ring decoders without the sigmoid) move these medians by up to 51 degrees, so the angles depend on details of the decoder.
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
print(f"Notebook ran in {(time.time() - t_start) / 60:.1f} min")
Notebook ran in 17.2 min
Summary¶
| claim | outcome |
|---|---|
| ext3ab-decoding | partially reproduced |
| ext3c-decoder-orthogonality | partially reproduced (assumed: the code's sigmoid output on the Ring decoders) |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The decoders are the code's models written in numpy and trained with the configs' settings, batched over splits; torch's random streams are not reproduced, so splits, initial weights and batches differ from the paper's. As in the code (dataset.py, whose per-unit normalisation is commented out), the decoders get rates in Hz, and the Ring decoders keep the sigmoid output of models.py. The units entering the decoders are selected with the Methods' AUC p-value (1,000 shuffles), which ranks the same three Triangle and three Ring sessions highest as the code's scripts name; the code's own selection used its resampling p-value. Ten seeds are used, as the legend says, where run_1_triangle.py lists twenty. The synthetic datasets follow dataset.py with numpy's generator. The code's simultaneity filter leaves 20 to 104 units per Triangle session and 44 to 90 per Ring session, below the "~100" simultaneously recorded neurons of the Methods in most sessions.
What was found about the dataset. It holds the trials and spike-sorted units the decoders need; see the main notebook's summary for the metadata findings.