import hashlib
import io
import json
import os
import re
import time
import urllib.request
import warnings

os.environ.setdefault("OMP_NUM_THREADS", "1")
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import sigprop as sp

warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 200)

started = time.time()
assets = sp.asset_list()
seg = sorted([a for a in assets if "desc-segmentation" in a["path"]], key=lambda a: sp.session_key(a["path"])[0])
img = sorted([a for a in assets if "desc-imaging" in a["path"]], key=lambda a: sp.session_key(a["path"])[0])
print(f"{len(assets)} assets in dandiset {sp.DANDISET_ID} version {sp.DANDISET_VERSION}: "
      f"{len(img)} imaging, {len(seg)} segmentation files")
depth = sp.load_many(sp.load_imaging_depth, [sp.asset_url(a["asset"]) for a in img])
print(f"read frame depths of the imaging files in {time.time() - started:.0f} s")
res = sp.analyse_all([sp.asset_url(a["asset"]) for a in seg])
print(f"loaded and processed {len(res)} recordings in {time.time() - started:.0f} s")
head = sp.head_neurons()
print("head neurons in the paper's code:", len(head))

223 assets in dandiset 001075 version 0.240930.1859: 113 imaging, 110 segmentation files

read frame depths of the imaging files in 170 s

loaded and processed 110 recordings in 629 s
head neurons in the paper's code: 188

sub_img = {sp.session_key(a["path"])[0]: sp.session_key(a["path"])[1] for a in img}
sub_seg = {sp.session_key(a["path"])[0]: sp.session_key(a["path"])[1] for a in seg}
print("Subjects with an imaging file:", len(sub_img), "; with a segmentation file:", len(sub_seg))
print("Subjects with imaging only:", sorted(set(sub_img) - set(sub_seg)))
print("Every subject has one session:", all(sub_img[s] == sub_seg[s] for s in sub_seg))
sizes = pd.Series([a["size"] for a in seg]) / 1e6
print(f"Segmentation file sizes {sizes.min():.1f} to {sizes.max():.1f} MB; imaging files "
      f"{min(a['size'] for a in img) / 1e9:.1f} to {max(a['size'] for a in img) / 1e9:.1f} GB")

# The two published versions: which files differ, and how.
from dandi.dandiapi import DandiAPIClient
client = DandiAPIClient()
old = {a.path: a for a in client.get_dandiset(sp.DANDISET_ID, "0.240920.1434").get_assets()}
new = {a["path"]: a["asset"] for a in assets}
print("\nVersions 0.240920.1434 and 0.240930.1859 list the same paths:", set(old) == set(new))
same_blob = {p: old[p].blob == new[p].blob for p in new}
same_size = {p: old[p].size == new[p].size for p in new}
print("Imaging files: same blob %d, same size %d, of %d" % (sum(same_blob[a["path"]] for a in img), sum(same_size[a["path"]] for a in img), len(img)))
print("Segmentation files: same blob %d, same size %d, of %d" % (sum(same_blob[a["path"]] for a in seg), sum(same_size[a["path"]] for a in seg), len(seg)))
import h5py, remfile
def h5_names(asset):
    with h5py.File(remfile.File(sp.asset_url(asset)), "r") as f:
        names = []
        f.visit(lambda n: names.append(n) if "OptogeneticStimulusTarget" not in n and not n.startswith("specifications") else None)
        return set(names)
p0 = seg[0]["path"]
n_old, n_new = h5_names(old[p0]), h5_names(new[p0])
print(f"\n{p0}: only in 0.240920.1434:", sorted(n_old - n_new))
print(f"{p0}: only in 0.240930.1859:", sorted(n_new - n_old))

Subjects with an imaging file: 113 ; with a segmentation file: 110
Subjects with imaging only: [20, 23, 33]
Every subject has one session: True
Segmentation file sizes 1.9 to 42.1 MB; imaging files 4.4 to 60.5 GB

Versions 0.240920.1434 and 0.240930.1859 list the same paths: True
Imaging files: same blob 0, same size 113, of 113
Segmentation files: same blob 0, same size 0, of 110

sub-0/sub-0_ses-20210824-104000_desc-segmentation_ophys+ogen.nwb: only in 0.240920.1434: ['processing/ophys/GreenSignals/GreenSignal', 'processing/ophys/GreenSignals/GreenSignal/data', 'processing/ophys/GreenSignals/GreenSignal/table_region', 'processing/ophys/GreenSignals/GreenSignal/timestamps', 'processing/ophys/RedSignals/RedSignal', 'processing/ophys/RedSignals/RedSignal/data', 'processing/ophys/RedSignals/RedSignal/table_region', 'processing/ophys/RedSignals/RedSignal/timestamps']
sub-0/sub-0_ses-20210824-104000_desc-segmentation_ophys+ogen.nwb: only in 0.240930.1859: ['intervals/OptogeneticStimulusTable/target_pumpprobe_id', 'processing/ophys/GreenSignals/BaseGreenSignal', 'processing/ophys/GreenSignals/BaseGreenSignal/data', 'processing/ophys/GreenSignals/BaseGreenSignal/table_region', 'processing/ophys/GreenSignals/BaseGreenSignal/timestamps', 'processing/ophys/GreenSignals/InterpolatedGreenSignal', 'processing/ophys/GreenSignals/InterpolatedGreenSignal/data', 'processing/ophys/GreenSignals/InterpolatedGreenSignal/timestamps', 'processing/ophys/RedSignals/BaseRedSignal', 'processing/ophys/RedSignals/BaseRedSignal/data', 'processing/ophys/RedSignals/BaseRedSignal/table_region', 'processing/ophys/RedSignals/BaseRedSignal/timestamps', 'processing/ophys/RedSignals/InterpolatedRedSignal', 'processing/ophys/RedSignals/InterpolatedRedSignal/data', 'processing/ophys/RedSignals/InterpolatedRedSignal/timestamps']

meta = pd.DataFrame([{
    "subject": int(r["subject_id"]), "strain": r["strain"], "genotype": r["genotype"],
    "stage": r["growth_stage"], "age": r["age"], "start": r["session_start"][:19],
    "minutes": (r["t1"] - r["t0"]) / 60, "volumes": r["n_vol"], "rois": r["n_roi"],
    "stims": len(r["stim_start"]), "dt": r["dt"], "dt_max": r["dt_max"],
    "nan_frac": r["nan_frac"]} for r in res])
print(meta[["strain", "genotype", "stage", "age"]].value_counts().to_string())
print("\nSubjects with strain 'none':", meta[meta.strain == "none"].subject.tolist())
print("\nRecording length (min): min %.1f, median %.1f, max %.1f; longer than 40 min: %d, longer than 44 min: %d"
      % (meta.minutes.min(), meta.minutes.median(), meta.minutes.max(), (meta.minutes > 40).sum(), (meta.minutes > 44).sum()))
print("Volume interval (s): median over recordings %.4f (2 volumes/s in the Methods); largest gap %.2f s"
      % (meta.dt.median(), meta.dt_max.max()))
print("Recordings with a gap above 1 s:", meta[meta.dt_max > 1].subject.tolist())
print("Stimulations per recording: min %d, median %d, max %d, total %d"
      % (meta.stims.min(), meta.stims.median(), meta.stims.max(), meta.stims.sum()))
print("Recordings with fewer than 10 stimulations:", meta[meta.stims < 10][["subject", "stims"]].values.tolist())
dep = pd.DataFrame(depth)
dep["subject"] = [sp.session_key(a["path"])[0] for a in img]
print("\nImaging files: frame interval %s s; depth per frame %.1f to %.1f um; plane step median %.3f um (range %.3f to %.3f)"
      % (np.unique(dep.frame_dt.round(4)).tolist(), dep.depth_min.min(), dep.depth_max.max(), dep.plane_step_um.median(),
         dep.plane_step_um.min(), dep.plane_step_um.max()))
m = meta.merge(dep, on="subject")
print("Frames per volume (imaging frames / segmentation volumes): median %.1f"
      % (m.n_frames / m.volumes).median())
print("Indicator stated for the red channel:", sorted(set(dep.red_indicator)))

strain  genotype  stage  age
AML462  WT        L4     P1D    62
                  YA     P1D    46
none    WT        L4     P1D     1
                  YA     P1D     1

Subjects with strain 'none': [24, 25]

Recording length (min): min 3.3, median 31.7, max 48.7; longer than 40 min: 13, longer than 44 min: 3
Volume interval (s): median over recordings 0.5011 (2 volumes/s in the Methods); largest gap 8.41 s
Recordings with a gap above 1 s: [2, 46, 102]
Stimulations per recording: min 3, median 54, max 92, total 5614
Recordings with fewer than 10 stimulations: [[24, 3], [25, 3], [29, 5], [52, 7], [54, 6], [55, 8], [56, 4], [57, 5], [59, 4], [62, 6], [103, 8]]

Imaging files: frame interval [0.02] s; depth per frame 2.5 to 37.5 um; plane step median 1.219 um (range 1.216 to 1.423)
Frames per volume (imaging frames / segmentation volumes): median 25.0
Indicator stated for the red channel: ['GCaMP6s']

st = pd.DataFrame([{
    "subject": int(r["subject_id"]), "dur": d, "power": p, "manual": np.isfinite(m_),
    "tx": x, "ty": y, "tz": z}
    for r in res for d, p, m_, x, y, z in zip(r["stim_stop"] - r["stim_start"], r["stim_power"],
                                               r["stim_manual_roi"], r["target_x"], r["target_y"], r["target_z"])])
print("Stimulation durations (s):", np.unique(st.dur.round(3)).tolist(), "; powers (W):", np.unique(st.power).tolist())
isi = np.concatenate([np.diff(r["stim_start"]) for r in res])
print("Interval between stimulations (s): 5th %.1f, median %.2f, 95th %.1f percentile"
      % tuple(np.percentile(isi, [5, 50, 95])))
print("Stimulations with a manually located target ROI:", int(st.manual.sum()), "of", len(st))
zero = (st.tx == 0) & (st.ty == 0)
print("Stimulations with target position (0, 0):", int(zero.sum()), "; of these with a manual ROI:", int((zero & st.manual).sum()))
print("Target 'depth_in_um' range: %.3f to %.3f (imaging depths %.1f to %.1f um)"
      % (st.tz.min(), st.tz.max(), dep.depth_min.min(), dep.depth_max.max()))
print("Target x range %.0f to %.0f, y %.0f to %.0f pixels; ROI centroid x %.0f to %.0f, y %.0f to %.0f, plane index %.0f to %.0f"
      % (st.tx[~zero].min(), st.tx.max(), st.ty[~zero].min(), st.ty.max(),
         min(r["centroids"][:, 0].min() for r in res), max(r["centroids"][:, 0].max() for r in res),
         min(r["centroids"][:, 1].min() for r in res), max(r["centroids"][:, 1].max() for r in res),
         min(r["centroids"][:, 2].min() for r in res), max(r["centroids"][:, 2].max() for r in res)))
print("Green and red centroids identical in every file:", all(r["centroids_red_equal"] for r in res))

Stimulation durations (s): [0.5] ; powers (W): [0.0012]
Interval between stimulations (s): 5th 30.8, median 30.81, 95th 38.0 percentile
Stimulations with a manually located target ROI: 2815 of 5614
Stimulations with target position (0, 0): 122 ; of these with a manual ROI: 0
Target 'depth_in_um' range: -0.875 to 0.875 (imaging depths 2.5 to 37.5 um)
Target x range 122 to 459, y 118 to 376 pixels; ROI centroid x 0 to 510, y 14 to 504, plane index 0 to 24
Green and red centroids identical in every file: True

units = sorted({r["green_unit"] for r in res})
print("Unit of the green and red series:", units)
neg = [(int(r["subject_id"]), round(r["green_range"][0])) for r in res if r["green_range"][0] < -100]
print("Recordings whose green trace goes below -100 (subject, minimum):", neg)
g = np.array([r["green_range"] for r in res]); rd = np.array([r["red_range"] for r in res])
print("Green (base) per recording: min %.1f, median of medians %.1f, max %.1f" % (g[:, 0].min(), np.median(g[:, 1]), g[:, 2].max()))
print("Red (base) per recording: min %.1f, median of medians %.1f, max %.1f" % (rd[:, 0].min(), np.median(rd[:, 1]), rd[:, 2].max()))
print("Fraction of missing (NaN) values in the base green traces: median %.3f, max %.3f" % (meta.nan_frac.median(), meta.nan_frac.max()))
print("Interpolated traces equal the base traces where the base is defined:", all(r["interp_matches_base"] for r in res))
print("All series of a file share one timestamp array:", all(r["same_timestamps"] for r in res))
print("Photobleaching fits that failed (trace left uncorrected): %d of %d traces"
      % (sum(r["photobleach_failed"] for r in res), sum(r["n_roi"] for r in res)))
pbs = pd.Series(np.concatenate([r["photobleach_status"] for r in res]))
print("  by kind:", pbs.value_counts().to_dict())
neg_recs = {k for k, r in enumerate(res) if r["green_range"][0] < -100}
pb_neg = sum(int(np.sum(res[k]["photobleach_status"] != "ok")) for k in neg_recs)
print("  of the failed fits, in the %d recordings with traces below -100: %d" % (len(neg_recs), pb_neg))
pbr = pd.Series({int(r["subject_id"]): int(r["photobleach_failed"]) for r in res})
print("  recordings with any failed fit: %d; most failures in one recording (subject: count):" % (pbr > 0).sum(),
      pbr.sort_values(ascending=False).head(12).to_dict())

