import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import warnings
import collections
import numpy as np
import pandas as pd
from scipy import stats, signal
import matplotlib.pyplot as plt
import neuropal as npl

T_START = time.time()
warnings.filterwarnings("ignore", category=RuntimeWarning)  # empty means over animals lacking a pair
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
print("dandiset", npl.DANDISET_ID, "version", npl.DANDISET_VERSION)

dandiset 000541 version 0.241009.1457

assets = npl.list_assets()
print(len(assets), "assets;", round(sum(a["size"] for a in assets) / 1e9, 2), "GB")
recs = npl.load_all()
recs = sorted(recs, key=lambda r: r["identifier"])
print("loaded in", round(time.time() - T_START), "s")

21 assets; 30.53 GB

loaded in 118 s

rows = []
for r in recs:
    named = [n for n in r["labels"] if n]
    allnan = int(np.sum(np.all(~np.isfinite(r["raw"]), axis=0)))
    order = sorted(r["stimuli"], key=lambda s: r["stimuli"][s][0])
    rows.append({
        "file": r["path"].split("/")[-1], "identifier": r["identifier"], "strain": r["strain"],
        "stage": r["growth_stage"], "description": r["session_description"],
        "frames": r["raw"].shape[0], "movie_frames": r["n_movie_frames"], "rate_Hz": r["rate"],
        "duration_s": r["raw"].shape[0] / r["rate"], "rois": r["raw"].shape[1], "named": len(named),
        "unnamed": r["raw"].shape[1] - len(named), "dup_names": len(named) - len(set(named)),
        "all_nan_traces": allnan, "two_series_equal": r["raw_equals_other_series"],
        "stim_order": "/".join(s[:4] for s in order),
        "stim_times": ";".join(f"{r['stimuli'][s][0]}-{r['stimuli'][s][1]}" for s in order),
        "raw_min": float(np.nanmin(r["raw"])), "raw_max": float(np.nanmax(r["raw"])),
    })
sanity = pd.DataFrame(rows)
print(sanity.drop(columns=["file"]).to_string())
print()
print("files:", len(recs), "| strains:", list(sanity.strain.unique()), "| stages:", list(sanity.stage.unique()),
      "| descriptions:", list(sanity.description.unique()))
print("trace unit:", {r["unit"] for r in recs}, "| rate:", {r["rate"] for r in recs},
      "| starting time:", {r["t0"] for r in recs})
print("grid spacing (um):", {tuple(np.round(r["grid_spacing_um"], 4)) for r in recs})
print("frames per trace: %d to %d (%.1f to %.1f s at 4 Hz)" % (sanity.frames.min(), sanity.frames.max(),
      sanity.frames.min() / 4, sanity.frames.max() / 4))
print("trace frames = movie frames - 1 in every file:", bool((sanity.frames == sanity.movie_frames - 1).all()))
print("stimulus orders:", dict(collections.Counter(sanity.stim_order)))
print("stimulus windows:", list(sanity.stim_times.unique()))
print("stimulus table description:", {r["stim_description"] for r in recs})
print("the dNMF series equals the 'raw' series (NaN equal to NaN) in", int(sanity.two_series_equal.sum()),
      "of", len(recs), "files")
print("raw trace values: %.3f to %.3f (unitless)" % (sanity.raw_min.min(), sanity.raw_max.max()))

     identifier   strain stage                                   description  frames  movie_frames  rate_Hz  duration_s  rois  named  unnamed  dup_names  all_nan_traces  two_series_equal      stim_order                         stim_times   raw_min   raw_max
0   20190924_01  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   181    180        1          0               0              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.042480  3.887675
1   20190924_02  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   178    176        2          0               0              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.018610  4.197489
2   20190924_03  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     962           963      4.0      240.50   188    184        4          0               0              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.052671  3.450128
3   20190925_01  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     962           963      4.0      240.50   184    180        4          0               1              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.016917  3.346579
4   20190925_04  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     960           961      4.0      240.00   177    176        1          0               2              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.066597  2.707961
5   20190928_01  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   172    169        3          0               0              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.024956  2.610775
6   20190928_03  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   168    167        1          0               0              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.022305  3.683224
7   20190928_05  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   172    164        8          0               0              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.013182  2.592344
8   20190928_07  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     959           960      4.0      239.75   180    179        1          0               2              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.012988  3.904508
9   20190928_08  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     960           961      4.0      240.00   183    183        0          0               1              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.033154  2.316842
10  20190928_11  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     962           963      4.0      240.50   181    177        4          0               2              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.008192  4.099419
11  20190928_13  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   175    175        0          0               0              True  nacl/pent/buta  60.5-70.5;120.5-130.5;180.5-190.5  0.031965  2.536674
12  20190929_02  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   166    164        2          0               0              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.011046  3.383612
13  20190929_03  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   178    175        3          0               0              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.018072  3.224996
14  20190929_05  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   170    168        2          0               0              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.019185  3.251345
15  20190929_06  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   178    175        3          0               0              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.011302  3.321618
16  20190929_07  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     936           937      4.0      234.00   177    176        1          0               0              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.035262  2.745039
17  20191030_03  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   185    182        3          0               2              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.029313  3.400185
18  20191030_07  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     960           961      4.0      240.00   180    175        5          0               3              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.005829  2.209842
19  20191104_08  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   179    177        2          0               0              True  buta/nacl/pent  60.5-70.5;120.5-130.5;180.5-190.5  0.125432  3.199358
20  20191104_10  OH16230    YA  C. elegans head NeuroPAL and Calcium imaging     961           962      4.0      240.25   174    174        0          0               0              True  pent/buta/nacl  60.5-70.5;120.5-130.5;180.5-190.5  0.033234  2.487772

files: 21 | strains: ['OH16230'] | stages: ['YA'] | descriptions: ['C. elegans head NeuroPAL and Calcium imaging']
trace unit: {'unitless'} | rate: {4.0} | starting time: {0.0}
grid spacing (um): {(np.float64(0.2667), np.float64(0.2667), np.float64(1.5))}
frames per trace: 936 to 962 (234.0 to 240.5 s at 4 Hz)
trace frames = movie frames - 1 in every file: True
stimulus orders: {'pent/buta/nacl': 7, 'nacl/pent/buta': 7, 'buta/nacl/pent': 7}
stimulus windows: ['60.5-70.5;120.5-130.5;180.5-190.5']
stimulus table description: {'Chemical stimuli delivered to  the worm. Each worms gets 160 mM NaCl, 10^-4mM 2-butanone, and 10^-4 2,3-pentanedione in a random order. Start and stop time represent seconds from start of the experiment'}
the dNMF series equals the 'raw' series (NaN equal to NaN) in 21 of 21 files
raw trace values: 0.006 to 4.197 (unitless)

# The calcium volume's stated grid spacing against its coordinates: the same named neurons in the
# NeuroPAL volume and in the first frame of the calcium segmentation.
sc = pd.DataFrame([{"file": r["identifier"], "x": r["xyz_scale"][0], "y": r["xyz_scale"][1], "z": r["xyz_scale"][2]}
                   for r in recs])
print("NeuroPAL voxel coordinate per calcium voxel coordinate, same neurons (slope per axis), over 21 files:")
print(sc[["x", "y", "z"]].describe().loc[["min", "50%", "max"]].round(3).to_string())
print("NeuroPAL volume shape:", {r["neuropal_shape"] for r in recs}, "| calcium movie shape (t, x, y, z):",
      {r["movie_shape"][1:] for r in recs})
print("stated grid spacing (um), calcium:", {tuple(np.round(r["grid_spacing_um"], 4)) for r in recs},
      "| NeuroPAL:", {tuple(np.round(r["neuropal_grid_spacing_um"], 4)) for r in recs})
print("calcium spacing used here (NeuroPAL spacing times the rounded slope):",
      {tuple(np.round(r["neuropal_grid_spacing_um"] * np.round(r["xyz_scale"]), 4)) for r in recs})
print("extent of the named neurons (um, corrected), x/y/z median over files:",
      np.round(np.median([np.ptp(r["centre_um"], axis=0) for r in recs], axis=0), 1))

NeuroPAL voxel coordinate per calcium voxel coordinate, same neurons (slope per axis), over 21 files:
         x      y      z
min  1.951  1.937  0.970
50%  1.999  2.015  1.033
max  2.046  2.088  1.072
NeuroPAL volume shape: {(512, 256, 21, 5)} | calcium movie shape (t, x, y, z): {(256, 128, 21)}
stated grid spacing (um), calcium: {(np.float64(0.2667), np.float64(0.2667), np.float64(1.5))} | NeuroPAL: {(np.float64(0.2667), np.float64(0.2667), np.float64(1.5))}
calcium spacing used here (NeuroPAL spacing times the rounded slope): {(np.float64(0.5333), np.float64(0.5333), np.float64(1.5))}
extent of the named neurons (um, corrected), x/y/z median over files: [130.7  49.1  21. ]

# Recordings deposited more than once: identical trace matrices, or any trace identical across files.
fp = {}
for r in recs:
    x = np.nan_to_num(r["raw"])
    fp[r["identifier"]] = x
same = []
ids = list(fp)
for i in range(len(ids)):
    for j in range(i + 1, len(ids)):
        a, b = fp[ids[i]], fp[ids[j]]
        n = min(len(a), len(b))
        ca = a[:n].T[np.std(a[:n], axis=0) > 0]
        cb = b[:n].T[np.std(b[:n], axis=0) > 0]
        za = (ca - ca.mean(1, keepdims=True)) / ca.std(1, keepdims=True)
        zb = (cb - cb.mean(1, keepdims=True)) / cb.std(1, keepdims=True)
        cmax = (za @ zb.T / n).max()
        if cmax > 0.999:
            same.append((ids[i], ids[j], cmax))
print("pairs of files sharing a trace correlated above 0.999:", same if same else "none")

names = sorted({n for r in recs for n in r["labels"] if n})
print("distinct neuron names over the 21 heads:", len(names))
cnt = collections.Counter(n for r in recs for n in set(r["labels"]) if n)
print("names in 3 or fewer heads:", {n: c for n, c in cnt.items() if c <= 3})
nonneuronal = [n for n in names if n in ("AMSOL", "AMSOR", "HMC")]
print("names that are not neurons (amphid sheath glia AMSO, head mesodermal cell HMC):", nonneuronal,
      {n: cnt[n] for n in nonneuronal})
print("DD6 (a motor neuron of the posterior ventral cord) named in", cnt["DD6"], "head file(s):",
      [r["identifier"] for r in recs if "DD6" in r["labels"]])

pairs of files sharing a trace correlated above 0.999: none
distinct neuron names over the 21 heads: 193
names in 3 or fewer heads: {'SABD': 3, 'RIGR': 3, 'DA1': 2, 'VD2': 2, 'AS1': 2, 'DD6': 1, 'AMSOR': 3}
names that are not neurons (amphid sheath glia AMSO, head mesodermal cell HMC): ['AMSOL', 'AMSOR', 'HMC'] {'AMSOL': 4, 'AMSOR': 3, 'HMC': 6}
DD6 (a motor neuron of the posterior ventral cord) named in 1 head file(s): ['20190929_06']

Z = npl.load_zenodo_heads()          # the paper's deposit, for the "deposit" reading and the link check
DEP = npl.deposit_slots(Z)
DFF = {r["identifier"]: npl.dff(r["raw"], r["rate"]) for r in recs}
MODES = ("one", "deposit", "both", "none")   # "one" is the main reading
data = {mode: npl.prepare(recs, mode, DEP, DFF) for mode in MODES}
D = data["one"]
cls_of = D["class_of"]
for mode, d in data.items():
    kept = [set(a["traces"]) for a in d["animals"]]
    print(f"reading '{mode}': removed {sum(len(a['duplicates']) for a in d['animals'])} named traces;"
          f" kept {sum(len(k) for k in kept)} named traces;"
          f" traces dropped as all NaN: {sum(len(a['dropped_nan']) for a in d['animals'])}")
allvals = np.concatenate([np.concatenate([t[4:] for t in a["traces"].values()]) for a in D["animals"]])
print("dF/F0 over all kept traces (main reading): %.2f to %.2f, median %.3f; values beyond +-10: %d of %d"
      % (allvals.min(), allvals.max(), np.median(allvals), np.sum(np.abs(allvals) > 10), allvals.size))
print("classes (left and right pooled):", len(set(cls_of.values())))

reading 'one': removed 338 named traces; kept 3325 named traces; traces dropped as all NaN: 13
reading 'deposit': removed 126 named traces; kept 3550 named traces; traces dropped as all NaN: 0
reading 'both': removed 651 named traces; kept 3012 named traces; traces dropped as all NaN: 13
reading 'none': removed 0 named traces; kept 3663 named traces; traces dropped as all NaN: 13
dF/F0 over all kept traces (main reading): -41.59 to 89.48, median -0.010; values beyond +-10: 4955 of 3174325
classes (left and right pooled): 112

# Per file: nearest-neighbour distance between named neurons, pairs within 5 um, pairs within 5 um
# at r >= 0.95, traces each reading removes, and how many of those the deposit also lacks.
rows = []
for r in recs:
    d = DFF[r["identifier"]]
    idx = [i for i, n in enumerate(r["labels"]) if n and np.all(np.isfinite(d[4:, i]))]
    xyz = r["centre_um"][idx]
    dist = np.linalg.norm(xyz[:, None] - xyz[None], axis=-1)
    np.fill_diagonal(dist, np.inf)
    cc = np.corrcoef(d[4:, idx].T)
    iu = np.triu_indices(len(idx), 1)
    named = {i for i, n in enumerate(r["labels"]) if n}
    lacking = {i for i in named if (r["identifier"], r["labels"][i]) not in DEP}
    row = {"file": r["identifier"], "named": len(named), "nn_um_median": np.median(dist.min(1)),
           "nn_um_p5": np.percentile(dist.min(1), 5), "pairs_5um": int((dist[iu] <= 5).sum()),
           "pairs_5um_r95": int(((dist[iu] <= 5) & (cc[iu] >= 0.95)).sum()),
           "pairs_any_r95": int((cc[iu] >= 0.95).sum()), "deposit_lacks": len(lacking)}
    for mode in ("both", "one"):
        rm = npl.near_duplicates(r, d, mode)
        row[f"rm_{mode}"] = len(rm)
        row[f"rm_{mode}_and_lacking"] = len(rm & lacking)
    rows.append(row)
xt = pd.DataFrame(rows)
print(xt.round(2).to_string(index=False))
tot = xt.drop(columns=["file", "nn_um_median", "nn_um_p5"]).sum()
print("totals:", tot.to_dict())
lack_names = collections.Counter(n for r in recs for n in r["labels"] if n and (r["identifier"], n) not in DEP)
print("named traces the deposit lacks, by name (top 12):", dict(lack_names.most_common(12)))
print("of these, all NaN in the dandiset:", sum(1 for r in recs for i, n in enumerate(r["labels"])
      if n and (r["identifier"], n) not in DEP and not np.all(np.isfinite(DFF[r["identifier"]][4:, i]))))
# How close and how correlated are the traces the deposit lacks, against their nearest named neighbour?
near = []
for r in recs:
    d = DFF[r["identifier"]]
    idx = [i for i, n in enumerate(r["labels"]) if n and np.all(np.isfinite(d[4:, i]))]
    cc = np.corrcoef(d[4:, idx].T)
    xyz = r["centre_um"][idx]
    dist = np.linalg.norm(xyz[:, None] - xyz[None], axis=-1)
    np.fill_diagonal(dist, np.inf)
    for k, i in enumerate(idx):
        if (r["identifier"], r["labels"][i]) not in DEP:
            j = int(np.argmin(dist[k]))
            near.append({"dist_um": dist[k, j], "r": cc[k, j], "max_r_within_5um": np.max(np.where(dist[k] <= 5, cc[k], -1))})
near = pd.DataFrame(near)
print("finite traces the deposit lacks: %d; nearest neighbour within 5 um: %d; with r >= 0.95 to a neuron within 5 um: %d"
      % (len(near), int((near.dist_um <= 5).sum()), int((near.max_r_within_5um >= 0.95).sum())))

       file  named  nn_um_median  nn_um_p5  pairs_5um  pairs_5um_r95  pairs_any_r95  deposit_lacks  rm_both  rm_both_and_lacking  rm_one  rm_one_and_lacking
