Reproducing Scherrer, Lynch et al. (2023) from dandiset 000691¶
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. The claims file claims/000691/10.1038-s41592-023-01806-1.yaml was written on
2026-10-20 from the full text of the version of record (main text, figure legends, Methods,
Extended Data legends and reporting summary) and the paper's code, GitHub
FeeLab/Feescope-Paper-Code at commit d596646, by a reader who did not see the data, and checked
by a second reader. It decides what is tested and how. This is the second attempt on the dandiset
(000691-scherrer-2023-02); the first (-01) worked from claims written from the abstract and
the public legends, and is deprecated.
Dandiset. DANDI:000691, version 0.260827.1358. Why this version: the paper cites no DANDI version (its data statement points to an Academic Torrents deposit), and 0.260827.1358 is the only published version, the one the first attempt used. The draft was last modified when that version was published and holds the same two assets with the same blobs (printed in the data check), so the data have not changed since the first attempt.
What the dandiset holds. One session of one mouse (4 mm cranial window, jGCaMP7f) imaged with the Kiloscope while it explores a square loop maze: the imaging video (consistent with the raw video, see the motion notebook), 2948 EXTRACT footprints with traces and masks, and a behaviour video (a second asset, an AVI). It holds no position of the mouse, no curation of the footprints, no pixel size, and none of the data of the other three mice, the three birds, or the optical tests.
Outside inputs. The claims file places the optical measurements of five claims (bead images,
their region lists, and images of a thin layer of dye) in the paper's code repository at commit
d596646 (data: elsewhere); they are read from there. Every other result is computed from the
dandiset. Method details come from the claims file, which cites the Methods and the code; see the
helper module kiloscope.py.
Notebooks of this entry. This notebook holds the data check and every claim but four. The
decoding claims are in decoding_quadrant_direction.ipynb (fig2-quadrant-direction-decoding,
fig2k-decoder-size, fig2l-local-decoding-map) and decoding_octant.ipynb (ed9b-octant-decoding,
ed9c-local-decoding-map), and the motion claim in motion_scherrer_2023.ipynb
(ed7-motion-small), each to keep within the time limit. Outputs added after the first
verification are in decoding_null_n2_quadrant_direction.ipynb and decoding_null_n2_octant.ipynb
(more null decoders of reading n2) and decoding_reduced.ipynb (reduced runs of fig2k, fig2l and
ed9c, which decide no outcome).
Random numbers are drawn from a generator per section, each with its own fixed seed.
| Claim | Data | Notebook | Status |
|---|---|---|---|
| abstract-thousands-of-neurons | dandiset | this one | tested |
| ed6a-curated-footprints | dandiset | this one | blocked: the curation is not marked |
| ed6bcd-other-mice | dandiset | this one | blocked: data missing |
| ed5-songbirds | dandiset | this one | blocked: data missing |
| birds-singing-rate | dandiset | this one | blocked: data missing |
| fig2d-lman-sparse-bursting | dandiset | this one | blocked: no method fixes it (its data are also missing) |
| abstract-kiloscope-fov | dandiset | this one | tested |
| abstract-weights | dandiset | this one | blocked: data missing |
| fig1e-featherscope-psf | elsewhere (GitHub) | this one | tested |
| fig1f-kiloscope-psf | elsewhere (GitHub) | this one | tested |
| ed1cf-field-curvature | elsewhere (GitHub) | this one | blocked: no method fixes it; shown |
| ed1-ed4-resolution | dandiset | this one | blocked: data missing |
| ed2b-featherscope-flat-illumination | elsewhere (GitHub) | this one | tested |
| ed2c-kiloscope-flat-illumination | elsewhere (GitHub) | this one | tested |
| ed3-adc | dandiset | this one | blocked: data missing |
| ed7-motion-small | dandiset | motion | tested |
| fig2-quadrant-direction-decoding | dandiset | decoding (quadrant) | see that notebook |
| fig2k-decoder-size | dandiset | decoding (quadrant) | see that notebook |
| fig2l-local-decoding-map | dandiset | decoding (quadrant) | see that notebook |
| ed8-coactive-distributed | dandiset | this one | tested |
| ed9b-octant-decoding | dandiset | decoding (octant) | see that notebook |
| ed9c-local-decoding-map | dandiset | decoding (octant) | see that notebook |
| ed10-localized-activity | dandiset | this one | blocked: no method fixes it; shown |
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from dandi.dandiapi import DandiAPIClient
from scipy import ndimage
from scipy.optimize import linear_sum_assignment
import kiloscope as ks
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
urls = ks.asset_urls()
f = ks.open_nwb(urls)
S = ks.load_session(f)
M = ks.read_masks(urls[ks.NWB_PATH], f)
POS = ks.mask_positions(M)
print(f"NWB traces, timestamps and masks loaded in {time.time() - started:.0f} s")
t0 = time.time()
# Frames kept for display: four in the opening part the code leaves out, four after it.
SAMPLE_FRAMES = (0, 640, 1280, 1900, 4000, 8000, 12000, 15000)
track, B = ks.load_behaviour(urls, S, keep_frames=SAMPLE_FRAMES)
B_code = ks.behaviour(S["t_img"], S["t_video"], track, unwrap="code")
print(f"behaviour video decoded and tracked in {time.time() - t0:.0f} s")
NWB traces, timestamps and masks loaded in 31 s
behaviour video decoded and tracked in 157 s
Data check¶
Each check of the metadata against the values is stated under its output. The last part of the
check gives, for each finding of the first attempt (000691-scherrer-2023-01), whether it is fixed.
with DandiAPIClient() as client:
dandiset = client.get_dandiset(ks.DANDISET_ID, ks.DANDISET_VERSION)
raw = dandiset.get_raw_metadata()
versions = [(v.identifier, str(v.created)[:19], str(v.modified)[:19])
for v in client.get_dandiset(ks.DANDISET_ID).get_versions()]
draft = client.get_dandiset(ks.DANDISET_ID, "draft")
assets = {v: {a.path: (a.size, a.blob) for a in client.get_dandiset(ks.DANDISET_ID, v).get_assets()}
for v in (ks.DANDISET_VERSION, "draft")}
print("Versions (identifier, created, modified):")
for v in versions:
print(" ", v)
print("Draft last modified:", str(draft.version.modified)[:19])
summary = raw["assetsSummary"]
print("Assets of the pinned version (size in bytes, blob id):")
for p in urls:
print(" ", p, assets[ks.DANDISET_VERSION][p])
print("Draft holds the same assets with the same blobs:", assets["draft"] == assets[ks.DANDISET_VERSION])
print("Dandiset summary: subjects", summary.get("numberOfSubjects"), "; species",
[s["name"] for s in summary.get("species", [])])
print(" measurement techniques:", [m["name"] for m in summary.get("measurementTechnique", [])])
print(" variables measured:", summary.get("variableMeasured"))
print("Dandiset description says the data come from:",
raw["description"].split("This dataset includes")[1].split(".")[0].strip())
print()
m = S["meta"]
for k in ["session_description", "experiment_description", "device_description", "imaging_series",
"imaging_description", "video_description", "video_external_file", "indicator", "location",
"excitation_lambda", "imaging_plane_fields", "grid_spacing_present", "virus", "surgery", "subject",
"trace_description", "trace_unit", "segmentation_type", "cellfind_min_snr",
"plane_segmentation_columns", "movie_shape", "movie_dtype", "movie_chunks", "mask_shape"]:
print(f"{k}: {m[k]}")
print("processing modules:", list(f["processing"].keys()), "; acquisition:", list(f["acquisition"].keys()))
print("intervals or epochs in the file:", [k for k in f.keys() if k in ("intervals", "epochs")])
Versions (identifier, created, modified):
('draft', '2023-10-30 15:17:41', '2026-08-27 13:58:02')
('0.260827.1358', '2026-08-27 13:58:02', '2026-08-27 13:58:03')
Draft last modified: 2026-08-27 13:58:02
Assets of the pinned version (size in bytes, blob id):
sub-4mm-mouse-1/sub-4mm-mouse-1_image+ophys/60d277d6-24fd-4817-bb02-41382db59172_external_file_0.avi (105329192, '1da756f3-ab4d-4796-bf27-ad8f9054eaf0')
sub-4mm-mouse-1/sub-4mm-mouse-1_image+ophys.nwb (15232760735, 'c3d0b627-6abb-4b8e-a44b-e51fe53cccfc')
Draft holds the same assets with the same blobs: True
Dandiset summary: subjects 1 ; species ['Mus musculus - House mouse']
measurement techniques: ['analytical technique', 'two-photon microscopy technique', 'surgical technique']
variables measured: ['TwoPhotonSeries', 'ImagingPlane', 'ProcessingModule', 'PlaneSegmentation', 'OpticalChannel']
Dandiset description says the data come from: the calcium imaging data from our microscope optimized for weight
session_description: This session includes calcium imaging recorded from a head-mounted microscope in a freely moving mouse while simultaneously recording more than a thousand neurons in cortex.
experiment_description: Here we present a fluorescence microscope light path that enables imaging, during free behavior, of thousands of neurons in mice and hundreds of neurons in juvenile songbirds. The light path eliminates traditional illumination optics, allowing for head-mounted microscopes that have both a lower weight and a larger field of view (FOV) than previously possible. Using this light path, we designed two microscopes: one optimized for FOV (~4 mm FOV; 1.4 g), and the other optimized for weight (1.0 mm FOV; 1.0 g). This dataset includes the calcium imaging data from our microscope optimized for weight.For this experiment, we stereotactically located a region of the mouse brain at coordinates associatied with primary visual and somatosensory regions, performed approximately 20 viral injections of GCaMP7f across a 4-mm-diameter region of cortex, and then implanted a 4-mm-diameter glass window in the skull. After letting each mouse recover, we installed our head-mounted microscope and allowed it to freely explore a circular maze while recording calcium activity for several minutes.
device_description: Head-mounted microscope (Kiloscope) optimized for weight (1.0 mm FOV; 1.0 g).
imaging_series: TwoPhotonSeries
imaging_description: Calcium imaging data at ~30 fps in grayscale format
video_description: Behavior video of animal moving in environment at ~30 fps
video_external_file: sub-4mm-mouse-1_image+ophys/60d277d6-24fd-4817-bb02-41382db59172_external_file_0.avi
indicator: unknown
location: unknown
excitation_lambda: nan
imaging_plane_fields: ['OpticalChannel', 'channel_0', 'description', 'device', 'excitation_lambda', 'indicator', 'location']
grid_spacing_present: False
virus: GCaMP7f
surgery: 4mm glass window implant over right visual and somatosensory cortex
subject: {'age': 'P3M', 'description': '~20 virus injections in layer 2/3 of cortex', 'genotype': 'C57/B6', 'sex': 'M', 'species': 'Mus musculus', 'subject_id': '4mm-mouse-1'}
trace_description: Array of df/F traces.
trace_unit: n.a.
segmentation_type: EXTRACTSegmentation
cellfind_min_snr: 2.0
plane_segmentation_columns: ['image_mask', 'ROICentroids']
movie_shape: (36596, 808, 608)
movie_dtype: uint16
movie_chunks: (1200, 20, 20)
mask_shape: (2948, 808, 608)
processing modules: ['ophys'] ; acquisition: ['TwoPhotonSeries', 'Video: home_arena_2021-07-26T13_50_50']
intervals or epochs in the file: []
The dandiset has one published version, 0.260827.1358, and the draft holds the same assets with the
same blobs, so nothing has changed since the first attempt. The device description names the
large-FOV microscope (Kiloscope) but gives the size and weight of the other one ("optimized for
weight (1.0 mm FOV; 1.0 g)"), and the dandiset description also says the data come from "our
microscope optimized for weight". The data are from the Kiloscope: a 4 mm window (the surgery
field), thousands of footprints, and frames of 808 x 608 pixels, the size the paper's code uses for
the Kiloscope (sf3/widefield/analysis.m). The imaging is stored as a TwoPhotonSeries and the
dandiset lists two-photon microscopy, but the microscope is a one-photon head-mounted
epifluorescence microscope. The imaging plane gives the indicator and location as "unknown", an
excitation wavelength of NaN and no grid spacing (pixel size), although general/virus says
GCaMP7f. The plane segmentation has only image_mask and ROICentroids: no column marks the
footprints kept by the paper's curation. The only processing module is ophys: there is no
position of the mouse, and no interval marks the opening of the behaviour video. The
segmentation parameters record cellfind_min_snr = 2, as the Methods state.
def rate_summary(t):
dt = np.diff(t)
return {"n": len(t), "first_s": round(t[0], 3), "last_s": round(t[-1], 3),
"median_dt_ms": round(1000 * np.median(dt), 2), "rate_hz": round(1 / np.median(dt), 2),
"min_dt_ms": round(1000 * dt.min(), 2), "max_dt_ms": round(1000 * dt.max(), 2),
"within_10pct": round(float(np.mean(np.abs(dt / np.median(dt) - 1) < 0.1)), 3)}
print(pd.DataFrame({"imaging": rate_summary(S["t_img"]), "traces": rate_summary(S["t_trace"]),
"behaviour video": rate_summary(S["t_video"])}).to_string())
print("Trace timestamps equal the imaging timestamps:", bool(np.array_equal(S["t_img"], S["t_trace"])))
print(f"Imaging frames: {m['movie_shape'][0]} in the movie, {len(S['t_img'])} timestamps, {S['traces'].shape[0]} trace samples")
print(f"\nAVI: {track['n_frames']} frames of {track['frame_shape']}, nominal {track['nominal_rate']:.0f} frames/s, "
f"so {track['n_frames'] / track['nominal_rate']:.0f} s at that rate; its NWB timestamps span "
f"{S['t_video'][-1] - S['t_video'][0]:.0f} s")
print("AVI frames equal NWB video timestamps:", track["n_frames"] == len(S["t_video"]))
T = S["traces"]
print(f"\nTraces: {T.shape[0]} frames x {T.shape[1]} footprints; values {T.min():.3f} to {T.max():.3f}; "
f"fraction exactly zero {np.mean(T == 0):.3f}; footprints never above zero {(T.max(0) == 0).sum()}")
print("Per-footprint maximum, percentiles 0/5/50/95/100:", np.round(np.percentile(T.max(0), [0, 5, 50, 95, 100]), 4))
print("ROI ids are 0..n-1 and the trace table references them in order:",
bool(np.array_equal(S["ids"], np.arange(T.shape[1])) and np.array_equal(S["rois"], S["ids"])))
c = S["centroids"]
print(f"ROICentroids: first column {c[:, 0].min()} to {c[:, 0].max()} (movie axis 1 has {m['movie_shape'][1]} pixels), "
f"second column {c[:, 1].min()} to {c[:, 1].max()} (movie axis 2 has {m['movie_shape'][2]})")
print(f"\nMasks: {M['n_chunks']} chunks, {M['compressed_bytes'] / 1e6:.1f} MB compressed; non-zero pixels per footprint, "
f"percentiles 0/50/100: {np.percentile(POS['n_pixels'], [0, 50, 100])}")
print(f" largest weight per footprint {POS['max_weight'].min()} to {POS['max_weight'].max()}; footprints with no weight "
f"equal to 1: {int(np.isnan(POS['peak'][:, 0]).sum())}")
pk = POS["peak"]
d_row, d_col = pk[:, 0] - (c[:, 1] + 1), pk[:, 1] - (c[:, 0] + 1)
print(f" the code's position (first pixel of weight 1) against ROICentroids (axis 2 + 1, axis 1 + 1): equal for "
f"{np.mean((d_row == 0) & (d_col == 0)):.3f} of footprints, largest difference {max(abs(d_row).max(), abs(d_col).max()):.0f} px")
cw = POS["centroid"]
print(f" ROICentroids equal the rounded weighted centroid (axis 1, axis 2) for "
f"{np.mean((np.round(cw[:, 1] - 1) == c[:, 0]) & (np.round(cw[:, 0] - 1) == c[:, 1])):.3f} of footprints")
print(f" weighted centroid against the code's position: median distance {np.median(np.hypot(*(cw - pk).T)):.2f} px, "
f"largest {np.hypot(*(cw - pk).T).max():.2f} px")
print(f" positions span rows {np.nanmin(pk[:, 0]):.0f} to {np.nanmax(pk[:, 0]):.0f} and columns {np.nanmin(pk[:, 1]):.0f} "
f"to {np.nanmax(pk[:, 1]):.0f} (MATLAB orientation: 608 rows, 808 columns)")
imaging traces behaviour video n 36596.000 36596.000 15521.000 first_s 0.000 0.000 0.259 last_s 1234.515 1234.515 1234.473 median_dt_ms 33.720 33.720 79.910 rate_hz 29.660 29.660 12.510 min_dt_ms 7.830 7.830 42.050 max_dt_ms 62.460 62.460 111.680 within_10pct 0.944 0.944 0.926 Trace timestamps equal the imaging timestamps: True Imaging frames: 36596 in the movie, 36596 timestamps, 36596 trace samples AVI: 15521 frames of (480, 640, 3), nominal 30 frames/s, so 517 s at that rate; its NWB timestamps span 1234 s AVI frames equal NWB video timestamps: True Traces: 36596 frames x 2948 footprints; values 0.000 to 0.348; fraction exactly zero 0.752; footprints never above zero 0 Per-footprint maximum, percentiles 0/5/50/95/100: [0.0142 0.0307 0.0639 0.157 0.3481] ROI ids are 0..n-1 and the trace table references them in order: True ROICentroids: first column 164 to 685 (movie axis 1 has 808 pixels), second column 30 to 587 (movie axis 2 has 608) Masks: 32768 chunks, 9.2 MB compressed; non-zero pixels per footprint, percentiles 0/50/100: [ 4. 32. 164.] largest weight per footprint 1.0 to 1.0; footprints with no weight equal to 1: 0 the code's position (first pixel of weight 1) against ROICentroids (axis 2 + 1, axis 1 + 1): equal for 0.844 of footprints, largest difference 1 px ROICentroids equal the rounded weighted centroid (axis 1, axis 2) for 0.425 of footprints weighted centroid against the code's position: median distance 0.73 px, largest 3.27 px positions span rows 31 to 588 and columns 165 to 686 (MATLAB orientation: 608 rows, 808 columns)
The imaging runs at a median 29.7 frames/s over 1234 s, as the series description says (about 30
frames/s), with irregular intervals. The behaviour video's timestamps run at a median of 12.5
frames/s over the same 1234 s, while its description says about 30 frames/s and the AVI header says
30 frames/s. The AVI has as many frames as the NWB file has timestamps for it, so the timestamps
are taken as the true frame times: at 30 frames/s the AVI would last only 517 s. The traces are
stored as "df/F" with unit "n.a."; three quarters of their values are exactly zero, as EXTRACT's
non-negative output gives. Every mask has a largest weight of exactly 1, so the code's position (the
first pixel of weight 1) exists for every footprint; it equals the stored ROICentroids for 84% of
footprints and is within one pixel for the rest; the stored values are described as "the x, y, (z) centroids of
each ROI" but equal the rounded weighted centroid for only 43% of footprints (printed above), and
how they were computed is not said. The masks are stored transposed relative to the code's arrays (axis 1 is the
code's column, of 808).
# Which movie axis each centroid column indexes, and whether traces, masks and movie share a clock.
movie = f[ks.MOVIE]
rng_check = np.random.default_rng(691_00)
active = np.argsort(T[:2400].std(0))[::-1]
picked = active[:6]
controls = rng_check.choice(active[200:], size=6, replace=False)
n_check = 2400
lags = np.arange(-15, 16)
def hp(v):
return v - pd.Series(v).rolling(301, center=True, min_periods=1).median().to_numpy()
def xcorr(a, b):
out = []
for L in lags:
if L >= 0:
out.append(np.corrcoef(a[L:], b[:len(b) - L])[0, 1])
else:
out.append(np.corrcoef(a[:L], b[-L:])[0, 1])
return np.array(out)
rows = []
for i, j in zip(picked, controls):
c0, c1 = S["centroids"][i]
patch = hp(movie[:n_check, c0 - 2:c0 + 3, c1 - 2:c1 + 3].astype(float).mean((1, 2)))
own, other = hp(T[:n_check, i].astype(float)), hp(T[:n_check, j].astype(float))
xo, xc = xcorr(patch, own), xcorr(patch, other)
rows.append({"footprint": i, "r_own_lag0": round(xo[15], 3), "best_lag_frames": int(lags[np.nanargmax(xo)]),
"r_own_best": round(np.nanmax(xo), 3), "r_other_footprint_lag0": round(xc[15], 3)})
print("Movie pixels (5 x 5 at each footprint's ROICentroid, first 2400 frames, slow drift removed) against traces:")
print(pd.DataFrame(rows).to_string(index=False))
Movie pixels (5 x 5 at each footprint's ROICentroid, first 2400 frames, slow drift removed) against traces:
footprint r_own_lag0 best_lag_frames r_own_best r_other_footprint_lag0
188 0.601 -1 0.602 0.052
69 0.696 1 0.696 -0.061
577 0.749 2 0.753 0.072
86 0.300 1 0.302 0.043
17 0.579 -9 0.611 -0.150
1627 0.869 1 0.869 -0.074
The movie pixels at a footprint's stored position follow its own trace at lag 0 and do not follow another footprint's trace (printed above), so the traces come from this movie on its clock and the positions index the movie as (axis 1, axis 2). These are the six footprints with the most variable traces in the first 2400 frames (80 s); the movie is chunked 1200 frames at a time, and the whole movie is read in the motion notebook.
x, y, area = track["x"], track["y"], track["area"]
tracked = np.arange(len(x)) >= ks.FIRST_VIDEO_FRAME
# The rotated crop against a full rotation of one frame (scipy.ndimage.rotate, bilinear).
rgb0, crop0 = track["frames"][4000]
full = ndimage.rotate(rgb0.astype(np.float64), ks.ROTATION_DEG, axes=(1, 0), reshape=True, order=1)
r0, r1 = ks.CROP_ROWS
c0_, c1_ = ks.CROP_COLS
print(f"Rotated frame {full.shape[:2]}; the crop against a full rotation of frame 4000: largest difference "
f"{np.abs(np.round(full[r0 - 1:r1, c0_ - 1:c1_]) - crop0).max():.0f} grey level")
print(f"Frames before frame {ks.FIRST_VIDEO_FRAME} (the code's 64 s at a nominal 30 frames/s), left out: "
f"{(~tracked).sum()}, up to t = {S['t_video'][ks.FIRST_VIDEO_FRAME]:.1f} s")
print("Frames from then on with no dark blob:", int(np.isnan(x[tracked]).sum()))
print("Blob area (pixels) from then on, percentiles 0/1/50/99/100:", np.percentile(area[tracked], [0, 1, 50, 99, 100]))
print(f"\nImaging frames within the tracked video: {len(B['sel'])} of {len(S['t_img'])}, t = {B['t'][0]:.1f} to {B['t'][-1]:.1f} s")
xmin, xmax, ymin, ymax = ks.RECT
dist_side = np.stack([abs(B["x"] - xmin), abs(B["y"] - ymin), abs(B["x"] - xmax), abs(B["y"] - ymax)], 1).min(1)
print(f"Distance from the nearest side of the code's rectangle (x {xmin} to {xmax}, y {ymin} to {ymax}), pixels, "
f"percentiles 50/95/99/100: {np.round(np.percentile(dist_side, [50, 95, 99, 100]), 1)}")
print("Fraction of imaging frames per quadrant (0 to 3):", np.round(np.bincount(B["quadrant"], minlength=4) / len(B["quadrant"]), 3))
print("Fraction per octant:", np.round(np.bincount(B["octant"], minlength=8) / len(B["octant"]), 3))
step = np.diff(B["sunwrap"])
print(f"Distance travelled: {np.abs(step).sum() / B['perimeter']:.1f} loop lengths ({np.abs(step).sum():.1f} m); "
f"largest step between imaging frames {np.abs(step).max() * 100:.1f} cm (loop {B['perimeter'] * 100:.1f} cm)")
vstep = np.hypot(np.diff(x[tracked]), np.diff(y[tracked]))
print("\nTracking quality, from frame 1920 on:")
print(" step between consecutive video frames (pixels), percentiles 50/90/99/99.9/100:",
np.round(np.percentile(vstep, [50, 90, 99, 99.9, 100]), 1))
print(" steps above 20, 40, 80 pixels:", [int((vstep > k).sum()) for k in (20, 40, 80)])
for j in np.flatnonzero(vstep > 40):
i = ks.FIRST_VIDEO_FRAME + j
print(f" frame {i} to {i + 1} (t = {S['t_video'][i]:.1f} s): {vstep[j]:.0f} px; blob areas {area[i]:.0f}, {area[i + 1]:.0f}")
small = (area < 50) & tracked
print(f" frames with blob area below 50 pixels: {small.sum()} of {tracked.sum()}")
bad = small.copy()
jump = np.flatnonzero(vstep > 40) + ks.FIRST_VIDEO_FRAME
bad[jump] = bad[jump + 1] = True
B_clean = ks.behaviour(S["t_img"], S["t_video"], dict(track, x=np.where(bad, np.nan, x), y=np.where(bad, np.nan, y)))
same_frames = np.array_equal(B_clean["sel"], B["sel"])
print(f" without those {bad.sum()} video frames: same imaging frames {same_frames}; imaging frames whose quadrant changes "
f"{np.mean(B_clean['quadrant'] != B['quadrant']):.4f}, octant {np.mean(B_clean['octant'] != B['octant']):.4f}, "
f"running {np.mean(B_clean['running'] != B['running']):.4f}, direction {np.mean(B_clean['direction'] != B['direction']):.4f}")
print(f"Running frames (speed above median): {B['running'].sum()}; direction +1 in {np.mean(B['direction'][B['running']] > 0):.3f} of them")
laps = np.flatnonzero(np.abs(np.diff(B["s"])) > B["perimeter"] / 2)
resid = np.abs(np.diff(B_code["sunwrap"]))[laps] * 100
print(f"\nThe code's unwrap (MATLAB unwrap of s * 100 with tolerance max(s) * 100 / 2): at the {len(laps)} wraps of linear "
f"position it leaves steps of {resid.min():.2f} to {resid.max():.2f} cm (median {np.median(resid):.2f}); the loop unwrapped "
f"leaves {np.abs(step[laps]).min() * 100:.2f} to {np.abs(step[laps]).max() * 100:.2f} cm; median step between frames "
f"{np.median(np.abs(step)) * 100:.2f} cm")
print(f" labels that differ with the code's unwrap: running {np.mean(B_code['running'] != B['running']):.4f}, "
f"direction {np.mean(B_code['direction'] != B['direction']):.4f}, running and direction (8 classes) among frames "
f"running in both {np.mean((ks.quad_dir_labels(B_code) != ks.quad_dir_labels(B))[B['running'] & B_code['running']]):.4f}")
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
for q in range(4):
k = B["quadrant"] == q
ax[0].plot(B["x"][k], B["y"][k], ",", alpha=0.3, label=f"quadrant {q}")
ax[0].plot([xmin, xmax, xmax, xmin, xmin], [ymin, ymin, ymax, ymax, ymin], "k", lw=0.8)
ax[0].invert_yaxis(); ax[0].set_aspect("equal"); ax[0].set_title("position at imaging frames (pixels of the crop)")
ax[0].legend(markerscale=20, fontsize=7)
ax[1].plot(B["t"], B["s"], lw=0.3); ax[1].set_xlabel("time (s)"); ax[1].set_ylabel("linear position (m)")
plt.tight_layout(); plt.show()
print(f"Plotted: x {B['x'].min():.0f} to {B['x'].max():.0f}, y {B['y'].min():.0f} to {B['y'].max():.0f}, "
f"linear position {B['s'].min():.2f} to {B['s'].max():.2f} m; nothing outside the axes")
Rotated frame (666, 515); the crop against a full rotation of frame 4000: largest difference 0 grey level Frames before frame 1920 (the code's 64 s at a nominal 30 frames/s), left out: 1920, up to t = 153.1 s Frames from then on with no dark blob: 0 Blob area (pixels) from then on, percentiles 0/1/50/99/100: [ 4. 68. 278. 405. 462.] Imaging frames within the tracked video: 32056 of 36596, t = 153.1 to 1234.4 s Distance from the nearest side of the code's rectangle (x 44 to 228, y 34 to 270), pixels, percentiles 50/95/99/100: [ 3.9 10.3 14.3 25.3] Fraction of imaging frames per quadrant (0 to 3): [0.312 0.229 0.251 0.208] Fraction per octant: [0.121 0.191 0.105 0.124 0.135 0.116 0.099 0.109] Distance travelled: 47.7 loop lengths (70.3 m); largest step between imaging frames 7.3 cm (loop 147.1 cm) Tracking quality, from frame 1920 on: step between consecutive video frames (pixels), percentiles 50/90/99/99.9/100: [ 2.1 7.6 13.2 18. 21.7] steps above 20, 40, 80 pixels: [2, 0, 0] frames with blob area below 50 pixels: 69 of 13601 without those 69 video frames: same imaging frames True; imaging frames whose quadrant changes 0.0003, octant 0.0003, running 0.0104, direction 0.0021 Running frames (speed above median): 16028; direction +1 in 0.520 of them The code's unwrap (MATLAB unwrap of s * 100 with tolerance max(s) * 100 / 2): at the 98 wraps of linear position it leaves steps of 0.01 to 3.08 cm (median 1.50); the loop unwrapped leaves 0.01 to 4.70 cm; median step between frames 0.12 cm labels that differ with the code's unwrap: running 0.0086, direction 0.0068, running and direction (8 classes) among frames running in both 0.0057
Plotted: x 29 to 241, y 24 to 287, linear position -0.03 to 1.50 m; nothing outside the axes
# Sample frames of the behaviour video: the code's rotated crop, with the tracked centroid (red cross).
fig, ax = plt.subplots(2, 4, figsize=(13, 7))
rows = []
for a_, i in zip(ax.ravel(), SAMPLE_FRAMES):
crop = track["frames"][i][1]
cx, cy, ca = ks._largest_dark_blob(crop)
a_.imshow(crop)
a_.plot([xmin, xmax, xmax, xmin, xmin], [ymin, ymin, ymax, ymax, ymin], color="w", lw=0.6)
if np.isfinite(cx):
a_.plot(cx - 1, cy - 1, "r+", ms=14, mew=2)
a_.set_title(f"frame {i}, t = {S['t_video'][i]:.0f} s{' (left out)' if i < ks.FIRST_VIDEO_FRAME else ''}", fontsize=8)
a_.axis("off")
rows.append({"frame": i, "t_s": round(S["t_video"][i], 1), "left_out": i < ks.FIRST_VIDEO_FRAME,
"blob_area_px": ca, "x": round(cx, 1), "y": round(cy, 1)})
plt.tight_layout(); plt.show()
print(pd.DataFrame(rows).to_string(index=False))
frame t_s left_out blob_area_px x y
0 0.3 True 294 227.4 240.7
640 51.1 True 193 203.1 40.3
1280 102.0 True 329 141.8 35.1
1900 151.5 True 302 228.9 97.7
4000 318.4 False 279 163.1 274.5
8000 636.3 False 317 90.6 42.2
12000 954.4 False 317 46.2 78.9
15000 1193.1 False 308 231.2 194.4
The dandiset gives no position, so it is tracked here with the code's rule: frames from 1920 on, rotated by 86.79 degrees and cropped as the code does (the crop equals a full rotation to within one grey level, printed above), and the code's colour rule. The track lies along the code's rectangle (distances printed above), so the rectangle's fixed corners fit this video. The mouse runs about 48 loop lengths and visits every quadrant and octant. No step between consecutive video frames exceeds 40 pixels, and leaving out the 69 frames with a blob under 50 pixels changes the quadrant of 0.03% of imaging frames. The sample frames show the crop and the tracked centroid on the mouse. The 64 s the code skips, 1920 frames at the AVI's nominal rate, are 153 s on the timestamps. The code unwraps linear position with MATLAB's unwrap on centimetres, which corrects each lap's jump by a multiple of 2 pi cm; as the loop (147.1 cm) is not such a multiple, the position gains an error of up to pi cm at every lap. At the 98 wraps its steps are 0.01 to 3.08 cm where the position unwrapped by the loop length steps 0.01 to 4.70 cm (printed above). The claims file says position is unwrapped, which is done here by the loop length; with the code's version, about 1% of the running and direction labels differ (printed above).
# Shared clock of behaviour and imaging: decode the quadrant (4 classes) with labels shifted in time.
rng_clock = np.random.default_rng(691_01)
rate = 1 / np.median(np.diff(S["t_img"]))
act = ks.moving_max(T[B["sel"]])
run = B["running"]
rows = []
for lag_s in [-60, -20, -5, -2, -1, 0, 1, 2, 5, 20, 60]:
shifted = np.roll(B["quadrant"], int(round(lag_s * rate)))
dm = ks.decoding_matrix(act[run], shifted[run], 4, ks.LAMBDA_QUAD_DIR, 3, rng_clock, n_jobs=4)
rows.append({"label_shift_s": lag_s, "median_accuracy_3_fits": round(dm["median"], 3)})
print("Quadrant decoding (chance 0.25), the average matrix of 3 fits per shift:")
print(pd.DataFrame(rows).to_string(index=False))
Quadrant decoding (chance 0.25), the average matrix of 3 fits per shift:
label_shift_s median_accuracy_3_fits
-60 0.293
-20 0.175
-5 0.286
-2 0.403
-1 0.457
0 0.495
1 0.493
2 0.519
5 0.420
20 0.275
60 0.221
Quadrant decoding is best with the labels shifted by 0 to +2 s (0.493 to 0.519), falls at 5 s on either side (0.286 and 0.420) and lies around chance, 0.175 to 0.293, at 20 and 60 s (printed above). The behaviour timestamps and the imaging share a clock to within a second or two; offsets that small cannot be told apart by this check.
# Footprints deposited twice: pairs with near-identical traces and nearby positions.
Z = (T - T.mean(0)) / T.std(0)
C = (Z.T @ Z) / len(Z)
np.fill_diagonal(C, 0)
iu = np.triu_indices_from(C, 1)
corr = C[iu]
dist = np.hypot(*(S["centroids"][iu[0]] - S["centroids"][iu[1]]).T)
for thr in [0.7, 0.8, 0.9]:
k = corr > thr
print(f"pairs with trace correlation > {thr}: {k.sum()}; of them within 10 px: {(k & (dist <= 10)).sum()}")
print("largest trace correlation between two footprints:", round(float(corr.max()), 3))
print("NWB files in the dandiset (sessions):", sum(p.endswith(".nwb") for p in urls))
del Z, C
pairs with trace correlation > 0.7: 1; of them within 10 px: 0 pairs with trace correlation > 0.8: 0; of them within 10 px: 0 pairs with trace correlation > 0.9: 0; of them within 10 px: 0 largest trace correlation between two footprints: 0.77 NWB files in the dandiset (sessions): 1
There is one session, so no recording can be deposited twice, and within it no two footprints have traces correlated above 0.77.
# The findings of the first attempt (000691-scherrer-2023-01), checked again on the same version.
sub = requests.get(ks.code_file_url("f3+sf8/subregion_map_octant.m"), timeout=60).text
sub_q = requests.get(ks.code_file_url("f3+sf8/subregion_map_quaddir.m"), timeout=60).text
def code_value(text, name):
return next(l.split("=")[1].strip(" ;") for l in text.splitlines() if l.strip().startswith(name + " ="))
earlier = [
("data_missing: one mouse of four, no bird or optical data",
f"subjects {summary.get('numberOfSubjects')}, species {[s['name'] for s in summary.get('species', [])]}, assets {len(urls)}"),
("metadata_error: device described as the weight-optimized microscope", m["device_description"]),
("metadata_error: one-photon imaging stored as two-photon", f"{m['imaging_series']}; techniques "
f"{[t['name'] for t in summary.get('measurementTechnique', []) if 'two-photon' in t['name']]}"),
("metadata_error: behaviour video described as 30 fps", f"description '{m['video_description']}', AVI header "
f"{track['nominal_rate']:.0f} frames/s, timestamps {1 / np.median(np.diff(S['t_video'])):.2f} frames/s"),
("metadata_insufficient: no pixel size; indicator and location unknown",
f"grid spacing present {m['grid_spacing_present']}, indicator '{m['indicator']}', location '{m['location']}', "
f"excitation {m['excitation_lambda']}"),
("selection_undocumented: curated footprints not marked",
f"{T.shape[1]} footprints, columns {m['plane_segmentation_columns']}"),
("data_missing: no position of the mouse; opening of the video not marked",
f"processing modules {list(f['processing'].keys())}; no intervals"),
("code_mismatch: legend 32 footprints within 270 um, code 64 within 50 px",
f"subregion_map_octant.m nsize {code_value(sub, 'nsize')}, rsize {code_value(sub, 'rsize')}; "
f"subregion_map_quaddir.m nsize {code_value(sub_q, 'nsize')}, rsize {code_value(sub_q, 'rsize')}; 50 px = "
f"{50 * ks.UM_PER_PIXEL_ILLUMINATION:.0f} um at 5.94, {50 * ks.UM_PER_PIXEL_BEADS:.1f} um at 5.35"),
("paper_inconsistency: Extended Data Fig. 8 refers to Figure 3 and Fig. 2d", "in the paper's text, not the data"),
]
print(pd.DataFrame(earlier, columns=["finding of -01", "what the data and code show now"]).to_string(index=False))
finding of -01 what the data and code show now
data_missing: one mouse of four, no bird or optical data subjects 1, species ['Mus musculus - House mouse'], assets 2
metadata_error: device described as the weight-optimized microscope Head-mounted microscope (Kiloscope) optimized for weight (1.0 mm FOV; 1.0 g).
metadata_error: one-photon imaging stored as two-photon TwoPhotonSeries; techniques ['two-photon microscopy technique']
metadata_error: behaviour video described as 30 fps description 'Behavior video of animal moving in environment at ~30 fps', AVI header 30 frames/s, timestamps 12.51 frames/s
metadata_insufficient: no pixel size; indicator and location unknown grid spacing present False, indicator 'unknown', location 'unknown', excitation nan
selection_undocumented: curated footprints not marked 2948 footprints, columns ['image_mask', 'ROICentroids']
data_missing: no position of the mouse; opening of the video not marked processing modules ['ophys']; no intervals
code_mismatch: legend 32 footprints within 270 um, code 64 within 50 px subregion_map_octant.m nsize 64, rsize 50; subregion_map_quaddir.m nsize 64, rsize 50; 50 px = 297 um at 5.94, 267.5 um at 5.35
paper_inconsistency: Extended Data Fig. 8 refers to Figure 3 and Fig. 2d in the paper's text, not the data
Status of the findings of 000691-scherrer-2023-01, from the output above:
- One mouse of four, no bird or optical data: still present (one subject, two assets).
- The device description gives the weight-optimized microscope's FOV and weight: still present.
- One-photon imaging stored as a
TwoPhotonSeriesand listed as two-photon microscopy: still present. - The behaviour video described as about 30 fps, with timestamps at 12.5 frames/s: still present.
- No pixel size; indicator and location "unknown", no excitation wavelength: still present.
- The curated footprints not marked: still present.
- No position of the mouse; the opening of the video not marked: still present.
- The code's local maps use 64 footprints within 50 pixels where the legend says 32 within 270 um: still present for the number of footprints in both map scripts. The radius differs only at the illumination script's 5.94 um per pixel (297 um); at the bead script's 5.35 um per pixel, 50 pixels are 267.5 um, close to 270 um, so that part of the finding depends on the pixel size.
- Extended Data Fig. 8's legend refers to "Figure 3" and to a ridgeline plot in "Fig. 2d" (the ridgeline is Fig. 2h): still present in the version of record's full text, which the first attempt had not read.
New findings of this attempt, shown above and in the claim sections: the stored ROICentroids are
closer to the peak pixel than to a weighted centroid, and how they were computed is not said; the
code's unwrap gains an error of up to pi cm at each lap; the bead scripts order their images by file date, which the
repository does not keep (the Featherscope regions confirm the file-number order); and, in the
motion notebook, the movie's rigid shifts against a template are about a quarter of those the paper reports,
while nothing says whether the movie is raw or motion-corrected (its values are consistent with the raw video).
abstract-thousands-of-neurons¶
The Kiloscope records thousands of neurons at once in a freely behaving mouse; in the mouse with the best expression EXTRACT gave 2,948 footprints.
Criterion: the number of footprints EXTRACT returned for the one session is within 2948 +/- 295, and the footprints were recorded together in a session of free exploration of the maze.
n_fp = T.shape[1]
print(f"EXTRACT footprints in the session: {n_fp} (ids {len(S['ids'])}, masks {M['shape'][0]}, trace columns {T.shape[1]})")
print(f"expected 2948 +/- 295: within tolerance {abs(n_fp - 2948) <= 295}")
print(f"recorded together: one movie of {m['movie_shape'][0]} frames, traces of all footprints on its {len(S['t_img'])} timestamps")
print(f"session description: '{m['session_description']}'")
print(f"free exploration: {np.abs(np.diff(B['sunwrap'])).sum() / B['perimeter']:.1f} loop lengths in "
f"{B['t'][-1] - B['t'][0]:.0f} s of tracked video, quadrants visited {int((np.bincount(B['quadrant'], minlength=4) > 0).sum())} of 4, "
f"direction +1 in {np.mean(B['direction'] > 0):.2f} of frames")
EXTRACT footprints in the session: 2948 (ids 2948, masks 2948, trace columns 2948) expected 2948 +/- 295: within tolerance True recorded together: one movie of 36596 frames, traces of all footprints on its 36596 timestamps session description: 'This session includes calcium imaging recorded from a head-mounted microscope in a freely moving mouse while simultaneously recording more than a thousand neurons in cortex.' free exploration: 47.7 loop lengths in 1081 s of tracked video, quadrants visited 4 of 4, direction +1 in 0.53 of frames
Result: reproduced. The session holds 2948 EXTRACT footprints (masks, ids and traces agree), within 2948 +/- 295, recorded together in one movie of 36596 frames while the mouse ran 47.7 loop lengths and visited all four quadrants. These are the footprints before the paper's curation, which kept 2,382; the dandiset does not mark which.
Blocked claims: data not in the dandiset¶
ed6a-curated-footprints (2,382 after curation), ed6bcd-other-mice (1,029, 401 and 112 footprints in three other mice), ed5-songbirds and birds-singing-rate and fig2d-lman-sparse-bursting (the birds), abstract-weights (hardware weights), ed1-ed4-resolution (USAF target images) and ed3-adc (test tones on the converter). The output below shows what the dandiset and the code repository hold.
print("Dandiset: subjects", summary.get("numberOfSubjects"), "; species", [s["name"] for s in summary.get("species", [])],
"; assets", list(urls))
print("Plane segmentation columns:", m["plane_segmentation_columns"], "(no curation flag)")
top = ks.list_code_dir("")
print("Code repository at", ks.CODE_COMMIT[:7], "top level:", top)
for d in ["f2+sf4", "sf3/feather", "sf3/widefield", "sf1", "sf2"]:
print(f" {d}:", ks.list_code_dir(d)[:12])
Dandiset: subjects 1 ; species ['Mus musculus - House mouse'] ; assets ['sub-4mm-mouse-1/sub-4mm-mouse-1_image+ophys/60d277d6-24fd-4817-bb02-41382db59172_external_file_0.avi', 'sub-4mm-mouse-1/sub-4mm-mouse-1_image+ophys.nwb'] Plane segmentation columns: ['image_mask', 'ROICentroids'] (no curation flag)
Code repository at d596646 top level: ['LICENSE.md', 'README.md', 'f2+sf4', 'f3+sf8', 'sf1', 'sf10', 'sf2', 'sf3', 'sf9']
f2+sf4: ['feather', 'feather_analysis_secondorder.m', 'feather_roi.csv', 'wf', 'wf_analysis_secondorder.m', 'wf_roi.csv'] sf3/feather: ['analysis.m', 'axicons', 'imImport.m', 'led', 'reference']
sf3/widefield: ['analysis.m', 'axicons', 'imImport.m', 'led', 'reference', 'scalebar'] sf1: ['featherscope.ZDA', 'featherscope.zmx', 'kiloscope.ZDA', 'kiloscope.zmx']
sf2: ['axicons.nb']
Result: blocked. The dandiset holds one mouse and no curation flag, no bird data, and no optical test data (printed above). The code repository holds bead images and dye images (used below), Zemax files and an axicon notebook, but no USAF target images and no ADC recordings. ed6a-curated-footprints is blocked as metadata_insufficient (the data are there, the selection is not marked), fig2d-lman-sparse-bursting as methods_underspecified (no measure of sparse bursting is fixed; its data are also missing), and the others as data_missing. The paper's Academic Torrents deposit lists only this mouse's session in its file list (read outside this notebook, not used for any outcome).
abstract-kiloscope-fov¶
The Kiloscope's field of view is 4.8 x 3.6 mm (about 4 mm in the abstract).
Criterion: the frame's width and height in pixels times the object-plane pixel size; a pixel size stored with the data is used when there is one; otherwise the code's 5.94 and 5.35 um per pixel are the two readings, and the claim is reproduced when both sides are within 4.8 +/- 0.48 and 3.6 +/- 0.36 mm under a reading.
h_px, w_px = m["movie_shape"][1], m["movie_shape"][2]
print(f"frame: {h_px} x {w_px} pixels (movie axes 1 and 2); pixel size stored in the file: {m['grid_spacing_present']}")
rows = []
for name, um in [("5.94 um/px (sf3/widefield/analysis.m)", ks.UM_PER_PIXEL_ILLUMINATION), ("5.35 um/px (f2+sf4/wf_analysis_secondorder.m)", ks.UM_PER_PIXEL_BEADS)]:
long_mm, short_mm = h_px * um / 1000, w_px * um / 1000
rows.append({"reading": name, "longer_mm": round(long_mm, 3), "shorter_mm": round(short_mm, 3),
"longer_within": abs(long_mm - 4.8) <= 0.48, "shorter_within": abs(short_mm - 3.6) <= 0.36})
print(pd.DataFrame(rows).to_string(index=False))
ext = np.nanmax(POS["peak"], 0) - np.nanmin(POS["peak"], 0)
print(f"for reference, footprint positions span {ext[1]:.0f} x {ext[0]:.0f} px = {ext[1] * 5.94 / 1000:.2f} x {ext[0] * 5.94 / 1000:.2f} mm at 5.94, "
f"{ext[1] * 5.35 / 1000:.2f} x {ext[0] * 5.35 / 1000:.2f} mm at 5.35")
frame: 808 x 608 pixels (movie axes 1 and 2); pixel size stored in the file: False
reading longer_mm shorter_mm longer_within shorter_within
5.94 um/px (sf3/widefield/analysis.m) 4.800 3.612 True True
5.35 um/px (f2+sf4/wf_analysis_secondorder.m) 4.323 3.253 True True
for reference, footprint positions span 521 x 557 px = 3.09 x 3.31 mm at 5.94, 2.79 x 2.98 mm at 5.35
Result: reproduced. The file stores no pixel size, so the criterion's two readings from the code apply. Frames of 808 x 608 pixels are 4.800 x 3.612 mm at 5.94 um per pixel and 4.323 x 3.253 mm at 5.35 um per pixel; both sides are within tolerance under both readings, the second only just (4.323 against a lower bound of 4.32, 3.253 against 3.24). The pixel size comes from the code, not the dandiset.
fig1e-featherscope-psf and fig1f-kiloscope-psf¶
The Featherscope's resolution (FWHM of the bead PSF) is about 4 um at the centre of the field of view and 7 um at the edges. The Kiloscope's is about 4 um at the centre and about 5 um at the edges.
Criterion: as in f2+sf4/*_analysis_secondorder.m, each bead region at its listed slice is fitted
with the code's rotated 2D Gaussian with second-order pixel correction (1.68 and 5.35 um per
pixel); FWHM is 2.355 times each axis's width; radius from the code's centre, beads grouped to the
nearest of 0, 100, ..., 500 um (Featherscope) or 0, 250, ..., 1500 um (Kiloscope); centre is the 0
um group and edge the last. Readings: (a) median major-axis FWHM, (b) median minor-axis FWHM, (c)
mean major-axis, (d) mean minor-axis. Reproduced when centre and edge are both within 10% of 4 and
7 um (Featherscope) or 4 and 5 um (Kiloscope) under a reading.
The images and region lists are read from the code repository at commit d596646.
psf_setups = {
"fig1e-featherscope-psf": dict(folder="f2+sf4/feather", csv="f2+sf4/feather_roi.csv", dp=1.68, centre=(273, 400),
groups=np.arange(0, 501, 100), start_at_max=False, ftol=1e-6, max_nfev=10000, expect=(4, 7)),
"fig1f-kiloscope-psf": dict(folder="f2+sf4/wf/beads", csv="f2+sf4/wf_roi.csv", dp=5.35, centre=(259, 420),
groups=np.arange(0, 1501, 250), start_at_max=True, ftol=1e-8, max_nfev=100000, expect=(4, 5)),
}
stacks, rois = {}, {}
for name, s in psf_setups.items():
files = sorted((n for n in ks.list_code_dir(s["folder"]) if n.startswith("beads") and n.endswith(".png")),
key=lambda n: int(n[5:-4]))
with ThreadPoolExecutor(8) as ex:
stacks[name] = np.array(list(ex.map(lambda n: ks.read_png(ks.code_file_url(f"{s['folder']}/{n}")), files)))
tab = pd.read_csv(ks.code_file_url(s["csv"]))
tab.columns = ["row"] + list(tab.columns[1:])
rois[name] = tab
st = stacks[name]
print(f"{name}: {len(files)} images {files[0]} to {files[-1]}, each {st.shape[1:]} {st.dtype}, values {st.min()} to {st.max()}; "
f"{len(tab)} bead regions, slices {tab.Slice.min()} to {tab.Slice.max()}")
# The order of the images: the regions name their image (Featherscope), and the listed slice should be
# near the slice where the bead is brightest.
if tab.Label.str.contains("beads").all():
named = tab.Label.str.extract(r"beads(\d+)")[0].astype(int) + 1
print(f" regions whose label (beadsN) is slice N + 1 in file-number order: {(named == tab.Slice).sum()} of {len(tab)}")
peak = np.array([[img[r.BY:r.BY + r.Height, r.BX:r.BX + r.Width].max() for img in st] for _, r in tab.iterrows()])
off = peak.argmax(1) + 1 - tab.Slice.to_numpy()
print(f" brightest slice minus listed slice, percentiles 10/50/90: {np.percentile(off, [10, 50, 90])}; within 5 slices: {np.mean(abs(off) <= 5):.2f}")
fig1e-featherscope-psf: 161 images beads0.png to beads160.png, each (608, 808) uint8, values 6 to 255; 155 bead regions, slices 1 to 161 regions whose label (beadsN) is slice N + 1 in file-number order: 155 of 155
brightest slice minus listed slice, percentiles 10/50/90: [-9. 0. 8.]; within 5 slices: 0.71
fig1f-kiloscope-psf: 91 images beads0.png to beads90.png, each (608, 808) uint8, values 14 to 255; 118 bead regions, slices 13 to 63
brightest slice minus listed slice, percentiles 10/50/90: [-12. 0. 6.6]; within 5 slices: 0.72
def psf_fits(name, x_scale="code"):
s, tab, st = psf_setups[name], rois[name], stacks[name]
rows = []
for _, r in tab.iterrows():
roi = st[int(r.Slice) - 1, r.BY:r.BY + r.Height, r.BX:r.BX + r.Width]
fit = ks.fit_bead(roi, s["dp"], s["start_at_max"], s["ftol"], s["max_nfev"], x_scale)
hh, hw = (r.Height - 1) / 2, (r.Width - 1) / 2
rad = s["dp"] * np.hypot(s["centre"][0] - (r.BY + hh + 1), s["centre"][1] - (r.BX + hw + 1))
rows.append({"r_um": rad, "slice": r.Slice, **fit})
out = pd.DataFrame(rows)
g = s["groups"]
out["group"] = g[np.abs(out.r_um.to_numpy()[:, None] - g[None]).argmin(1)]
return out
def readings(fits, name):
s = psf_setups[name]
by = fits.groupby("group")
table = pd.DataFrame({"n": by.size(), "median_major": by.fwhm_major.median(), "median_minor": by.fwhm_minor.median(),
"mean_major": by.fwhm_major.mean(), "mean_minor": by.fwhm_minor.mean()}).round(3)
centre_g, edge_g = s["groups"][0], s["groups"][-1]
rows = []
for label, col in [("(a) median major", "median_major"), ("(b) median minor", "median_minor"),
("(c) mean major", "mean_major"), ("(d) mean minor", "mean_minor")]:
cv, ev = table.loc[centre_g, col], table.loc[edge_g, col]
rows.append({"reading": label, "centre_um": cv, "edge_um": ev,
"centre_within": abs(cv - s["expect"][0]) <= 0.1 * s["expect"][0],
"edge_within": abs(ev - s["expect"][1]) <= 0.1 * s["expect"][1]})
return table, pd.DataFrame(rows)
psf = {}
for name in psf_setups:
t0 = time.time()
psf[name] = psf_fits(name)
fits = psf[name]
table, read = readings(fits, name)
print(f"\n{name}: {len(fits)} fits in {time.time() - t0:.0f} s; stopped by the evaluation limit: {(fits.status == 0).sum()}; "
f"FWHM major {fits.fwhm_major.min():.2f} to {fits.fwhm_major.max():.2f} um, minor {fits.fwhm_minor.min():.2f} to {fits.fwhm_minor.max():.2f} um")
print(table.to_string())
print(read.to_string(index=False))
print(" a reading with both within tolerance:", bool((read.centre_within & read.edge_within).any()))
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for a_, name in zip(ax, psf):
fits = psf[name]
a_.scatter(fits.r_um, fits.fwhm_major, s=8, label="major axis")
a_.scatter(fits.r_um, fits.fwhm_minor, s=8, label="minor axis")
a_.set_ylim(0, 12); a_.set_xlabel("radius (um)"); a_.set_ylabel("FWHM (um)"); a_.set_title(name, fontsize=9); a_.legend(fontsize=7)
plt.tight_layout(); plt.show()
for name in psf:
print(f"{name}: points above the plotted 12 um: {(psf[name].fwhm_major > 12).sum()}")
fig1e-featherscope-psf: 155 fits in 1 s; stopped by the evaluation limit: 0; FWHM major 2.48 to 9.89 um, minor 2.18 to 7.94 um
n median_major median_minor mean_major mean_minor
group
0 11 3.919 2.609 3.887 2.690
100 27 3.975 2.879 3.932 2.907
200 34 3.902 3.037 3.875 3.063
300 33 4.548 3.887 4.900 3.962
400 31 5.972 4.562 6.030 4.636
500 19 6.433 5.203 6.883 5.244
reading centre_um edge_um centre_within edge_within
(a) median major 3.919 6.433 True True
(b) median minor 2.609 5.203 False False
(c) mean major 3.887 6.883 True True
(d) mean minor 2.690 5.244 False False
a reading with both within tolerance: True
fig1f-kiloscope-psf: 118 fits in 204 s; stopped by the evaluation limit: 1; FWHM major 2.05 to 27.77 um, minor 1.05 to 5.64 um
n median_major median_minor mean_major mean_minor
group
0 11 4.467 2.722 4.537 2.669
250 14 4.361 2.535 4.793 2.546
500 27 4.712 3.955 5.092 3.798
750 26 5.480 3.473 5.473 3.242
1000 19 6.018 3.594 6.072 3.440
1250 10 5.126 3.230 5.410 3.200
1500 11 5.441 3.306 7.373 3.094
reading centre_um edge_um centre_within edge_within
(a) median major 4.467 5.441 False True
(b) median minor 2.722 3.306 False False
(c) mean major 4.537 7.373 False False
(d) mean minor 2.669 3.094 False False
a reading with both within tolerance: False
fig1e-featherscope-psf: points above the plotted 12 um: 0 fig1f-kiloscope-psf: points above the plotted 12 um: 1
# Variant: the same fits with the optimizer's steps scaled by the Jacobian (scipy's x_scale="jac"), to show how
# much the fits depend on the optimizer.
psf_jac = {}
for name in psf_setups:
fits_j = psf_jac[name] = psf_fits(name, x_scale="jac")
_, read_j = readings(fits_j, name)
dif = np.abs(fits_j.fwhm_major - psf[name].fwhm_major)
print(f"\n{name}, Jacobian scaling: beads whose major FWHM moves by more than 0.5 um: {(dif > 0.5).sum()} of {len(dif)}; "
f"fits with a lower cost than the main fit: {(fits_j.cost < psf[name].cost - 1e-6).sum()}")
print(read_j.to_string(index=False))
fig1e-featherscope-psf, Jacobian scaling: beads whose major FWHM moves by more than 0.5 um: 0 of 155; fits with a lower cost than the main fit: 60
reading centre_um edge_um centre_within edge_within
(a) median major 3.919 6.433 True True
(b) median minor 2.609 5.203 False False
(c) mean major 3.887 6.884 True True
(d) mean minor 2.690 5.244 False False
fig1f-kiloscope-psf, Jacobian scaling: beads whose major FWHM moves by more than 0.5 um: 15 of 118; fits with a lower cost than the main fit: 46
reading centre_um edge_um centre_within edge_within
(a) median major 3.639 5.892 True False
(b) median minor 2.387 2.634 False False
(c) mean major 4.261 6.135 True False
(d) mean minor 2.513 2.878 False False
# For each bead, the lower-cost of the two fits (the main fit and the Jacobian-scaled one), the one nearer the
# least-squares minimum that lsqcurvefit seeks.
psf_best = {}
for name in psf_setups:
main_, jac_ = psf[name], psf_jac[name]
use_j = (jac_.cost < main_.cost).to_numpy()
best = psf_best[name] = main_.copy()
best.loc[use_j] = jac_.loc[use_j]
_, read_b = readings(best, name)
print(f"\n{name}, the lower-cost fit for each bead (the Jacobian-scaled fit for {use_j.sum()} of {len(best)} beads):")
print(read_b.to_string(index=False))
print(" a reading with both within tolerance:", bool((read_b.centre_within & read_b.edge_within).any()))
print(f" largest major FWHM: main fits {np.round(np.sort(main_.fwhm_major)[-3:], 2).tolist()}, lower-cost fits "
f"{np.round(np.sort(best.fwhm_major)[-3:], 2).tolist()} um")
# The Kiloscope's centre and edge beads, each with both fits.
main_, jac_ = psf["fig1f-kiloscope-psf"], psf_jac["fig1f-kiloscope-psf"]
both = pd.DataFrame({"group": main_.group, "r_um": main_.r_um, "major_main": main_.fwhm_major, "minor_main": main_.fwhm_minor,
"cost_main": main_.cost, "major_jac": jac_.fwhm_major, "minor_jac": jac_.fwhm_minor, "cost_jac": jac_.cost})
print("\nfig1f-kiloscope-psf, beads of the 0 and 1500 um groups:")
print(both[both.group.isin([0, 1500])].round(3).to_string())
fig1e-featherscope-psf, the lower-cost fit for each bead (the Jacobian-scaled fit for 69 of 155 beads):
reading centre_um edge_um centre_within edge_within
(a) median major 3.919 6.433 True True
(b) median minor 2.609 5.203 False False
(c) mean major 3.887 6.884 True True
(d) mean minor 2.690 5.244 False False
a reading with both within tolerance: True
largest major FWHM: main fits [9.15, 9.26, 9.89], lower-cost fits [9.15, 9.26, 9.89] um
fig1f-kiloscope-psf, the lower-cost fit for each bead (the Jacobian-scaled fit for 66 of 118 beads):
reading centre_um edge_um centre_within edge_within
(a) median major 4.467 5.892 False False
(b) median minor 2.387 2.634 False False
(c) mean major 4.685 7.944 False False
(d) mean minor 2.509 2.904 False False
a reading with both within tolerance: False
largest major FWHM: main fits [8.81, 9.51, 27.77], lower-cost fits [9.51, 14.2, 27.77] um
fig1f-kiloscope-psf, beads of the 0 and 1500 um groups:
group r_um major_main minor_main cost_main major_jac minor_jac cost_jac
0 0 93.434 4.084 1.680 461.402 8.842 1.838 451.498
3 0 32.100 4.467 4.017 166.698 4.466 4.018 166.698
4 0 92.858 3.375 2.722 232.050 2.775 2.311 232.021
30 1500 1536.577 4.755 3.343 132.365 3.139 2.069 131.398
31 1500 1521.652 27.768 2.697 60.092 5.160 1.369 79.942
32 1500 1438.782 8.387 3.514 195.227 8.386 3.514 195.227
33 1500 1479.858 6.566 3.362 228.591 6.567 3.362 228.591
34 1500 1405.455 3.235 1.811 123.280 5.938 2.850 132.656
50 0 32.100 4.779 2.948 134.915 3.631 2.582 134.649
60 0 93.434 6.421 1.782 241.677 3.167 1.869 253.870
61 0 32.100 4.467 4.017 166.698 4.466 4.018 166.698
62 0 57.621 5.727 2.387 352.206 5.728 2.387 352.206
63 0 45.710 3.982 3.403 146.270 3.067 2.473 143.310
64 0 81.664 3.696 2.206 172.476 3.639 2.175 172.476
65 0 94.651 5.476 1.345 208.288 4.068 1.300 221.333
71 0 108.329 3.433 2.850 165.002 3.024 2.669 165.000
111 1500 1624.897 6.313 4.978 324.973 6.313 4.978 324.973
112 1500 1490.590 3.853 2.359 402.624 3.176 2.377 402.605
113 1500 1632.802 3.687 2.729 240.910 3.267 2.634 240.909
115 1500 1453.008 5.892 3.478 342.750 5.892 3.478 342.750
116 1500 1439.441 5.441 2.456 162.527 5.441 2.456 162.527
117 1500 1407.330 5.205 3.306 158.925 14.202 2.569 129.867
Result: reproduced (fig1e-featherscope-psf). Under reading (a) the median major-axis FWHM is 3.919 um at the centre (11 beads) and 6.433 um in the 500 um group (19 beads), within 10% of 4 and 7 um; reading (c) gives 3.887 and 6.883 um (6.884 with Jacobian scaling), also within. The minor axis is narrower (2.609 and 5.203 um under b) and fails at both positions. All 155 regions name their image in file-number order. With Jacobian scaling no bead's major FWHM moves by more than 0.5 um, and taking for each bead the lower-cost of the two fits gives the same readings (3.919 and 6.433 um under a).
Result: partially reproduced (fig1f-kiloscope-psf). No reading puts both values within tolerance, with either optimizer or with the lower-cost fit for each bead. Under (a) the edge (1500 um group, 11 beads) is 5.441 um, within 5 +/- 0.5, and the centre (11 beads) is 4.467 um, 0.067 um above the bound of 4.4. The mean major axis (c) is 4.537 um at the centre and 7.373 um at the edge, where one bead is fitted at 27.77 um; the minor axis gives 2.722 and 3.306 um (b). The fits depend on the optimizer, as the PSF is narrower than a pixel (5.35 um): with Jacobian scaling 15 of 118 beads' major FWHM move by more than 0.5 um and (a) gives 3.639 um at the centre (within) and 5.892 um at the edge (0.392 um above the bound). Taking for each bead the lower-cost fit, the one nearer the least-squares minimum, gives 4.467 and 5.892 um under (a), both outside, and 4.685 and 7.944 um under (c); the 27.77 um edge bead keeps its main fit, which has the lower cost (60.092 against 79.942). These fits are not the same in every execution: the seeds are fixed, but the values printed above differ between runs on different machines (the attempt record gives those of the run reviewed in the second verification, where the lower-cost fits put the centre within tolerance and the edge outside). The listed slice is within 5 slices of the bead's brightest slice for 72% of the Kiloscope's regions, and the image order of that stack is not confirmed by its labels. Under every reading and every fit the FWHM is larger at the edge than at the centre, and under (a) one of the two positions is within tolerance with the main fits (the edge) and with the Jacobian-scaled fits (the centre), so the claim holds in part.
ed1cf-field-curvature¶
Both microscopes have field curvature small enough for typical imaging depths.
The paper gives no bound on field curvature and no typical imaging depth, so no method it fixes
decides the claim (methods_underspecified). As the criterion asks, the slice of best focus (the
Slice column of the region lists, over the 100 um stack of the Methods) is reported against
radius for both microscopes.
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for a_, name in zip(ax, psf):
fits = psf[name]
by = fits.groupby("group").slice
print(f"{name}: images in the stack {len(stacks[name])}; slice of best focus by radius group:")
print(pd.DataFrame({"n": by.size(), "median": by.median(), "min": by.min(), "max": by.max()}).to_string())
print(f" Spearman correlation of slice with radius: {pd.Series(fits.slice.to_numpy()).corr(pd.Series(fits.r_um.to_numpy()), method='spearman'):.2f}")
a_.scatter(fits.r_um, fits.slice, s=8); a_.set_xlabel("radius (um)"); a_.set_ylabel("slice of best focus"); a_.set_title(name.split("-")[1], fontsize=9)
plt.tight_layout(); plt.show()
fig1e-featherscope-psf: images in the stack 161; slice of best focus by radius group:
n median min max
group
0 11 6.0 1 7
100 27 8.0 1 13
200 34 13.0 8 21
300 33 20.0 9 29
400 31 50.0 20 76
500 19 89.0 38 161
Spearman correlation of slice with radius: 0.96
fig1f-kiloscope-psf: images in the stack 91; slice of best focus by radius group:
n median min max
group
0 11 48.0 43 52
250 14 47.5 39 54
500 27 45.0 26 63
750 26 37.0 25 53
1000 19 32.0 22 52
1250 10 26.0 21 37
1500 11 29.0 13 47
Spearman correlation of slice with radius: -0.66
Result: blocked (methods_underspecified). The paper sets no bound on field curvature, so nothing decides the claim. As the criterion asks, the slice of best focus is shown: for the Featherscope it rises from a median of 6 at the centre to 89 in the 500 um group (of 161 slices; Spearman 0.96 with radius), and for the Kiloscope it falls from 48 at the centre to 29 in the 1500 um group (of 91 slices; Spearman -0.66). The Methods say the stack spans 100 um; if its slices are evenly spaced, these are about 52 um and 21 um of focal shift across the field.
ed2b-featherscope-flat-illumination and ed2c-kiloscope-flat-illumination¶
With the axicon correction optics the signal profile is flat across the field of view, where total vignetting makes the signal fall toward the edges (Extended Data Fig. 2b for the Featherscope, 2c for the Kiloscope).
Criterion: radial profiles as sf3/*/analysis.m computes them; the corrected profile is the
axicons image, total vignetting the reference image, optical vignetting the led image. Edge radius
readings 0.5 and 0.45 mm (Featherscope), 1.8 and 2.0 mm (Kiloscope). Reproduced when, under both,
the corrected profile at the edge radius is higher than the total-vignetting profile and its fall
from the centre is smaller.
The code's profiles run from 0 to rLength - 1 um (499 and 1999 um), so the 0.5 mm and 2.0 mm readings are taken at the last bin, 1 um short; a profile extended to include 500 and 2000 um is shown as a check (a standard choice that cannot change the outcome).
setups = {
"ed2b-featherscope": dict(folder="sf3/feather", um=1.68, r_length=500, med=3, radii=(500, 450),
files={"corrected": "axicons/rgb2021-01-29T16_31_30.png", "optical vignetting": "led/feather.png",
"total vignetting": "reference/rgb2021-01-29T16_54_50.png"}),
"ed2c-kiloscope": dict(folder="sf3/widefield", um=5.94, r_length=2000, med=10, radii=(1800, 2000),
files={"corrected": "axicons/rgb2021-01-08T14_40_22.png", "optical vignetting": "led/rgb2021-01-29T15_25_25.png",
"total vignetting": "reference/rgb2021-01-08T15_23_50.png"}),
}
images, profiles, extended = {}, {}, {}
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for a_, (name, s) in zip(ax, setups.items()):
images[name], profiles[name], extended[name] = {}, {}, {}
for label, fn in s["files"].items():
img = ks.read_png(ks.code_file_url(f"{s['folder']}/{fn}"))
img = ks.unpack_rgb565(img) if img.ndim == 3 else img.astype(float)
images[name][label] = img
profiles[name][label] = ks.radial_profile(img, s["um"], s["r_length"], s["med"])
extended[name][label] = ks.radial_profile(img, s["um"], max(s["radii"]) + 1, s["med"])
a_.plot(profiles[name][label], label=label)
print(f"{name} {label}: image {img.shape}, values {img.min():.0f} to {img.max():.0f}, weighted centre "
f"{tuple(round(v, 1) for v in ks.weighted_centre(img))}")
a_.set_xlim(0, s["r_length"]); a_.set_ylim(0, 1.1)
for R in s["radii"]:
a_.axvline(min(R, s["r_length"] - 1), color="k", ls=":")
a_.set_xlabel("radius (um)"); a_.set_title(name); a_.legend(fontsize=7)
plt.show()
for name, s in setups.items():
print(f"\n{name}: radii with values above the plotted 1.1:", {k: int(np.sum(p > 1.1)) for k, p in profiles[name].items()})
rows = []
for R in s["radii"]:
for kind, prof_set in [("code's profile", profiles[name]), ("extended profile", extended[name])]:
Rb = min(R, len(prof_set["corrected"]) - 1)
corr_, tot_, opt_ = (prof_set[k] for k in ("corrected", "total vignetting", "optical vignetting"))
fall_c, fall_t = corr_[0] - corr_[Rb], tot_[0] - tot_[Rb]
rows.append({"edge_reading_um": R, "profile": kind, "bin_um": Rb, "corrected": round(corr_[Rb], 3),
"total": round(tot_[Rb], 3), "optical": round(opt_[Rb], 3), "fall_corrected": round(fall_c, 3),
"fall_total": round(fall_t, 3), "corrected_higher": corr_[Rb] > tot_[Rb], "fall_smaller": fall_c < fall_t})
print(pd.DataFrame(rows).to_string(index=False))
ed2b-featherscope corrected: image (608, 808), values 23 to 716, weighted centre (np.float64(390.2), np.float64(324.8))
ed2b-featherscope optical vignetting: image (608, 808), values 18 to 195, weighted centre (np.float64(398.4), np.float64(322.9))
ed2b-featherscope total vignetting: image (608, 808), values 17 to 778, weighted centre (np.float64(379.8), np.float64(299.6))
ed2c-kiloscope corrected: image (608, 808), values 23 to 444, weighted centre (np.float64(439.0), np.float64(299.5))
ed2c-kiloscope optical vignetting: image (608, 808), values 37 to 489, weighted centre (np.float64(417.4), np.float64(285.2))
ed2c-kiloscope total vignetting: image (608, 808), values 27 to 888, weighted centre (np.float64(426.2), np.float64(288.9))
ed2b-featherscope: radii with values above the plotted 1.1: {'corrected': 0, 'optical vignetting': 0, 'total vignetting': 0}
edge_reading_um profile bin_um corrected total optical fall_corrected fall_total corrected_higher fall_smaller
500 code's profile 499 0.511 0.193 0.494 0.470 0.814 True True
500 extended profile 500 0.507 0.192 0.494 0.475 0.815 True True
450 code's profile 450 0.628 0.257 0.561 0.353 0.750 True True
450 extended profile 450 0.628 0.257 0.561 0.353 0.750 True True
ed2c-kiloscope: radii with values above the plotted 1.1: {'corrected': 0, 'optical vignetting': 0, 'total vignetting': 0}
edge_reading_um profile bin_um corrected total optical fall_corrected fall_total corrected_higher fall_smaller
1800 code's profile 1800 0.557 0.180 0.309 0.437 0.843 True True
1800 extended profile 1800 0.557 0.180 0.309 0.437 0.843 True True
2000 code's profile 1999 0.420 0.148 0.272 0.573 0.875 True True
2000 extended profile 2000 0.419 0.148 0.272 0.574 0.875 True True
# Sensitivity of the comparison to the centre: the corrected and total-vignetting images about other centres.
for name, s in setups.items():
rows = []
for label, shift in [("intensity-weighted (main)", (0, 0)), ("+100 um in x", (1, 0)), ("-100 um in x", (-1, 0)),
("+100 um in y", (0, 1)), ("-100 um in y", (0, -1)), ("image centre", None)]:
row = {"centre": label}
for R in s["radii"]:
Rb = min(R, s["r_length"] - 1)
vals = {}
for k in ("corrected", "total vignetting"):
img = images[name][k]
xc, yc = ks.weighted_centre(img)
step = 100 / s["um"]
centre = ((img.shape[1] + 1) / 2, (img.shape[0] + 1) / 2) if shift is None else (xc + shift[0] * step, yc + shift[1] * step)
p = ks.radial_profile(img, s["um"], s["r_length"], s["med"], centre=centre)
vals[k] = (p[Rb], p[0] - p[Rb])
row[f"holds_at_{R}um"] = vals["corrected"][0] > vals["total vignetting"][0] and vals["corrected"][1] < vals["total vignetting"][1]
row[f"corrected_{R}"] = round(vals["corrected"][0], 3)
row[f"total_{R}"] = round(vals["total vignetting"][0], 3)
rows.append(row)
print(f"\n{name}:")
print(pd.DataFrame(rows).to_string(index=False))
ed2b-featherscope:
centre holds_at_500um corrected_500 total_500 holds_at_450um corrected_450 total_450
intensity-weighted (main) True 0.511 0.193 True 0.628 0.257
+100 um in x True 0.453 0.207 True 0.547 0.268
-100 um in x True 0.581 0.220 True 0.700 0.287
+100 um in y True 0.475 0.234 True 0.552 0.293
-100 um in y True 0.680 0.253 True 0.779 0.323
image centre True 0.526 0.193 True 0.651 0.256
ed2c-kiloscope:
centre holds_at_1800um corrected_1800 total_1800 holds_at_2000um corrected_2000 total_2000
intensity-weighted (main) True 0.557 0.180 True 0.420 0.148
+100 um in x True 0.584 0.190 True 0.435 0.156
-100 um in x True 0.543 0.175 True 0.416 0.145
+100 um in y True 0.563 0.191 True 0.419 0.154
-100 um in y True 0.562 0.177 True 0.420 0.146
image centre True 0.510 0.177 True 0.402 0.145
# Which images sf3/*/analysis.m loads: its lines that read an image, against the files used above.
import re
for name, s in setups.items():
script = requests.get(ks.code_file_url(f"{s['folder']}/analysis.m"), timeout=60).text
lines = [(i + 1, line.strip()) for i, line in enumerate(script.splitlines()) if re.search(r"imImport\(|imread\(", line)]
print(f"{s['folder']}/analysis.m:")
for i, line in lines:
print(f" line {i}: {line}")
loaded = [re.search(r"'([^']+)'", line).group(1) for _, line in lines]
print(f" the files used here, in the order corrected, optical vignetting, total vignetting, are the ones it loads as axicon, "
f"led and ref: {loaded == list(s['files'].values())}")
sf3/feather/analysis.m:
line 20: axicon = double(imImport('axicons/rgb2021-01-29T16_31_30.png'));
line 21: led = double(imread('led/feather.png'));
line 22: ref = double(imImport('reference/rgb2021-01-29T16_54_50.png'));
the files used here, in the order corrected, optical vignetting, total vignetting, are the ones it loads as axicon, led and ref: True
sf3/widefield/analysis.m:
line 18: axicon = double(imImport('axicons/rgb2021-01-08T14_40_22.png'));
line 19: led = double(imImport('led/rgb2021-01-29T15_25_25.png'));
line 20: ref = double(imImport('reference/rgb2021-01-08T15_23_50.png'));
the files used here, in the order corrected, optical vignetting, total vignetting, are the ones it loads as axicon, led and ref: True
Result: reproduced (ed2b-featherscope-flat-illumination). At both edge readings the corrected
profile is higher than the total-vignetting profile and falls less from the centre: at 0.5 mm (the
code's last bin, 499 um) 0.511 against 0.193, falling by 0.470 against 0.814; at 0.45 mm 0.628
against 0.257, falling by 0.353 against 0.750. The bin at 500 um of an extended profile differs by at
most 0.005. With the centre moved by 100 um in either axis, or at the image centre, the comparison
holds at both radii. The three images are the ones sf3/feather/analysis.m loads, in the same
roles (its lines 20 to 22, printed above). The corrected profile still falls to about half its peak at 0.5 mm; the
criterion judges it against total vignetting, as the legend does.
Result: reproduced (ed2c-kiloscope-flat-illumination). At 1.8 mm the corrected profile is 0.557
against 0.180 for total vignetting (falling by 0.437 against 0.843), and at 2.0 mm (the code's
last bin, 1999 um) 0.420 against 0.148 (falling by 0.573 against 0.875). The comparison holds at
both radii for all six centres tried. The three images are the ones sf3/widefield/analysis.m loads, in the same roles
(its lines 18 to 20).
ed8-coactive-distributed¶
The correlated activity in each NMF component is driven by footprints spread across the field of view.
Criterion: NMF of the traces into eight components (multiplicative updates, 4 replicates, tolerances 1e-7); each footprint to the component in which its weight ranks highest; footprints of a component are active in its event when their maximum within 25 s of the event time exceeds half their maximum over the recording. Event times: (a) each component's largest temporal weight; (b) the code's 820, 770, 990, 700, 400, 420, 1000 and 410 s, each assigned one to one to the component with the largest temporal weights at those times. Reproduced when every component's event has active footprints in all four quarters of the field of view (split at the midlines of the extent of all footprint positions, positions as for fig2l: the first pixel of weight 1) under both readings.
The curated set of 2,382 is not marked in the dandiset, so all 2,948 footprints are used, the variant the claims file names.
t0 = time.time()
nmf = ks.nmf_replicates(T, 8, 4, seed=691_08)
W, H = nmf["W"], nmf["H"]
group = ks.nmf_groups(H)
t_tr = S["t_trace"]
print(f"NMF in {time.time() - t0:.0f} s; errors of the 4 replicates {np.round(nmf['errors'], 3)} (iterations {nmf['n_iter']})")
print("footprints per component:", np.bincount(group, minlength=8).tolist())
ev_a = t_tr[W.argmax(0)]
code_times = np.array([820, 770, 990, 700, 400, 420, 1000, 410])
W_at = W[np.searchsorted(t_tr, code_times)] # rows: the code's times; columns: components
rows_, cols_ = linear_sum_assignment(-W_at)
ev_b = np.empty(8); ev_b[cols_] = code_times[rows_]
rows_n, cols_n = linear_sum_assignment(-(W_at / W.max(0)))
ev_b_norm = np.empty(8); ev_b_norm[cols_n] = code_times[rows_n]
print("event times (a), largest temporal weight (s):", np.round(ev_a, 1).tolist())
print("event times (b), the code's times assigned to components (s):", ev_b.tolist())
print(" assignment with each component's weights scaled to their maximum gives the same times:", bool(np.array_equal(ev_b, ev_b_norm)))
print(" temporal weight at the assigned time / component maximum:", np.round(W[np.searchsorted(t_tr, ev_b), np.arange(8)] / W.max(0), 2).tolist())
peak_all = T.max(0)
def events(times, positions):
pos = positions
rmid = (np.nanmin(pos[:, 0]) + np.nanmax(pos[:, 0])) / 2
cmid = (np.nanmin(pos[:, 1]) + np.nanmax(pos[:, 1])) / 2
quarter = 2 * (pos[:, 0] > rmid) + (pos[:, 1] > cmid)
rows, actives = [], []
for c in range(8):
members = np.flatnonzero(group == c)
win = np.abs(t_tr - times[c]) <= 25
act_ = members[T[win][:, members].max(0) > 0.5 * peak_all[members]]
actives.append(act_)
q = np.bincount(quarter[act_], minlength=4)
rows.append({"component": c, "event_t_s": round(times[c], 1), "members": len(members), "active": len(act_),
"upper_left": q[0], "upper_right": q[1], "lower_left": q[2], "lower_right": q[3], "all_four": bool((q > 0).all())})
return pd.DataFrame(rows), actives, quarter
results = {}
for label, times in [("(a) largest weight", ev_a), ("(b) the code's times", ev_b)]:
for pos_label, pos in [("code's position", POS["peak"]), ("centroid (variant)", POS["centroid"])]:
tab, actives, quarter = events(times, pos)
results[(label, pos_label)] = (tab, actives)
print(f"\nEvent times {label}, positions: {pos_label}; events with active footprints in all four quarters: {tab.all_four.sum()} of 8")
if pos_label == "code's position":
print(tab.to_string(index=False))
print(" all footprints per quarter:", np.bincount(quarter, minlength=4).tolist())
else:
print(" active per quarter, smallest count over events and quarters:", int(tab[["upper_left", "upper_right", "lower_left", "lower_right"]].to_numpy().min()))
NMF in 35 s; errors of the 4 replicates [97.604 97.562 97.565 97.603] (iterations [100, 100, 100, 100])
footprints per component: [389, 388, 329, 381, 333, 385, 405, 338]
event times (a), largest temporal weight (s): [286.4, 953.0, 791.5, 1150.7, 518.1, 1165.2, 122.9, 80.4]
event times (b), the code's times assigned to components (s): [770.0, 1000.0, 700.0, 410.0, 990.0, 420.0, 820.0, 400.0]
assignment with each component's weights scaled to their maximum gives the same times: True
temporal weight at the assigned time / component maximum: [0.4699999988079071, 0.10999999940395355, 0.03999999910593033, 0.05999999865889549, 0.550000011920929, 0.05000000074505806, 0.12999999523162842, 0.6800000071525574]
Event times (a) largest weight, positions: code's position; events with active footprints in all four quarters: 8 of 8
component event_t_s members active upper_left upper_right lower_left lower_right all_four
0 286.4 389 223 80 103 6 34 True
1 953.0 388 250 103 120 1 26 True
2 791.5 329 215 72 112 4 27 True
3 1150.7 381 326 122 156 16 32 True
4 518.1 333 192 72 85 3 32 True
5 1165.2 385 308 137 130 10 31 True
6 122.9 405 280 115 119 11 35 True
7 80.4 338 260 90 116 7 47 True
all footprints per quarter: [1114, 1378, 85, 371]
Event times (a) largest weight, positions: centroid (variant); events with active footprints in all four quarters: 8 of 8
active per quarter, smallest count over events and quarters: 1
Event times (b) the code's times, positions: code's position; events with active footprints in all four quarters: 7 of 8
component event_t_s members active upper_left upper_right lower_left lower_right all_four
0 770.0 389 219 79 99 8 33 True
1 1000.0 388 151 55 79 0 17 False
2 700.0 329 125 46 64 2 13 True
3 410.0 381 136 47 71 10 8 True
4 990.0 333 208 94 76 5 33 True
5 420.0 385 134 67 48 4 15 True
6 820.0 405 222 91 93 11 27 True
7 400.0 338 158 42 92 3 21 True
all footprints per quarter: [1114, 1378, 85, 371]
Event times (b) the code's times, positions: centroid (variant); events with active footprints in all four quarters: 7 of 8
active per quarter, smallest count over events and quarters: 0
# Variant: the code's window, frames lb * 30 to ub * 30 (single_component_raster.m), which at the imaging's 29.7
# frames/s is not exactly 25 s on either side of the time.
pos = POS["peak"]
rmid = (np.nanmin(pos[:, 0]) + np.nanmax(pos[:, 0])) / 2
cmid = (np.nanmin(pos[:, 1]) + np.nanmax(pos[:, 1])) / 2
quarter = 2 * (pos[:, 0] > rmid) + (pos[:, 1] > cmid)
rows = []
for c in range(8):
members = np.flatnonzero(group == c)
lo, hi = int((ev_b[c] - 25) * 30) - 1, int((ev_b[c] + 25) * 30)
act_ = members[T[lo:hi][:, members].max(0) > 0.5 * peak_all[members]]
q = np.bincount(quarter[act_], minlength=4)
rows.append({"component": c, "time_s": ev_b[c], "window_s": f"{t_tr[lo]:.1f} to {t_tr[hi - 1]:.1f}", "active": len(act_),
"per_quarter": q.tolist(), "all_four": bool((q > 0).all())})
print("Variant of reading (b) with the code's frame window:")
print(pd.DataFrame(rows).to_string(index=False))
# Maps of the active footprints (the code's Extended Data Fig. 8b): footprints thresholded at 10% of their
# maximum; pixels with one active footprint blue, two or more red, inactive footprints grey.
roi, a1, a2, w = M["roi"], M["a1"], M["a2"], M["w"]
above = w > 0.1 * POS["max_weight"][roi]
outline = np.zeros((M["shape"][2], M["shape"][1]), bool)
outline[a2[above], a1[above]] = True
for label in ["(a) largest weight", "(b) the code's times"]:
tab, actives = results[(label, "code's position")]
fig, ax = plt.subplots(2, 4, figsize=(14, 5.5))
for c, a_ in enumerate(ax.ravel()):
count = np.zeros_like(outline, dtype=int)
k = above & np.isin(roi, actives[c])
np.add.at(count, (a2[k], a1[k]), 1)
rgb = np.ones((*outline.shape, 3))
rgb[outline] = 0.6
rgb[count == 1] = (0, 0.5, 1)
rgb[count >= 2] = (1, 0, 0)
a_.imshow(rgb); a_.axis("off")
a_.set_title(f"component {c}: {len(actives[c])} of {tab.members[c]} active, t = {tab.event_t_s[c]:.0f} s", fontsize=8)
fig.suptitle(f"Event times {label}", fontsize=10)
plt.tight_layout(); plt.show()
Variant of reading (b) with the code's frame window:
component time_s window_s active per_quarter all_four
0 770.0 753.9 to 804.5 207 [82, 90, 7, 28] True
1 1000.0 986.7 to 1037.3 154 [60, 80, 0, 14] False
2 700.0 683.1 to 733.7 127 [49, 65, 2, 11] True
3 410.0 389.6 to 440.2 123 [42, 67, 9, 5] True
4 990.0 976.6 to 1027.2 193 [88, 68, 4, 33] True
5 420.0 399.7 to 450.3 133 [66, 47, 5, 15] True
6 820.0 804.5 to 855.1 247 [101, 106, 10, 30] True
7 400.0 379.5 to 430.1 160 [42, 94, 2, 22] True
# Variant: quarters split at the image's midlines, not at the midlines of the extent of the footprints' positions.
n_rows, n_cols = M["shape"][2], M["shape"][1] # MATLAB rows and columns of the masks
pos = POS["peak"]
has_pos = ~np.isnan(pos).any(1)
quarter_img = 2 * (pos[:, 0] > (n_rows + 1) / 2) + (pos[:, 1] > (n_cols + 1) / 2)
print(f"image of {n_rows} rows x {n_cols} columns, midlines at row {(n_rows + 1) / 2} and column {(n_cols + 1) / 2} "
f"(the extent's midlines: row {rmid:.1f}, column {cmid:.1f}); footprints without a position: {(~has_pos).sum()}")
print("all footprints per image quarter (upper left, upper right, lower left, lower right):", np.bincount(quarter_img[has_pos], minlength=4).tolist())
for label in ["(a) largest weight", "(b) the code's times"]:
tab, actives = results[(label, "code's position")]
counts = np.array([np.bincount(quarter_img[a_[has_pos[a_]]], minlength=4) for a_ in actives])
print(f"\nEvent times {label}: events with active footprints in all four image quarters {int((counts > 0).all(1).sum())} of 8")
print(pd.DataFrame(counts, columns=["upper_left", "upper_right", "lower_left", "lower_right"]).assign(
component=range(8), event_t_s=tab.event_t_s).to_string(index=False))
image of 608 rows x 808 columns, midlines at row 304.5 and column 404.5 (the extent's midlines: row 309.5, column 425.5); footprints without a position: 0
all footprints per image quarter (upper left, upper right, lower left, lower right): [1001, 1465, 86, 396]
Event times (a) largest weight: events with active footprints in all four image quarters 8 of 8
upper_left upper_right lower_left lower_right component event_t_s
71 110 6 36 0 286.4
96 126 2 26 1 953.0
63 120 4 28 2 791.5
109 164 16 37 3 1150.7
68 89 3 32 4 518.1
127 139 10 32 5 1165.2
107 125 9 39 6 122.9
81 121 7 51 7 80.4
Event times (b) the code's times: events with active footprints in all four image quarters 7 of 8
upper_left upper_right lower_left lower_right component event_t_s
67 110 8 34 0 770.0
50 84 0 17 1 1000.0
37 71 3 14 2 700.0
43 74 11 8 3 410.0
90 78 5 35 4 990.0
63 52 3 16 5 420.0
80 102 10 30 6 820.0
39 93 3 23 7 400.0
Result: blocked (methods_underspecified). Under event reading (a) every component's event has active footprints in all four quarters of the field of view (8 of 8, with 192 to 326 active footprints). Under reading (b), the code's times, 7 of 8 do: component 1's event at 1000 s has 151 active footprints and none in the lower-left quarter, which holds only 85 of all 2948 footprints. The two readings give different outcomes, so by the criterion the claim is methods_underspecified. The code's times fall where this NMF's components have weights of 0.04 to 0.68 of their maximum, as a new NMF run does not reproduce the code's components. With weighted centroids as positions the counts are the same (8 of 8 and 7 of 8), and with the code's frame window (frames 30 times the times) component 1's event again has none in the lower-left quarter (7 of 8). With the quarters split at the image's midlines (row 304.5, column 404.5) in place of the midlines of the footprints' extent, the lower-left quarter holds 86 of 2948 footprints and the counts are again 8 of 8 under (a) and 7 of 8 under (b), component 1 at 1000 s having none there. Under reading (b), five components' assigned times fall where their temporal weights are 0.04 to 0.13 of their maximum (components 1, 2, 3, 5 and 6 are at 0.11, 0.04, 0.06, 0.05 and 0.13), so those times are not events of this NMF run.
ed10-localized-activity¶
Some footprints have activity localized to subregions of the maze.
The paper selected footprints by eye; neither the paper nor the code sets a test for "localized",
so the claim cannot be judged by a method the paper fixes (methods_underspecified). As the
criterion asks, the maps of place_cells.m and each footprint's entropy and map peak are reported.
fast = B["speed"] > np.percentile(B["speed"], 80, method="hazen")
t0 = time.time()
PM = ks.place_measures(T[B["sel"]], B, fast)
print(f"computed in {time.time() - t0:.0f} s; frames above the 80th percentile of speed: {PM['n_frames']}; "
f"position bins with no such frame: {PM['empty_bins']} of 100 (occupied bins hold {PM['occupancy'][PM['occupancy'] > 0].min():.0f} "
f"to {PM['occupancy'].max():.0f} frames)")
e, pk_ = PM["entropy"], PM["peak"]
print(f"entropy of mean activity over linear position, bits (largest possible {PM['max_entropy']:.2f}): percentiles "
f"0/5/50/95/100 {np.round(np.percentile(e, [0, 5, 50, 95, 100]), 3)}")
print(f"peak of the smoothed map of the fraction of time active: percentiles 0/5/50/95/100 {np.round(np.percentile(pk_, [0, 5, 50, 95, 100]), 3)}")
print("Spearman correlation of entropy with peak:", round(pd.Series(e).corr(pd.Series(pk_), method="spearman"), 2))
order = np.argsort(pk_)[::-1]
print("footprints with the 12 highest peaks:", order[:12].tolist(), "; their entropies", np.round(e[order[:12]], 2).tolist())
print("The code's 12 selected footprints (ifig, place_cells.m) index the curated set of 2,382, which the dandiset does not mark, "
"so they cannot be identified here.")
fig, ax = plt.subplots(2, 6, figsize=(15, 5))
for a_, i in zip(ax.ravel(), PM["maps"]):
a_.imshow(PM["maps"][i], cmap="viridis"); a_.set_title(f"footprint {i}: peak {pk_[i]:.2f}", fontsize=8); a_.axis("off")
plt.tight_layout(); plt.show()
print("map values plotted:", round(min(v.min() for v in PM["maps"].values()), 3), "to", round(max(v.max() for v in PM["maps"].values()), 3),
"(each panel on its own colour scale)")
computed in 15 s; frames above the 80th percentile of speed: 6411; position bins with no such frame: 1 of 100 (occupied bins hold 1 to 153 frames) entropy of mean activity over linear position, bits (largest possible 6.63): percentiles 0/5/50/95/100 [3.402 5.814 6.24 6.441 6.541] peak of the smoothed map of the fraction of time active: percentiles 0/5/50/95/100 [0.157 0.359 0.567 0.928 0.957] Spearman correlation of entropy with peak: 0.23 footprints with the 12 highest peaks: [183, 648, 483, 2888, 1414, 2809, 358, 2701, 303, 2470, 262, 1412] ; their entropies [6.34, 5.52, 6.22, 6.46, 6.28, 5.82, 6.34, 6.36, 6.48, 5.86, 6.4, 5.99] The code's 12 selected footprints (ifig, place_cells.m) index the curated set of 2,382, which the dandiset does not mark, so they cannot be identified here.
map values plotted: 0.0 to 0.957 (each panel on its own colour scale)
Result: blocked (methods_underspecified). No test of "localized" is set by the paper or the code, so nothing decides the claim. As the criterion asks, the measures are shown: over 6411 frames above the 80th percentile of speed, the entropy of mean activity over linear position runs from 3.40 to 6.54 bits (median 6.24; largest possible 6.63 over the 99 occupied bins), and the peak of the smoothed map of the fraction of time active from 0.157 to 0.957 (median 0.567). The code's 12 examples index the curated set and cannot be identified here.
print(f"total run time {(time.time() - started) / 60:.1f} min")
total run time 8.8 min
Summary¶
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.
| Claim | Data | Outcome (this attempt) | Outcome of -01 |
|---|---|---|---|
| abstract-thousands-of-neurons | dandiset | reproduced | reproduced |
| abstract-kiloscope-fov | dandiset | reproduced | blocked: metadata insufficient |
| ed7-motion-small (motion notebook) | dandiset | not reproduced | (new) |
| fig2-quadrant-direction-decoding (decoding notebook) | dandiset | blocked: not runnable | (dropped in -01) |
| fig2k-decoder-size (decoding notebook) | dandiset | blocked: not runnable | (new) |
| fig2l-local-decoding-map (decoding notebook) | dandiset | blocked: not runnable | (new) |
| ed9b-octant-decoding (decoding notebook) | dandiset | blocked: not runnable | reproduced |
| ed9c-local-decoding-map (decoding notebook) | dandiset | blocked: not runnable | reproduced |
| ed8-coactive-distributed | dandiset | blocked: methods underspecified | reproduced |
| ed10-localized-activity | dandiset | blocked: methods underspecified | reproduced |
| ed6a-curated-footprints | dandiset | blocked: metadata insufficient | blocked: metadata insufficient |
| ed6bcd-other-mice, ed5-songbirds, abstract-weights, ed1-ed4-resolution, ed3-adc | dandiset | blocked: data missing | blocked: data missing |
| birds-singing-rate | dandiset | blocked: data missing | (new) |
| fig2d-lman-sparse-bursting | dandiset | blocked: methods underspecified (data also missing) | (new) |
| fig1e-featherscope-psf | elsewhere (GitHub) | reproduced | (new) |
| fig1f-kiloscope-psf | elsewhere (GitHub) | partially reproduced | (new) |
| ed1cf-field-curvature | elsewhere (GitHub) | blocked: methods underspecified | (new) |
| ed2b-featherscope-flat-illumination | elsewhere (GitHub) | reproduced | partially reproduced |
| ed2c-kiloscope-flat-illumination | elsewhere (GitHub) | reproduced | partially reproduced |
Where an outcome differs from 000691-scherrer-2023-01, the claim's criterion changed: the claims
file was rewritten from the full text. ed2b and ed2c are now judged against total vignetting, not
against a level of 0.8 chosen by the earlier author; abstract-kiloscope-fov now admits the code's
pixel sizes; ed8 now asks for active footprints in all four quarters under two readings of the
event times; ed10 is now methods_underspecified; ed9b and ed9c now use the paper's null, a
thousand decoders on shuffled activity, whose reading with 200 fits per decoder cannot run here.
How the analysis differs from the paper's.
- All 2948 EXTRACT footprints are used where the claims name the 2,382 curated ones (all 2948 is the claims file's variant); the dandiset does not mark the curated set.
- Position is tracked from the behaviour video with the code's rule; the paper's tracking output is not in the dandiset. Linear position is unwrapped by the loop length, as the claims file says; the code's unwrap differs in about 1% of labels.
- MATLAB's fits are replaced: linear SVMs by scikit-learn (SGD with 10 passes for all footprints, as MATLAB's default for more than 100 predictors), lsqcurvefit by scipy's trust-region-reflective least squares, nnmf by scikit-learn's multiplicative-update NMF, NoRMCorre by a rigid registration written after it.
- The decoding claims' null with 200 fits per shuffled decoder, and the decoder-size and map analyses, need 64,000 to 276,000 fits each and were not run; their cost is measured in the decoding notebooks. Eight more (n2) null decoders per claim, and reduced runs of the decoder-size and map analyses, are in the added notebooks; they decide no outcome.
What was found about the dandiset (with the status of each finding of -01 in the data check): it holds one of four mice and no bird or optical data; the device is described as the weight-optimized microscope; one-photon imaging is stored as two-photon; the behaviour video is described as 30 frames/s but runs at 12.5; there is no pixel size, indicator, location, curation or position; and nothing says whether the movie is raw or motion-corrected. Its values are consistent with the raw video (see the motion notebook), and its shifts are about a quarter of those the paper reports.