Unit of the green and red series: ['n.a.']
Recordings whose green trace goes below -100 (subject, minimum): [(2, -196646), (3, -2549), (4, -1135), (34, -8622), (35, -4304), (36, -8201), (38, -562), (39, -965), (40, -1489)]
Green (base) per recording: min -196646.1, median of medians 15.0, max 31263.6
Red (base) per recording: min 0.0, median of medians 48.6, max 3074.5
Fraction of missing (NaN) values in the base green traces: median 0.420, max 0.682
Interpolated traces equal the base traces where the base is defined: True
All series of a file share one timestamp array: True
Photobleaching fits that failed (trace left uncorrected): 897 of 14039 traces
  by kind: {'ok': 13142, 'missing': 729, 'near zero': 138, 'failed': 28, 'nonpositive': 2}
  of the failed fits, in the 9 recordings with traces below -100: 4
  recordings with any failed fit: 76; most failures in one recording (subject: count): {88: 77, 92: 67, 73: 55, 79: 41, 45: 33, 75: 33, 85: 33, 69: 28, 70: 24, 93: 23, 48: 22, 89: 21}

# Check of the shortened photobleaching fit: re-run 12 recordings, drawn with a fixed seed, with
# the fit on every volume and no limit on evaluations, as in the code, and compare the outputs.
rng_pb = np.random.default_rng(1075)
pick = sorted(rng_pb.choice(len(res), 12, replace=False).tolist())
full = sp.analyse_all([res[k]["url"] for k in pick], pb_points=None, pb_nfev=None)
rows = []
for k, f in zip(pick, full):
    r = res[k]
    ok = np.isfinite(r["dff"]) & np.isfinite(f["dff"])
    rows.append({"subject": int(r["subject_id"]), "traces": r["n_roi"],
                 "pb_not_ok_short": int(np.sum(r["photobleach_status"] != "ok")),
                 "pb_not_ok_full": int(np.sum(f["photobleach_status"] != "ok")),
                 "dff_pearson": np.corrcoef(r["dff"][ok], f["dff"][ok])[0, 1],
                 "dff_median_abs_diff": np.median(np.abs(r["dff"][ok] - f["dff"][ok])),
                 "sign_agree": np.mean(np.sign(r["dff"][ok]) == np.sign(f["dff"][ok])),
                 "auto_agree": np.mean(r["auto"] == f["auto"]),
                 "responder_calls_agree": np.mean([np.array_equal(np.isin(np.arange(r["n_roi"]), a), np.isin(np.arange(r["n_roi"]), b))
                                                   for a, b in zip(r["responders"], f["responders"])]),
                 "roi_calls_agree": np.mean(np.concatenate([np.isin(np.arange(r["n_roi"]), a) == np.isin(np.arange(r["n_roi"]), b)
                                                             for a, b in zip(r["responders"], f["responders"])])),
                 "responders_short": sum(len(a) for a in r["responders"]), "responders_full": sum(len(b) for b in f["responders"])})
pbc = pd.DataFrame(rows)
print(pbc.round(4).to_string(index=False))
print("Over the 12 recordings: failed fits %d (short) and %d (full); median Pearson of <dF/F0>_t %.4f; sign agreement %.4f; "
      "autoresponse agreement %.4f; classifier calls per (stimulation, ROI) agreeing %.4f; responses detected %d (short) and %d (full)"
      % (pbc.pb_not_ok_short.sum(), pbc.pb_not_ok_full.sum(), pbc.dff_pearson.median(),
         np.average(pbc.sign_agree, weights=pbc.traces), pbc.auto_agree.mean(), np.average(pbc.roi_calls_agree, weights=pbc.traces),
         pbc.responders_short.sum(), pbc.responders_full.sum()))

 subject  traces  pb_not_ok_short  pb_not_ok_full  dff_pearson  dff_median_abs_diff  sign_agree  auto_agree  responder_calls_agree  roi_calls_agree  responders_short  responders_full
      13      90                0               0       0.9958               0.0008      0.9826      1.0000                 0.6769           0.9952               835              825
      22     169                4               3       0.9846               0.0005      0.9852      1.0000                 0.6533           0.9962              1212             1206
      27      91                1               1       0.9963               0.0010      0.9868      0.9538                 0.5385           0.9902              1203             1183
      29      80                1               1       1.0000               0.0000      1.0000      1.0000                 1.0000           1.0000                98               98
      40     127                2               2       1.0000               0.0004      0.9931      1.0000                 0.7231           0.9973              1749             1747
      59      83                2               2       1.0000               0.0000      1.0000      1.0000                 1.0000           1.0000                48               48
      81     135                4               5       0.9979               0.0008      0.9900      1.0000                 0.4783           0.9950               877              876
      89     140               21              21       0.9969               0.0007      0.9881      1.0000                 0.3750           0.9938               629              630
      91     147                7               6       0.9215               0.0006      0.9909      0.9808                 0.5769           0.9945              1088             1100
      96     119                1               1       0.9932               0.0008      0.9845      1.0000                 0.4878           0.9947               445              427
     102     123                9               5       0.9992               0.0007      0.9887      1.0000                 0.2642           0.9903               896              869
     112     110                7               7       0.9993               0.0008      0.9907      1.0000                 0.3667           0.9911              1081             1094
Over the 12 recordings: failed fits 59 (short) and 54 (full); median Pearson of <dF/F0>_t 0.9974; sign agreement 0.9896; autoresponse agreement 0.9946; classifier calls per (stimulation, ROI) agreeing 0.9947; responses detected 10161 (short) and 10103 (full)

lab = []
for r in res:
    names = [sp.atlas_name(x) for x in r["roi_label"]]
    named = [n for n in names if n]
    in_head = [n for n in named if n in head]
    dup = pd.Series(in_head).value_counts()
    lab.append({"subject": int(r["subject_id"]), "rois": r["n_roi"], "labelled": len(named), "head_labelled": len(in_head),
                "head_distinct": len(set(in_head)), "dup_labels": int((dup > 1).sum()),
                "other_ids": ",".join(r["neuropal_ids_other"])})
lab = pd.DataFrame(lab)
print("ROIs per recording: min %d, median %d, max %d" % (lab.rois.min(), lab.rois.median(), lab.rois.max()))
print("ROIs with a NeuroPAL label: median %d per recording; with a head-neuron label: median %d (%d distinct)"
      % (lab.labelled.median(), lab.head_labelled.median(), lab.head_distinct.median()))
print("Recordings in which a head label is used for more than one ROI:", int((lab.dup_labels > 0).sum()),
      "; such labels per recording, max", lab.dup_labels.max())
print("Head neurons labelled in at least one recording:", len({sp.atlas_name(x) for r in res for x in r["roi_label"]} & set(head)), "of", len(head))
print("Recordings whose neuropal_ids hold text other than a number or blank:", int((lab.other_ids != "").sum()),
      "; values:", sorted({v for s in lab.other_ids for v in s.split(",") if v}))
conf = np.concatenate([r["roi_conf"][np.array([bool(x) for x in r["roi_label"]])] for r in res])
print("Label confidences of labelled ROIs:", pd.Series(conf).round(2).value_counts().sort_index().to_dict())
print("Recordings with ROIs labelled 'target':", sum(any(x == "target" for x in r["roi_label"]) for r in res),
      "; such ROIs:", sum(int(np.sum(r["roi_label"] == "target")) for r in res))
never = sorted(set(head) - {sp.atlas_name(x) for r in res for x in r["roi_label"]})
print("Head neurons never labelled:", never)
print("ROIs labelled RID: %d, in %d recordings" % (sum(int(np.sum(r["roi_label"] == "RID")) for r in res),
                                               sum(any(x == "RID" for x in r["roi_label"]) for r in res)))
man_rid = sum(int(np.sum([r["roi_label"][int(m)] == "RID" for m in r["stim_manual_roi"] if np.isfinite(m)])) for r in res)
print("Manually located targets whose ROI is labelled RID:", man_rid)
non_head = pd.Series([sp.atlas_name(x) for r in res for x in r["roi_label"] if x and sp.atlas_name(x) not in head]).value_counts()
print("Most frequent labels that are not head-neuron names:", non_head.head(15).to_dict())

ROIs per recording: min 80, median 125, max 281
ROIs with a NeuroPAL label: median 67 per recording; with a head-neuron label: median 60 (60 distinct)
Recordings in which a head label is used for more than one ROI: 1 ; such labels per recording, max 1
Head neurons labelled in at least one recording: 186 of 188
Recordings whose neuropal_ids hold text other than a number or blank: 63 ; values: ['I4', 'RIAR', 'check', 'merge']
Label confidences of labelled ROIs: {-1.0: 356, 0.4: 2, 0.5: 16, 0.6: 121, 0.7: 566, 0.8: 1148, 0.9: 394, 1.0: 4849, 11.0: 1}
Recordings with ROIs labelled 'target': 21 ; such ROIs: 49
Head neurons never labelled: ['AIAR', 'I5']
ROIs labelled RID: 11, in 11 recordings
Manually located targets whose ROI is labelled RID: 5
Most frequent labels that are not head-neuron names: {'RIG': 51, 'target': 49, 'AMso': 34, 'RMDD': 30, 'IL1V': 26, 'CEPV': 25, 'SMDD': 21, 'SAAD': 19, 'IL2D': 18, 'AMSoL': 17, 'OLQV': 16, 'IL2V': 15, 'URX': 15, 'AMSoR': 15, 'SMB': 14}

tg = pd.DataFrame([{"subject": int(r["subject_id"]), "z_slope": r["z_slope"], "z_offset": r["z_offset"],
                    "med_dist_px": np.median(r["stim_dist_px"][r["stim_how"] == "nearest"]) if (r["stim_how"] == "nearest").any() else np.nan,
                    "nearest": int((r["stim_how"] == "nearest").sum()), "manual": int((r["stim_how"] == "manual").sum()),
                    "none": int((r["stim_how"] == "no target").sum()),
                    "xy_differs": int(np.sum((r["stim_roi"] != r["stim_roi_xy"]) & (r["stim_how"] == "nearest")))} for r in res])
print("Stimulated ROI assigned: manual %d, nearest centroid %d, no target %d"
      % (tg.manual.sum(), tg.nearest.sum(), tg.none.sum()))
d_all = np.concatenate([r["stim_dist_px"][r["stim_how"] == "nearest"] for r in res])
print("Distance from target to nearest centroid (pixels, 0.42 um each): median %.1f, 90th %.1f, 99th %.1f percentile"
      % tuple(np.percentile(d_all, [50, 90, 99])))
print("Nearest-centroid assignments that change when z is ignored: %d of %d" % (tg.xy_differs.sum(), tg.nearest.sum()))
print("Fitted z map slope (planes per unit): 5th %.1f, median %.1f, 95th %.1f percentile" % tuple(np.percentile(tg.z_slope, [5, 50, 95])))
man_check = []
for r in res:
    ok = (r["stim_how"] == "manual") & ((r["target_x"] > 0) | (r["target_y"] > 0))
    if ok.any():
        man_check.append(np.mean(r["stim_roi"][ok] == r["stim_roi_xy"][ok]))
print("Among manual targets with a stored position, the share whose nearest xy centroid is the manual ROI: median over recordings %.2f"
      % np.median(man_check))

Stimulated ROI assigned: manual 2815, nearest centroid 2677, no target 122
Distance from target to nearest centroid (pixels, 0.42 um each): median 4.1, 90th 8.2, 99th 19.1 percentile
Nearest-centroid assignments that change when z is ignored: 460 of 2677
Fitted z map slope (planes per unit): 5th -13.8, median 0.0, 95th 16.0 percentile
Among manual targets with a stored position, the share whose nearest xy centroid is the manual ROI: median over recordings 0.60

clock = pd.DataFrame([{"subject": int(r["subject_id"]), "stims": len(r["auto"]),
                       "auto_true": r["auto"].mean(), "auto_shift15": r["auto_shift"].mean(),
                       "auto_xy": r["auto_xy"].mean()} for r in res])
a_true = np.concatenate([r["auto"] for r in res]); a_shift = np.concatenate([r["auto_shift"] for r in res])
print("Autoresponse rate, true stimulation times: %.3f; times shifted by +15 s: %.3f (of %d stimulations)"
      % (a_true.mean(), a_shift.mean(), len(a_true)))
print("Autoresponse rate with the nearest centroid in x and y only: %.3f" % np.concatenate([r["auto_xy"] for r in res]).mean())
diff = clock.auto_true - clock.auto_shift15
print("Recordings where the true rate exceeds the shifted rate: %d of %d" % ((diff > 0).sum(), len(clock)))
print("Recordings where it does not (clock not established):")
print(clock[diff <= 0].to_string(index=False))
dff_t = np.concatenate([r["dff"][np.arange(len(r["auto"])), np.clip(r["stim_roi"], 0, None)][r["stim_roi"] >= 0] for r in res])
dff_s = np.concatenate([r["dff_shift"][np.arange(len(r["auto"])), np.clip(r["stim_roi"], 0, None)][r["stim_roi"] >= 0] for r in res])
print("Stimulated ROI <dF/F0>_t: median %.3f at true times, %.3f at shifted times" % (np.nanmedian(dff_t), np.nanmedian(dff_s)))

