import collections
import hashlib
import time
import warnings

import numpy as np
import pandas as pd
import scipy.stats as st
import yaml
import matplotlib.pyplot as plt
from lifelines import KaplanMeierFitter
from statsmodels.stats.multitest import multipletests

import dandi_io
import mouse_eye as me

T0 = time.time()
warnings.filterwarnings("ignore", message="Mean of empty slice")   # NaN frames in the DeepLabCut traces
pd.set_option("display.width", 200)
np.set_printoptions(precision=3, suppress=True, linewidth=160)

CLAIMS = yaml.safe_load(open("../../claims/001347/10.1126-science.adt3971.yaml"))
CLAIM = {c["id"]: c for c in CLAIMS["claims"]}


def expected(cid):
    # Expected values of a claim, read from the claims file, for comparison only.
    return [(e["quantity"], e["value"], e["tolerance"]) for e in CLAIM[cid].get("expected", [])]


def bh(p):
    return multipletests(p, method="fdr_bh")[1]


def within(x, v, tol):
    return abs(x - v) <= tol + 1e-12


print("dandiset", dandi_io.DANDISET, dandi_io.VERSION)

dandiset 001347 0.250528.0702

zf = me.fetch_code_zip()
LISTS = me.code_lists(zf)
for script, v in LISTS.items():
    keys = {k: v[k] for k in ("ntrials",) if k in v}
    if "experiments" in v:
        print(f"{script}: ntrials={keys.get('ntrials')}")
        for cond, ids in v["experiments"].items():
            print(f"    {cond:9s} {ids}")
    if "SKIP_PUFFS" in v:
        print(f"    SKIP_PUFFS {v['SKIP_PUFFS']}")
    if "expids" in v:
        print(f"{script}: expids {v['expids']}")

fig1_drugs: ntrials=40
    ket50     [275, 276, 277, 280, 281]
    pcp20     [287, 288, 289, 290, 291]
    sal       [282, 283, 284, 285, 286]
    kx        [316, 338, 339, 343, 344]
    dzp       [325, 326, 327, 333, 334]
    bup       [340, 341, 342, 347, 348]
fig1_drugs_first_puff: ntrials=40
    ket50     [275, 276, 277, 280, 281]
    pcp20     [287, 288, 289, 290, 291]
    sal       [282, 283, 284, 285, 286]
    kx        [316, 338, 339, 343, 344]
    dzp       [325, 326, 327, 333, 334]
    bup       [340, 341, 342, 347, 348]
fig1_dose: ntrials=40
    ket50     [275, 276, 277, 280, 281]
    pcp20     [287, 288, 289, 290, 291]
    sal       [282, 283, 284, 285, 286]
    kx        [316, 338, 339, 343, 344]
    ket25     [355, 356, 357, 362, 363]
    pcp10     [358, 359, 360, 361, 364]
    ket13     [401, 400, 399, 398, 397]
    ket6      [403, 405, 404, 406, 426]
    pcp5      [370, 368, 366, 365, 434]
    pcp1      [430, 432, 431, 435, 436]
    ket35     [1020, 1021, 1022, 1023, 1024, 1025, 1026, 1027]
fig1_psi: ntrials=20
    18 PSI    [1004, 1001, 1002, 1003, 1030, 1031, 1029, 1028]
    8 PSI     [1034, 1033, 1032, 1035, 1008, 1005, 1007, 1006]
    5 PSI     [1014, 1016, 1019, 1017, 1011, 1012, 1013, 1009]
    28 PSI    [1020, 1021, 1022, 1023, 1024, 1025, 1026, 1027]
fig1_days: ntrials=40
    day1      [668, 669, 670, 671, 674]
    day2      [672, 673, 677, 678, 679]
    day3      [675, 676, 682, 683, 684]
    day4      [680, 681, 687, 688, 689]
    day5      [685, 686, 690, 691, 692]
fig1_sex: ntrials=40
    female    [282, 283, 284, 285, 286]
    male      [700, 701, 702, 708, 709]
fig1_muscimol: ntrials=20
    intact    [975, 971, 972, 973, 974]
    sal       [984, 985, 980, 981, 982]
    vpm_musc  [989, 990, 986, 987, 988]
fig1_persistence: ntrials=40
    ket50     [275, 276, 277, 280, 281]
    pcp20     [287, 288, 289, 290, 291]
    sal       [282, 283, 284, 285, 286]
    kx        [316, 338, 339, 343, 344]
    dzp       [325, 326, 327, 333, 334]
    bup       [340, 341, 342, 347, 348]
fig1_skippuff: ntrials=40
    ket50     [411, 412, 413, 420, 421]
    SKIP_PUFFS [6, 6, 5, 6, 5, 5, 5, 6, 6, 5, 5, 6, 5, 6, 5, 6, 6, 5, 6, 5]
fig1_npx_behavior: expids ['m94-549', 'm304-876', 'm11208-959', 'm11208-962', 'm11210-960', 'm11219-1108', 'm11221-1119', 'm11223-1125', 'm11225-1131', 'm11227-1132']
convert_npx_mouse_airpuff_to_nwb: expids ['20220228_m94_549', '20220831_m304_876', '20230218_m11208_959', '20230226_m11208_962', '20230218_m11210_960', '20230823_m11219_1108', '20230901_m11221_1119', '20230902_m11223_1125', '20230907_m11225_1131', '20230908_m11227_1132']

IDS = me.mouse_ids()
DATA, META = me.load_all(IDS)
meta = pd.DataFrame(META).T
meta["n_frames"] = [len(DATA[e]) for e in meta.index]
meta["n_trials"] = meta["n_frames"] / me.FRAMES_PER_TRIAL
meta["min"] = [DATA[e].min() for e in meta.index]
meta["max"] = [DATA[e].max() for e in meta.index]
meta["n_nan"] = [int(np.isnan(DATA[e]).sum()) for e in meta.index]
print("files:", len(IDS))
print("identifier equals file session id:", all(meta.loc[e, "identifier"] == str(e) for e in IDS))
print("subject_id equals session id:", all(meta.loc[e, "subject_id"] == str(e) for e in IDS))
for col in ["rate", "unit", "sex", "age", "subject_description"]:
    print(f"{col}: {dict(collections.Counter(meta[col]))}")
print("session_description examples:", list(meta["session_description"][:3]))
print("trials per file:", dict(collections.Counter(meta["n_trials"])))
print("whole trials only:", bool((meta["n_frames"] % me.FRAMES_PER_TRIAL == 0).all()))
print("20-trial sessions:", [e for e in IDS if meta.loc[e, "n_trials"] == 20])
print(f"eye closure range {meta['min'].min():.3f} to {meta['max'].max():.3f} (correlation with a "
      f"closed-eye template, no unit); missing values {meta['n_nan'].sum()}")

files: 142
identifier equals file session id: True
subject_id equals session id: True
rate: {200.0: 142}
unit: {'Eye closure, sampled at 200.0Hz': 142}
sex: {'O': 142}
age: {'P8W/': 142}
subject_description: {'mouse_airpuff': 142}
session_description examples: ['Airpuff session 275', 'Airpuff session 276', 'Airpuff session 277']
trials per file: {40.0: 103, 20.0: 39}
whole trials only: True
20-trial sessions: [971, 972, 973, 974, 975, 980, 981, 982, 984, 985, 986, 987, 988, 989, 990, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1011, 1012, 1013, 1014, 1016, 1017, 1019, 1028, 1029, 1030, 1031, 1032, 1033, 1034, 1035]
eye closure range -0.665 to 0.997 (correlation with a closed-eye template, no unit); missing values 0

