import time
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from dandi.dandiapi import DandiAPIClient
from joblib import Parallel, delayed
from scipy import stats
from sklearn.decomposition import NMF

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)
print(f"NWB loaded in {time.time() - started:.0f} s")
video_url = next(u for p, u in urls.items() if p.endswith(".avi"))
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)
# Decoders are fit in parallel; their random numbers are drawn in order first, so results do not
# depend on the number of workers.
track = ks.track_mouse(video_url, keep_frames=SAMPLE_FRAMES)
B = ks.behaviour(S["t_img"], S["t_video"], track)
print(f"behaviour video decoded and tracked in {time.time() - t0:.0f} s")

NWB loaded in 8 s

behaviour video decoded and tracked in 52 s

with DandiAPIClient() as client:
    raw = client.get_dandiset(ks.DANDISET_ID, ks.DANDISET_VERSION).get_raw_metadata()
    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")}
summary = raw["assetsSummary"]
print("Assets (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", "device_description", "imaging_series", "imaging_description",
          "video_description", "video_external_file", "indicator", "location", "excitation_lambda",
          "grid_spacing_present", "virus", "surgery", "subject", "trace_description", "trace_unit",
          "segmentation_type", "cellfind_min_snr", "plane_segmentation_columns", "movie_shape",
          "movie_dtype", "movie_chunks"]:
    print(f"{k}: {m[k]}")

Assets (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.
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
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)

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"])))
n_frames = len(track["x"])
print(f"\nAVI: {n_frames} frames of {track['frame_shape']}, nominal {track['nominal_rate']:.0f} frames/s, "
      f"so {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:", 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"Centroids: 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]})")
um = ks.UM_PER_PIXEL
print(f"At the code's {um} um per pixel (the file gives none): frame {m['movie_shape'][1] * um / 1000:.2f} x "
      f"{m['movie_shape'][2] * um / 1000:.2f} mm; centroid extent {np.ptp(c[:, 0]) * um / 1000:.2f} x {np.ptp(c[:, 1]) * um / 1000:.2f} mm")

                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

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
Centroids: first column 164 to 685 (movie axis 1 has 808 pixels), second column 30 to 587 (movie axis 2 has 608)
At the code's 5.94 um per pixel (the file gives none): frame 4.80 x 3.61 mm; centroid extent 3.09 x 3.31 mm

# Which movie axis each centroid column indexes, and whether traces, masks and movie share a clock.
movie = f["acquisition/TwoPhotonSeries/data"]
masks = f["processing/ophys/ImageSegmentation/PlaneSegmentation/image_mask"]
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]
    win = masks[i, max(c0 - 8, 0):c0 + 9, max(c1 - 8, 0):c1 + 9]
    peak = np.unravel_index(np.argmax(win), win.shape)
    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, "mask_peak_offset": (peak[0] - min(c0, 8), peak[1] - min(c1, 8)),
                 "mask_max": round(float(win.max()), 3), "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 centroid, first 2400 frames, slow drift removed) against traces:")
print(pd.DataFrame(rows).to_string(index=False))

Movie pixels (5 x 5 at each centroid, first 2400 frames, slow drift removed) against traces:
 footprint mask_peak_offset  mask_max  r_own_lag0  best_lag_frames  r_own_best  r_other_footprint_lag0
       188           (0, 0)       1.0       0.601               -1       0.602                   0.052
        69          (0, -1)       1.0       0.696                1       0.696                  -0.061
       577           (0, 0)       1.0       0.749                2       0.753                   0.072
        86           (0, 0)       1.0       0.300                1       0.302                   0.043
        17           (0, 0)       1.0       0.579               -9       0.611                  -0.150
      1627          (0, -1)       1.0       0.869                1       0.869                  -0.074

x, y, area = track["x"], track["y"], track["area"]
skipped = np.arange(len(x)) < ks.FIRST_VIDEO_FRAME
print(f"Frames before frame {ks.FIRST_VIDEO_FRAME} (the code's 64 s at a nominal 30 frames/s), left out: "
      f"{skipped.sum()}, up to t = {S['t_video'][ks.FIRST_VIDEO_FRAME]:.1f} s")
print("Frames from then on with no dark blob in the maze:", int(np.isnan(x[~skipped]).sum()))
print("Blob area (pixels) from then on, percentiles 0/1/50/99/100:", np.percentile(area[~skipped], [0, 1, 50, 99, 100]))
print(f"\nImaging frames overlapping the tracked video: {len(B['sel'])} of {len(S['t_img'])}, "
      f"t = {B['t'][0]:.1f} to {B['t'][-1]:.1f} s")
xmin, xmax, ymin, ymax = B["rect"]
print(f"Loop rectangle (2nd to 98th percentiles of position): x {xmin:.0f} to {xmax:.0f}, y {ymin:.0f} to {ymax:.0f} pixels "
      f"of the {track['frame_shape'][1]} x {track['frame_shape'][0]} frame")
print("Fraction of imaging frames per arm (quadrant 0 to 3):", np.round(np.bincount(B["quadrant"]) / 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; net {abs(B['sunwrap'][-1] - B['sunwrap'][0]) / B['perimeter']:.1f}")
print(f"Largest step between imaging frames: {np.abs(step).max():.1f} pixels (loop length {B['perimeter']:.0f})")
vstep = np.hypot(np.diff(x[~skipped]), np.diff(y[~skipped]))
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]}, {area[i + 1]}; "
          f"frames {i + 1} to {i + 2}: {np.hypot(x[i + 2] - x[i + 1], y[i + 2] - y[i + 1]):.0f} px")
small = (area < 50) & ~skipped
print(f"  frames with blob area below 50 pixels: {small.sum()}; below 100: {((area < 100) & ~skipped).sum()} (of {(~skipped).sum()})")
# The labels again, with the small blobs and both ends of every step above 40 pixels left out.
bad = small.copy()
jump = np.flatnonzero(vstep > 40) + ks.FIRST_VIDEO_FRAME
bad[jump] = bad[jump + 1] = True
clean = dict(track, x=np.where(bad, np.nan, x), y=np.where(bad, np.nan, y))
B_clean = ks.behaviour(S["t_img"], S["t_video"], clean)
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 arm 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")
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"arm {q}")
ax[0].invert_yaxis(); ax[0].set_aspect("equal"); ax[0].set_title("position at imaging frames (pixels)")
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 (pixels)")
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():.0f} to {B['s'].max():.0f}; nothing outside the axes")

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 in the maze: 0
Blob area (pixels) from then on, percentiles 0/1/50/99/100: [  3.  74. 288. 419. 478.]