Autoresponse rate, true stimulation times: 0.552; times shifted by +15 s: 0.103 (of 5614 stimulations)
Autoresponse rate with the nearest centroid in x and y only: 0.490
Recordings where the true rate exceeds the shifted rate: 110 of 110
Recordings where it does not (clock not established):
Empty DataFrame
Columns: [subject, stims, auto_true, auto_shift15, auto_xy]
Index: []
Stimulated ROI <dF/F0>_t: median 0.212 at true times, 0.035 at shifted times

fp = []
for k, r in enumerate(res):
    fp.append((r["n_vol"], r["n_roi"], tuple(np.round(r["stim_start"][:5] - r["t0"], 2)),
               tuple(np.round(r["centroids"][:3].ravel(), 1))))
pairs = [(int(res[a]["subject_id"]), int(res[b]["subject_id"]), what)
         for a in range(len(fp)) for b in range(a + 1, len(fp))
         for what in ["volumes and ROIs", "stimulation times", "centroids"]
         if (what == "volumes and ROIs" and fp[a][:2] == fp[b][:2])
         or (what == "stimulation times" and len(fp[a][2]) > 2 and fp[a][2] == fp[b][2])
         or (what == "centroids" and fp[a][3] == fp[b][3])]
print("Pairs of recordings sharing a fingerprint:", pairs if pairs else "none")
starts = pd.Series([r["session_start"] for r in res])
print("Session start times repeated:", starts[starts.duplicated()].tolist())

Pairs of recordings sharing a fingerprint: [(3, 109, 'volumes and ROIs')]
Session start times repeated: []

obs = sp.observations(res, head)
print("Observation rows:", len(obs), "from", obs.rec.nunique(), "recordings")
print("Stimulations with an assigned target ROI:", obs.groupby(["rec", "stim"]).ngroups)
st_auto = obs.groupby(["rec", "stim"]).agg(auto=("auto", "first"), j=("j", "first"), j_head=("j_head", "first"))
print("Of these with an autoresponse: %d; with a head-neuron label on the target: %d; both: %d"
      % (st_auto.auto.sum(), st_auto.j_head.sum(), (st_auto.auto & st_auto.j_head).sum()))
print("Rows with F0 <= 0: %d (%.4f); with F0 <= 0 among head pairs with an autoresponse: %d"
      % ((obs.f0 <= 0).sum(), (obs.f0 <= 0).mean(), ((obs.f0 <= 0) & obs.i_head & obs.j_head & obs.auto).sum()))
print("Rows from the %d recordings with traces below -100: %d; rows whose trace was not corrected for photobleaching: %d"
      % (len(neg_recs), obs.rec.isin(neg_recs).sum(), (~obs.pb_ok).sum()))
big = obs.dff.abs() > 2
print("Rows beyond +-2: %d; of these with F0 <= 0: %d, from the recordings with negative traces: %d, without photobleaching correction: %d"
      % (big.sum(), (big & (obs.f0 <= 0)).sum(), (big & obs.rec.isin(neg_recs)).sum(), (big & ~obs.pb_ok).sum()))
print("<dF/F0>_t of usable observed traces: 1st %.3f, median %.4f, 99th %.3f percentile; beyond +-2: %d rows"
      % (np.nanpercentile(obs.dff, 1), np.nanmedian(obs.dff), np.nanpercentile(obs.dff, 99), int((obs.dff.abs() > 2).sum())))

Observation rows: 433216 from 110 recordings
Stimulations with an assigned target ROI: 5492
Of these with an autoresponse: 3101; with a head-neuron label on the target: 3509; both: 2034
Rows with F0 <= 0: 3172 (0.0073); with F0 <= 0 among head pairs with an autoresponse: 761
Rows from the 9 recordings with traces below -100: 49753; rows whose trace was not corrected for photobleaching: 2039
Rows beyond +-2: 178; of these with F0 <= 0: 20, from the recordings with negative traces: 42, without photobleaching correction: 78
<dF/F0>_t of usable observed traces: 1st -0.304, median -0.0012, 99th 0.552 percentile; beyond +-2: 178 rows

n_both = sum(1 for r in res if r["n_vol"] > 0 and len(r["stim_start"]) > 0)
print("Subjects in the dandiset:", len(sub_img), "; recordings with traces and a stimulation table:", n_both)
print("WT genotype among them:", int((meta.genotype == "WT").sum()))
print("Median interval between stimulations: %.2f s (criterion: within 10%% of 30 s, 27 to 33)" % np.median(isi))
share = (meta.minutes <= 44).mean()
print("Recordings lasting at most 44 min: %d of %d (%.3f; criterion at least 0.9)" % ((meta.minutes <= 44).sum(), len(meta), share))
print("Stored stimulus durations:", np.unique(st.dur.round(3)).tolist())
print("Recording length per animal (min):", meta.sort_values("minutes").minutes.round(1).tolist())

Subjects in the dandiset: 113 ; recordings with traces and a stimulation table: 110
WT genotype among them: 110
Median interval between stimulations: 30.81 s (criterion: within 10% of 30 s, 27 to 33)
Recordings lasting at most 44 min: 107 of 110 (0.973; criterion at least 0.9)
Stored stimulus durations: [0.5]
Recording length per animal (min): [3.3, 3.5, 5.6, 5.8, 6.3, 6.5, 6.9, 7.4, 7.6, 7.8, 9.8, 10.3, 11.1, 11.7, 13.3, 13.6, 15.0, 16.6, 20.9, 21.8, 22.5, 22.8, 22.9, 23.2, 23.8, 24.1, 24.3, 25.5, 25.5, 25.9, 26.1, 26.2, 26.2, 26.4, 26.5, 26.6, 27.0, 27.4, 27.6, 27.7, 27.7, 27.8, 27.9, 28.0, 28.0, 28.1, 28.4, 28.5, 29.9, 30.7, 30.8, 31.2, 31.5, 31.6, 31.7, 31.7, 31.7, 31.8, 31.9, 32.0, 32.1, 32.1, 32.2, 32.3, 32.9, 33.3, 33.4, 33.4, 33.6, 33.7, 33.7, 33.9, 34.0, 34.3, 34.5, 34.6, 34.7, 34.7, 34.8, 35.1, 35.1, 36.7, 36.7, 36.7, 36.8, 37.1, 37.1, 37.2, 37.2, 37.2, 37.3, 37.4, 37.7, 37.8, 38.5, 38.9, 39.4, 40.3, 41.5, 42.1, 42.1, 42.3, 42.4, 42.4, 42.5, 42.5, 43.1, 45.0, 45.2, 48.7]

def pair_counts(o, require_auto=True, diag=True):
    h = o[o.i_head & o.j_head]
    if require_auto:
        h = h[h.auto]
    if not diag:
        h = h[h.i != h.j]
    n = h.groupby(["i", "j"]).size()
    return n, h

rows = []
for name, ra, dg in [("(a) with i = j, autoresponse", True, True), ("(b) i != j, autoresponse", True, False),
                     ("(c) with i = j, no autoresponse required", False, True)]:
    n, h = pair_counts(obs, ra, dg)
    rows.append({"reading": name, "pairs": len(n), "share_of_188^2": len(n) / 188 ** 2,
                 "neurons_recorded": h[h.i != h.j].i.nunique(), "max_obs": int(n.max()),
                 "animals": h.rec.nunique()})
pm = pd.DataFrame(rows)
print(pm.to_string(index=False))
print("Expected: 23,433 pairs (tolerance 2,343), 186 neurons (18.6), 59 observations (5.9)")
n_a, h_a = pair_counts(obs)
print("\nMost observed pairs, reading (a):", n_a.sort_values(ascending=False).head(5).to_dict())
per_animal = h_a.groupby("rec").apply(lambda x: x.groupby(["i", "j"]).ngroups)
print("Pairs contributed per animal (reading a): min %d, median %d, max %d; animals contributing %d"
      % (per_animal.min(), per_animal.median(), per_animal.max(), len(per_animal)))
stim_n = h_a.groupby("j").apply(lambda x: x.groupby(["rec", "stim"]).ngroups)
print("Head neurons stimulated with an autoresponse: %d (Source Data Extended Data Fig. 3b lists 168 stimulated)" % len(stim_n))

                                 reading  pairs  share_of_188^2  neurons_recorded  max_obs  animals
            (a) with i = j, autoresponse  20587        0.582475               184       50      107
                (b) i != j, autoresponse  20422        0.577807               184       50      107
(c) with i = j, no autoresponse required  23813        0.673749               184       86      108
Expected: 23,433 pairs (tolerance 2,343), 186 neurons (18.6), 59 observations (5.9)

Most observed pairs, reading (a): {('M3R', np.str_('RMDVR')): 50, ('RMDVR', np.str_('RMDVR')): 49, ('AVEL', np.str_('AVEL')): 49, ('I1R', np.str_('RMDVR')): 48, ('M3R', np.str_('AVEL')): 48}
Pairs contributed per animal (reading a): min 22, median 696, max 2281; animals contributing 107

Head neurons stimulated with an autoresponse: 172 (Source Data Extended Data Fig. 3b lists 168 stimulated)

SOURCE = "https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-023-06683-4/MediaObjects/"
def source_sheet(fname, sheet):
    with urllib.request.urlopen(SOURCE + fname, timeout=120) as r:
        raw = r.read()
    print(f"{fname}: {len(raw)} bytes, md5 {hashlib.md5(raw).hexdigest()}")
    x = pd.read_excel(io.BytesIO(raw), sheet_name=sheet, header=None)
    m = x.iloc[1:, 1:].astype(float)
    m.index = x.iloc[1:, 0].astype(str).values
    m.columns = x.iloc[0, 1:].astype(str).values
    return m.stack(future_stack=True).rename_axis(["i", "j"])

sd_obs = source_sheet("41586_2023_6683_MOESM12_ESM.xlsx", "EDFig4a")
sd_q = source_sheet("41586_2023_6683_MOESM12_ESM.xlsx", "EDFig4b")
sd_dff = source_sheet("41586_2023_6683_MOESM11_ESM.xlsx", "EDFig3a")
sd = pd.DataFrame({"n_paper": sd_obs, "q": sd_q, "dff_paper": sd_dff})
sd = sd[sd.n_paper > 0]
print("Source Data: pairs observed %d, with q %d, q < 0.05 %d, largest count %d"
      % (len(sd), sd.q.notna().sum(), (sd.q < 0.05).sum(), sd.n_paper.max()))

def link_stats(here, label):
    j = sd.join(here, how="left")
    found = j.n.notna()
    both = j[found]
    rho_n = stats.spearmanr(both.n_paper, both.n).statistic
    ok = both.dff_paper.notna() & both["mean"].notna()
    r_dff = stats.pearsonr(both.dff_paper[ok], both["mean"][ok]).statistic
    rho_dff = stats.spearmanr(both.dff_paper[ok], both["mean"][ok]).statistic
    extra = len(here.index.difference(sd.index))
    return {"reading": label, "paper_pairs_found": found.mean(), "pairs_here_not_in_paper": extra,
            "spearman_counts": rho_n, "pearson_dff": r_dff, "spearman_dff": rho_dff,
            "median_count_ratio": np.median(both.n / both.n_paper), "pairs_compared": int(ok.sum())}

41586_2023_6683_MOESM12_ESM.xlsx: 378622 bytes, md5 8d55d8113826b809c9610e46accdf093

41586_2023_6683_MOESM12_ESM.xlsx: 378622 bytes, md5 8d55d8113826b809c9610e46accdf093

41586_2023_6683_MOESM11_ESM.xlsx: 282522 bytes, md5 02ca4ec308dc5e398f059c4f9a2c63d1
Source Data: pairs observed 23433, with q 23433, q < 0.05 1310, largest count 59

# The paper's export on OSF (doi 10.17605/OSF.IO/E2SYT), file exported_data.tar.gz (OSF guid 9mecf,
# version 1, md5 7ae518fea414f664e2493238afe209bb). It is streamed once and only its small per-recording
# files are kept: the folder name, the stimulated ROI of every stimulation and the ROI labels.
import tarfile
OSF_URL, OSF_MD5 = "https://osf.io/download/9mecf/", "7ae518fea414f664e2493238afe209bb"

class _Hashed:
    def __init__(self, f):
        self.f, self.md5, self.n = f, hashlib.md5(), 0
    def read(self, n=-1):
        b = self.f.read(n)
        self.md5.update(b)
        self.n += len(b)
        return b

import zlib
texts, big_txt = {}, {}
def _has_rid(k):
    if (k, "labels") not in texts or (k, "stim_neurons") not in texts:
        return None
    le = [x.strip() for x in texts[(k, "labels")].split("\n")]
    return any(0 <= int(s) < len(le) and le[int(s)] == "RID" for s in texts[(k, "stim_neurons")].split())
with urllib.request.urlopen(OSF_URL, timeout=600) as r:
    hr = _Hashed(r)
    with tarfile.open(fileobj=hr, mode="r|gz") as tf:
        for m in tf:
            mm = re.search(r"/(\d+)_(ds_name|stim_neurons|labels|gcamp|t|stim_volume_i)\.txt$", m.name)
            if mm:
                k, what = int(mm.group(1)), mm.group(2)
                if what in ("gcamp", "t", "stim_volume_i"):
                    # the traces are kept (compressed) only for recordings that stimulate an ROI labelled RID
                    if _has_rid(k) is not False:
                        big_txt[(k, what)] = zlib.compress(tf.extractfile(m).read())
                else:
                    texts[(k, what)] = tf.extractfile(m).read().decode()
                    if _has_rid(k) is False:
                        for w in ("gcamp", "t", "stim_volume_i"):
                            big_txt.pop((k, w), None)
    while hr.read(1 << 20):
        pass