20190924_01    180          3.06      1.51        246             27             38              3       33                    1      19                   0
20190924_02    176          3.06      1.19        219             20             35              6       34                    1      18                   1
20190924_03    184          2.72      1.19        279             35             91             10       53                    2      28                   2
20190925_01    180          2.87      1.19        210             29             72              6       50                    2      24                   2
20190925_04    176          2.67      1.19        325             34             73              5       37                    0      20                   0
20190928_01    169          3.24      1.19        239             19             32              4       36                    0      18                   0
20190928_03    167          3.45      1.51        193             18             47              4       26                    2      14                   1
20190928_05    164          3.43      1.51        171             11             24              3       20                    1      10                   1
20190928_07    179          3.24      1.50        200             19             42              6       32                    1      17                   1
20190928_08    183          2.77      1.19        320             26             65              4       35                    1      18                   0
20190928_11    177          3.06      1.60        210             11             23              4       16                    1       9                   1
20190928_13    175          3.06      1.07        249             19             38              3       26                    1      13                   0
20190929_02    164          3.45      1.59        155             12             12             18       22                    3      12                   3
20190929_03    175          3.11      1.19        195             19             40              4       29                    0      15                   0
20190929_05    168          2.87      1.50        211             17             22              5       34                    3      17                   2
20190929_06    175          3.11      1.59        172             14             16              6       27                    1      14                   1
20190929_07    176          2.82      1.19        264             15             16              3       29                    2      14                   2
20191030_03    182          3.08      1.50        301              9             22              9       17                    3       8                   1
20191030_07    175          3.44      1.19        147              9             24             10       17                    0       8                   0
20191104_08    177          2.82      0.53        285             37             69              7       48                    1      26                   1
20191104_10    174          3.06      1.07        250             19             33              6       30                    1      16                   1
totals: {'named': 3676, 'pairs_5um': 4841, 'pairs_5um_r95': 419, 'pairs_any_r95': 834, 'deposit_lacks': 126, 'rm_both': 651, 'rm_both_and_lacking': 27, 'rm_one': 338, 'rm_one_and_lacking': 20}
named traces the deposit lacks, by name (top 12): {'AVL': 11, 'I5': 8, 'HMC': 6, 'AMSOL': 4, 'RIMR': 4, 'DB2': 4, 'SABVL': 4, 'RMHR': 3, 'VA1': 3, 'SIBDL': 3, 'I1L': 3, 'I1R': 3}
of these, all NaN in the dandiset: 13
finite traces the deposit lacks: 113; nearest neighbour within 5 um: 79; with r >= 0.95 to a neuron within 5 um: 27

def resp_at(a, n, s0, s1):
    dur, pre = npl.window_means(a["traces"][n], s0, s1, a["rate"])
    return dur - pre

clock = []
for cls, stims, sign in (("ASE", ["nacl"], 1), ("AWC", ["butanone", "pentanedione"], -1)):
    for s in stims:
        true, shifted = [], []
        for a in D["animals"]:
            for n in (cls + "L", cls + "R"):
                if n not in a["traces"]:
                    continue
                true.append(resp_at(a, n, *a["stimuli"][s]))
                for o in npl.STIMULI:
                    if o != s:
                        shifted.append(resp_at(a, n, *a["stimuli"][o]))
        clock.append({"neuron": cls, "stimulus": s, "true_mean": np.mean(true), "n_true": len(true),
                      "other_windows_mean": np.mean(shifted), "n_other": len(shifted),
                      "mannwhitney_p": stats.mannwhitneyu(true, shifted).pvalue})
print(pd.DataFrame(clock).round(4).to_string())
# The ASE NaCl response with the stimulus table shifted in time: a response locked to the table peaks at 0 s.
shift_rows = []
for sh in np.arange(-20, 20.1, 2.5):
    v = [resp_at(a, n, a["stimuli"]["nacl"][0] + sh, a["stimuli"]["nacl"][1] + sh)
         for a in D["animals"] for n in ("ASEL", "ASER") if n in a["traces"]]
    shift_rows.append({"shift_s": sh, "ASE_nacl_mean_response": np.mean(v), "n": len(v)})
print(pd.DataFrame(shift_rows).round(3).to_string(index=False))

  neuron      stimulus  true_mean  n_true  other_windows_mean  n_other  mannwhitney_p
0    ASE          nacl     1.1937      40              0.3628       80         0.0000
1    AWC      butanone    -0.5924      42             -0.4794       84         0.2899
2    AWC  pentanedione    -0.4772      42             -0.5370       84         0.7075
 shift_s  ASE_nacl_mean_response  n
   -20.0                  -0.043 40
   -17.5                  -0.044 40
   -15.0                  -0.038 40
   -12.5                  -0.021 40
   -10.0                   0.007 40
    -7.5                   0.176 40
    -5.0                   0.580 40
    -2.5                   0.940 40
     0.0                   1.194 40
     2.5                   1.069 40
     5.0                   0.363 40
     7.5                  -0.309 40
    10.0                  -0.814 40
    12.5                  -0.999 40
    15.0                  -0.771 40
    17.5                  -0.515 40
    20.0                  -0.313 40

zfiles = ["2019" + str(f) for f in Z["files"]]
our_ids = [a["id"] for a in data["none"]["animals"]]
print("deposit files:", len(zfiles), "| matched to dandiset identifiers:", sum(f in our_ids for f in zfiles),
      "| same order:", zfiles == our_ids)
znames = [str(s) for s in Z["stim_names"]]
zorder = {f: [znames[k - 1] for k in np.atleast_1d(Z["stims"][i])] for i, f in enumerate(zfiles)}
ourorder = {a["id"]: sorted(a["stimuli"], key=lambda s: a["stimuli"][s][0]) for a in data["none"]["animals"]}
print("stimulus order agrees in", sum(zorder[f] == ourorder[f] for f in zfiles), "of", len(zfiles), "files")
print("deposit stimulus windows:", np.asarray(Z["stim_times"]).tolist(), "| fps:", Z["fps"])

zt = npl.zenodo_traces_by_file(Z)
zt_slot = DEP
NONE = {a["id"]: a for a in data["none"]["animals"]}
link_rows, count_mismatch = [], []
for cls, (trs, is_l) in zt.items():
    if np.all(np.isnan(is_l)):
        groups = [(cls, list(range(len(trs))))]
    else:
        groups = [(cls + "L", [k for k in range(len(trs)) if is_l[k] == 1]),
                  (cls + "R", [k for k in range(len(trs)) if is_l[k] == 0])]
    for name, ks in groups:
        if len(ks) != len(our_ids):
            count_mismatch.append((name, len(ks)))
            continue
        for f, k in zip(our_ids, ks):
            z = np.asarray(trs[k], float)
            in_dep = bool(np.isfinite(z).any())
            in_dan = name in NONE[f]["traces"]
            row = {"file": f, "neuron": name, "in_deposit": in_dep, "in_dandiset": in_dan, "r": np.nan,
                   "slope": np.nan, "deposit_len": len(z), "nan_head": -1}
            if in_dep and in_dan:
                ours = NONE[f]["traces"][name]
                nn = min(len(z), len(ours))
                z, ours = z[:nn], ours[:nn]
                ok = np.isfinite(z) & np.isfinite(ours)
                row.update(r=np.corrcoef(z[ok], ours[ok])[0, 1], slope=np.polyfit(ours[ok], z[ok], 1)[0],
                           nan_head=int(np.argmax(np.isfinite(z))))
            link_rows.append(row)
link = pd.DataFrame(link_rows)
deposit_names = sum(1 if np.all(np.isnan(il)) else int((il == 1).any()) + int((il == 0).any()) for _, il in zt.values())
print("deposit neuron names:", deposit_names, "| classes:", len(zt))
print("names without one slot per file in the deposit:", count_mismatch if count_mismatch else "none")
allnamed = {(a["id"], n) for a in NONE.values() for n in a["traces"]}
print("named dandiset traces:", len(allnamed), "| deposit slots with a trace:", int(link.in_deposit.sum()))
print("in both:", int((link.in_deposit & link.in_dandiset).sum()),
      "| in the deposit only:", int((link.in_deposit & ~link.in_dandiset).sum()),
      "| in the dandiset only:", int((~link.in_deposit & link.in_dandiset).sum()),
      "| dandiset traces with a name the deposit lacks:", sorted({n for f, n in allnamed} - set(link.neuron)))
only_dan = link[~link.in_deposit & link.in_dandiset]
print("dandiset-only traces by name:", only_dan.neuron.value_counts().to_dict())
link = link[link.in_deposit & link.in_dandiset].copy()
print("Pearson r between dandiset dF/F0 and deposit trace: median %.4f, 5th percentile %.4f, min %.4f"
      % (link.r.median(), link.r.quantile(0.05), link.r.min()))
print("fraction with r >= 0.9: %.3f; r >= 0.95: %.3f; r >= 0.99: %.3f"
      % ((link.r >= 0.9).mean(), (link.r >= 0.95).mean(), (link.r >= 0.99).mean()))
print("regression slope deposit on dandiset: median %.3f (IQR %.3f to %.3f)"
      % (link.slope.median(), link.slope.quantile(0.25), link.slope.quantile(0.75)))
link["len_diff"] = [row.deposit_len - len(NONE[row.file]["traces"][row.neuron]) for row in link.itertuples()]
print("deposit trace length minus dandiset frames (after dropping the last frame):",
      link.groupby("file").len_diff.first().value_counts().to_dict(),
    "| leading NaN frames in the deposit:", sorted(link.nan_head.unique()))
# Frame alignment: the deposit's trace against the dandiset's, shifted by -1 to +1 frames.
lag_r = collections.defaultdict(list)
for row in link.itertuples():
    z = np.asarray(zt_slot[(row.file, row.neuron)], float)
    ours = NONE[row.file]["traces"][row.neuron]
    for lag in (-1, 0, 1):
        a, b = z[max(lag, 0):], ours[max(-lag, 0):]
        nn = min(len(a), len(b))
        ok = np.isfinite(a[:nn]) & np.isfinite(b[:nn])
        lag_r[(row.file, lag)].append(np.corrcoef(a[:nn][ok], b[:nn][ok])[0, 1])
lag_tab = pd.DataFrame({lag: {f: np.median(lag_r[(f, lag)]) for f in our_ids} for lag in (-1, 0, 1)})
print("median r per file with the deposit shifted by -1, 0, +1 frames; best lag per file:",
      dict(collections.Counter(lag_tab.idxmax(axis=1))))
r1 = np.concatenate([lag_r[(f, 1)] for f in our_ids])
print("with deposit frame k + 1 against dandiset frame k: median r %.4f, fraction with r >= 0.9: %.3f"
      % (np.median(r1), np.mean(r1 >= 0.9)))
print(link.groupby("file").r.agg(["count", "median", "min"]).join(lag_tab.rename(columns=lambda c: f"r_lag{c:+d}")).round(4).to_string())

deposit files: 21 | matched to dandiset identifiers: 21 | same order: True
stimulus order agrees in 21 of 21 files
deposit stimulus windows: [[60.5, 70.5], [120.5, 130.5], [180.5, 190.5]] | fps: 4

deposit neuron names: 189 | classes: 109
names without one slot per file in the deposit: none
named dandiset traces: 3663 | deposit slots with a trace: 3550
in both: 3550 | in the deposit only: 0 | in the dandiset only: 101 | dandiset traces with a name the deposit lacks: ['AMSOL', 'AMSOR', 'DD6', 'HMC']
dandiset-only traces by name: {'AVL': 11, 'I5': 8, 'DB2': 4, 'RIMR': 4, 'SABVL': 4, 'AIAR': 3, 'RMHR': 3, 'SIBDL': 3, 'VA1': 3, 'AIAL': 2, 'ASJL': 2, 'AVBR': 2, 'RIML': 2, 'RIR': 2, 'RMGR': 2, 'SAADR': 2, 'VD1': 2, 'ADAR': 1, 'AIBL': 1, 'AINL': 1, 'AIYL': 1, 'AIYR': 1, 'AIZL': 1, 'ASHR': 1, 'ASIL': 1, 'ASIR': 1, 'ASJR': 1, 'AVBL': 1, 'AVDL': 1, 'AVFL': 1, 'AVFR': 1, 'AVG': 1, 'AVJL': 1, 'CEPDL': 1, 'FLPL': 1, 'FLPR': 1, 'I1L': 1, 'I4': 1, 'M3L': 1, 'M3R': 1, 'RIAR': 1, 'RIBR': 1, 'RIFL': 1, 'RIFR': 1, 'RIH': 1, 'RIPL': 1, 'RIPR': 1, 'RIS': 1, 'RMEL': 1, 'RMHL': 1, 'SAADL': 1, 'SABVR': 1, 'SIADR': 1, 'SIBDR': 1, 'SIBVL': 1, 'SIBVR': 1, 'SMDDL': 1, 'SMDDR': 1, 'VB2': 1}
Pearson r between dandiset dF/F0 and deposit trace: median 0.9718, 5th percentile 0.7582, min 0.2573
fraction with r >= 0.9: 0.817; r >= 0.95: 0.665; r >= 0.99: 0.191
regression slope deposit on dandiset: median 0.985 (IQR 0.948 to 0.998)
deposit trace length minus dandiset frames (after dropping the last frame): {0: 11, -1: 8, 1: 2} | leading NaN frames in the deposit: [np.int64(4)]

median r per file with the deposit shifted by -1, 0, +1 frames; best lag per file: {1: 21}
with deposit frame k + 1 against dandiset frame k: median r 0.9803, fraction with r >= 0.9: 0.858
             count  median     min  r_lag-1  r_lag+0  r_lag+1
file                                                         
20190924_01    177  0.9820  0.5925   0.9719   0.9820   0.9868
20190924_02    170  0.9855  0.6681   0.9764   0.9855   0.9890
20190924_03    174  0.9857  0.4538   0.9776   0.9857   0.9883
20190925_01    174  0.9736  0.4257   0.9606   0.9736   0.9807
20190925_04    171  0.9772  0.5673   0.9612   0.9772   0.9814
20190928_01    165  0.9696  0.5740   0.9462   0.9696   0.9772
20190928_03    163  0.9772  0.7462   0.9647   0.9772   0.9824
20190928_05    161  0.9659  0.4538   0.9471   0.9659   0.9757
20190928_07    173  0.9763  0.5164   0.9636   0.9763   0.9829
20190928_08    179  0.9518  0.2573   0.9279   0.9518   0.9794
20190928_11    173  0.9795  0.4005   0.9683   0.9795   0.9843
20190928_13    172  0.9723  0.4827   0.9548   0.9723   0.9798
20190929_02    146  0.9651  0.6987   0.9526   0.9651   0.9731
20190929_03    171  0.9636  0.5493   0.9497   0.9636   0.9718
20190929_05    163  0.9718  0.6229   0.9583   0.9718   0.9792
20190929_06    169  0.9695  0.6611   0.9570   0.9695   0.9772
20190929_07    173  0.9722  0.4762   0.9565   0.9722   0.9806
20191030_03    173  0.9703  0.4734   0.9470   0.9703   0.9797
20191030_07    165  0.9470  0.4938   0.9288   0.9470   0.9798
20191104_08    170  0.9612  0.5264   0.9233   0.9612   0.9706
20191104_10    168  0.9741  0.2803   0.9538   0.9741   0.9818

ss = []
for mode in MODES:
    d = data[mode]
    kept = [set(a["traces"]) for a in d["animals"]]
    by_name = collections.Counter(n for k in kept for n in k)
    by_cls = collections.Counter()
    for k in kept:
        by_cls.update({cls_of[n] for n in k})
    named_any = {n for r in recs for n in r["labels"] if n}
    ss.append({"duplicates": mode, "heads": len(d["animals"]),
               "names_any_head": len(named_any),
               "names_any_head_neurons_only": len(named_any - {"AMSOL", "AMSOR", "HMC"}),
               "names_after_removals": sum(v > npl.MIN_N for v in by_name.values()),
               "median_animals_per_name": float(np.median(list(by_name.values()))),
               "median_animals_per_name_N>3": float(np.median([v for v in by_name.values() if v > npl.MIN_N])),
               "median_animals_per_class": float(np.median(list(by_cls.values()))),
               "median_animals_per_class_N>3": float(np.median([v for v in by_cls.values() if v > npl.MIN_N]))})
ss = pd.DataFrame(ss)
print(ss.to_string())
dep_name = collections.Counter(n for f, n in DEP)
dep_cls = collections.Counter(c for f, c in {(f, npl.class_maps({m for _, m in DEP})[n]) for f, n in DEP})
print("the deposit itself: names", len(dep_name), "| median animals per name", np.median(list(dep_name.values())),
      "| classes", len(dep_cls), "| median animals per class", np.median(list(dep_cls.values())))
print("tail recordings in the dandiset:", sum("tail" in r["session_description"].lower() for r in recs))
print("expected: heads 21 (tol 0), head neurons 189 (170.1 to 207.9), tails 21, tail neurons 42, median 18 (16.2 to 19.8)")

  duplicates  heads  names_any_head  names_any_head_neurons_only  names_after_removals  median_animals_per_name  median_animals_per_name_N>3  median_animals_per_class  median_animals_per_class_N>3
0        one     21             193                          190                   184                     19.0                         19.0                      20.0                          20.0
1    deposit     21             193                          190                   184                     21.0                         21.0                      21.0                          21.0
2       both     21             193                          190                   184                     17.0                         17.0                      19.5                          20.0
3       none     21             193                          190                   185                     21.0                         21.0                      21.0                          21.0

the deposit itself: names 189 | median animals per name 21.0 | classes 109 | median animals per class 21.0
tail recordings in the dandiset: 0
expected: heads 21 (tol 0), head neurons 189 (170.1 to 207.9), tails 21, tail neurons 42, median 18 (16.2 to 19.8)