Imaging frames overlapping the tracked video: 32056 of 36596, t = 153.1 to 1234.4 s
Loop rectangle (2nd to 98th percentiles of position): x 143 to 395, y 119 to 320 pixels of the 640 x 480 frame
Fraction of imaging frames per arm (quadrant 0 to 3): [0.203 0.315 0.232 0.251]
Fraction per octant: [0.087 0.115 0.113 0.202 0.098 0.133 0.129 0.122]
Distance travelled: 48.5 loop lengths; net 2.3
Largest step between imaging frames: 62.2 pixels (loop length 905)

Tracking quality, from frame 1920 on:
  step between consecutive video frames (pixels), percentiles 50/90/99/99.9/100: [  2.2   7.7  13.2  18.2 113.5]
  steps above 20, 40, 80 pixels: [5, 2, 2]
    frame 13388 to 13389 (t = 1064.8 s): 112 px; blob areas 8, 3; frames 13389 to 13390: 113 px
    frame 13389 to 13390 (t = 1064.9 s): 113 px; blob areas 3, 7; frames 13390 to 13391: 6 px
  frames with blob area below 50 pixels: 55; below 100: 358 (of 13601)
  without those 55 video frames: same imaging frames True; imaging frames whose arm changes 0.0005, octant 0.0006, running 0.0107, direction 0.0018
Running frames (speed above median): 16028; direction +1 in 0.514 of them

Plotted: x 133 to 405, y 111 to 328, linear position 10 to 896; nothing outside the axes

# Sample frames of the behaviour video, with the tracked centroid (red cross) and the area searched (white box).
fig, ax = plt.subplots(2, 4, figsize=(14, 5.5))
r0, r1 = ks.ARENA_ROWS
c0, c1 = ks.ARENA_COLS
for a_, i in zip(ax.ravel(), SAMPLE_FRAMES):
    a_.imshow(track["frames"][i])
    a_.plot([c0, c1, c1, c0, c0], [r0, r0, r1, r1, r0], color="w", lw=0.8)
    if np.isfinite(x[i]):
        a_.plot(x[i], y[i], "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")
plt.tight_layout(); plt.show()
print(pd.DataFrame({"frame": SAMPLE_FRAMES, "t_s": [round(S["t_video"][i], 1) for i in SAMPLE_FRAMES],
                    "left_out": [i < ks.FIRST_VIDEO_FRAME for i in SAMPLE_FRAMES],
                    "blob_area_px": [int(area[i]) for i in SAMPLE_FRAMES],
                    "x": np.round(x[list(SAMPLE_FRAMES)], 0), "y": np.round(y[list(SAMPLE_FRAMES)], 0)}).to_string(index=False))
print("Blob area (pixels), percentiles 1/50/99, in the frames left out:", np.percentile(area[skipped], [1, 50, 99]),
      "; from frame 1920 on:", np.percentile(area[~skipped], [1, 50, 99]))

 frame    t_s  left_out  blob_area_px     x     y
     0    0.3      True           304 188.0 316.0
   640   51.1      True           206 387.0 281.0
  1280  102.0      True           344 389.0 219.0
  1900  151.5      True           311 331.0 310.0
  4000  318.4     False           296 151.0 254.0
  8000  636.3     False           340 379.0 168.0
 12000  954.4     False           320 340.0 126.0
 15000 1193.1     False           316 235.0 318.0
Blob area (pixels), percentiles 1/50/99, in the frames left out: [ 56.   320.   516.81] ; from frame 1920 on: [ 74. 288. 419.]

# Shared clock of behaviour and imaging: decode the arm (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)))
    acc = ks.decode_many(act[run], shifted[run], 4, ks.LAMBDA_QUAD_DIR, 3, rng_clock, n_jobs=-1)
    rows.append({"label_shift_s": lag_s, "median_class_accuracy": round(float(np.median(np.nanmedian(acc, 1))), 3)})
clock = pd.DataFrame(rows)
print("Arm decoding (chance 0.25), 3 decoders per shift:")
print(clock.to_string(index=False))

Arm decoding (chance 0.25), 3 decoders per shift:
 label_shift_s  median_class_accuracy
           -60                  0.218
           -20                  0.244
            -5                  0.333
            -2                  0.447
            -1                  0.442
             0                  0.506
             1                  0.549
             2                  0.532
             5                  0.413
            20                  0.203
            60                  0.226

# Footprints deposited twice: pairs with near-identical traces and nearby centroids.
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, 0.95]:
    k = corr > thr
    print(f"pairs with trace correlation > {thr}: {k.sum()}; of them within 10 px: {(k & (dist <= 10)).sum()}; "
          f"median centroid distance {np.median(dist[k]) if k.any() else float('nan'):.0f} px")
print("largest trace correlation between two footprints:", round(float(corr.max()), 3))
dup = (corr > 0.9) & (dist <= 10)
parent = np.arange(T.shape[1])
def find(a):
    while parent[a] != a:
        parent[a] = parent[parent[a]]
        a = parent[a]
    return a
for a, b in zip(iu[0][dup], iu[1][dup]):
    parent[find(a)] = find(b)
n_distinct = len({find(a) for a in range(T.shape[1])})
print("footprints after merging near-duplicates (r > 0.9, within 10 px):", n_distinct)
del Z, C

pairs with trace correlation > 0.7: 1; of them within 10 px: 0; median centroid distance 29 px
pairs with trace correlation > 0.8: 0; of them within 10 px: 0; median centroid distance nan px
pairs with trace correlation > 0.9: 0; of them within 10 px: 0; median centroid distance nan px
pairs with trace correlation > 0.95: 0; of them within 10 px: 0; median centroid distance nan px
largest trace correlation between two footprints: 0.77
footprints after merging near-duplicates (r > 0.9, within 10 px): 2948

nonzero = int((T.max(0) > 0).sum())
arms = int((np.bincount(B["quadrant"], minlength=4) > 0).sum())
print(f"footprints in the file: {T.shape[1]}; with a trace above zero at some time: {nonzero}")
print(f"after merging near-duplicates: {n_distinct}")
print(f"arms of the maze visited during imaging: {arms} of 4; distance run {np.abs(np.diff(B['sunwrap'])).sum() / B['perimeter']:.1f} loop lengths")
events = (T > 0.5 * T.max(0)).sum(0)
print("frames above half its own maximum, per footprint, percentiles 5/50/95:", np.percentile(events, [5, 50, 95]))
print("Variant thresholds: footprints whose maximum exceeds 1, 2, 5 times the median of all maxima:",
      [int((T.max(0) > k * np.median(T.max(0))).sum()) for k in (1, 2, 5)])