print("exported_data.tar.gz: %d bytes, md5 %s, matches the OSF record: %s" % (hr.n, hr.md5.hexdigest(), hr.md5.hexdigest() == OSF_MD5))
n_exp = len({k for k, _ in texts})
exp = []
for k in range(n_exp):
    folder = texts[(k, "ds_name")].strip()
    m = re.search(r"pumpprobe_(\d{8})_(\d{6})", folder)
    exp.append({"k": k, "key": f"{m.group(1)}-{m.group(2)}",
                "stim": np.array([int(x) for x in texts[(k, "stim_neurons")].split()]),
                "labels": texts[(k, "labels")].split("\n")})
exp = {e["key"]: e for e in exp}
print("Exported recordings with distinct session keys:", len(exp))
rid_k = sorted(k for k in range(n_exp) if _has_rid(k))
print("Exported WT recordings that stimulate an ROI labelled RID:", len(rid_k), "; traces kept for", len({k for k, _ in big_txt}))
osf_roi = []
for r in res:
    key = r["session_start"][:10].replace("-", "") + "-" + r["session_start"][11:19].replace(":", "")
    e = exp.get(key)
    osf_roi.append(e["stim"] if e is not None and len(e["stim"]) == len(r["stim_start"]) else None)
excluded = {(k, e) for k, roi in enumerate(osf_roi) if roi is not None for e in np.where(roi == -2)[0]}
print("Stimulations the export marks excluded (-2), over the recordings here:", len(excluded))

exported_data.tar.gz: 523093816 bytes, md5 7ae518fea414f664e2493238afe209bb, matches the OSF record: True
Exported recordings with distinct session keys: 113
Exported WT recordings that stimulate an ROI labelled RID: 7 ; traces kept for 7
Stimulations the export marks excluded (-2), over the recordings here: 950

rows, assign = [], []
for r in res:
    key = r["session_start"][:10].replace("-", "") + "-" + r["session_start"][11:19].replace(":", "")
    e = exp.get(key)
    row = {"subject": int(r["subject_id"]), "key": key, "matched": e is not None}
    if e is not None:
        row["same_stims"] = len(e["stim"]) == len(r["stim_start"])
        labs_e = e["labels"][:-1] if e["labels"] and e["labels"][-1] == "" else e["labels"]
        n_lab = len(labs_e)
        row["rois_here"] = r["n_roi"]
        row["labels_exported"] = n_lab
        ok = row["same_stims"]
        if ok:
            both = (e["stim"] >= 0) & (r["stim_roi"] >= 0)
            row["both_assigned"] = int(both.sum())
            row["agree"] = int(np.sum(e["stim"][both] == r["stim_roi"][both]))
            row["agree_nearest"] = int(np.sum((e["stim"] == r["stim_roi"])[both & (r["stim_how"] == "nearest")]))
            row["n_nearest"] = int(np.sum(both & (r["stim_how"] == "nearest")))
            row["agree_manual"] = int(np.sum((e["stim"] == r["stim_roi"])[both & (r["stim_how"] == "manual")]))
            row["n_manual"] = int(np.sum(both & (r["stim_how"] == "manual")))
            row["exported_excluded"] = int(np.sum(e["stim"] == -2))
            row["exported_unassigned"] = int(np.sum(e["stim"] == -1))
            row["exported_minus3"] = int(np.sum(e["stim"] == -3))
            labs = labs_e[:r["n_roi"]]
            same = [a.strip() == b.strip() for a, b in zip(labs, r["roi_label"])]
            row["labels_equal"] = float(np.mean(same))
            named = [(a.strip(), b.strip()) for a, b in zip(labs, r["roi_label"]) if a.strip() not in ("", "merge") or b.strip()]
            row["labels_equal_named"] = float(np.mean([a == b for a, b in named])) if named else np.nan
            row["export_merge"] = sum(a.strip() == "merge" for a in labs)
    rows.append(row)
ol = pd.DataFrame(rows)
mt = ol[ol.matched & ol.same_stims.fillna(False).astype(bool)]
print("Recordings with traces: %d; matched to an export by session key: %d; with the same number of stimulations: %d"
      % (len(ol), ol.matched.sum(), len(mt)))
print("Exported recordings with no recording here:", sorted(set(exp) - set(ol.key)))
print("Session ids of the subjects with an imaging file only:", {s_: sub_img[s_] for s_ in sorted(set(sub_img) - set(sub_seg))})
print("ROI count equal to the exported label count in %d of %d matched recordings; differing:"
      % ((mt.rois_here == mt.labels_exported).sum(), len(mt)), mt[mt.rois_here != mt.labels_exported][["subject", "rois_here", "labels_exported"]].values.tolist())
print("Stimulations assigned by both: %d; same ROI: %d (%.3f)" % (mt.both_assigned.sum(), mt.agree.sum(), mt.agree.sum() / mt.both_assigned.sum()))
print("  manual targets: %d of %d agree; nearest-centroid targets: %d of %d agree (%.3f)"
      % (mt.agree_manual.sum(), mt.n_manual.sum(), mt.agree_nearest.sum(), mt.n_nearest.sum(), mt.agree_nearest.sum() / mt.n_nearest.sum()))
print("Exported stimulations marked excluded (-2): %d, unassigned (-1): %d, labelled outside the reference volume (-3): %d"
      % (mt.exported_excluded.sum(), mt.exported_unassigned.sum(), mt.exported_minus3.sum()))
print("Share of ROI labels equal to the exported labels: median %.3f; recordings below 0.9: %d"
      % (mt.labels_equal.median(), (mt.labels_equal < 0.9).sum()))
print("Exported labels 'merge' (blank here): %d ROIs in %d recordings" % (mt.export_merge.sum(), (mt.export_merge > 0).sum()))
print("Over ROIs named on either side, leaving out 'merge': median agreement %.3f; recordings below 0.9 (subject, agreement):"
      % mt.labels_equal_named.median(), mt[mt.labels_equal_named < 0.9][["subject", "labels_equal_named"]].round(3).values.tolist())
per = (mt.agree / mt.both_assigned)
print("Agreement per recording: min %.2f, median %.2f; recordings below 0.6: %s"
      % (per.min(), per.median(), mt[per < 0.6].subject.tolist()))
low_recs = set(mt[per < 0.6].subject.tolist())
mt = mt.merge(tg[["subject", "z_slope", "z_offset"]], on="subject", how="left")
mt["pos_agree"] = mt.agree_nearest / mt.n_nearest
print("\nPosition-found agreement per recording against the fitted z slope:")
print("  z slope exactly 0: %d recordings, median agreement %.3f; nonzero: %d, median %.3f"
      % ((mt.z_slope == 0).sum(), mt.pos_agree[mt.z_slope == 0].median(), (mt.z_slope != 0).sum(), mt.pos_agree[mt.z_slope != 0].median()))
print("  Spearman of agreement with |slope|: %.3f" % stats.spearmanr(mt.pos_agree, mt.z_slope.abs(), nan_policy="omit").statistic)
print("  Recordings with overall agreement below 0.6:")
print(mt[mt.subject.isin(low_recs)][["subject", "both_assigned", "agree", "n_nearest", "agree_nearest", "pos_agree", "z_slope", "z_offset"]].round(3).to_string(index=False))
# Recordings whose labels disagree with the export: the first named ROIs side by side.
for s_ in mt[mt.labels_equal_named < 0.9].subject:
    r = next(r for r in res if int(r["subject_id"]) == s_)
    key = r["session_start"][:10].replace("-", "") + "-" + r["session_start"][11:19].replace(":", "")
    le = [x.strip() for x in exp[key]["labels"]]
    le = le[:-1] if le and le[-1] == "" else le   # trailing empty entry, as above
    lf = [x.strip() for x in r["roi_label"]]
    pairs_ = [(k, le[k] if k < len(le) else "-", lf[k]) for k in range(len(lf)) if (k < len(le) and le[k] not in ("", "merge")) or lf[k]][:8]
    se, sf = {x for x in le if x not in ("", "merge")}, {x for x in lf if x}
    print("\nsubject %d: ROIs %d here, %d exported; distinct names exported %d, here %d, shared %d" % (s_, len(lf), len(le), len(se), len(sf), len(se & sf)))
    print("  first named ROIs (index, exported, here):", pairs_)
# RID among the exported labels and stimulated ROIs.
rid_lab = rid_stim = 0
for r in res:
    key = r["session_start"][:10].replace("-", "") + "-" + r["session_start"][11:19].replace(":", "")
    e = exp.get(key)
    if e is None:
        continue
    le = [x.strip() for x in e["labels"]]
    rid_lab += sum(x == "RID" for x in le)
    rid_stim += sum(1 for s in e["stim"] if 0 <= s < len(le) and le[s] == "RID")
print("\nExported ROI labels RID: %d; exported stimulations whose ROI is labelled RID: %d" % (rid_lab, rid_stim))

Recordings with traces: 110; matched to an export by session key: 110; with the same number of stimulations: 110
Exported recordings with no recording here: ['20210917-132736', '20211008-115929', '20220330-110253']
Session ids of the subjects with an imaging file only: {20: '20210917-132736', 23: '20211008-115929', 33: '20220330-110253'}
ROI count equal to the exported label count in 109 of 110 matched recordings; differing: [[11, 161, 200]]
Stimulations assigned by both: 4479; same ROI: 3913 (0.874)
  manual targets: 2812 of 2812 agree; nearest-centroid targets: 1101 of 1667 agree (0.660)
Exported stimulations marked excluded (-2): 950, unassigned (-1): 156, labelled outside the reference volume (-3): 28
Share of ROI labels equal to the exported labels: median 0.975; recordings below 0.9: 23
Exported labels 'merge' (blank here): 732 ROIs in 63 recordings
Over ROIs named on either side, leaving out 'merge': median agreement 1.000; recordings below 0.9 (subject, agreement): [[11.0, 0.0], [64.0, 0.0], [76.0, 0.0], [77.0, 0.0], [86.0, 0.0], [91.0, 0.0]]
Agreement per recording: min 0.24, median 0.97; recordings below 0.6: [27, 28, 30, 35, 42, 50, 51, 71, 86]

Position-found agreement per recording against the fitted z slope:
  z slope exactly 0: 3 recordings, median agreement 0.579; nonzero: 107, median 0.739
  Spearman of agreement with |slope|: 0.478
  Recordings with overall agreement below 0.6:
 subject  both_assigned  agree  n_nearest  agree_nearest  pos_agree  z_slope  z_offset
      27             27     13         17              3      0.176     -2.0      13.5
      28             33      8         31              6      0.194      2.5      13.0
      30             47     20         35              8      0.229      0.5      10.5
      35             51     26         39             14      0.359     -8.5      13.0
      42             37     19         20              2      0.100      5.0      14.0
      50             34     18         22              6      0.273     -5.5      14.0
      51             37     22         23              8      0.348     -7.0      11.5
      71             38     18         23              3      0.130     -8.0      11.0
      86             43     24         27              8      0.296     15.5      13.0

subject 11: ROIs 161 here, 200 exported; distinct names exported 0, here 41, shared 0
  first named ROIs (index, exported, here): [(1, '', 'I2L'), (3, '', 'MCL'), (4, '', 'RMEL'), (6, '', 'I1L'), (7, '', 'AMso'), (8, '', 'M3L'), (9, '', 'OLLL'), (11, '', 'AVEL')]

subject 64: ROIs 110 here, 110 exported; distinct names exported 70, here 4, shared 0
  first named ROIs (index, exported, here): [(2, 'DB1', ''), (3, 'SABD', 'AIY3'), (4, 'RIAR', ''), (5, 'RMDVR', 'AFDR'), (6, 'BAGR', ''), (8, 'DA1', ''), (14, 'DB2', ''), (16, 'SABVL', '')]

subject 76: ROIs 162 here, 162 exported; distinct names exported 102, here 0, shared 0
  first named ROIs (index, exported, here): [(2, 'RMDVL', ''), (3, 'AWCON', ''), (5, 'I2L', ''), (7, 'RMEL', ''), (8, 'OLLL', ''), (9, 'SAAVL', ''), (11, 'AVDL', ''), (12, 'AWAL', '')]

subject 77: ROIs 140 here, 140 exported; distinct names exported 87, here 15, shared 2
  first named ROIs (index, exported, here): [(0, 'RMDVR', ''), (1, 'OLQDR', ''), (2, 'FLPR', ''), (4, 'I1R', ''), (6, '', 'target'), (7, 'M3R', ''), (8, 'OLQVR', 'target'), (9, 'URYVR', '')]

subject 86: ROIs 119 here, 119 exported; distinct names exported 76, here 69, shared 35
  first named ROIs (index, exported, here): [(8, 'ADAL', 'AVHL'), (9, 'RMDVL', 'ADAL'), (10, 'AINL', 'AWCON'), (12, 'RIAL', ''), (13, 'AWBL', 'FLPL'), (14, 'M5', 'URBR'), (15, 'FLPL', ''), (16, '', 'ASJR')]