h_all = collections.defaultdict(list); h_20 = collections.defaultdict(list); h_trial = collections.defaultdict(list)
for e in IDS:
    x = DATA[e]
    h_all[hashlib.md5(x.tobytes()).hexdigest()].append(e)
    h_20[hashlib.md5(x[:20 * me.FRAMES_PER_TRIAL].tobytes()).hexdigest()].append(e)
    for i in range(len(x) // me.FRAMES_PER_TRIAL):
        h_trial[hashlib.md5(x[i * me.FRAMES_PER_TRIAL:(i + 1) * me.FRAMES_PER_TRIAL].tobytes()).hexdigest()].append((e, i))
print("identical files:", [v for v in h_all.values() if len(v) > 1])
print("identical first 20 trials:", [v for v in h_20.values() if len(v) > 1])
print("identical single trials across or within files:", [v for v in h_trial.values() if len(v) > 1])
X = np.stack([DATA[e][:20 * me.FRAMES_PER_TRIAL] for e in IDS])
X = X - X.mean(1, keepdims=True); X /= np.linalg.norm(X, axis=1, keepdims=True)
C = X @ X.T; np.fill_diagonal(C, -1)
order = np.dstack(np.unravel_index(np.argsort(C, axis=None)[::-1], C.shape))[0]
top = [(IDS[a], IDS[b], round(float(C[a, b]), 3)) for a, b in order[:10] if a < b]
print("highest correlations between sessions (first 20 trials):", top)
del X, C

identical files: []
identical first 20 trials: []
identical single trials across or within files: []

highest correlations between sessions (first 20 trials): [(1023, 1026, 0.912), (1026, 1029, 0.889), (681, 1023, 0.874), (343, 1023, 0.871), (343, 681, 0.864)]

def jump(tr, t0):
    return tr[:, me.window(t0, (0.1, 0.2))].mean(1) - tr[:, me.window(t0, (-0.1, 0.0))].mean(1)

rows = []
for e in IDS:
    nt = int(meta.loc[e, "n_trials"])
    tr = me.prox_norm(me.trials(DATA[e], nt))[:20]
    at = np.mean([jump(tr, 3 * k).mean() for k in range(8)])
    mid = np.mean([jump(tr, 3 * k + 1.5).mean() for k in range(7)])
    rows.append((e, at, mid))
clock = pd.DataFrame(rows, columns=["session", "jump_at_puffs", "jump_between_puffs"]).set_index("session")
print(clock.describe().round(3))
print("sessions where the response at puff times is not larger than between puffs:",
      list(clock.index[clock.jump_at_puffs <= clock.jump_between_puffs]))

       jump_at_puffs  jump_between_puffs
count        142.000             142.000
mean           0.447              -0.005
std            0.168               0.004
min            0.058              -0.015
25%            0.324              -0.007
50%            0.428              -0.004
75%            0.558              -0.002
max            1.041               0.006
sessions where the response at puff times is not larger than between puffs: []

NPX = {p: me.npx_load(p) for p in me.npx_paths()}
b0 = NPX[me.npx_paths()[0]]
for k, v in b0.items():
    if isinstance(v, np.ndarray):
        print(f"  {k:22s} shape {v.shape}" + (f" value {v}" if v.size <= 8 else ""))
rows = []
for p, b in NPX.items():
    fps = int(b["fps"][0]); on = b["events-on"]; bt = b["eyeclosure_by_trial"].T
    start = on[0] - int(b["w_pre"][0]) * fps
    seg = b["eyeclosure"][start:start + bt.shape[1]]
    ok = np.allclose(seg, bt[0], equal_nan=True)
    gaps = np.diff(on[:40]) / fps
    rows.append(dict(file=p.split("/")[0], subject=b["subject_id"], n_onsets=len(on),
                     n_puff_trials=b["eyeclosure_by_puff"].shape[0] // 8, ket_trial=int(b["ket_trial"][0]),
                     pre=f"{b['trial_range_preket'][0]}-{b['trial_range_preket'][-1]}",
                     ket=f"{b['trial_range_ket'][0]}-{b['trial_range_ket'][-1]}", fps=fps,
                     ref_w=tuple(b["ref_w"]), aff_w=tuple(b["aff_w"]), led_delay=int(b["led_delay"][0]),
                     by_trial_matches=ok, t_trial_first=b["t_trial_eyeclosure"][0],
                     median_gap_s=np.median(gaps), longest_gap_s=gaps.max(), longest_gap_after=int(gaps.argmax()) + 1,
                     nan_frac=np.isnan(b["eyeclosure"]).mean(),
                     ec_min=np.nanmin(b["eyeclosure"]), ec_max=np.nanmax(b["eyeclosure"])))
npx_tab = pd.DataFrame(rows)
print(npx_tab.to_string())
print("identical eye traces between files:",
      [v for v in collections.Counter(hashlib.md5(b["eyeclosure"].tobytes()).hexdigest() for b in NPX.values()).items() if v[1] > 1])

  aff_w                  shape (2,) value [0.3 0.8]
  events-off             shape (48,)
  events-on              shape (48,)
  eyeclosure             shape (558508,)
  eyeclosure_by_puff     shape (384, 400)
  eyeclosure_by_trial    shape (7300, 48)
  eyesize                shape (558508,)
  eyesize_by_puff        shape (384, 400)
  eyesize_by_trial       shape (7300, 48)
  fps                    shape (1,) value [100]
  ket_puff               shape (1,) value [160]
  ket_trial              shape (1,) value [20]
  led                    shape (558508,)
  led_delay              shape (1,) value [5]
  npuffs_per_trial       shape (1,) value [8]
  p_post                 shape (1,) value [2]
  p_pre                  shape (1,) value [2]
  puff_aff               shape (384,)
  puff_ref               shape (384,)
  ref_w                  shape (2,) value [0.  0.3]
  t_eyeclosure           shape (400,)
  t_trial_eyeclosure     shape (7300,)
  trial_range_ket        shape (8,) value [22 23 24 25 26 27 28 29]
  trial_range_preket     shape (20,)
  w_post                 shape (1,) value [70]
  w_pre                  shape (1,) value [3]
              file      subject  n_onsets  n_puff_trials  ket_trial   pre    ket  fps       ref_w       aff_w  led_delay  by_trial_matches  t_trial_first  median_gap_s  longest_gap_s  longest_gap_after  nan_frac  ec_min  ec_max
0   sub-m11208-959   m11208-959        48             48         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        110.61         154.82                 20  0.002414     0.0     1.0
1   sub-m11208-962   m11208-962        59             59         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        124.48         434.09                 20  0.008956     0.0     1.0
2   sub-m11210-960   m11210-960        48             48         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        110.60         180.81                 20  0.000025     0.0     1.0
3  sub-m11219-1108  m11219-1108        58             58         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        125.65         165.93                 20  0.011032     0.0     1.0
4  sub-m11221-1119  m11221-1119        58             58         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        125.66         150.00                 20  0.008903     0.0     1.0
5  sub-m11223-1125  m11223-1125        58             58         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        125.66         154.01                 20  0.011217     0.0     1.0
6  sub-m11225-1131  m11225-1131        58             58         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        125.67         149.91                 20  0.008292     0.0     1.0
7  sub-m11227-1132  m11227-1132        58             58         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        125.66         298.04                 20  0.011044     0.0     1.0
8     sub-m304-876     m304-876        50             50         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        110.69         475.72                 20  0.000034     0.0     1.0
9      sub-m94-549      m94-549        79             40         20  0-19  22-29  100  (0.0, 0.3)  (0.3, 0.8)          5              True           -8.0        110.66         216.39                 20  0.000221     0.0     1.0
identical eye traces between files: []

conv = LISTS["convert_npx_mouse_airpuff_to_nwb"]["expids"]
conv_dates = collections.defaultdict(list)
for x in conv:
    conv_dates[x.split("_")[0]].append(f"{x.split('_')[1]}-{x.split('_')[2]}")
npx_sessions = sorted(dandi_io.list_assets(dandi_io.NPX_DANDISET, dandi_io.NPX_VERSION))
print(f"001326 sessions: {len(npx_sessions)}")
have = set(p.split("/")[0].replace("sub-", "") for p in NPX)
for p in npx_sessions:
    date = p.split("ses-")[1].split("_")[0]
    beh = conv_dates.get(date, [])
    print(f"   {p:45s} behaviour session(s) in the converter for this date: {beh}; "
          f"in 001347: {bool(beh) and all(x in have for x in beh)}")
print("converter sessions:", conv, "; all in 001347:", all(f"{x.split('_')[1]}-{x.split('_')[2]}" in have for x in conv))
conv_src = zf.read(me.SCRIPTS + "convert_npx_mouse_airpuff_to_nwb.py").decode()
print("converter lines that name sessions it does not convert:")
print(me.grep_lines(conv_src, ["# extra", "#   "]))

001326 sessions: 13
   sub-0/sub-0_ses-20230901_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11221-1119']; in 001347: True
   sub-1/sub-1_ses-20230218_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11208-959', 'm11210-960']; in 001347: True
   sub-1/sub-1_ses-20230226_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11208-962']; in 001347: True
   sub-2/sub-2_ses-20230902_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11223-1125']; in 001347: True
   sub-3/sub-3_ses-20230218_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11208-959', 'm11210-960']; in 001347: True
   sub-3/sub-3_ses-20230225_ecephys.nwb          behaviour session(s) in the converter for this date: []; in 001347: False
   sub-4/sub-4_ses-20230907_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11225-1131']; in 001347: True
   sub-5/sub-5_ses-20230706_ecephys.nwb          behaviour session(s) in the converter for this date: []; in 001347: False
   sub-5/sub-5_ses-20230822_ecephys.nwb          behaviour session(s) in the converter for this date: []; in 001347: False
   sub-6/sub-6_ses-20230823_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11219-1108']; in 001347: True
   sub-7/sub-7_ses-20220831_ecephys.nwb          behaviour session(s) in the converter for this date: ['m304-876']; in 001347: True
   sub-8/sub-8_ses-20220228_ecephys.nwb          behaviour session(s) in the converter for this date: ['m94-549']; in 001347: True
   sub-9/sub-9_ses-20230908_ecephys.nwb          behaviour session(s) in the converter for this date: ['m11227-1132']; in 001347: True
converter sessions: ['20220228_m94_549', '20220831_m304_876', '20230218_m11208_959', '20230226_m11208_962', '20230218_m11210_960', '20230823_m11219_1108', '20230901_m11221_1119', '20230902_m11223_1125', '20230907_m11225_1131', '20230908_m11227_1132'] ; all in 001347: True
converter lines that name sessions it does not convert:
   24  # extra = ['20220223_m95_544', '20230225_m11210_961',
   25  #          '20230706_m11217_1104', '20230822_m11217_1107',]

for p, b in NPX.items():
    fps = int(b["fps"][0]); ec = b["eyeclosure"]
    on = b["events-on"][:20] + 5 * fps
    d = []
    for shift, puffs in ((0.0, range(8)), (1.5, range(7)), (0.0, [0]), (1.5, [0])):
        vals = []
        for o in on:
            for k in puffs:
                q = int(o + (3 * k + shift) * fps)
                vals.append(np.nanmean(ec[q + int(0.1 * fps):q + int(0.3 * fps)]) - np.nanmean(ec[q - int(0.2 * fps):q]))
        d.append(np.nanmean(vals))
    print(f"{p.split('/')[0]:18s} response at puff times {d[0]:.3f}, at mid-points {d[1]:.3f}; "
          f"first puff only {d[2]:.3f}, 1.5 s after it {d[3]:.3f}")

sub-m11208-959     response at puff times 0.195, at mid-points -0.009; first puff only 0.208, 1.5 s after it -0.006
sub-m11208-962     response at puff times 0.124, at mid-points -0.011; first puff only 0.125, 1.5 s after it -0.008
sub-m11210-960     response at puff times 0.286, at mid-points -0.009; first puff only 0.331, 1.5 s after it 0.001
sub-m11219-1108    response at puff times 0.098, at mid-points -0.009; first puff only 0.258, 1.5 s after it -0.017
sub-m11221-1119    response at puff times -0.005, at mid-points -0.004; first puff only 0.069, 1.5 s after it -0.009
sub-m11223-1125    response at puff times 0.251, at mid-points -0.026; first puff only 0.250, 1.5 s after it -0.039
sub-m11225-1131    response at puff times 0.180, at mid-points -0.007; first puff only 0.211, 1.5 s after it -0.006
sub-m11227-1132    response at puff times 0.225, at mid-points -0.021; first puff only 0.269, 1.5 s after it -0.002
sub-m304-876       response at puff times 0.297, at mid-points -0.033; first puff only 0.518, 1.5 s after it -0.031
sub-m94-549        response at puff times 0.212, at mid-points -0.012; first puff only 0.300, 1.5 s after it -0.035

listed = {}
for script, v in LISTS.items():
    for cond, ids in v.get("experiments", {}).items():
        for e in ids:
            listed.setdefault(e, []).append((script, cond, v["ntrials"]))
in_data = set(IDS)
missing = sorted(set(listed) - in_data)
unlisted = sorted(in_data - set(listed))
print(f"(1) listed ids: {len(listed)}; in the dataset: {len(set(listed) & in_data)}; listed but missing: {missing}")
print(f"    sessions in the dataset in no list: {unlisted}")
short = [(e, s, c, n, int(meta.loc[e, 'n_trials'])) for e, uses in listed.items() for s, c, n in uses
         if e in in_data and meta.loc[e, "n_trials"] < n]
print(f"(2) sessions shorter than their script reads: {short}")
multi = {e: sorted({(s.replace('fig1_', ''), c) for s, c, _ in u}) for e, u in listed.items()
         if len({c for _, c, _ in u}) > 1}
print("    ids under more than one condition name:")
for e, u in sorted(multi.items()):
    print(f"      {e}: {u}")
twenty = sorted(e for e in IDS if meta.loc[e, "n_trials"] == 20)
twenty_scripts = {s for e in twenty for s, _, _ in listed[e]}
print(f"    the 20-trial sessions are listed only by: {sorted(twenty_scripts)}")

(1) listed ids: 142; in the dataset: 142; listed but missing: []
    sessions in the dataset in no list: []
(2) sessions shorter than their script reads: []
    ids under more than one condition name:
      282: [('dose', 'sal'), ('drugs', 'sal'), ('drugs_first_puff', 'sal'), ('persistence', 'sal'), ('sex', 'female')]
      283: [('dose', 'sal'), ('drugs', 'sal'), ('drugs_first_puff', 'sal'), ('persistence', 'sal'), ('sex', 'female')]
      284: [('dose', 'sal'), ('drugs', 'sal'), ('drugs_first_puff', 'sal'), ('persistence', 'sal'), ('sex', 'female')]
      285: [('dose', 'sal'), ('drugs', 'sal'), ('drugs_first_puff', 'sal'), ('persistence', 'sal'), ('sex', 'female')]
      286: [('dose', 'sal'), ('drugs', 'sal'), ('drugs_first_puff', 'sal'), ('persistence', 'sal'), ('sex', 'female')]
      1020: [('dose', 'ket35'), ('psi', '28 PSI')]
      1021: [('dose', 'ket35'), ('psi', '28 PSI')]
      1022: [('dose', 'ket35'), ('psi', '28 PSI')]
      1023: [('dose', 'ket35'), ('psi', '28 PSI')]
      1024: [('dose', 'ket35'), ('psi', '28 PSI')]
      1025: [('dose', 'ket35'), ('psi', '28 PSI')]
      1026: [('dose', 'ket35'), ('psi', '28 PSI')]
      1027: [('dose', 'ket35'), ('psi', '28 PSI')]
    the 20-trial sessions are listed only by: ['fig1_muscimol', 'fig1_psi']

E6 = LISTS["fig1_drugs"]["experiments"]
link_rows = []
for cond in ["sal", "ket50", "pcp20", "bup", "dzp", "kx"]:
    R = me.condition_values(DATA, E6[cond], 40, me.PRE_CODE, me.POST_CODE, 21)
    for i, e in enumerate(E6[cond]):
        link_rows.append(dict(cond=cond, session=e, pre_ratio=R["pre_ratio"][i], post_ratio=R["post_ratio"][i],
                              change=R["post_ratio"][i] - R["pre_ratio"][i],
                              post_early=R["post_early"][i], post_late=R["post_late"][i]))
link = pd.DataFrame(link_rows)
print(link.round(3).to_string(index=False))
med = link.groupby("cond")["change"].median()
a1 = bool((link[link.cond.isin(["ket50", "pcp20"])].change < 0).all())
a2 = bool(abs(med["sal"]) < min(abs(med["ket50"]), abs(med["pcp20"])))
kx = link[link.cond == "kx"]
b_ok = bool((kx.post_early < 0.1).all() and (kx.post_late < 0.1).all())
print(f"(3a) every ketamine and PCP session drops: {a1}; median change saline {med['sal']:.3f}, "
      f"ketamine {med['ket50']:.3f}, PCP {med['pcp20']:.3f}; saline smaller: {a2}")
print(f"(3b) K/X infusion early max {kx.post_early.max():.3f}, late max {kx.post_late.max():.3f}; both below 0.1: {b_ok}")

skip = np.array(LISTS["fig1_skippuff"]["SKIP_PUFFS"])
match_post, match_pre = [], []
for e in LISTS["fig1_skippuff"]["experiments"]["ket50"]:
    tr = me.prox_norm(me.trials(DATA[e], 40))
    J = np.stack([jump(tr, 3 * k) for k in range(8)], 1)
    inferred = J.argmin(1)
    match_pre.append((inferred[:20] == skip).mean()); match_post.append((inferred[20:] == skip).mean())
    print(f"     skip session {e}: SKIP_PUFFS matched in {match_pre[-1]:.0%} of preinfusion and {match_post[-1]:.0%} of infusion trials")
c_ok = bool(np.mean(match_post) >= 0.9)
print(f"(3c) skip-puff labels, infusion trials: {np.mean(match_post):.0%} matched; preinfusion {np.mean(match_pre):.0%}; holds: {c_ok}")
for e in [282, 316]:
    tr = me.prox_norm(me.trials(DATA[e], 40))
    J = np.stack([jump(tr, 3 * k) for k in range(8)], 1)
    print(f"     session {e} (also listed as sal_skip / kx_skip by the converter): smallest early response "
          f"over puffs in preinfusion trials, min over trials {J[:20].min(1).min():.3f}")
LINK_HOLDS = (not missing) and (not short) and a1 and a2 and b_ok and c_ok
print("link holds:", LINK_HOLDS)

 cond  session  pre_ratio  post_ratio  change  post_early  post_late
  sal      282      0.535       0.473  -0.062       0.795      0.376
  sal      283      0.594       0.632   0.038       0.796      0.503
  sal      284      0.617       0.622   0.006       0.624      0.388
  sal      285      0.621       0.849   0.228       0.918      0.780
  sal      286      0.920       0.947   0.027       0.761      0.720
ket50      275      0.678       0.212  -0.467       0.398      0.084
ket50      276      0.672      -0.130  -0.802       0.474     -0.062
ket50      277      0.666      -0.267  -0.933       0.250     -0.067
ket50      280      0.542      -0.035  -0.578       0.467     -0.017
ket50      281      0.734       0.039  -0.695       0.593      0.023
pcp20      287      0.702       0.244  -0.458       0.466      0.114
pcp20      288      0.704      -0.094  -0.798       0.393     -0.037
pcp20      289      0.616      -0.052  -0.668       0.410     -0.021
pcp20      290      0.591       0.060  -0.530       0.446      0.027
pcp20      291      0.935      -0.190  -1.125       0.427     -0.081
  bup      340      0.637       0.846   0.209       0.348      0.294
  bup      341      0.818       0.513  -0.306       0.743      0.381
  bup      342      0.839       0.784  -0.055       0.463      0.363
  bup      347      0.315       0.430   0.115       0.379      0.163
  bup      348      0.624       0.622  -0.002       0.420      0.261
  dzp      325      0.746       0.602  -0.144       0.520      0.313
  dzp      326      0.725       0.696  -0.029       0.592      0.412
  dzp      327      0.832       0.440  -0.392       0.407      0.179
  dzp      333      0.614       0.711   0.097       0.473      0.336
  dzp      334      0.822       0.631  -0.191       0.406      0.256
   kx      316      0.787       0.225  -0.562       0.041      0.009
   kx      338      0.563       0.818   0.255       0.062      0.051
   kx      339      0.586      -1.660  -2.246      -0.020      0.034
   kx      343      0.473       0.313  -0.160       0.046      0.014
   kx      344      0.728     -11.006 -11.735      -0.001      0.013
(3a) every ketamine and PCP session drops: True; median change saline 0.027, ketamine -0.695, PCP -0.668; saline smaller: True
(3b) K/X infusion early max 0.062, late max 0.051; both below 0.1: True
     skip session 411: SKIP_PUFFS matched in 100% of preinfusion and 100% of infusion trials
     skip session 412: SKIP_PUFFS matched in 100% of preinfusion and 100% of infusion trials
     skip session 413: SKIP_PUFFS matched in 95% of preinfusion and 100% of infusion trials
     skip session 420: SKIP_PUFFS matched in 80% of preinfusion and 100% of infusion trials
     skip session 421: SKIP_PUFFS matched in 100% of preinfusion and 100% of infusion trials
(3c) skip-puff labels, infusion trials: 100% matched; preinfusion 95%; holds: True
     session 282 (also listed as sal_skip / kx_skip by the converter): smallest early response over puffs in preinfusion trials, min over trials 0.421
     session 316 (also listed as sal_skip / kx_skip by the converter): smallest early response over puffs in preinfusion trials, min over trials 0.122
link holds: True

READINGS = {"a (code, 6-20 / 23-37)": (me.PRE_CODE, me.POST_CODE), "b (3-20 / 23-40)": (me.PRE_B, me.POST_B)}
FP = {}
for rname, (sp, sq) in READINGS.items():
    for scale in ("code", "methods"):
        FP[(rname, scale)] = {d: me.condition_values(DATA, E6[d], 40, sp, sq, 0, scale) for d in ("sal", "ket50", "pcp20")}
for rname in READINGS:
    R = FP[(rname, "code")]
    t = st.ttest_rel(R["sal"]["post_ratio"], R["ket50"]["post_ratio"])
    print(f"reading {rname}")
    for d in ("sal", "ket50"):
        print(f"   {d:6s} per mouse {R[d]['post_ratio'].round(3)}  mean {R[d]['post_ratio'].mean():.3f} "
              f"s.e.m. {st.sem(R[d]['post_ratio']):.3f}")
    print(f"   paired t = {t.statistic:.3f}, p = {t.pvalue:.4f}")
exp = expected("fig1k-mouse-late-eye-closure-ketamine")
R = FP[(list(READINGS)[0], "code")]
fig1k_vals = [R["sal"]["post_ratio"].mean(), R["ket50"]["post_ratio"].mean()]
for (q, v, tol), x in zip(exp, fig1k_vals):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.3f}; within: {within(x, v, tol)}")