print("For reference, the paper's curated count (Extended Data Fig. 6a) is 2382; the dandiset does not mark which footprints it kept.")

footprints in the file: 2948; with a trace above zero at some time: 2948
after merging near-duplicates: 2948
arms of the maze visited during imaging: 4 of 4; distance run 48.5 loop lengths

frames above half its own maximum, per footprint, percentiles 5/50/95: [ 117.35  387.5  1217.3 ]

Variant thresholds: footprints whose maximum exceeds 1, 2, 5 times the median of all maxima: [1474, 298, 3]
For reference, the paper's curated count (Extended Data Fig. 6a) is 2382; the dandiset does not mark which footprints it kept.

rng_qd = np.random.default_rng(691_02)
labels_qd = ks.quad_dir_labels(B)[run]
X_run = act[run]
print("running frames:", len(labels_qd), "; frames per class:", np.bincount(labels_qd, minlength=8).tolist())
t0 = time.time()
real_qd = ks.decode_many(X_run, labels_qd, 8, ks.LAMBDA_QUAD_DIR, 50, rng_qd, n_jobs=-1)
null_qd = ks.decode_many(X_run, labels_qd, 8, ks.LAMBDA_QUAD_DIR, 50, rng_qd, null=True, n_jobs=-1)
print(f"100 decoders in {time.time() - t0:.0f} s")

def decoding_report(real, null, names):
    per_real, per_null = np.nanmedian(real, 1), np.nanmedian(null, 1)
    print("classes absent from a test half (counted out of", real.size, "):", int(np.isnan(real).sum()))
    print(pd.DataFrame({"class": names, "mean_accuracy": np.nanmean(real, 0).round(3),
                        "null_mean": np.nanmean(null, 0).round(3),
                        "null_95th": np.nanpercentile(null, 95, axis=0).round(3)}).to_string(index=False))
    med, q95 = np.median(per_real), np.percentile(per_null, 95)
    print(f"median class accuracy per decoder: median {med:.3f} (range {per_real.min():.3f} to {per_real.max():.3f}); "
          f"null median {np.median(per_null):.3f}, 95th percentile {q95:.3f}; chance {1 / len(names):.3f}")
    print(f"ratio to chance {med * len(names):.2f}; real decoders above the null's 95th percentile: {(per_real > q95).sum()} of {len(per_real)}")
    print(f"code's summary (median of the averaged matrix's diagonal): {np.median(np.nanmean(real, 0)):.3f}")
    return med, q95

names_qd = [f"arm{q}{'+' if d else '-'}" for q in range(4) for d in (0, 1)]
med_qd, q95_qd = decoding_report(real_qd, null_qd, names_qd)
fig, ax = plt.subplots(figsize=(6, 3))
ax.boxplot([np.nanmedian(real_qd, 1), np.nanmedian(null_qd, 1)]); ax.set_xticks([1, 2], ["decoders", "shifted labels"])
ax.axhline(1 / 8, color="k", ls=":", label="chance"); ax.set_ylabel("median class accuracy"); ax.legend()
plt.show()

running frames: 16028 ; frames per class: [1527, 1537, 2882, 2401, 1715, 2056, 1669, 2241]

100 decoders in 69 s
classes absent from a test half (counted out of 400 ): 0
class  mean_accuracy  null_mean  null_95th
arm0-          0.489      0.058      0.176
arm0+          0.464      0.097      0.309
arm1-          0.506      0.194      0.507
arm1+          0.594      0.147      0.425
arm2-          0.239      0.125      0.408
arm2+          0.443      0.133      0.356
arm3-          0.236      0.112      0.334
arm3+          0.459      0.162      0.406
median class accuracy per decoder: median 0.415 (range 0.304 to 0.604); null median 0.087, 95th percentile 0.145; chance 0.125
ratio to chance 3.32; real decoders above the null's 95th percentile: 50 of 50
code's summary (median of the averaged matrix's diagonal): 0.461

# Variant: all frames (no running selection), 20 decoders each.
rng_qd_var = np.random.default_rng(691_03)
labels_all = ks.quad_dir_labels(B)
r_all = ks.decode_many(act, labels_all, 8, ks.LAMBDA_QUAD_DIR, 20, rng_qd_var, n_jobs=-1)
n_all = ks.decode_many(act, labels_all, 8, ks.LAMBDA_QUAD_DIR, 20, rng_qd_var, null=True, n_jobs=-1)
print(f"all frames: median class accuracy {np.median(np.nanmedian(r_all, 1)):.3f}; "
      f"null 95th percentile {np.percentile(np.nanmedian(n_all, 1), 95):.3f}")

all frames: median class accuracy 0.348; null 95th percentile 0.127

# Check of the SGD solver: the same splits fit with an exact solver of the same objective
# (liblinear, up to 10,000 iterations), 3 decoders and 3 with shifted labels.
seeds = np.random.SeedSequence(691_09).spawn(6)
checks = Parallel(n_jobs=-1)(delayed(ks.solver_check)(X_run, labels_qd, 8, ks.LAMBDA_QUAD_DIR, np.random.default_rng(s), null=k >= 3)
                             for k, s in enumerate(seeds))
check_qd = pd.DataFrame(checks).assign(labels=["real"] * 3 + ["shifted"] * 3)
print(check_qd.round(4).to_string(index=False))

 median_acc_sgd  objective_sgd  median_acc_exact  objective_exact  exact_converged  labels
         0.5354         2.3405            0.2615           0.1184            False    real
         0.4427         2.0439            0.2991           0.0897            False    real
         0.4420         2.5547            0.3181           0.1122            False    real
         0.1521         3.2066            0.1309           0.1086            False shifted
         0.0525         3.0650            0.1310           0.0863            False shifted
         0.0980         2.7171            0.1641           0.1017            False shifted