subject 91: ROIs 147 here, 147 exported; distinct names exported 97, here 81, shared 59
  first named ROIs (index, exported, here): [(0, 'RMDVL', 'AVEL'), (1, 'URBL', 'RIFL'), (2, 'I2L', 'DB2'), (4, 'M5', 'AIYR'), (5, 'ASHL', 'target'), (7, 'AVEL', 'AVDL'), (8, 'CEPVL', ''), (9, 'URYDL', 'M4')]

Exported ROI labels RID: 12; exported stimulations whose ROI is labelled RID: 12

# Variant: pairs-measured with the exported assignment of the stimulated ROI (and its exclusions).
# The autoresponse is re-evaluated for the exported ROI.
res_osf = []
for r, roi in zip(res, osf_roi):
    res_osf.append(roi)
obs_osf = sp.observations(res, head, stim_roi=res_osf)
for name, dg in [("(a) with i = j", True), ("(b) i != j", False)]:
    n, hh = pair_counts(obs_osf, True, dg)
    print("Exported assignment, reading %s: %d pairs, %d neurons recorded, max %d observations"
          % (name, len(n), hh[hh.i != hh.j].i.nunique(), n.max()))
h_osf = obs_osf[obs_osf.i_head & obs_osf.j_head & obs_osf.auto]
main_osf = h_osf.groupby(["i", "j"]).dff.agg(["size", "mean"]).rename(columns={"size": "n"})
# Variant: exported assignment and exported ROI labels.
osf_lab = []
for r in res:
    key = r["session_start"][:10].replace("-", "") + "-" + r["session_start"][11:19].replace(":", "")
    e = exp.get(key)
    le = None
    if e is not None:
        le = e["labels"][:-1] if e["labels"] and e["labels"][-1] == "" else e["labels"]
        le = np.array([x.strip() for x in le]) if len(le) == r["n_roi"] else None
    osf_lab.append(le)
obs_osf2 = sp.observations(res, head, stim_roi=res_osf, labels=osf_lab)
n2, hh2 = pair_counts(obs_osf2, True, True)
print("Exported assignment and labels, reading (a): %d pairs, %d neurons recorded, max %d observations (recordings used: %d)"
      % (len(n2), hh2[hh2.i != hh2.j].i.nunique(), n2.max(), obs_osf2.rec.nunique()))
h_osf2 = obs_osf2[obs_osf2.i_head & obs_osf2.j_head & obs_osf2.auto]
main_osf2 = h_osf2.groupby(["i", "j"]).dff.agg(["size", "mean"]).rename(columns={"size": "n"})
print(pd.DataFrame([link_stats(main_osf, "exported assignment"),
                    link_stats(main_osf2, "exported assignment and labels")]).round(3).to_string(index=False))

Exported assignment, reading (a) with i = j: 20351 pairs, 183 neurons recorded, max 55 observations
Exported assignment, reading (b) i != j: 20189 pairs, 183 neurons recorded, max 55 observations

Exported assignment and labels, reading (a): 20564 pairs, 184 neurons recorded, max 56 observations (recordings used: 107)
                       reading  paper_pairs_found  pairs_here_not_in_paper  spearman_counts  pearson_dff  spearman_dff  median_count_ratio  pairs_compared
           exported assignment              0.833                      824            0.820        0.675         0.573                 1.0           19527
exported assignment and labels              0.845                      757            0.834        0.684         0.586                 1.0           19807

h = obs[obs.i_head & obs.j_head & obs.auto]
main = h.groupby(["i", "j"]).dff.agg(["size", "mean"]).rename(columns={"size": "n"})
dec = h.groupby(["i", "j"]).dff_decay.agg(["size", "mean"]).rename(columns={"size": "n"})
rsp = h[h.responded | (h.i == h.j)].groupby(["i", "j"]).dff.agg(["size", "mean"]).rename(columns={"size": "n"})
link = pd.DataFrame([link_stats(main, "main: all observations, Methods dF/F0"),
                     link_stats(dec, "variant: code's decaying correction"),
                     link_stats(rsp, "variant: detected responses only")])
print(link.round(3).to_string(index=False))
j = sd.join(main, how="inner")
print("\nCount differences (here minus paper) over pairs found: median %.0f, 5th %.0f, 95th %.0f percentile; equal in %.3f"
      % (np.median(j.n - j.n_paper), np.percentile(j.n - j.n_paper, 5), np.percentile(j.n - j.n_paper, 95), np.mean(j.n == j.n_paper)))
print("Sign of the mean response agrees in %.3f of pairs compared" % np.mean(np.sign(j["mean"]) == np.sign(j.dff_paper)))
# Variants that remove rows a non-positive F0, a failed photobleaching fit, the export's exclusions or
# a poorly linked recording could distort.
low_rec_idx = {k for k, r in enumerate(res) if int(r["subject_id"]) in low_recs}
def filtered(o):
    h = o[o.i_head & o.j_head & o.auto]
    ex = pd.Series([(a, b) in excluded for a, b in zip(h.rec, h.stim)], index=h.index)
    f = {"main": pd.Series(True, index=h.index),
         "F0 > 0 only": h.f0 > 0,
         "without the 9 recordings with negative traces": ~h.rec.isin(neg_recs),
         "photobleaching-corrected traces only": h.pb_ok,
         "without the export's excluded stimulations": ~ex,
         "without the 9 recordings with agreement < 0.6": ~h.rec.isin(low_rec_idx)}
    f["all of these filters"] = np.logical_and.reduce([v for k, v in f.items() if k != "main"])
    return {k: h[v].groupby(["i", "j"]).dff.agg(["size", "mean"]).rename(columns={"size": "n"}) for k, v in f.items()}
var_tabs = filtered(obs)
print("\nVariants (main reading otherwise):")
print(pd.DataFrame([link_stats(t, k) for k, t in var_tabs.items()]).round(3).to_string(index=False))

                              reading  paper_pairs_found  pairs_here_not_in_paper  spearman_counts  pearson_dff  spearman_dff  median_count_ratio  pairs_compared
main: all observations, Methods dF/F0              0.817                     1440            0.768        0.551         0.493               1.000           19147
  variant: code's decaying correction              0.817                     1440            0.768        0.616         0.620               1.000           19147
     variant: detected responses only              0.360                      305            0.471        0.458         0.517               0.286            8443

Count differences (here minus paper) over pairs found: median 0, 5th -5, 95th 3 percentile; equal in 0.274
Sign of the mean response agrees in 0.652 of pairs compared

Variants (main reading otherwise):

                                      reading  paper_pairs_found  pairs_here_not_in_paper  spearman_counts  pearson_dff  spearman_dff  median_count_ratio  pairs_compared
                                         main              0.817                     1440            0.768        0.551         0.493               1.000           19147
                                  F0 > 0 only              0.814                     1425            0.766        0.552         0.493               1.000           19069
without the 9 recordings with negative traces              0.799                     1369            0.752        0.534         0.477               0.846           18724
         photobleaching-corrected traces only              0.816                     1438            0.767        0.551         0.493               1.000           19126
   without the export's excluded stimulations              0.808                     1197            0.780        0.601         0.511               0.900           18936
without the 9 recordings with agreement < 0.6              0.799                     1246            0.770        0.543         0.494               0.900           18727
                         all of these filters              0.768                      994            0.765        0.569         0.489               0.731           17985

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].scatter(j.n_paper, j.n, s=3, alpha=0.3)
ax[0].plot([0, 60], [0, 60], "k--", lw=0.8)
ax[0].set_xlabel("observations, Source Data"); ax[0].set_ylabel("observations, this analysis")
lim = 1.0
inside = (j.dff_paper.abs() <= lim) & (j["mean"].abs() <= lim)
ax[1].scatter(j.dff_paper[inside], j["mean"][inside], s=3, alpha=0.3)
ax[1].plot([-lim, lim], [-lim, lim], "k--", lw=0.8)
ax[1].set_xlim(-lim, lim); ax[1].set_ylim(-lim, lim)
ax[1].set_xlabel("mean dF/F0, Source Data"); ax[1].set_ylabel("mean dF/F0, this analysis")
plt.tight_layout(); plt.show()
print("Counts: paper 1 to %d, here %d to %d. dF/F0 panel: %d of %d pairs outside +-%.0f not shown"
      % (j.n_paper.max(), j.n.min(), j.n.max(), (~inside).sum(), len(j), lim))

Counts: paper 1 to 59, here 1 to 50. dF/F0 panel: 22 of 19147 pairs outside +-1 not shown

x = obs[(obs.j == "AVJR") & (obs.i == "AVDR") & obs.auto]
x_all = obs[(obs.j == "AVJR") & (obs.i == "AVDR")]
print("AVJR stimulations with an autoresponse and AVDR observed: %d trials in %d animals (%d without the autoresponse requirement)"
      % (len(x), x.rec.nunique(), len(x_all)))
print("Mean <dF/F0>_t of AVDR: %.4f (median %.4f, range %.3f to %.3f); trials with a detected AVDR response: %d"
      % (x.dff.mean(), x.dff.median(), x.dff.min(), x.dff.max(), x.responded.sum()))
print("Source Data value for this pair: mean %.4f over %d observations, q = %.3g"
      % (sd.dff_paper.get(("AVDR", "AVJR"), np.nan), sd.n_paper.get(("AVDR", "AVJR"), 0), sd.q.get(("AVDR", "AVJR"), np.nan)))
print("Variants: code's decaying correction, mean %.4f; detected responses only, mean %.4f (n = %d); without the autoresponse requirement, mean %.4f"
      % (x.dff_decay.mean(), x[x.responded].dff.mean(), x.responded.sum(), x_all.dff.mean()))
if len(x) > 1:
    w = stats.wilcoxon(x.dff)
    print("Wilcoxon signed-rank test of the trials against 0: p = %.3g" % w.pvalue)
print("Per animal:")
per_an = x.groupby("subject").agg(trials=("dff", "size"), mean_dff=("dff", "mean"), responses=("responded", "sum"))
print("Animals whose mean has the expected sign (positive): %d of %d" % ((np.sign(per_an.mean_dff) == 1).sum(), len(per_an)))
print(per_an.round(4).to_string())

AVJR stimulations with an autoresponse and AVDR observed: 27 trials in 17 animals (37 without the autoresponse requirement)
Mean <dF/F0>_t of AVDR: 0.2364 (median 0.1689, range -0.136 to 1.156); trials with a detected AVDR response: 12
Source Data value for this pair: mean 0.2461 over 24 observations, q = 6.07e-07
Variants: code's decaying correction, mean 0.1971; detected responses only, mean 0.3476 (n = 12); without the autoresponse requirement, mean 0.1787
Wilcoxon signed-rank test of the trials against 0: p = 5.49e-05
Per animal:
Animals whose mean has the expected sign (positive): 16 of 17
         trials  mean_dff  responses
subject                             
0             3    0.1071          0
107           1    0.4532          1
108           2    0.1089          0
110           1    0.4715          1
111           1    0.0482          0
14            1    0.4495          1
30            2    0.1368          1
34            4    0.3289          4
36            3    0.1130          2
4             1    0.1018          0
40            1    0.1321          1
45            1    1.1556          0
63            1    0.2825          0
67            2   -0.0907          0
69            1    0.1774          0
75            1    0.1029          0
98            1    0.7217          1

x = obs[(obs.j == "AVER") & (obs.i == "AVAR") & obs.auto]
x_all = obs[(obs.j == "AVER") & (obs.i == "AVAR")]
print("AVER stimulations with an autoresponse and AVAR observed: %d trials in %d animals (%d without the autoresponse requirement)"
      % (len(x), x.rec.nunique(), len(x_all)))
print("Mean <dF/F0>_t of AVAR: %.4f (median %.4f, range %.3f to %.3f); trials with a detected AVAR response: %d"
      % (x.dff.mean(), x.dff.median(), x.dff.min(), x.dff.max(), x.responded.sum()))
print("Source Data value for this pair: mean %.4f over %d observations, q = %.3g"
      % (sd.dff_paper.get(("AVAR", "AVER"), np.nan), sd.n_paper.get(("AVAR", "AVER"), 0), sd.q.get(("AVAR", "AVER"), np.nan)))
print("Variants: code's decaying correction, mean %.4f; detected responses only, mean %.4f (n = %d); without the autoresponse requirement, mean %.4f"
      % (x.dff_decay.mean(), x[x.responded].dff.mean(), x.responded.sum(), x_all.dff.mean()))
if len(x) > 1:
    w = stats.wilcoxon(x.dff)
    print("Wilcoxon signed-rank test of the trials against 0: p = %.3g" % w.pvalue)
print("Per animal:")
per_an = x.groupby("subject").agg(trials=("dff", "size"), mean_dff=("dff", "mean"), responses=("responded", "sum"))
print("Animals whose mean has the expected sign (positive): %d of %d" % ((np.sign(per_an.mean_dff) == 1).sum(), len(per_an)))
print(per_an.round(4).to_string())

AVER stimulations with an autoresponse and AVAR observed: 7 trials in 5 animals (13 without the autoresponse requirement)
Mean <dF/F0>_t of AVAR: 0.1589 (median 0.1204, range -0.065 to 0.495); trials with a detected AVAR response: 4
Source Data value for this pair: mean 0.2044 over 11 observations, q = 0.00104
Variants: code's decaying correction, mean 0.1632; detected responses only, mean 0.3015 (n = 4); without the autoresponse requirement, mean 0.0921
Wilcoxon signed-rank test of the trials against 0: p = 0.219
Per animal:
Animals whose mean has the expected sign (positive): 3 of 5
         trials  mean_dff  responses