reading a (code, 6-20 / 23-37)
   sal    per mouse [0.225 0.217 0.144 0.208 0.515]  mean 0.262 s.e.m. 0.065
   ket50  per mouse [ 0.122 -0.161 -0.191 -0.088 -0.019]  mean -0.068 s.e.m. 0.056
   paired t = 4.737, p = 0.0091
reading b (3-20 / 23-40)
   sal    per mouse [0.205 0.234 0.161 0.22  0.539]  mean 0.272 s.e.m. 0.068
   ket50  per mouse [ 0.087 -0.163 -0.172 -0.061 -0.011]  mean -0.064 s.e.m. 0.049
   paired t = 4.751, p = 0.0090
expected late/early first puff, saline: 0.26 +/- 0.06; measured 0.262; within: True
expected late/early first puff, ketamine: -0.07 +/- 0.06; measured -0.068; within: True

R = FP[(list(READINGS)[0], "code")]
t = np.arange(me.FRAMES_PER_TRIAL) / me.FPS - me.PRE_S
fig, ax = plt.subplots(1, 2, figsize=(8, 2.6), sharey=True)
rng_txt = []
for a, d in zip(ax, ("sal", "ket50")):
    pre = np.mean([me.mean_trace(DATA, e, 40, me.PRE_CODE) for e in E6[d]], 0)
    post = np.mean([me.mean_trace(DATA, e, 40, me.POST_CODE) for e in E6[d]], 0)
    a.plot(t, pre, lw=0.7, label="preinfusion"); a.plot(t, post, lw=0.7, label="infusion")
    a.set_xlim(-1, 30); a.set_title(d); a.set_xlabel("time from first puff (s)")
    rng_txt.append(f"{d}: pre {pre.min():.2f} to {pre.max():.2f}, infusion {post.min():.2f} to {post.max():.2f}")
ax[0].set_ylabel("eye closure (baseline subtracted)"); ax[0].legend(frameon=False, fontsize=7)
plt.tight_layout(); plt.show()
print("plotted ranges over the whole trial:", "; ".join(rng_txt))

plotted ranges over the whole trial: sal: pre -0.02 to 0.91, infusion -0.03 to 0.89; ket50: pre -0.02 to 0.86, infusion -0.04 to 0.64

exp = expected("fig1l-mouse-dissociatives-reduce-late")
for rname in READINGS:
    R = FP[(rname, "code")]
    tests = [st.ttest_rel(R[d]["pre_ratio"], R[d]["post_ratio"]) for d in ("sal", "ket50", "pcp20")]
    p = bh([x.pvalue for x in tests])
    print(f"reading {rname}")
    for d, x, pc in zip(("sal", "ket50", "pcp20"), tests, p):
        print(f"   {d:6s} pre {R[d]['pre_ratio'].mean():.3f}  infusion {R[d]['post_ratio'].mean():.3f}  "
              f"t = {x.statistic:.3f}  BH p = {pc:.4f}  mice lower on infusion: {(R[d]['post_ratio'] < R[d]['pre_ratio']).sum()}/5")
R = FP[(list(READINGS)[0], "code")]
fig1l_vals = [R["sal"]["pre_ratio"].mean(), R["sal"]["post_ratio"].mean(), R["ket50"]["pre_ratio"].mean(),
              R["ket50"]["post_ratio"].mean(), R["pcp20"]["pre_ratio"].mean(), R["pcp20"]["post_ratio"].mean()]
for (q, v, tol), x in zip(exp, fig1l_vals):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.3f}; within: {within(x, v, tol)}")

reading a (code, 6-20 / 23-37)
   sal    pre 0.286  infusion 0.262  t = 0.733  BH p = 0.5041  mice lower on infusion: 3/5
   ket50  pre 0.371  infusion -0.068  t = 9.525  BH p = 0.0020  mice lower on infusion: 5/5
   pcp20  pre 0.349  infusion -0.033  t = 3.409  BH p = 0.0406  mice lower on infusion: 5/5
reading b (3-20 / 23-40)
   sal    pre 0.293  infusion 0.272  t = 0.527  BH p = 0.6259  mice lower on infusion: 2/5
   ket50  pre 0.390  infusion -0.064  t = 8.440  BH p = 0.0032  mice lower on infusion: 5/5
   pcp20  pre 0.372  infusion -0.035  t = 3.532  BH p = 0.0363  mice lower on infusion: 5/5
expected pre-saline: 0.29 +/- 0.08; measured 0.286; within: True
expected saline: 0.26 +/- 0.06; measured 0.262; within: True
expected pre-ketamine: 0.37 +/- 0.06; measured 0.371; within: True
expected ketamine: -0.07 +/- 0.06; measured -0.068; within: True
expected pre-PCP: 0.35 +/- 0.1; measured 0.349; within: True
expected PCP: -0.03 +/- 0.05; measured -0.033; within: True

exp = expected("figs2bc-mouse-ketamine-early-present-late-abolished")
s2bc = {}
for (rname, scale), R in FP.items():
    te = st.ttest_rel(R["sal"]["post_early"], R["ket50"]["post_early"])
    tl = st.ttest_rel(R["sal"]["post_late"], R["ket50"]["post_late"])
    vals = [R["sal"]["post_early"].mean(), R["ket50"]["post_early"].mean(), R["sal"]["post_late"].mean(), R["ket50"]["post_late"].mean()]
    s2bc[(rname, scale)] = vals
    print(f"reading {rname}, scale {scale}: early saline {vals[0]:.3f} ketamine {vals[1]:.3f} (t = {te.statistic:.3f}, "
          f"p = {te.pvalue:.4f}); late saline {vals[2]:.3f} ketamine {vals[3]:.3f} (t = {tl.statistic:.3f}, p = {tl.pvalue:.4f})")
R = FP[(list(READINGS)[0], "code")]
print("late, saline infusion trials, per mouse:", R["sal"]["post_late"].round(3), "s.e.m.", round(st.sem(R["sal"]["post_late"]), 3),
      "median", round(np.median(R["sal"]["post_late"]), 3))
print("late, saline preinfusion trials of the same sessions, mean:", round(R["sal"]["pre_late"].mean(), 3))
print("ketamine early above zero in every mouse:", bool((R["ket50"]["post_early"] > 0).all()))
for (q, v, tol), x in zip(exp, s2bc[(list(READINGS)[0], "code")]):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.3f}; within: {within(x, v, tol)}")

reading a (code, 6-20 / 23-37), scale code: early saline 0.778 ketamine 0.451 (t = 6.854, p = 0.0024); late saline 0.203 ketamine -0.025 (t = 4.951, p = 0.0078)
reading a (code, 6-20 / 23-37), scale methods: early saline 0.722 ketamine 0.463 (t = 2.873, p = 0.0453); late saline 0.186 ketamine -0.020 (t = 5.155, p = 0.0067)
reading b (3-20 / 23-40), scale code: early saline 0.770 ketamine 0.430 (t = 6.065, p = 0.0037); late saline 0.208 ketamine -0.024 (t = 4.908, p = 0.0080)
reading b (3-20 / 23-40), scale methods: early saline 0.715 ketamine 0.441 (t = 2.993, p = 0.0402); late saline 0.191 ketamine -0.020 (t = 5.164, p = 0.0067)
late, saline infusion trials, per mouse: [0.191 0.17  0.089 0.19  0.374] s.e.m. 0.047 median 0.19
late, saline preinfusion trials of the same sessions, mean: 0.24
ketamine early above zero in every mouse: True
expected early, saline: 0.78 +/- 0.05; measured 0.778; within: True
expected early, ketamine: 0.45 +/- 0.06; measured 0.451; within: True
expected late, saline: 0.34 +/- 0.07; measured 0.203; within: False
expected late, ketamine: -0.02 +/- 0.02; measured -0.025; within: True

src_p = zf.read(me.SCRIPTS + "fig1_persistence.py").decode()
print("fig1_persistence.py, the lines that choose the trials and the start frame:")
print(me.grep_lines(src_p, ["ntrials =", "fps =", "aff_w =", "load_files_nwb(", "corr.reshape", "for exp in experiments",
                            "last_puff_ind =", "baseline_vals =", "b = m >", "first_inds = b", "no_first =", "first_inds[no_first]"]))