def head_checks(d, k_sd):
    out = []
    for a in d["animals"]:
        row = {"animal": a["id"]}
        for cls, stim, sign in (("ASE", "nacl", 1), ("AWC", "butanone", -1), ("AWC", "pentanedione", -1)):
            ok = False
            for n in (cls + "L", cls + "R"):
                if n not in a["traces"]:
                    continue
                tr = a["traces"][n]
                s0, s1 = a["stimuli"][stim]
                t = np.arange(len(tr)) / a["rate"]
                pre = tr[(t >= s0 - 10) & (t < s0)]
                diff = tr[(t >= s0) & (t < s1)].mean() - pre.mean()
                ok |= sign * diff > k_sd * pre.std()
            row[f"{cls}_{stim}"] = ok
        row["all_three"] = all(v for k, v in row.items() if k != "animal")
        out.append(row)
    return pd.DataFrame(out)

for label, mode, k in [(f"({'a' if k == 0 else 'b'}) {m}", m, k) for m in ("none", "one", "deposit", "both") for k in (0, 2)]:
    hc = head_checks(data[mode], k)
    print(label, "heads passing all three:", int(hc.all_three.sum()), "of", len(hc),
          "| per stimulus:", hc.drop(columns=["animal", "all_three"]).sum().to_dict())
hc_a = head_checks(data["none"], 0)
print(hc_a.to_string())

(a) none heads passing all three: 19 of 21 | per stimulus: {'ASE_nacl': 21, 'AWC_butanone': 20, 'AWC_pentanedione': 19}
(b) none heads passing all three: 16 of 21 | per stimulus: {'ASE_nacl': 20, 'AWC_butanone': 19, 'AWC_pentanedione': 18}
(a) one heads passing all three: 19 of 21 | per stimulus: {'ASE_nacl': 21, 'AWC_butanone': 20, 'AWC_pentanedione': 19}
(b) one heads passing all three: 16 of 21 | per stimulus: {'ASE_nacl': 20, 'AWC_butanone': 19, 'AWC_pentanedione': 18}
(a) deposit heads passing all three: 19 of 21 | per stimulus: {'ASE_nacl': 21, 'AWC_butanone': 20, 'AWC_pentanedione': 19}
(b) deposit heads passing all three: 16 of 21 | per stimulus: {'ASE_nacl': 20, 'AWC_butanone': 19, 'AWC_pentanedione': 18}
(a) both heads passing all three: 18 of 21 | per stimulus: {'ASE_nacl': 19, 'AWC_butanone': 20, 'AWC_pentanedione': 19}
(b) both heads passing all three: 15 of 21 | per stimulus: {'ASE_nacl': 18, 'AWC_butanone': 19, 'AWC_pentanedione': 18}
         animal  ASE_nacl  AWC_butanone  AWC_pentanedione  all_three
0   20190924_01      True          True              True       True
1   20190924_02      True          True              True       True
2   20190924_03      True          True              True       True
3   20190925_01      True          True              True       True
4   20190925_04      True          True              True       True
5   20190928_01      True          True              True       True
6   20190928_03      True          True              True       True
7   20190928_05      True          True              True       True
8   20190928_07      True          True              True       True
9   20190928_08      True          True              True       True
10  20190928_11      True          True              True       True
11  20190928_13      True          True              True       True
12  20190929_02      True          True              True       True
13  20190929_03      True          True              True       True
14  20190929_05      True          True              True       True
15  20190929_06      True          True              True       True
16  20190929_07      True          True              True       True
17  20191030_03      True          True              True       True
18  20191030_07      True          True             False      False
19  20191104_08      True          True              True       True
20  20191104_10      True         False             False      False

def test_table(d, norm="zscore"):
    rows = []
    for s in npl.STIMULI:
        for c, v in npl.class_tests(npl.responses(d, s, norm)).items():
            rows.append({"stimulus": s, "cls": c, **v})
    t = pd.DataFrame(rows)
    t["q_all"] = npl.storey_q(t.p.values)
    t["q_stim"] = np.nan
    for s in npl.STIMULI:
        m = t.stimulus == s
        t.loc[m, "q_stim"] = npl.storey_q(t.loc[m, "p"].values)
    t["q_bh_all"] = npl.bh_q(t.p.values)
    return t

TT = {mode: test_table(data[mode]) for mode in MODES}
TT["one, dF/F0"] = test_table(D, "raw")
TT["one, min-max"] = test_table(D, "minmax")
T = TT["one"]
print("tests per reading:", {m: len(t) for m, t in TT.items()}, "| classes tested:", T.cls.nunique())
print("non-neurons tested:", sorted(set(T.cls) & npl.NON_NEURONS))
print("Storey pi0, all tests:", round(npl.storey_pi0(T.p.values), 3),
      "| per stimulus:", {s: round(npl.storey_pi0(T[T.stimulus == s].p.values), 3) for s in npl.STIMULI})
for mode, t in TT.items():
    print(f"[{mode}] significant tests: q over all {int((t.q_all <= 0.05).sum())}, q per stimulus "
          f"{int((t.q_stim <= 0.05).sum())}, Benjamini-Hochberg over all {int((t.q_bh_all <= 0.05).sum())}, "
          f"p {int((t.p <= 0.05).sum())} of {len(t)}")
print("Storey pi0 over all tests, bootstrap seeds 0 to 9 (main reading):", [round(npl.storey_pi0(T.p.values, seed=k), 3) for k in range(10)])
READINGS = {"q_all": "(a) q over all tests", "q_stim": "(b) q per stimulus", "p": "(c) uncorrected p",
            "q_bh_all": "(a') Benjamini-Hochberg over all tests, a conservative check"}

def sig(t, cls, stim, key):
    r = t[(t.cls == cls) & (t.stimulus == stim)]
    if r.empty:
        return None
    return bool(r[key].iloc[0] <= 0.05), float(r.mean_diff.iloc[0])

def show(t, classes):
    sub = t[t.cls.isin(classes)].set_index(["cls", "stimulus"])[["n_traces", "n_animals", "mean_diff", "t", "p", "q_all", "q_stim"]]
    print(sub.round(5).to_string())

tests per reading: {'one': 282, 'deposit': 282, 'both': 282, 'none': 282, 'one, dF/F0': 282, 'one, min-max': 282} | classes tested: 94
non-neurons tested: []

Storey pi0, all tests: 0.222 | per stimulus: {'nacl': 0.236, 'butanone': 0.177, 'pentanedione': 0.411}
[one] significant tests: q over all 141, q per stimulus 128, Benjamini-Hochberg over all 64, p 106 of 282
[deposit] significant tests: q over all 106, q per stimulus 113, Benjamini-Hochberg over all 77, p 109 of 282
[both] significant tests: q over all 93, q per stimulus 114, Benjamini-Hochberg over all 61, p 105 of 282
[none] significant tests: q over all 103, q per stimulus 112, Benjamini-Hochberg over all 76, p 111 of 282
[one, dF/F0] significant tests: q over all 61, q per stimulus 72, Benjamini-Hochberg over all 50, p 82 of 282
[one, min-max] significant tests: q over all 111, q per stimulus 118, Benjamini-Hochberg over all 67, p 112 of 282

Storey pi0 over all tests, bootstrap seeds 0 to 9 (main reading): [0.222, 0.222, 0.222, 0.222, 0.222, 0.222, 0.222, 0.222, 0.222, 0.222]

S6 = npl.load_table_s6()
pcol = {"nacl": "NaCl p", "butanone": "2-butanone p", "pentanedione": "2,3-pentanedione p"}
cicol = {"nacl": "NaCl CI (ΔF/F0)", "butanone": "2-butanone CI (ΔF/F0)", "pentanedione": "2,3-pentanedione CI (ΔF/F0)"}
print("Table S6 classes:", len(S6), "| classes tested here and in Table S6:", len(set(T.cls) & set(S6.index)))
# The same tests on the deposit's own traces (deposit frame k + 1 is dandiset frame k, link-zenodo-traces;
# the deposit's last frame is NaN, the frame the STAR Methods remove).
stim_of = {a["id"]: a["stimuli"] for a in D["animals"]}
dep_data = {"animals": [{"id": f, "stimuli": stim_of[f], "rate": npl.RATE_HZ,
                         "traces": {n: z[1:-1] for (ff, n), z in DEP.items() if ff == f}} for f in stim_of],
            "class_of": npl.class_maps({n for _, n in DEP})}
print("deposit traces: NaN in the first 4 frames and the last frame, and nowhere else, in",
      sum(1 for z in DEP.values() if np.isfinite(z[4:-1]).all() and not np.isfinite(z[:4]).any() and not np.isfinite(z[-1])), "of", len(DEP))
TT_dep = {"deposit traces, dF/F0": test_table(dep_data, "raw"), "deposit traces, z-scored": test_table(dep_data)}
cmp_rows = []
for key, t in {**TT, **TT_dep}.items():
    row = {"reading": key}
    nt = t[t.stimulus == "nacl"].set_index("cls").n_traces
    com = [c for c in nt.index if c in S6.index]
    diff = nt[com] - S6.loc[com, "Samples"].astype(float)
    row.update(traces_equal=int((diff == 0).sum()), traces_minus_paper_median=float(np.median(diff)),
               traces_minus_paper_mean=round(float(diff.mean()), 2))
    for s in npl.STIMULI:
        ts = t[t.stimulus == s].set_index("cls")
        pp = pd.to_numeric(S6.loc[com, pcol[s]], errors="coerce")
        ok = pp.notna()
        lp, lq = np.log10(ts.loc[com, "p"][ok]), np.log10(pp[ok].astype(float))
        lohi = S6.loc[com, cicol[s]].astype(str).str.split(",", expand=True).apply(pd.to_numeric, errors="coerce")
        mid = lohi.mean(axis=1)[ok]
        row.update({f"{s[:4]}_rho": round(stats.spearmanr(lp, lq)[0], 3), f"{s[:4]}_dlog10p": round(float(np.median(np.abs(lp - lq))), 3),
                    f"{s[:4]}_sig_agree": int(((ts.loc[com, "p"][ok] <= 0.05) == (pp[ok] <= 0.05)).sum()),
                    f"{s[:4]}_effect_ratio": round(float(np.median(mid / ts.loc[com, "mean_diff"][ok])), 3)})
    row["n_classes"] = int(ok.sum())
    cmp_rows.append(row)
print(pd.DataFrame(cmp_rows).to_string(index=False))

# The q values of this attempt on the paper's own p values, against the paper's q values.
p6 = np.concatenate([pd.to_numeric(S6[pcol[s]], errors="coerce").dropna().values for s in npl.STIMULI])
q6 = np.concatenate([pd.to_numeric(S6[pcol[s].replace(" p", " q")], errors="coerce")[pd.to_numeric(S6[pcol[s]], errors="coerce").notna()].values for s in npl.STIMULI])
rank6 = np.empty(len(p6)); rank6[np.argsort(p6)] = np.arange(1, len(p6) + 1)
pos = p6 > 0
print("Table S6: %d stimulus tests, p <= 0.05 in %d, q <= 0.05 in %d; pi0 implied by its q values (q * rank / (m p)): median %.4f"
      % (len(p6), (p6 <= 0.05).sum(), (q6 <= 0.05).sum(), np.median((q6 * rank6 / (len(p6) * p6))[pos])))
print("storey_pi0 on the paper's p values, seeds 0 to 9:", [round(npl.storey_pi0(p6, seed=k), 3) for k in range(10)])
q6_ours = npl.storey_q(p6)
print("our q on the paper's p against the paper's q: q <= 0.05 agrees in %d of %d tests; largest |log10 ratio| %.3f"
      % (((q6_ours <= 0.05) == (q6 <= 0.05)).sum(), len(p6), np.max(np.abs(np.log10(q6_ours[pos] / q6[pos])))))

Table S6 classes: 109 | classes tested here and in Table S6: 94
deposit traces: NaN in the first 4 frames and the last frame, and nowhere else, in 3550 of 3550

                 reading  traces_equal  traces_minus_paper_median  traces_minus_paper_mean  nacl_rho  nacl_dlog10p  nacl_sig_agree  nacl_effect_ratio  buta_rho  buta_dlog10p  buta_sig_agree  buta_effect_ratio  pent_rho  pent_dlog10p  pent_sig_agree  pent_effect_ratio  n_classes
                     one            18                        0.0                     0.60     0.960         0.187              84              0.127     0.902         0.145              87              0.142     0.933         0.153              89              0.178         94
                 deposit            18                        2.0                     2.91     0.959         0.157              87              0.123     0.869         0.140              87              0.144     0.938         0.180              87              0.139         94
                    both            12                       -2.0                    -2.30     0.944         0.204              87              0.129     0.929         0.127              87              0.124     0.913         0.208              85              0.128         94
                    none            10                        3.0                     3.84     0.928         0.157              83              0.123     0.869         0.140              88              0.167     0.916         0.184              86              0.164         94
              one, dF/F0            18                        0.0                     0.60     0.728         0.400              74              0.271     0.637         0.382              80              0.448     0.621         0.385              80              0.421         94
            one, min-max            18                        0.0                     0.60     0.944         0.207              82              0.769     0.865         0.157              85              0.882     0.918         0.231              89              0.980         94
   deposit traces, dF/F0            18                        2.0                     2.91     0.748         0.361              72              0.329     0.569         0.395              80              0.445     0.654         0.448              78              0.389         94
deposit traces, z-scored            18                        2.0                     2.91     0.960         0.123              87              0.133     0.869         0.152              85              0.162     0.938         0.167              88              0.153         94
Table S6: 285 stimulus tests, p <= 0.05 in 106, q <= 0.05 in 102; pi0 implied by its q values (q * rank / (m p)): median 0.4236

storey_pi0 on the paper's p values, seeds 0 to 9: [0.478, 0.478, 0.421, 0.478, 0.478, 0.478, 0.421, 0.421, 0.478, 0.478]
our q on the paper's p against the paper's q: q <= 0.05 agrees in 273 of 285 tests; largest |log10 ratio| 0.060

def side_values(d, name, stim, norm="zscore"):
    out = []
    for a in d["animals"]:
        if name in a["traces"]:
            dur, pre = npl.window_means(npl.normalise(a["traces"][name], norm), *a["stimuli"][stim], a["rate"])
            out.append((a["id"], dur, pre))
    return out

for mode, norm in [(m, "zscore") for m in MODES] + [("one", "raw"), ("one", "minmax")]:
    d, t = data[mode], TT[mode if norm == "zscore" else ("one, dF/F0" if norm == "raw" else "one, min-max")]
    L = side_values(d, "ASEL", "nacl", norm)
    R = side_values(d, "ASER", "nacl", norm)
    mode = mode if norm == "zscore" else f"{mode}, {norm}"
    rd = np.array([x[1] for x in R]); rb = np.array([x[2] for x in R])
    pr = stats.ttest_rel(rd, rb).pvalue
    q_all = npl.storey_q(np.r_[t.p.values, pr])[-1]
    q_st = npl.storey_q(np.r_[t[t.stimulus == "nacl"].p.values, pr])[-1]
    ld = np.array([x[1] for x in L]); lb = np.array([x[2] for x in L])
    asym_resp = stats.ttest_ind(ld - lb, rd - rb)
    asym_mean = stats.ttest_ind(ld, rd)
    print(f"[{mode}] ASER NaCl: n={len(R)}, mean during-before={np.mean(rd - rb):+.3f}, "
          f"rising in {np.sum(rd > rb)} of {len(R)}, p={pr:.3g}, q(a)={q_all:.3g}, q(b)={q_st:.3g}")
    print(f"[{mode}] ASEL NaCl: n={len(L)}, mean during-before={np.mean(ld - lb):+.3f}, rising in {np.sum(ld > lb)} of {len(L)}")
    print(f"[{mode}] asymmetry (a) responses: t={asym_resp.statistic:.3f}, p={asym_resp.pvalue:.3g};"
          f" (b) means during: t={asym_mean.statistic:.3f}, p={asym_mean.pvalue:.3g}")

[one] ASER NaCl: n=21, mean during-before=+1.338, rising in 14 of 21, p=0.00282, q(a)=0.0036, q(b)=0.00313
[one] ASEL NaCl: n=19, mean during-before=+3.450, rising in 19 of 19
[one] asymmetry (a) responses: t=4.587, p=4.79e-05; (b) means during: t=4.955, p=1.53e-05

[deposit] ASER NaCl: n=21, mean during-before=+1.338, rising in 14 of 21, p=0.00282, q(a)=0.00709, q(b)=0.00575
[deposit] ASEL NaCl: n=21, mean during-before=+3.357, rising in 21 of 21
[deposit] asymmetry (a) responses: t=4.456, p=6.57e-05; (b) means during: t=4.475, p=6.2e-05

[both] ASER NaCl: n=15, mean during-before=+1.384, rising in 10 of 15, p=0.0151, q(a)=0.0272, q(b)=0.0124
[both] ASEL NaCl: n=17, mean during-before=+3.522, rising in 17 of 17
[both] asymmetry (a) responses: t=4.075, p=0.000311; (b) means during: t=4.845, p=3.61e-05

[none] ASER NaCl: n=21, mean during-before=+1.338, rising in 14 of 21, p=0.00282, q(a)=0.00753, q(b)=0.00578
[none] ASEL NaCl: n=21, mean during-before=+3.357, rising in 21 of 21
[none] asymmetry (a) responses: t=4.456, p=6.57e-05; (b) means during: t=4.475, p=6.2e-05