subject                             
26            1   -0.0179          0
4             1    0.1204          1
53            2    0.3153          2
70            2   -0.0376          0
73            1    0.4548          1

x = obs[(obs.j == "SAADL") & (obs.i == "OLLR") & obs.auto]
x_all = obs[(obs.j == "SAADL") & (obs.i == "OLLR")]
print("SAADL stimulations with an autoresponse and OLLR observed: %d trials in %d animals (%d without the autoresponse requirement)"
      % (len(x), x.rec.nunique(), len(x_all)))
print("Mean <dF/F0>_t of OLLR: %.4f (median %.4f, range %.3f to %.3f); trials with a detected OLLR response: %d"
      % (x.dff.mean(), x.dff.median(), x.dff.min(), x.dff.max(), x.responded.sum()))
print("Source Data value for this pair: mean %.4f over %d observations, q = %.3g"
      % (sd.dff_paper.get(("OLLR", "SAADL"), np.nan), sd.n_paper.get(("OLLR", "SAADL"), 0), sd.q.get(("OLLR", "SAADL"), np.nan)))
print("Variants: code's decaying correction, mean %.4f; detected responses only, mean %.4f (n = %d); without the autoresponse requirement, mean %.4f"
      % (x.dff_decay.mean(), x[x.responded].dff.mean(), x.responded.sum(), x_all.dff.mean()))
if len(x) > 1:
    w = stats.wilcoxon(x.dff)
    print("Wilcoxon signed-rank test of the trials against 0: p = %.3g" % w.pvalue)
print("Per animal:")
per_an = x.groupby("subject").agg(trials=("dff", "size"), mean_dff=("dff", "mean"), responses=("responded", "sum"))
print("Animals whose mean has the expected sign (negative): %d of %d" % ((np.sign(per_an.mean_dff) == -1).sum(), len(per_an)))
print(per_an.round(4).to_string())

SAADL stimulations with an autoresponse and OLLR observed: 5 trials in 4 animals (9 without the autoresponse requirement)
Mean <dF/F0>_t of OLLR: -0.0591 (median -0.0940, range -0.233 to 0.202); trials with a detected OLLR response: 2
Source Data value for this pair: mean -0.0404 over 5 observations, q = 0.00326
Variants: code's decaying correction, mean 0.0136; detected responses only, mean -0.0154 (n = 2); without the autoresponse requirement, mean -0.0249
Wilcoxon signed-rank test of the trials against 0: p = 0.438
Per animal:
Animals whose mean has the expected sign (negative): 3 of 4
         trials  mean_dff  responses
subject                             
1             2   -0.1170          0
101           1   -0.2330          1
50            1    0.2023          1
94            1   -0.0307          0

step = dict(zip(dep.subject, dep.plane_step_um))
rows = []
for r in res:
    s = step[int(r["subject_id"])]
    for e in np.where(r["auto"])[0]:
        jj = r["stim_roi"][e]
        dist = sp.roi_distance_um(r, np.arange(r["n_roi"]), jj, s)
        v = r["dff"][e]
        ok = np.isfinite(v)
        rows.append(pd.DataFrame({"rec": int(r["subject_id"]), "dist": dist[ok], "dff": v[ok], "f0": r["f0"][e][ok]}))
sdist = pd.concat(rows, ignore_index=True)
bands = {"0 um (target)": sdist.dist == 0, "4-10 um": (sdist.dist >= 4) & (sdist.dist < 10),
         "10-16 um": (sdist.dist >= 10) & (sdist.dist < 16), "4-16 um": (sdist.dist >= 4) & (sdist.dist < 16)}
tab = pd.DataFrame([{"band": b, "n": int(m.sum()), "median": sdist.dff[m].median(),
                     "p75": sdist.dff[m].quantile(0.75), "frac_below_0.1": (sdist.dff[m] < 0.1).mean()}
                    for b, m in bands.items()])
print(tab.round(3).to_string(index=False))
ratio = tab.loc[0, "median"] / tab.loc[1, "median"]
print("Target median over 4-10 um median: %.1f (criterion at least 3)" % ratio)
print("Expected fraction below 0.1: 0.75 (tolerance 0.075); the Source Data give 0.616 (4-10 um), 0.823 (10-16 um), 0.749 pooled")
for thr in (0.05, 0.1, 0.2):
    print("Variant threshold %.2f: fraction of 4-16 um neighbours below it %.3f" % (thr, (sdist.dff[bands['4-16 um']] < thr).mean()))
per = sdist[bands["4-16 um"]].groupby("rec").dff.apply(lambda v: (v < 0.1).mean())
print("Per animal, fraction of 4-16 um neighbours below 0.1: min %.2f, median %.2f, max %.2f (%d animals)"
      % (per.min(), per.median(), per.max(), len(per)))
pos = sdist.f0 > 0
for b in ["4-10 um", "10-16 um", "4-16 um"]:
    print("Variant F0 > 0 only, %s: fraction below 0.1 %.3f (n = %d)" % (b, (sdist.dff[bands[b] & pos] < 0.1).mean(), (bands[b] & pos).sum()))

         band     n  median   p75  frac_below_0.1
0 um (target)  2448   0.376 0.601           0.092
      4-10 um 20358   0.061 0.183           0.604
     10-16 um 36157   0.006 0.071           0.817
      4-16 um 56515   0.022 0.105           0.741
Target median over 4-10 um median: 6.2 (criterion at least 3)
Expected fraction below 0.1: 0.75 (tolerance 0.075); the Source Data give 0.616 (4-10 um), 0.823 (10-16 um), 0.749 pooled
Variant threshold 0.05: fraction of 4-16 um neighbours below it 0.608
Variant threshold 0.10: fraction of 4-16 um neighbours below it 0.741
Variant threshold 0.20: fraction of 4-16 um neighbours below it 0.869
Per animal, fraction of 4-16 um neighbours below 0.1: min 0.50, median 0.77, max 0.97 (110 animals)
Variant F0 > 0 only, 4-10 um: fraction below 0.1 0.605 (n = 20180)
Variant F0 > 0 only, 10-16 um: fraction below 0.1 0.817 (n = 35913)
Variant F0 > 0 only, 4-16 um: fraction below 0.1 0.741 (n = 56093)

fig, ax = plt.subplots(figsize=(6, 4))
for b in ["0 um (target)", "4-10 um", "10-16 um"]:
    v = np.sort(sdist.dff[bands[b]].values)
    ax.plot(v, np.arange(1, len(v) + 1) / len(v), label=b)
ax.axvline(0.1, color="k", lw=0.8, ls="--"); ax.axhline(0.75, color="k", lw=0.8, ls=":")
ax.set_xlim(-0.5, 1.5); ax.set_xlabel("<dF/F0>_t"); ax.set_ylabel("cumulative fraction"); ax.legend()
plt.show()
print("Values outside the x axis (-0.5 to 1.5): %d of %d" % (((sdist.dff < -0.5) | (sdist.dff > 1.5)).sum(), len(sdist)))

Values outside the x axis (-0.5 to 1.5): 508 of 244219

sig_pairs = sd[sd.q < 0.05]
jq = sig_pairs.join(main, how="left")
for name, sel in [("(a) including i = j", np.ones(len(jq), bool)),
                  ("(b) excluding i = j", jq.index.get_level_values(0) != jq.index.get_level_values(1))]:
    q = jq[sel]
    found = q["mean"].notna()
    print("%s: q < 0.05 pairs %d, observed here %d; negative mean here %d (%.4f); negative in the Source Data %d (%.4f)"
          % (name, len(q), found.sum(), (q["mean"][found] < 0).sum(), (q["mean"][found] < 0).mean(),
             (q.dff_paper < 0).sum(), (q.dff_paper < 0).mean()))
print("Expected 0.11 (tolerance 0.011)")
qd = jq[jq["mean"].notna()]
print("Sign agreement with the Source Data over these pairs: %.3f" % np.mean(np.sign(qd["mean"]) == np.sign(qd.dff_paper)))
dq = sig_pairs.join(dec, how="left")
print("Variant, code's decaying correction, reading (a): negative %.4f" % (dq["mean"].dropna() < 0).mean())
for thr in (0.01, 0.02, 0.1):
    qq = sd[sd.q < thr].join(main, how="left")["mean"].dropna()
    print("Variant q < %.2f: %d pairs, negative %.4f" % (thr, len(qq), (qq < 0).mean()))
def neg_frac(t):
    qq = sig_pairs.join(t, how="left")["mean"]
    off = qq.index.get_level_values(0) != qq.index.get_level_values(1)
    return {"observed (a)": int(qq.notna().sum()), "negative (a)": (qq.dropna() < 0).mean(),
            "negative (b)": (qq[off].dropna() < 0).mean()}
print("\nVariants, fraction of q < 0.05 pairs with a negative mean (expected 0.11, tolerance 0.011):")
vt = {**var_tabs, "exported assignment (its exclusions applied)": main_osf, "exported assignment and labels": main_osf2,
      "exported assignment, F0 > 0 and corrected traces only": h_osf[(h_osf.f0 > 0) & h_osf.pb_ok].groupby(["i", "j"]).dff.agg(["size", "mean"]).rename(columns={"size": "n"})}
print(pd.DataFrame([{"variant": k, **neg_frac(t)} for k, t in vt.items()]).round(4).to_string(index=False))

(a) including i = j: q < 0.05 pairs 1310, observed here 1161; negative mean here 175 (0.1507); negative in the Source Data 151 (0.1153)
(b) excluding i = j: q < 0.05 pairs 1151, observed here 1005; negative mean here 175 (0.1741); negative in the Source Data 151 (0.1312)
Expected 0.11 (tolerance 0.011)
Sign agreement with the Source Data over these pairs: 0.891
Variant, code's decaying correction, reading (a): negative 0.0956
Variant q < 0.01: 689 pairs, negative 0.0726
Variant q < 0.02: 837 pairs, negative 0.0920
Variant q < 0.10: 1681 pairs, negative 0.2124

Variants, fraction of q < 0.05 pairs with a negative mean (expected 0.11, tolerance 0.011):
                                              variant  observed (a)  negative (a)  negative (b)
                                                 main          1161        0.1507        0.1741
                                          F0 > 0 only          1160        0.1509        0.1743
        without the 9 recordings with negative traces          1137        0.1530        0.1772
                 photobleaching-corrected traces only          1161        0.1507        0.1741
           without the export's excluded stimulations          1157        0.1521        0.1758
        without the 9 recordings with agreement < 0.6          1141        0.1560        0.1803
                                 all of these filters          1108        0.1579        0.1832
         exported assignment (its exclusions applied)          1179        0.1399        0.1611
                       exported assignment and labels          1190        0.1336        0.1535
exported assignment, F0 > 0 and corrected traces only          1179        0.1416        0.1631

FIG4 = {"URX": "URXL", "ADL": "ADLR", "AWB": "AWBR"}
def rid_rows(o, label):
    x = o[o.j == "RID"]
    n_stim = x.groupby(["rec", "stim"]).ngroups
    x = x[x.dff.abs() <= 2]
    rows = []
    for cls, shown in FIG4.items():
        for reading, names in [("(a) %s, L and R" % cls, [cls + "L", cls + "R"]), ("(b) %s" % shown, [shown])]:
            y = x[x.i.isin(names)]
            rows.append({"source": label, "reading": reading, "neuron": cls, "rid_stims": n_stim, "trials": len(y),
                         "animals": y.rec.nunique(), "mean": y.dff.mean() if len(y) else np.nan,
                         "median": y.dff.median() if len(y) else np.nan, "detected": int(y.responded.sum()),
                         "positive_trials": int((y.dff > 0).sum()),
                         "pctile_among_head_pairs": stats.percentileofscore(main["mean"].dropna(), y.dff.mean()) if len(y) else np.nan})
    return rows

rid_stims = obs[obs.j == "RID"].groupby(["rec", "stim"]).agg(subject=("subject", "first"), how=("how", "first"), auto=("auto", "first"))
print("RID stimulations here: %d in %d animals (manual targets %d, nearest centroid %d); with a RID autoresponse: %d"
      % (len(rid_stims), rid_stims.subject.nunique(), (rid_stims.how == "manual").sum(), (rid_stims.how == "nearest").sum(), rid_stims.auto.sum()))
xr = obs[(obs.j == "RID")]
print("Observations of URX, ADL or AWB after a RID stimulation: %d; beyond +-2 and excluded: %d"
      % (xr.i.str[:3].isin(list(FIG4)).sum(), (xr.i.str[:3].isin(list(FIG4)) & (xr.dff.abs() > 2)).sum()))
rid_wt = pd.DataFrame(rid_rows(obs, "dandiset, assignment here")
                      + rid_rows(obs[obs.auto], "dandiset, autoresponse required")
                      + rid_rows(obs_osf, "dandiset traces, exported assignment")
                      + rid_rows(obs_osf2, "dandiset traces, exported assignment and labels"))
print(rid_wt.round(3).to_string(index=False))
print("\nPer trial, main source (subject, stimulation, observed neuron, <dF/F0>_t, detected):")
xx = xr[xr.i.str[:3].isin(list(FIG4)) & (xr.dff.abs() <= 2)]
print(xx[["subject", "stim", "how", "i", "dff", "responded"]].round(3).to_string(index=False))