i0 = int(200 * (6 + 8 * 3 + 0.3))
print(f"the code's start frame: {i0} = {i0 / 200:.1f} s into the 81-s window, {i0 / 200 - 6:.1f} s after the first puff, "
      f"{i0 / 200 - 6 - 21:.1f} s after the eighth (at 21 s); the trials: all ntrials = 40 rows of each saline session\n")
sal = LISTS["fig1_persistence"]["experiments"]["sal"]
km = {}
for ntr, tlab in ((20, "20 preinfusion trials"), (40, "all 40 trials (code)")):
    for start, slab in ((30.3, "a: 3.3 s after 8th puff (code)"), (27.3, "b: 0.3 s after 8th puff")):
        T, Ev, per = [], [], []
        for e in sal:
            tt, ev = me.survival_times(me.trials(DATA[e], 40)[:ntr], start)
            T += list(tt); Ev += list(ev)
            per.append(KaplanMeierFitter().fit(tt, ev).median_survival_time_)
        k = KaplanMeierFitter().fit(T, Ev)
        km[(ntr, slab)] = k.median_survival_time_
        print(f"{tlab:22s} {slab:32s} n = {len(T)} trials, {int(np.sum(Ev))} reopen; pooled median {k.median_survival_time_:.2f} s; "
              f"per mouse {np.round(per, 2)}")
for th in (1.1, 1.2, 1.3):
    T, Ev = [], []
    for e in sal:
        tt, ev = me.survival_times(me.trials(DATA[e], 40)[:20], 27.3, thresh=th)
        T += list(tt); Ev += list(ev)
    print(f"threshold {th} x baseline, 20 trials, start b: median {KaplanMeierFitter().fit(T, Ev).median_survival_time_:.2f} s")
R = me.condition_values(DATA, sal, 40, np.arange(0, 20), np.arange(0, 20), 0)
R8 = me.condition_values(DATA, sal, 40, np.arange(0, 20), np.arange(0, 20), 21)
print("late eye closure, trials 1-20, per mouse: first puff", R["pre_late"].round(3), "eighth puff", R8["pre_late"].round(3))
print(f"mean first {R['pre_late'].mean():.3f}, eighth {R8['pre_late'].mean():.3f}; higher on eighth in "
      f"{(R8['pre_late'] > R['pre_late']).sum()}/5 mice")
exp = expected("fig1m-mouse-late-closure-accumulates-and-persists")
print(f"\nexpected median {exp[0][1]} +/- {exp[0][2]} s")
for ntr, tlab in ((20, "20 preinfusion (criterion)"), (40, "all 40 (code, not the criterion)")):
    for slab in ("a: 3.3 s after 8th puff (code)", "b: 0.3 s after 8th puff"):
        x = km[(ntr, slab)]
        print(f"   {tlab:33s} {slab:32s} median {x:.2f} s, minus 15.16: {x - exp[0][1]:+.2f}; within: {within(x, exp[0][1], exp[0][2])}")

fig1_persistence.py, the lines that choose the trials and the start frame:
   50      ntrials = 40
   51      fps = 200
   54      aff_w = [0.3, 0.8]
   59      corrs = etu.load_files_nwb(experiment_ids, datapath, ntrials)
   63          by_trial[key] = corr.reshape(ntrials, -1)
   68      for exp in experiments['sal']:
   70          last_puff_ind = int(fps * (6 + 8 * 3 + aff_w[0]))
   71          baseline_vals = m[:, 0:fps * 5].mean(1)
   73          b = m > 1.2 * bmat
   74          first_inds = b[:, last_puff_ind:].argmin(1)
   75          no_first = np.where(b[:, last_puff_ind:].min(1) == 1)[0]
   76          first_inds[no_first] = fps * 50
the code's start frame: 6060 = 30.3 s into the 81-s window, 24.3 s after the first puff, 3.3 s after the eighth (at 21 s); the trials: all ntrials = 40 rows of each saline session

20 preinfusion trials  a: 3.3 s after 8th puff (code)   n = 100 trials, 86 reopen; pooled median 13.05 s; per mouse [22.54  4.3   5.33 23.33 31.43]
20 preinfusion trials  b: 0.3 s after 8th puff          n = 100 trials, 86 reopen; pooled median 15.45 s; per mouse [25.54  7.3   5.34 26.33 34.43]
all 40 trials (code)   a: 3.3 s after 8th puff (code)   n = 200 trials, 166 reopen; pooled median 15.16 s; per mouse [22.54  5.44  5.33 24.24 33.37]
all 40 trials (code)   b: 0.3 s after 8th puff          n = 200 trials, 166 reopen; pooled median 18.16 s; per mouse [25.54  8.44  5.34 27.24 36.37]
threshold 1.1 x baseline, 20 trials, start b: median 14.92 s
threshold 1.2 x baseline, 20 trials, start b: median 15.45 s
threshold 1.3 x baseline, 20 trials, start b: median 16.84 s
late eye closure, trials 1-20, per mouse: first puff [0.321 0.208 0.072 0.249 0.369] eighth puff [0.449 0.463 0.379 0.578 0.716]
mean first 0.244, eighth 0.517; higher on eighth in 5/5 mice

expected median 15.16 +/- 1.52 s
   20 preinfusion (criterion)        a: 3.3 s after 8th puff (code)   median 13.05 s, minus 15.16: -2.11; within: False
   20 preinfusion (criterion)        b: 0.3 s after 8th puff          median 15.45 s, minus 15.16: +0.29; within: True
   all 40 (code, not the criterion)  a: 3.3 s after 8th puff (code)   median 15.16 s, minus 15.16: +0.00; within: True
   all 40 (code, not the criterion)  b: 0.3 s after 8th puff          median 18.16 s, minus 15.16: +3.00; within: False

P = LISTS["fig1_psi"]["experiments"]
psi = {}
for ps in ("8 PSI", "18 PSI", "28 PSI"):
    psi[ps] = np.array([[me.per_mouse(me.prox_norm(me.trials(DATA[e], 20))[me.TRIALS_20], 3 * k)["ratio"] for k in range(8)]
                        for e in P[ps]])
    print(f"{ps}: mean late/early by puff {psi[ps].mean(0).round(3)}")
tests, labels = [], []
for ps in ("8 PSI", "18 PSI", "28 PSI"):
    tests.append(st.ttest_rel(psi[ps][:, 0], psi[ps][:, 7])); labels.append(f"{ps} puff 1 vs 8")
for a, b in (("8 PSI", "18 PSI"), ("8 PSI", "28 PSI"), ("18 PSI", "28 PSI")):
    tests.append(st.ttest_rel(psi[a][:, 7], psi[b][:, 7])); labels.append(f"puff 8, {a} vs {b}")
pc = bh([x.pvalue for x in tests])
for l, x, p in zip(labels, tests, pc):
    print(f"   {l:26s} t = {x.statistic:7.3f}  BH p = {p:.4f}")
exp = expected("figs2g-mouse-scaling-with-puff-number-not-pressure")
s2g_vals = [psi["8 PSI"][:, 0].mean(), psi["8 PSI"][:, 7].mean(), psi["18 PSI"][:, 0].mean(),
            psi["18 PSI"][:, 7].mean(), psi["28 PSI"][:, 0].mean(), psi["28 PSI"][:, 7].mean()]
for (q, v, tol), x in zip(exp, s2g_vals):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.3f}; within: {within(x, v, tol)}")

8 PSI: mean late/early by puff [0.256 0.327 0.423 0.491 0.518 0.517 0.525 0.517]

18 PSI: mean late/early by puff [0.239 0.33  0.408 0.453 0.488 0.476 0.513 0.515]

28 PSI: mean late/early by puff [0.286 0.318 0.384 0.419 0.444 0.444 0.451 0.447]
   8 PSI puff 1 vs 8          t =  -4.834  BH p = 0.0038
   18 PSI puff 1 vs 8         t =  -9.243  BH p = 0.0001
   28 PSI puff 1 vs 8         t = -11.633  BH p = 0.0000
   puff 8, 8 PSI vs 18 PSI    t =   0.024  BH p = 0.9815
   puff 8, 8 PSI vs 28 PSI    t =   1.106  BH p = 0.4578
   puff 8, 18 PSI vs 28 PSI   t =   0.718  BH p = 0.5955
expected 8 PSI puff 1: 0.256 +/- 0.067; measured 0.256; within: True
expected 8 PSI puff 8: 0.517 +/- 0.058; measured 0.517; within: True
expected 18 PSI puff 1: 0.239 +/- 0.032; measured 0.239; within: True
expected 18 PSI puff 8: 0.515 +/- 0.037; measured 0.515; within: True
expected 28 PSI puff 1: 0.286 +/- 0.068; measured 0.286; within: True
expected 28 PSI puff 8: 0.447 +/- 0.061; measured 0.447; within: True

DOSE = LISTS["fig1_dose"]["experiments"]
od = ["sal", "ket50", "ket35", "ket25", "ket13", "ket6", "pcp20", "pcp10", "pcp5", "pcp1"]
should = {"ket50", "ket35", "pcp20", "pcp10"}
dose_out = {}
for rname, (sp, sq) in READINGS.items():
    R = {d: me.condition_values(DATA, DOSE[d], 40, sp, sq, 21) for d in od}
    tests = [st.ttest_rel(R[d]["pre_ratio"], R[d]["post_ratio"]) for d in od]
    pc = bh([x.pvalue for x in tests])
    agree = 0
    print(f"reading {rname}")
    for d, x, p in zip(od, tests, pc):
        sig = (p < 0.05) and R[d]["post_ratio"].mean() < R[d]["pre_ratio"].mean()
        ok = sig == (d in should); agree += ok
        print(f"   {d:6s} n={len(DOSE[d])} pre {R[d]['pre_ratio'].mean():.3f} infusion {R[d]['post_ratio'].mean():.3f} "
              f"t = {x.statistic:6.3f} BH p = {p:.4f} significant drop: {sig} as claimed: {ok}")
    dose_out[rname] = agree
    print(f"   conditions agreeing with the claim: {agree}/10")

reading a (code, 6-20 / 23-37)
   sal    n=5 pre 0.657 infusion 0.705 t = -0.980 BH p = 0.4252 significant drop: False as claimed: True
   ket50  n=5 pre 0.659 infusion -0.036 t =  8.488 BH p = 0.0106 significant drop: True as claimed: True
   ket35  n=8 pre 0.447 infusion 0.169 t =  3.843 BH p = 0.0159 significant drop: True as claimed: True
   ket25  n=5 pre 0.614 infusion 0.439 t =  1.232 BH p = 0.3566 significant drop: False as claimed: True
   ket13  n=5 pre 0.573 infusion 1.000 t = -1.270 BH p = 0.3566 significant drop: False as claimed: True
   ket6   n=5 pre 0.702 infusion 0.678 t =  0.329 BH p = 0.7585 significant drop: False as claimed: True
   pcp20  n=5 pre 0.709 infusion -0.006 t =  6.080 BH p = 0.0159 significant drop: True as claimed: True
   pcp10  n=5 pre 0.719 infusion 0.205 t =  5.404 BH p = 0.0159 significant drop: True as claimed: True
   pcp5   n=5 pre 0.644 infusion 0.505 t =  2.682 BH p = 0.1103 significant drop: False as claimed: True
   pcp1   n=5 pre 0.610 infusion 0.395 t =  2.064 BH p = 0.1799 significant drop: False as claimed: True
   conditions agreeing with the claim: 10/10
reading b (3-20 / 23-40)
   sal    n=5 pre 0.651 infusion 0.696 t = -0.807 BH p = 0.5164 significant drop: False as claimed: True
   ket50  n=5 pre 0.655 infusion -0.025 t =  8.807 BH p = 0.0092 significant drop: True as claimed: True
   ket35  n=8 pre 0.454 infusion 0.175 t =  3.735 BH p = 0.0192 significant drop: True as claimed: True
   ket25  n=5 pre 0.598 infusion 0.445 t =  1.049 BH p = 0.4416 significant drop: False as claimed: True
   ket13  n=5 pre 0.570 infusion 1.122 t = -1.234 BH p = 0.4066 significant drop: False as claimed: True
   ket6   n=5 pre 0.682 infusion 0.665 t =  0.261 BH p = 0.8068 significant drop: False as claimed: True
   pcp20  n=5 pre 0.698 infusion -0.007 t =  5.992 BH p = 0.0192 significant drop: True as claimed: True
   pcp10  n=5 pre 0.712 infusion 0.219 t =  4.964 BH p = 0.0192 significant drop: True as claimed: True
   pcp5   n=5 pre 0.641 infusion 0.497 t =  2.530 BH p = 0.1253 significant drop: False as claimed: True
   pcp1   n=5 pre 0.606 infusion 0.399 t =  2.390 BH p = 0.1253 significant drop: False as claimed: True
   conditions agreeing with the claim: 10/10