[one, raw] ASER NaCl: n=21, mean during-before=+0.594, rising in 14 of 21, p=0.00547, q(a)=0.0177, q(b)=0.00945
[one, raw] ASEL NaCl: n=19, mean during-before=+1.856, rising in 19 of 19
[one, raw] asymmetry (a) responses: t=4.565, p=5.11e-05; (b) means during: t=4.587, p=4.78e-05

[one, minmax] ASER NaCl: n=21, mean during-before=+0.208, rising in 14 of 21, p=0.00258, q(a)=0.00594, q(b)=0.00443
[one, minmax] ASEL NaCl: n=19, mean during-before=+0.553, rising in 19 of 19
[one, minmax] asymmetry (a) responses: t=5.014, p=1.27e-05; (b) means during: t=5.031, p=1.21e-05

SIX = ["AIN", "AVF", "AVH", "AVJ", "I2", "MI"]
show(T[T.stimulus == "nacl"], SIX)
for mode, t in TT.items():
    for key, lab in READINGS.items():
        res = {c: sig(t, c, "nacl", key) for c in SIX}
        print(f"[{mode}] {lab}: significant {sum(v[0] for v in res.values())} of 6:",
              [c for c, v in res.items() if v[0]])

              n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                    
AVH nacl            31         19    0.85212  5.32426  0.00001  0.00004  0.00004
AVJ nacl            35         21    0.38690  2.40719  0.02165  0.01574  0.01574
AIN nacl            41         21    0.41730  3.20669  0.00264  0.00344  0.00309
I2  nacl            39         20    0.39842  3.90890  0.00037  0.00068  0.00063
MI  nacl            17         17    0.29184  3.86399  0.00137  0.00192  0.00171
AVF nacl            32         19    0.24723  2.24774  0.03185  0.02074  0.01913
[one] (a) q over all tests: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[one] (b) q per stimulus: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[one] (c) uncorrected p: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[one] (a') Benjamini-Hochberg over all tests, a conservative check: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[deposit] (a) q over all tests: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[deposit] (b) q per stimulus: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[deposit] (c) uncorrected p: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[deposit] (a') Benjamini-Hochberg over all tests, a conservative check: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[both] (a) q over all tests: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[both] (b) q per stimulus: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[both] (c) uncorrected p: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[both] (a') Benjamini-Hochberg over all tests, a conservative check: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[none] (a) q over all tests: significant 5 of 6: ['AIN', 'AVH', 'AVJ', 'I2', 'MI']
[none] (b) q per stimulus: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[none] (c) uncorrected p: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[none] (a') Benjamini-Hochberg over all tests, a conservative check: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[one, dF/F0] (a) q over all tests: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[one, dF/F0] (b) q per stimulus: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[one, dF/F0] (c) uncorrected p: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']
[one, dF/F0] (a') Benjamini-Hochberg over all tests, a conservative check: significant 3 of 6: ['AVH', 'I2', 'MI']
[one, min-max] (a) q over all tests: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[one, min-max] (b) q per stimulus: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[one, min-max] (c) uncorrected p: significant 6 of 6: ['AIN', 'AVF', 'AVH', 'AVJ', 'I2', 'MI']
[one, min-max] (a') Benjamini-Hochberg over all tests, a conservative check: significant 4 of 6: ['AIN', 'AVH', 'I2', 'MI']

for mode, t in TT.items():
    rows = []
    for s in npl.STIMULI:
        ts = t[t.stimulus == s]
        row = {"stimulus": s, "classes": len(ts), "frac_negative_mean": round((ts.mean_diff < 0).mean(), 3)}
        for key in READINGS:
            row[f"down_{key}"] = int(((ts[key] <= 0.05) & (ts.mean_diff < 0)).sum())
            row[f"up_{key}"] = int(((ts[key] <= 0.05) & (ts.mean_diff > 0)).sum())
        rows.append(row)
    print(f"[{mode}]")
    print(pd.DataFrame(rows).to_string(index=False))

[one]
    stimulus  classes  frac_negative_mean  down_q_all  up_q_all  down_q_stim  up_q_stim  down_p  up_p  down_q_bh_all  up_q_bh_all
        nacl       94               0.234           6        48            6         49       3    41              2           22
    butanone       94               0.245           6        34            6         34       5    24              4           16
pentanedione       94               0.362          17        30           11         22      11    22              7           13
[deposit]
    stimulus  classes  frac_negative_mean  down_q_all  up_q_all  down_q_stim  up_q_stim  down_p  up_p  down_q_bh_all  up_q_bh_all
        nacl       94               0.234           2        36            3         36       3    36              2           22
    butanone       94               0.234           5        23            5         25       5    24              5           18
pentanedione       94               0.372          14        26           15         29      14    27             12           18
[both]
    stimulus  classes  frac_negative_mean  down_q_all  up_q_all  down_q_stim  up_q_stim  down_p  up_p  down_q_bh_all  up_q_bh_all
        nacl       94               0.223           2        36            5         49       4    39              2           21
    butanone       94               0.245           4        23            4         23       4    25              4           15
pentanedione       94               0.362           8        20           10         23      10    23              5           14
[none]
    stimulus  classes  frac_negative_mean  down_q_all  up_q_all  down_q_stim  up_q_stim  down_p  up_p  down_q_bh_all  up_q_bh_all
        nacl       94               0.223           2        36            3         38       3    38              2           23
    butanone       94               0.234           5        23            5         25       5    25              5           17
pentanedione       94               0.362          14        23           15         26      14    26             11           18
[one, dF/F0]
    stimulus  classes  frac_negative_mean  down_q_all  up_q_all  down_q_stim  up_q_stim  down_p  up_p  down_q_bh_all  up_q_bh_all
        nacl       94               0.266           3        15            4         26       4    26              3           13
    butanone       94               0.277           5        14            4         14       5    17              4           12
pentanedione       94               0.457          10        14           10         14      13    17              7           11
[one, min-max]
    stimulus  classes  frac_negative_mean  down_q_all  up_q_all  down_q_stim  up_q_stim  down_p  up_p  down_q_bh_all  up_q_bh_all
        nacl       94               0.234           5        39            5         40       5    39              2           21
    butanone       94               0.255           6        23            6         22       6    23              4           17
pentanedione       94               0.372          14        24           17         28      15    24             10           13

text = {r["identifier"]: [*r["labels"], *r["neuropal_labels"], r["subject_description"], r["session_description"],
                         r["stim_description"]] for r in recs}
import re
onoff = sorted({s for v in text.values() for s in v if s and re.search(r"(^|[^A-Z])(ON|OFF)([^A-Z]|$)|AWC_?(ON|OFF)", s.upper())})
print("AWC names in the files:", sorted({n for r in recs for n in r["labels"] if n.startswith("AWC")}),
      "| strings naming ON or OFF as a word, or AWCON/AWCOFF, among the trace labels, NeuroPAL labels and descriptions of all 21 files:", onoff)
is_off = dict(zip(zfiles, np.asarray(Z["is_AWCL_OFF"], float)))
print("deposit: AWCL is OFF in", int(np.nansum(list(is_off.values()))), "of", len(is_off), "heads")

# The NeuroPAL volume: srsx-3 colors AWC OFF (Results). Per head, each channel at AWCL over AWCR,
# divided by the same ratio for the panneuronal TagRFP-T.
col = []
for r in recs:
    c = r["awc_color"]
    row = {"file": r["identifier"], "deposit_AWCL_OFF": int(is_off[r["identifier"]])}
    for ch in c["AWCL"]:
        if ch != "Tag RFP-T":
            row[ch] = (c["AWCL"][ch] / c["AWCL"]["Tag RFP-T"]) / (c["AWCR"][ch] / c["AWCR"]["Tag RFP-T"])
    col.append(row)
col = pd.DataFrame(col)
print(col.round(2).to_string(index=False))
for ch in col.columns[2:]:
    agree = int(((col[ch] > 1) == (col.deposit_AWCL_OFF == 1)).sum())
    print(f"{ch}: AWCL brighter (relative to TagRFP-T) where the deposit has AWCL OFF in {agree} of 21 heads, "
          f"binomial p = {stats.binomtest(max(agree, 21 - agree), 21).pvalue:.2g}")
img_off = dict(zip(col.file, (col["CyOFP1"] > 1).astype(float)))

print("NeuroPAL fluorophore of srsx-3 in the paper's Table S1:", npl.reporter_fluorophore("srsx-3"))
for label, side in (("deposit label", is_off), ("NeuroPAL CyOFP1 label", img_off)):
    for mode, norm in [(m, "zscore") for m in MODES] + [("one", "raw")]:
        d = data[mode]
        for s in ("butanone", "pentanedione"):
            on, off = {"resp": [], "mean": []}, {"resp": [], "mean": []}
            for a in d["animals"]:
                for sd in "LR":
                    n = "AWC" + sd
                    if n not in a["traces"]:
                        continue
                    dur, pre = npl.window_means(npl.normalise(a["traces"][n], norm), *a["stimuli"][s], a["rate"])
                    tgt = off if (side[a["id"]] == 1) == (sd == "L") else on
                    tgt["resp"].append(dur - pre)
                    tgt["mean"].append(dur)
            for meas in ("resp", "mean"):
                r = stats.ttest_ind(on[meas], off[meas])
                print(f"[{label}, {mode}{'' if norm == 'zscore' else ', ' + norm}] {s}, {meas}: ON {np.mean(on[meas]):+.3f} (n={len(on[meas])}) vs OFF "
                      f"{np.mean(off[meas]):+.3f} (n={len(off[meas])}), t={r.statistic:.2f}, p={r.pvalue:.3g}")

AWC names in the files: ['AWCL', 'AWCR'] | strings naming ON or OFF as a word, or AWCON/AWCOFF, among the trace labels, NeuroPAL labels and descriptions of all 21 files: []
deposit: AWCL is OFF in 10 of 21 heads
       file  deposit_AWCL_OFF  mTagBFP2  CyOFP1  mNeptune 2.5  GFP-GCaMP
20190924_01                 1      1.30    3.75          1.02       1.50
20190924_02                 0      0.48    0.00          0.41       0.81
20190924_03                 1      1.05    4.72          1.24       2.31
20190925_01                 1      1.04    1.90          1.76       1.70
20190925_04                 0      0.83    0.34          1.68       0.63
20190928_01                 0      1.21    1.43          0.79       0.87
20190928_03                 0      0.61    0.27          0.66       0.40
20190928_05                 1      0.78    1.42          0.80       0.29
20190928_07                 0      0.63    0.63          0.44       0.54
20190928_08                 0      1.31    2.18          1.28       1.33
20190928_11                 1      1.05    1.87          1.30       1.37
20190928_13                 0      1.61    0.49          0.94       0.73
20190929_02                 1      1.20    7.99          1.54       3.67
20190929_03                 1      0.55    1.22          0.68       1.01
20190929_05                 0      1.79    0.21          1.26       1.99
20190929_06                 1      2.00    2.16          2.28       2.39
20190929_07                 1      0.55    1.04          0.62       0.75
20191030_03                 0      1.11    0.77          1.32       1.00
20191030_07                 1      0.89    1.02          1.00       1.31
20191104_08                 0      1.21    0.78          0.82       0.96
20191104_10                 0      0.98    0.55          0.90       0.82
mTagBFP2: AWCL brighter (relative to TagRFP-T) where the deposit has AWCL OFF in 11 of 21 heads, binomial p = 1
CyOFP1: AWCL brighter (relative to TagRFP-T) where the deposit has AWCL OFF in 19 of 21 heads, binomial p = 0.00022
mNeptune 2.5: AWCL brighter (relative to TagRFP-T) where the deposit has AWCL OFF in 14 of 21 heads, binomial p = 0.19
GFP-GCaMP: AWCL brighter (relative to TagRFP-T) where the deposit has AWCL OFF in 17 of 21 heads, binomial p = 0.0072

NeuroPAL fluorophore of srsx-3 in the paper's Table S1: CyOFP1
[deposit label, one] butanone, resp: ON -1.015 (n=21) vs OFF -1.168 (n=21), t=0.50, p=0.62
[deposit label, one] butanone, mean: ON -1.079 (n=21) vs OFF -0.989 (n=21), t=-0.29, p=0.771
[deposit label, one] pentanedione, resp: ON -1.044 (n=21) vs OFF -0.685 (n=21), t=-1.50, p=0.143
[deposit label, one] pentanedione, mean: ON -1.083 (n=21) vs OFF -0.800 (n=21), t=-1.18, p=0.246
[deposit label, deposit] butanone, resp: ON -1.015 (n=21) vs OFF -1.168 (n=21), t=0.50, p=0.62
[deposit label, deposit] butanone, mean: ON -1.079 (n=21) vs OFF -0.989 (n=21), t=-0.29, p=0.771
[deposit label, deposit] pentanedione, resp: ON -1.044 (n=21) vs OFF -0.685 (n=21), t=-1.50, p=0.143
[deposit label, deposit] pentanedione, mean: ON -1.083 (n=21) vs OFF -0.800 (n=21), t=-1.18, p=0.246

[deposit label, both] butanone, resp: ON -0.953 (n=19) vs OFF -1.139 (n=19), t=0.57, p=0.575
[deposit label, both] butanone, mean: ON -1.051 (n=19) vs OFF -0.932 (n=19), t=-0.37, p=0.717
[deposit label, both] pentanedione, resp: ON -1.091 (n=19) vs OFF -0.698 (n=19), t=-1.52, p=0.138
[deposit label, both] pentanedione, mean: ON -1.133 (n=19) vs OFF -0.772 (n=19), t=-1.40, p=0.171
[deposit label, none] butanone, resp: ON -1.015 (n=21) vs OFF -1.168 (n=21), t=0.50, p=0.62
[deposit label, none] butanone, mean: ON -1.079 (n=21) vs OFF -0.989 (n=21), t=-0.29, p=0.771
[deposit label, none] pentanedione, resp: ON -1.044 (n=21) vs OFF -0.685 (n=21), t=-1.50, p=0.143
[deposit label, none] pentanedione, mean: ON -1.083 (n=21) vs OFF -0.800 (n=21), t=-1.18, p=0.246
[deposit label, one, raw] butanone, resp: ON -0.648 (n=21) vs OFF -0.536 (n=21), t=-0.58, p=0.563
[deposit label, one, raw] butanone, mean: ON -0.632 (n=21) vs OFF -0.488 (n=21), t=-0.67, p=0.508
[deposit label, one, raw] pentanedione, resp: ON -0.599 (n=21) vs OFF -0.356 (n=21), t=-1.75, p=0.0873
[deposit label, one, raw] pentanedione, mean: ON -0.598 (n=21) vs OFF -0.382 (n=21), t=-1.70, p=0.0961
[NeuroPAL CyOFP1 label, one] butanone, resp: ON -1.155 (n=21) vs OFF -1.027 (n=21), t=-0.42, p=0.68
[NeuroPAL CyOFP1 label, one] butanone, mean: ON -1.100 (n=21) vs OFF -0.967 (n=21), t=-0.44, p=0.666

[NeuroPAL CyOFP1 label, one] pentanedione, resp: ON -1.092 (n=21) vs OFF -0.636 (n=21), t=-1.94, p=0.06
[NeuroPAL CyOFP1 label, one] pentanedione, mean: ON -1.109 (n=21) vs OFF -0.775 (n=21), t=-1.40, p=0.17
[NeuroPAL CyOFP1 label, deposit] butanone, resp: ON -1.155 (n=21) vs OFF -1.027 (n=21), t=-0.42, p=0.68
[NeuroPAL CyOFP1 label, deposit] butanone, mean: ON -1.100 (n=21) vs OFF -0.967 (n=21), t=-0.44, p=0.666
[NeuroPAL CyOFP1 label, deposit] pentanedione, resp: ON -1.092 (n=21) vs OFF -0.636 (n=21), t=-1.94, p=0.06
[NeuroPAL CyOFP1 label, deposit] pentanedione, mean: ON -1.109 (n=21) vs OFF -0.775 (n=21), t=-1.40, p=0.17
[NeuroPAL CyOFP1 label, both] butanone, resp: ON -1.108 (n=19) vs OFF -0.984 (n=19), t=-0.38, p=0.708
[NeuroPAL CyOFP1 label, both] butanone, mean: ON -1.075 (n=19) vs OFF -0.908 (n=19), t=-0.51, p=0.611
[NeuroPAL CyOFP1 label, both] pentanedione, resp: ON -1.145 (n=19) vs OFF -0.644 (n=19), t=-1.97, p=0.0563
[NeuroPAL CyOFP1 label, both] pentanedione, mean: ON -1.161 (n=19) vs OFF -0.743 (n=19), t=-1.63, p=0.112
[NeuroPAL CyOFP1 label, none] butanone, resp: ON -1.155 (n=21) vs OFF -1.027 (n=21), t=-0.42, p=0.68
[NeuroPAL CyOFP1 label, none] butanone, mean: ON -1.100 (n=21) vs OFF -0.967 (n=21), t=-0.44, p=0.666
[NeuroPAL CyOFP1 label, none] pentanedione, resp: ON -1.092 (n=21) vs OFF -0.636 (n=21), t=-1.94, p=0.06
[NeuroPAL CyOFP1 label, none] pentanedione, mean: ON -1.109 (n=21) vs OFF -0.775 (n=21), t=-1.40, p=0.17
[NeuroPAL CyOFP1 label, one, raw] butanone, resp: ON -0.694 (n=21) vs OFF -0.491 (n=21), t=-1.07, p=0.293
[NeuroPAL CyOFP1 label, one, raw] butanone, mean: ON -0.624 (n=21) vs OFF -0.496 (n=21), t=-0.59, p=0.555
[NeuroPAL CyOFP1 label, one, raw] pentanedione, resp: ON -0.596 (n=21) vs OFF -0.358 (n=21), t=-1.71, p=0.0941
[NeuroPAL CyOFP1 label, one, raw] pentanedione, mean: ON -0.586 (n=21) vs OFF -0.393 (n=21), t=-1.51, p=0.139

for mode, t in TT.items():
    print(f"[{mode}]")
    show(t, ["ASJ", "RIC"])
for mode, t in TT.items():
    for key, lab in READINGS.items():
        def part(cls, stim, want):
            s, m = sig(t, cls, stim, key)
            if want == "up":
                return s and m > 0
            if want == "down":
                return s and m < 0
            return not (s and m > 0)
        parts = {"ASJ+NaCl": part("ASJ", "nacl", "up"), "ASJ-but": part("ASJ", "butanone", "down"),
                 "ASJ-pent": part("ASJ", "pentanedione", "down"), "RIC+pent": part("RIC", "pentanedione", "up"),
                 "RIC not+NaCl": part("RIC", "nacl", "not_up"), "RIC not+but": part("RIC", "butanone", "not_up")}
        print(f"[{mode}] {lab}: {sum(parts.values())} of 6 hold:", parts)

[one]
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
ASJ nacl                37         21    0.79770  4.13358  0.00020  0.00042  0.00043
RIC nacl                40         21    0.16765  2.20836  0.03317  0.02137  0.01930
ASJ butanone            37         21   -0.51942 -4.18396  0.00018  0.00038  0.00033
RIC butanone            40         21   -0.02232 -0.25536  0.79979  0.19262  0.15006
ASJ pentanedione        37         21   -0.24859 -1.78883  0.08206  0.04007  0.07555
RIC pentanedione        40         21    0.28813  2.96973  0.00508  0.00557  0.01034
[deposit]
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
ASJ nacl                37         21    0.78267  3.99527  0.00031  0.00112  0.00081
RIC nacl                42         21    0.14539  1.96412  0.05632  0.06534  0.05321
ASJ butanone            37         21   -0.56025 -4.54764  0.00006  0.00034  0.00019
RIC butanone            42         21    0.00025  0.00294  0.99767  0.46791  0.26534
ASJ pentanedione        37         21   -0.37866 -2.72692  0.00982  0.01881  0.01201
RIC pentanedione        42         21    0.32572  3.39271  0.00154  0.00445  0.00307
[both]
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
ASJ nacl                33         20    0.82111  3.93333  0.00042  0.00172  0.00094
RIC nacl                37         21    0.16770  2.24313  0.03114  0.04415  0.01821
ASJ butanone            33         20   -0.56218 -4.24290  0.00018  0.00086  0.00073
RIC butanone            37         21   -0.03129 -0.33705  0.73804  0.38466  0.33456
ASJ pentanedione        33         20   -0.26510 -1.75423  0.08896  0.08840  0.08155
RIC pentanedione        37         21    0.31931  3.24987  0.00251  0.00703  0.00674
[none]
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
ASJ nacl                40         21    0.73535  4.01898  0.00026  0.00114  0.00082
RIC nacl                42         21    0.14539  1.96412  0.05632  0.06965  0.05368
ASJ butanone            40         21   -0.52642 -4.51461  0.00006  0.00037  0.00020
RIC butanone            42         21    0.00025  0.00294  0.99767  0.49878  0.26534
ASJ pentanedione        40         21   -0.31768 -2.35636  0.02357  0.03693  0.02510
RIC pentanedione        42         21    0.32572  3.39271  0.00154  0.00475  0.00339
[one, dF/F0]
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
ASJ nacl                37         21    0.04841  0.22056  0.82668  0.49336  0.26539
RIC nacl                40         21   -0.00747 -0.23295  0.81702  0.49094  0.26527
ASJ butanone            37         21   -0.17191 -3.96714  0.00033  0.00260  0.00383
RIC butanone            40         21    0.01731  0.74737  0.45933  0.34393  0.43252
ASJ pentanedione        37         21   -0.05524 -0.85163  0.40005  0.31880  0.40972
RIC pentanedione        40         21    0.03060  1.01310  0.31726  0.27328  0.33083
[one, min-max]
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
ASJ nacl                37         21    0.12054  3.95560  0.00034  0.00128  0.00098
RIC nacl                40         21    0.02331  1.87023  0.06897  0.06981  0.05162
ASJ butanone            37         21   -0.08713 -4.40902  0.00009  0.00048  0.00044
RIC butanone            40         21   -0.00172 -0.12177  0.90371  0.38478  0.32275
ASJ pentanedione        37         21   -0.04887 -2.19761  0.03449  0.04190  0.02112
RIC pentanedione        40         21    0.04286  2.82293  0.00745  0.01450  0.00808
[one] (a) q over all tests: 5 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': False, 'RIC not+but': True}
[one] (b) q per stimulus: 4 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': False, 'RIC not+but': True}
[one] (c) uncorrected p: 4 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': False, 'RIC not+but': True}
[one] (a') Benjamini-Hochberg over all tests, a conservative check: 5 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[deposit] (a) q over all tests: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[deposit] (b) q per stimulus: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[deposit] (c) uncorrected p: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[deposit] (a') Benjamini-Hochberg over all tests, a conservative check: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[both] (a) q over all tests: 4 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': False, 'RIC not+but': True}
[both] (b) q per stimulus: 4 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': False, 'RIC not+but': True}
[both] (c) uncorrected p: 4 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': False, 'RIC not+but': True}
[both] (a') Benjamini-Hochberg over all tests, a conservative check: 5 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[none] (a) q over all tests: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[none] (b) q per stimulus: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[none] (c) uncorrected p: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[none] (a') Benjamini-Hochberg over all tests, a conservative check: 5 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, dF/F0] (a) q over all tests: 3 of 6 hold: {'ASJ+NaCl': False, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': False, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, dF/F0] (b) q per stimulus: 3 of 6 hold: {'ASJ+NaCl': False, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': False, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, dF/F0] (c) uncorrected p: 3 of 6 hold: {'ASJ+NaCl': False, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': False, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, dF/F0] (a') Benjamini-Hochberg over all tests, a conservative check: 3 of 6 hold: {'ASJ+NaCl': False, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': False, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, min-max] (a) q over all tests: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, min-max] (b) q per stimulus: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, min-max] (c) uncorrected p: 6 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': True, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}
[one, min-max] (a') Benjamini-Hochberg over all tests, a conservative check: 5 of 6 hold: {'ASJ+NaCl': True, 'ASJ-but': True, 'ASJ-pent': False, 'RIC+pent': True, 'RIC not+NaCl': True, 'RIC not+but': True}

FIVE = ["AIN", "AVH", "RIF", "RIG", "RIR"]
show(T, FIVE)
for mode, t in TT.items():
    for key, lab in READINGS.items():
        res = {c: [s[:4] for s in npl.STIMULI if sig(t, c, s, key) and sig(t, c, s, key)[0]] for c in FIVE}
        print(f"[{mode}] {lab}: {sum(bool(v) for v in res.values())} of 5 respond:", res)

                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
AVH nacl                31         19    0.85212  5.32426  0.00001  0.00004  0.00004
AIN nacl                41         21    0.41730  3.20669  0.00264  0.00344  0.00309
RIR nacl                17         17    0.21230  0.94101  0.36069  0.11865  0.11963
RIG nacl                11          9    0.14383  0.58616  0.57076  0.15925  0.16662
RIF nacl                20         11   -0.17708 -0.89512  0.38192  0.12354  0.12300
AVH butanone            31         19    0.40669  2.60357  0.01421  0.01206  0.01076
AIN butanone            41         21    0.21074  2.12263  0.04003  0.02429  0.02382
RIR butanone            17         17    0.86443  3.28451  0.00467  0.00540  0.00474
RIG butanone            11          9   -0.09540 -0.34168  0.73966  0.18497  0.14703
RIF butanone            20         11   -0.07642 -0.71213  0.48504  0.14716  0.12631
AVH pentanedione        31         19    0.92494  4.82456  0.00004  0.00014  0.00018
AIN pentanedione        41         21    0.14416  1.48531  0.14530  0.06054  0.10710
RIR pentanedione        17         17    1.68121  6.71619  0.00000  0.00002  0.00003
RIG pentanedione        11          9   -0.54415 -3.03571  0.01255  0.01137  0.01936
RIF pentanedione        20         11   -0.34784 -3.13495  0.00545  0.00588  0.01054
[one] (a) q over all tests: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[one] (b) q per stimulus: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[one] (c) uncorrected p: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[one] (a') Benjamini-Hochberg over all tests, a conservative check: 4 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'pent'], 'RIF': ['pent'], 'RIG': [], 'RIR': ['buta', 'pent']}

[deposit] (a) q over all tests: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[deposit] (b) q per stimulus: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[deposit] (c) uncorrected p: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[deposit] (a') Benjamini-Hochberg over all tests, a conservative check: 4 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': [], 'RIR': ['buta', 'pent']}
[both] (a) q over all tests: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}

[both] (b) q per stimulus: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[both] (c) uncorrected p: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[both] (a') Benjamini-Hochberg over all tests, a conservative check: 4 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'pent'], 'RIF': ['pent'], 'RIG': [], 'RIR': ['buta', 'pent']}
[none] (a) q over all tests: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[none] (b) q per stimulus: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[none] (c) uncorrected p: 5 of 5 respond: {'AIN': ['nacl', 'buta'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[none] (a') Benjamini-Hochberg over all tests, a conservative check: 4 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': [], 'RIR': ['buta', 'pent']}
[one, dF/F0] (a) q over all tests: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta']}
[one, dF/F0] (b) q per stimulus: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta']}
[one, dF/F0] (c) uncorrected p: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta']}

[one, dF/F0] (a') Benjamini-Hochberg over all tests, a conservative check: 3 of 5 respond: {'AIN': [], 'AVH': ['nacl', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': []}
[one, min-max] (a) q over all tests: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[one, min-max] (b) q per stimulus: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[one, min-max] (c) uncorrected p: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}
[one, min-max] (a') Benjamini-Hochberg over all tests, a conservative check: 5 of 5 respond: {'AIN': ['nacl'], 'AVH': ['nacl', 'buta', 'pent'], 'RIF': ['pent'], 'RIG': ['pent'], 'RIR': ['buta', 'pent']}

tail_names = [n for r in recs for n in r["labels"] if n[:3] in ("LUA", "PVN", "PVQ", "PVR", "PVW")]
print("tail neurons named in the dandiset:", tail_names if tail_names else "none",
      "| session descriptions:", sorted({r["session_description"] for r in recs}))

tail neurons named in the dandiset: none | session descriptions: ['C. elegans head NeuroPAL and Calcium imaging']

PHAR = ["I1", "I2", "I3", "I4", "I5", "I6", "M1", "M2", "M3", "M4", "M5", "MC", "MI", "NSM"]
print("pharyngeal classes tested:", sorted(set(T.cls) & set(PHAR)), "| not tested:", sorted(set(PHAR) - set(T.cls)))
show(T, PHAR)
for mode, t in TT.items():
    for key, lab in READINGS.items():
        res = {s: [c for c in PHAR if sig(t, c, s, key) and sig(t, c, s, key)[0]] for s in npl.STIMULI}
        print(f"[{mode}] {lab}:", {s: (len(v), v) for s, v in res.items()})

pharyngeal classes tested: ['I1', 'I2', 'I3', 'I4', 'I5', 'I6', 'M1', 'M2', 'M3', 'M4', 'M5', 'MC', 'MI', 'NSM'] | not tested: []
                  n_traces  n_animals  mean_diff        t        p    q_all   q_stim
cls stimulus                                                                        
M1  nacl                20         20    0.10669  1.96166  0.06462  0.03484  0.02992
I5  nacl                17         17    0.06933  0.47295  0.64264  0.16767  0.17702
M5  nacl                21         21    0.10821  1.28716  0.21274  0.07812  0.07786
M2  nacl                42         21    0.11934  2.24518  0.03021  0.01988  0.01865
I4  nacl                20         20   -0.05871 -0.62598  0.53878  0.15662  0.16180
I6  nacl                21         21    0.01952  0.11575  0.90900  0.20659  0.21957
I2  nacl                39         20    0.39842  3.90890  0.00037  0.00068  0.00063
NSM nacl                42         21    0.06078  1.04341  0.30287  0.10288  0.10355
MC  nacl                35         21    0.25755  4.58954  0.00006  0.00018  0.00022
M3  nacl                41         21    0.16858  2.14084  0.03843  0.02378  0.02082
M4  nacl                21         21    0.16300  2.58336  0.01775  0.01356  0.01364
MI  nacl                17         17    0.29184  3.86399  0.00137  0.00192  0.00171
I3  nacl                17         17    0.68431  3.00213  0.00844  0.00879  0.00816
I1  nacl                18         10    0.00320  0.02293  0.98198  0.21850  0.23215
M1  butanone            20         20    0.10684  2.48328  0.02252  0.01592  0.01501
I5  butanone            17         17    0.05927  0.32435  0.74987  0.18525  0.14703
M5  butanone            21         21    0.18559  1.95533  0.06466  0.03484  0.03476
M2  butanone            42         21    0.12626  2.61147  0.01254  0.01137  0.00995
I4  butanone            20         20    0.20140  1.75819  0.09481  0.04294  0.03951
I6  butanone            21         21    0.34721  2.88889  0.00908  0.00920  0.00800
I2  butanone            39         20    0.23049  1.85844  0.07087  0.03691  0.03579
NSM butanone            42         21    0.28666  3.74039  0.00056  0.00093  0.00078
MC  butanone            35         21    0.22767  3.03489  0.00459  0.00540  0.00474
M3  butanone            41         21    0.22630  4.13878  0.00017  0.00038  0.00033
M4  butanone            21         21    0.13871  1.53189  0.14122  0.05923  0.05604
MI  butanone            17         17    0.13888  2.17342  0.04511  0.02660  0.02592
I3  butanone            17         17    0.18742  2.64726  0.01757  0.01356  0.01220
I1  butanone            18         10    0.05963  0.60344  0.55418  0.15793  0.12667
M1  pentanedione        20         20    0.11161  1.77085  0.09263  0.04243  0.08140
I5  pentanedione        17         17    0.03608  0.38653  0.70419  0.17819  0.32415
M5  pentanedione        21         21    0.00143  0.01858  0.98536  0.21850  0.40553
M2  pentanedione        42         21    0.05991  0.70905  0.48231  0.14705  0.26267
I4  pentanedione        20         20   -0.02909 -0.32821  0.74634  0.18525  0.33951
I6  pentanedione        21         21    0.06399  0.64090  0.52887  0.15578  0.28013
I2  pentanedione        39         20    0.32530  2.64399  0.01184  0.01121  0.01936
NSM pentanedione        42         21    0.08518  1.18213  0.24396  0.08609  0.14513
MC  pentanedione        35         21    0.22014  3.30680  0.00223  0.00297  0.00576
M3  pentanedione        41         21    0.02746  0.57301  0.56985  0.15925  0.29379
M4  pentanedione        21         21    0.13537  1.47640  0.15541  0.06348  0.10926
MI  pentanedione        17         17    0.11460  2.58436  0.01996  0.01471  0.02417
I3  pentanedione        17         17    0.06050  1.62312  0.12410  0.05380  0.09793
I1  pentanedione        18         10    0.05738  1.32107  0.20399  0.07680  0.12722

[one] (a) q over all tests: {'nacl': (8, ['I2', 'I3', 'M1', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (11, ['I2', 'I3', 'I4', 'I6', 'M1', 'M2', 'M3', 'M5', 'MC', 'MI', 'NSM']), 'pentanedione': (4, ['I2', 'M1', 'MC', 'MI'])}
[one] (b) q per stimulus: {'nacl': (8, ['I2', 'I3', 'M1', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (11, ['I2', 'I3', 'I4', 'I6', 'M1', 'M2', 'M3', 'M5', 'MC', 'MI', 'NSM']), 'pentanedione': (3, ['I2', 'MC', 'MI'])}
[one] (c) uncorrected p: {'nacl': (7, ['I2', 'I3', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (3, ['I2', 'MC', 'MI'])}

[one] (a') Benjamini-Hochberg over all tests, a conservative check: {'nacl': (4, ['I2', 'I3', 'MC', 'MI']), 'butanone': (4, ['I6', 'M3', 'MC', 'NSM']), 'pentanedione': (1, ['MC'])}

[deposit] (a) q over all tests: {'nacl': (7, ['I2', 'I3', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (4, ['I2', 'I3', 'MC', 'MI'])}

[deposit] (b) q per stimulus:

 {'nacl': (7, ['I2', 'I3', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (4, ['I2', 'I3', 'MC', 'MI'])}
[deposit] (c) uncorrected p: {'nacl': (7, ['I2', 'I3', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (4, ['I2', 'I3', 'MC', 'MI'])}
[deposit] (a') Benjamini-Hochberg over all tests, a conservative check: {'nacl': (4, ['I2', 'I3', 'MC', 'MI']), 'butanone': (6, ['I3', 'I6', 'M2', 'M3', 'MC', 'NSM']), 'pentanedione': (1, ['MC'])}
[both] (a) q over all tests: {'nacl': (6, ['I2', 'I3', 'M2', 'M4', 'MC', 'MI']), 'butanone': (7, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'NSM']), 'pentanedione': (2, ['I2', 'MC'])}

[both] (b) q per stimulus: {'nacl': (8, ['I2', 'I3', 'M1', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (7, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'NSM']), 'pentanedione': (3, ['I2', 'I3', 'MC'])}

[both] (c) uncorrected p: {'nacl': (7, ['I2', 'I3', 'M2', 'M3', 'M4', 'MC', 'MI']), 'butanone': (7, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'NSM']), 'pentanedione': (3, ['I2', 'I3', 'MC'])}

[both] (a') Benjamini-Hochberg over all tests, a conservative check: {'nacl': (3, ['I2', 'MC', 'MI']), 'butanone': (4, ['I6', 'M3', 'MC', 'NSM']), 'pentanedione': (1, ['MC'])}
[none] (a) q over all tests: {'nacl': (6, ['I2', 'I3', 'M2', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (3, ['I2', 'I3', 'MC'])}
[none] (b) q per stimulus: {'nacl': (6, ['I2', 'I3', 'M2', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (4, ['I2', 'I3', 'MC', 'MI'])}
[none] (c) uncorrected p: {'nacl': (6, ['I2', 'I3', 'M2', 'M4', 'MC', 'MI']), 'butanone': (8, ['I3', 'I6', 'M1', 'M2', 'M3', 'MC', 'MI', 'NSM']), 'pentanedione': (4, ['I2', 'I3', 'MC', 'MI'])}

[none] (a') Benjamini-Hochberg over all tests, a conservative check: {'nacl': (4, ['I2', 'I3', 'MC', 'MI']), 'butanone': (6, ['I3', 'I6', 'M2', 'M3', 'MC', 'NSM']), 'pentanedione': (1, ['MC'])}

for mode, t in TT.items():
    for key, lab in READINGS.items():
        sets = {s: set(t[(t.stimulus == s) & (t[key] <= 0.05)].cls) for s in npl.STIMULI}
        pairs = {f"{a[:4]}&{b[:4]}": (len(sets[a] & sets[b]), sets[a] == sets[b])
                 for i, a in enumerate(npl.STIMULI) for b in npl.STIMULI[i + 1:]}
        print(f"[{mode}] {lab}: sizes", {s[:4]: len(v) for s, v in sets.items()},
              "| shared, identical:", pairs, "| in all three:", len(set.intersection(*sets.values())))

[one] (a) q over all tests: sizes {'nacl': 54, 'buta': 40, 'pent': 47} | shared, identical: {'nacl&buta': (25, False), 'nacl&pent': (28, False), 'buta&pent': (26, False)} | in all three: 18
[one] (b) q per stimulus: sizes {'nacl': 55, 'buta': 40, 'pent': 33} | shared, identical: {'nacl&buta': (26, False), 'nacl&pent': (21, False), 'buta&pent': (19, False)} | in all three: 14
[one] (c) uncorrected p: sizes {'nacl': 44, 'buta': 29, 'pent': 33} | shared, identical: {'nacl&buta': (16, False), 'nacl&pent': (17, False), 'buta&pent': (17, False)} | in all three: 9
[one] (a') Benjamini-Hochberg over all tests, a conservative check: sizes {'nacl': 24, 'buta': 20, 'pent': 20} | shared, identical: {'nacl&buta': (8, False), 'nacl&pent': (6, False), 'buta&pent': (12, False)} | in all three: 5
[deposit] (a) q over all tests: sizes {'nacl': 38, 'buta': 28, 'pent': 40} | shared, identical: {'nacl&buta': (15, False), 'nacl&pent': (20, False), 'buta&pent': (19, False)} | in all three: 11
[deposit] (b) q per stimulus: sizes {'nacl': 39, 'buta': 30, 'pent': 44} | shared, identical: {'nacl&buta': (17, False), 'nacl&pent': (22, False), 'buta&pent': (20, False)} | in all three: 13
[deposit] (c) uncorrected p: sizes {'nacl': 39, 'buta': 29, 'pent': 41} | shared, identical: {'nacl&buta': (16, False), 'nacl&pent': (21, False), 'buta&pent': (19, False)} | in all three: 12
[deposit] (a') Benjamini-Hochberg over all tests, a conservative check: sizes {'nacl': 24, 'buta': 23, 'pent': 30} | shared, identical: {'nacl&buta': (10, False), 'nacl&pent': (10, False), 'buta&pent': (16, False)} | in all three: 8
[both] (a) q over all tests: sizes {'nacl': 38, 'buta': 27, 'pent': 28} | shared, identical: {'nacl&buta': (14, False), 'nacl&pent': (12, False), 'buta&pent': (14, False)} | in all three: 7
[both] (b) q per stimulus: sizes {'nacl': 54, 'buta': 27, 'pent': 33} | shared, identical: {'nacl&buta': (20, False), 'nacl&pent': (22, False), 'buta&pent': (16, False)} | in all three: 12
[both] (c) uncorrected p: sizes {'nacl': 43, 'buta': 29, 'pent': 33} | shared, identical: {'nacl&buta': (16, False), 'nacl&pent': (16, False), 'buta&pent': (16, False)} | in all three: 10
[both] (a') Benjamini-Hochberg over all tests, a conservative check: sizes {'nacl': 23, 'buta': 19, 'pent': 19} | shared, identical: {'nacl&buta': (7, False), 'nacl&pent': (6, False), 'buta&pent': (13, False)} | in all three: 5
[none] (a) q over all tests: sizes {'nacl': 38, 'buta': 28, 'pent': 37} | shared, identical: {'nacl&buta': (14, False), 'nacl&pent': (17, False), 'buta&pent': (18, False)} | in all three: 10
[none] (b) q per stimulus: sizes {'nacl': 41, 'buta': 30, 'pent': 41} | shared, identical: {'nacl&buta': (17, False), 'nacl&pent': (21, False), 'buta&pent': (20, False)} | in all three: 13
[none] (c) uncorrected p: sizes {'nacl': 41, 'buta': 30, 'pent': 40} | shared, identical: {'nacl&buta': (17, False), 'nacl&pent': (21, False), 'buta&pent': (20, False)} | in all three: 13
[none] (a') Benjamini-Hochberg over all tests, a conservative check: sizes {'nacl': 25, 'buta': 22, 'pent': 29} | shared, identical: {'nacl&buta': (10, False), 'nacl&pent': (9, False), 'buta&pent': (14, False)} | in all three: 7
[one, dF/F0] (a) q over all tests: sizes {'nacl': 18, 'buta': 19, 'pent': 24} | shared, identical: {'nacl&buta': (9, False), 'nacl&pent': (10, False), 'buta&pent': (12, False)} | in all three: 7
[one, dF/F0] (b) q per stimulus: sizes {'nacl': 30, 'buta': 18, 'pent': 24} | shared, identical: {'nacl&buta': (10, False), 'nacl&pent': (10, False), 'buta&pent': (11, False)} | in all three: 7
[one, dF/F0] (c) uncorrected p: sizes {'nacl': 30, 'buta': 22, 'pent': 30} | shared, identical: {'nacl&buta': (11, False), 'nacl&pent': (11, False), 'buta&pent': (14, False)} | in all three: 8
[one, dF/F0] (a') Benjamini-Hochberg over all tests, a conservative check: sizes {'nacl': 16, 'buta': 16, 'pent': 18} | shared, identical: {'nacl&buta': (7, False), 'nacl&pent': (8, False), 'buta&pent': (9, False)} | in all three: 5
[one, min-max] (a) q over all tests: sizes {'nacl': 44, 'buta': 29, 'pent': 38} | shared, identical: {'nacl&buta': (15, False), 'nacl&pent': (19, False), 'buta&pent': (19, False)} | in all three: 11
[one, min-max] (b) q per stimulus: sizes {'nacl': 45, 'buta': 28, 'pent': 45} | shared, identical: {'nacl&buta': (16, False), 'nacl&pent': (21, False), 'buta&pent': (22, False)} | in all three: 12
[one, min-max] (c) uncorrected p: sizes {'nacl': 44, 'buta': 29, 'pent': 39} | shared, identical: {'nacl&buta': (15, False), 'nacl&pent': (19, False), 'buta&pent': (19, False)} | in all three: 11
[one, min-max] (a') Benjamini-Hochberg over all tests, a conservative check: sizes {'nacl': 23, 'buta': 21, 'pent': 23} | shared, identical: {'nacl&buta': (10, False), 'nacl&pent': (7, False), 'buta&pent': (12, False)} | in all three: 6

def avf_rows(d):
    rows = []
    for a in d["animals"]:
        for n in ("AVFL", "AVFR"):
            if n not in a["traces"]:
                continue
            x = a["traces"][n][4:]
            y = npl.medfilt_matlab(x, int(round(0.5 * a["rate"])))
            dur = len(y) / a["rate"]
            f, pxx = signal.welch(y - y.mean(), fs=a["rate"], nperseg=256)
            fpk = f[1:][np.argmax(pxx[1:])]
            row = {"animal": a["id"], "neuron": n, "duration_s": dur, "spectral_peak_Hz": fpk,
                   "spectral_ok": 0.1 <= fpk <= 1.0}
            for k in (1.0, 0.5):
                pk, _ = signal.find_peaks(y, prominence=k * y.std())
                row[f"freq_{k}sd"] = len(pk) / dur
            rows.append(row)
    return pd.DataFrame(rows)

for mode in MODES:
    av = avf_rows(data[mode])
    print(f"[{mode}] AVF traces: {len(av)} from {av.animal.nunique()} heads; spectral peak in 0.1-1 Hz: {int(av.spectral_ok.sum())}")
    for k in ("1.0", "0.5"):
        col = f"freq_{k}sd"
        print(f"[{mode}] peaks at {k} SD: median over all {av[col].median():.3f} Hz "
              f"(range {av[col].min():.3f} to {av[col].max():.3f});"
              f" over spectral-peak traces {av[av.spectral_ok][col].median():.3f} Hz (n={int(av.spectral_ok.sum())})")
    print(f"[{mode}] median spectral peak {av.spectral_peak_Hz.median():.3f} Hz")
av = avf_rows(D)
print(av.round(3).to_string())
fig, axs = plt.subplots(5, 1, figsize=(9, 8), sharex=True)
for ax, a in zip(axs, [a for a in D["animals"] if "AVFL" in a["traces"]][:5]):
    for n, col in (("AVFL", "tab:blue"), ("AVFR", "tab:red")):
        if n in a["traces"]:
            x = a["traces"][n][4:]
            ax.plot(np.arange(len(x)) / 4 + 1, npl.medfilt_matlab(x, 2), color=col, lw=0.7)
    ax.set_ylabel(a["id"][-7:], fontsize=7)
axs[-1].set_xlabel("time (s)")
plt.tight_layout()
plt.show()

[one] AVF traces: 32 from 19 heads; spectral peak in 0.1-1 Hz: 3
[one] peaks at 1.0 SD: median over all 0.121 Hz (range 0.046 to 0.268); over spectral-peak traces 0.172 Hz (n=3)
[one] peaks at 0.5 SD: median over all 0.199 Hz (range 0.104 to 0.301); over spectral-peak traces 0.180 Hz (n=3)
[one] median spectral peak 0.016 Hz

[deposit] AVF traces: 31 from 18 heads; spectral peak in 0.1-1 Hz: 3
[deposit] peaks at 1.0 SD: median over all 0.125 Hz (range 0.046 to 0.268); over spectral-peak traces 0.172 Hz (n=3)
[deposit] peaks at 0.5 SD: median over all 0.205 Hz (range 0.109 to 0.301); over spectral-peak traces 0.180 Hz (n=3)
[deposit] median spectral peak 0.016 Hz

[both] AVF traces: 29 from 18 heads; spectral peak in 0.1-1 Hz: 3
[both] peaks at 1.0 SD: median over all 0.121 Hz (range 0.046 to 0.268); over spectral-peak traces 0.172 Hz (n=3)
[both] peaks at 0.5 SD: median over all 0.192 Hz (range 0.104 to 0.301); over spectral-peak traces 0.180 Hz (n=3)
[both] median spectral peak 0.016 Hz

[none] AVF traces: 33 from 19 heads; spectral peak in 0.1-1 Hz: 3
[none] peaks at 1.0 SD: median over all 0.121 Hz (range 0.046 to 0.268); over spectral-peak traces 0.172 Hz (n=3)
[none] peaks at 0.5 SD: median over all 0.201 Hz (range 0.104 to 0.301); over spectral-peak traces 0.180 Hz (n=3)
[none] median spectral peak 0.016 Hz

         animal neuron  duration_s  spectral_peak_Hz  spectral_ok  freq_1.0sd  freq_0.5sd
0   20190924_01   AVFL      239.00             0.188         True       0.172       0.180
1   20190924_01   AVFR      239.00             0.016        False       0.138       0.197
2   20190924_02   AVFL      239.00             0.016        False       0.176       0.226
3   20190924_02   AVFR      239.00             0.031        False       0.146       0.201
4   20190924_03   AVFL      239.25             0.016        False       0.205       0.247
5   20190924_03   AVFR      239.25             0.016        False       0.067       0.104
6   20190925_01   AVFL      239.25             0.016        False       0.109       0.238
7   20190925_01   AVFR      239.25             0.016        False       0.125       0.247
8   20190925_04   AVFL      238.75             0.016        False       0.121       0.138
9   20190925_04   AVFR      238.75             0.016        False       0.117       0.184
10  20190928_03   AVFR      239.00             0.250         True       0.268       0.289
11  20190928_05   AVFL      239.00             0.031        False       0.142       0.222
12  20190928_07   AVFL      238.50             0.016        False       0.164       0.252
13  20190928_07   AVFR      238.50             0.016        False       0.151       0.226
14  20190928_08   AVFL      238.75             0.234         True       0.113       0.155
15  20190928_08   AVFR      238.75             0.016        False       0.109       0.184
16  20190928_11   AVFL      239.25             0.016        False       0.113       0.184
17  20190928_11   AVFR      239.25             0.016        False       0.117       0.142
18  20190928_13   AVFL      239.00             0.016        False       0.155       0.301
19  20190928_13   AVFR      239.00             0.016        False       0.050       0.146
20  20190929_03   AVFL      239.00             0.016        False       0.079       0.163
21  20190929_03   AVFR      239.00             0.016        False       0.121       0.205
22  20190929_05   AVFL      239.00             0.016        False       0.105       0.167
23  20190929_06   AVFL      239.00             0.016        False       0.046       0.109
24  20190929_06   AVFR      239.00             0.016        False       0.142       0.192
25  20190929_07   AVFR      232.75             0.016        False       0.052       0.133
26  20191030_03   AVFL      239.00             0.016        False       0.121       0.238
27  20191030_07   AVFL      238.75             0.016        False       0.071       0.239
28  20191030_07   AVFR      238.75             0.016        False       0.230       0.272
29  20191104_08   AVFL      239.00             0.016        False       0.084       0.184
30  20191104_10   AVFL      239.00             0.016        False       0.180       0.280
31  20191104_10   AVFR      239.00             0.016        False       0.205       0.276

def corr_mats(d, window=30.0):
    allnames = sorted({n for a in d["animals"] for n in a["traces"]})
    ix = {n: i for i, n in enumerate(allnames)}
    M = {}
    for a in d["animals"]:
        names = sorted(a["traces"])
        X = np.array([a["traces"][n] for n in names]).T
        t = np.arange(len(X)) / a["rate"]
        for s in npl.STIMULI:
            s0 = a["stimuli"][s][0]
            C = np.corrcoef(X[(t >= s0) & (t < s0 + window)].T)
            full = np.full((len(allnames), len(allnames)), np.nan)
            ii = [ix[n] for n in names]
            full[np.ix_(ii, ii)] = C
            M[(a["id"], s)] = full
    return M, allnames

def mat_sim(A, B):
    iu = np.triu_indices_from(A, 1)
    a, b = A[iu], B[iu]
    ok = np.isfinite(a) & np.isfinite(b)
    return np.corrcoef(a[ok], b[ok])[0, 1]

SPAIRS = [("nacl", "butanone"), ("nacl", "pentanedione"), ("butanone", "pentanedione")]
for mode in MODES:
    M, allnames = corr_mats(data[mode])
    ids = [a["id"] for a in data[mode]["animals"]]
    sims = {}
    for i in range(len(ids)):
        for j in range(i + 1, len(ids)):
            for x in npl.STIMULI:
                for y in npl.STIMULI:
                    sims[(i, j, x, y)] = mat_sim(M[(ids[i], x)], M[(ids[j], y)])
    for x, y in SPAIRS:
        same = np.array([(sims[(i, j, x, x)] + sims[(i, j, y, y)]) / 2 for i in range(len(ids)) for j in range(i + 1, len(ids))])
        diff = np.array([(sims[(i, j, x, y)] + sims[(i, j, y, x)]) / 2 for i in range(len(ids)) for j in range(i + 1, len(ids))])
        w = stats.wilcoxon(same, diff, alternative="greater")
        print(f"[{mode}] (a) {x[:4]}/{y[:4]}: same {np.mean(same):.4f}, different {np.mean(diff):.4f}, "
              f"same > different in {np.sum(same > diff)} of {len(same)} animal pairs, Wilcoxon p={w.pvalue:.3g}")
    print(f"[{mode}] (a) mean similarity between animals, same stimulus:",
          {x[:4]: round(np.mean([sims[(i, j, x, x)] for i in range(len(ids)) for j in range(i + 1, len(ids))]), 4) for x in npl.STIMULI})
    # (b) random halves
    rng = np.random.default_rng(611)
    wins = collections.Counter()
    gaps = collections.defaultdict(list)
    for rep in range(100):
        perm = rng.permutation(len(ids))
        h1, h2 = perm[: len(ids) // 2], perm[len(ids) // 2:]
        with np.errstate(all="ignore"):
            mean1 = {s: np.nanmean([M[(ids[k], s)] for k in h1], axis=0) for s in npl.STIMULI}
            mean2 = {s: np.nanmean([M[(ids[k], s)] for k in h2], axis=0) for s in npl.STIMULI}
        for x, y in SPAIRS:
            same = (mat_sim(mean1[x], mean2[x]) + mat_sim(mean1[y], mean2[y])) / 2
            diff = (mat_sim(mean1[x], mean2[y]) + mat_sim(mean1[y], mean2[x])) / 2
            wins[(x, y)] += same > diff
            gaps[(x, y)].append(same - diff)
    for x, y in SPAIRS:
        print(f"[{mode}] (b) {x[:4]}/{y[:4]}: same > different in {wins[(x, y)]} of 100 splits, "
              f"mean gap {np.mean(gaps[(x, y)]):.4f} (range {np.min(gaps[(x, y)]):.4f} to {np.max(gaps[(x, y)]):.4f})")

[one] (a) nacl/buta: same 0.0558, different 0.0481, same > different in 126 of 210 animal pairs, Wilcoxon p=0.000634
[one] (a) nacl/pent: same 0.0487, different 0.0410, same > different in 131 of 210 animal pairs, Wilcoxon p=0.000138
[one] (a) buta/pent: same 0.0497, different 0.0467, same > different in 107 of 210 animal pairs, Wilcoxon p=0.172
[one] (a) mean similarity between animals, same stimulus: {'nacl': np.float64(0.0548), 'buta': np.float64(0.0567), 'pent': np.float64(0.0427)}

[one] (b) nacl/buta: same > different in 91 of 100 splits, mean gap 0.0222 (range -0.0217 to 0.0479)
[one] (b) nacl/pent: same > different in 96 of 100 splits, mean gap 0.0288 (range -0.0128 to 0.0564)
[one] (b) buta/pent: same > different in 95 of 100 splits, mean gap 0.0212 (range -0.0077 to 0.0526)

[deposit] (a) nacl/buta: same 0.0671, different 0.0614, same > different in 116 of 210 animal pairs, Wilcoxon p=0.00783
[deposit] (a) nacl/pent: same 0.0579, different 0.0514, same > different in 130 of 210 animal pairs, Wilcoxon p=0.00177
[deposit] (a) buta/pent: same 0.0596, different 0.0573, same > different in 102 of 210 animal pairs, Wilcoxon p=0.269
[deposit] (a) mean similarity between animals, same stimulus: {'nacl': np.float64(0.0654), 'buta': np.float64(0.0688), 'pent': np.float64(0.0505)}

[deposit] (b) nacl/buta: same > different in 50 of 100 splits, mean gap -0.0009 (range -0.0481 to 0.0575)
[deposit] (b) nacl/pent: same > different in 71 of 100 splits, mean gap 0.0087 (range -0.0377 to 0.0440)
[deposit] (b) buta/pent: same > different in 16 of 100 splits, mean gap -0.0084 (range -0.0329 to 0.0279)

[both] (a) nacl/buta: same 0.0507, different 0.0402, same > different in 131 of 210 animal pairs, Wilcoxon p=1.24e-05
[both] (a) nacl/pent: same 0.0452, different 0.0357, same > different in 128 of 210 animal pairs, Wilcoxon p=2.61e-05
[both] (a) buta/pent: same 0.0463, different 0.0443, same > different in 101 of 210 animal pairs, Wilcoxon p=0.337
[both] (a) mean similarity between animals, same stimulus: {'nacl': np.float64(0.0496), 'buta': np.float64(0.0519), 'pent': np.float64(0.0408)}

[both] (b) nacl/buta: same > different in 98 of 100 splits, mean gap 0.0299 (range -0.0094 to 0.0617)
[both] (b) nacl/pent: same > different in 98 of 100 splits, mean gap 0.0382 (range -0.0167 to 0.0677)
[both] (b) buta/pent: same > different in 99 of 100 splits, mean gap 0.0249 (range -0.0012 to 0.0529)

[none] (a) nacl/buta: same 0.0627, different 0.0573, same > different in 123 of 210 animal pairs, Wilcoxon p=0.00762
[none] (a) nacl/pent: same 0.0535, different 0.0473, same > different in 130 of 210 animal pairs, Wilcoxon p=0.00135
[none] (a) buta/pent: same 0.0553, different 0.0531, same > different in 110 of 210 animal pairs, Wilcoxon p=0.196
[none] (a) mean similarity between animals, same stimulus: {'nacl': np.float64(0.0608), 'buta': np.float64(0.0645), 'pent': np.float64(0.0462)}

[none] (b) nacl/buta: same > different in 56 of 100 splits, mean gap 0.0075 (range -0.0508 to 0.0663)
[none] (b) nacl/pent: same > different in 82 of 100 splits, mean gap 0.0175 (range -0.0315 to 0.0546)
[none] (b) buta/pent: same > different in 71 of 100 splits, mean gap 0.0103 (range -0.0277 to 0.0506)

def shifted_mats(d, offset=30.0, window=30.0):
    allnames = sorted({n for a in d["animals"] for n in a["traces"]})
    ix = {n: i for i, n in enumerate(allnames)}
    M = {}
    for a in d["animals"]:
        names = sorted(a["traces"])
        X = np.array([a["traces"][n] for n in names]).T
        t = np.arange(len(X)) / a["rate"]
        for s in npl.STIMULI:
            s0 = a["stimuli"][s][0] + offset
            C = np.corrcoef(X[(t >= s0) & (t < s0 + window)].T)
            full = np.full((len(allnames), len(allnames)), np.nan)
            ii = [ix[n] for n in names]
            full[np.ix_(ii, ii)] = C
            M[(a["id"], s)] = full
    return M

Mn = shifted_mats(D)
ids = [a["id"] for a in D["animals"]]
for x, y in SPAIRS:
    same, diff = [], []
    for i in range(len(ids)):
        for j in range(i + 1, len(ids)):
            same.append((mat_sim(Mn[(ids[i], x)], Mn[(ids[j], x)]) + mat_sim(Mn[(ids[i], y)], Mn[(ids[j], y)])) / 2)
            diff.append((mat_sim(Mn[(ids[i], x)], Mn[(ids[j], y)]) + mat_sim(Mn[(ids[i], y)], Mn[(ids[j], x)])) / 2)
    same, diff = np.array(same), np.array(diff)
    print(f"[one] null window {x[:4]}/{y[:4]}: same {same.mean():.4f}, different {diff.mean():.4f}, "
          f"Wilcoxon p={stats.wilcoxon(same, diff, alternative='greater').pvalue:.3g}")
rng = np.random.default_rng(612)
wins = collections.Counter()
for rep in range(100):
    perm = rng.permutation(len(ids))
    h1, h2 = perm[: len(ids) // 2], perm[len(ids) // 2:]
    with np.errstate(all="ignore"):
        m1 = {s: np.nanmean([Mn[(ids[k], s)] for k in h1], axis=0) for s in npl.STIMULI}
        m2 = {s: np.nanmean([Mn[(ids[k], s)] for k in h2], axis=0) for s in npl.STIMULI}
    for x, y in SPAIRS:
        wins[(x, y)] += (mat_sim(m1[x], m2[x]) + mat_sim(m1[y], m2[y])) / 2 > (mat_sim(m1[x], m2[y]) + mat_sim(m1[y], m2[x])) / 2
print("[one] null window (b): same > different in", {f"{x[:4]}/{y[:4]}": wins[(x, y)] for x, y in SPAIRS}, "of 100 splits")

[one] null window nacl/buta: same 0.0548, different 0.0552, Wilcoxon p=0.531

[one] null window nacl/pent: same 0.0592, different 0.0502, Wilcoxon p=7.49e-06

[one] null window buta/pent: same 0.0575, different 0.0535, Wilcoxon p=0.0134

[one] null window (b): same > different in {'nacl/buta': np.int64(40), 'nacl/pent': np.int64(100), 'buta/pent': np.int64(100)} of 100 splits

def class_traces(d, norm="zscore"):
    out = {}
    for a in d["animals"]:
        g = collections.defaultdict(list)
        for n, tr in a["traces"].items():
            g[cls_of[n]].append(npl.normalise(tr, norm))
        out[a["id"]] = {c: np.mean(v, axis=0) for c, v in g.items()}
    return out

def class_scale(CT):
    """Each class's traces divided by the SD of that class's values pooled over animals."""
    sd = collections.defaultdict(list)
    for d_ in CT.values():
        for c, v in d_.items():
            sd[c].append(v[4:])
    sd = {c: np.nanstd(np.concatenate(v)) for c, v in sd.items()}
    return {k: {c: v / sd[c] for c, v in d_.items()} for k, d_ in CT.items()}

for mode, norm in [(m, "z-scored") for m in MODES] + [("one", "dF/F0"), ("one", "class-scaled dF/F0"),
                                                     ("none", "dF/F0"), ("none", "class-scaled dF/F0")]:
    CT = class_traces(data[mode], "zscore" if norm == "z-scored" else "raw")
    if norm == "class-scaled dF/F0":
        CT = class_scale(CT)
    common = sorted(set.intersection(*[set(v) for v in CT.values()]))
    nfr = min(len(next(iter(v.values()))) for v in CT.values())
    A = {k: np.array([v[c][:nfr] for c in common]).T for k, v in CT.items()}
    stim_win = {}
    for a in data[mode]["animals"]:
        t = np.arange(nfr) / a["rate"]
        for s in npl.STIMULI:
            s0 = a["stimuli"][s][0]
            stim_win[(a["id"], s)] = A[a["id"]][(t >= s0) & (t < s0 + 30)][:120]
    for reading in ("a", "b"):
        if reading == "a":
            X = np.vstack([A[k][4:] for k in A])
            X = X[np.all(np.isfinite(X), axis=1)]
        else:
            X = np.vstack(list(stim_win.values()))
            X = X[np.all(np.isfinite(X), axis=1)]
        mu = X.mean(0)
        U, S, Vt = np.linalg.svd(X - mu, full_matrices=False)
        ev = S ** 2 / np.sum(S ** 2)
        P = {k: (w - mu) @ Vt[:3].T for k, w in stim_win.items()}
        mean = {s: np.mean([P[(a, s)] for a in A], axis=0) for s in npl.STIMULI}
        sem = {s: np.linalg.norm(stats.sem([P[(a, s)] for a in A], axis=0), axis=1) for s in npl.STIMULI}
        res = []
        for x, y in SPAIRS:
            dist = np.linalg.norm(mean[x] - mean[y], axis=1).mean()
            thr = 2 * np.mean([sem[x].mean(), sem[y].mean()])
            res.append(f"{x[:4]}/{y[:4]} dist {dist:.3f} vs 2xSEM {thr:.3f} -> {'distinguishable' if dist > thr else 'not'}")
        top = [common[i] for i in np.argsort(-np.abs(Vt[0]))[:3]]
        print(f"[{mode}, {norm}] ({reading}) {len(common)} classes, variance explained by PC1-3: {np.round(ev[:3], 3)}, "
              f"largest PC1 loadings {top}; " + "; ".join(res))

[one, z-scored] (a) 48 classes, variance explained by PC1-3: [0.276 0.114 0.046], largest PC1 loadings ['NSM', 'URAD', 'M3']; nacl/buta dist 1.049 vs 2xSEM 1.479 -> not; nacl/pent dist 1.701 vs 2xSEM 1.328 -> distinguishable; buta/pent dist 1.339 vs 2xSEM 1.445 -> not
[one, z-scored] (b) 48 classes, variance explained by PC1-3: [0.173 0.121 0.067], largest PC1 loadings ['AIB', 'AFD', 'URYV']; nacl/buta dist 1.067 vs 2xSEM 1.473 -> not; nacl/pent dist 1.724 vs 2xSEM 1.316 -> distinguishable; buta/pent dist 1.302 vs 2xSEM 1.432 -> not

[deposit, z-scored] (a) 71 classes, variance explained by PC1-3: [0.278 0.164 0.045], largest PC1 loadings ['URB', 'ASK', 'OLL']; nacl/buta dist 1.314 vs 2xSEM 1.941 -> not; nacl/pent dist 2.143 vs 2xSEM 1.745 -> distinguishable; buta/pent dist 1.869 vs 2xSEM 1.845 -> distinguishable
[deposit, z-scored] (b) 71 classes, variance explained by PC1-3: [0.209 0.143 0.07 ], largest PC1 loadings ['AVA', 'AVE', 'RMDV']; nacl/buta dist 1.385 vs 2xSEM 1.916 -> not; nacl/pent dist 2.730 vs 2xSEM 1.718 -> distinguishable; buta/pent dist 2.175 vs 2xSEM 1.816 -> distinguishable
[both, z-scored] (a) 29 classes, variance explained by PC1-3: [0.306 0.09  0.056], largest PC1 loadings ['NSM', 'M3', 'AUA']; nacl/buta dist 0.787 vs 2xSEM 1.155 -> not; nacl/pent dist 1.364 vs 2xSEM 1.072 -> distinguishable; buta/pent dist 0.958 vs 2xSEM 1.163 -> not
[both, z-scored] (b) 29 classes, variance explained by PC1-3: [0.203 0.102 0.096], largest PC1 loadings ['M5', 'AUA', 'NSM']; nacl/buta dist 1.059 vs 2xSEM 1.156 -> not; nacl/pent dist 1.600 vs 2xSEM 1.066 -> distinguishable; buta/pent dist 0.902 vs 2xSEM 1.160 -> not

[none, z-scored] (a) 81 classes, variance explained by PC1-3: [0.274 0.159 0.042], largest PC1 loadings ['URB', 'ASK', 'RIP']; nacl/buta dist 1.337 vs 2xSEM 2.049 -> not; nacl/pent dist 2.211 vs 2xSEM 1.834 -> distinguishable; buta/pent dist 1.942 vs 2xSEM 1.951 -> not
[none, z-scored] (b) 81 classes, variance explained by PC1-3: [0.199 0.148 0.062], largest PC1 loadings ['AVA', 'AVE', 'RMDV']; nacl/buta dist 1.397 vs 2xSEM 2.024 -> not; nacl/pent dist 2.716 vs 2xSEM 1.811 -> distinguishable; buta/pent dist 2.246 vs 2xSEM 1.921 -> distinguishable
[one, dF/F0] (a) 48 classes, variance explained by PC1-3: [0.757 0.054 0.036], largest PC1 loadings ['RIM', 'ASJ', 'RID']; nacl/buta dist 2.220 vs 2xSEM 2.562 -> not; nacl/pent dist 2.697 vs 2xSEM 2.725 -> not; buta/pent dist 1.117 vs 2xSEM 1.440 -> not
[one, dF/F0] (b) 48 classes, variance explained by PC1-3: [0.751 0.056 0.031], largest PC1 loadings ['RIM', 'AIY', 'ASJ']; nacl/buta dist 2.234 vs 2xSEM 2.567 -> not; nacl/pent dist 2.794 vs 2xSEM 2.747 -> distinguishable; buta/pent dist 1.145 vs 2xSEM 1.458 -> not

[one, class-scaled dF/F0] (a) 48 classes, variance explained by PC1-3: [0.176 0.117 0.094], largest PC1 loadings ['URAV', 'M3', 'NSM']; nacl/buta dist 1.609 vs 2xSEM 1.868 -> not; nacl/pent dist 1.404 vs 2xSEM 1.597 -> not; buta/pent dist 1.229 vs 2xSEM 1.474 -> not
[one, class-scaled dF/F0] (b) 48 classes, variance explained by PC1-3: [0.173 0.109 0.09 ], largest PC1 loadings ['CEPV', 'IL1V', 'AIY']; nacl/buta dist 1.651 vs 2xSEM 1.894 -> not; nacl/pent dist 1.632 vs 2xSEM 1.712 -> not; buta/pent dist 1.265 vs 2xSEM 1.581 -> not
[none, dF/F0] (a) 81 classes, variance explained by PC1-3: [0.493 0.176 0.069], largest PC1 loadings ['SIBV', 'RIM', 'RMEV']; nacl/buta dist 2.162 vs 2xSEM 2.887 -> not; nacl/pent dist 2.692 vs 2xSEM 2.851 -> not; buta/pent dist 2.550 vs 2xSEM 2.200 -> distinguishable

[none, dF/F0] (b) 81 classes, variance explained by PC1-3: [0.464 0.215 0.068], largest PC1 loadings ['SIBV', 'RIM', 'RMEV']; nacl/buta dist 2.194 vs 2xSEM 2.909 -> not; nacl/pent dist 2.809 vs 2xSEM 2.905 -> not; buta/pent dist 2.605 vs 2xSEM 2.237 -> distinguishable
[none, class-scaled dF/F0] (a) 81 classes, variance explained by PC1-3: [0.178 0.151 0.084], largest PC1 loadings ['M3', 'RIP', 'URB']; nacl/buta dist 1.914 vs 2xSEM 2.440 -> not; nacl/pent dist 2.045 vs 2xSEM 2.106 -> not; buta/pent dist 2.290 vs 2xSEM 2.104 -> distinguishable

[none, class-scaled dF/F0] (b) 81 classes, variance explained by PC1-3: [0.171 0.148 0.082], largest PC1 loadings ['URYV', 'RMDV', 'CEPV']; nacl/buta dist 2.053 vs 2xSEM 2.462 -> not; nacl/pent dist 2.358 vs 2xSEM 2.234 -> distinguishable; buta/pent dist 2.360 vs 2xSEM 2.193 -> distinguishable

chem, gap = npl.load_cook()
print("Cook et al. synapse-count matrices: chemical", chem.shape, "total", chem.values.sum(), "| gap junction", gap.shape,
      "total", gap.values.sum(), "| gap matrix symmetric:", bool(np.allclose(gap.values, gap.values.T)))
dnames = sorted({n for r in recs for n in r["labels"] if n})
unmatched = [n for n in dnames if n not in chem.index or n not in gap.index]
print("dandiset names:", len(dnames), "| matched:", len(dnames) - len(unmatched), "| unmatched:", unmatched)

Cook et al. synapse-count matrices: chemical (468, 468) total 11760.0 | gap junction (466, 466) total 4094.0 | gap matrix symmetric: True
dandiset names: 193 | matched: 190 | unmatched: ['AMSOL', 'AMSOR', 'HMC']

def animal_corr(d, which="all"):
    """Per-animal correlation matrices on a common name index; which = all, stim or nonstim frames."""
    names = sorted({n for a in d["animals"] for n in a["traces"]} & set(chem.index) & set(gap.index))
    ix = {n: i for i, n in enumerate(names)}
    mats = []
    for a in d["animals"]:
        nm = sorted(set(a["traces"]) & set(ix))
        X = np.array([a["traces"][n] for n in nm]).T[4:]
        t = (np.arange(len(X)) + 4) / a["rate"]
        st = np.zeros(len(X), bool)
        for s0, s1 in a["stimuli"].values():
            st |= (t >= s0) & (t < s1)
        if which == "stim":
            X = X[st]
        elif which == "nonstim":
            X = X[~st]
        full = np.full((len(names), len(names)), np.nan)
        ii = [ix[n] for n in nm]
        full[np.ix_(ii, ii)] = np.corrcoef(X.T)
        mats.append(full)
    return names, np.array(mats)

def r2_table(names, mats, chem_m, gap_m, corr="pearson", logc=False, min_animals=0):
    iu = np.triu_indices(len(names), 1)
    G = gap_m.loc[names, names].values
    Cm = chem_m.loc[names, names].values
    syn = {"electrical": G[iu], "chemical (a) sum": (Cm + Cm.T)[iu], "chemical (b) max": np.maximum(Cm, Cm.T)[iu]}
    present = np.isfinite(mats[:, iu[0], iu[1]])
    navail = present.sum(0)
    with np.errstate(all="ignore"):
        fa = np.abs(np.nanmean(mats[:, iu[0], iu[1]], axis=0))
    out = {}
    for kind, cnt in syn.items():
        cnt = np.log1p(cnt) if logc else cnt
        for fread in ("a", "b"):
            for pread in ("a", "b"):
                keep = (navail > min_animals) & (navail > 0)
                if pread == "b":
                    keep &= cnt > 0
                if fread == "a":
                    f, c = fa[keep], cnt[keep]
                else:
                    sel = present[:, keep]
                    f = np.abs(mats[:, iu[0][keep], iu[1][keep]][sel])
                    c = np.broadcast_to(cnt[keep], sel.shape)[sel]
                r = stats.pearsonr(f, c)[0] if corr == "pearson" else stats.spearmanr(f, c)[0]
                out[(kind, f"functional ({fread})", f"pairs ({pread})")] = (100 * r ** 2, len(f))
    return out

cn_names, cn_mats = animal_corr(D, "all")
print("neurons in the comparison:", len(cn_names))
base = r2_table(cn_names, cn_mats, chem, gap)
tab = pd.DataFrame([{"synapses": k[0], "functional": k[1], "pairs": k[2], "R2_percent": v[0], "n": v[1]} for k, v in base.items()])
print(tab.round(3).to_string(index=False))
for mode in ("deposit", "both", "none"):
    nm, mt = animal_corr(data[mode], "all")
    alt = r2_table(nm, mt, chem, gap)
    print(f"[{mode}]", {f"{k[0][:12]} {k[1][-3:]} {k[2][-3:]}": round(v[0], 3) for k, v in alt.items()})
chem5, gap5 = npl.load_cook(npl.COOK_ADJACENCY, "hermaphrodite chemical", "herm gap jn symmetric")
alt5 = r2_table(cn_names, cn_mats, chem5, gap5)
print("[SI 5 weighted adjacency, a variant]", {f"{k[0][:12]} {k[1][-3:]} {k[2][-3:]}": round(v[0], 3) for k, v in alt5.items()})

neurons in the comparison: 190
        synapses     functional     pairs  R2_percent      n
      electrical functional (a) pairs (a)       1.282  17869
      electrical functional (a) pairs (b)       3.465    478
      electrical functional (b) pairs (a)       0.342 260172
      electrical functional (b) pairs (b)       1.078   6966
chemical (a) sum functional (a) pairs (a)       0.231  17869
chemical (a) sum functional (a) pairs (b)       0.203   1688
chemical (a) sum functional (b) pairs (a)       0.117 260172
chemical (a) sum functional (b) pairs (b)       0.127  25707
chemical (b) max functional (a) pairs (a)       0.173  17869
chemical (b) max functional (a) pairs (b)       0.037   1688
chemical (b) max functional (b) pairs (a)       0.102 260172
chemical (b) max functional (b) pairs (b)       0.067  25707
[deposit] {'electrical (a) (a)': np.float64(1.457), 'electrical (a) (b)': np.float64(3.68), 'electrical (b) (a)': np.float64(0.451), 'electrical (b) (b)': np.float64(1.221), 'chemical (a) (a) (a)': np.float64(0.314), 'chemical (a) (a) (b)': np.float64(0.248), 'chemical (a) (b) (a)': np.float64(0.165), 'chemical (a) (b) (b)': np.float64(0.178), 'chemical (b) (a) (a)': np.float64(0.251), 'chemical (b) (a) (b)': np.float64(0.07), 'chemical (b) (b) (a)': np.float64(0.147), 'chemical (b) (b) (b)': np.float64(0.104)}

[both] {'electrical (a) (a)': np.float64(1.375), 'electrical (a) (b)': np.float64(3.538), 'electrical (b) (a)': np.float64(0.312), 'electrical (b) (b)': np.float64(0.829), 'chemical (a) (a) (a)': np.float64(0.212), 'chemical (a) (a) (b)': np.float64(0.182), 'chemical (a) (b) (a)': np.float64(0.098), 'chemical (a) (b) (b)': np.float64(0.088), 'chemical (b) (a) (a)': np.float64(0.162), 'chemical (b) (a) (b)': np.float64(0.038), 'chemical (b) (b) (a)': np.float64(0.087), 'chemical (b) (b) (b)': np.float64(0.045)}

[none] {'electrical (a) (a)': np.float64(1.354), 'electrical (a) (b)': np.float64(3.367), 'electrical (b) (a)': np.float64(0.405), 'electrical (b) (b)': np.float64(1.155), 'chemical (a) (a) (a)': np.float64(0.294), 'chemical (a) (a) (b)': np.float64(0.225), 'chemical (a) (b) (a)': np.float64(0.152), 'chemical (a) (b) (b)': np.float64(0.172), 'chemical (b) (a) (a)': np.float64(0.234), 'chemical (b) (a) (b)': np.float64(0.058), 'chemical (b) (b) (a)': np.float64(0.136), 'chemical (b) (b) (b)': np.float64(0.103)}

[SI 5 weighted adjacency, a variant] {'electrical (a) (a)': np.float64(1.093), 'electrical (a) (b)': np.float64(3.483), 'electrical (b) (a)': np.float64(0.376), 'electrical (b) (b)': np.float64(1.807), 'chemical (a) (a) (a)': np.float64(0.192), 'chemical (a) (a) (b)': np.float64(0.061), 'chemical (a) (b) (a)': np.float64(0.113), 'chemical (a) (b) (b)': np.float64(0.089), 'chemical (b) (a) (a)': np.float64(0.14), 'chemical (b) (a) (b)': np.float64(0.001), 'chemical (b) (b) (a)': np.float64(0.097), 'chemical (b) (b) (b)': np.float64(0.044)}

var = {}
for th in (5, 10, 15):
    var[f"pairs in > {th} animals"] = r2_table(cn_names, cn_mats, chem, gap, min_animals=th)
var["Spearman"] = r2_table(cn_names, cn_mats, chem, gap, corr="spearman")
var["log(1 + count)"] = r2_table(cn_names, cn_mats, chem, gap, logc=True)
for which in ("stim", "nonstim"):
    nm, mt = animal_corr(D, which)
    var[f"{which} frames only"] = r2_table(nm, mt, chem, gap)
vt = pd.DataFrame([{"variation": v, "synapses": k[0], "functional": k[1], "pairs": k[2], "R2_percent": x[0], "n": x[1]}
                   for v, tbl in var.items() for k, x in tbl.items()])
print(vt.pivot_table(index=["synapses", "functional", "pairs"], columns="variation", values="R2_percent").round(3).to_string())
first = vt[(vt.functional == "functional (a)") & (vt.pairs == "pairs (a)") & (vt.synapses != "chemical (b) max")]
print("largest R2 under the first-listed readings: %.3f%% (%s)" % (first.R2_percent.max(),
      first.loc[first.R2_percent.idxmax(), ["variation", "synapses"]].tolist()))
print("largest R2 over all readings and variations: %.3f%%" % vt.R2_percent.max())

variation                                  Spearman  log(1 + count)  nonstim frames only  pairs in > 10 animals  pairs in > 15 animals  pairs in > 5 animals  stim frames only
synapses         functional     pairs                                                                                                                                         
chemical (a) sum functional (a) pairs (a)     0.202           0.324                0.178                  0.309                  0.300                 0.312             0.191
                                pairs (b)     0.030           0.133                0.099                  0.256                  0.260                 0.256             0.520
                 functional (b) pairs (a)     0.123           0.165                0.095                  0.117                  0.080                 0.118             0.058
                                pairs (b)     0.061           0.108                0.095                  0.123                  0.068                 0.132             0.113
chemical (b) max functional (a) pairs (a)     0.198           0.276                0.129                  0.225                  0.208                 0.233             0.149
                                pairs (b)     0.011           0.010                0.005                  0.046                  0.026                 0.050             0.264
                 functional (b) pairs (a)     0.121           0.152                0.083                  0.101                  0.066                 0.103             0.050
                                pairs (b)     0.023           0.050                0.049                  0.063                  0.022                 0.070             0.069
electrical       functional (a) pairs (a)     0.677           1.493                1.200                  1.508                  1.217                 1.443             0.730
                                pairs (b)     2.434           3.743                3.486                  3.329                  2.325                 3.061             1.987
                 functional (b) pairs (a)     0.250           0.406                0.313                  0.348                  0.281                 0.342             0.081
                                pairs (b)     0.865           1.240                1.091                  1.120                  1.093                 1.061             0.148
largest R2 under the first-listed readings: 1.508% (['pairs in > 10 animals', 'electrical'])
largest R2 over all readings and variations: 3.743%

present = {"acquisition": sorted({k for r in recs[:1] for k in ("CalciumImageSeries", "NeuroPALImageRaw")}),
           "strains": sorted({r["strain"] for r in recs}), "descriptions": sorted({r["session_description"] for r in recs})}
print(present)
print("runtime so far: %.1f min" % ((time.time() - T_START) / 60))

{'acquisition': ['CalciumImageSeries', 'NeuroPALImageRaw'], 'strains': ['OH16230'], 'descriptions': ['C. elegans head NeuroPAL and Calcium imaging']}
runtime so far: 6.7 min

Claim	Tested here
link-zenodo-traces	tested (link)
fig6-sample-sizes	tested (heads only; tails not in the dandiset)
fig6-all-heads-responsive	tested
fig6d-ase-nacl	tested
fig6-nacl-unknown-interneurons	tested
fig6c-odors-more-inhibitory	tested
fig6d-awc-asymmetric	blocked: the dandiset does not label AWC ON and OFF (variants with the Zenodo labels and with labels read from the NeuroPAL image are shown)
fig6-asj-ric	tested
fig6-interneurons-downstream	tested
fig6-tail-responses	blocked: no tail recordings in the dandiset
fig6f-pharyngeal	tested
fig6-stimulus-sets-distinct	tested
figs1h-avf-cyclic	tested
fig6gi-stimulus-specific-correlations	tested
figs2a-pca-trajectories	tested
link-connectome-names	tested (link)
fig6jk-connectome-correlation	tested (heads only)
fig6-connectome-variations	tested (heads only)
fig1, figs1, fig2, fig4c, fig5 and fig7bc claims (8)	blocked: the data are not in the dandiset
fig4bd-gaba-network, fig4e-extrasynaptic, fig7de-optimal-coloring	blocked: outside inputs or software, not tested here (see the attempt record)

Claim	Outcome
link-zenodo-traces	partially reproduced: files and stimuli match, 81.7% of traces at r of 0.9 (85.8% with the one-frame offset)
link-connectome-names	reproduced: 190 of 193 names, the 3 others not neurons
fig6-sample-sizes	partially reproduced (assumed): head values match under the main selection reading, no tails
fig6-all-heads-responsive	not reproduced: 19 of 21 heads
fig6d-ase-nacl	reproduced
fig6-nacl-unknown-interneurons	reproduced on z-scored traces (assumed); 4 of 6 on dF/F0
fig6c-odors-more-inhibitory	reproduced, thin margin for 2-butanone
fig6d-awc-asymmetric	blocked: metadata insufficient (no ON/OFF label; the NeuroPAL image carries the side)
fig6-asj-ric	partially reproduced (assumed): 5 of 6 parts; 6 of 6 under two selection readings, 3 on dF/F0
fig6-interneurons-downstream	reproduced
fig6-tail-responses	blocked: data missing
fig6f-pharyngeal	reproduced
fig6-stimulus-sets-distinct	reproduced (a weak test)
figs1h-avf-cyclic	not reproduced: 0.12 to 0.20 Hz
fig6gi-stimulus-specific-correlations	partially reproduced (assumed): the odors are not told apart; the null window also passes in part
figs2a-pca-trajectories	not reproduced: at most two of three pairs
fig6jk-connectome-correlation	partially reproduced: all R2 below 5%, none in tolerance, heads only
fig6-connectome-variations	reproduced: largest R2 1.51%
fig1, figs1, fig2, fig4c, fig5, fig7bc claims (8)	blocked: data missing
fig4bd-gaba-network, fig4e-extrasynaptic	blocked (other): not tested in this attempt; the inputs the paper names (Table S4, the cited transmitter maps, the Cook et al. connectome) are public
fig7de-optimal-coloring	blocked: not runnable (MATLAB solver)

NeuroPAL whole-brain activity (Yemini et al. 2021) from dandiset 000541¶

Data check¶

Trace processing¶

Clock of the stimulus table¶

link-zenodo-traces¶

fig6-sample-sizes¶

fig6-all-heads-responsive¶

The response tests used by the claims below¶

The tests against the paper's Table S6¶

fig6d-ase-nacl¶

fig6-nacl-unknown-interneurons¶

fig6c-odors-more-inhibitory¶

fig6d-awc-asymmetric (blocked)¶

fig6-asj-ric¶

fig6-interneurons-downstream¶

fig6-tail-responses (blocked)¶

fig6f-pharyngeal¶

fig6-stimulus-sets-distinct¶

figs1h-avf-cyclic¶

fig6gi-stimulus-specific-correlations¶

figs2a-pca-trajectories¶

link-connectome-names¶

fig6jk-connectome-correlation¶

fig6-connectome-variations¶

Claims not testable on this dandiset¶

Summary¶