rng_oct = np.random.default_rng(691_04)
labels_oct = B["octant"][run]
print("frames per octant:", np.bincount(labels_oct, minlength=8).tolist())
real_oct = ks.decode_many(X_run, labels_oct, 8, ks.LAMBDA_OCTANT, 50, rng_oct, n_jobs=-1)
null_oct = ks.decode_many(X_run, labels_oct, 8, ks.LAMBDA_OCTANT, 50, rng_oct, null=True, n_jobs=-1)
med_oct, q95_oct = decoding_report(real_oct, null_oct, [f"octant{k}" for k in range(8)])
conf = np.nanmean(real_oct, 0)
fig, ax = plt.subplots(figsize=(6, 3))
ax.bar(range(8), conf); ax.axhline(1 / 8, color="k", ls=":"); ax.set_xlabel("octant"); ax.set_ylabel("mean accuracy")
plt.show()
print(f"Plotted accuracies {conf.min():.3f} to {conf.max():.3f}")

frames per octant: [1285, 1779, 2444, 2839, 1808, 1963, 1974, 1936]

classes absent from a test half (counted out of 400 ): 0
  class  mean_accuracy  null_mean  null_95th
octant0          0.280      0.058      0.210
octant1          0.318      0.103      0.391
octant2          0.356      0.190      0.431
octant3          0.499      0.180      0.379
octant4          0.355      0.112      0.307
octant5          0.260      0.117      0.359
octant6          0.236      0.128      0.331
octant7          0.337      0.116      0.301
median class accuracy per decoder: median 0.313 (range 0.211 to 0.416); null median 0.095, 95th percentile 0.144; chance 0.125
ratio to chance 2.50; real decoders above the null's 95th percentile: 50 of 50
code's summary (median of the averaged matrix's diagonal): 0.328

Plotted accuracies 0.236 to 0.499

# Check of the SGD solver, as for the previous claim: the same splits fit with the exact solver
# (liblinear, up to 10,000 iterations), 3 decoders and 3 with shifted labels.
seeds = np.random.SeedSequence(691_12).spawn(6)
checks = Parallel(n_jobs=-1)(delayed(ks.solver_check)(X_run, labels_oct, 8, ks.LAMBDA_OCTANT, np.random.default_rng(s), null=k >= 3)
                             for k, s in enumerate(seeds))
check_oct = pd.DataFrame(checks).assign(labels=["real"] * 3 + ["shifted"] * 3)
print(check_oct.round(4).to_string(index=False))

 median_acc_sgd  objective_sgd  median_acc_exact  objective_exact  exact_converged  labels
         0.3728         4.1232            0.2917           0.0882            False    real
         0.3012         3.8869            0.2204           0.0889            False    real
         0.3377         4.2872            0.2820           0.0746            False    real
         0.0937         4.5371            0.0939           0.0939            False shifted
         0.0630         6.1753            0.0919           0.0944            False shifted
         0.0768         5.5076            0.1451           0.0708            False shifted

def local_map(n_near, radius_px, spacing, n_dec, seed, solver=None):
    """Grid points with at least n_near centroids within radius_px; at each, n_dec decoders and
    n_dec shifted-label decoders on its n_near nearest footprints, each point with its own
    generator spawned from the seed, run in parallel."""
    cen = S["centroids"].astype(float)
    points = []
    for a in np.arange(cen[:, 0].min(), cen[:, 0].max() + 1, spacing):
        for b in np.arange(cen[:, 1].min(), cen[:, 1].max() + 1, spacing):
            d = np.hypot(cen[:, 0] - a, cen[:, 1] - b)
            if (d <= radius_px).sum() >= n_near:
                points.append((a, b, int((d <= radius_px).sum()), np.argsort(d)[:n_near]))
    seeds = np.random.SeedSequence(seed).spawn(len(points))
    fits = Parallel(n_jobs=-1)(delayed(ks.local_decoders)(X_run[:, idx], labels_oct, ks.LAMBDA_OCTANT, n_dec, s, solver)
                               for (_, _, _, idx), s in zip(points, seeds))
    rows = []
    for (a, b, n_within, _), fit in zip(points, fits):
        real, null = fit["real"], fit["null"]
        z = (real.mean() - null.mean()) / null.std(ddof=1)
        rows.append({"c0": a, "c1": b, "score": real.mean(), "null_mean": null.mean(), "null_sd": null.std(ddof=1),
                     "z": z, "p": stats.norm.sf(z), "n_within": n_within, "at_iteration_limit": fit["at_iteration_limit"],
                     # Variant test: Welch's t-test of the 10 real against the 10 shifted decoders, one-sided.
                     "p_t": stats.ttest_ind(real, null, equal_var=False, alternative="greater").pvalue})
    out = pd.DataFrame(rows)
    out["q_bh"] = stats.false_discovery_control(out.p.to_numpy())
    out["significant"] = out.q_bh < 0.05
    out["significant_t"] = stats.false_discovery_control(out.p_t.to_numpy()) < 0.05
    return out

def map_report(m, label):
    sig = m[m.significant]
    print(f"{label}: {len(m)} grid points used; significant at BH 5%: {len(sig)} ({len(sig) / len(m):.0%})")
    print(f"  score over all points {m.score.min():.3f} to {m.score.max():.3f} (median {m.score.median():.3f}); "
          f"null means {m.null_mean.min():.3f} to {m.null_mean.max():.3f}")
    print(f"  null standard deviation per point, percentiles 5/50/95: {np.round(np.percentile(m.null_sd, [5, 50, 95]), 4)}; "
          f"z percentiles 5/50/95: {np.round(np.percentile(m.z, [5, 50, 95]), 2)}; largest z {m.z.max():.2f}, smallest q {m.q_bh.min():.3f}")
    print(f"  z the best point needs to pass alone (p < 0.05 / {len(m)}): {stats.norm.isf(0.05 / len(m)):.2f}")
    print(f"  exact fits that reached the iteration limit: {m.at_iteration_limit.sum()} of {20 * len(m)}")
    if len(sig):
        print(f"  significant points' scores {sig.score.min():.3f} to {sig.score.max():.3f} (range {sig.score.max() - sig.score.min():.3f}); "
              f"smallest z among them {sig.z.min():.2f}")
    sig_t = m[m.significant_t]
    print(f"  variant, Welch t-test in place of the z-test: significant at BH 5%: {len(sig_t)}"
          + (f"; scores {sig_t.score.min():.3f} to {sig_t.score.max():.3f} (range {sig_t.score.max() - sig_t.score.min():.3f})" if len(sig_t) else ""))