o5 = ["sal", "ket50", "pcp20", "bup", "dzp"]
EIGHTH = {}
for rname, (sp, sq) in READINGS.items():
    EIGHTH[rname] = {d: me.condition_values(DATA, E6[d], 40, sp, sq, 21) for d in E6}
    R = EIGHTH[rname]
    tests = [st.ttest_rel(R[d]["pre_ratio"], R[d]["post_ratio"]) for d in o5]
    pc = bh([x.pvalue for x in tests])
    print(f"reading {rname}")
    for d, x, p in zip(o5, tests, pc):
        print(f"   {d:6s} pre {R[d]['pre_ratio'].mean():.3f} infusion {R[d]['post_ratio'].mean():.3f} "
              f"t = {x.statistic:7.4f} BH p = {p:.4f} mean change {np.mean(R[d]['post_ratio'] - R[d]['pre_ratio']):+.3f}")

reading a (code, 6-20 / 23-37)
   sal    pre 0.657 infusion 0.705 t = -0.9797 BH p = 0.4784 mean change +0.047
   ket50  pre 0.659 infusion -0.036 t =  8.4882 BH p = 0.0053 mean change -0.695
   pcp20  pre 0.709 infusion -0.006 t =  6.0804 BH p = 0.0092 mean change -0.716
   bup    pre 0.647 infusion 0.639 t =  0.0878 BH p = 0.9342 mean change -0.008
   dzp    pre 0.748 infusion 0.616 t =  1.6120 BH p = 0.3038 mean change -0.132
reading b (3-20 / 23-40)
   sal    pre 0.651 infusion 0.696 t = -0.8073 BH p = 0.5810 mean change +0.045
   ket50  pre 0.655 infusion -0.025 t =  8.8073 BH p = 0.0046 mean change -0.680
   pcp20  pre 0.698 infusion -0.007 t =  5.9920 BH p = 0.0098 mean change -0.705
   bup    pre 0.647 infusion 0.644 t =  0.0372 BH p = 0.9721 mean change -0.003
   dzp    pre 0.749 infusion 0.593 t =  2.1353 BH p = 0.1660 mean change -0.156

o6 = ["sal", "ket50", "pcp20", "bup", "dzp", "kx"]
exp = expected("fig1p-mouse-anesthesia-abolishes-both")
fig1p_vals = {}
for rname, (sp, sq) in READINGS.items():
    for scale in ("code", "methods"):
        R = {d: me.condition_values(DATA, E6[d], 40, sp, sq, 21, scale) for d in o6}
        out = []
        for m in ("early", "late"):
            tests = [st.ttest_rel(R[d]["pre_" + m], R[d]["post_" + m]) for d in o6]
            pc = bh([x.pvalue for x in tests])
            out.append(f"{m}: K/X pre {R['kx']['pre_' + m].mean():.3f} infusion {R['kx']['post_' + m].mean():.3f} "
                       f"t = {tests[-1].statistic:.4f} BH p = {pc[-1]:.4f}")
        fig1p_vals[(rname, scale)] = [R["kx"]["post_early"].mean(), R["kx"]["post_late"].mean()]
        print(f"reading {rname}, scale {scale}: " + "; ".join(out))
for (q, v, tol), x in zip(exp, fig1p_vals[(list(READINGS)[0], "code")]):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.3f}; within: {within(x, v, tol)}")

reading a (code, 6-20 / 23-37), scale code: early: K/X pre 0.773 infusion 0.025 t = 8.6017 BH p = 0.0060; late: K/X pre 0.470 infusion 0.024 t = 12.0602 BH p = 0.0016
reading a (code, 6-20 / 23-37), scale methods: early: K/X pre 0.777 infusion 0.023 t = 11.9189 BH p = 0.0017; late: K/X pre 0.479 infusion 0.025 t = 12.9131 BH p = 0.0006
reading b (3-20 / 23-40), scale code: early: K/X pre 0.783 infusion 0.025 t = 8.8342 BH p = 0.0030; late: K/X pre 0.466 infusion 0.026 t = 11.1774 BH p = 0.0022

reading b (3-20 / 23-40), scale methods: early: K/X pre 0.787 infusion 0.023 t = 12.6979 BH p = 0.0013; late: K/X pre 0.475 infusion 0.027 t = 11.7544 BH p = 0.0009
expected early eye closure, K/X: 0.03 +/- 0.02; measured 0.025; within: True
expected late eye closure, K/X: 0.02 +/- 0.01; measured 0.024; within: True

exp = expected("fig1no-mouse-dissociatives-no-accumulation")
for rname, (sp, sq) in READINGS.items():
    for t_end, wlab in ((23.2, "21.2-23.2 s (code)"), (23.0, "21.2-23.0 s (legend)")):
        bp = {d: np.array([me.fit_decay(me.mean_trace(DATA, e, 40, sp), 21.2, t_end) for e in E6[d]]) for d in o5}
        bq = {d: np.array([me.fit_decay(me.mean_trace(DATA, e, 40, sq), 21.2, t_end) for e in E6[d]]) for d in o5}
        fails = sum(int(np.isnan(bp[d]).sum() + np.isnan(bq[d]).sum()) for d in o5)
        tests = [st.ttest_rel(bp[d], bq[d]) for d in o5]
        pc = bh([x.pvalue for x in tests])
        print(f"reading {rname}, window {wlab}: failed fits {fails}")
        for d, x, p in zip(o5, tests, pc):
            print(f"   {d:6s} pre {bp[d].mean():6.2f} infusion {bq[d].mean():6.2f} t = {x.statistic:7.4f} BH p = {p:.4f}")
        if rname == list(READINGS)[0] and t_end == 23.2:
            fig1no_vals = [bp["ket50"].mean(), bq["ket50"].mean(), bp["pcp20"].mean(), bq["pcp20"].mean()]
            print("   per mouse, ketamine infusion:", bq["ket50"].round(2), " PCP infusion:", bq["pcp20"].round(2))
for (q, v, tol), x in zip(exp, fig1no_vals):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.2f}; within: {within(x, v, tol)}")
R = EIGHTH[list(READINGS)[0]]
R1 = {d: me.condition_values(DATA, E6[d], 40, me.PRE_CODE, me.POST_CODE, 0) for d in ("sal", "ket50", "pcp20")}
for d in ("sal", "ket50", "pcp20"):
    print(f"{d}: infusion late, first puff {R1[d]['post_late'].mean():.3f}, eighth puff {R[d]['post_late'].mean():.3f} "
          f"(eighth above first in {(R[d]['post_late'] > R1[d]['post_late']).sum()}/5 mice, paired t = "
          f"{st.ttest_rel(R[d]['post_late'], R1[d]['post_late']).statistic:.2f}, p = {st.ttest_rel(R[d]['post_late'], R1[d]['post_late']).pvalue:.3f}); "
          f"infusion early on first puff {R1[d]['post_early'].mean():.3f}")
# Effect sizes: infusion late eye closure, eighth minus first puff, mean and paired 95% CI per condition;
# and each drug's rise minus the saline rise (Welch 95% CI; the sessions are treated as independent,
# an assumption, since the files carry no mouse label).
rise = {}
for d in ("sal", "ket50", "pcp20"):
    x = st.ttest_rel(R[d]["post_late"], R1[d]["post_late"]); ci = x.confidence_interval()
    rise[d] = R[d]["post_late"] - R1[d]["post_late"]
    print(f"{d}: infusion late, eighth minus first puff {rise[d].mean():+.3f}, 95% CI [{ci.low:+.3f}, {ci.high:+.3f}] "
          f"(paired t = {x.statistic:.2f}, p = {x.pvalue:.3f}, n = {len(rise[d])})")
for d in ("ket50", "pcp20"):
    x = st.ttest_ind(rise[d], rise["sal"], equal_var=False); ci = x.confidence_interval()
    print(f"{d} rise minus saline rise {rise[d].mean() - rise['sal'].mean():+.3f}, 95% CI [{ci.low:+.3f}, {ci.high:+.3f}] "
          f"(Welch t = {x.statistic:.2f}, p = {x.pvalue:.4f})")

reading a (code, 6-20 / 23-37), window 21.2-23.2 s (code): failed fits 0
   sal    pre   3.33 infusion   3.96 t = -0.8204 BH p = 0.7634
   ket50  pre   3.53 infusion  11.89 t = -4.5198 BH p = 0.0266
   pcp20  pre   3.04 infusion  13.20 t = -5.3534 BH p = 0.0266
   bup    pre   4.15 infusion   3.60 t =  0.4065 BH p = 0.8243
   dzp    pre   3.29 infusion   3.03 t =  0.2370 BH p = 0.8243
   per mouse, ketamine infusion: [ 7.42 14.77 17.    9.78 10.47]  PCP infusion: [ 7.33 14.39 13.26 14.08 16.93]

reading a (code, 6-20 / 23-37), window 21.2-23.0 s (legend): failed fits 0
   sal    pre   3.36 infusion   4.04 t = -0.8772 BH p = 0.7165
   ket50  pre   3.64 infusion  11.82 t = -4.5097 BH p = 0.0269
   pcp20  pre   3.10 infusion  13.15 t = -5.5518 BH p = 0.0258
   bup    pre   4.16 infusion   3.72 t =  0.3332 BH p = 0.8668
   dzp    pre   3.29 infusion   3.09 t =  0.1788 BH p = 0.8668

reading b (3-20 / 23-40), window 21.2-23.2 s (code): failed fits 0
   sal    pre   3.43 infusion   3.74 t = -0.4731 BH p = 0.8260
   ket50  pre   3.35 infusion  11.77 t = -4.9200 BH p = 0.0198
   pcp20  pre   3.23 infusion  13.29 t = -5.5511 BH p = 0.0198
   bup    pre   4.49 infusion   3.45 t =  0.6955 BH p = 0.8260
   dzp    pre   3.23 infusion   3.09 t =  0.1465 BH p = 0.8906

reading b (3-20 / 23-40), window 21.2-23.0 s (legend): failed fits 0
   sal    pre   3.45 infusion   3.85 t = -0.5918 BH p = 0.7323
   ket50  pre   3.46 infusion  11.71 t = -4.9531 BH p = 0.0194
   pcp20  pre   3.34 infusion  13.25 t = -5.8062 BH p = 0.0194
   bup    pre   4.52 infusion   3.61 t =  0.6169 BH p = 0.7323
   dzp    pre   3.25 infusion   3.13 t =  0.1216 BH p = 0.9091
expected decay rate pre-ketamine: 3.53 +/- 0.28; measured 3.53; within: True
expected decay rate ketamine: 11.89 +/- 1.74; measured 11.89; within: True
expected decay rate pre-PCP: 3.04 +/- 0.46; measured 3.04; within: True
expected decay rate PCP: 13.2 +/- 1.59; measured 13.20; within: True
sal: infusion late, first puff 0.203, eighth puff 0.553 (eighth above first in 5/5 mice, paired t = 5.30, p = 0.006); infusion early on first puff 0.778
ket50: infusion late, first puff -0.025, eighth puff -0.007 (eighth above first in 4/5 mice, paired t = 1.65, p = 0.175); infusion early on first puff 0.451
pcp20: infusion late, first puff -0.010, eighth puff 0.000 (eighth above first in 4/5 mice, paired t = 0.93, p = 0.404); infusion early on first puff 0.427
sal: infusion late, eighth minus first puff +0.350, 95% CI [+0.167, +0.534] (paired t = 5.30, p = 0.006, n = 5)
ket50: infusion late, eighth minus first puff +0.017, 95% CI [-0.012, +0.046] (paired t = 1.65, p = 0.175, n = 5)
pcp20: infusion late, eighth minus first puff +0.010, 95% CI [-0.019, +0.039] (paired t = 0.93, p = 0.404, n = 5)
ket50 rise minus saline rise -0.333, 95% CI [-0.516, -0.151] (Welch t = -4.98, p = 0.0067)
pcp20 rise minus saline rise -0.341, 95% CI [-0.523, -0.158] (Welch t = -5.09, p = 0.0062)

DAYS = LISTS["fig1_days"]["experiments"]
days = ["day1", "day2", "day3", "day4", "day5"]
for rname, (sp, sq) in READINGS.items():
    R = {d: me.condition_values(DATA, DAYS[d], 40, sp, sq, 21) for d in days}
    n_sig = 0
    print(f"reading {rname}")
    for m, lab in (("ratio", "normalised"), ("early", "reflexive"), ("late", "affective")):
        tests = [st.ttest_rel(R["day1"]["pre_" + m], R[d]["pre_" + m]) for d in days[1:]]
        pc = bh([x.pvalue for x in tests]); n_sig += int((pc < 0.05).sum())
        print(f"   {lab:10s} means {[round(R[d]['pre_' + m].mean(), 3) for d in days]} |t| {[round(abs(x.statistic), 2) for x in tests]} BH p {pc.round(4)}")
    print(f"   tests with corrected p < 0.05: {n_sig}/12")