RID stimulations here: 8 in 5 animals (manual targets 5, nearest centroid 3); with a RID autoresponse: 5
Observations of URX, ADL or AWB after a RID stimulation: 14; beyond +-2 and excluded: 0

                                         source          reading neuron  rid_stims  trials  animals  mean  median  detected  positive_trials  pctile_among_head_pairs
                      dandiset, assignment here (a) URX, L and R    URX          8       1        1 0.032   0.032         0                1                   67.309
                      dandiset, assignment here         (b) URXL    URX          8       0        0   NaN     NaN         0                0                      NaN
                      dandiset, assignment here (a) ADL, L and R    ADL          8       0        0   NaN     NaN         0                0                      NaN
                      dandiset, assignment here         (b) ADLR    ADL          8       0        0   NaN     NaN         0                0                      NaN
                      dandiset, assignment here (a) AWB, L and R    AWB          8      12        5 0.076   0.052         3                7                   81.692
                      dandiset, assignment here         (b) AWBR    AWB          8       7        4 0.094   0.053         2                4                   85.233
                dandiset, autoresponse required (a) URX, L and R    URX          5       1        1 0.032   0.032         0                1                   67.309
                dandiset, autoresponse required         (b) URXL    URX          5       0        0   NaN     NaN         0                0                      NaN
                dandiset, autoresponse required (a) ADL, L and R    ADL          5       0        0   NaN     NaN         0                0                      NaN
                dandiset, autoresponse required         (b) ADLR    ADL          5       0        0   NaN     NaN         0                0                      NaN
                dandiset, autoresponse required (a) AWB, L and R    AWB          5       8        4 0.056   0.019         2                4                   76.922
                dandiset, autoresponse required         (b) AWBR    AWB          5       4        3 0.071   0.020         1                2                   80.585
           dandiset traces, exported assignment (a) URX, L and R    URX         11       1        1 0.032   0.032         0                1                   67.309
           dandiset traces, exported assignment         (b) URXL    URX         11       0        0   NaN     NaN         0                0                      NaN
           dandiset traces, exported assignment (a) ADL, L and R    ADL         11       0        0   NaN     NaN         0                0                      NaN
           dandiset traces, exported assignment         (b) ADLR    ADL         11       0        0   NaN     NaN         0                0                      NaN
           dandiset traces, exported assignment (a) AWB, L and R    AWB         11      14        5 0.100   0.071         5                9                   86.394
           dandiset traces, exported assignment         (b) AWBR    AWB         11       9        4 0.127   0.116         4                6                   89.853
dandiset traces, exported assignment and labels (a) URX, L and R    URX         12       1        1 0.032   0.032         0                1                   67.309
dandiset traces, exported assignment and labels         (b) URXL    URX         12       0        0   NaN     NaN         0                0                      NaN
dandiset traces, exported assignment and labels (a) ADL, L and R    ADL         12       0        0   NaN     NaN         0                0                      NaN
dandiset traces, exported assignment and labels         (b) ADLR    ADL         12       0        0   NaN     NaN         0                0                      NaN
dandiset traces, exported assignment and labels (a) AWB, L and R    AWB         12      16        6 0.087   0.052         5                9                   83.995
dandiset traces, exported assignment and labels         (b) AWBR    AWB         12      10        5 0.113   0.085         4                6                   88.255

Per trial, main source (subject, stimulation, observed neuron, <dF/F0>_t, detected):
subject  stim     how    i    dff  responded
     25     0  manual AWBL -0.074      False
     25     0  manual AWBR -0.054      False
     25     1  manual AWBL  0.090      False
     25     1  manual AWBR -0.008      False
     25     2  manual AWBL  0.051      False
     25     2  manual AWBR -0.013      False
     26    49 nearest AWBL -0.103      False
     26    49 nearest  ADL  0.011       True
     41    50  manual AWBR  0.298       True
     41    50  manual URXR  0.032      False
     41    50  manual AWBL  0.294       True
     42     0 nearest AWBR  0.264       True
     42     2  manual AWBR  0.053      False
     84    19 nearest AWBR  0.116      False

# The paper's unc-31 export on OSF (doi 10.17605/OSF.IO/E2SYT), file exported_data_unc31.tar.gz
# (guid ty6sp, version 1, md5 e5c83242baa131bca95ee58bc7b58bec), streamed once.
U_URL, U_MD5 = "https://osf.io/download/ty6sp/", "e5c83242baa131bca95ee58bc7b58bec"
utx = {}
with urllib.request.urlopen(U_URL, timeout=600) as r:
    hr = _Hashed(r)
    with tarfile.open(fileobj=hr, mode="r|gz") as tf:
        for m in tf:
            mm = re.search(r"/(\d+)_(ds_name|stim_neurons|labels|gcamp|t|stim_volume_i)\.txt$", m.name)
            if mm:
                utx[(int(mm.group(1)), mm.group(2))] = tf.extractfile(m).read().decode()
    while hr.read(1 << 20):
        pass
print("exported_data_unc31.tar.gz: %d bytes, md5 %s, matches the OSF record: %s" % (hr.n, hr.md5.hexdigest(), hr.md5.hexdigest() == U_MD5))

def export_recording(get, k, name):
    g = np.loadtxt(io.StringIO(get(k, "gcamp")), ndmin=2)
    t = np.loadtxt(io.StringIO(get(k, "t")))
    sv = np.loadtxt(io.StringIO(get(k, "stim_volume_i")), ndmin=1).astype(int)
    stim = np.array([int(x) for x in get(k, "stim_neurons").split()])
    lab = [x.strip() for x in get(k, "labels").split("\n")]
    lab = (lab + [""] * g.shape[1])[:g.shape[1]]
    same = len(sv) == len(stim)
    return sp.analyse_export(g, t, sv, lab, name), (stim if same else None), same, g.shape

u_keys = sorted({k for k, _ in utx})
u_res, u_stim, u_info = [], [], []
for k in u_keys:
    an, stim, same, shape = export_recording(lambda kk, w: utx[(kk, w)], k, "unc31-%d" % k)
    le = list(an["roi_label"])
    n_rid = sum(1 for s in (stim if stim is not None else []) if 0 <= s < len(le) and le[s] == "RID")
    u_info.append({"k": k, "folder": utx[(k, "ds_name")].strip().split("/")[-2], "volumes": shape[0], "rois": shape[1],
                   "stims": len(stim) if stim is not None else np.nan, "stim_counts_match": same, "rid_stims": n_rid,
                   "pb_not_ok": int(np.sum(an["photobleach_status"] != "ok"))})
    u_res.append(an); u_stim.append(stim)
print(pd.DataFrame(u_info).to_string(index=False))
obs_u = sp.observations(u_res, head, stim_roi=u_stim)

wt_res, wt_stim = [], []
for k in rid_k:
    an, stim, same, shape = export_recording(lambda kk, w: zlib.decompress(big_txt[(kk, w)]).decode() if w in ("gcamp", "t", "stim_volume_i") else texts[(kk, w)], k, "wt-export-%d" % k)
    wt_res.append(an); wt_stim.append(stim)
obs_wte = sp.observations(wt_res, head, stim_roi=wt_stim)
print("WT export recordings processed for the variant:", len(wt_res))

exported_data_unc31.tar.gz: 84924311 bytes, md5 e5c83242baa131bca95ee58bc7b58bec, matches the OSF record: True

 k                    folder  volumes  rois  stims  stim_counts_match  rid_stims  pb_not_ok
 0 pumpprobe_20211216_164305     4080    91     60               True          0          0
 1 pumpprobe_20220113_101730     5035   126     79               True          0          4
 2 pumpprobe_20220113_112922     3853    89     60               True          0          1
 3 pumpprobe_20220113_142740     3853    82     60               True          0          0
 4 pumpprobe_20220113_165125     3854   118     60               True          0          1
 5 pumpprobe_20220211_094256     1462   100     19               True          4          1
 6 pumpprobe_20220211_141039      683   116      5               True          5          0
 7 pumpprobe_20220211_151211      550   107      4               True          4         19
 8 pumpprobe_20220412_112143     4493   134     68               True          3         11
 9 pumpprobe_20220422_104335     3416   110     52               True          0          3
10 pumpprobe_20220415_104205     4599   150     69               True          4         12
11 pumpprobe_20220422_144057     4522   135     65               True          2          4
12 pumpprobe_20220505_111403     4571   131     68               True          0         11
13 pumpprobe_20220510_110807     4785   148     75               True          0         22
14 pumpprobe_20220510_162733     4748   160     75               True          0         37
15 pumpprobe_20220511_124404     4920   111     71               True          3          8
16 pumpprobe_20220511_150909     4448   135     67               True          3         23
17 pumpprobe_20220511_161734     4566   115     68               True          0          4

WT export recordings processed for the variant: 7

rid_u = pd.DataFrame(rid_rows(obs_u, "unc-31 export") + rid_rows(obs_wte, "WT export"))
print(rid_u.round(3).to_string(index=False))
wt_main = rid_wt[rid_wt.source == "dandiset, assignment here"].set_index("reading")
wt_exp = rid_u[rid_u.source == "WT export"].set_index("reading")
un = rid_u[rid_u.source == "unc-31 export"].set_index("reading")
cmp_ = pd.DataFrame({"WT dandiset mean": wt_main["mean"], "WT dandiset trials": wt_main.trials,
                     "WT export mean": wt_exp["mean"], "WT export trials": wt_exp.trials,
                     "unc-31 mean": un["mean"], "unc-31 trials": un.trials})
cmp_["reduced vs WT dandiset"] = np.where(cmp_[["WT dandiset mean", "unc-31 mean"]].notna().all(axis=1), cmp_["unc-31 mean"] < cmp_["WT dandiset mean"], np.nan)
cmp_["reduced vs WT export"] = np.where(cmp_[["WT export mean", "unc-31 mean"]].notna().all(axis=1), cmp_["unc-31 mean"] < cmp_["WT export mean"], np.nan)
def trials(o, names):
    y = o[(o.j == "RID") & o.i.isin(names) & (o.dff.abs() <= 2)]
    return y.dff.values
mw = {}
for cls, shown in FIG4.items():
    for reading, names in [("(a) %s, L and R" % cls, [cls + "L", cls + "R"]), ("(b) %s" % shown, [shown])]:
        a_, b_, e_ = trials(obs, names), trials(obs_u, names), trials(obs_wte, names)
        mw[reading] = {"p, unc-31 below WT dandiset": stats.mannwhitneyu(b_, a_, alternative="less").pvalue if len(a_) and len(b_) else np.nan,
                       "p, unc-31 below WT export": stats.mannwhitneyu(b_, e_, alternative="less").pvalue if len(e_) and len(b_) else np.nan}
cmp_ = cmp_.join(pd.DataFrame(mw).T)
print("\n" + cmp_.round(3).to_string())
print("(one-sided Mann-Whitney tests of the trials, for information; the criterion compares means)")
for (lab_, o_) in [("unc-31", obs_u), ("WT export", obs_wte)]:
    y = o_[(o_.j == "RID") & o_.i.str[:3].isin(list(FIG4)) & (o_.dff.abs() <= 2)]
    print("\n%s, per trial (recording, stimulation, observed neuron, <dF/F0>_t, detected):" % lab_)
    print(y[["subject", "stim", "i", "dff", "responded"]].round(3).to_string(index=False))

       source          reading neuron  rid_stims  trials  animals   mean  median  detected  positive_trials  pctile_among_head_pairs
unc-31 export (a) URX, L and R    URX         28      25        5  0.031   0.050         7               14                   66.819
unc-31 export         (b) URXL    URX         28      16        5  0.018   0.035         3                9                   60.116
unc-31 export (a) ADL, L and R    ADL         28      14        4  0.061   0.036         3                9                   78.205
unc-31 export         (b) ADLR    ADL         28       9        3  0.033   0.031         2                6                   67.980
unc-31 export (a) AWB, L and R    AWB         28      50        7  0.067   0.047        17               32                   79.783
unc-31 export         (b) AWBR    AWB         28      25        7  0.048   0.048         5               18                   74.328
    WT export (a) URX, L and R    URX         12       5        4  0.109   0.092         3                5                   87.745
    WT export         (b) URXL    URX         12       3        3  0.136   0.092         2                3                   90.810
    WT export (a) ADL, L and R    ADL         12       1        1 -0.062  -0.062         0                0                   16.656
    WT export         (b) ADLR    ADL         12       0        0    NaN     NaN         0                0                      NaN
    WT export (a) AWB, L and R    AWB         12      18        7  0.086   0.068         7               12                   83.849
    WT export         (b) AWBR    AWB         12      11        6  0.112   0.126         6                7                   88.051

                  WT dandiset mean  WT dandiset trials  WT export mean  WT export trials  unc-31 mean  unc-31 trials  reduced vs WT dandiset  reduced vs WT export  p, unc-31 below WT dandiset  p, unc-31 below WT export
reading                                                                                                                                                                                                                   
(a) URX, L and R             0.032                   1           0.109                 5        0.031             25                     1.0                   1.0                        0.538                      0.085
(b) URXL                       NaN                   0           0.136                 3        0.018             16                     NaN                   1.0                          NaN                      0.086
(a) ADL, L and R               NaN                   0          -0.062                 1        0.061             14                     NaN                   0.0                          NaN                      1.000
(b) ADLR                       NaN                   0             NaN                 0        0.033              9                     NaN                   NaN                          NaN                        NaN
(a) AWB, L and R             0.076                  12           0.086                18        0.067             50                     1.0                   1.0                        0.440                      0.268
(b) AWBR                     0.094                   7           0.112                11        0.048             25                     1.0                   1.0                        0.312                      0.129
(one-sided Mann-Whitney tests of the trials, for information; the criterion compares means)