t0 = time.time()
radius = 45       # 270 um at 5.94 um per pixel, rounded to 45 pixels as in the criterion
amap = local_map(32, radius, 30, 10, 691_05)
print(f"computed in {time.time() - t0:.0f} s; radius {radius} px; exact solver (32 predictors)")
map_report(amap, "32 footprints within 270 um")
best = amap.sort_values("score", ascending=False).head(5)
print("best points (pixels):"); print(best.round(3).to_string(index=False))
fig, ax = plt.subplots(figsize=(5, 4.5))
sc = ax.scatter(amap.c1, amap.c0, c=amap.score, s=40, marker="s", cmap="viridis")
ax.scatter(amap.c1[~amap.significant], amap.c0[~amap.significant], s=40, marker="s", facecolors="none", edgecolors="k")
ax.invert_yaxis(); ax.set_aspect("equal"); plt.colorbar(sc, label="score (black outline: not significant)")
ax.set_xlabel("second centroid column (px)"); ax.set_ylabel("first centroid column (px)")
plt.show()
print(f"Plotted scores {amap.score.min():.3f} to {amap.score.max():.3f}; all points inside the axes")

computed in 625 s; radius 45 px; exact solver (32 predictors)
32 footprints within 270 um: 150 grid points used; significant at BH 5%: 19 (13%)
  score over all points 0.068 to 0.187 (median 0.125); null means 0.065 to 0.123
  null standard deviation per point, percentiles 5/50/95: [0.0182 0.03   0.046 ]; z percentiles 5/50/95: [-0.17  1.04  3.07]; largest z 5.81, smallest q 0.000
  z the best point needs to pass alone (p < 0.05 / 150): 3.40
  exact fits that reached the iteration limit: 2264 of 3000
  significant points' scores 0.112 to 0.187 (range 0.075); smallest z among them 2.53
  variant, Welch t-test in place of the z-test: significant at BH 5%: 80; scores 0.108 to 0.187 (range 0.079)
best points (pixels):
   c0    c1  score  null_mean  null_sd     z     p  n_within  at_iteration_limit  p_t  q_bh  significant  significant_t
314.0 240.0  0.187      0.085    0.035 2.896 0.002       256                  17  0.0 0.024         True           True
344.0 180.0  0.186      0.089    0.032 3.048 0.001       187                  12  0.0 0.019         True           True
284.0 210.0  0.184      0.092    0.041 2.248 0.012       203                  18  0.0 0.072        False           True
314.0 180.0  0.177      0.077    0.030 3.297 0.000       115                  14  0.0 0.012         True           True
284.0 180.0  0.174      0.089    0.027 3.090 0.001       136                  13  0.0 0.019         True           True

Plotted scores 0.068 to 0.187; all points inside the axes

# Variants. (a) The solver of the first round of this attempt: SGD with 10 passes on the same
# points and seeds. (b) The code's choice, 64 footprints within 50 pixels, on a 60 pixel grid
# (the criterion's 40 pixel grid puts the notebook over the 30 minute limit on four cores)
# (exact solver, 64 predictors).
t0 = time.time()
amap_sgd = local_map(32, radius, 30, 10, 691_05, solver="sgd")
map_report(amap_sgd, "(a) 32 footprints within 270 um, SGD solver")
amap64 = local_map(64, 50, 60, 10, 691_06)
map_report(amap64, "(b) 64 footprints within 50 px (code), 60 px grid")
print(f"variants computed in {time.time() - t0:.0f} s")

# The two solvers on the same splits at four points of the map (one decoder each, real labels).
cen = S["centroids"].astype(float)
picked_points = amap.iloc[np.linspace(0, len(amap) - 1, 4).astype(int)]
seeds = np.random.SeedSequence(691_11).spawn(4)
rows = []
for (_, pt), s in zip(picked_points.iterrows(), seeds):
    idx = np.argsort(np.hypot(cen[:, 0] - pt.c0, cen[:, 1] - pt.c1))[:32]
    rows.append({"c0": pt.c0, "c1": pt.c1, **ks.solver_check(X_run[:, idx], labels_oct, 8, ks.LAMBDA_OCTANT, np.random.default_rng(s), max_iter=100000)})
print("SGD and exact solver on the same split (exact up to 100,000 iterations); objective = summed mean hinge + penalty:")
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("For reference, predicting every class as absent (w = 0, b = -1) gives an objective of 2 * (sum of class fractions) = 2.0")

(a) 32 footprints within 270 um, SGD solver: 150 grid points used; significant at BH 5%: 1 (1%)
  score over all points 0.021 to 0.158 (median 0.068); null means 0.012 to 0.096
  null standard deviation per point, percentiles 5/50/95: [0.0344 0.0526 0.07  ]; z percentiles 5/50/95: [-0.39  0.33  1.77]; largest z 3.41, smallest q 0.048
  z the best point needs to pass alone (p < 0.05 / 150): 3.40
  exact fits that reached the iteration limit: 0 of 3000
  significant points' scores 0.137 to 0.137 (range 0.000); smallest z among them 3.41
  variant, Welch t-test in place of the z-test: significant at BH 5%: 5; scores 0.106 to 0.158 (range 0.052)

(b) 64 footprints within 50 px (code), 60 px grid: 34 grid points used; significant at BH 5%: 22 (65%)
  score over all points 0.130 to 0.211 (median 0.163); null means 0.098 to 0.131
  null standard deviation per point, percentiles 5/50/95: [0.0168 0.0235 0.0321]; z percentiles 5/50/95: [0.88 2.   3.99]; largest z 5.15, smallest q 0.000
  z the best point needs to pass alone (p < 0.05 / 34): 2.97
  exact fits that reached the iteration limit: 680 of 680
  significant points' scores 0.141 to 0.211 (range 0.070); smallest z among them 1.89
  variant, Welch t-test in place of the z-test: significant at BH 5%: 33; scores 0.130 to 0.211 (range 0.081)
variants computed in 387 s

SGD and exact solver on the same split (exact up to 100,000 iterations); objective = summed mean hinge + penalty:
   c0    c1  median_acc_sgd  objective_sgd  median_acc_exact  objective_exact  exact_converged
164.0 210.0          0.0895        24.6867            0.1745           2.0000             True
374.0 240.0          0.0000        68.9142            0.1156           1.9764            False
554.0 270.0          0.1592        19.8103            0.1428           1.9918            False
674.0 420.0          0.0497        30.9869            0.1141           2.0000             True
For reference, predicting every class as absent (w = 0, b = -1) gives an objective of 2 * (sum of class fractions) = 2.0