reading a (code, 6-20 / 23-37)
   normalised means [np.float64(0.422), np.float64(0.466), np.float64(0.305), np.float64(0.502), np.float64(0.478)] |t| [np.float64(0.74), np.float64(1.38), np.float64(1.24), np.float64(0.48)] BH p [0.657 0.564 0.564 0.657]
   reflexive  means [np.float64(0.392), np.float64(0.359), np.float64(0.41), np.float64(0.677), np.float64(0.428)] |t| [np.float64(0.91), np.float64(0.35), np.float64(3.31), np.float64(0.63)] BH p [0.743 0.743 0.119 0.743]
   affective  means [np.float64(0.166), np.float64(0.166), np.float64(0.124), np.float64(0.341), np.float64(0.208)] |t| [np.float64(0.02), np.float64(1.13), np.float64(3.19), np.float64(0.63)] BH p [0.988 0.642 0.132 0.748]
   tests with corrected p < 0.05: 0/12
reading b (3-20 / 23-40)
   normalised means [np.float64(0.427), np.float64(0.466), np.float64(0.307), np.float64(0.484), np.float64(0.48)] |t| [np.float64(0.65), np.float64(1.37), np.float64(0.86), np.float64(0.47)] BH p [0.663 0.663 0.663 0.663]
   reflexive  means [np.float64(0.4), np.float64(0.368), np.float64(0.416), np.float64(0.681), np.float64(0.442)] |t| [np.float64(0.85), np.float64(0.33), np.float64(3.32), np.float64(0.72)] BH p [0.679 0.758 0.118 0.679]
   affective  means [np.float64(0.171), np.float64(0.171), np.float64(0.128), np.float64(0.333), np.float64(0.214)] |t| [np.float64(0.02), np.float64(1.12), np.float64(2.82), np.float64(0.66)] BH p [0.985 0.65  0.191 0.724]
   tests with corrected p < 0.05: 0/12

SEX = LISTS["fig1_sex"]["experiments"]
for rname, (sp, sq) in READINGS.items():
    R = {s: me.condition_values(DATA, SEX[s], 40, sp, sq, 21) for s in ("male", "female")}
    tests = [st.ttest_ind(R["male"]["pre_ratio"], R["female"]["pre_ratio"]),
             st.ttest_ind(R["male"]["post_ratio"], R["female"]["post_ratio"]),
             st.ttest_rel(R["male"]["pre_ratio"], R["male"]["post_ratio"]),
             st.ttest_rel(R["female"]["pre_ratio"], R["female"]["post_ratio"])]
    pc = bh([x.pvalue for x in tests])
    print(f"reading {rname}: means male pre {R['male']['pre_ratio'].mean():.3f} saline {R['male']['post_ratio'].mean():.3f}; "
          f"female pre {R['female']['pre_ratio'].mean():.3f} saline {R['female']['post_ratio'].mean():.3f}")
    for l, x, p in zip(["pre male vs female", "saline male vs female", "male pre vs saline", "female pre vs saline"], tests, pc):
        print(f"   {l:24s} t = {x.statistic:6.3f} BH p = {p:.4f}")

reading a (code, 6-20 / 23-37): means male pre 0.445 saline 0.567; female pre 0.657 saline 0.705
   pre male vs female       t = -2.079 BH p = 0.2851
   saline male vs female    t = -0.983 BH p = 0.3827
   male pre vs saline       t = -1.080 BH p = 0.3827
   female pre vs saline     t = -0.980 BH p = 0.3827
reading b (3-20 / 23-40): means male pre 0.463 saline 0.557; female pre 0.651 saline 0.696
   pre male vs female       t = -1.945 BH p = 0.3506
   saline male vs female    t = -1.084 BH p = 0.4648
   male pre vs saline       t = -0.927 BH p = 0.4648
   female pre vs saline     t = -0.807 BH p = 0.4648

SK = LISTS["fig1_skippuff"]["experiments"]["ket50"]
skip_ok = []
for k in (5, 6):
    rows = []
    for e in SK:
        m = me.prox_norm(me.trials(DATA[e], 40))[np.where(skip == k)[0]].mean(0)
        ts = 3 * k
        rows.append(dict(mouse=e, n_trials=int((skip == k).sum()),
                         rise_at_skip=m[me.window(ts, (0.1, 0.2))].mean() - m[me.window(ts, (-0.1, 0.0))].mean(),
                         level_at_skip=m[me.window(ts, (0.1, 0.2))].mean(),
                         late_at_skip=m[me.window(ts, (0.3, 0.8))].mean(),
                         rise_prev_puff=m[me.window(3 * (k - 1), (0.1, 0.2))].mean() - m[me.window(3 * (k - 1), (-0.1, 0.0))].mean()))
    df = pd.DataFrame(rows)
    print(f"skipped puff {k + 1}:")
    print(df.round(3).to_string(index=False))
    rise, sem = df.rise_at_skip.mean(), st.sem(df.rise_at_skip)
    ok = (rise <= sem) and (df.level_at_skip.mean() > 0)
    skip_ok.append(ok)
    print(f"   mean rise {rise:.4f}, s.e.m. {sem:.4f}; mean level {df.level_at_skip.mean():.3f}; "
          f"rise at the previous (delivered) puff {df.rise_prev_puff.mean():.3f}; criterion met: {ok}")
print("both skipped positions meet the criterion:", all(skip_ok))

skipped puff 6:
 mouse  n_trials  rise_at_skip  level_at_skip  late_at_skip  rise_prev_puff
   411        10        -0.016          0.374         0.379           0.570
   412        10        -0.005          0.301         0.299           0.336
   413        10         0.008          0.057         0.054           0.223
   420        10         0.028          0.116         0.091           0.160
   421        10        -0.011          0.171         0.159           0.558
   mean rise 0.0008, s.e.m. 0.0079; mean level 0.204; rise at the previous (delivered) puff 0.369; criterion met: True
skipped puff 7:
 mouse  n_trials  rise_at_skip  level_at_skip  late_at_skip  rise_prev_puff
   411        10        -0.004          0.297         0.297           0.516
   412        10        -0.004          0.253         0.246           0.397
   413        10         0.013          0.047         0.043           0.222
   420        10         0.003          0.106         0.099           0.177
   421        10        -0.005          0.207         0.212           0.555
   mean rise 0.0004, s.e.m. 0.0035; mean level 0.182; rise at the previous (delivered) puff 0.373; criterion met: True
both skipped positions meet the criterion: True

assets = dandi_io.list_assets()
print("assets with 'lick' or 'reward' in the path:", [p for p in assets if "lick" in p.lower() or "reward" in p.lower()])
print("acquisitions in a sub-mouse file:", ["eye_closure"])
print("acquisitions in the Neuropixels-rig files:", sorted(k for k in b0 if k != "subject_id"))
print("code record files mentioning lick:", [n for n in zf.namelist() if n.endswith(".py") and b"lick" in zf.read(n).lower()])
print("sub-mouse sessions not used by any listed condition:", unlisted)

assets with 'lick' or 'reward' in the path: []
acquisitions in a sub-mouse file: ['eye_closure']
acquisitions in the Neuropixels-rig files: ['aff_w', 'events-off', 'events-on', 'eyeclosure', 'eyeclosure_by_puff', 'eyeclosure_by_trial', 'eyesize', 'eyesize_by_puff', 'eyesize_by_trial', 'fps', 'ket_puff', 'ket_trial', 'led', 'led_delay', 'npuffs_per_trial', 'p_post', 'p_pre', 'puff_aff', 'puff_ref', 'ref_w', 't_eyeclosure', 't_trial_eyeclosure', 'trial_range_ket', 'trial_range_preket', 'w_post', 'w_pre']
code record files mentioning lick: []
sub-mouse sessions not used by any listed condition: []

MU = LISTS["fig1_muscimol"]["experiments"]
R = {d: me.condition_values(DATA, MU[d], 20, me.TRIALS_20, me.TRIALS_20, 21) for d in MU}
tests = [st.ttest_rel(R["sal"]["pre_ratio"], R["intact"]["pre_ratio"]),
         st.ttest_rel(R["sal"]["pre_ratio"], R["vpm_musc"]["pre_ratio"]),
         st.ttest_rel(R["intact"]["pre_ratio"], R["vpm_musc"]["pre_ratio"])]
pc = bh([x.pvalue for x in tests])
for l, x, p in zip(["saline vs no injection", "saline vs muscimol", "no injection vs muscimol"], tests, pc):
    print(f"late/early {l:26s} t = {x.statistic:6.3f} BH p = {p:.4f}")
print({d: (R[d]["pre_ratio"].round(3), round(R[d]["pre_ratio"].mean(), 3)) for d in R})
F = {d: me.condition_values(DATA, MU[d], 20, me.TRIALS_20, me.TRIALS_20, 0, w_early=(0.0, 0.04)) for d in MU}
for test in (st.ttest_ind, st.ttest_rel):
    tt = [test(F["sal"]["pre_early"], F["intact"]["pre_early"]), test(F["sal"]["pre_early"], F["vpm_musc"]["pre_early"]),
          test(F["vpm_musc"]["pre_early"], F["intact"]["pre_early"])]
    pp = bh([x.pvalue for x in tt])
    print(f"fast early 0-40 ms ({test.__name__}): means {[round(F[d]['pre_early'].mean(), 4) for d in ('intact', 'sal', 'vpm_musc')]} "
          f"saline vs muscimol t = {tt[1].statistic:.3f} BH p = {pp[1]:.4f}; all BH p {pp.round(4)}")
exp = expected("figs6-mouse-vpm-inhibition-reduces-affective")
s6_vals = [R["sal"]["pre_ratio"].mean(), R["vpm_musc"]["pre_ratio"].mean()]
for (q, v, tol), x in zip(exp, s6_vals):
    print(f"expected {q}: {v} +/- {tol}; measured {x:.3f}; within: {within(x, v, tol)}")

late/early saline vs no injection     t = -0.455 BH p = 0.6727

late/early saline vs muscimol         t =  3.845 BH p = 0.0289
late/early no injection vs muscimol   t =  3.789 BH p = 0.0289
{'intact': (array([0.613, 0.749, 0.542, 0.657, 0.603]), np.float64(0.633)), 'sal': (array([0.724, 0.586, 0.618, 0.434, 0.649]), np.float64(0.602)), 'vpm_musc': (array([ 0.339, -0.284, -0.078,  0.303,  0.247]), np.float64(0.105))}
fast early 0-40 ms (ttest_ind): means [np.float64(0.0801), np.float64(0.0755), np.float64(0.0887)] saline vs muscimol t = -0.367 BH p = 0.8898; all BH p [0.89 0.89 0.89]
fast early 0-40 ms (ttest_rel): means [np.float64(0.0801), np.float64(0.0755), np.float64(0.0887)] saline vs muscimol t = -0.379 BH p = 0.9013; all BH p [0.901 0.901 0.901]
expected late/early, saline injection: 0.6 +/- 0.05; measured 0.602; within: True
expected late/early, muscimol: 0.11 +/- 0.12; measured 0.105; within: True

src_an = me.fetch_analysis_py()
ROOT = me.SCRIPTS.split("mouse_eyepuff/")[0]
src_npx = zf.read(me.SCRIPTS + "fig1_npx_behavior.py").decode()
src_etu = zf.read(ROOT + "ketmouse/eye_trace_utils.py").decode()
print("fig1_npx_behavior.py (the script the criterion names), main():")
print(me.show_lines(src_npx, *me.def_lines(src_npx, "main")))
print("\nketmouse/eye_trace_utils.py, plot_npx_eyeclosure, which main() calls (whole function):")
print(me.show_lines(src_etu, *me.def_lines(src_etu, "plot_npx_eyeclosure")))
print("\nanalysis.py (Zenodo 14876694, md5 checked), make_puffseries_behavior_plot, the lines that cut, baseline and normalise:")
print(me.def_grep(src_an, "make_puffseries_behavior_plot",
                  ["trial_events =", "pertrial.append", "peranimal.append", "preket_behavior.append",
                   "postket_behavior.append", "ste ="]))
MAXKEYS = (".max()", "nanmax", "np.max(", "amax")
for name, s in (("fig1_npx_behavior.py main + plot_npx_eyeclosure",
                 me.def_source(src_npx, "main") + me.def_source(src_etu, "plot_npx_eyeclosure")),
                ("analysis.py make_puffseries_behavior_plot", me.def_source(src_an, "make_puffseries_behavior_plot"))):
    print(f"{name}: divides by a maximum: {any(k in s for k in MAXKEYS)}")

fig1_npx_behavior.py (the script the criterion names), main():
   18  def main():
   19  
   20      for expid in expids:
   21          datadir, plotdir = io.get_dirs(pt=expid)
   22          savepath = os.path.join(datadir, f'sub-{expid}_ses-puff_behavior.nwb')
   23          nwb_loaded = etu.npx_load_nwb(expid, savepath)
   24  
   25          etu.plot_npx_eyeclosure(expid, nwb_loaded, plotdir)
   26          plt.show()

