Scherrer, Lynch et al. (2023), dandiset 000691: decoding quadrant and direction

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 quadrant and direction 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
fig2-quadrant-direction-decoding	dandiset	run under null (n1); (n2) not run
fig2k-decoder-size	dandiset	not run; cost measured
fig2l-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 263 s

Data check

The full check of the data is in the main notebook. Here: what this notebook decodes.

LABELS = ks.quad_dir_labels(B)[run]
LAMBDA = ks.LAMBDA_QUAD_DIR
NAMES = [f"{q + 1}{'+' if d else '-'}" for q in range(4) for d in (0, 1)]
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: {'1-': 2845, '1+': 2458, '2-': 1758, '2+': 2025, '3-': 1617, '3+': 2263, '4-': 1474, '4+': 1588}
each fit trains on 7854 frames and tests on 7854; lambda 3.2e-07; 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).

fig2-quadrant-direction-decoding

Linear SVM decoders trained on half of the data decode the mouse's combined quadrant and movement direction on the other half at an accuracy at least three times chance, above that of decoders trained on shuffled data.

Criterion: eight classes (four quadrants by two directions), lambda 3.2e-7, the average matrix of 200 fits. The null is 1,000 decoders on their own shuffles (each footprint's activity circularly shifted by its own offset); p is the fraction of null median accuracies at or above the real one. Null readings: (n1) one fit per shuffle, (n2) the average matrix of 200 fits per shuffle; when they disagree on p < 0.001 the claim is methods_underspecified. Accuracy readings: (a) the median of the diagonal at least 0.375, (b) every diagonal element at least 0.375. Reproduced when p < 0.001 and the accuracy holds under a reading; partially reproduced when p < 0.001 and the median exceeds 1/8 but neither reading holds.

rng_real = np.random.default_rng(69120)
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=69121, 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 115 s (0.57 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.437; chance 0.125; ratio 3.49
class  accuracy  per_fit_5th  per_fit_95th
   1-     0.540        0.157         0.899
   1+     0.587        0.259         0.861
   2-     0.294        0.002         0.643
   2+     0.473        0.153         0.780
   3-     0.239        0.000         0.628
   3+     0.426        0.097         0.785
   4-     0.431        0.044         0.855
   4+     0.443        0.070         0.788

null (n1): 1000 shuffles, one fit each, in 603 s; NaN medians 0
  null median accuracies: percentiles 0/50/95/99/100 [0.    0.074 0.147 0.175 0.227]
  null values at or above the real 0.437: 0; p = 0.0000 (fraction), 0.0010 ((k + 1) / 1001)

matrix values 0.011 to 0.587; null values 0.000 to 0.227, all inside the axes

# The accuracy readings.
print(f"(a) median of the diagonal {real['median']:.3f} >= 0.375: {real['median'] >= 0.375}")
print(f"(b) every diagonal element >= 0.375: {bool((real['diag'] >= 0.375).all())} (smallest {real['diag'].min():.3f}, classes below: "
      f"{[n_ for n_, v in zip(NAMES, real['diag']) if v < 0.375]})")
print(f"median above 1/8: {real['median'] > 0.125}")

(a) median of the diagonal 0.437 >= 0.375: True
(b) every diagonal element >= 0.375: False (smallest 0.239, classes below: ['2-', '3-'])
median above 1/8: True

# 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.138 0.299 0.43  0.55  0.604]
(n1) null: largest 0.227; real single fits at or below it: 2; for the weakest real single fit, null values at or above it: 78 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(69122).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 116 s: median accuracy [0.136] (the real decoder 0.437)
(n2) needs 200000 fits; at the 0.60 s per fit measured above with 4 fits in parallel, 33.5 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.437, 3.49 times chance; none of the 1000 shuffled decoders reaches it (the largest is 0.227), so p < 0.001; and accuracy reading (a) holds (0.437 >= 0.375) while (b) does not (classes 2- and 3- at 0.294 and 0.239). The (n1) null is of single fits and the real matrix an average of 200: on equal terms, the real single fits have median accuracies of 0.138 to 0.604 (median 0.430), 2 of the 200 lie at or below the largest null value, and the weakest has 78 of the 1000 null values at or above it. Reading (n2), 1000 shuffles of 200 fits each, needs 200,000 fits, 33.3 hours on this machine at the measured rate; the one (n2) null decoder run here gives 0.136 (eight more are in decoding_null_n2_quadrant_direction.ipynb). The criterion judges p under both readings, so the claim cannot be decided within the time limit.

fig2k-decoder-size

Decoding accuracy depends on the number of footprints used, with P at or below 2.7 x 10^-4 (Mann-Whitney U) for every decoder size from 2 to 2,048 footprints.

Criterion: for each size 2^1 to 2^11, 64 random subsets; for each, the average matrix of 100 fits and the same subset shuffled, again 100 fits; Mann-Whitney U under two readings. It depends on fig2-quadrant-direction-decoding. Its cost is measured here: one fit at each size, timed.

rng_k = np.random.default_rng(691_22)
rows = []
for size in 2 ** np.arange(1, 12):
    cols = rng_k.choice(X_run.shape[1], size, replace=False)
    t0 = time.time()
    dm = ks.decoding_matrix(X_run[:, cols], LABELS, 8, ks.LAMBDA_OTHER, 1, rng_k)
    rows.append({"size": size, "solver": ks._solver(X_run[:, cols], None), "seconds_per_fit": round(time.time() - t0, 2),
                 "fits_needed": 64 * 2 * 100})
cost = pd.DataFrame(rows)
print(cost.to_string(index=False))
total = (cost.seconds_per_fit * cost.fits_needed).sum()
print(f"total {cost.fits_needed.sum()} fits, {total / 3600:.1f} hours on one core, {total / 4 / 3600:.1f} hours on four")

 size solver  seconds_per_fit  fits_needed
    2  exact             0.03        12800
    4  exact             0.04        12800
    8  exact             0.10        12800
   16  exact             0.41        12800
   32  exact             0.63        12800
   64  exact             0.89        12800
  128    sgd             0.16        12800
  256    sgd             0.25        12800
  512    sgd             0.44        12800
 1024    sgd             0.88        12800
 2048    sgd             1.70        12800
total 140800 fits, 19.7 hours on one core, 4.9 hours on four

Result: blocked (not_runnable). It depends on fig2-quadrant-direction-decoding, which is blocked, and it needs 140,800 fits, measured here at 4.8 hours on four cores. A reduced run is in decoding_reduced.ipynb.

fig2l-local-decoding-map

Decoders built from the 32 footprints closest to a point decode quadrant and direction with accuracies that differ substantially across the field of view.

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 10 fits on its 32 nearest footprints; the null is 1,000 decoders on shuffled data from 32 random footprints (one fit or 10 each); Benjamini-Hochberg 5%. It depends on fig2-quadrant-direction-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 * 10 + 1000 * 10
    print(f"{g.positions}, {g.radius_px} px: {g.points_used} points x 10 fits + 1000 null decoders x 10 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): 1.01 s
code's position, 45 px: 5411 points x 10 fits + 1000 null decoders x 10 fits = 64110 fits, 18.1 hours on one core, 4.5 on four
code's position, 50 px: 5905 points x 10 fits + 1000 null decoders x 10 fits = 69050 fits, 19.4 hours on one core, 4.9 on four
centroid (variant), 45 px: 5409 points x 10 fits + 1000 null decoders x 10 fits = 64090 fits, 18.1 hours on one core, 4.5 on four
centroid (variant), 50 px: 5897 points x 10 fits + 1000 null decoders x 10 fits = 68970 fits, 19.4 hours on one core, 4.9 on four

Result: blocked (not_runnable). It depends on fig2-quadrant-direction-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 10 fits per point and per null decoder that is 64,090 to 69,050 fits, 4.4 to 4.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.5 min

Summary

Claim	Outcome (this attempt)	Outcome of -01
fig2-quadrant-direction-decoding	blocked: not runnable	(dropped in -01)
fig2k-decoder-size	blocked: not runnable	(new)
fig2l-local-decoding-map	blocked: not runnable	(new)

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) did not test fig2-quadrant-direction-decoding: its claims file had dropped the claim, which the full text states.