rng_pf = np.random.default_rng(691_07)
fast = B["speed"] > np.percentile(B["speed"], 80)
s_run = B["s"][fast]
lo, hi = np.percentile(s_run, [1, 99])
bins = np.clip(np.floor(np.clip((s_run - lo) / (hi - lo), 0, 0.999) * 100), 0, 99).astype(int)
A = T[B["sel"]][fast].astype(np.float64)
n = len(bins)

def entropies(bin_labels):
    onehot = np.zeros((100, n)); onehot[bin_labels, np.arange(n)] = 1
    occ = onehot.sum(1)
    mean_act = (onehot @ A) / np.maximum(occ, 1)[:, None]
    p = mean_act / mean_act.sum(0)
    with np.errstate(divide="ignore", invalid="ignore"):
        return -np.nansum(np.where(p > 0, p * np.log2(p), 0), axis=0), occ

obs, occ = entropies(bins)
print(f"frames above the 80th percentile of speed: {n}; frames per bin {occ.min():.0f} to {occ.max():.0f}")
nulls = np.array([entropies(np.roll(bins, int(rng_pf.integers(n // 10, n - n // 10))))[0] for _ in range(100)])
below = (nulls <= obs).sum(0)
n_fp = A.shape[1]
expected = n_fp / 101
for k, label in [(0, "below all 100 (criterion)"), (5, "at most 5 null values below (p <= 0.06)")]:
    count = int((below <= k).sum())
    fp_rate = (k + 1) / 101
    p_binom = stats.binomtest(count, n_fp, fp_rate, alternative="greater").pvalue
    print(f"{label}: {count} of {n_fp} footprints ({count / n_fp:.1%}); expected by chance {n_fp * fp_rate:.0f}; binomial p = {p_binom:.2g}")
print(f"entropy, bits (maximum log2(100) = {np.log2(100):.2f}): observed median {np.median(obs):.3f}, "
      f"range {obs.min():.3f} to {obs.max():.3f}; null median {np.median(nulls):.3f}")
# Empirical false-positive count: each shifted copy in turn treated as the observed data and
# compared with the other 99 (the same shifts serve every footprint, so counts are correlated).
loo = np.array([int((np.delete(nulls, j, 0) > nulls[j]).all(0).sum()) for j in range(100)])
print(f"shifted copy as observed, below all 99 others: count per copy median {np.median(loo):.0f}, 95th percentile "
      f"{np.percentile(loo, 95):.0f}, maximum {loo.max()} (expected {n_fp / 100:.0f}); copies with a count of at least "
      f"{int((below == 0).sum())}: {(loo >= (below == 0).sum()).sum()} of 100")
gain = (nulls.mean(0) - obs)
print("entropy below the null mean, localized footprints, percentiles 5/50/95 (bits):",
      np.round(np.percentile(gain[below == 0], [5, 50, 95]), 3))
top = np.argsort(gain)[::-1][:6]
fig, ax = plt.subplots(1, 6, figsize=(14, 2), sharey=True)
for a_, i in zip(ax, top):
    prof = np.bincount(bins, weights=A[:, i], minlength=100) / np.maximum(occ, 1)
    a_.plot(prof / prof.max()); a_.set_title(f"footprint {i}", fontsize=8)
ax[0].set_ylabel("activity / max"); fig.supxlabel("linear position bin", fontsize=8)
plt.show()
print("Profiles are each scaled to their maximum (0 to 1).")

frames above the 80th percentile of speed: 6411; frames per bin 0 to 224

below all 100 (criterion): 100 of 2948 footprints (3.4%); expected by chance 29; binomial p = 4.5e-25
at most 5 null values below (p <= 0.06): 448 of 2948 footprints (15.2%); expected by chance 175; binomial p = 1.7e-72
entropy, bits (maximum log2(100) = 6.64): observed median 6.128, range 3.642 to 6.513; null median 6.183
shifted copy as observed, below all 99 others: count per copy median 28, 95th percentile 54, maximum 62 (expected 29); copies with a count of at least 100: 0 of 100
entropy below the null mean, localized footprints, percentiles 5/50/95 (bits): [0.201 0.448 0.899]

Profiles are each scaled to their maximum (0 to 1).

rng_nmf = np.random.default_rng(691_08)
t0 = time.time()
fits = []
for rep in range(4):
    model = NMF(8, solver="mu", init="random", max_iter=200, random_state=int(rng_nmf.integers(2**31)))
    W = model.fit_transform(T)
    fits.append((model.reconstruction_err_, W, model.components_))
err, W, H = min(fits, key=lambda z: z[0])
print(f"NMF in {time.time() - t0:.0f} s; reconstruction errors {[round(z[0], 2) for z in fits]}")
rank = np.argsort(np.argsort(H, axis=1), axis=1)       # rank of each footprint within each component
group = rank.argmax(0)
cen_um = S["centroids"] * ks.UM_PER_PIXEL
def spread(p):
    return np.sqrt(((p - p.mean(0)) ** 2).sum(1).mean())
all_spread = spread(cen_um)
peak_all = T.max(0)
def events_table(threshold, t_from=0.0):
    rows = []
    for c in range(8):
        members = np.flatnonzero(group == c)
        later = S["t_trace"] >= t_from
        tpk = S["t_trace"][later][np.argmax(W[later, c])]
        win = np.abs(S["t_trace"] - tpk) <= 25
        act_ = members[T[win][:, members].max(0) > threshold * peak_all[members]]
        d = np.hypot(*(cen_um[act_][:, None, :] - cen_um[act_][None, :, :]).transpose(2, 0, 1)) if len(act_) > 1 else np.zeros((1, 1))
        rows.append({"component": c, "members": len(members), "event_t_s": round(tpk, 1), "active": len(act_),
                     "spread_ratio": round(spread(cen_um[act_]) / all_spread, 3) if len(act_) > 1 else np.nan,
                     "pairwise_90th_mm": round(np.percentile(d[np.triu_indices_from(d, 1)], 90) / 1000, 2) if len(act_) > 1 else np.nan})
    return pd.DataFrame(rows)
ev = events_table(0.5)
print(f"spread of all footprints: {all_spread / 1000:.2f} mm (RMS distance from their mean, at {ks.UM_PER_PIXEL} um/px)")
print(ev.to_string(index=False))
passing = int(((ev.active >= 10) & (ev.spread_ratio >= 0.5)).sum())
print(f"events with at least 10 active footprints and spread ratio >= 0.5: {passing} of 8")
for thr in (0.3, 0.7):
    e2 = events_table(thr)
    print(f"variant, active above {thr} of maximum: {int(((e2.active >= 10) & (e2.spread_ratio >= 0.5)).sum())} of 8; "
          f"spread ratios {e2.spread_ratio.min():.2f} to {e2.spread_ratio.max():.2f}")
ev_late = events_table(0.5, t_from=B["t"][0])
print(f"variant, event times from t = {B['t'][0]:.0f} s on (after the opening part the behaviour analysis leaves out): "
      f"{int(((ev_late.active >= 10) & (ev_late.spread_ratio >= 0.5)).sum())} of 8; event times {ev_late.event_t_s.tolist()}; "
      f"spread ratios {ev_late.spread_ratio.min():.2f} to {ev_late.spread_ratio.max():.2f}")

# Baselines: members active in randomly placed 50 s windows, and the spread of random sets of footprints.
rng_ed8 = np.random.default_rng(691_10)
t_tr = S["t_trace"]
rows = []
for c in range(8):
    members = np.flatnonzero(group == c)
    fracs = []
    for tc in rng_ed8.uniform(t_tr[0] + 25, t_tr[-1] - 25, 100):
        win = np.abs(t_tr - tc) <= 25
        fracs.append(np.mean(T[win][:, members].max(0) > 0.5 * peak_all[members]))
    n_act = int(ev.active[c])
    rand_spread = [spread(cen_um[rng_ed8.choice(len(cen_um), n_act, replace=False)]) / all_spread for _ in range(200)]
    rows.append({"component": c, "event_fraction_active": round(n_act / ev.members[c], 3),
                 "random_windows_median": round(float(np.median(fracs)), 3), "random_windows_95th": round(float(np.percentile(fracs, 95)), 3),
                 "event_spread_ratio": ev.spread_ratio[c], "random_set_spread_5th": round(float(np.percentile(rand_spread, 5)), 3),
                 "random_set_spread_95th": round(float(np.percentile(rand_spread, 95)), 3)})
print("\nBaselines (100 random 50 s windows per component; 200 random sets of footprints of the event's size):")
print(pd.DataFrame(rows).to_string(index=False))
fig, ax = plt.subplots(2, 4, figsize=(12, 5))
for c, a_ in enumerate(ax.ravel()):
    members = np.flatnonzero(group == c)
    tpk = S["t_trace"][np.argmax(W[:, c])]
    win = np.abs(S["t_trace"] - tpk) <= 25
    act_ = members[T[win][:, members].max(0) > 0.5 * peak_all[members]]
    a_.scatter(S["centroids"][:, 1], S["centroids"][:, 0], s=1, color="0.8")
    a_.scatter(S["centroids"][act_, 1], S["centroids"][act_, 0], s=3, color="r")
    a_.set_title(f"component {c}: {len(act_)} active", fontsize=8); a_.invert_yaxis(); a_.set_aspect("equal"); a_.axis("off")
plt.show()

NMF in 54 s; reconstruction errors [np.float32(97.79), np.float32(97.78), np.float32(97.79), np.float32(97.8)]
spread of all footprints: 0.98 mm (RMS distance from their mean, at 5.94 um/px)
 component  members  event_t_s  active  spread_ratio  pairwise_90th_mm
         0      391      286.4     222         1.042              2.29
         1      397      953.0     257         0.945              2.05
         2      322      791.5     208         1.007              2.18
         3      376     1150.7     320         0.999              2.16
         4      337      518.1     196         1.055              2.27
         5      380     1165.2     305         0.956              2.06
         6      407      122.9     281         0.992              2.13
         7      338       80.4     262         1.050              2.26
events with at least 10 active footprints and spread ratio >= 0.5: 8 of 8
variant, active above 0.3 of maximum: 8 of 8; spread ratios 0.95 to 1.05
variant, active above 0.7 of maximum: 8 of 8; spread ratios 0.92 to 1.07

variant, event times from t = 153 s on (after the opening part the behaviour analysis leaves out): 8 of 8; event times [286.4, 953.0, 791.5, 1150.7, 518.1, 1165.2, 843.6, 535.3]; spread ratios 0.94 to 1.05

Baselines (100 random 50 s windows per component; 200 random sets of footprints of the event's size):
 component  event_fraction_active  random_windows_median  random_windows_95th  event_spread_ratio  random_set_spread_5th  random_set_spread_95th
         0                  0.568                  0.450                0.568               1.042                  0.956                   1.037
         1                  0.647                  0.418                0.645               0.945                  0.962                   1.040
         2                  0.646                  0.399                0.574               1.007                  0.953                   1.040
         3                  0.851                  0.435                0.581               0.999                  0.965                   1.033
         4                  0.582                  0.519                0.632               1.055                  0.951                   1.042
         5                  0.803                  0.420                0.580               0.956                  0.962                   1.034
         6                  0.690                  0.409                0.611               0.992                  0.954                   1.032
         7                  0.775                  0.475                0.705               1.050                  0.963                   1.036

setups = {
    "ed2b-featherscope": dict(folder="sf3/feather", um=1.68, max_um=500, med=3, R=400,
                              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, max_um=2000, med=10, R=1500,
                           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"}),
}
profiles = {}
fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))
for a_, (name, s) in zip(ax, setups.items()):
    profiles[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)
        prof = ks.radial_profile(img, s["um"], s["max_um"], s["med"])
        profiles[name][label] = prof
        a_.plot(prof, label=label)
        print(f"{name} {label}: image {img.shape}, values {img.min():.0f} to {img.max():.0f}")
    a_.set_xlim(0, s["max_um"]); a_.set_ylim(0, 1.1); a_.axvline(s["R"], color="k", ls=":")
    a_.set_xlabel("radius (um)"); a_.set_title(name); a_.legend(fontsize=7)
plt.show()
for name, s in setups.items():
    p = profiles[name]
    R = s["R"]
    rows = []
    for label, prof in p.items():
        rows.append({"profile": label, "min_0_to_R": round(np.nanmin(prof[:R + 1]), 3),
                     **{f"at_{r}um": round(prof[r], 3) for r in (R // 4, R // 2, 3 * R // 4, R)}})
    print(f"\n{name} (R = {R} um):")
    print(pd.DataFrame(rows).to_string(index=False))
    n_out = {label: int(np.sum(prof[:s['max_um']] > 1.1)) for label, prof in p.items()}
    print("radii with values above the plotted 1.1:", n_out)

# Sensitivity of the corrected profiles to the radius and the centre.
for name, s in setups.items():
    img = ks.read_png(ks.code_file_url(f"{s['folder']}/{s['files']['corrected']}"))
    img = ks.unpack_rgb565(img) if img.ndim == 3 else img.astype(float)
    xc, yc = ks.weighted_centre(img)
    h, w = img.shape
    step = round(100 / s["um"])       # 100 um in pixels
    centres = {"intensity-weighted (main)": (xc, yc), "image centre": ((w + 1) / 2, (h + 1) / 2),
               "weighted + 100 um in x": (xc + step, yc), "weighted - 100 um in x": (xc - step, yc),
               "weighted + 100 um in y": (xc, yc + step), "weighted - 100 um in y": (xc, yc - step)}
    rows = []
    for label, centre in centres.items():
        prof = ks.radial_profile(img, s["um"], s["max_um"], s["med"], centre=centre)
        row = {"centre": label, "x_px": round(centre[0], 1), "y_px": round(centre[1], 1)}
        for f_ in (0.95, 1.0, 1.05):
            r = round(f_ * s["R"])
            row[f"min_0_to_{r}um"] = round(np.nanmin(prof[:r + 1]), 3)
        rows.append(row)
    print(f"\n{name}, corrected profile, minimum out to R - 5%, R, R + 5%:")
    print(pd.DataFrame(rows).to_string(index=False))

ed2b-featherscope corrected: image (608, 808), values 23 to 716

ed2b-featherscope optical vignetting: image (608, 808), values 18 to 195

ed2b-featherscope total vignetting: image (608, 808), values 17 to 778

ed2c-kiloscope corrected: image (608, 808), values 23 to 444

ed2c-kiloscope optical vignetting: image (608, 808), values 37 to 489

ed2c-kiloscope total vignetting: image (608, 808), values 27 to 888

ed2b-featherscope (R = 400 um):
           profile  min_0_to_R  at_100um  at_200um  at_300um  at_400um
         corrected       0.742     0.973     0.953     0.901     0.742
optical vignetting       0.628     0.966     0.875     0.762     0.628
  total vignetting       0.342     0.945     0.785     0.556     0.342
radii with values above the plotted 1.1: {'corrected': 0, 'optical vignetting': 0, 'total vignetting': 0}

ed2c-kiloscope (R = 1500 um):
           profile  min_0_to_R  at_375um  at_750um  at_1125um  at_1500um
         corrected       0.797     0.940     0.944      0.922      0.797
optical vignetting       0.389     0.886     0.696      0.525      0.389
  total vignetting       0.253     0.865     0.599      0.386      0.253
radii with values above the plotted 1.1: {'corrected': 0, 'optical vignetting': 0, 'total vignetting': 0}

ed2b-featherscope, corrected profile, minimum out to R - 5%, R, R + 5%:
                   centre  x_px  y_px  min_0_to_380um  min_0_to_400um  min_0_to_420um
intensity-weighted (main) 390.2 324.8           0.785           0.742           0.698
             image centre 404.5 304.5           0.805           0.766           0.723
   weighted + 100 um in x 450.2 324.8           0.692           0.650           0.608
   weighted - 100 um in x 330.2 324.8           0.853           0.811           0.768
   weighted + 100 um in y 390.2 384.8           0.653           0.626           0.601
   weighted - 100 um in y 390.2 264.8           0.899           0.867           0.829

ed2c-kiloscope, corrected profile, minimum out to R - 5%, R, R + 5%:
                   centre  x_px  y_px  min_0_to_1425um  min_0_to_1500um  min_0_to_1575um
intensity-weighted (main) 439.0 299.5            0.846            0.797            0.742
             image centre 404.5 304.5            0.728            0.695            0.654
   weighted + 100 um in x 456.0 299.5            0.880            0.831            0.776
   weighted - 100 um in x 422.0 299.5            0.807            0.765            0.713
   weighted + 100 um in y 439.0 316.5            0.810            0.766            0.716
   weighted - 100 um in y 439.0 282.5            0.845            0.798            0.745

print(f"total run time {(time.time() - started) / 60:.1f} min")

total run time 25.9 min

Claim	Data	Status
abstract-thousands-of-neurons	dandiset	tested
ed6a-curated-footprints	dandiset	blocked: the curation is not in the dandiset
ed6bcd-other-mice	dandiset	blocked: data missing
ed5-songbirds	dandiset	blocked: data missing
abstract-kiloscope-fov	dandiset	blocked: no pixel size in the dandiset
abstract-weights	dandiset	blocked: data missing
ed1-ed4-resolution	dandiset	blocked: data missing
ed3-adc	dandiset	blocked: data missing
ed8-coactive-distributed	dandiset	tested
fig2-quadrant-direction-decoding	dandiset	tested
ed9b-octant-decoding	dandiset	tested
ed9c-local-decoding-map	dandiset	tested
ed10-localized-activity	dandiset	tested
ed2b-featherscope-flat-illumination	elsewhere (GitHub)	tested
ed2c-kiloscope-flat-illumination	elsewhere (GitHub)	tested

Claim	Data	Outcome
abstract-thousands-of-neurons	dandiset	reproduced
fig2-quadrant-direction-decoding	dandiset	reproduced
ed9b-octant-decoding	dandiset	reproduced
ed9c-local-decoding-map	dandiset	reproduced
ed10-localized-activity	dandiset	reproduced
ed8-coactive-distributed	dandiset	reproduced
ed6a-curated-footprints	dandiset	blocked: metadata insufficient
abstract-kiloscope-fov	dandiset	blocked: metadata insufficient
ed6bcd-other-mice, ed5-songbirds, abstract-weights, ed1-ed4-resolution, ed3-adc	dandiset	blocked: data missing
ed2b-featherscope-flat-illumination	elsewhere (GitHub)	partially reproduced
ed2c-kiloscope-flat-illumination	elsewhere (GitHub)	partially reproduced

Reproducing Scherrer, Lynch et al. (2023) from dandiset 000691¶

Data check¶

abstract-thousands-of-neurons¶

fig2-quadrant-direction-decoding¶

ed9b-octant-decoding¶

ed9c-local-decoding-map¶

ed10-localized-activity¶

ed8-coactive-distributed¶

ed2b-featherscope-flat-illumination and ed2c-kiloscope-flat-illumination¶

Summary¶