ketmouse/eye_trace_utils.py, plot_npx_eyeclosure, which main() calls (whole function):
  130  def plot_npx_eyeclosure(expid, data, savedir=None):
  131      by_trial_closure = data['eyeclosure_by_trial']
  132      trial_range_preket = data['trial_range_preket']
  133      trial_range_ket = data['trial_range_ket']
  134      ket_trial = data['ket_trial']
  135      w_pre = data['w_pre']
  136      w_post = data['w_post']
  137      led_delay = data['led_delay']
  138      t_trial = data['t_trial_eyeclosure']
  139  
  140      fontsize = 8
  141      preket_closure = np.nanmedian(by_trial_closure[trial_range_preket, :],
  142                                    axis=0)
  143      ket_closure = np.nanmedian(by_trial_closure[trial_range_ket, :],
  144                                 axis=0)
  145      plt.figure(figsize=(3.4, 1.75))
  146      plt.subplot(1, 2, 1)
  147      plt.imshow(by_trial_closure[:, :], aspect='auto',
  148                 interpolation='nearest',
  149                 extent=[-w_pre - led_delay, w_post - led_delay,
  150                         by_trial_closure.shape[0], 0],
  151                 clim=[0, 1])
  152  
  153      plt.axvline((0), linestyle='--', color='k', linewidth=0.5)
  154      plt.axhline(ket_trial, linestyle='--', color='k', linewidth=0.5)
  155      # plt.title(f'{expid}, eye closure')
  156      plt.ylabel('Trial #', fontsize=fontsize)
  157      plt.xlabel('Time (s) ', fontsize=fontsize)
  158  
  159      plt.subplot(1, 2, 2)
  160      plt.plot(t_trial, preket_closure, '#6AADD6', linewidth=1)
  161      plt.plot(t_trial, ket_closure, '#CD1D58', linewidth=1)
  162      plt.axvline((0), linestyle='--', color='k', linewidth=0.5)
  163      plt.ylabel('Eye closure', fontsize=fontsize)
  164      plt.xlabel('Time (s)', fontsize=fontsize)
  165      plt.xticks([0, 25, 50], fontsize=fontsize)
  166      # plt.legend(['Pre', 'Ket.'], frameon=False, handlelength=0.5, handletextpad=0.3, loc='right')
  167      # plt.title('Average across trials')
  168      simpleaxis(plt.gca())
  169      plt.tight_layout(rect=(0, 0, 1, 0.97), w_pad=2)
  170      # plt.suptitle(f'{expid.split("_")[-2]}-{expid.split("_")[-1]}',
  171      #              fontsize=fontsize)
  172      plt.suptitle(f'{expid}',
  173                   fontsize=fontsize)
  174      if savedir is not None:
  175          plt.savefig(os.path.join(savedir, f'{expid}_eyeclosure.png'),
  176                      dpi=300)
  177      plt.show()

analysis.py (Zenodo 14876694, md5 checked), make_puffseries_behavior_plot, the lines that cut, baseline and normalise:
   38  def make_puffseries_behavior_plot(save_path, all_data, do_save=True, save_dpi=300, display_dpi=100, first_only=False, plot_late=True):
   55              trial_events = b['events']['on'] + 5*b['fps']
   68                              pertrial.append(b['eyeclosure'][p-int(1*b['fps']): p+int(2*b['fps'])+extra])
   72                          pertrial.append(b['eyeclosure'][p-int(1*b['fps']): p+int(2*b['fps'])+extra])
   73                  peranimal.append(np.concatenate(pertrial) - np.mean(pertrial[0][:100]))
   78              preket_behavior.append(np.nanmean(peranimal[:20,:],0).flatten()/np.nanmean(peranimal[:20,:],0).flatten().max())
   80              postket_behavior.append(np.nanmean(peranimal[22:32,:],0).flatten()/np.nanmean(peranimal[22:32,:],0).flatten().max())
  103      ste = np.std(np.stack(preket_behavior),0)/np.sqrt(np.stack(preket_behavior).shape[0])
  108      ste = np.std(np.stack(postket_behavior),0)/np.sqrt(np.stack(postket_behavior).shape[0])
  114          ste = np.std(np.stack(lateket_behavior),0)/np.sqrt(np.stack(lateket_behavior).shape[0])
fig1_npx_behavior.py main + plot_npx_eyeclosure: divides by a maximum: False
analysis.py make_puffseries_behavior_plot: divides by a maximum: True

def ci95(x, ddof=1):
    return 1.96 * np.std(x, ddof=ddof) / np.sqrt(len(x))


def fig4b_table(max_norm, ket_kind, agg):
    # Per session: first-puff late (stored aff_w) pre and ketamine, and accumulation (8th - 1st).
    rows = []
    for b in NPX.values():
        pre = b["trial_range_preket"]
        ket = b["trial_range_ket"] if ket_kind == "stored" else slice(22, 32)
        aw = b["aff_w"]
        t, mp, mk, per = me.npx_series(b, first_only=True, max_norm=max_norm, pre=pre, ket=ket, agg=agg)
        w = (t >= aw[0]) & (t < aw[1])
        t2, sp, sk, _ = me.npx_series(b, first_only=False, max_norm=max_norm, pre=pre, ket=ket, agg=agg)
        w1 = (t2 >= aw[0]) & (t2 < aw[1]); w8 = (t2 >= 21 + aw[0]) & (t2 < 21 + aw[1])
        rows.append(dict(session=b["subject_id"], late_pre=np.nanmean(mp[w]), late_ket=np.nanmean(mk[w]),
                         acc_pre=np.nanmean(sp[w8]) - np.nanmean(sp[w1]), acc_ket=np.nanmean(sk[w8]) - np.nanmean(sk[w1]),
                         nan_baseline_trials=int(np.isnan(per[:, :int(b["fps"][0])]).any(1).sum())))
    return pd.DataFrame(rows)


fig4b_res = {}
MAIN4B = ("not normalised (criterion)", "stored 23-30", "mean")
for mx, mlab in ((False, "not normalised (criterion)"), (True, "max-normalised (analysis.py)")):
    for kr, klab in (("stored", "stored 23-30"), ("code", "23-32 (analysis.py, Table S3)")):
        for agg in ("mean", "median"):
            df = fig4b_table(mx, kr, agg)
            n = len(df)
            lp, lk = df.late_pre.mean(), df.late_ket.mean()
            r = dict(late_pre=lp, late_ket=lk, acc_pre=df.acc_pre.mean(), acc_ket=df.acc_ket.mean(), n=n,
                     n_lower=int((df.late_ket < df.late_pre).sum()), n_acc_lower=int((df.acc_ket < df.acc_pre).sum()))
            for dd in (1, 0):
                r[f"ci_pre_{dd}"], r[f"ci_ket_{dd}"] = ci95(df.late_pre, dd), ci95(df.late_ket, dd)
                r[f"sep_{dd}"] = bool((lp - r[f"ci_pre_{dd}"]) > (lk + r[f"ci_ket_{dd}"]))
            dci = st.ttest_rel(df.late_pre, df.late_ket).confidence_interval()
            r["diff"], r["diff_ci"] = lp - lk, (dci.low, dci.high)
            fig4b_res[(mlab, klab, agg)] = r
            tag = "  <- main reading" if (mlab, klab, agg) == MAIN4B else ""
            print(f"{mlab}; ketamine trials {klab}; trials aggregated by {agg}; n = {n} sessions{tag}")
            print(f"   late, first puff: pre {lp:.3f} +/- {r['ci_pre_1']:.3f}, ketamine {lk:.3f} +/- {r['ci_ket_1']:.3f} "
                  f"(1.96 s.e.m.); intervals separate: {r['sep_1']}. With ddof 0 (analysis.py): +/- {r['ci_pre_0']:.3f}, "
                  f"+/- {r['ci_ket_0']:.3f}, separate: {r['sep_0']}")
            print(f"   pre minus ketamine {r['diff']:.3f}, paired 95% CI [{dci.low:.3f}, {dci.high:.3f}] (effect size; the paper "
                  f"reports no test); lower on ketamine in {r['n_lower']}/{n} sessions")
            print(f"   accumulation (late after 8th minus after 1st): pre {r['acc_pre']:.3f} +/- {ci95(df.acc_pre):.3f}, "
                  f"ketamine {r['acc_ket']:.3f} +/- {ci95(df.acc_ket):.3f}; lower on ketamine: {r['acc_ket'] < r['acc_pre']}, "
                  f"in {r['n_acc_lower']}/{n} sessions")
            if (mlab, klab, agg) == MAIN4B:
                print(df.round(3).to_string(index=False))
                print(f"   trials (of 40 per session) with a NaN frame in the baseline second, which np.mean in analysis.py "
                      f"would drop: {df.nan_baseline_trials.sum()} in total")

not normalised (criterion); ketamine trials stored 23-30; trials aggregated by mean; n = 10 sessions  <- main reading
   late, first puff: pre 0.174 +/- 0.061, ketamine 0.147 +/- 0.034 (1.96 s.e.m.); intervals separate: False. With ddof 0 (analysis.py): +/- 0.058, +/- 0.032, separate: False
   pre minus ketamine 0.027, paired 95% CI [-0.039, 0.094] (effect size; the paper reports no test); lower on ketamine in 6/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.058 +/- 0.028, ketamine 0.046 +/- 0.031; lower on ketamine: True, in 5/10 sessions
    session  late_pre  late_ket  acc_pre  acc_ket  nan_baseline_trials
 m11208-959     0.117     0.169    0.076    0.010                    1
 m11208-962     0.074     0.170    0.044    0.004                    0
 m11210-960     0.142     0.053    0.065    0.017                    0
m11219-1108     0.206     0.155    0.080    0.081                    0
m11221-1119     0.119     0.189    0.031    0.041                    0
m11223-1125     0.146     0.093    0.029    0.027                    0
m11225-1131     0.127     0.116    0.001    0.006                    0
m11227-1132     0.204     0.096    0.001    0.013                    1
   m304-876     0.430     0.224    0.121    0.131                    0
    m94-549     0.175     0.203    0.132    0.132                    0
   trials (of 40 per session) with a NaN frame in the baseline second, which np.mean in analysis.py would drop: 2 in total
not normalised (criterion); ketamine trials stored 23-30; trials aggregated by median; n = 10 sessions
   late, first puff: pre 0.169 +/- 0.063, ketamine 0.144 +/- 0.034 (1.96 s.e.m.); intervals separate: False. With ddof 0 (analysis.py): +/- 0.060, +/- 0.032, separate: False
   pre minus ketamine 0.026, paired 95% CI [-0.045, 0.096] (effect size; the paper reports no test); lower on ketamine in 6/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.050 +/- 0.025, ketamine 0.043 +/- 0.025; lower on ketamine: True, in 5/10 sessions
not normalised (criterion); ketamine trials 23-32 (analysis.py, Table S3); trials aggregated by mean; n = 10 sessions
   late, first puff: pre 0.174 +/- 0.061, ketamine 0.148 +/- 0.037 (1.96 s.e.m.); intervals separate: False. With ddof 0 (analysis.py): +/- 0.058, +/- 0.035, separate: False
   pre minus ketamine 0.026, paired 95% CI [-0.042, 0.094] (effect size; the paper reports no test); lower on ketamine in 6/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.058 +/- 0.028, ketamine 0.040 +/- 0.029; lower on ketamine: True, in 8/10 sessions

not normalised (criterion); ketamine trials 23-32 (analysis.py, Table S3); trials aggregated by median; n = 10 sessions
   late, first puff: pre 0.169 +/- 0.063, ketamine 0.143 +/- 0.037 (1.96 s.e.m.); intervals separate: False. With ddof 0 (analysis.py): +/- 0.060, +/- 0.035, separate: False
   pre minus ketamine 0.026, paired 95% CI [-0.048, 0.101] (effect size; the paper reports no test); lower on ketamine in 5/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.050 +/- 0.025, ketamine 0.037 +/- 0.028; lower on ketamine: True, in 6/10 sessions
max-normalised (analysis.py); ketamine trials stored 23-30; trials aggregated by mean; n = 10 sessions
   late, first puff: pre 0.630 +/- 0.080, ketamine 0.318 +/- 0.061 (1.96 s.e.m.); intervals separate: True. With ddof 0 (analysis.py): +/- 0.076, +/- 0.058, separate: True
   pre minus ketamine 0.312, paired 95% CI [0.228, 0.396] (effect size; the paper reports no test); lower on ketamine in 10/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.200 +/- 0.083, ketamine 0.087 +/- 0.045; lower on ketamine: True, in 8/10 sessions