unc-31, per trial (recording, stimulation, observed neuron, <dF/F0>_t, detected):
 subject  stim    i    dff  responded
 unc31-5     0 AWBR  0.005      False
 unc31-5     0 ADLR -0.058      False
 unc31-5     0 URXL -0.099      False
 unc31-5     0 AWBL -0.053      False
 unc31-5     1 AWBR -0.029      False
 unc31-5     1 ADLR  0.041      False
 unc31-5     1 URXL  0.122      False
 unc31-5     1 AWBL -0.008      False
 unc31-5     2 AWBR  0.081      False
 unc31-5     2 ADLR  0.097      False
 unc31-5     2 URXL -0.056      False
 unc31-5     2 AWBL  0.088      False
 unc31-5     3 AWBR -0.091      False
 unc31-5     3 ADLR -0.037      False
 unc31-5     3 URXL  0.020      False
 unc31-5     3 AWBL -0.066      False
 unc31-6     0 AWBR -0.031      False
 unc31-6     0  URX -0.070      False
 unc31-6     0 AWBL -0.060      False
 unc31-6     1 AWBR  0.011      False
 unc31-6     1  URX  0.050       True
 unc31-6     1 AWBL -0.027       True
 unc31-6     2 AWBR -0.034      False
 unc31-6     2  URX -0.033      False
 unc31-6     2 AWBL  0.084       True
 unc31-6     3 AWBR  0.008      False
 unc31-6     3  URX -0.007      False
 unc31-6     3 AWBL  0.011       True
 unc31-6     4 AWBR -0.081      False
 unc31-6     4  URX -0.079      False
 unc31-6     4 AWBL -0.053      False
 unc31-7     0 AWBR -0.092      False
 unc31-7     0 AWBL -0.022      False
 unc31-7     1 AWBR  0.069      False
 unc31-7     1 AWBL  0.140      False
 unc31-7     2 AWBR  0.078      False
 unc31-7     2 AWBL  0.261       True
 unc31-7     3 AWBR  0.082       True
 unc31-7     3 AWBL  0.212       True
 unc31-8     0 AWBL  0.005      False
 unc31-8     0 URXL -0.106      False
 unc31-8     0 ADLR -0.025      False
 unc31-8     0 AWBR  0.012      False
 unc31-8     1 AWBL -0.018      False
 unc31-8     1 URXL  0.051      False
 unc31-8     1 ADLR  0.031      False
 unc31-8     1 AWBR  0.045      False
 unc31-8     2 AWBL -0.021      False
 unc31-8     2 URXL -0.031      False
 unc31-8     2 ADLR  0.013      False
 unc31-8     2 AWBR  0.048      False
unc31-10     0 AWBR  0.065      False
unc31-10     0 URXR -0.001      False
unc31-10     0 URXL -0.044      False
unc31-10     0 AWBL  0.130       True
unc31-10     1 AWBR  0.091       True
unc31-10     1 URXR -0.044      False
unc31-10     1 URXL  0.059      False
unc31-10     1 AWBL  0.328       True
unc31-10     2 AWBR  0.165      False
unc31-10     2 URXR -0.016      False
unc31-10     2 URXL -0.026      False
unc31-10     2 AWBL  0.270       True
unc31-10     3 AWBR  0.260       True
unc31-10     3 URXR  0.080       True
unc31-10     3 URXL  0.050      False
unc31-10     3 AWBL  0.537      False
unc31-11     3 AWBL -0.118       True
unc31-11     3 ADLL -0.051      False
unc31-11     3 URXL  0.113       True
unc31-11     3 URXR  0.065       True
unc31-11     3 ADLR  0.042       True
unc31-11     3 AWBR -0.013      False
unc31-11     4 AWBL -0.028       True
unc31-11     4 ADLL  0.437       True
unc31-11     4 URXL  0.074       True
unc31-11     4 URXR  0.288       True
unc31-11     4 ADLR  0.194       True
unc31-11     4 AWBR  0.189      False
unc31-15     0 URXR -0.079      False
unc31-15     0 URXL -0.097      False
unc31-15     0 ADLL  0.109      False
unc31-15     1 URXR  0.099       True
unc31-15     1 URXL  0.079      False
unc31-15     1 ADLL  0.073      False
unc31-15     2 URXR  0.096      False
unc31-15     2 URXL  0.179       True
unc31-15     2 ADLL -0.010      False
unc31-16     0 AWBR  0.080      False
unc31-16     0 AWBL  0.182       True
unc31-16     1 AWBR  0.139       True
unc31-16     1 AWBL  0.194       True
unc31-16     2 AWBR  0.150       True
unc31-16     2 AWBL  0.168      False

WT export, per trial (recording, stimulation, observed neuron, <dF/F0>_t, detected):
     subject  stim    i    dff  responded
wt-export-25     0 AWBL -0.107      False
wt-export-25     0 AWBR -0.076      False
wt-export-25     1 AWBL  0.093      False
wt-export-25     1 AWBR -0.007      False
wt-export-25     2 AWBL  0.044      False
wt-export-25     2 AWBR -0.018      False
wt-export-26    49 AWBL -0.104      False
wt-export-26    49  ADL  0.009       True
wt-export-41    50 AWBR  0.299       True
wt-export-41    50 URXR  0.033      False
wt-export-41    50 AWBL  0.291       True
wt-export-42     0 AWBR  0.262       True
wt-export-42     1 AWBR  0.236       True
wt-export-42     2 AWBR  0.043       True
wt-export-42     3 AWBR  0.243       True
wt-export-76    22 URXL  0.048      False
wt-export-76    22 AWBL  0.003      False
wt-export-76    22 ADLL -0.062      False
wt-export-76    22 AWBR -0.020      False
wt-export-84    19 AWBR  0.126      False
wt-export-84    19 URXR  0.104       True
wt-export-84    19 AWBL  0.099      False
wt-export-84    19  ADL  0.131      False
wt-export-84    19 URXL  0.269       True
wt-export-89    19 URXL  0.092       True
wt-export-89    19 AWBR  0.139       True

# Build the paper's kernel-fitting package and the packages it imports, at fixed commits, in a
# temporary environment with this kernel's Python and numpy (own uv cache, so nothing is reused).
import shutil, subprocess, sys, tempfile
UV = shutil.which("uv")
assert UV, "uv is needed for this check"
PKGS = [("pumpprobe", "https://github.com/leiferlab/pumpprobe", "1dbc5e0a2b609d54bc9b1c90c73d4e3bf183d3c7"),
        ("mistofrutta", "https://github.com/francescorandi/mistofrutta", "11ac8c3e40743b6e35f54568e2de36b95466bbc4"),
        ("wormbrain", "https://github.com/leiferlab/wormbrain", "d94f7a1456622f77eebcc181de1d8504618764c4"),
        ("wormdatamodel", "https://github.com/leiferlab/wormdatamodel", "2ab956199e3931de41a190d2b9985e961df3810c")]
with tempfile.TemporaryDirectory() as tmp:
    env0 = {**os.environ, "UV_CACHE_DIR": tmp + "/cache"}
    envpy = tmp + "/env/bin/python"
    subprocess.run([UV, "venv", "-q", "-p", "%d.%d" % sys.version_info[:2], tmp + "/env"], check=True, env=env0)
    subprocess.run([UV, "pip", "install", "-q", "-p", envpy, "numpy==" + np.__version__, "setuptools", "gitpython"],
                   check=True, env=env0)
    inc = subprocess.run([envpy, "-c", "import numpy; print(numpy.get_include())"],
                         capture_output=True, text=True, check=True).stdout.strip()
    print("Python %d.%d, numpy %s" % (*sys.version_info[:2], np.__version__))
    def build(name, url, commit, env):
        p = subprocess.run([UV, "pip", "install", "-p", envpy, "--no-build-isolation", "--no-deps",
                            "%s @ git+%s@%s" % (name, url, commit)], capture_output=True, text=True, env=env)
        if p.returncode == 0:
            return "built"
        err = [l.strip() for l in p.stderr.splitlines()
               if re.search(r"error:|Invalid version", l) and not re.search(r"Failed to build|build backend|error: Command", l)]
        return "failed: " + (err[0] if err else p.stderr.strip().splitlines()[-1])
    print("pumpprobe %s, as published: %s" % (PKGS[0][2][:7], build(*PKGS[0], env0)))
    with_inc = {**env0, "CPATH": inc}
    for name, url, commit in PKGS:
        print("%s %s, numpy headers added by hand: %s" % (name, commit[:7], build(name, url, commit, with_inc)))
    p = subprocess.run([envpy, "-c", "import pumpprobe"], capture_output=True, text=True)
    print("import pumpprobe:", "ok" if p.returncode == 0 else p.stderr.strip().splitlines()[-1])

Python 3.13, numpy 2.5.4

pumpprobe 1dbc5e0, as published: failed: pumpprobe/_integration.cpp:3:10: fatal error: numpy/arrayobject.h: No such file or directory

pumpprobe 1dbc5e0, numpy headers added by hand: built

mistofrutta 11ac8c3, numpy headers added by hand: failed: raise InvalidVersion(f"Invalid version: {version!r}")

wormbrain d94f7a1, numpy headers added by hand: failed: dsmm/dsmm.cpp:6:10: fatal error: Eigen/Dense: No such file or directory

wormdatamodel 2ab9561, numpy headers added by hand: failed: wormdatamodel/data/_legacy_c.cpp:63:65: error: ‘NPY_IN_ARRAY’ was not declared in this scope
import pumpprobe: ModuleNotFoundError: No module named 'mistofrutta'

print(f"Total run time {time.time() - started:.0f} s")

Total run time 1045 s

Claim	Data	Tested
recordings-protocol: 113 WT animals, up to 40 min, a stimulation every 30 s	dandiset	yes
pairs-measured: 23,433 pairs, 186 of 188 neurons, up to 59 observations	dandiset	yes
fig1e-avjr-avdr: AVJR stimulation evokes AVDR activity	dandiset	yes
fig1f-aver-avar: AVER stimulation evokes AVAR activity	dandiset	yes
fig1g-saadl-ollr: SAADL stimulation inhibits OLLR	dandiset	yes
edfig2ij-spatial-specificity: neighbours of the target mostly do not respond	dandiset	yes
sourcedata-link: per-pair counts and means agree with the Source Data	with inputs	yes
fig2d-inhibitory-fraction: 11% of q < 0.05 connections are inhibitory	with inputs	yes
osf-link: recordings and stimulated neurons link to the OSF export	with labels	yes
fig2c-1310-connected, fig2c-bilateral, fig2efg-path-length	dandiset	blocked: data missing (control recordings)
fig3-anatomy-predictions	dandiset	blocked: data missing (model predictions)
edfig6-kernel-stereotypy	dandiset	blocked: not runnable
edfig7b-unc31-fewer, fig5-extrasynaptic-screen	dandiset	blocked: data missing
fig4-rid-wt: URX, ADL and AWB respond to RID stimulation	dandiset	yes
fig4-rid-unc31: these responses are reduced in unc-31	elsewhere (OSF)	yes
fig6-spontaneous	elsewhere (OSF)	blocked: not runnable (needs the kernels of edfig6)

Claim	Data	Outcome
recordings-protocol	dandiset	partially reproduced (110 of 113 animals have traces)
pairs-measured	dandiset	partially reproduced (20,587 pairs, 184 neurons, at most 50 observations)
fig1e-avjr-avdr	dandiset	reproduced
fig1f-aver-avar	dandiset	reproduced (7 trials, 5 animals)
fig1g-saadl-ollr	dandiset	reproduced (5 trials, 4 animals; sign not robust)
edfig2ij-spatial-specificity	dandiset	reproduced
sourcedata-link	with inputs	partially reproduced
fig2d-inhibitory-fraction	with inputs	not reproduced (15.1% against 11%)
osf-link	with labels	partially reproduced (66% of position-found targets agree)
fig2c-1310-connected, fig2c-bilateral, fig2efg-path-length	dandiset	blocked: data missing
fig3-anatomy-predictions	dandiset	blocked: data missing
edfig6-kernel-stereotypy	dandiset	blocked: not runnable
edfig7b-unc31-fewer, fig5-extrasynaptic-screen	dandiset	blocked: data missing
fig4-rid-wt	dandiset	partially reproduced (AWB only; URX seen once, ADL never)
fig4-rid-unc31	elsewhere	reproduced (small differences, few WT trials)
fig6-spontaneous	elsewhere	blocked: not runnable (needs the kernels of edfig6)

Reproducing Randi et al. (2023) from dandiset 001075¶

Data check¶

Observations¶

recordings-protocol¶

pairs-measured¶

Source Data and the per-pair link¶

osf-link¶

sourcedata-link¶

fig1e-avjr-avdr¶

fig1f-aver-avar¶

fig1g-saadl-ollr¶

edfig2ij-spatial-specificity¶

fig2d-inhibitory-fraction¶

fig4-rid-wt¶

fig4-rid-unc31¶

edfig6-kernel-stereotypy (blocked)¶

Summary¶