Scherrer, Lynch et al. (2023), dandiset 000691: decoding octant¶
Paper. Scherrer JR, Lynch GF, Zhang JJ, Fee MS (2023). An optical design enabling lightweight and large field-of-view head-mounted microscopes. Nature Methods 20, 546-549. doi:10.1038/s41592-023-01806-1.
Claims. claims/000691/10.1038-s41592-023-01806-1.yaml, written from the full text of the
version of record and the paper's code (GitHub FeeLab/Feescope-Paper-Code, commit d596646). It
decides what is tested and how.
Dandiset. DANDI:000691, version
0.260827.1358, the only published version (the paper cites none); see the main notebook,
reproduce_scherrer_2023.ipynb, for the reasons and the full data check. This notebook holds the
octant decoding claims, apart to keep each notebook within the time limit.
Footprints. The claims use the 2,382 curated footprints, with all 2,948 as a named variant. The dandiset does not mark the curated set, so all 2,948 are used.
| Claim | Data | Status |
|---|---|---|
| ed9b-octant-decoding | dandiset | run under null (n1); (n2) not run |
| ed9c-local-decoding-map | dandiset | not run; cost measured |
More (n2) null decoders are run in decoding_null_n2_quadrant_direction.ipynb and
decoding_null_n2_octant.ipynb, and reduced runs of the map and decoder-size analyses in
decoding_reduced.ipynb; the reduced runs decide no outcome.
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import kiloscope as ks
pd.set_option("display.width", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = ks.asset_urls()
f = ks.open_nwb(urls)
S = ks.load_session(f)
POS = ks.mask_positions(ks.read_masks(urls[ks.NWB_PATH], f))
track, B = ks.load_behaviour(urls, S)
T = S["traces"]
run = B["running"]
X_run = ks.moving_max(T[B["sel"]])[run]
print(f"loaded and tracked in {time.time() - started:.0f} s")
loaded and tracked in 258 s
Data check¶
The full check of the data is in the main notebook. Here: what this notebook decodes.
LABELS = B["octant"][run]
LAMBDA = ks.LAMBDA_OTHER
NAMES = [f"octant {k + 1}" for k in range(8)]
print(f"imaging frames within the tracked video: {len(B['sel'])}; running (speed above its median): {run.sum()}")
print(f"activity: moving maximum over 30 frames of {X_run.shape[1]} footprints; values {X_run.min():.3f} to {X_run.max():.3f}")
print("frames per class:", dict(zip(NAMES, np.bincount(LABELS, minlength=8).tolist())))
n = len(LABELS)
print(f"each fit trains on {round(n * ks.HOLDOUT) - round(n * ks.BUFFER)} frames and tests on {n - round(n * ks.HOLDOUT) - round(n * ks.BUFFER)}; "
f"lambda {LAMBDA}; solver for {X_run.shape[1]} predictors: {ks._solver(X_run, None)}")
imaging frames within the tracked video: 32056; running (speed above its median): 16028
activity: moving maximum over 30 frames of 2948 footprints; values 0.000 to 0.348
frames per class: {'octant 1': 2557, 'octant 2': 2746, 'octant 3': 1858, 'octant 4': 1925, 'octant 5': 2044, 'octant 6': 1836, 'octant 7': 1372, 'octant 8': 1690}
each fit trains on 7854 frames and tests on 7854; lambda 6.75e-08; solver for 2948 predictors: sgd
The labels come from the position tracked in the behaviour video with the code's rule (checked in the main notebook). Every class occurs in the running frames (printed above).
ed9b-octant-decoding¶
Position within the maze can be decoded with octant resolution from the activity, above decoders trained on shuffled data.
Criterion: eight octants (direction not coded), lambda 6.75e-8, the average matrix of 200 fits; the null is 1,000 decoders on their own shuffles, under the readings (n1) one fit per shuffle and (n2) the average of 200 fits per shuffle. Reproduced when the real median accuracy exceeds every null value (p < 0.001) and exceeds chance (1/8) under both; not reproduced when it fails under both; otherwise methods_underspecified.
rng_real = np.random.default_rng(69130)
t0 = time.time()
real = ks.decoding_matrix(X_run, LABELS, 8, LAMBDA, 200, rng_real, n_jobs=4)
sec_real = (time.time() - t0) / 200
print(f"200 fits in {time.time() - t0:.0f} s ({sec_real:.2f} s per fit, 4 in parallel); classes absent from a test part: {real['absent']} of 1600")
print(f"median accuracy (median of the diagonal of the average matrix): {real['median']:.3f}; chance 0.125; ratio {real['median'] / 0.125:.2f}")
print(pd.DataFrame({"class": NAMES, "accuracy": real["diag"].round(3),
"per_fit_5th": np.nanpercentile(real["fit_diag"], 5, axis=0).round(3),
"per_fit_95th": np.nanpercentile(real["fit_diag"], 95, axis=0).round(3)}).to_string(index=False))
t0 = time.time()
X_F = np.asfortranarray(X_run) # columns contiguous, for fast shuffles
null1 = ks.null_medians(X_F, LABELS, 8, LAMBDA, 1000, 1, seed=69131, n_jobs=4)
sec_null = (time.time() - t0) / 1000
print(f"\nnull (n1): 1000 shuffles, one fit each, in {time.time() - t0:.0f} s; NaN medians {np.isnan(null1).sum()}")
print(f" null median accuracies: percentiles 0/50/95/99/100 {np.round(np.nanpercentile(null1, [0, 50, 95, 99, 100]), 3)}")
k = int((null1 >= real["median"]).sum())
print(f" null values at or above the real {real['median']:.3f}: {k}; p = {k / 1000:.4f} (fraction), {(k + 1) / 1001:.4f} ((k + 1) / 1001)")
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
im = ax[0].imshow(real["matrix"], cmap="RdBu_r", vmin=0, vmax=2 * 0.125 + max(0, np.nanmax(real["matrix"]) - 0.25))
ax[0].set_xticks(range(8), NAMES, rotation=45, fontsize=7); ax[0].set_yticks(range(8), NAMES, fontsize=7)
ax[0].set_xlabel("actual"); ax[0].set_ylabel("predicted"); plt.colorbar(im, ax=ax[0])
ax[1].hist(null1, 40, color="0.6", label="null (n1)"); ax[1].axvline(real["median"], color="r", label="decoders")
ax[1].axvline(0.125, color="k", ls=":", label="chance"); ax[1].set_xlabel("median accuracy"); ax[1].legend(fontsize=7)
plt.tight_layout(); plt.show()
print(f"matrix values {np.nanmin(real['matrix']):.3f} to {np.nanmax(real['matrix']):.3f}; null values {np.nanmin(null1):.3f} to {np.nanmax(null1):.3f}, all inside the axes")
200 fits in 116 s (0.58 s per fit, 4 in parallel); classes absent from a test part: 0 of 1600 median accuracy (median of the diagonal of the average matrix): 0.326; chance 0.125; ratio 2.61 class accuracy per_fit_5th per_fit_95th octant 1 0.336 0.078 0.611 octant 2 0.473 0.134 0.785 octant 3 0.353 0.115 0.675 octant 4 0.302 0.049 0.621 octant 5 0.252 0.039 0.607 octant 6 0.362 0.135 0.609 octant 7 0.317 0.072 0.649 octant 8 0.271 0.016 0.560
null (n1): 1000 shuffles, one fit each, in 607 s; NaN medians 0 null median accuracies: percentiles 0/50/95/99/100 [0.004 0.093 0.153 0.182 0.229] null values at or above the real 0.326: 0; p = 0.0000 (fraction), 0.0010 ((k + 1) / 1001)
matrix values 0.018 to 0.473; null values 0.004 to 0.229, all inside the axes
# On equal terms with the (n1) null, which is single fits: the median accuracy of each of the 200 real fits on its
# own (NaN when a class is absent from its test part, as for a null decoder).
fit_med = np.median(real["fit_diag"], axis=1)
ok = ~np.isnan(fit_med)
print(f"real single fits with every class in the test part: {ok.sum()} of 200; their median accuracies, percentiles 0/5/50/95/100 "
f"{np.round(np.percentile(fit_med[ok], [0, 5, 50, 95, 100]), 3)}")
print(f"(n1) null: largest {np.nanmax(null1):.3f}; real single fits at or below it: {(fit_med[ok] <= np.nanmax(null1)).sum()}; "
f"for the weakest real single fit, null values at or above it: {(null1 >= fit_med[ok].min()).sum()} of 1000")
real single fits with every class in the test part: 200 of 200; their median accuracies, percentiles 0/5/50/95/100 [0.179 0.227 0.32 0.379 0.417] (n1) null: largest 0.229; real single fits at or below it: 11; for the weakest real single fit, null values at or above it: 12 of 1000
# Null (n2): each null decoder is the average matrix of 200 fits on its shuffle. One such decoder is run to show
# where these values lie; the 1,000 the criterion asks for are not run (cost below).
t0 = time.time()
null2_sample = []
for s_ in np.random.SeedSequence(69132).spawn(1):
rng_s = np.random.default_rng(s_)
null2_sample.append(ks.decoding_matrix(ks.shuffle_columns(X_F, rng_s), LABELS, 8, LAMBDA, 200, rng_s, n_jobs=4)["median"])
print(f"one (n2) null decoder, 200 fits, in {time.time() - t0:.0f} s: median accuracy {np.round(null2_sample, 3).tolist()} "
f"(the real decoder {real['median']:.3f})")
fits_n2 = 1000 * 200
print(f"(n2) needs {fits_n2} fits; at the {sec_null:.2f} s per fit measured above with 4 fits in parallel, {fits_n2 * sec_null / 3600:.1f} hours "
f"on this machine, where a notebook may take 30 minutes")
one (n2) null decoder, 200 fits, in 119 s: median accuracy [0.11] (the real decoder 0.326) (n2) needs 200000 fits; at the 0.61 s per fit measured above with 4 fits in parallel, 33.7 hours on this machine, where a notebook may take 30 minutes
Result: blocked (not_runnable). Under null reading (n1) the claim holds: the average matrix of
200 fits has a median accuracy of 0.326, 2.61 times chance; none of the 1000 shuffled decoders
reaches it (the largest is 0.229), so p < 0.001, and it exceeds chance. On equal terms with the
single-fit null, the real single fits have median accuracies of 0.179 to 0.417 (median 0.320); 11
of the 200 lie at or below the largest null value, and the weakest has 12 of the 1000 null values
at or above it. Reading (n2), 1000 shuffles of 200 fits each, needs 200,000 fits, 33.7 hours
on this machine at the measured rate; the one (n2) null decoder run here gives 0.110 (eight more are in
decoding_null_n2_octant.ipynb). The criterion asks for both readings, so the claim cannot be
decided within the time limit.
ed9c-local-decoding-map¶
Decoders built from the 32 footprints closest to a point decode octant position significantly in parts of the field of view, with accuracy that varies across it.
Criterion: map points every 5 image pixels over rows 1 to 608 and columns 140 to 700; a point is used when at least 32 footprints lie within 270 um (45 pixels at 5.94 um, or 50 at 5.35); at each, the average matrix of 40 fits on its 32 nearest footprints; the null is 1,000 decoders on shuffled data from 32 random footprints (one fit or 40 each); Benjamini-Hochberg 5%. It depends on ed9b-octant-decoding. Its cost is measured here.
pos = POS["peak"] # (MATLAB row, column), the code's position
rows = []
for pos_label, p in [("code's position", POS["peak"]), ("centroid (variant)", POS["centroid"])]:
for radius in (45, 50):
n_used = 0
for r in range(1, round((608 - 1) / 5) + 1):
for c in range(1, round((700 - 140) / 5) + 1):
coord = (r * 5 + 1, c * 5 + 140)
if (np.hypot(p[:, 0] - coord[0], p[:, 1] - coord[1]) <= radius).sum() >= 32:
n_used += 1
rows.append({"positions": pos_label, "radius_px": radius, "points_used": n_used})
grid = pd.DataFrame(rows)
print(grid.to_string(index=False))
rng_m = np.random.default_rng(691_23)
t0 = time.time()
for _ in range(4):
cols = rng_m.choice(X_run.shape[1], 32, replace=False)
ks.decoding_matrix(X_run[:, cols], LABELS, 8, ks.LAMBDA_OTHER, 1, rng_m)
sec32 = (time.time() - t0) / 4
print(f"one fit on 32 footprints (exact solver, as MATLAB's default for 100 or fewer predictors): {sec32:.2f} s")
for _, g in grid.iterrows():
fits = g.points_used * 40 + 1000 * 40
print(f"{g.positions}, {g.radius_px} px: {g.points_used} points x 40 fits + 1000 null decoders x 40 fits = {fits} fits, "
f"{fits * sec32 / 3600:.1f} hours on one core, {fits * sec32 / 4 / 3600:.1f} on four")
positions radius_px points_used code's position 45 5411 code's position 50 5905 centroid (variant) 45 5409 centroid (variant) 50 5897
one fit on 32 footprints (exact solver, as MATLAB's default for 100 or fewer predictors): 0.66 s code's position, 45 px: 5411 points x 40 fits + 1000 null decoders x 40 fits = 256440 fits, 47.1 hours on one core, 11.8 on four code's position, 50 px: 5905 points x 40 fits + 1000 null decoders x 40 fits = 276200 fits, 50.7 hours on one core, 12.7 on four centroid (variant), 45 px: 5409 points x 40 fits + 1000 null decoders x 40 fits = 256360 fits, 47.1 hours on one core, 11.8 on four centroid (variant), 50 px: 5897 points x 40 fits + 1000 null decoders x 40 fits = 275880 fits, 50.7 hours on one core, 12.7 on four
Result: blocked (not_runnable). It depends on ed9b-octant-decoding, which is blocked. With the
code's positions, 5411 grid points are used at 45 pixels and 5905 at 50 pixels (5409 and 5897 with
centroids); with 40 fits per point and per null decoder that is 256,360 to 276,200 fits, about 11.9
to 12.8 hours on four cores. A reduced run is in decoding_reduced.ipynb.
print(f"total run time {(time.time() - started) / 60:.1f} min")
total run time 18.4 min
Summary¶
| Claim | Outcome (this attempt) | Outcome of -01 |
|---|---|---|
| ed9b-octant-decoding | blocked: not runnable | reproduced |
| ed9c-local-decoding-map | blocked: not runnable | reproduced |
How the analysis differs from the paper's. All 2948 footprints are used where the claims name the 2,382 curated ones, which the dandiset does not mark; position is tracked from the behaviour video with the code's rule; MATLAB's linear SVMs are replaced by scikit-learn's SGD (10 passes, for more than 100 predictors) and LinearSVC (exact, for 100 or fewer), as MATLAB's defaults choose. Only null reading (n1) is run in full; (n2) and the analyses that depend on these claims are costed, not run (more n2 decoders and reduced runs are in the added notebooks). The first attempt (-01) judged ed9b and ed9c against decoders with shifted labels; the claims file now uses the paper's null, shuffled activity, with two readings of a null decoder.