max-normalised (analysis.py); ketamine trials stored 23-30; trials aggregated by median; n = 10 sessions
   late, first puff: pre 0.656 +/- 0.080, ketamine 0.312 +/- 0.064 (1.96 s.e.m.); intervals separate: True. With ddof 0 (analysis.py): +/- 0.076, +/- 0.061, separate: True
   pre minus ketamine 0.343, paired 95% CI [0.224, 0.463] (effect size; the paper reports no test); lower on ketamine in 10/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.169 +/- 0.055, ketamine 0.086 +/- 0.040; lower on ketamine: True, in 8/10 sessions
max-normalised (analysis.py); ketamine trials 23-32 (analysis.py, Table S3); trials aggregated by mean; n = 10 sessions
   late, first puff: pre 0.630 +/- 0.080, ketamine 0.322 +/- 0.054 (1.96 s.e.m.); intervals separate: True. With ddof 0 (analysis.py): +/- 0.076, +/- 0.051, separate: True
   pre minus ketamine 0.308, paired 95% CI [0.229, 0.387] (effect size; the paper reports no test); lower on ketamine in 10/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.200 +/- 0.083, ketamine 0.074 +/- 0.044; lower on ketamine: True, in 9/10 sessions

max-normalised (analysis.py); ketamine trials 23-32 (analysis.py, Table S3); trials aggregated by median; n = 10 sessions
   late, first puff: pre 0.656 +/- 0.080, ketamine 0.314 +/- 0.057 (1.96 s.e.m.); intervals separate: True. With ddof 0 (analysis.py): +/- 0.076, +/- 0.054, separate: True
   pre minus ketamine 0.342, paired 95% CI [0.231, 0.453] (effect size; the paper reports no test); lower on ketamine in 10/10 sessions
   accumulation (late after 8th minus after 1st): pre 0.169 +/- 0.055, ketamine 0.073 +/- 0.045; lower on ketamine: True, in 9/10 sessions

# Stored per-puff values: puff_aff (late, aff_w) and puff_ref (early, ref_w) on the first puff of each
# series, averaged over the stored trial ranges per session; across sessions mean and 1.96 s.e.m.
rows = []
for b in NPX.values():
    pa, pr = b["puff_aff"].reshape(-1, 8), b["puff_ref"].reshape(-1, 8)
    pre, ket = b["trial_range_preket"], b["trial_range_ket"]
    bp, tb = b["eyeclosure_by_puff"], b["t_eyeclosure"]
    raw_late = np.nanmean(bp[:, (tb >= 0.3 - 1e-9) & (tb < 0.8 - 1e-9)], 1)
    raw_early = np.nanmean(bp[:, (tb >= -1e-9) & (tb < 0.3 - 1e-9)], 1)
    rows.append(dict(session=b["subject_id"], aff_pre=np.nanmean(pa[pre, 0]), aff_ket=np.nanmean(pa[ket, 0]),
                     ref_pre=np.nanmean(pr[pre, 0]), ref_ket=np.nanmean(pr[ket, 0]),
                     ratio_pre=np.nanmean(pa[pre, 0] / pr[pre, 0]), ratio_ket=np.nanmean(pa[ket, 0] / pr[ket, 0]),
                     base_pre=np.nanmean(bp.reshape(-1, 8, bp.shape[1])[pre, 0][:, (tb >= -1) & (tb < 0)]),
                     base_ket=np.nanmean(bp.reshape(-1, 8, bp.shape[1])[ket, 0][:, (tb >= -1) & (tb < 0)]),
                     d_aff=np.nanmax(np.abs(raw_late - b["puff_aff"])), d_ref=np.nanmax(np.abs(raw_early - b["puff_ref"]))))
sdf = pd.DataFrame(rows)
print(sdf.drop(columns=["d_aff", "d_ref"]).round(3).to_string(index=False))
print(f"stored values against the raw (not baseline-subtracted) mean of eyeclosure_by_puff: puff_ref over 0-0.3 s, max "
      f"difference {sdf.d_ref.max():.4f}; puff_aff over 0.3-0.8 s, max difference {sdf.d_aff.max():.4f}")
for q, lab in (("aff", "puff_aff"), ("ref", "puff_ref"), ("ratio", "puff_aff/puff_ref"), ("base", "raw 1 s before the first puff")):
    a, k = sdf[f"{q}_pre"], sdf[f"{q}_ket"]
    print(f"{lab:30s} pre {a.mean():.3f} +/- {ci95(a):.3f}, ketamine {k.mean():.3f} +/- {ci95(k):.3f} (1.96 s.e.m.); "
          f"intervals separate: {(a.mean() - ci95(a)) > (k.mean() + ci95(k))}; lower on ketamine in {(k < a).sum()}/{len(a)}")

    session  aff_pre  aff_ket  ref_pre  ref_ket  ratio_pre  ratio_ket  base_pre  base_ket
 m11208-959    0.616    0.303    0.684    0.528      0.914      0.575     0.499     0.129
 m11208-962    0.830    0.339    0.851    0.634      0.973      0.519     0.756     0.163
 m11210-960    0.733    0.336    0.851    0.438      0.862      0.769     0.588     0.282
m11219-1108    0.704    0.341    0.726    0.546      0.966      0.631     0.497     0.186
m11221-1119    0.799    0.414    0.735    0.472      1.099      0.881     0.679     0.223
m11223-1125    0.787    0.283    0.850    0.439      0.925      0.654     0.639     0.186
m11225-1131    0.762    0.429    0.825    0.714      0.924      0.602     0.634     0.307
m11227-1132    0.638    0.283    0.661    0.508      0.965      0.571     0.436     0.182
   m304-876    0.754    0.440    0.764    0.723      0.985      0.609     0.324     0.210
    m94-549    0.749    0.413    0.866    0.782      0.865      0.532     0.581     0.222
stored values against the raw (not baseline-subtracted) mean of eyeclosure_by_puff: puff_ref over 0-0.3 s, max difference 0.0000; puff_aff over 0.3-0.8 s, max difference 0.0164
puff_aff                       pre 0.737 +/- 0.042, ketamine 0.358 +/- 0.038 (1.96 s.e.m.); intervals separate: True; lower on ketamine in 10/10
puff_ref                       pre 0.781 +/- 0.047, ketamine 0.578 +/- 0.078 (1.96 s.e.m.); intervals separate: True; lower on ketamine in 10/10
puff_aff/puff_ref              pre 0.948 +/- 0.042, ketamine 0.634 +/- 0.069 (1.96 s.e.m.); intervals separate: True; lower on ketamine in 10/10
raw 1 s before the first puff  pre 0.563 +/- 0.079, ketamine 0.209 +/- 0.033 (1.96 s.e.m.); intervals separate: True; lower on ketamine in 10/10

fig, ax = plt.subplots(1, 2, figsize=(8, 2.6))
rng = []
for a, mx in zip(ax, (False, True)):
    P_, K_ = [], []
    for b in NPX.values():
        t, mp, mk, _ = me.npx_series(b, first_only=True, max_norm=mx, pre=b["trial_range_preket"], ket=b["trial_range_ket"])
        P_.append(mp); K_.append(mk)
    P_, K_ = np.array(P_), np.array(K_)
    for X, lab in ((P_, "preinfusion"), (K_, "ketamine")):
        m, h = np.nanmean(X, 0), 1.96 * np.nanstd(X, 0, ddof=1) / np.sqrt(len(X))
        a.plot(t, m, label=lab); a.fill_between(t, m - h, m + h, alpha=0.25, linewidth=0)
        rng.append(f"{'max-normalised' if mx else 'not normalised'} {lab}: {np.nanmin(m - h):.2f} to {np.nanmax(m + h):.2f}")
    a.axvspan(0.3, 0.8, color="0.9", zorder=0)
    a.set_title("max-normalised (analysis.py)" if mx else "not normalised (criterion)"); a.set_xlabel("time from first puff (s)")
ax[0].legend(frameon=False, fontsize=7); ax[0].set_ylabel("eye closure")
plt.tight_layout(); plt.show()
print("plotted ranges (mean +/- 1.96 s.e.m.), all inside the axes:", "; ".join(rng))

plotted ranges (mean +/- 1.96 s.e.m.), all inside the axes: not normalised preinfusion: -0.01 to 0.33; not normalised ketamine: -0.00 to 0.55; max-normalised preinfusion: -0.07 to 1.05; max-normalised ketamine: -0.01 to 1.00

print(f"runtime {(time.time() - T0) / 60:.1f} min")

runtime 2.5 min

Claim	Data	Status
mouse-session-labels-link	with labels	tested
fig1k-mouse-late-eye-closure-ketamine	with labels	tested
fig1l-mouse-dissociatives-reduce-late	with labels	tested
figs2bc-mouse-ketamine-early-present-late-abolished	with labels	tested
fig1m-mouse-late-closure-accumulates-and-persists	with labels	tested
figs2g-mouse-scaling-with-puff-number-not-pressure	with labels	tested
figs2h-mouse-dose-threshold	with labels	tested
figs5c-mouse-other-drugs-no-disconnection	with labels	tested
fig1p-mouse-anesthesia-abolishes-both	with labels	tested
fig1no-mouse-dissociatives-no-accumulation	with labels	tested
figs2l-mouse-repeatable-across-days	with labels	tested
figs2n-mouse-no-sex-difference	with labels	tested
figs2o-mouse-no-anticipation	with labels	tested
figs4-mouse-reward-seeking-reduced	dataset	blocked: lick data not in the dataset
figs6-mouse-vpm-inhibition-reduces-affective	with labels	tested
fig4b-mouse-recorded-affective-reduced	dataset	tested (10 of 13 sessions)

Claim	Outcome	Key numbers
mouse-session-labels-link	holds	all listed ids present, lengths as read, ketamine/PCP drops, K/X near zero, skip labels match
fig1k-mouse-late-eye-closure-ketamine	reproduced	saline 0.262, ketamine -0.068, t = 4.737, p = 0.009
fig1l-mouse-dissociatives-reduce-late	reproduced	t = 0.73 / 9.53 / 3.41, BH p = 0.504 / 0.002 / 0.041
figs2bc-mouse-ketamine-early-present-late-abolished	partially reproduced	tests as Table S3; saline late 0.203 against 0.34
fig1m-mouse-late-closure-accumulates-and-persists	reproduced (reading b)	20 trials: median 15.45 s (b), 13.05 s (a); the code's 40 trials and start frame give 15.16 s, the paper's value
figs2g-mouse-scaling-with-puff-number-not-pressure	reproduced	statistics as Table S3
figs2h-mouse-dose-threshold	reproduced	10 of 10 conditions as claimed
figs5c-mouse-other-drugs-no-disconnection	reproduced	ketamine 0.0053, PCP 0.0092, bup 0.93, dzp 0.30
fig1p-mouse-anesthesia-abolishes-both	reproduced	K/X early 0.025, late 0.024
fig1no-mouse-dissociatives-no-accumulation	reproduced (assumed, test reading of "not above")	decay rates as Table S3; late on 8th minus 1st puff +0.017 [-0.012, +0.046] (ketamine), +0.010 [-0.019, +0.039] (PCP), saline +0.350 [+0.167, +0.534]
figs2l-mouse-repeatable-across-days	reproduced	0 of 12 tests significant
figs2n-mouse-no-sex-difference	reproduced	0 of 4 tests significant
figs2o-mouse-no-anticipation	reproduced	no rise at the skipped puff, closure above baseline
figs4-mouse-reward-seeking-reduced	blocked (data missing)	no lick data in the dataset
figs6-mouse-vpm-inhibition-reduces-affective	reproduced	saline 0.602, muscimol 0.105
fig4b-mouse-recorded-affective-reduced	not reproduced (10 of 13 sessions)	criterion's method: 0.174 +/- 0.061 against 0.147 +/- 0.034, intervals overlap; separate only with analysis.py's max normalisation (0.630 against 0.318)

Mouse eyepuff behaviour: Kauvar, Richman, Liu et al. (2025), dandiset 001347¶

Data check¶

The code record and its experiment lists¶

The head-fixed mouse files¶

The Neuropixels-rig files¶

Link claim: mouse-session-labels-link¶

fig1k-mouse-late-eye-closure-ketamine¶

fig1l-mouse-dissociatives-reduce-late¶

figs2bc-mouse-ketamine-early-present-late-abolished¶

fig1m-mouse-late-closure-accumulates-and-persists¶

figs2g-mouse-scaling-with-puff-number-not-pressure¶

figs2h-mouse-dose-threshold¶

figs5c-mouse-other-drugs-no-disconnection¶

fig1p-mouse-anesthesia-abolishes-both¶

fig1no-mouse-dissociatives-no-accumulation¶

figs2l-mouse-repeatable-across-days¶

figs2n-mouse-no-sex-difference¶

figs2o-mouse-no-anticipation¶

figs4-mouse-reward-seeking-reduced (blocked)¶

figs6-mouse-vpm-inhibition-reduces-affective¶

fig4b-mouse-recorded-affective-reduced¶

Summary¶