import time
import warnings
import numpy as np
import pandas as pd
import scipy
import scipy.cluster
import scipy.spatial
import scipy.stats
import kauvar as K

warnings.filterwarnings("ignore", category=RuntimeWarning)
print(f"dandiset {K.DANDISET}, version {K.VERSION}")
urls = K.asset_urls()
paths = K.allen_tree()
res, meta = {}, {}
t0 = time.time()
for p, u in urls:
    s = K.load_session(p, u)
    meta[s.name] = dict(subject=s.subject, units=len(s.spikes), unsorted_units=s.n_unsorted_units, spikes=s.n_spikes_total,
                        n_spikes_column=int(s.n_spikes.sum()), snr_min=s.snr.min(), n_spikes_min=int(s.n_spikes.min()),
                        quality=sorted(set(s.quality)), duration_s=s.maxt, stim=s.stim.stimulus.value_counts().to_dict(),
                        stim_dur=s.stim.groupby("stimulus").apply(lambda g: (g.stop_time - g.start_time).agg(["min", "max"]).round(4).tolist()).to_dict(),
                        puff_isi=np.round(np.percentile(np.diff(s.onsets("puff").reshape(-1, 8), axis=1), [0, 100]), 3).tolist(),
                        series_gap=np.round(np.percentile(np.diff(s.onsets("puff").reshape(-1, 8)[:, 0]), [0, 100]), 1).tolist(),
                        **s.meta)
    res[s.name] = K.process_session(s, paths)
    del s
names = list(res)
print(f"{len(names)} sessions loaded and reduced in {time.time() - t0:.0f} s")

dandiset 001326, version 0.250528.1957

13 sessions loaded and reduced in 200 s

m = pd.DataFrame(meta).T
print("Units table: every unit is labelled", sorted({q for v in m.quality for q in v}), "; smallest SNR",
      round(m.snr_min.min(), 2), "(paper: SNR above 2); smallest spike count", m.n_spikes_min.min(), "(paper: more than 500)")
print("n_spikes column equals the stored spike times in every file:", bool((m.spikes == m.n_spikes_column).all()))
print("units whose spike times are not in ascending order:", int(m.unsorted_units.sum()), "of", int(m.units.sum()))
print(m[["subject", "units", "spikes", "duration_s", "devices", "age", "sex"]].to_string())
print("\nStimulus presentations per session (count) and duration range (s):")
print(pd.DataFrame({k: v for k, v in m.stim.items()}).T.fillna(0).astype(int).to_string())
print(pd.Series(m.stim_dur).to_string())
print("\nPuff onset intervals within a series (s, min and max over sessions):", [min(v[0] for v in m.puff_isi), max(v[1] for v in m.puff_isi)])
print("Intervals between series starts (s, min and max):", [min(v[0] for v in m.series_gap), max(v[1] for v in m.series_gap)])
rows = []
for k in names:
    o = res[k]
    rows.append({"session": k, "subject": o["subject"], "puff series": o["n_series"], "series of 8": set(o["series_sizes"]) == {8},
                 "code ketamine onset (s)": round(o["ket_start_code"], 1), "infusion in file (s)": round(o["ket_infusion_nwb"], 1),
                 "tone series": len(o.get("tone_series_sizes", [])), "regions (finest, no NA)": len(set(o["regions_all"]) - {"NA"}),
                 "units": o["n_units"], "puff-modulated": o["n_mod"],
                 "clock puff true": round(o["clock_puff"]["true"], 3), "clock puff shifted": round(o["clock_puff"]["shifted 1.5 s"], 3),
                 "clock tone true": round(o["clock_tone"]["true"], 3) if "clock_tone" in o else np.nan,
                 "clock tone shifted": round(o["clock_tone"]["shifted 1.5 s"], 3) if "clock_tone" in o else np.nan})
sess = pd.DataFrame(rows).set_index("session")
print(sess.to_string())
print(f"\nsessions {len(names)}, mice {sess.subject.nunique()}, units {sess.units.sum()}, puff-modulated {sess['puff-modulated'].sum()}")
print("every session has at least 40 puff series:", bool((sess["puff series"] >= 40).all()),
      "; the infusion in the file falls after series 20 in every session:", bool((sess["infusion in file (s)"] > sess["code ketamine onset (s)"] - 30).all()))
print("rates: z-scored per unit; z values of the stored windows range",
      round(min(res[k]["pre_avg"].min() for k in names), 2), "to", round(max(res[k]["pre_avg"].max() for k in names), 2), "(Pre averages)")

Units table: every unit is labelled [np.str_('good')] ; smallest SNR 2.0 (paper: SNR above 2); smallest spike count 517 (paper: more than 500)
n_spikes column equals the stored spike times in every file: True
units whose spike times are not in ascending order: 9554 of 9554
                   subject units    spikes   duration_s                                                                               devices       age sex
sub-0_ses-20230901       0   942  30275497  7860.708509  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-1_ses-20230218       1   583  17014893  5881.223584  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-1_ses-20230226       1  1020  45462064  8090.584986  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-2_ses-20230902       2   708  16019860  7920.803845  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-3_ses-20230218       3   431  10144317  5881.294849  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-3_ses-20230225       3   917  35381600  7947.431481  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-4_ses-20230907       4   885  29845185  8103.080661  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-5_ses-20230706       5   154   2753026  8101.440109                                            [Neuropixels_Probe_0, Neuropixels_Probe_3]  P8W/P12W   F
sub-5_ses-20230822       5   664  16205211   7921.21314  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-6_ses-20230823       6   591  17459504  7920.848412  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-7_ses-20220831       7   928  37991946  7537.026668  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-8_ses-20220228       8   851  42459218  9174.911337  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F
sub-9_ses-20230908       9   880  23843952   8402.29841  [Neuropixels_Probe_0, Neuropixels_Probe_1, Neuropixels_Probe_2, Neuropixels_Probe_3]  P8W/P12W   F

Stimulus presentations per session (count) and duration range (s):
                    puff  light  tone
sub-0_ses-20230901   480    240   232
sub-1_ses-20230218   400      0     0
sub-1_ses-20230226   488    248   240
sub-2_ses-20230902   480    240   232
sub-3_ses-20230218   400      0     0
sub-3_ses-20230225   496    248   248
sub-4_ses-20230907   480    240   240
sub-5_ses-20230706   480    240   240
sub-5_ses-20230822   480    240   232
sub-6_ses-20230823   480    240   232
sub-7_ses-20220831   416      0     0
sub-8_ses-20220228   320      0     0
sub-9_ses-20230908   480    240   240
sub-0_ses-20230901    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-1_ses-20230218                           {'puff': ['0.25', '0.25']}
sub-1_ses-20230226    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-2_ses-20230902    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-3_ses-20230218                           {'puff': ['0.25', '0.25']}
sub-3_ses-20230225    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-4_ses-20230907    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-5_ses-20230706    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-5_ses-20230822    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-6_ses-20230823    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...
sub-7_ses-20220831                           {'puff': ['0.25', '0.25']}
sub-8_ses-20220228                           {'puff': ['0.25', '0.25']}
sub-9_ses-20230908    {'light': ['0.25', '0.25'], 'puff': ['0.25', '...

Puff onset intervals within a series (s, min and max over sessions): [2.63, 3.38]
Intervals between series starts (s, min and max): [106.4, 475.7]
                   subject  puff series  series of 8  code ketamine onset (s)  infusion in file (s)  tone series  regions (finest, no NA)  units  puff-modulated  clock puff true  clock puff shifted  clock tone true  clock tone shifted
session                                                                                                                                                                                                                                   
sub-0_ses-20230901       0           60         True                   2747.6                2841.6           29                       31    942             767            0.171              -0.015            0.128              -0.015
sub-1_ses-20230218       1           50         True                   2401.4                2500.2            0                       32    583             408            0.141              -0.003              NaN                 NaN
sub-1_ses-20230226       1           61         True                   2645.0                3023.1           30                       33   1020             641            0.092              -0.012            0.098              -0.002
sub-2_ses-20230902       2           60         True                   2748.4                2846.4           29                       33    708             461            0.063               0.003            0.050              -0.019
sub-3_ses-20230218       3           50         True                   2403.4                2528.2            0                       24    431             316            0.075              -0.027              NaN                 NaN
sub-3_ses-20230225       3           62         True                   2636.8                2764.4           31                       33    917             589            0.057              -0.018            0.043              -0.005
sub-4_ses-20230907       4           60         True                   2751.0                2844.9           30                       33    885             661            0.120              -0.016            0.057               0.004
sub-5_ses-20230706       5           60         True                   2747.3                2861.7           30                       11    154              71            0.033              -0.013            0.015              -0.003
sub-5_ses-20230822       5           60         True                   2751.7                2886.6           29                       20    664             405            0.068              -0.016            0.041              -0.020
sub-6_ses-20230823       6           60         True                   2750.7                2860.6           29                       22    591             429            0.129               0.002            0.156              -0.004
sub-7_ses-20220831       7           52         True                   2283.8                2703.4            0                       29    928             645            0.089              -0.002              NaN                 NaN
sub-8_ses-20220228       8           40         True                   2285.6                2445.9            0                       37    851             675            0.194              -0.025              NaN                 NaN
sub-9_ses-20230908       9           60         True                   2747.5                2989.5           30                       35    880             630            0.104              -0.010            0.051              -0.011

sessions 13, mice 10, units 9554, puff-modulated 6698
every session has at least 40 puff series: True ; the infusion in the file falls after series 20 in every session: True
rates: z-scored per unit; z values of the stored windows range -1.83 to 18.91 (Pre averages)

# Added in revision. Epochs against the file's infusion time: the first puff of series 20 to 24 relative to
# the "Ketamine infusion start" of each file (negative: before it). Ketamine is series 23 to 32 (code), the
# Methods-literal reading 22 to 32.
ep = pd.DataFrame({k: {f"series {s + 1}": res[k]["series_rel_infusion"][s, 0] for s in range(19, 24)} for k in names}).T
ep["code onset - infusion"] = [res[k]["ket_start_code"] - res[k]["ket_infusion_nwb"] for k in names]
print(ep.round(1).to_string())
print("sessions where series 23 starts before the file's infusion time:", list(ep.index[ep["series 23"] < 0]))
print("sessions where series 22 starts before the file's infusion time:", list(ep.index[ep["series 22"] < 0]))
# tone and light onsets inside the windows of the analyses: the 30 s before puff series (fig6h, fig6j, Pre 1-20
# and Ketamine 23-32) and the 50 s after each puff series (fig5d, fig5g, figs33)
rows = []
for k in names:
    for kind in ("tone", "light"):
        w = res[k]["stim_in_windows"][kind]
        for wn, v in w.items():
            used = np.r_[0:20, 22:32]
            rows.append({"session": k, "stimulus": kind, "window": wn, "series with an onset inside (of 30 used)": int((v[used] > 0).sum()),
                         "onsets inside": int(v[used].sum())})
sw = pd.DataFrame(rows)
print("\ntone and light onsets inside the windows of Pre series 1-20 and Ketamine series 23-32:")
print(sw.pivot_table(index="session", columns=["stimulus", "window"], values="series with an onset inside (of 30 used)").to_string())
print("total onsets inside, by stimulus and window:", sw.groupby(["stimulus", "window"])["onsets inside"].sum().to_dict())
# added in revision 2: series with a tone or a light onset in the 30 s before them, per session
either = {k: int(((res[k]["stim_in_windows"]["tone"]["30 s before series"][used] > 0)
                  | (res[k]["stim_in_windows"]["light"]["30 s before series"][used] > 0)).sum()) for k in names}
print("series (of 30 used) with a tone or light onset in the 30 s before them:", either)
pattern = {kind: sorted({tuple(int(x) + 1 for x in used if res[k]["stim_in_windows"][kind]["30 s before series"][x] > 0) for k in names} - {()})
           for kind in ("tone", "light")}
print("series numbers with an onset in the 30 s before them, in the sessions that have any:", pattern)

                    series 20  series 21  series 22  series 23  series 24  code onset - infusion
sub-0_ses-20230901     -145.0        5.0      134.5      265.6      390.3                  -94.0
sub-1_ses-20230218     -149.8        5.0      119.4      235.5      345.1                  -98.8
sub-1_ses-20230226     -429.1        5.0      129.5      260.6      380.2                 -378.1
sub-2_ses-20230902     -149.0        5.0      134.5      265.6      390.3                  -98.0
sub-3_ses-20230218     -175.8        5.0      119.4      235.5      345.1                 -124.8
sub-3_ses-20230225     -178.6        5.0      129.5      260.6      380.2                 -127.6
sub-4_ses-20230907     -144.9        5.0      134.5      265.7      390.3                  -93.9
sub-5_ses-20230706     -165.4        5.0      134.5      265.6      390.3                 -114.4
sub-5_ses-20230822     -185.9        5.0      134.5      265.6      390.3                 -134.9
sub-6_ses-20230823     -160.9        5.0      134.5      265.6      390.2                 -109.9
sub-7_ses-20220831     -470.7        5.0      119.5      235.7      345.3                 -419.6
sub-8_ses-20220228     -211.4        5.0      119.5      235.7      345.3                 -160.3
sub-9_ses-20230908     -293.0        5.0      134.5      265.6      390.3                 -242.0
sessions where series 23 starts before the file's infusion time: []
sessions where series 22 starts before the file's infusion time: []

tone and light onsets inside the windows of Pre series 1-20 and Ketamine series 23-32:
stimulus                        light                                 tone                  
window             30 s before series 50 s after series 30 s before series 50 s after series
session                                                                                     
sub-0_ses-20230901               10.0               0.0                4.0               0.0
sub-1_ses-20230218                0.0               0.0                0.0               0.0
sub-1_ses-20230226               10.0              15.0                4.0               0.0
sub-2_ses-20230902               10.0               0.0                4.0               0.0
sub-3_ses-20230218                0.0               0.0                0.0               0.0
sub-3_ses-20230225               10.0              15.0                4.0               0.0
sub-4_ses-20230907               10.0               0.0                4.0               0.0
sub-5_ses-20230706               10.0               0.0                4.0               0.0
sub-5_ses-20230822               10.0               0.0                4.0               0.0
sub-6_ses-20230823               10.0               0.0                4.0               0.0
sub-7_ses-20220831                0.0               0.0                0.0               0.0
sub-8_ses-20220228                0.0               0.0                0.0               0.0
sub-9_ses-20230908               10.0               0.0                4.0               0.0
total onsets inside, by stimulus and window: {('light', '30 s before series'): 108, ('light', '50 s after series'): 30, ('tone', '30 s before series'): 45, ('tone', '50 s after series'): 0}
series (of 30 used) with a tone or light onset in the 30 s before them: {'sub-0_ses-20230901': 14, 'sub-1_ses-20230218': 0, 'sub-1_ses-20230226': 14, 'sub-2_ses-20230902': 14, 'sub-3_ses-20230218': 0, 'sub-3_ses-20230225': 14, 'sub-4_ses-20230907': 14, 'sub-5_ses-20230706': 14, 'sub-5_ses-20230822': 14, 'sub-6_ses-20230823': 14, 'sub-7_ses-20220831': 0, 'sub-8_ses-20220228': 0, 'sub-9_ses-20230908': 14}
series numbers with an onset in the 30 s before them, in the sessions that have any: {'tone': [(5, 7, 9, 29)], 'light': [(4, 6, 8, 14, 16, 20, 24, 28, 30, 32)]}

n_sessions, n_mice = len(names), sess.subject.nunique()
n_mod = int(sess["puff-modulated"].sum())
print(f"sessions: {n_sessions} (expected 13); mice: {n_mice} (expected 10)")
print(f"units after the quality criteria (all units in the files): {sess.units.sum()} (range per session {sess.units.min()}-{sess.units.max()})")
print(f"puff-modulated units: {n_mod} (expected 6969, tolerance 697); with the spike-count reading of the 250 minimum: "
      f"{sum(res[k]['sel_count_reading'] for k in names)}")
print("per test (units passing, summed over sessions):", pd.DataFrame({k: res[k]["n_test"] for k in names}).sum(1).to_dict())
reg = sess["regions (finest, no NA)"]
print(f"regions per session at the finest Allen level: median {reg.median():.0f}, range {reg.min()}-{reg.max()}; "
      f"sessions with at least 24: {(reg >= 24).sum()} of {len(reg)}")

sessions: 13 (expected 13); mice: 10 (expected 10)
units after the quality criteria (all units in the files): 9554 (range per session 154-1020)
puff-modulated units: 6698 (expected 6969, tolerance 697); with the spike-count reading of the 250 minimum: 6491
per test (units passing, summed over sessions): {'pre vs 0-200 ms': 4931, 'pre vs 300-800 ms': 4533, 'pre vs 1-2 s': 1257, 'baseline first vs last puff': 310, 'response first vs last puff': 1345}
regions per session at the finest Allen level: median 32, range 11-37; sessions with at least 24: 10 of 13

/tmp/ipykernel_2286/626211601.py:7: Pandas4Warning: Starting with pandas version 4.0 all arguments of sum will be keyword-only.
  print("per test (units passing, summed over sessions):", pd.DataFrame({k: res[k]["n_test"] for k in names}).sum(1).to_dict())

import scanpy as sc

pre_avg = np.vstack([res[k]["pre_avg"] for k in names])
ket_avg = np.vstack([res[k]["ket_avg"] for k in names])
mod_regions = np.hstack([res[k]["regions"] for k in names])
mod_session = np.hstack([[k] * res[k]["n_mod"] for k in names])
X = np.nan_to_num(pre_avg - pre_avg[:, :100].mean(1, keepdims=True)).astype(np.float32)
# Added in revision 2. The float32 PCA rounds differently with the number of BLAS threads, and that is enough to
# change the neighbor graph and the partition (on the same data, 1, 2 and 4 threads gave three different
# partitions). BLAS is limited to one thread here, so the partition does not depend on the machine; this is a
# standard choice for a detail of computation, not a method choice. The PCA in float64 is shown in the stability
# section as a rounding-level perturbation.
from threadpoolctl import threadpool_limits
with threadpool_limits(1):
    np.random.seed(42)
    ann = sc.AnnData(X)
    sc.tl.pca(ann, svd_solver="arpack")
    sc.pp.neighbors(ann, n_neighbors=30, n_pcs=20)
    sc.tl.umap(ann)
    sc.tl.leiden(ann, flavor="leidenalg")
labels = ann.obs["leiden"].astype(int).to_numpy()
nclus = labels.max() + 1
sizes = np.bincount(labels)
expected_sizes = [1070, 832, 803, 663, 639, 597, 547, 508, 473, 344, 317, 176]
print(f"puff-modulated units clustered: {len(labels)}; clusters: {nclus} (expected 12, tolerance 1)")
srt = np.sort(sizes)[::-1]
rows = []
for i, e in enumerate(expected_sizes):
    got = srt[i] if i < len(srt) else 0
    rows.append({"rank": i, "size": got, "expected": e, "tolerance": round(0.1 * e), "within": abs(got - e) <= round(0.1 * e)})
print(pd.DataFrame(rows).to_string(index=False))

def cluster_means(avg, labels):
    """analysis.py get_cell_clusters: per cell z-scored over -1..2 s, averaged, minus the mean of -1..0 s."""
    out = []
    for c in range(labels.max() + 1):
        zz = np.nanmean(scipy.stats.zscore(avg[labels == c], axis=1), 0)
        out.append(zz - zz[:100].mean())
    return np.array(out)

def peak_decay(trace, frac=0.25):
    """analysis.py find_decay: |trace| normalized by its largest |value| after onset; the peak is its argmax and
    the decay the time from onset until it falls below frac after the peak; a peak whose z-score (over the
    trace) is below 1 in absolute value gets decay 2 s, as the code does. Returns peak value (signed),
    peak time, decay."""
    zt = scipy.stats.zscore(trace)
    a = np.abs(trace / np.nanmax(np.abs(trace[100:])))
    pk = int(np.argmax(a[100:]))
    below = np.where(a[100:][pk:] < frac)[0]
    decay = (pk + below[0]) * K.BIN if len(below) else np.nan
    if abs(zt[100:][pk]) < 1:
        decay = 2.0
    return trace[100 + pk], pk * K.BIN, decay

cm = cluster_means(X, labels)
# hierarchical reordering as get_cell_clusters (moving average of 20 bins, correlation, cosine, average linkage)
cp = np.array([K.moving_avg(c, 20) for c in cm])
link = scipy.cluster.hierarchy.linkage(scipy.spatial.distance.pdist(np.nan_to_num(np.corrcoef(cp)), metric="cosine"), method="average")
reordering = scipy.cluster.hierarchy.leaves_list(link)
disp = {c: i for i, c in enumerate(reordering)}       # figure number of each size-ranked cluster
shape = []
for c in range(nclus):
    pk, tpk, dec = peak_decay(cm[c])
    late = np.abs(cm[c][200:300]).mean()
    shape.append({"cluster (size rank)": c, "figure number": disp[c], "cells": sizes[c], "peak": round(pk, 2),
                  "peak time s": tpk, "decay to 25% s": dec, "mean |1-2 s| / |peak|": round(late / abs(pk), 2), "min |1-2 s| / |peak|": round(np.abs(cm[c][200:300]).min() / abs(pk), 2),
                  "rapid and brief": (dec <= 0.25) if not np.isnan(dec) else False,
                  "persistent": late > 0.25 * abs(pk), "sign": "up" if pk > 0 else "down"})
shape = pd.DataFrame(shape)
print(shape.to_string(index=False))
rb = shape["rapid and brief"].sum()
pers = shape[shape.persistent]
print(f"rapid and brief clusters: {rb}; persistent clusters: {len(pers)}, of which peak within 250 ms: "
      f"{(pers['peak time s'] < 0.25).sum()}, later: {(pers['peak time s'] >= 0.25).sum()}")
for numbering in ("cluster (size rank)", "figure number"):
    sub = shape.set_index(numbering).loc[[0, 1, 2]]
    print(f"clusters 0, 1, 2 by {numbering}: persistent {list(sub.persistent)}, peak times {list(sub['peak time s'])}")
print("cells per cluster by session:")
print(pd.crosstab(mod_session, labels).to_string())

/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

/tmp/ipykernel_2286/535742623.py:20: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg")

puff-modulated units clustered: 6698; clusters: 11 (expected 12, tolerance 1)
 rank  size  expected  tolerance  within
    0   963      1070        107    True
    1   772       832         83    True
    2   664       803         80   False
    3   657       663         66    True
    4   651       639         64    True
    5   604       597         60    True
    6   590       547         55    True
    7   579       508         51   False
    8   554       473         47   False
    9   491       344         34   False
   10   173       317         32   False
   11     0       176         18   False
 cluster (size rank)  figure number  cells  peak  peak time s  decay to 25% s  mean |1-2 s| / |peak|  min |1-2 s| / |peak|  rapid and brief  persistent sign
                   0              0    963 -1.06         0.11            1.49                   0.26                  0.07            False        True down
                   1              3    772  1.40         0.52            1.35                   0.22                  0.08            False       False   up
                   2              5    664  3.65         0.09            0.49                   0.03                  0.01            False       False   up
                   3              6    657  2.10         0.19            0.58                   0.09                  0.04            False       False   up
                   4              1    651 -1.90         0.35            1.16                   0.14                  0.00            False       False down
                   5              2    604  2.26         0.50            1.24                   0.18                  0.03            False       False   up
                   6              9    590  3.28         0.07            0.18                   0.02                  0.00             True       False   up
                   7              8    579  3.99         0.10            0.32                   0.01                  0.00            False       False   up
                   8              7    554  2.96         0.29            0.55                   0.02                  0.00            False       False   up
                   9             10    491  6.29         0.04            0.10                   0.01                  0.00             True       False   up
                  10              4    173  0.49         0.77            0.86                   0.42                  0.12            False        True   up
rapid and brief clusters: 2; persistent clusters: 2, of which peak within 250 ms: 1, later: 1
clusters 0, 1, 2 by cluster (size rank): persistent [True, False, False], peak times [0.11, 0.52, 0.09]
clusters 0, 1, 2 by figure number: persistent [True, False, False], peak times [0.11, 0.35000000000000003, 0.5]
cells per cluster by session:
col_0                0   1    2    3   4   5   6   7   8   9   10
row_0                                                            
sub-0_ses-20230901   78  62  116   48  91  53  67  93  78  65  16
sub-1_ses-20230218   50  30   42   43  37  21  56  54  26  42   7
sub-1_ses-20230226  106  74   28  102  33  15  67  76  56  60  24
sub-2_ses-20230902  104  61   19   58  36  23  48  27  40  30  15
sub-3_ses-20230218   39  65   10   41  13  57  21  24  28   9   9
sub-3_ses-20230225  116  76   34   53  70  74  38  43  40  32  13
sub-4_ses-20230907   74  87   61   50  74  54  60  69  70  47  15
sub-5_ses-20230706   13  13    9    3   8   4   5   2   3   5   6
sub-5_ses-20230822   69  49   24   48  35  53  29  32  32  13  21
sub-6_ses-20230823   49  25   78   20  67  26  45  35  36  37  11
sub-7_ses-20220831  120  98   67   58  64  80  52  33  33  26  14
sub-8_ses-20220228   75  58  125   52  59  90  49  46  45  68   8
sub-9_ses-20230908   70  74   51   81  64  54  53  45  67  57  14

def cluster_trace(avg, mask):
    """Table S3, Fig. 4G: cells of the cluster pooled, baseline-subtracted (1 s before the puff), z-scored over -1..2 s."""
    m = np.nanmean(avg[mask] - avg[mask][:, :100].mean(1, keepdims=True), 0)
    return scipy.stats.zscore(m)

FRONTAL = ["ORB", "ACA", "PL", "ILA", "MOs", "AI"]
def is_frontal(r):
    return any(K.under(r, f, paths) for f in FRONTAL)

def overrep(labels, nclus, seed=2024):
    """analysis.py get_cluster_null_distribution: per region with more than 30 clustered cells, the fraction
    in each cluster against 20,000 uniformly random labelings (label 0 of the draws), two-sided; BH-FDR
    across regions. Returns the corrected p (regions x clusters), the fractions and the regions."""
    rng = np.random.default_rng(seed)
    valid = mod_regions != "NA"
    ureg, ucnt = np.unique(mod_regions[valid], return_counts=True)
    big = ureg[ucnt > 30]
    pvals, fracs = {}, {}
    for r in big:
        ncell = (mod_regions == r).sum()
        null = (rng.integers(nclus, size=(20000, ncell)) == 0).mean(1)
        obs = np.bincount(labels[mod_regions == r], minlength=nclus) / ncell
        pvals[r] = np.array([np.mean(np.abs(null - 1 / nclus) >= abs(o - 1 / nclus)) for o in obs])
        fracs[r] = obs
    P = pd.DataFrame(pvals).T
    F = pd.DataFrame(fracs).T
    Q = P.apply(lambda col: pd.Series(scipy.stats.false_discovery_control(np.maximum(col.values, 1 / 20000)), index=col.index))
    return Q, F, big

def figs31b_test(Q, F, big, nclus, dec, thr):
    """PAG and VPM in a faster-half cluster and a frontal/limbic region in a slower-half cluster, with clusters
    ranked by `dec` (decay to 25% of the peak). With an odd count the faster half is the nclus // 2 fastest."""
    fast_half = set(np.argsort(np.nan_to_num(dec, nan=99), kind="stable")[: nclus // 2])
    over = (Q < thr) & (F > 1 / nclus)
    pag = [int(c) for c in np.where(over.loc["PAG"])[0] if c in fast_half] if "PAG" in big else []
    vpm = [int(c) for c in np.where(over.loc["VPM"])[0] if c in fast_half] if "VPM" in big else []
    fr = [(r, int(c)) for r in big if is_frontal(r) for c in np.where(over.loc[r])[0] if c not in fast_half]
    return over, fast_half, pag, vpm, fr

Q, F, big = overrep(labels, nclus)
# two rankings of cluster speed: find_decay on the cluster mean of z-scored cells (get_cell_clusters, as
# fig4de) and on the pooled baseline-subtracted trace that table S3 gives for Fig. 4G (verification variant)
rankings = {"cluster mean (fig4de)": np.array([peak_decay(cm[c])[2] for c in range(nclus)]),
            "pooled trace (Fig. 4G)": np.array([peak_decay(cluster_trace(pre_avg, labels == c))[2] for c in range(nclus)])}
print(f"regions with more than 30 clustered cells: {len(big)}; chance fraction 1/{nclus} = {1 / nclus:.3f} (paper 1/12 = 0.083)")
for rk, dec in rankings.items():
    order = np.argsort(np.nan_to_num(dec, nan=99), kind="stable")
    print(f"\nranking by {rk}, fastest first:", [(int(c), round(float(dec[c]), 2)) for c in order])
f31 = []
for thr in (0.05, 0.01):
    over = (Q < thr) & (F > 1 / nclus)
    print(f"\nFDR threshold {thr}: significantly overrepresented (region: clusters)")
    for r in big:
        cl = [int(c) for c in np.where(over.loc[r])[0]]
        if cl:
            print(f"  {r}: {cl}" + ("  frontal/limbic" if is_frontal(r) else ""))
    for rk, dec in rankings.items():
        _, fast_half, pag, vpm, fr = figs31b_test(Q, F, big, nclus, dec, thr)
        print(f"  [{rk}] faster half {sorted(int(c) for c in fast_half)}: PAG in a faster-half cluster {pag}; VPM {vpm}; "
              f"frontal/limbic in a slower-half cluster {fr}")
        f31.append({"FDR": thr, "ranking": rk, "PAG fast": bool(pag), "VPM fast": bool(vpm), "frontal slow": bool(fr),
                    "holds": bool(pag) and bool(vpm) and bool(fr)})
print()
print(pd.DataFrame(f31).to_string(index=False))

regions with more than 30 clustered cells: 54; chance fraction 1/11 = 0.091 (paper 1/12 = 0.083)

ranking by cluster mean (fig4de), fastest first: [(9, 0.1), (6, 0.18), (7, 0.32), (2, 0.49), (8, 0.55), (3, 0.58), (10, 0.86), (4, 1.16), (5, 1.24), (1, 1.35), (0, 1.49)]

ranking by pooled trace (Fig. 4G), fastest first: [(9, 0.1), (6, 0.16), (10, 0.18), (1, 0.27), (7, 0.29), (0, 0.31), (2, 0.4), (3, 0.48), (8, 0.49), (4, 0.88), (5, 0.96)]

FDR threshold 0.05: significantly overrepresented (region: clusters)
  ACAd5: [2]  frontal/limbic
  ACAv1: [1]  frontal/limbic
  ACAv2/3: [1]  frontal/limbic
  ACAv5: [0]  frontal/limbic
  ACB: [2]
  AId5: [0]  frontal/limbic
  AON: [5]
  CA1: [2]
  CP: [0, 10]
  EPd: [0]
  FF: [9]
  HY: [4]
  ILA1: [5]  frontal/limbic
  LSr: [0, 4]
  LSv: [0]
  MB: [2, 9]
  MD: [5]
  MOp5: [0]
  MOp6a: [0]
  MOs6a: [0]  frontal/limbic
  MRN: [2]
  NPC: [9]
  OLF: [0]
  ORBl5: [0]  frontal/limbic
  ORBl6a: [1]  frontal/limbic
  ORBm5: [5]  frontal/limbic
  ORBvl2/3: [1, 5]  frontal/limbic
  ORBvl5: [3, 5]  frontal/limbic
  OT: [0]
  PAG: [9]
  PIR: [0, 1]
  PL1: [0, 5]  frontal/limbic
  PL2/3: [0]  frontal/limbic
  PL5: [2]  frontal/limbic
  PO: [3, 8]
  PoT: [0, 1]
  RSPv5: [0, 3]
  RSPv6a: [3]
  RT: [9]
  SI: [5]
  SUB: [3]
  VPM: [2, 8]
  VPMpc: [1, 3]
  ZI: [9]
  [cluster mean (fig4de)] faster half [2, 6, 7, 8, 9]: PAG in a faster-half cluster [9]; VPM [2, 8]; frontal/limbic in a slower-half cluster [(np.str_('ACAv1'), 1), (np.str_('ACAv2/3'), 1), (np.str_('ACAv5'), 0), (np.str_('AId5'), 0), (np.str_('ILA1'), 5), (np.str_('MOs6a'), 0), (np.str_('ORBl5'), 0), (np.str_('ORBl6a'), 1), (np.str_('ORBm5'), 5), (np.str_('ORBvl2/3'), 1), (np.str_('ORBvl2/3'), 5), (np.str_('ORBvl5'), 3), (np.str_('ORBvl5'), 5), (np.str_('PL1'), 0), (np.str_('PL1'), 5), (np.str_('PL2/3'), 0)]
  [pooled trace (Fig. 4G)] faster half [1, 6, 7, 9, 10]: PAG in a faster-half cluster [9]; VPM []; frontal/limbic in a slower-half cluster [(np.str_('ACAd5'), 2), (np.str_('ACAv5'), 0), (np.str_('AId5'), 0), (np.str_('ILA1'), 5), (np.str_('MOs6a'), 0), (np.str_('ORBl5'), 0), (np.str_('ORBm5'), 5), (np.str_('ORBvl2/3'), 5), (np.str_('ORBvl5'), 3), (np.str_('ORBvl5'), 5), (np.str_('PL1'), 0), (np.str_('PL1'), 5), (np.str_('PL2/3'), 0), (np.str_('PL5'), 2)]

FDR threshold 0.01: significantly overrepresented (region: clusters)
  ACAd5: [2]  frontal/limbic
  ACAv2/3: [1]  frontal/limbic
  ACAv5: [0]  frontal/limbic
  AId5: [0]  frontal/limbic
  AON: [5]
  CA1: [2]
  CP: [0]
  EPd: [0]
  ILA1: [5]  frontal/limbic
  LSr: [0, 4]
  LSv: [0]
  MB: [2, 9]
  MD: [5]
  MOp5: [0]
  MOp6a: [0]
  MRN: [2]
  NPC: [9]
  OLF: [0]
  ORBl5: [0]  frontal/limbic
  ORBm5: [5]  frontal/limbic
  ORBvl2/3: [1, 5]  frontal/limbic
  ORBvl5: [3, 5]  frontal/limbic
  OT: [0]
  PAG: [9]
  PIR: [0, 1]
  PL2/3: [0]  frontal/limbic
  PO: [3, 8]
  RSPv5: [0]
  RSPv6a: [3]
  RT: [9]
  SI: [5]
  SUB: [3]
  VPM: [8]
  VPMpc: [1, 3]
  ZI: [9]
  [cluster mean (fig4de)] faster half [2, 6, 7, 8, 9]: PAG in a faster-half cluster [9]; VPM [8]; frontal/limbic in a slower-half cluster [(np.str_('ACAv2/3'), 1), (np.str_('ACAv5'), 0), (np.str_('AId5'), 0), (np.str_('ILA1'), 5), (np.str_('ORBl5'), 0), (np.str_('ORBm5'), 5), (np.str_('ORBvl2/3'), 1), (np.str_('ORBvl2/3'), 5), (np.str_('ORBvl5'), 3), (np.str_('ORBvl5'), 5), (np.str_('PL2/3'), 0)]
  [pooled trace (Fig. 4G)] faster half [1, 6, 7, 9, 10]: PAG in a faster-half cluster [9]; VPM []; frontal/limbic in a slower-half cluster [(np.str_('ACAd5'), 2), (np.str_('ACAv5'), 0), (np.str_('AId5'), 0), (np.str_('ILA1'), 5), (np.str_('ORBl5'), 0), (np.str_('ORBm5'), 5), (np.str_('ORBvl2/3'), 5), (np.str_('ORBvl5'), 3), (np.str_('ORBvl5'), 5), (np.str_('PL2/3'), 0)]

 FDR                ranking  PAG fast  VPM fast  frontal slow  holds
0.05  cluster mean (fig4de)      True      True          True   True
0.05 pooled trace (Fig. 4G)      True     False          True  False
0.01  cluster mean (fig4de)      True      True          True   True
0.01 pooled trace (Fig. 4G)      True     False          True  False

def fig4g_table(labels, nclus):
    rows = []
    for c in range(nclus):
        tp = cluster_trace(pre_avg, labels == c)
        tk = cluster_trace(ket_avg, labels == c)
        rows.append({"cluster": c, "cells": int((labels == c).sum()), "pre decay s": peak_decay(tp)[2],
                     "cosine distance": scipy.spatial.distance.cosine(tp, tk),
                     "|diff| 1-2 s": abs(tp[200:300].mean() - tk[200:300].mean())})
    g = pd.DataFrame(rows)
    rho = {meas: scipy.stats.spearmanr(g["pre decay s"], g[meas]) for meas in ("cosine distance", "|diff| 1-2 s")}
    return g, rho

g, rho = fig4g_table(labels, nclus)
print(g.round(3).to_string(index=False))
for meas, r in rho.items():
    print(f"{meas}: Spearman rho = {r.statistic:.3f}, p = {r.pvalue:.4f} (n = {nclus} clusters)")

 cluster  cells  pre decay s  cosine distance  |diff| 1-2 s
       0    963         0.31            1.171         0.464
       1    772         0.27            0.637         0.617
       2    664         0.40            0.060         0.089
       3    657         0.48            0.172         0.145
       4    651         0.88            0.622         0.363
       5    604         0.96            0.139         0.228
       6    590         0.16            0.583         0.420
       7    579         0.29            0.233         0.153
       8    554         0.49            0.024         0.084
       9    491         0.10            0.081         0.165
      10    173         0.18            0.754         0.831
cosine distance: Spearman rho = -0.273, p = 0.4171 (n = 11 clusters)
|diff| 1-2 s: Spearman rho = -0.409, p = 0.2115 (n = 11 clusters)

# Stability of the clustering (added in revision). Leiden's partition depends on its random seed; scanpy's
# default (0) is the main result above, as the code calls sc.tl.leiden without a seed. The neighbor graph is
# the same for every seed. For each Leiden seed: the fig4de counts, the figs31b test under both rankings and
# both FDR thresholds, and the fig4g correlations.
stab = []
for seed in range(10):
    sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")
    lab = ann.obs[f"leiden_{seed}"].astype(int).to_numpy()
    nc = lab.max() + 1
    srt_s = np.sort(np.bincount(lab))[::-1]
    nwithin = sum(abs((srt_s[i] if i < len(srt_s) else 0) - e) <= round(0.1 * e) for i, e in enumerate(expected_sizes))
    cm_s = cluster_means(X, lab)
    pk_dec = [peak_decay(cm_s[c]) for c in range(nc)]
    late = [np.abs(cm_s[c][200:300]).mean() > 0.25 * abs(pk_dec[c][0]) for c in range(nc)]
    rapid = sum((d <= 0.25) for _, _, d in pk_dec if not np.isnan(d))
    pers_early = sum(l and t < 0.25 for l, (_, t, _) in zip(late, pk_dec))
    Qs, Fs, bigs = overrep(lab, nc)
    row = {"seed": seed, "clusters": nc, "sizes within (of 12)": nwithin, "rapid and brief": rapid,
           "persistent": sum(late), "persistent, early peak": pers_early,
           "fig4de holds": abs(nc - 12) <= 1 and nwithin == 12 and rapid >= 1 and sum(late) >= 3 and pers_early >= 2 and sum(late) - pers_early >= 1}
    for rk, dec in (("mean", np.array([d for _, _, d in pk_dec])),
                    ("pooled", np.array([peak_decay(cluster_trace(pre_avg, lab == c))[2] for c in range(nc)]))):
        for thr in (0.05, 0.01):
            _, _, pag, vpm, fr = figs31b_test(Qs, Fs, bigs, nc, dec, thr)
            miss = [n for n, v in (("PAG", pag), ("VPM", vpm), ("frontal", fr)) if not v]
            row[f"S31B {rk} {thr}"] = "holds" if not miss else "no " + "+".join(miss)
    _, rho_s = fig4g_table(lab, nc)
    row["4G rho cosine"] = rho_s["cosine distance"].statistic
    row["4G p cosine"] = rho_s["cosine distance"].pvalue
    row["4G rho |diff|"] = rho_s["|diff| 1-2 s"].statistic
    row["4G p |diff|"] = rho_s["|diff| 1-2 s"].pvalue
    stab.append(row)
stab = pd.DataFrame(stab)
print(stab.round(3).to_string(index=False))
print("\nseeds (of 10) where each test holds:")
print(f"fig4de: {int(stab['fig4de holds'].sum())}; cluster count within 12 +- 1: {int((abs(stab.clusters - 12) <= 1).sum())}")
for c in [c for c in stab.columns if c.startswith("S31B")]:
    print(f"{c}: {int((stab[c] == 'holds').sum())}")
print(f"fig4g: rho > 0 with p < 0.05 in {((stab['4G rho cosine'] > 0) & (stab['4G p cosine'] < 0.05)).sum()} seeds (cosine) and "
      f"{((stab['4G rho |diff|'] > 0) & (stab['4G p |diff|'] < 0.05)).sum()} (|diff|)")

# Added in revision 2: the same pipeline with the PCA run in float64, a change at the level of rounding, and
# the partition from 4 BLAS threads (what an unpinned run gave on a 4-core machine), against the main partition.
from sklearn.metrics import adjusted_rand_score
for lbl, dtype, nthr in (("float64 PCA, 1 thread", np.float64, 1), ("float32 PCA, 4 threads", np.float32, 4)):
    with threadpool_limits(nthr):
        a2 = sc.AnnData(X.astype(dtype))
        sc.tl.pca(a2, svd_solver="arpack")
        sc.pp.neighbors(a2, n_neighbors=30, n_pcs=20)
        sc.tl.leiden(a2, flavor="leidenalg")
    l2 = a2.obs["leiden"].astype(int).to_numpy()
    s2 = np.sort(np.bincount(l2))[::-1]
    nw2 = sum(abs((s2[i] if i < len(s2) else 0) - e) <= round(0.1 * e) for i, e in enumerate(expected_sizes))
    print(f"{lbl}: {l2.max() + 1} clusters, sizes {s2.tolist()}, {nw2} of 12 sizes within tolerance, "
          f"adjusted Rand index against the main partition {adjusted_rand_score(labels, l2):.3f}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

/tmp/ipykernel_2286/2756076075.py:7: UserWarning: The `igraph` implementation of leiden clustering is *orders of magnitude faster*. Set the flavor argument to (and install if needed) 'igraph' to use it.
  sc.tl.leiden(ann, flavor="leidenalg", random_state=seed, key_added=f"leiden_{seed}")

first2_pre = np.concatenate([res[k]["first2_pre"] for k in names])      # (cells, 20, 2, 300)
first2_ket = np.concatenate([res[k]["first2_ket"] for k in names])      # (cells, 10, 2, 300), series 23-32
first2_ket_lit = np.concatenate([res[k]["first2_ket_lit"] for k in names])  # series 22-32
red_regions = np.array([K.reduce_region(r) for r in mod_regions])

def dprime_abs(p):
    """region_time_course calculate_dprime on (..., trials, 300): |mean - baseline mean| / pooled sd."""
    bm = p[..., :100].mean((-1, -2))
    bs = np.nanstd(p[..., :100].reshape(p.shape[:-2] + (-1,)), -1)
    rm = p.mean(-2)
    rs = np.nanstd(p, -2)
    return np.abs(rm - bm[..., None]) / np.sqrt(0.5 * (bs[..., None] ** 2 + rs ** 2))

def rise_decay(p, thr=0.25):
    """region_time_course _find_rise_decay, vectorized over leading axes of p (..., trials, 300)."""
    d = dprime_abs(p)
    d = d - d[..., :100].mean(-1, keepdims=True)
    d = d / d.max(-1, keepdims=True)
    post = d[..., 100:]
    pk = post.argmax(-1)
    above = post > thr
    rise = np.where(above.any(-1), above.argmax(-1) * K.BIN, np.nan)
    idx = np.arange(post.shape[-1])
    after = (post < thr) & (idx >= pk[..., None])
    decay = np.where(after.any(-1), after.argmax(-1) * K.BIN, np.nan)
    bad = pk > 75
    rise[bad] = np.nan
    decay[bad] = np.nan
    return rise, decay

def region_times(first2_p, first2_k, min_cells=30, strict=True, dthr=1.3, nreps=1000, seed=0):
    rng = np.random.default_rng(seed)
    out = {}
    reg_u, reg_n = np.unique(red_regions[red_regions != "NA"], return_counts=True)
    for r, nr in zip(reg_u, reg_n):
        if (nr <= min_cells) if strict else (nr < min_cells):
            continue
        f = red_regions == r
        allp = first2_p[f].mean(0).reshape(-1, 300)
        allp = allp - allp[:, :100].mean(1, keepdims=True)
        peak = np.max(dprime_abs(allp[None])[0][100:])
        if peak < dthr:
            continue
        w = np.array([(f & (labels == c)).sum() for c in range(nclus)], float)
        w[w <= 1] = 0
        w /= w.sum()
        res_r = {"n": int(nr), "peak_dprime": peak}
        for cond, arr in (("pre", first2_p), ("ket", first2_k)):
            data = []
            for c in range(nclus):
                m = f & (labels == c)
                if w[c] == 0:
                    data.append(np.zeros((arr.shape[1] * 2, 300)))
                    continue
                d = np.nanmean(arr[m], 0).reshape(-1, 300)
                data.append(d - d[:, :100].mean(1, keepdims=True))
            data = np.stack(data)                                   # (clusters, trials, 300)
            nt = data.shape[1]
            ix = rng.integers(0, nt, size=(nreps, nt))
            ri = np.empty((nclus, nreps)); de = np.empty((nclus, nreps))
            for c in range(nclus):                                  # (reps, trials, 300) per cluster
                ri[c], de[c] = rise_decay(data[c][ix].astype(np.float32))
            res_r[cond + "_rise"] = np.nansum(w[:, None] * np.where(w[:, None] > 0, ri, 0), 0)
            res_r[cond + "_decay"] = np.nansum(w[:, None] * np.where(w[:, None] > 0, de, 0), 0)
            # verification variant: a cluster whose rise or decay is NaN (peak after 750 ms, or no fall below
            # 25%) adds 0 in the code's nansum; here the weights are renormalized over the clusters with a
            # value, and a replicate with none is NaN. Also the weighted fraction of NaN cluster values.
            for nm, v in (("rise", ri), ("decay", de)):
                ok = np.isfinite(v) & (w[:, None] > 0)
                wsum = (w[:, None] * ok).sum(0)
                res_r[f"{cond}_{nm}_rn"] = np.where(wsum > 0, (w[:, None] * np.where(ok, v, 0)).sum(0) / np.where(wsum > 0, wsum, 1), np.nan)
                res_r[f"{cond}_{nm}_nanfrac"] = 1 - wsum.mean()
        out[r] = res_r
    for catch in ("NA", "STR", "MB", "HY"):          # region_time_course drops these catch-all labels
        out.pop(catch, None)
    return out

# One pass per Ketamine reading at the widest settings (>= 30 cells, peak d' >= 1.3); the code's rule
# (> 30 cells) and the d' threshold of 2 (table S3) are subsets of it, applied below.
rt_all = {kname: region_times(first2_pre, karr, min_cells=30, strict=False, dthr=1.3)
          for kname, karr in (("ket 23-32", first2_ket), ("ket 22-32", first2_ket_lit))}
rt = {}
for dthr in (1.3, 2.0):
    for label, keep in (("code >30", lambda o: o["n"] > 30), ("Methods >=30", lambda o: o["n"] >= 30)):
        for kname, out in rt_all.items():
            rt[(dthr, label, kname)] = {r: o for r, o in out.items() if keep(o) and o["peak_dprime"] >= dthr}
# verification control: Pre has 20 series (40 replicate puffs) and Ketamine 10 (20); Pre from series 11 to 20
# only gives both 20, at the main settings
rt_sub = region_times(first2_pre[:, 10:20], first2_ket, min_cells=30, strict=True, dthr=1.3)
rows = []
for key, out in rt.items():
    med_change = np.array([np.nanmedian(o["ket_rise"] - o["pre_rise"]) for o in out.values()])
    pre_rise = np.array([np.nanmedian(o["pre_rise"]) for o in out.values()])
    rows.append({"d' threshold": key[0], "region rule": key[1], "ketamine series": key[2], "regions": len(out),
                 "median over regions of median rise change (s)": np.median(med_change),
                 "regions with Pre median rise < 0.1 s": int((pre_rise < 0.1).sum()),
                 "within one bin": abs(np.median(med_change)) <= 0.01 + 1e-9})
print(pd.DataFrame(rows).to_string(index=False))
main = rt[(1.3, "code >30", "ket 23-32")]
print(pd.DataFrame({r: {"n": o["n"], "peak d'": round(o["peak_dprime"], 2), "Pre rise": np.nanmedian(o["pre_rise"]),
                        "Ket rise": np.nanmedian(o["ket_rise"]), "Pre decay": np.nanmedian(o["pre_decay"]),
                        "Ket decay": np.nanmedian(o["ket_decay"])} for r, o in main.items()}).T.sort_values("Pre decay").to_string())
print("\nweighted fraction of NaN cluster values (cluster weight given 0 by nansum), main settings, over regions:")
nanf = pd.DataFrame({r: {f"{c} {m}": o[f"{c}_{m}_nanfrac"] for c in ("pre", "ket") for m in ("rise", "decay")} for r, o in main.items()}).T
print(nanf.describe().loc[["mean", "50%", "max"]].round(3).to_string())
print("regions where any NaN fraction exceeds 0.1:", list(nanf.index[(nanf > 0.1).any(axis=1)]))
rows = []
for lab, out, sfx in (("code nansum (main)", main, ""), ("renormalized weights", main, "_rn"), ("Pre series 11-20 only", rt_sub, ""), ("Pre series 11-20, renormalized", rt_sub, "_rn")):
    med_change = np.array([np.nanmedian(o["ket_rise" + sfx] - o["pre_rise" + sfx]) for o in out.values()])
    pre_rise = np.array([np.nanmedian(o["pre_rise" + sfx]) for o in out.values()])
    rows.append({"variant": lab, "regions": len(out), "median over regions of median rise change (s)": np.median(med_change),
                 "regions with Pre median rise < 0.1 s": int((pre_rise < 0.1).sum()), "within one bin": abs(np.median(med_change)) <= 0.01 + 1e-9})
print("\nrise-time variants (threshold 1.3, code rule, Ketamine 23-32):")
print(pd.DataFrame(rows).to_string(index=False))

 d' threshold  region rule ketamine series  regions  median over regions of median rise change (s)  regions with Pre median rise < 0.1 s  within one bin
          1.3     code >30       ket 23-32       35                                       0.000122                                    32            True
          1.3     code >30       ket 22-32       35                                      -0.003015                                    32            True
          1.3 Methods >=30       ket 23-32       35                                       0.000122                                    32            True
          1.3 Methods >=30       ket 22-32       35                                      -0.003015                                    32            True
          2.0     code >30       ket 23-32       31                                       0.000122                                    30            True
          2.0     code >30       ket 22-32       31                                      -0.003015                                    30            True
          2.0 Methods >=30       ket 23-32       31                                       0.000122                                    30            True
          2.0 Methods >=30       ket 22-32       31                                      -0.003015                                    30            True
           n  peak d'  Pre rise  Ket rise  Pre decay  Ket decay
LSv     44.0     1.34  0.032386  0.052273   0.202955   0.298636
RT      37.0     3.71  0.025588  0.025000   0.212941   0.149412
CA3     35.0     2.79  0.031471  0.010735   0.227647   0.143676
MOp     96.0     1.71  0.040691  0.143564   0.257926   0.460426
PAG    138.0     3.45  0.010882  0.008676   0.268529   0.145515
FF      51.0     3.03  0.059216  0.045294   0.281176   0.248431
NPC     84.0     3.50  0.010476  0.013571   0.284107   0.154524
TH      53.0     3.68  0.055577  0.030769   0.300865   0.227308
SCig    42.0     2.35  0.090714  0.051905   0.316905   0.194643
MRN     79.0     3.13  0.021795  0.009936   0.328910   0.195128
PoT     68.0     2.59  0.029104  0.028358   0.347612   0.236940
MOs     74.0     2.25  0.043425  0.052397   0.363151   0.263082
ZI     171.0     5.24  0.021023  0.017281   0.379971   0.174532
ILA     63.0     2.40  0.027581  0.029677   0.390081   0.162661
LP      82.0     2.54  0.082284  0.043272   0.393704   0.246173
CA1     99.0     3.76  0.030102  0.035204   0.401888   0.214745
ORBm    45.0     2.01  0.041136  0.050909   0.420795   0.233182
ACAv   154.0     3.29  0.051071  0.053279   0.433182   0.244383
SUB     47.0     2.64  0.059022  0.090870   0.437283   0.311413
DG-mo   57.0     2.97  0.038333  0.033333   0.458684   0.204035
CP     274.0     3.19  0.048175  0.077591   0.495109   0.460821
RSPv   297.0     3.18  0.049360  0.080067   0.512643   0.382626
AIv    132.0     2.10  0.062423  0.070577   0.513269   0.287269
PO     277.0     4.58  0.027762  0.028809   0.526264   0.376986
ACB    205.0     3.20  0.041317  0.041707   0.527439   0.341780
ORBl   390.0     3.92  0.053987  0.082038   0.532231   0.389192
ACAd   153.0     3.72  0.046732  0.038562   0.542353   0.261242
PL     193.0     3.38  0.034119  0.030570   0.549896   0.225648
VPM    131.0     3.20  0.062385  0.041538   0.557423   0.266962
MD     185.0     3.08  0.053777  0.038723   0.570815   0.246141
LSr    403.0     4.00  0.021216  0.032022   0.585695   0.220583
OLF    181.0     2.37  0.094006  0.089972   0.618011   0.387348
ORBvl  127.0     2.07  0.125952  0.054127   0.643016   0.286627
AON     62.0     1.48  0.123115  0.046885   0.646557   0.429836
PIR    154.0     1.76  0.147143  0.109610   0.751396   0.438766

weighted fraction of NaN cluster values (cluster weight given 0 by nansum), main settings, over regions:
      pre rise  pre decay  ket rise  ket decay
mean     0.089      0.087     0.246      0.246
50%      0.082      0.082     0.201      0.201
max      0.259      0.259     0.580      0.577
regions where any NaN fraction exceeds 0.1: [np.str_('ACAd'), np.str_('ACAv'), np.str_('ACB'), np.str_('AIv'), np.str_('AON'), np.str_('CA1'), np.str_('CA3'), np.str_('CP'), np.str_('DG-mo'), np.str_('FF'), np.str_('ILA'), np.str_('LP'), np.str_('LSr'), np.str_('LSv'), np.str_('MD'), np.str_('MOp'), np.str_('MOs'), np.str_('MRN'), np.str_('NPC'), np.str_('OLF'), np.str_('ORBl'), np.str_('ORBm'), np.str_('ORBvl'), np.str_('PAG'), np.str_('PL'), np.str_('RSPv'), np.str_('RT'), np.str_('SCig'), np.str_('SUB'), np.str_('TH'), np.str_('VPM'), np.str_('ZI')]

rise-time variants (threshold 1.3, code rule, Ketamine 23-32):
                       variant  regions  median over regions of median rise change (s)  regions with Pre median rise < 0.1 s  within one bin
            code nansum (main)       35                                       0.000122                                    32            True
          renormalized weights       35                                       0.006471                                    31            True
         Pre series 11-20 only       36                                      -0.003976                                    35            True
Pre series 11-20, renormalized       36                                       0.004043                                    31            True

rows = []
for key, out in rt.items():
    ch = np.concatenate([o["ket_decay"] - o["pre_decay"] for o in out.values()])
    shorter = sum(np.nanmedian(o["ket_decay"]) < np.nanmedian(o["pre_decay"]) for o in out.values())
    rows.append({"d' threshold": key[0], "region rule": key[1], "ketamine series": key[2], "regions": len(out),
                 "median of pooled bootstrap changes (s)": np.nanmedian(ch),
                 "regions with Ket median decay < Pre": shorter})
print(pd.DataFrame(rows).to_string(index=False))
rows = []
for lab, out, sfx in (("code nansum (main)", main, ""), ("renormalized weights", main, "_rn"), ("Pre series 11-20 only", rt_sub, ""), ("Pre series 11-20, renormalized", rt_sub, "_rn")):
    ch = np.concatenate([o["ket_decay" + sfx] - o["pre_decay" + sfx] for o in out.values()])
    shorter = sum(np.nanmedian(o["ket_decay" + sfx]) < np.nanmedian(o["pre_decay" + sfx]) for o in out.values())
    rows.append({"variant": lab, "regions": len(out), "median of pooled bootstrap changes (s)": np.nanmedian(ch),
                 "regions with Ket median decay < Pre": shorter})
print("\ndecay-time variants (threshold 1.3, code rule, Ketamine 23-32):")
print(pd.DataFrame(rows).to_string(index=False))

 d' threshold  region rule ketamine series  regions  median of pooled bootstrap changes (s)  regions with Ket median decay < Pre
          1.3     code >30       ket 23-32       35                               -0.162383                                   33
          1.3     code >30       ket 22-32       35                               -0.158051                                   33
          1.3 Methods >=30       ket 23-32       35                               -0.162383                                   33
          1.3 Methods >=30       ket 22-32       35                               -0.158051                                   33
          2.0     code >30       ket 23-32       31                               -0.165027                                   31
          2.0     code >30       ket 22-32       31                               -0.161015                                   31
          2.0 Methods >=30       ket 23-32       31                               -0.165027                                   31
          2.0 Methods >=30       ket 22-32       31                               -0.161015                                   31

decay-time variants (threshold 1.3, code rule, Ketamine 23-32):
                       variant  regions  median of pooled bootstrap changes (s)  regions with Ket median decay < Pre
            code nansum (main)       35                               -0.162383                                   33
          renormalized weights       35                               -0.124458                                   29
         Pre series 11-20 only       36                               -0.083488                                   31
Pre series 11-20, renormalized       36                               -0.058447                                   26

FRONTAL_LIMBIC = ["ORB", "ACA", "PL", "ILA", "MOs", "MOp", "AI"]
rows = []
for key, out in rt.items():
    if key[1] != "code >30" or key[2] != "ket 23-32":
        continue
    dec = pd.Series({r: np.nanmedian(o["pre_decay"]) for r, o in out.items()})
    fl = [r for r in dec.index if any(K.under(r, f, paths) for f in FRONTAL_LIMBIC)]
    mb = [r for r in dec.index if K.under(r, "MB", paths)]
    med_all = dec.median()
    rows.append({"d' threshold": key[0], "frontal/limbic regions": ", ".join(fl), "median decay frontal/limbic": dec[fl].median(),
                 "midbrain regions": ", ".join(mb), "median decay midbrain": dec[mb].median(), "median all": med_all,
                 "VPM": dec.get("VPM", np.nan), "MD": dec.get("MD", np.nan)})
print(pd.DataFrame(rows).T.to_string())

                                                                                 0                                                 1
d' threshold                                                                   1.3                                               2.0
frontal/limbic regions       ACAd, ACAv, AIv, ILA, MOp, MOs, ORBl, ORBm, ORBvl, PL  ACAd, ACAv, AIv, ILA, MOs, ORBl, ORBm, ORBvl, PL
median decay frontal/limbic                                               0.473226                                          0.513269
midbrain regions                                               MRN, NPC, PAG, SCig                               MRN, NPC, PAG, SCig
median decay midbrain                                                     0.300506                                          0.300506
median all                                                                0.433182                                          0.433182
VPM                                                                       0.557423                                          0.557423
MD                                                                        0.570815                                          0.570815

# Session readings for the cross-puff and population analyses (Figs. 5 and 6): the paper excludes one
# unnamed animal; each mouse with two sessions is a reading (fig5d-emotion-dimension-accumulates).
multi = sorted({res[k]["subject"] for k in names if sum(res[j]["subject"] == res[k]["subject"] for j in names) > 1})
READINGS = {f"without sub-{m}": [k for k in names if res[k]["subject"] != m and res[k]["n_mod"] >= 20] for m in multi}
for r, ks in READINGS.items():
    print(f"{r}: {len(ks)} sessions from {len({res[k]['subject'] for k in ks})} mice")
TONE_SESSIONS = [k for k in names if "tone" in res[k]]
for r, ks in READINGS.items():
    t = [k for k in ks if k in TONE_SESSIONS]
    print(f"  tone sessions {r}: {len(t)} sessions from {len({res[k]['subject'] for k in t})} mice")

def bh(p):
    return scipy.stats.false_discovery_control(np.asarray(p))

def ci95(x):
    x = np.asarray(x, float)
    m = x.mean(); h = 1.96 * x.std(ddof=0) / np.sqrt(len(x))
    return m, m - h, m + h

without sub-1: 11 sessions from 9 mice
without sub-3: 11 sessions from 9 mice
without sub-5: 11 sessions from 9 mice
  tone sessions without sub-1: 8 sessions from 7 mice
  tone sessions without sub-3: 8 sessions from 7 mice
  tone sessions without sub-5: 7 sessions from 7 mice

def zproj(trace, base="own"):
    z = scipy.stats.zscore(trace)
    return z - z[:100].mean()

rows = []
per_session = {}
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        for puffs in ("first", "all"):
            for ep in ("code", "Methods-literal"):
                kk = ("ket", "post") if ep == "code" else ("ket_lit", "post_lit")
                vals = {c: np.array([zproj(res[k]["cd_proj"][("persist", var, cc)][puffs])[200:300].mean() for k in ks])
                        for c, cc in (("Pre", "pre"), ("Ketamine", kk[0]), ("Post", kk[1]))}
                t1 = scipy.stats.ttest_rel(vals["Pre"], vals["Ketamine"])
                t2 = scipy.stats.ttest_rel(vals["Post"], vals["Ketamine"])
                q = bh([t1.pvalue, t2.pvalue])
                means = [vals[c].mean() for c in ("Pre", "Ketamine", "Post")]
                within = [abs(means[0] - 0.86) <= 0.1, abs(means[1] - 0.02) <= 0.25, abs(means[2] - 0.52) <= 0.41]
                rows.append({"reading": rd, "variance": var, "puffs": puffs, "epochs": ep, "Pre": means[0], "Ket": means[1], "Post": means[2],
                             "t Pre-Ket": t1.statistic, "q Pre-Ket": q[0], "t Post-Ket": t2.statistic, "q Post-Ket": q[1],
                             "means within": all(within), "tests hold": q[0] < 0.05 and q[1] < 0.05 and means[1] < min(means[0], means[2])})
                if var == "unit" and puffs == "first" and ep == "code":
                    per_session[rd] = pd.DataFrame(vals, index=ks)
f6d = pd.DataFrame(rows)
print(f6d.round(3).to_string(index=False))
print("\nper session, Methods variance, first puff, code epochs (reading without sub-5):")
print(per_session["without sub-5"].round(3).to_string())
# verification variant (leakage): weights from the odd Pre series only; Pre is then scored on the even Pre
# series, which the weights never saw (first puff, code epochs)
rows = []
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        vals = {c: np.array([zproj(res[k]["cd_proj_heldout"][(var, cc)])[200:300].mean() for k in ks])
                for c, cc in (("Pre (training series)", "pre_train"), ("Pre (held-out series)", "pre_test"), ("Ketamine", "ket"), ("Post", "post"))}
        t1 = scipy.stats.ttest_rel(vals["Pre (held-out series)"], vals["Ketamine"]); t2 = scipy.stats.ttest_rel(vals["Post"], vals["Ketamine"])
        q = bh([t1.pvalue, t2.pvalue])
        rows.append({"reading": rd, "variance": var, **{c: v.mean() for c, v in vals.items()}, "q held-out Pre-Ket": q[0], "q Post-Ket": q[1]})
print("\nheld-out check, first puff, code epochs:")
print(pd.DataFrame(rows).round(3).to_string(index=False))

      reading variance puffs          epochs   Pre   Ket  Post  t Pre-Ket  q Pre-Ket  t Post-Ket  q Post-Ket  means within  tests hold
without sub-1     unit first            code 0.827 0.055 0.501      5.301      0.001       2.213       0.051          True       False
without sub-1     unit first Methods-literal 0.827 0.156 0.453      5.266      0.001       1.436       0.182          True       False
without sub-1     unit   all            code 0.455 0.076 0.419      5.889      0.000       2.998       0.013         False        True
without sub-1     unit   all Methods-literal 0.455 0.126 0.368      5.090      0.001       2.460       0.034         False        True
without sub-1   pooled first            code 0.818 0.066 0.424      4.907      0.001       1.906       0.086          True       False
without sub-1   pooled first Methods-literal 0.818 0.185 0.394      4.686      0.002       1.207       0.255          True       False
without sub-1   pooled   all            code 0.445 0.089 0.352      4.961      0.001       2.021       0.071         False       False
without sub-1   pooled   all Methods-literal 0.445 0.155 0.319      4.831      0.001       1.628       0.134         False       False
without sub-3     unit first            code 0.813 0.000 0.514      5.320      0.001       3.045       0.012          True        True
without sub-3     unit first Methods-literal 0.813 0.091 0.488      5.413      0.001       2.636       0.025          True        True
without sub-3     unit   all            code 0.450 0.069 0.348      5.766      0.000       2.719       0.022         False        True
without sub-3     unit   all Methods-literal 0.450 0.114 0.304      5.605      0.000       2.124       0.060         False       False
without sub-3   pooled first            code 0.794 0.045 0.397      5.084      0.001       2.585       0.027          True        True
without sub-3   pooled first Methods-literal 0.794 0.144 0.383      4.985      0.001       2.174       0.055          True       False
without sub-3   pooled   all            code 0.435 0.066 0.241      5.383      0.001       1.431       0.183         False       False
without sub-3   pooled   all Methods-literal 0.435 0.124 0.216      6.065      0.000       0.966       0.357         False       False
without sub-5     unit first            code 0.867 0.007 0.511      5.527      0.001       2.490       0.032          True        True
without sub-5     unit first Methods-literal 0.867 0.141 0.423      5.218      0.001       1.419       0.186          True       False
without sub-5     unit   all            code 0.454 0.087 0.460      5.605      0.000       3.538       0.005         False        True
without sub-5     unit   all Methods-literal 0.454 0.156 0.413      4.780      0.001       2.601       0.026         False        True
without sub-5   pooled first            code 0.855 0.019 0.512      5.643      0.000       2.613       0.026          True        True
without sub-5   pooled first Methods-literal 0.855 0.172 0.423      5.095      0.001       1.337       0.211          True       False
without sub-5   pooled   all            code 0.441 0.052 0.394      5.891      0.000       2.737       0.021         False        True
without sub-5   pooled   all Methods-literal 0.441 0.143 0.347      5.063      0.001       1.840       0.096         False       False

per session, Methods variance, first puff, code epochs (reading without sub-5):
                      Pre  Ketamine   Post
sub-0_ses-20230901  1.206     0.554  1.331
sub-1_ses-20230218  0.909    -0.534  0.657
sub-1_ses-20230226  0.809    -0.002  0.979
sub-2_ses-20230902  0.608     0.274  0.457
sub-3_ses-20230218  1.071    -0.116  1.599
sub-3_ses-20230225  0.796     0.181 -0.101
sub-4_ses-20230907  0.777     0.471  1.031
sub-6_ses-20230823  0.734    -0.062  0.445
sub-7_ses-20220831  0.870     0.439  0.034
sub-8_ses-20220228  1.163    -0.854 -1.184
sub-9_ses-20230908  0.599    -0.269  0.378

held-out check, first puff, code epochs:
      reading variance  Pre (training series)  Pre (held-out series)  Ketamine  Post  q held-out Pre-Ket  q Post-Ket
without sub-1     unit                  0.797                  0.814     0.050 0.510               0.002       0.054
without sub-1   pooled                  0.795                  0.796     0.067 0.411               0.005       0.101
without sub-3     unit                  0.805                  0.753    -0.016 0.532               0.002       0.012
without sub-3   pooled                  0.810                  0.715     0.029 0.376               0.006       0.033
without sub-5     unit                  0.827                  0.847    -0.008 0.518               0.001       0.032
without sub-5   pooled                  0.833                  0.824     0.009 0.493               0.001       0.033

rows = []
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        for ref in ("1s", "70ms"):
            tr = {c: np.array([scipy.stats.zscore(res[k]["cd_proj"][("fast", var, ref, cc)]["all"]) for k in ks])
                  for c, cc in (("Pre", "pre"), ("Ketamine", "ket"), ("Post", "post"))}
            vals = {}
            for c, a in tr.items():
                a = a - a[:, :100].mean()           # analysis.py subtracts one mean over sessions and baseline bins
                vals[c] = a[:, 100:107].mean(1)
            t1 = scipy.stats.ttest_rel(vals["Pre"], vals["Ketamine"])
            t2 = scipy.stats.ttest_rel(vals["Post"], vals["Ketamine"])
            q = bh([t1.pvalue, t2.pvalue])
            rows.append({"reading": rd, "variance": var, "reference": ref, "Pre": vals["Pre"].mean(), "Ket": vals["Ketamine"].mean(),
                         "Post": vals["Post"].mean(), "t Pre-Ket": t1.statistic, "q Pre-Ket": q[0], "t Post-Ket": t2.statistic,
                         "q Post-Ket": q[1], "Pre-Ket diff": (vals["Pre"] - vals["Ketamine"]).mean()})
print(pd.DataFrame(rows).round(4).to_string(index=False))

      reading variance reference    Pre    Ket   Post  t Pre-Ket  q Pre-Ket  t Post-Ket  q Post-Ket  Pre-Ket diff
without sub-1     unit        1s 4.0067 3.7224 3.5657     1.1732     0.2679     -1.3252      0.2679        0.2843
without sub-1     unit      70ms 3.9596 3.8488 3.5535     0.4379     0.6708     -2.0175      0.1426        0.1108
without sub-1   pooled        1s 4.2494 4.0557 3.7165     0.6466     0.5324     -2.5450      0.0582        0.1937
without sub-1   pooled      70ms 4.2026 4.0441 3.7140     0.5488     0.5952     -2.3094      0.0871        0.1585
without sub-3     unit        1s 4.0907 3.9166 3.7594     0.7513     0.4698     -1.4016      0.3826        0.1741
without sub-3     unit      70ms 4.0716 4.0552 3.7340     0.0701     0.9455     -2.1949      0.1058        0.0165
without sub-3   pooled        1s 4.3152 4.2962 3.9496     0.0729     0.9433     -2.5210      0.0607        0.0190
without sub-3   pooled      70ms 4.2855 4.2910 3.9382    -0.0220     0.9829     -2.4271      0.0713       -0.0055
without sub-5     unit        1s 4.2251 4.1753 3.9770     0.4045     0.6944     -2.0697      0.1306        0.0498
without sub-5     unit      70ms 4.1955 4.2598 3.9654    -0.4602     0.6552     -3.0308      0.0253       -0.0643
without sub-5   pooled        1s 4.4700 4.4988 4.1785    -0.1413     0.8904     -2.9660      0.0283       -0.0288
without sub-5   pooled      70ms 4.4334 4.4961 4.1614    -0.3263     0.7509     -3.1028      0.0224       -0.0627

for var in ("unit", "pooled"):
    for rd, ks in READINGS.items():
        tops = np.hstack([res[k]["top_units_regions"][var] for k in ks])
        incl = np.hstack([res[k]["regions"] for k in ks])
        top_tl = pd.Series([K.top_level(r, paths) for r in tops]).value_counts()
        inc_tl = pd.Series([K.top_level(r, paths) for r in incl]).value_counts()
        tab = pd.DataFrame({"included units": inc_tl, "top-10% units": top_tl}).fillna(0).astype(int)
        tab["share of top units"] = (tab["top-10% units"] / tab["top-10% units"].sum()).round(3)
        div = [d for d in tab.index if d != "NA" and tab.loc[d, "included units"] >= 30]
        missing = [d for d in div if tab.loc[d, "top-10% units"] == 0]
        print(f"{var} variance, {rd}: divisions with >= 30 included units: {div}; without top units: {missing}; "
              f"largest share {tab.drop('NA', errors='ignore')['share of top units'].max():.3f}")
        if var == "unit" and rd == "without sub-5":
            print(tab.to_string())

unit variance, without sub-1: divisions with >= 30 included units: ['CTXsp', 'HPF', 'HY', 'Isocortex', 'MB', 'OLF', 'PAL', 'STR', 'TH']; without top units: []; largest share 0.254
unit variance, without sub-3: divisions with >= 30 included units: ['HPF', 'HY', 'Isocortex', 'MB', 'OLF', 'PAL', 'STR', 'TH']; without top units: []; largest share 0.241
unit variance, without sub-5: divisions with >= 30 included units: ['CTXsp', 'HPF', 'HY', 'Isocortex', 'MB', 'OLF', 'PAL', 'STR', 'TH']; without top units: []; largest share 0.209
           included units  top-10% units  share of top units
CTXsp                  96             15               0.024
HPF                   263             40               0.064
HY                    334             33               0.053
Isocortex            1778            131               0.209
MB                    711             47               0.075
NA                    507             34               0.054
OLF                   448             91               0.145
PAL                    85             13               0.021
STR                   975            110               0.175
TH                   1025            113               0.180
pooled variance, without sub-1: divisions with >= 30 included units: ['CTXsp', 'HPF', 'HY', 'Isocortex', 'MB', 'OLF', 'PAL', 'STR', 'TH']; without top units: []; largest share 0.318
pooled variance, without sub-3: divisions with >= 30 included units: ['HPF', 'HY', 'Isocortex', 'MB', 'OLF', 'PAL', 'STR', 'TH']; without top units: []; largest share 0.289
pooled variance, without sub-5: divisions with >= 30 included units: ['CTXsp', 'HPF', 'HY', 'Isocortex', 'MB', 'OLF', 'PAL', 'STR', 'TH']; without top units: []; largest share 0.266

FIT_VARIANTS = {"code": {}, "bounded b": {"bounded": True}, "demeaned": {"demean": True}, "demeaned, bounded b": {"demean": True, "bounded": True}}

def degenerate(b):
    """A fit whose timescale is not finite, not positive, or beyond the 2.5 s fit window."""
    b = np.asarray(b, float)
    return ~np.isfinite(b) | (b <= 0) | (b > 2.5)

def timescale_rows(get_sigs, ks, exp_pre, exp_ket, labels):
    """Per-series fits averaged (mean, as the Methods, and median) and one concatenated fit, under each fit
    variant; paired t-test across sessions. get_sigs(k) returns (Pre signals, Ketamine signals)."""
    rows = []
    for fv, kw in FIT_VARIANTS.items():
        a_pre, a_ket, m_pre, m_ket, b_pre, b_ket, ndeg, nfit = [], [], [], [], [], [], 0, 0
        for k in ks:
            ps, kt = get_sigs(k)
            tp = np.array([K.extract_timescale(s, **kw) for s in ps]); tk = np.array([K.extract_timescale(s, **kw) for s in kt])
            ndeg += degenerate(tp).sum() + degenerate(tk).sum(); nfit += len(tp) + len(tk)
            a_pre.append(np.nanmean(tp)); a_ket.append(np.nanmean(tk))
            m_pre.append(np.nanmedian(tp)); m_ket.append(np.nanmedian(tk))
            b_pre.append(K.extract_timescale(ps.reshape(-1), **kw)); b_ket.append(K.extract_timescale(kt.reshape(-1), **kw))
        for lab, p_, k_ in (("(a) per series, mean", a_pre, a_ket), ("(a') per series, median", m_pre, m_ket), ("(b) concatenated", b_pre, b_ket)):
            p_, k_ = np.array(p_), np.array(k_)
            t = scipy.stats.ttest_rel(p_, k_)
            rows.append(dict(labels, fit=fv, combine=lab, Pre=p_.mean(), Ket=k_.mean(), t=t.statistic, p=t.pvalue,
                             holds=bool(t.pvalue < 0.05 and k_.mean() < p_.mean()),
                             within=bool(abs(p_.mean() - exp_pre[0]) <= exp_pre[1] and abs(k_.mean() - exp_ket[0]) <= exp_ket[1]),
                             **{"degenerate per-series fits": f"{ndeg} of {nfit}" if lab.startswith("(a") else
                                f"{int(degenerate(p_).sum() + degenerate(k_).sum())} of {2 * len(p_)}"},
                             **{"Ket range": f"{k_.min():.3f}-{k_.max():.3f}"}))
    return rows

rows = []
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        get = lambda k, var=var: (res[k]["base_sigs"][var]["pre_sigs"][10:20], res[k]["base_sigs"][var]["ket_sigs"])
        rows += timescale_rows(get, ks, (0.45, 0.14), (0.14, 0.12), {"reading": rd, "variance": var})
f6h = pd.DataFrame(rows)
print("Pre series 11-20 against Ketamine 23-32. 'code' is analysis.py extract_timescale; the other fits are verification variants.")
print(f6h[f6h.fit == "code"].drop(columns="fit").round(4).to_string(index=False))
print("\nvariants: does the test hold (p < 0.05, Ketamine shorter), and are the means within tolerance?")
piv = f6h.assign(res=lambda d: d.holds.map({True: "holds", False: "fails"}) + d.within.map({True: "", False: " (out)"}) + " p=" + d.p.round(3).astype(str))
print(piv.pivot_table(index=["reading", "variance", "combine"], columns="fit", values="res", aggfunc="first").to_string())

# verification variant: only the series whose 30 s window holds no tone or light onset (cell above lists them)
rows = []
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        def get(k, var=var):
            w = res[k]["stim_in_windows"]
            okp = np.array([(w["tone"]["30 s before series"][s] + w["light"]["30 s before series"][s]) == 0 for s in range(10, 20)])
            okk = np.array([(w["tone"]["30 s before series"][s] + w["light"]["30 s before series"][s]) == 0 for s in range(22, 32)])
            return res[k]["base_sigs"][var]["pre_sigs"][10:20][okp], res[k]["base_sigs"][var]["ket_sigs"][okk]
        rows += [r for r in timescale_rows(get, ks, (0.45, 0.14), (0.14, 0.12), {"reading": rd, "variance": var}) if r["fit"] == "code"]
print("\nwindows without tone or light onsets only (code fit):")
print(pd.DataFrame(rows).drop(columns="fit").round(4).to_string(index=False))

Pre series 11-20 against Ketamine 23-32. 'code' is analysis.py extract_timescale; the other fits are verification variants.
      reading variance                 combine    Pre     Ket       t      p  holds  within degenerate per-series fits     Ket range
without sub-1     unit    (a) per series, mean 0.4515  0.3439  0.6040 0.5593  False   False                   4 of 220   0.019-1.710
without sub-1     unit (a') per series, median 0.2955  0.1267  2.2226 0.0505  False   False                   4 of 220   0.020-0.842
without sub-1     unit        (b) concatenated 0.3445  0.1202  3.0708 0.0118   True    True                    0 of 22   0.019-0.794
without sub-1   pooled    (a) per series, mean 0.5244 77.9188 -0.9978 0.3419  False   False                   5 of 220 0.031-853.339
without sub-1   pooled (a') per series, median 0.3012  0.2585  0.3284 0.7494  False   False                   5 of 220   0.028-1.024
without sub-1   pooled        (b) concatenated 0.3925  0.2083  1.8098 0.1004  False    True                    0 of 22   0.019-0.766
without sub-3     unit    (a) per series, mean 0.4614  0.3286  0.7413 0.4756  False   False                   4 of 220   0.019-1.710
without sub-3     unit (a') per series, median 0.3171  0.1268  2.3977 0.0375   True    True                   4 of 220   0.020-0.842
without sub-3     unit        (b) concatenated 0.3579  0.1207  3.1494 0.0103   True    True                    0 of 22   0.019-0.794
without sub-3   pooled    (a) per series, mean 0.5436 77.9087 -0.9974 0.3421  False   False                   5 of 220 0.035-853.339
without sub-3   pooled (a') per series, median 0.3352  0.2520  0.6012 0.5611  False    True                   5 of 220   0.029-1.024
without sub-3   pooled        (b) concatenated 0.4107  0.1876  2.1865 0.0536  False    True                    0 of 22   0.019-0.766
without sub-5     unit    (a) per series, mean 0.4837  0.2003  3.1881 0.0097   True    True                   3 of 220   0.026-1.049
without sub-5     unit (a') per series, median 0.3391  0.1253  2.7348 0.0210   True    True                   3 of 220   0.025-0.842
without sub-5     unit        (b) concatenated 0.3882  0.1188  3.6955 0.0041   True    True                    0 of 22   0.025-0.794
without sub-5   pooled    (a) per series, mean 0.5571  0.2213  2.9419 0.0147   True    True                   2 of 220   0.031-0.945
without sub-5   pooled (a') per series, median 0.3575  0.1455  2.4990 0.0315   True    True                   2 of 220   0.028-0.761
without sub-5   pooled        (b) concatenated 0.4170  0.1608  2.8857 0.0162   True    True                    0 of 22   0.030-0.766

variants: does the test hold (p < 0.05, Ketamine shorter), and are the means within tolerance?
fit                                                       bounded b                 code             demeaned  demeaned, bounded b
reading       variance combine                                                                                                    
without sub-1 pooled   (a') per series, median  fails (out) p=0.749  fails (out) p=0.749  fails (out) p=0.091  fails (out) p=0.091
                       (a) per series, mean     fails (out) p=0.542  fails (out) p=0.342  fails (out) p=0.342        fails p=0.126
                       (b) concatenated                 fails p=0.1          fails p=0.1        fails p=0.098        fails p=0.098
              unit     (a') per series, median   fails (out) p=0.05   fails (out) p=0.05  holds (out) p=0.047  holds (out) p=0.047
                       (a) per series, mean           fails p=0.077  fails (out) p=0.559  fails (out) p=0.425        fails p=0.078
                       (b) concatenated               holds p=0.012        holds p=0.012        holds p=0.012        holds p=0.012
without sub-3 pooled   (a') per series, median        fails p=0.561        fails p=0.561        fails p=0.061        fails p=0.061
                       (a) per series, mean     fails (out) p=0.437  fails (out) p=0.342  fails (out) p=0.342        fails p=0.096
                       (b) concatenated               fails p=0.054        fails p=0.054        fails p=0.054        fails p=0.054
              unit     (a') per series, median        holds p=0.037        holds p=0.037        holds p=0.028        holds p=0.028
                       (a) per series, mean            holds p=0.05  fails (out) p=0.476  fails (out) p=0.436         holds p=0.05
                       (b) concatenated                holds p=0.01         holds p=0.01        holds p=0.011        holds p=0.011
without sub-5 pooled   (a') per series, median        holds p=0.032        holds p=0.031        holds p=0.022        holds p=0.022
                       (a) per series, mean           holds p=0.029        holds p=0.015        holds p=0.012        holds p=0.024
                       (b) concatenated               holds p=0.016        holds p=0.016        holds p=0.015        holds p=0.015
              unit     (a') per series, median        holds p=0.021        holds p=0.021        holds p=0.019        holds p=0.019
                       (a) per series, mean           holds p=0.024         holds p=0.01        holds p=0.008        holds p=0.024
                       (b) concatenated               holds p=0.004        holds p=0.004        holds p=0.004        holds p=0.004

windows without tone or light onsets only (code fit):
      reading variance                 combine    Pre    Ket      t      p  holds  within degenerate per-series fits   Ket range
without sub-1     unit    (a) per series, mean 0.4707 0.1875 2.7550 0.0203   True    True                   3 of 156 0.019-1.049
without sub-1     unit (a') per series, median 0.3494 0.1247 2.0337 0.0694  False    True                   3 of 156 0.019-0.842
without sub-1     unit        (b) concatenated 0.3653 0.1200 3.0105 0.0131   True    True                    0 of 22 0.019-0.794
without sub-1   pooled    (a) per series, mean 0.5618 0.3059 1.4828 0.1689  False   False                   3 of 156 0.032-0.945
without sub-1   pooled (a') per series, median 0.3580 0.1614 1.5818 0.1448  False    True                   3 of 156 0.024-0.761
without sub-1   pooled        (b) concatenated 0.3960 0.1810 2.1246 0.0596  False    True                    0 of 22 0.020-0.766
without sub-3     unit    (a) per series, mean 0.4764 0.1720 3.0575 0.0121   True    True                   3 of 156 0.019-1.049
without sub-3     unit (a') per series, median 0.3591 0.1263 2.1070 0.0614  False    True                   3 of 156 0.019-0.842
without sub-3     unit        (b) concatenated 0.3704 0.1206 3.0589 0.0121   True    True                    0 of 22 0.019-0.794
without sub-3   pooled    (a) per series, mean 0.5771 0.2969 1.6200 0.1363  False   False                   3 of 156 0.032-0.945
without sub-3   pooled (a') per series, median 0.3825 0.1580 1.7602 0.1089  False    True                   3 of 156 0.024-0.761
without sub-3   pooled        (b) concatenated 0.4059 0.1604 2.4800 0.0325   True    True                    0 of 22 0.020-0.766
without sub-5     unit    (a) per series, mean 0.5106 0.1870 3.2207 0.0092   True    True                   3 of 164 0.026-1.049
without sub-5     unit (a') per series, median 0.3918 0.1260 2.4168 0.0363   True    True                   3 of 164 0.023-0.842
without sub-5     unit        (b) concatenated 0.4090 0.1204 3.6030 0.0048   True    True                    0 of 22 0.025-0.794
without sub-5   pooled    (a) per series, mean 0.6100 0.2075 2.9219 0.0153   True   False                   2 of 164 0.032-0.945
without sub-5   pooled (a') per series, median 0.4122 0.1407 2.1942 0.0530  False    True                   2 of 164 0.030-0.761
without sub-5   pooled        (b) concatenated 0.4438 0.1587 2.9232 0.0152   True    True                    0 of 22 0.032-0.766

rows = []
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        for lab, sname, meas in (("(a) Methods", "methods", "per_series_mean_r"), ("(b) code", "code", "concat_median_abs_r")):
            pre_v, ket_v, nu = [], [], []
            for k in ks:
                c = res[k]["coupling"][(var, sname)]
                f = (lambda d: np.nanmean(d[meas])) if meas == "per_series_mean_r" else (lambda d: d[meas])
                pre_v.append(f(c["pre"])); ket_v.append(f(c["ket"])); nu.append(c["n_units"])
            pre_v, ket_v = np.array(pre_v), np.array(ket_v)
            ok = ~(np.isnan(pre_v) | np.isnan(ket_v))
            t = scipy.stats.ttest_rel(pre_v[ok], ket_v[ok])
            ratio = ket_v[ok] / pre_v[ok]
            rows.append({"reading": rd, "variance": var, "measure": lab, "Pre": pre_v.mean(), "Ket": ket_v.mean(),
                         "Ket/Pre mean": ratio.mean(), "ratio range": f"{ratio.min():.2f}-{ratio.max():.2f}", "t": t.statistic, "p": t.pvalue,
                         "neurons per session": f"{min(nu)}-{max(nu)}", "sessions without a value": int((~ok).sum())})
print(pd.DataFrame(rows).round(4).to_string(index=False))
# verification variant: reading (a) with only the series whose 30 s window holds no tone or light onset
def clean(k, series):
    w = res[k]["stim_in_windows"]
    return np.array([(w["tone"]["30 s before series"][s] + w["light"]["30 s before series"][s]) == 0 for s in series])

rows = []
for rd, ks in READINGS.items():
    for var in ("unit", "pooled"):
        pre_v, ket_v = [], []
        for k in ks:
            c = res[k]["coupling"][(var, "methods")]
            pre_v.append(np.nanmean(c["pre"]["per_series_mean_r"][clean(k, range(0, 10))]))
            ket_v.append(np.nanmean(c["ket"]["per_series_mean_r"][clean(k, range(22, 32))]))
        pre_v, ket_v = np.array(pre_v), np.array(ket_v)
        t = scipy.stats.ttest_rel(pre_v, ket_v)
        rows.append({"reading": rd, "variance": var, "measure": "(a) Methods, clean windows", "Ket/Pre mean": (ket_v / pre_v).mean(), "t": t.statistic, "p": t.pvalue})
print("\nwindows without tone or light onsets only (Pre series 1, 2, 3, 10 and Ketamine 23, 25, 26, 27, 31 in the sessions with those stimuli):")
print(pd.DataFrame(rows).round(4).to_string(index=False))

      reading variance     measure    Pre    Ket  Ket/Pre mean ratio range      t      p neurons per session  sessions without a value
without sub-1     unit (a) Methods 0.0236 0.0173        0.7457   0.12-1.69 2.3962 0.0376                8-77                         0
without sub-1     unit    (b) code 0.0349    NaN        0.7372   0.36-1.28 2.5557 0.0309                3-51                         1
without sub-1   pooled (a) Methods 0.0405 0.0202        0.5559  -0.06-1.40 3.2463 0.0088                8-77                         0
without sub-1   pooled    (b) code 0.0428    NaN        0.5711   0.32-0.89 4.3856 0.0018                7-65                         1
without sub-3     unit (a) Methods 0.0158 0.0115        0.7491   0.12-1.69 1.9056 0.0858                8-77                         0
without sub-3     unit    (b) code 0.0273    NaN        0.7564   0.36-1.28 2.3289 0.0448                3-51                         1
without sub-3   pooled (a) Methods 0.0310 0.0122        0.5115  -0.06-1.40 3.0119 0.0131                8-77                         0
without sub-3   pooled    (b) code 0.0369    NaN        0.5864   0.27-0.89 4.2114 0.0023                7-65                         1
without sub-5     unit (a) Methods 0.0236 0.0180        0.8170   0.35-1.69 2.1435 0.0577               32-77                         0
without sub-5     unit    (b) code 0.0382 0.0247        0.7422   0.36-1.28 2.7360 0.0210               17-51                         0
without sub-5   pooled (a) Methods 0.0352 0.0198        0.5861   0.14-1.40 3.7374 0.0039               32-77                         0
without sub-5   pooled    (b) code 0.0369 0.0188        0.5644   0.27-0.89 4.4947 0.0012               23-65                         0

windows without tone or light onsets only (Pre series 1, 2, 3, 10 and Ketamine 23, 25, 26, 27, 31 in the sessions with those stimuli):
      reading variance                    measure  Ket/Pre mean      t      p
without sub-1     unit (a) Methods, clean windows        0.7301 2.5771 0.0276
without sub-1   pooled (a) Methods, clean windows        0.5310 3.0504 0.0122
without sub-3     unit (a) Methods, clean windows        0.7184 2.2402 0.0490
without sub-3   pooled (a) Methods, clean windows        0.4806 2.8557 0.0171
without sub-5     unit (a) Methods, clean windows        0.7858 2.4353 0.0351
without sub-5   pooled (a) Methods, clean windows        0.5520 4.1778 0.0019

def emo_signal(o, cond, scale):
    """Projection over the series (8 puff windows, -1 to 23 s after the first puff) and the post-series window
    (from 1 s before the last puff, 50 s), baseline-subtracted with the 250 ms before the first puff."""
    f = o["l1_factor"] if scale == "l1" else 1.0
    s8 = o[cond]["s8"] * f
    post = o[cond]["post"] * f
    b = s8[75:100].mean()
    return s8 - b, post - b

def emo_checks(ks, group="all", tw="methods", var="unit", scale="unit"):
    d = {}
    for cond in ("pre", "ket"):
        s8 = []; post = []
        for k in ks:
            o = res[k]["emotion"].get((group, tw, var))
            if o is None:
                continue
            a, b = emo_signal(o, cond, scale)
            s8.append(a); post.append(b)
        d[cond] = (np.array(s8), np.array(post))
    s8p, postp = d["pre"]; s8k, postk = d["ket"]
    n = len(s8p)
    v8 = s8p[:, 7 * 300 + 75:7 * 300 + 100].mean(1)
    v2 = s8p[:, 300 + 75:300 + 100].mean(1)
    v10 = postp[:, 1100:1200].mean(1)          # 10 to 11 s after the last puff onset
    vend = postp[:, 4500:5000].mean(1)         # last 5 s of the 45 s after the series
    k8 = s8k[:, 7 * 300 + 75:7 * 300 + 100].mean(1)
    c8, c2, c10, cend, ck8 = ci95(v8), ci95(v2), ci95(v10), ci95(vend), ci95(k8)
    cd = ci95(v8 - v2)
    # "pre8>pre2" is the criterion's wording, the interval before puff 8 above the value before puff 2 (revised
    # after verification; the first version compared the means, kept as "pre8>pre2 (means)"); the paired
    # difference interval is a further variant. Only the criterion's wording enters "all hold".
    ok = {"pre8>0": c8[1] > 0, "pre8>pre2": c8[1] > c2[0], "10-11s>0": c10[1] > 0,
          "end includes 0": cend[1] <= 0 <= cend[2], "ket8 includes 0 or below pre": (ck8[1] <= 0 <= ck8[2]) or ck8[2] < c8[0]}
    extra = {"pre8>pre2 (means)": c8[0] > c2[0], "pre8-pre2 paired > 0": cd[1] > 0}
    return n, {"before puff 2": c2, "before puff 8": c8, "10-11 s after": c10, "last 5 s": cend, "ket before puff 8": ck8,
               "puff 8 minus puff 2": cd}, dict(ok, **extra)

CRITERION_CHECKS = ("pre8>0", "pre8>pre2", "10-11s>0", "end includes 0", "ket8 includes 0 or below pre")
rows = []
for rd, ks in READINGS.items():
    for tw in ("methods", "code"):
        for var in ("unit", "pooled"):
            for scale in ("unit", "l1"):
                n, c, ok = emo_checks(ks, "all", tw, var, scale)
                row = {"reading": rd, "target": tw, "variance": var, "scale": scale, "n": n}
                row.update({k: f"{v[0]:.3f} [{v[1]:.3f}, {v[2]:.3f}]" for k, v in c.items()})
                row.update(ok); row["all hold"] = all(v for k, v in ok.items() if k in CRITERION_CHECKS)
                rows.append(row)
f5d = pd.DataFrame(rows)
print(f5d.to_string(index=False))
print("\nchecks that hold, of", len(f5d), "readings:")
print(f5d[list(CRITERION_CHECKS) + ["pre8>pre2 (means)", "pre8-pre2 paired > 0", "all hold"]].sum().to_string())

# verification variant: without the two sessions whose 50 s post-series windows hold light onsets (sanity cell)
LIGHT_POST = [k for k in names if res[k]["stim_in_windows"]["light"]["50 s after series"][:20].sum() > 0]
rows = []
for rd, ks in READINGS.items():
    n, c, ok = emo_checks([k for k in ks if k not in LIGHT_POST], "all", "methods", "unit", "unit")
    rows.append({"reading": rd, "sessions": n, "last 5 s": f"{c['last 5 s'][0]:.3f} [{c['last 5 s'][1]:.3f}, {c['last 5 s'][2]:.3f}]",
                 **{k: v for k, v in ok.items() if k in CRITERION_CHECKS}})
print(f"\nwithout {LIGHT_POST} (Methods window, per-unit variance, unit vector):")
print(pd.DataFrame(rows).to_string(index=False))

      reading  target variance scale  n        before puff 2        before puff 8        10-11 s after             last 5 s      ket before puff 8  puff 8 minus puff 2  pre8>0  pre8>pre2  10-11s>0  end includes 0  ket8 includes 0 or below pre  pre8>pre2 (means)  pre8-pre2 paired > 0  all hold
without sub-1 methods     unit  unit 11 0.297 [0.080, 0.514] 0.928 [0.621, 1.235] 0.687 [0.470, 0.904] 0.251 [0.155, 0.347]  0.020 [-0.154, 0.194] 0.632 [0.427, 0.836]    True       True      True           False                          True               True                  True     False
without sub-1 methods     unit    l1 11 0.019 [0.005, 0.032] 0.059 [0.043, 0.075] 0.043 [0.031, 0.055] 0.016 [0.011, 0.021]  0.002 [-0.010, 0.014] 0.041 [0.030, 0.051]    True       True      True           False                          True               True                  True     False
without sub-1 methods   pooled  unit 11 0.461 [0.159, 0.763] 1.410 [0.974, 1.846] 1.057 [0.729, 1.385] 0.384 [0.240, 0.528]  0.040 [-0.074, 0.155] 0.949 [0.671, 1.227]    True       True      True           False                          True               True                  True     False
without sub-1 methods   pooled    l1 11 0.028 [0.009, 0.046] 0.086 [0.064, 0.108] 0.063 [0.046, 0.080] 0.023 [0.016, 0.030]  0.003 [-0.004, 0.010] 0.058 [0.044, 0.072]    True       True      True           False                          True               True                  True     False
without sub-1    code     unit  unit 11 0.525 [0.261, 0.789] 1.077 [0.737, 1.418] 0.638 [0.441, 0.835] 0.235 [0.137, 0.333] -0.054 [-0.198, 0.091] 0.553 [0.322, 0.783]    True       True      True           False                          True               True                  True     False
without sub-1    code     unit    l1 11 0.033 [0.018, 0.049] 0.068 [0.050, 0.086] 0.040 [0.029, 0.050] 0.014 [0.009, 0.020] -0.003 [-0.013, 0.007] 0.035 [0.022, 0.047]    True       True      True           False                          True               True                  True     False
without sub-1    code   pooled  unit 11 0.727 [0.386, 1.068] 1.570 [1.104, 2.035] 0.937 [0.656, 1.218] 0.346 [0.210, 0.482]  0.004 [-0.076, 0.084] 0.842 [0.558, 1.127]    True       True      True           False                          True               True                  True     False
without sub-1    code   pooled    l1 11 0.045 [0.025, 0.066] 0.098 [0.073, 0.123] 0.057 [0.042, 0.072] 0.020 [0.012, 0.028]  0.001 [-0.004, 0.007] 0.053 [0.037, 0.068]    True       True      True           False                          True               True                  True     False
without sub-3 methods     unit  unit 11 0.313 [0.106, 0.521] 1.054 [0.782, 1.326] 0.786 [0.589, 0.983] 0.302 [0.199, 0.405] -0.121 [-0.295, 0.052] 0.741 [0.520, 0.962]    True       True      True           False                          True               True                  True     False
without sub-3 methods     unit    l1 11 0.020 [0.007, 0.033] 0.067 [0.055, 0.080] 0.050 [0.040, 0.060] 0.020 [0.014, 0.025] -0.008 [-0.020, 0.004] 0.048 [0.036, 0.060]    True       True      True           False                          True               True                  True     False
without sub-3 methods   pooled  unit 11 0.465 [0.165, 0.764] 1.620 [1.270, 1.969] 1.208 [0.928, 1.488] 0.474 [0.331, 0.618] -0.078 [-0.225, 0.070] 1.155 [0.871, 1.439]    True       True      True           False                          True               True                  True     False
without sub-3 methods   pooled    l1 11 0.028 [0.010, 0.046] 0.098 [0.083, 0.113] 0.072 [0.058, 0.086] 0.028 [0.020, 0.036] -0.004 [-0.013, 0.005] 0.070 [0.057, 0.084]    True       True      True           False                          True               True                  True     False
without sub-3    code     unit  unit 11 0.517 [0.244, 0.789] 1.216 [0.940, 1.492] 0.709 [0.539, 0.878] 0.261 [0.163, 0.360] -0.112 [-0.242, 0.018] 0.699 [0.448, 0.950]    True       True      True           False                          True               True                  True     False
without sub-3    code     unit    l1 11 0.033 [0.017, 0.049] 0.078 [0.066, 0.090] 0.045 [0.036, 0.053] 0.017 [0.011, 0.022] -0.007 [-0.016, 0.002] 0.045 [0.031, 0.058]    True       True      True           False                          True               True                  True     False
without sub-3    code   pooled  unit 11 0.709 [0.349, 1.069] 1.783 [1.427, 2.138] 1.040 [0.801, 1.279] 0.386 [0.249, 0.523] -0.034 [-0.118, 0.049] 1.074 [0.744, 1.403]    True       True      True           False                          True               True                  True     False
without sub-3    code   pooled    l1 11 0.044 [0.023, 0.066] 0.111 [0.095, 0.127] 0.063 [0.051, 0.075] 0.023 [0.015, 0.031] -0.001 [-0.007, 0.005] 0.066 [0.050, 0.082]    True       True      True           False                          True               True                  True     False
without sub-5 methods     unit  unit 11 0.255 [0.037, 0.473] 1.033 [0.767, 1.300] 0.790 [0.609, 0.971] 0.303 [0.203, 0.403] -0.018 [-0.203, 0.167] 0.778 [0.580, 0.977]    True       True      True           False                          True               True                  True     False
without sub-5 methods     unit    l1 11 0.015 [0.002, 0.028] 0.063 [0.049, 0.077] 0.048 [0.039, 0.058] 0.019 [0.013, 0.025] -0.001 [-0.013, 0.011] 0.048 [0.036, 0.060]    True       True      True           False                          True               True                  True     False
without sub-5 methods   pooled  unit 11 0.368 [0.070, 0.666] 1.508 [1.115, 1.901] 1.178 [0.901, 1.456] 0.450 [0.296, 0.603] -0.040 [-0.201, 0.122] 1.140 [0.842, 1.439]    True       True      True           False                          True               True                  True     False
without sub-5 methods   pooled    l1 11 0.021 [0.004, 0.038] 0.087 [0.066, 0.108] 0.068 [0.053, 0.083] 0.026 [0.017, 0.035] -0.002 [-0.012, 0.008] 0.066 [0.049, 0.082]    True       True      True           False                          True               True                  True     False
without sub-5    code     unit  unit 11 0.444 [0.171, 0.717] 1.155 [0.856, 1.454] 0.695 [0.530, 0.861] 0.260 [0.160, 0.360] -0.041 [-0.160, 0.078] 0.711 [0.457, 0.966]    True       True      True           False                          True               True                  True     False
without sub-5    code     unit    l1 11 0.026 [0.011, 0.042] 0.069 [0.053, 0.085] 0.042 [0.033, 0.051] 0.016 [0.010, 0.022] -0.003 [-0.010, 0.005] 0.043 [0.028, 0.058]    True       True      True           False                          True               True                  True     False
without sub-5    code   pooled  unit 11 0.568 [0.223, 0.913] 1.616 [1.187, 2.045] 0.991 [0.748, 1.234] 0.368 [0.237, 0.500] -0.020 [-0.106, 0.066] 1.048 [0.693, 1.402]    True       True      True           False                          True               True                  True     False
without sub-5    code   pooled    l1 11 0.033 [0.013, 0.052] 0.093 [0.069, 0.116] 0.057 [0.044, 0.071] 0.021 [0.014, 0.029] -0.001 [-0.006, 0.004] 0.060 [0.041, 0.080]    True       True      True           False                          True               True                  True     False

checks that hold, of 24 readings:
pre8>0                          24
pre8>pre2                       24
10-11s>0                        24
end includes 0                   0
ket8 includes 0 or below pre    24
pre8>pre2 (means)               24
pre8-pre2 paired > 0            24
all hold                         0

without ['sub-1_ses-20230226', 'sub-3_ses-20230225'] (Methods window, per-unit variance, unit vector):
      reading  sessions             last 5 s  pre8>0  pre8>pre2  10-11s>0  end includes 0  ket8 includes 0 or below pre
without sub-1        10 0.265 [0.163, 0.367]    True       True      True           False                          True
without sub-3        10 0.305 [0.191, 0.419]    True       True      True           False                          True
without sub-5         9 0.327 [0.213, 0.441]    True       True      True           False                          True

def first_puff(ks, tw, var, scale):
    """analysis.py plot_affective_firstpuff_persist: signal from 0.25 s before the first puff, minus its first 250 ms."""
    out = {}
    for cond in ("pre", "ket"):
        e, l = [], []
        for k in ks:
            o = res[k]["emotion"][("all", tw, var)]
            f = o["l1_factor"] if scale == "l1" else 1.0
            sig = np.concatenate([o[cond]["s8"][75:], o[cond]["post"][400:5000]]) * f
            sig = sig - sig[:25].mean()
            e.append(sig[25:50].mean()); l.append(sig[225:325].mean())
        out[cond] = (np.array(e), np.array(l))
    return out

rows5e = []
for rd, ks in READINGS.items():
    for tw in ("methods", "code"):
        for var in ("unit", "pooled"):
            for scale in ("unit", "l1"):
                fp = first_puff(ks, tw, var, scale)
                te = scipy.stats.ttest_rel(fp["pre"][0], fp["ket"][0])
                tl = scipy.stats.ttest_rel(fp["pre"][1], fp["ket"][1])
                q = bh([te.pvalue, tl.pvalue])
                rows5e.append({"reading": rd, "target": tw, "variance": var, "scale": scale,
                               "early Pre": fp["pre"][0].mean(), "early Ket": fp["ket"][0].mean(), "t early": te.statistic, "q early": q[0],
                               "late Pre": fp["pre"][1].mean(), "late Ket": fp["ket"][1].mean(), "t late": tl.statistic, "q late": q[1],
                               "late within": abs(fp["pre"][1].mean() - 0.546) <= 0.34 and abs(fp["ket"][1].mean() - 0.005) <= 0.084,
                               "early within": abs(fp["pre"][0].mean() - 0.601) <= 0.55 and abs(fp["ket"][0].mean() - 0.332) <= 0.27})
f5e = pd.DataFrame(rows5e)
print(f5e.round(4).to_string(index=False))

      reading  target variance scale  early Pre  early Ket  t early  q early  late Pre  late Ket  t late  q late  late within  early within
without sub-1 methods     unit  unit     0.8706     0.4406   1.9219   0.1094    0.3108    0.0403  1.7573  0.1094         True          True
without sub-1 methods     unit    l1     0.0624     0.0412   1.3841   0.1964    0.0206    0.0032  1.9896  0.1493        False         False
without sub-1 methods   pooled  unit     0.4869     0.3133   1.1119   0.2922    0.4886   -0.0011  2.7199  0.0431         True          True
without sub-1 methods   pooled    l1     0.0346     0.0189   1.3887   0.1951    0.0303   -0.0001  3.0130  0.0261        False         False
without sub-1    code     unit  unit     1.0307     0.4790   2.5120   0.0308    0.4948    0.0451  3.1728  0.0199         True          True
without sub-1    code     unit    l1     0.0669     0.0409   1.4784   0.1701    0.0323    0.0035  3.4743  0.0120        False         False
without sub-1    code   pooled  unit     0.7657     0.3199   2.1299   0.0590    0.7094    0.0124  3.9334  0.0056         True          True
without sub-1    code   pooled    l1     0.0491     0.0185   2.4361   0.0351    0.0454    0.0011  4.1750  0.0038        False         False
without sub-3 methods     unit  unit     0.7985     0.3627   1.7385   0.1128    0.3585    0.0406  2.1257  0.1128         True          True
without sub-3 methods     unit    l1     0.0581     0.0356   1.3079   0.2202    0.0235    0.0030  2.4141  0.0729        False         False
without sub-3 methods   pooled  unit     0.3669     0.2380   0.6086   0.5564    0.5365    0.0279  2.8493  0.0345         True          True
without sub-3 methods   pooled    l1     0.0271     0.0137   0.9449   0.3670    0.0329    0.0015  3.1318  0.0213        False         False
without sub-3    code     unit  unit     0.9851     0.4617   2.2504   0.0481    0.5072    0.0420  3.3253  0.0154         True          True
without sub-3    code     unit    l1     0.0639     0.0394   1.3536   0.2057    0.0334    0.0032  3.7195  0.0080        False         False
without sub-3    code   pooled  unit     0.7056     0.3154   1.6266   0.1349    0.7073    0.0431  3.3396  0.0150         True          True
without sub-3    code   pooled    l1     0.0457     0.0178   1.9467   0.0802    0.0457    0.0028  3.6568  0.0088        False         False
without sub-5 methods     unit  unit     0.8040     0.3293   1.9800   0.1049    0.3015    0.0345  1.7833  0.1049         True          True
without sub-5 methods     unit    l1     0.0487     0.0193   2.1116   0.0830    0.0185    0.0026  1.9262  0.0830        False         False
without sub-5 methods   pooled  unit     0.3275     0.2580   0.3370   0.7431    0.4284   -0.0088  2.4170  0.0725         True          True
without sub-5 methods   pooled    l1     0.0199     0.0148   0.4203   0.6832    0.0247   -0.0008  2.5209  0.0607        False         False
without sub-5    code     unit  unit     0.9589     0.3981   2.6898   0.0227    0.4107    0.0275  2.7859  0.0227         True          True
without sub-5    code     unit    l1     0.0560     0.0224   2.9820   0.0138    0.0243    0.0018  2.9915  0.0138        False         False
without sub-5    code   pooled  unit     0.6081     0.3078   1.3079   0.2202    0.5336    0.0052  2.7636  0.0400         True          True
without sub-5    code   pooled    l1     0.0347     0.0170   1.4005   0.1916    0.0305   -0.0000  2.8427  0.0349        False         False

print(f5e[["reading", "target", "variance", "scale", "early Pre", "early Ket", "t early", "q early", "early within"]].round(4).to_string(index=False))

      reading  target variance scale  early Pre  early Ket  t early  q early  early within
without sub-1 methods     unit  unit     0.8706     0.4406   1.9219   0.1094          True
without sub-1 methods     unit    l1     0.0624     0.0412   1.3841   0.1964         False
without sub-1 methods   pooled  unit     0.4869     0.3133   1.1119   0.2922          True
without sub-1 methods   pooled    l1     0.0346     0.0189   1.3887   0.1951         False
without sub-1    code     unit  unit     1.0307     0.4790   2.5120   0.0308          True
without sub-1    code     unit    l1     0.0669     0.0409   1.4784   0.1701         False
without sub-1    code   pooled  unit     0.7657     0.3199   2.1299   0.0590          True
without sub-1    code   pooled    l1     0.0491     0.0185   2.4361   0.0351         False
without sub-3 methods     unit  unit     0.7985     0.3627   1.7385   0.1128          True
without sub-3 methods     unit    l1     0.0581     0.0356   1.3079   0.2202         False
without sub-3 methods   pooled  unit     0.3669     0.2380   0.6086   0.5564          True
without sub-3 methods   pooled    l1     0.0271     0.0137   0.9449   0.3670         False
without sub-3    code     unit  unit     0.9851     0.4617   2.2504   0.0481          True
without sub-3    code     unit    l1     0.0639     0.0394   1.3536   0.2057         False
without sub-3    code   pooled  unit     0.7056     0.3154   1.6266   0.1349          True
without sub-3    code   pooled    l1     0.0457     0.0178   1.9467   0.0802         False
without sub-5 methods     unit  unit     0.8040     0.3293   1.9800   0.1049          True
without sub-5 methods     unit    l1     0.0487     0.0193   2.1116   0.0830         False
without sub-5 methods   pooled  unit     0.3275     0.2580   0.3370   0.7431          True
without sub-5 methods   pooled    l1     0.0199     0.0148   0.4203   0.6832         False
without sub-5    code     unit  unit     0.9589     0.3981   2.6898   0.0227          True
without sub-5    code     unit    l1     0.0560     0.0224   2.9820   0.0138         False
without sub-5    code   pooled  unit     0.6081     0.3078   1.3079   0.2202          True
without sub-5    code   pooled    l1     0.0347     0.0170   1.4005   0.1916         False

rows_r = []
for rd, ks in READINGS.items():
    for tw in ("methods", "code"):
        for area in K.TOP_LEVEL:
            qual = [k for k in ks if (area, tw, "unit") in res[k]["emotion"]]
            if len(qual) < 2:
                continue
            n, c, ok = emo_checks(qual, area, tw, "unit", "unit")
            rows_r.append({"reading": rd, "target": tw, "area": area, "sessions": n,
                           "pre before puff 8": f"{c['before puff 8'][0]:.3f} [{c['before puff 8'][1]:.3f}, {c['before puff 8'][2]:.3f}]",
                           "pre before puff 2": f"{c['before puff 2'][0]:.3f}",
                           "ket before puff 8": f"{c['ket before puff 8'][0]:.3f} [{c['ket before puff 8'][1]:.3f}, {c['ket before puff 8'][2]:.3f}]",
                           "meets": ok["pre8>0"] and ok["pre8>pre2"] and ok["ket8 includes 0 or below pre"],
                           "meets (means)": ok["pre8>0"] and ok["pre8>pre2 (means)"] and ok["ket8 includes 0 or below pre"],
                           "meets (paired)": ok["pre8>0"] and ok["pre8-pre2 paired > 0"] and ok["ket8 includes 0 or below pre"]})
f32 = pd.DataFrame(rows_r)
print(f32.to_string(index=False))
print("\nareas meeting the conditions, under the criterion's wording (interval before puff 8 above the value before puff 2),")
print("comparing the means (the first version), and with the paired puff 8 minus puff 2 interval:")
print(f32.groupby(["reading", "target"])[["meets", "meets (means)", "meets (paired)"]].agg(lambda x: f"{x.sum()} of {len(x)}").to_string())
print("areas failing under the criterion's wording:", f32[~f32.meets].groupby("area").size().to_dict())

      reading  target      area  sessions    pre before puff 8 pre before puff 2      ket before puff 8  meets  meets (means)  meets (paired)
without sub-1 methods Isocortex        10 0.607 [0.400, 0.814]             0.211 -0.020 [-0.173, 0.133]   True           True            True
without sub-1 methods       OLF         8 0.189 [0.049, 0.330]             0.051 -0.077 [-0.219, 0.065]  False           True            True
without sub-1 methods       HPF         4 0.234 [0.016, 0.452]             0.119  0.093 [-0.026, 0.212]  False           True           False
without sub-1 methods       STR        10 0.226 [0.075, 0.377]            -0.007  0.042 [-0.027, 0.112]   True           True            True
without sub-1 methods        TH         6 0.735 [0.331, 1.139]             0.320  0.002 [-0.164, 0.169]   True           True            True
without sub-1 methods        HY         6 0.411 [0.153, 0.669]            -0.035 -0.004 [-0.078, 0.071]   True           True            True
without sub-1 methods        MB         5 0.444 [0.166, 0.723]             0.179  0.020 [-0.133, 0.172]  False           True            True
without sub-1    code Isocortex        10 0.720 [0.472, 0.967]             0.376 -0.103 [-0.241, 0.036]   True           True            True
without sub-1    code       OLF         8 0.204 [0.048, 0.360]             0.164 -0.081 [-0.246, 0.083]  False           True           False
without sub-1    code       HPF         4 0.251 [0.156, 0.347]             0.099  0.047 [-0.167, 0.261]   True           True           False
without sub-1    code       STR        10 0.287 [0.121, 0.454]             0.157 -0.018 [-0.100, 0.065]  False           True           False
without sub-1    code        TH         6 0.909 [0.512, 1.307]             0.439  0.052 [-0.122, 0.226]   True           True            True
without sub-1    code        HY         6 0.551 [0.273, 0.829]             0.122 -0.027 [-0.154, 0.101]   True           True            True
without sub-1    code        MB         5 0.500 [0.139, 0.861]             0.193  0.036 [-0.154, 0.227]  False           True            True
without sub-3 methods Isocortex        10 0.641 [0.437, 0.844]             0.268 -0.074 [-0.223, 0.074]   True           True            True
without sub-3 methods       OLF         6 0.262 [0.115, 0.409]             0.099 -0.036 [-0.214, 0.141]   True           True            True
without sub-3 methods       HPF         5 0.208 [0.028, 0.388]             0.077  0.049 [-0.073, 0.171]  False           True           False
without sub-3 methods       STR         9 0.280 [0.116, 0.443]            -0.009 -0.020 [-0.085, 0.046]   True           True            True
without sub-3 methods        TH         7 0.723 [0.368, 1.077]             0.334 -0.134 [-0.328, 0.060]   True           True            True
without sub-3 methods        HY         7 0.381 [0.147, 0.615]            -0.054 -0.014 [-0.081, 0.053]   True           True            True
without sub-3 methods        MB         6 0.527 [0.238, 0.816]             0.139  0.046 [-0.104, 0.196]   True           True            True
without sub-3    code Isocortex        10 0.747 [0.514, 0.981]             0.401 -0.137 [-0.273, 0.000]   True           True            True
without sub-3    code       OLF         6 0.290 [0.135, 0.444]             0.194 -0.129 [-0.331, 0.072]  False           True           False
without sub-3    code       HPF         5 0.241 [0.163, 0.320]             0.080  0.028 [-0.146, 0.202]   True           True           False
without sub-3    code       STR         9 0.333 [0.187, 0.479]             0.185 -0.028 [-0.110, 0.054]   True           True           False
without sub-3    code        TH         7 0.856 [0.507, 1.204]             0.344 -0.033 [-0.185, 0.118]   True           True            True
without sub-3    code        HY         7 0.502 [0.242, 0.762]             0.067  0.022 [-0.091, 0.135]   True           True            True
without sub-3    code        MB         6 0.629 [0.302, 0.957]             0.176  0.028 [-0.141, 0.197]   True           True            True
without sub-5 methods Isocortex        11 0.533 [0.333, 0.733]             0.180  0.008 [-0.096, 0.113]   True           True            True
without sub-5 methods       OLF         8 0.189 [0.049, 0.330]             0.051 -0.077 [-0.219, 0.065]  False           True            True
without sub-5 methods       HPF         5 0.208 [0.028, 0.388]             0.077  0.049 [-0.073, 0.171]  False           True           False
without sub-5 methods       STR        10 0.251 [0.096, 0.407]            -0.038  0.028 [-0.053, 0.109]   True           True            True
without sub-5 methods        TH         8 0.691 [0.375, 1.007]             0.292 -0.104 [-0.283, 0.074]   True           True            True
without sub-5 methods        HY         8 0.356 [0.146, 0.565]            -0.060 -0.022 [-0.082, 0.038]   True           True            True
without sub-5 methods        MB         7 0.503 [0.251, 0.754]             0.125  0.045 [-0.084, 0.174]   True           True            True
without sub-5    code Isocortex        11 0.608 [0.385, 0.831]             0.288 -0.071 [-0.169, 0.027]   True           True            True
without sub-5    code       OLF         8 0.204 [0.048, 0.360]             0.164 -0.081 [-0.246, 0.083]  False           True           False
without sub-5    code       HPF         5 0.241 [0.163, 0.320]             0.080  0.028 [-0.146, 0.202]   True           True           False
without sub-5    code       STR        10 0.282 [0.117, 0.447]             0.115 -0.008 [-0.088, 0.073]   True           True            True
without sub-5    code        TH         8 0.855 [0.550, 1.160]             0.317  0.002 [-0.145, 0.149]   True           True            True
without sub-5    code        HY         8 0.471 [0.237, 0.706]             0.069  0.001 [-0.105, 0.107]   True           True            True
without sub-5    code        MB         7 0.578 [0.283, 0.874]             0.156  0.037 [-0.109, 0.183]   True           True            True

areas meeting the conditions, under the criterion's wording (interval before puff 8 above the value before puff 2),
comparing the means (the first version), and with the paired puff 8 minus puff 2 interval:
                        meets meets (means) meets (paired)
reading       target                                      
without sub-1 code     4 of 7        7 of 7         4 of 7
              methods  4 of 7        7 of 7         6 of 7
without sub-3 code     6 of 7        7 of 7         4 of 7
              methods  6 of 7        7 of 7         6 of 7
without sub-5 code     6 of 7        7 of 7         5 of 7
              methods  5 of 7        7 of 7         6 of 7
areas failing under the criterion's wording: {'HPF': 3, 'MB': 2, 'OLF': 5, 'STR': 1}

def model_data(ks, tw, var, scale):
    segs = []
    for cond in ("pre", "ket"):
        a = []
        for k in ks:
            s8, post = emo_signal(res[k]["emotion"][("all", tw, var)], cond, scale)
            a.append(np.concatenate([s8[75:2400], post[500:5000]]))     # -0.25 to 23 s and 25 to 70 s after the first puff
        segs.append(np.mean(a, 0))
    return np.concatenate(segs)

import time as _time
rows_m = []
fits = {}
t0 = _time.time()
for rd, ks in READINGS.items():
    for tw in ("methods", "code"):
        for var in ("unit", "pooled"):
            for scale in ("unit", "l1"):
                data = model_data(ks, tw, var, scale)
                r = {}
                for kind in ("persist", "saturation", "broadcast", "symmetric"):
                    try:
                        r[kind] = K.compare_model(kind, data)
                    except Exception as e:
                        r[kind] = {"varexp": np.nan, "fit": [np.nan] * 5, "pred": np.full_like(data, np.nan), "null": np.full_like(data, np.nan), "error": str(e)}
                fits[(rd, tw, var, scale)] = (data, r)
                pm = r["persist"]
                rows_m.append({"reading": rd, "target": tw, "variance": var, "scale": scale, "VE persist": pm["varexp"],
                               "S": pm["fit"][0], "tauB": pm["fit"][1], "tauP": pm["fit"][2], "factor Pre": pm["fit"][3], "factor Ket": pm["fit"][4],
                               "tauP Pre": pm["fit"][2] * pm["fit"][3], "tauP Ket": pm["fit"][2] * pm["fit"][4]})
print(f"{len(fits)} fits in {_time.time() - t0:.0f} s")
print(pd.DataFrame(rows_m).round(4).to_string(index=False))

24 fits in 18 s
      reading  target variance scale  VE persist      S   tauB    tauP  factor Pre  factor Ket  tauP Pre  tauP Ket
without sub-1 methods     unit  unit      0.8379 1.4593 0.7842 27.8453      1.0000      0.0124   27.8453    0.3459
without sub-1 methods     unit    l1      0.8125 0.0905 0.7225 27.3776      0.9998      0.0140   27.3719    0.3829
without sub-1 methods   pooled  unit      0.9196 2.6684 1.2785 28.5822      0.9935      0.0084   28.3964    0.2399
without sub-1 methods   pooled    l1      0.9062 0.1636 1.2891 27.5047      0.9999      0.0087   27.5031    0.2390
without sub-1    code     unit  unit      0.8663 1.4480 0.4002 20.8507      1.0000      0.0130   20.8507    0.2701
without sub-1    code     unit    l1      0.8502 0.0904 0.3755 20.4583      0.9999      0.0134   20.4554    0.2744
without sub-1    code   pooled  unit      0.9377 2.1250 0.5497 22.5010      1.0000      0.0084   22.5010    0.1891
without sub-1    code   pooled    l1      0.9253 0.1313 0.5296 21.5266      1.0000      0.0079   21.5266    0.1691
without sub-3 methods     unit  unit      0.8863 1.6526 0.9555 30.4425      1.0000      0.0065   30.4425    0.1964
without sub-3 methods     unit    l1      0.8695 0.1038 0.9138 30.1858      0.9994      0.0068   30.1684    0.2058
without sub-3 methods   pooled  unit      0.9399 3.1155 1.5449 31.6297      0.9992      0.0046   31.6060    0.1449
without sub-3 methods   pooled    l1      0.9303 0.1940 1.6131 30.7823      1.0000      0.0044   30.7816    0.1340
without sub-3    code     unit  unit      0.8997 1.5095 0.4799 23.6423      1.0000      0.0095   23.6423    0.2250
without sub-3    code     unit    l1      0.8891 0.0947 0.4572 23.6213      0.9999      0.0094   23.6187    0.2212
without sub-3    code   pooled  unit      0.9479 2.2766 0.6480 25.1294      1.0000      0.0060   25.1294    0.1507
without sub-3    code   pooled    l1      0.9392 0.1388 0.6124 24.3038      1.0000      0.0050   24.3038    0.1223
without sub-5 methods     unit  unit      0.8655 1.7350 1.0958 32.1677      1.0000      0.0078   32.1677    0.2517
without sub-5 methods     unit    l1      0.8629 0.1097 1.1733 32.3865      0.9992      0.0068   32.3591    0.2201
without sub-5 methods   pooled  unit      0.9286 3.5042 2.0460 33.4923      0.9985      0.0054   33.4408    0.1818
without sub-5 methods   pooled    l1      0.9285 0.2073 2.1135 33.6526      0.9995      0.0054   33.6370    0.1802
without sub-5    code     unit  unit      0.8860 1.4408 0.5075 25.7084      1.0000      0.0092   25.7084    0.2353
without sub-5    code     unit    l1      0.8813 0.0873 0.5383 26.0900      0.9997      0.0081   26.0813    0.2126
without sub-5    code   pooled  unit      0.9468 2.1496 0.7620 27.7854      1.0000      0.0067   27.7850    0.1869
without sub-5    code   pooled    l1      0.9460 0.1255 0.7829 27.8718      1.0000      0.0065   27.8710    0.1801

# The fitted trace of each condition is the series (-0.25 to 23 s after the first puff, 2325 bins) followed by
# 25 to 70 s after it (4500 bins). Variance explained on the Ketamine puff-series segment (criterion) is taken
# over the Ketamine series bins only; the first version used the whole Ketamine half, shown as a variant.
N_SERIES_BINS = 2400 - 75
rows_a = []
for key, (data, r) in fits.items():
    half = len(data) // 2
    for seg, sl in (("Ketamine series", slice(half, half + N_SERIES_BINS)), ("Ketamine half (series and after)", slice(half, None))):
        ve = {}
        for kind, m in r.items():
            ve[kind] = 1 - np.mean((data[sl] - m["pred"][sl]) ** 2) / np.mean((data[sl] - m["null"][sl]) ** 2)
        best = max(ve, key=lambda k: -np.inf if np.isnan(ve[k]) else ve[k])
        srt_ve = sorted(v for v in ve.values() if np.isfinite(v))
        rows_a.append(dict(zip(("reading", "target", "variance", "scale"), key), segment=seg, **{f"VE {k}": v for k, v in ve.items()},
                           best=best, margin=srt_ve[-1] - srt_ve[-2]))
f33 = pd.DataFrame(rows_a)
for seg in ("Ketamine series", "Ketamine half (series and after)"):
    s = f33[f33.segment == seg]
    print(f"\n{seg}:")
    print(s.drop(columns="segment").round(4).to_string(index=False))
    print(f"tauP model best in {(s.best == 'persist').sum()} of {len(s)} readings; best otherwise: {s.loc[s.best != 'persist', 'best'].value_counts().to_dict()}")

Ketamine series:
      reading  target variance scale  VE persist  VE saturation  VE broadcast  VE symmetric       best  margin
without sub-1 methods     unit  unit      0.9030         0.8622        0.8586        0.7616    persist  0.0408
without sub-1 methods     unit    l1      0.8723         0.8207        0.8148        0.7356    persist  0.0516
without sub-1 methods   pooled  unit      0.9665         0.9495        0.9491        0.7696    persist  0.0170
without sub-1 methods   pooled    l1      0.9658         0.9504        0.9500        0.7772    persist  0.0154
without sub-1    code     unit  unit      0.9265         0.8706        0.8622        0.8206    persist  0.0559
without sub-1    code     unit    l1      0.9110         0.8492        0.8387        0.8059    persist  0.0618
without sub-1    code   pooled  unit      0.9635         0.9520        0.9506        0.8457    persist  0.0115
without sub-1    code   pooled    l1      0.9603         0.9582        0.9569        0.8563    persist  0.0021
without sub-3 methods     unit  unit      0.9196         0.8922        0.8914        0.7279    persist  0.0274
without sub-3 methods     unit    l1      0.8972         0.8613        0.8596        0.7122    persist  0.0359
without sub-3 methods   pooled  unit      0.9700         0.9669        0.9667        0.7396    persist  0.0032
without sub-3 methods   pooled    l1      0.9701         0.9692        0.9690        0.7447    persist  0.0009
without sub-3    code     unit  unit      0.9280         0.8862        0.8801        0.8076    persist  0.0418
without sub-3    code     unit    l1      0.9152         0.8703        0.8631        0.7972    persist  0.0449
without sub-3    code   pooled  unit      0.9644         0.9633        0.9618        0.8350    persist  0.0011
without sub-3    code   pooled    l1      0.9617         0.9697        0.9684        0.8437 saturation  0.0012
without sub-5 methods     unit  unit      0.9099         0.8759        0.8751        0.7125    persist  0.0339
without sub-5 methods     unit    l1      0.8992         0.8727        0.8722        0.6994    persist  0.0264
without sub-5 methods   pooled  unit      0.9651         0.9552        0.9545        0.7006    persist  0.0099
without sub-5 methods   pooled    l1      0.9639         0.9542        0.9530        0.6951    persist  0.0097
without sub-5    code     unit  unit      0.9259         0.8659        0.8620        0.7805    persist  0.0601
without sub-5    code     unit    l1      0.9180         0.8666        0.8639        0.7717    persist  0.0514
without sub-5    code   pooled  unit      0.9674         0.9469        0.9461        0.7946    persist  0.0205
without sub-5    code   pooled    l1      0.9672         0.9480        0.9477        0.7901    persist  0.0193
tauP model best in 23 of 24 readings; best otherwise: {'saturation': 1}

Ketamine half (series and after):
      reading  target variance scale  VE persist  VE saturation  VE broadcast  VE symmetric       best  margin
without sub-1 methods     unit  unit      0.8529         0.8385        0.8363        0.5077    persist  0.0144
without sub-1 methods     unit    l1      0.8228         0.8028        0.7993        0.4939    persist  0.0200
without sub-1 methods   pooled  unit      0.9542         0.9468        0.9466        0.5115    persist  0.0075
without sub-1 methods   pooled    l1      0.9515         0.9451        0.9449        0.5185    persist  0.0064
without sub-1    code     unit  unit      0.8769         0.8425        0.8366        0.5453    persist  0.0344
without sub-1    code     unit    l1      0.8595         0.8205        0.8132        0.5330    persist  0.0390
without sub-1    code   pooled  unit      0.9520         0.9453        0.9443        0.5729    persist  0.0068
without sub-1    code   pooled    l1      0.9492         0.9487        0.9478        0.5855    persist  0.0005
without sub-3 methods     unit  unit      0.8946         0.8802        0.8798        0.4493    persist  0.0144
without sub-3 methods     unit    l1      0.8753         0.8562        0.8555        0.4394    persist  0.0192
without sub-3 methods   pooled  unit      0.9655         0.9642        0.9638        0.4786    persist  0.0013
without sub-3 methods   pooled    l1      0.9650         0.9648        0.9644        0.4819    persist  0.0002
without sub-3    code     unit  unit      0.8955         0.8685        0.8649        0.5143    persist  0.0270
without sub-3    code     unit    l1      0.8836         0.8540        0.8499        0.5038    persist  0.0296
without sub-3    code   pooled  unit      0.9565         0.9573        0.9563        0.5638 saturation  0.0008
without sub-3    code   pooled    l1      0.9556         0.9619        0.9611        0.5723 saturation  0.0008
without sub-5 methods     unit  unit      0.8842         0.8686        0.8681        0.4557    persist  0.0156
without sub-5 methods     unit    l1      0.8759         0.8633        0.8630        0.4405    persist  0.0126
without sub-5 methods   pooled  unit      0.9620         0.9575        0.9551        0.4444    persist  0.0045
without sub-5 methods   pooled    l1      0.9608         0.9564        0.9532        0.4372    persist  0.0044
without sub-5    code     unit  unit      0.8883         0.8518        0.8496        0.4892    persist  0.0365
without sub-5    code     unit    l1      0.8798         0.8485        0.8473        0.4682    persist  0.0313
without sub-5    code   pooled  unit      0.9574         0.9457        0.9452        0.5104    persist  0.0117
without sub-5    code   pooled    l1      0.9572         0.9461        0.9457        0.5015    persist  0.0111
tauP model best in 22 of 24 readings; best otherwise: {'saturation': 2}

def tone_vals(ks, split, key, puffs, win, base="global"):
    out = {}
    for c in ("pre", "ket", "post"):
        a = np.array([scipy.stats.zscore(res[k]["tone"][split]["proj"][key + (c,)][puffs]) for k in ks])
        a = a - (a[:, :100].mean() if base == "global" else a[:, :100].mean(1, keepdims=True))
        out[c] = a[:, win].mean(1)
    return out

def tone_tests(v):
    t1 = scipy.stats.ttest_rel(v["pre"], v["ket"]); t2 = scipy.stats.ttest_rel(v["post"], v["ket"])
    q = bh([t1.pvalue, t2.pvalue])
    return t1.statistic, q[0], t2.statistic, q[1]

rows = []
for rd, ks0 in READINGS.items():
    ks = [k for k in ks0 if k in TONE_SESSIONS]
    for split in ("infusion", "code"):
        for var in ("unit", "pooled"):
            for ref in ("1s", "70ms"):
                v = tone_vals(ks, split, ("fast", var, ref), "all", slice(100, 107))
                t1, q1, t2, q2 = tone_tests(v)
                rows.append({"reading": rd, "sessions": len(ks), "split": split, "variance": var, "reference": ref,
                             "Pre": v["pre"].mean(), "Ket": v["ket"].mean(), "Post": v["post"].mean(),
                             "t Pre-Ket": t1, "q Pre-Ket": q1, "t Post-Ket": t2, "q Post-Ket": q2})
print(pd.DataFrame(rows).round(4).to_string(index=False))

      reading  sessions    split variance reference    Pre    Ket   Post  t Pre-Ket  q Pre-Ket  t Post-Ket  q Post-Ket
without sub-1         8 infusion     unit        1s 3.0267 0.9802 0.9053     3.2650     0.0275     -0.2734      0.7924
without sub-1         8 infusion     unit      70ms 2.5895 1.3429 0.9986     2.4572     0.0436     -2.4862      0.0436
without sub-1         8 infusion   pooled        1s 3.3142 1.3907 1.1113     3.0796     0.0357     -1.2871      0.2390
without sub-1         8 infusion   pooled      70ms 2.8318 1.4534 1.0376     2.5473     0.0383     -2.6085      0.0383
without sub-1         8     code     unit        1s 3.0599 1.0477 0.9762     2.8492     0.0494     -0.2563      0.8051
without sub-1         8     code     unit      70ms 2.5776 1.3931 1.0330     1.9813     0.0880     -2.0249      0.0880
without sub-1         8     code   pooled        1s 3.3372 1.4533 1.0795     2.8221     0.0514     -1.6569      0.1415
without sub-1         8     code   pooled      70ms 2.8356 1.5531 1.0433     2.2229     0.0616     -3.4461      0.0215
without sub-3         8 infusion     unit        1s 2.7692 1.2297 1.1575     3.1418     0.0327     -0.3115      0.7645
without sub-3         8 infusion     unit      70ms 2.3507 1.4457 1.1862     2.2142     0.1112     -1.8225      0.1112
without sub-3         8 infusion   pooled        1s 3.0290 1.5701 1.3056     2.6930     0.0619     -1.5260      0.1708
without sub-3         8 infusion   pooled      70ms 2.6166 1.5456 1.1968     2.1477     0.0688     -2.3386      0.0688
without sub-3         8     code     unit        1s 2.7998 1.3238 1.2250     2.5451     0.0768     -0.4103      0.6938
without sub-3         8     code     unit      70ms 2.3494 1.6219 1.2790     1.8468     0.1073     -2.4664      0.0861
without sub-3         8     code   pooled        1s 3.0485 1.7065 1.2724     2.5013     0.0409     -3.3438      0.0247
without sub-3         8     code   pooled      70ms 2.6253 1.7565 1.2938     1.9276     0.0953     -3.8361      0.0128
without sub-5         7 infusion     unit        1s 2.7134 1.3743 1.2160     2.2001     0.1402     -0.5634      0.5936
without sub-5         7 infusion     unit      70ms 2.2838 1.5737 1.2247     1.4766     0.1902     -2.2999      0.1222
without sub-5         7 infusion   pooled        1s 2.7808 1.6709 1.2970     1.9846     0.1564     -1.6199      0.1564
without sub-5         7 infusion   pooled      70ms 2.4107 1.7000 1.2713     1.6407     0.1520     -2.5288      0.0895
without sub-5         7     code     unit        1s 2.7725 1.5570 1.3635     1.8866     0.2163     -0.7203      0.4984
without sub-5         7     code     unit      70ms 2.2689 1.5125 1.3289     1.2128     0.3920     -0.9222      0.3920
without sub-5         7     code   pooled        1s 2.8121 1.6911 1.3852     1.7546     0.2597     -1.1365      0.2991
without sub-5         7     code   pooled      70ms 2.3935 1.6481 1.3056     1.3413     0.2284     -1.9034      0.2113

rows = []
for rd, ks0 in READINGS.items():
    ks = [k for k in ks0 if k in TONE_SESSIONS]
    for split in ("infusion", "code"):
        for var in ("unit", "pooled"):
            v = tone_vals(ks, split, ("persist", var), "all", slice(135, 150), base="own")
            t1, q1, t2, q2 = tone_tests(v)
            m = [v[c].mean() for c in ("pre", "ket", "post")]
            rows.append({"reading": rd, "sessions": len(ks), "split": split, "variance": var, "Pre": m[0], "Ket": m[1], "Post": m[2],
                         "t Pre-Ket": t1, "q Pre-Ket": q1, "t Post-Ket": t2, "q Post-Ket": q2,
                         "holds": q1 < 0.05 and q2 < 0.05 and m[1] < min(m[0], m[2]),
                         "within": abs(m[0] - 1.531) <= 0.16 and abs(m[1] - 0.392) <= 0.35 and abs(m[2] - 1.385) <= 0.66})
print(pd.DataFrame(rows).round(4).to_string(index=False))

      reading  sessions    split variance    Pre    Ket   Post  t Pre-Ket  q Pre-Ket  t Post-Ket  q Post-Ket  holds  within
without sub-1         8 infusion     unit 1.6784 0.3893 1.2816     5.6083     0.0016      2.6757      0.0317   True    True
without sub-1         8 infusion   pooled 1.5766 0.3501 1.1579     4.9560     0.0033      2.4550      0.0438   True    True
without sub-1         8     code     unit 1.6772 0.4851 1.2676     6.9246     0.0005      2.7975      0.0266   True    True
without sub-1         8     code   pooled 1.5867 0.5038 1.1784     6.2918     0.0008      2.3153      0.0538  False    True
without sub-3         8 infusion     unit 1.6389 0.3268 1.0882     5.7476     0.0014      2.7257      0.0295   True    True
without sub-3         8 infusion   pooled 1.5121 0.3131 0.9640     4.7916     0.0040      2.4170      0.0463   True    True
without sub-3         8     code     unit 1.6324 0.3554 1.0817     7.5289     0.0003      2.8879      0.0234   True    True
without sub-3         8     code   pooled 1.5201 0.3875 1.0052     6.7447     0.0005      2.3044      0.0546  False    True
without sub-5         7 infusion     unit 1.6388 0.4673 1.4996     7.8578     0.0004      2.8705      0.0284   True    True
without sub-5         7 infusion   pooled 1.5123 0.3565 1.3650     6.0215     0.0019      3.0380      0.0229   True    True
without sub-5         7     code     unit 1.6389 0.4919 1.4799     8.7250     0.0003      3.1089      0.0209   True    True
without sub-5         7     code   pooled 1.5221 0.4031 1.3760     7.1014     0.0008      3.5287      0.0124   True    True

def half_fall(tr):
    post = tr[100:]
    pk = int(np.argmax(post))
    below = np.where(post[pk:] < 0.5 * post[pk])[0]
    return (pk + below[0]) * K.BIN if len(below) else np.nan

rows = []
for rd, ks0 in READINGS.items():
    ks = [k for k in ks0 if k in TONE_SESSIONS]
    for split in ("infusion", "code"):
        for var in ("unit", "pooled"):
            dp = np.array([res[k]["tone"][split]["dprime_puff"][var] for k in ks])
            dt = np.array([res[k]["tone"][split]["dprime_tone"][var] for k in ks])
            cp, ct = ci95(dp[:, 200:300].mean(1)), ci95(dt[:, 200:300].mean(1))
            rows.append({"reading": rd, "sessions": len(ks), "split": split, "variance": var,
                         "puff half-peak fall (s)": half_fall(dp.mean(0)), "tone half-peak fall (s)": half_fall(dt.mean(0)),
                         "puff d' 1-2 s": f"{cp[0]:.2f} [{cp[1]:.2f}, {cp[2]:.2f}]", "tone d' 1-2 s": f"{ct[0]:.2f} [{ct[1]:.2f}, {ct[2]:.2f}]",
                         "holds": half_fall(dp.mean(0)) > half_fall(dt.mean(0)) and cp[1] > ct[2]})
print(pd.DataFrame(rows).to_string(index=False))

      reading  sessions    split variance  puff half-peak fall (s)  tone half-peak fall (s)     puff d' 1-2 s      tone d' 1-2 s  holds
without sub-1         8 infusion     unit                     0.95                     0.17 1.10 [0.76, 1.44] 0.10 [-0.05, 0.25]   True
without sub-1         8 infusion   pooled                     0.99                     0.17 1.06 [0.74, 1.38] 0.09 [-0.06, 0.24]   True
without sub-1         8     code     unit                     0.95                     0.17 1.10 [0.76, 1.44] 0.13 [-0.02, 0.27]   True
without sub-1         8     code   pooled                     0.99                     0.17 1.06 [0.74, 1.38] 0.12 [-0.02, 0.25]   True
without sub-3         8 infusion     unit                     0.92                     0.16 1.10 [0.76, 1.44]  0.17 [0.05, 0.30]   True
without sub-3         8 infusion   pooled                     0.98                     0.16 1.06 [0.74, 1.38]  0.19 [0.04, 0.33]   True
without sub-3         8     code     unit                     0.92                     0.16 1.10 [0.76, 1.44]  0.19 [0.06, 0.31]   True
without sub-3         8     code   pooled                     0.98                     0.16 1.06 [0.74, 1.38]  0.20 [0.06, 0.34]   True
without sub-5         7 infusion     unit                     0.90                     0.16 1.06 [0.73, 1.39]  0.17 [0.01, 0.34]   True
without sub-5         7 infusion   pooled                     0.97                     0.16 1.04 [0.73, 1.34] 0.19 [-0.00, 0.37]   True
without sub-5         7     code     unit                     0.90                     0.16 1.06 [0.73, 1.39]  0.19 [0.04, 0.35]   True
without sub-5         7     code   pooled                     0.97                     0.16 1.04 [0.73, 1.34]  0.21 [0.04, 0.38]   True

rows = []
for rd, ks0 in READINGS.items():
    ks = [k for k in ks0 if k in TONE_SESSIONS]
    for split in ("infusion", "code"):
        for var in ("unit", "pooled"):
            for pre_lab, pre_sl in (("Pre 5-10 (code)", slice(4, 10)), ("Pre all 10 (table S3)", slice(0, 10))):
                get = lambda k, split=split, var=var, pre_sl=pre_sl: (res[k]["tone"][split]["base_sigs"][var]["pre_sigs"][pre_sl],
                                                                      res[k]["tone"][split]["base_sigs"][var]["post_sigs"][1:6])
                rows += timescale_rows(get, ks, (0.268, 0.073), (0.089, 0.024),
                                       {"reading": rd, "sessions": len(ks), "split": split, "variance": var, "pre series": pre_lab})
f34g = pd.DataFrame(rows)
print("Ketamine is tone series 2 to 6 after the split. 'code' is analysis.py extract_timescale; the other fits are verification variants.")
print(f34g[f34g.fit == "code"].drop(columns="fit").round(4).to_string(index=False))
print("\nvariants: does the test hold (p < 0.05, Ketamine shorter), and are the means within tolerance?")
piv = f34g.assign(res=lambda d: d.holds.map({True: "holds", False: "fails"}) + d.within.map({True: "", False: " (out)"}) + " p=" + d.p.round(3).astype(str))
print(piv.pivot_table(index=["reading", "split", "variance", "pre series", "combine"], columns="fit", values="res", aggfunc="first").to_string())

Ketamine is tone series 2 to 6 after the split. 'code' is analysis.py extract_timescale; the other fits are verification variants.
      reading  sessions    split variance            pre series                 combine    Pre    Ket       t      p  holds  within degenerate per-series fits   Ket range
without sub-1         8 infusion     unit       Pre 5-10 (code)    (a) per series, mean 0.2795 0.1862  1.2379 0.2557  False   False                    0 of 88 0.018-0.557
without sub-1         8 infusion     unit       Pre 5-10 (code) (a') per series, median 0.2342 0.0795  3.6226 0.0085   True    True                    0 of 88 0.019-0.169
without sub-1         8 infusion     unit       Pre 5-10 (code)        (b) concatenated 0.2543 0.1000  2.6439 0.0332   True    True                    0 of 16 0.018-0.279
without sub-1         8 infusion     unit Pre all 10 (table S3)    (a) per series, mean 0.2757 0.1862  1.1878 0.2737  False   False                   0 of 120 0.018-0.557
without sub-1         8 infusion     unit Pre all 10 (table S3) (a') per series, median 0.2181 0.0795  4.6358 0.0024   True    True                   0 of 120 0.019-0.169
without sub-1         8 infusion     unit Pre all 10 (table S3)        (b) concatenated 0.2195 0.1000  2.6828 0.0314   True    True                    0 of 16 0.018-0.279
without sub-1         8 infusion   pooled       Pre 5-10 (code)    (a) per series, mean 0.2813 0.2515  0.2358 0.8203  False   False                    1 of 88 0.023-1.047
without sub-1         8 infusion   pooled       Pre 5-10 (code) (a') per series, median 0.2301 0.1150  2.0933 0.0746  False   False                    1 of 88 0.022-0.373
without sub-1         8 infusion   pooled       Pre 5-10 (code)        (b) concatenated 0.2564 0.1248  1.8006 0.1148  False   False                    0 of 16 0.021-0.476
without sub-1         8 infusion   pooled Pre all 10 (table S3)    (a) per series, mean 0.3059 0.2515  0.4023 0.6995  False   False                   2 of 120 0.023-1.047
without sub-1         8 infusion   pooled Pre all 10 (table S3) (a') per series, median 0.2124 0.1150  2.1349 0.0702  False   False                   2 of 120 0.022-0.373
without sub-1         8 infusion   pooled Pre all 10 (table S3)        (b) concatenated 0.2166 0.1248  1.5829 0.1575  False   False                    0 of 16 0.021-0.476
without sub-1         8     code     unit       Pre 5-10 (code)    (a) per series, mean 0.2556 0.1827  1.2594 0.2483  False   False                    0 of 80 0.018-0.418
without sub-1         8     code     unit       Pre 5-10 (code) (a') per series, median 0.2090 0.0845  4.5937 0.0025   True    True                    0 of 80 0.019-0.170
without sub-1         8     code     unit       Pre 5-10 (code)        (b) concatenated 0.2358 0.0778  4.0015 0.0052   True    True                    0 of 16 0.019-0.134
without sub-1         8     code     unit Pre all 10 (table S3)    (a) per series, mean 0.2634 0.1827  1.2784 0.2418  False   False                   0 of 112 0.018-0.418
without sub-1         8     code     unit Pre all 10 (table S3) (a') per series, median 0.2023 0.0845  4.2730 0.0037   True    True                   0 of 112 0.019-0.170
without sub-1         8     code     unit Pre all 10 (table S3)        (b) concatenated 0.2096 0.0778  4.2591 0.0038   True    True                    0 of 16 0.019-0.134
without sub-1         8     code   pooled       Pre 5-10 (code)    (a) per series, mean 0.2582 0.2750 -0.1497 0.8852  False   False                    1 of 80 0.051-0.962
without sub-1         8     code   pooled       Pre 5-10 (code) (a') per series, median 0.2004 0.0838  4.1789 0.0041   True    True                    1 of 80 0.024-0.152
without sub-1         8     code   pooled       Pre 5-10 (code)        (b) concatenated 0.2392 0.1109  2.5844 0.0362   True    True                    0 of 16 0.021-0.321
without sub-1         8     code   pooled Pre all 10 (table S3)    (a) per series, mean 0.3026 0.2750  0.2172 0.8342  False   False                   2 of 112 0.051-0.962
without sub-1         8     code   pooled Pre all 10 (table S3) (a') per series, median 0.1991 0.0838  4.5898 0.0025   True    True                   2 of 112 0.024-0.152
without sub-1         8     code   pooled Pre all 10 (table S3)        (b) concatenated 0.2077 0.1109  2.4895 0.0416   True    True                    0 of 16 0.021-0.321
without sub-3         8 infusion     unit       Pre 5-10 (code)    (a) per series, mean 0.2659 0.1911  1.0651 0.3222  False   False                    0 of 88 0.018-0.557
without sub-3         8 infusion     unit       Pre 5-10 (code) (a') per series, median 0.2220 0.0841  3.3879 0.0116   True    True                    0 of 88 0.019-0.169
without sub-3         8 infusion     unit       Pre 5-10 (code)        (b) concatenated 0.2401 0.1052  2.4441 0.0445   True    True                    0 of 16 0.018-0.279
without sub-3         8 infusion     unit Pre all 10 (table S3)    (a) per series, mean 0.2707 0.1911  1.1082 0.3044  False   False                   0 of 120 0.018-0.557
without sub-3         8 infusion     unit Pre all 10 (table S3) (a') per series, median 0.2141 0.0841  4.8897 0.0018   True    True                   0 of 120 0.019-0.169
without sub-3         8 infusion     unit Pre all 10 (table S3)        (b) concatenated 0.2106 0.1052  2.6996 0.0307   True    True                    0 of 16 0.018-0.279
without sub-3         8 infusion   pooled       Pre 5-10 (code)    (a) per series, mean 0.2679 0.2541  0.1122 0.9138  False   False                    1 of 88 0.023-1.047
without sub-3         8 infusion   pooled       Pre 5-10 (code) (a') per series, median 0.2161 0.1169  1.8555 0.1059  False   False                    1 of 88 0.022-0.373
without sub-3         8 infusion   pooled       Pre 5-10 (code)        (b) concatenated 0.2441 0.1283  1.6229 0.1486  False   False                    0 of 16 0.021-0.476
without sub-3         8 infusion   pooled Pre all 10 (table S3)    (a) per series, mean 0.3063 0.2541  0.3869 0.7103  False   False                   2 of 120 0.023-1.047
without sub-3         8 infusion   pooled Pre all 10 (table S3) (a') per series, median 0.2082 0.1169  2.0670 0.0776  False   False                   2 of 120 0.022-0.373
without sub-3         8 infusion   pooled Pre all 10 (table S3)        (b) concatenated 0.2119 0.1283  1.5046 0.1761  False   False                    0 of 16 0.021-0.476
without sub-3         8     code     unit       Pre 5-10 (code)    (a) per series, mean 0.2425 0.1877  1.0945 0.3100  False   False                    0 of 80 0.018-0.418
without sub-3         8     code     unit       Pre 5-10 (code) (a') per series, median 0.2012 0.0893  4.4994 0.0028   True    True                    0 of 80 0.019-0.170
without sub-3         8     code     unit       Pre 5-10 (code)        (b) concatenated 0.2219 0.0828  3.9438 0.0056   True    True                    0 of 16 0.019-0.134
without sub-3         8     code     unit Pre all 10 (table S3)    (a) per series, mean 0.2603 0.1877  1.2223 0.2611  False   False                   0 of 112 0.018-0.418
without sub-3         8     code     unit Pre all 10 (table S3) (a') per series, median 0.2106 0.0893  4.1460 0.0043   True    True                   0 of 112 0.019-0.170
without sub-3         8     code     unit Pre all 10 (table S3)        (b) concatenated 0.2019 0.0828  4.9049 0.0017   True    True                    0 of 16 0.019-0.134
without sub-3         8     code   pooled       Pre 5-10 (code)    (a) per series, mean 0.2456 0.2780 -0.3002 0.7728  False   False                    1 of 80 0.075-0.962
without sub-3         8     code   pooled       Pre 5-10 (code) (a') per series, median 0.1918 0.0865  3.7987 0.0067   True   False                    1 of 80 0.024-0.152
without sub-3         8     code   pooled       Pre 5-10 (code)        (b) concatenated 0.2284 0.1146  2.4188 0.0462   True   False                    0 of 16 0.021-0.321
without sub-3         8     code   pooled Pre all 10 (table S3)    (a) per series, mean 0.3052 0.2780  0.2144 0.8364  False   False                   2 of 112 0.075-0.962
without sub-3         8     code   pooled Pre all 10 (table S3) (a') per series, median 0.2067 0.0865  4.4884 0.0028   True    True                   2 of 112 0.024-0.152
without sub-3         8     code   pooled Pre all 10 (table S3)        (b) concatenated 0.2053 0.1146  2.5138 0.0402   True   False                    0 of 16 0.021-0.321
without sub-5         7 infusion     unit       Pre 5-10 (code)    (a) per series, mean 0.3134 0.1404  3.8005 0.0090   True   False                    0 of 77 0.029-0.382
without sub-5         7 infusion     unit       Pre 5-10 (code) (a') per series, median 0.2537 0.0902  3.6170 0.0111   True    True                    0 of 77 0.029-0.169
without sub-5         7 infusion     unit       Pre 5-10 (code)        (b) concatenated 0.2834 0.0818  3.9344 0.0077   True    True                    0 of 14 0.028-0.167
without sub-5         7 infusion     unit Pre all 10 (table S3)    (a) per series, mean 0.3209 0.1404  3.8783 0.0082   True   False                   0 of 105 0.029-0.382
without sub-5         7 infusion     unit Pre all 10 (table S3) (a') per series, median 0.2463 0.0902  5.4483 0.0016   True    True                   0 of 105 0.029-0.169
without sub-5         7 infusion     unit Pre all 10 (table S3)        (b) concatenated 0.2523 0.0818  5.2857 0.0019   True    True                    0 of 14 0.028-0.167
without sub-5         7 infusion   pooled       Pre 5-10 (code)    (a) per series, mean 0.3075 0.1454  4.1451 0.0060   True   False                    0 of 77 0.055-0.399
without sub-5         7 infusion   pooled       Pre 5-10 (code) (a') per series, median 0.2410 0.0855  3.3078 0.0163   True    True                    0 of 77 0.058-0.128
without sub-5         7 infusion   pooled       Pre 5-10 (code)        (b) concatenated 0.2809 0.0826  3.7990 0.0090   True    True                    0 of 14 0.049-0.142
without sub-5         7 infusion   pooled Pre all 10 (table S3)    (a) per series, mean 0.3511 0.1454  3.2372 0.0177   True   False                   1 of 105 0.055-0.399
without sub-5         7 infusion   pooled Pre all 10 (table S3) (a') per series, median 0.2344 0.0855  5.2703 0.0019   True    True                   1 of 105 0.058-0.128
without sub-5         7 infusion   pooled Pre all 10 (table S3)        (b) concatenated 0.2457 0.0826  6.1216 0.0009   True    True                    0 of 14 0.049-0.142
without sub-5         7     code     unit       Pre 5-10 (code)    (a) per series, mean 0.2904 0.1565  3.1042 0.0210   True   False                    0 of 70 0.030-0.406
without sub-5         7     code     unit       Pre 5-10 (code) (a') per series, median 0.2234 0.0915  4.7699 0.0031   True    True                    0 of 70 0.030-0.170
without sub-5         7     code     unit       Pre 5-10 (code)        (b) concatenated 0.2626 0.0857  4.4812 0.0042   True    True                    0 of 14 0.031-0.134
without sub-5         7     code     unit Pre all 10 (table S3)    (a) per series, mean 0.3098 0.1565  2.7990 0.0312   True   False                    0 of 98 0.030-0.406
without sub-5         7     code     unit Pre all 10 (table S3) (a') per series, median 0.2488 0.0915  5.8721 0.0011   True    True                    0 of 98 0.030-0.170
without sub-5         7     code     unit Pre all 10 (table S3)        (b) concatenated 0.2436 0.0857  5.2846 0.0019   True    True                    0 of 14 0.031-0.134
without sub-5         7     code   pooled       Pre 5-10 (code)    (a) per series, mean 0.2857 0.1614  3.6355 0.0109   True   False                    0 of 70 0.051-0.373
without sub-5         7     code   pooled       Pre 5-10 (code) (a') per series, median 0.2085 0.0891  3.9462 0.0076   True    True                    0 of 70 0.051-0.152
without sub-5         7     code   pooled       Pre 5-10 (code)        (b) concatenated 0.2630 0.0885  4.5187 0.0040   True    True                    0 of 14 0.044-0.141
without sub-5         7     code   pooled Pre all 10 (table S3)    (a) per series, mean 0.3503 0.1614  2.6635 0.0373   True   False                    1 of 98 0.051-0.373
without sub-5         7     code   pooled Pre all 10 (table S3) (a') per series, median 0.2340 0.0891  5.6493 0.0013   True    True                    1 of 98 0.051-0.152
without sub-5         7     code   pooled Pre all 10 (table S3)        (b) concatenated 0.2387 0.0885  6.2473 0.0008   True    True                    0 of 14 0.044-0.141

variants: does the test hold (p < 0.05, Ketamine shorter), and are the means within tolerance?
fit                                                                                      bounded b                 code             demeaned  demeaned, bounded b
reading       split    variance pre series            combine                                                                                                    
without sub-1 code     pooled   Pre 5-10 (code)       (a') per series, median        holds p=0.004        holds p=0.004        holds p=0.006        holds p=0.006
                                                      (a) per series, mean     fails (out) p=0.539  fails (out) p=0.885  fails (out) p=0.779  fails (out) p=0.358
                                                      (b) concatenated               holds p=0.036        holds p=0.036        holds p=0.036        holds p=0.036
                                Pre all 10 (table S3) (a') per series, median        holds p=0.002        holds p=0.003        holds p=0.004        holds p=0.004
                                                      (a) per series, mean     fails (out) p=0.426  fails (out) p=0.834  fails (out) p=0.596  fails (out) p=0.268
                                                      (b) concatenated               holds p=0.042        holds p=0.042        holds p=0.039        holds p=0.039
                       unit     Pre 5-10 (code)       (a') per series, median        holds p=0.003        holds p=0.003        holds p=0.006        holds p=0.006
                                                      (a) per series, mean     fails (out) p=0.248  fails (out) p=0.248  fails (out) p=0.114  fails (out) p=0.114
                                                      (b) concatenated               holds p=0.005        holds p=0.005        holds p=0.006        holds p=0.006
                                Pre all 10 (table S3) (a') per series, median        holds p=0.004        holds p=0.004        holds p=0.006        holds p=0.006
                                                      (a) per series, mean     fails (out) p=0.242  fails (out) p=0.242  fails (out) p=0.088  fails (out) p=0.088
                                                      (b) concatenated               holds p=0.004        holds p=0.004        holds p=0.004        holds p=0.004
              infusion pooled   Pre 5-10 (code)       (a') per series, median  fails (out) p=0.075  fails (out) p=0.075        holds p=0.007        holds p=0.007
                                                      (a) per series, mean     fails (out) p=0.419   fails (out) p=0.82  fails (out) p=0.702  fails (out) p=0.282
                                                      (b) concatenated         fails (out) p=0.115  fails (out) p=0.115  fails (out) p=0.119  fails (out) p=0.119
                                Pre all 10 (table S3) (a') per series, median   fails (out) p=0.07   fails (out) p=0.07        holds p=0.002        holds p=0.002
                                                      (a) per series, mean     fails (out) p=0.399  fails (out) p=0.699  fails (out) p=0.609  fails (out) p=0.253
                                                      (b) concatenated         fails (out) p=0.157  fails (out) p=0.157  fails (out) p=0.163  fails (out) p=0.163
                       unit     Pre 5-10 (code)       (a') per series, median        holds p=0.008        holds p=0.008  holds (out) p=0.009  holds (out) p=0.009
                                                      (a) per series, mean     fails (out) p=0.256  fails (out) p=0.256  fails (out) p=0.165  fails (out) p=0.165
                                                      (b) concatenated               holds p=0.033        holds p=0.033        holds p=0.034        holds p=0.034
                                Pre all 10 (table S3) (a') per series, median        holds p=0.002        holds p=0.002  holds (out) p=0.003  holds (out) p=0.003
                                                      (a) per series, mean     fails (out) p=0.274  fails (out) p=0.274  fails (out) p=0.148  fails (out) p=0.148
                                                      (b) concatenated               holds p=0.031        holds p=0.031        holds p=0.033        holds p=0.033
without sub-3 code     pooled   Pre 5-10 (code)       (a') per series, median  holds (out) p=0.007  holds (out) p=0.007   holds (out) p=0.01   holds (out) p=0.01
                                                      (a) per series, mean     fails (out) p=0.659  fails (out) p=0.773  fails (out) p=0.871  fails (out) p=0.434
                                                      (b) concatenated         holds (out) p=0.046  holds (out) p=0.046        holds p=0.046        holds p=0.046
                                Pre all 10 (table S3) (a') per series, median        holds p=0.003        holds p=0.003        holds p=0.004        holds p=0.004
                                                      (a) per series, mean     fails (out) p=0.428  fails (out) p=0.836    fails (out) p=0.6   fails (out) p=0.27
                                                      (b) concatenated          holds (out) p=0.04   holds (out) p=0.04        holds p=0.038        holds p=0.038
                       unit     Pre 5-10 (code)       (a') per series, median        holds p=0.003        holds p=0.003  holds (out) p=0.008  holds (out) p=0.008
                                                      (a) per series, mean      fails (out) p=0.31   fails (out) p=0.31  fails (out) p=0.142  fails (out) p=0.142
                                                      (b) concatenated               holds p=0.006        holds p=0.006        holds p=0.006        holds p=0.006
                                Pre all 10 (table S3) (a') per series, median        holds p=0.004        holds p=0.004        holds p=0.006        holds p=0.006
                                                      (a) per series, mean     fails (out) p=0.261  fails (out) p=0.261  fails (out) p=0.093  fails (out) p=0.093
                                                      (b) concatenated               holds p=0.002        holds p=0.002        holds p=0.002        holds p=0.002
              infusion pooled   Pre 5-10 (code)       (a') per series, median  fails (out) p=0.106  fails (out) p=0.106        holds p=0.012        holds p=0.012
                                                      (a) per series, mean     fails (out) p=0.506  fails (out) p=0.914  fails (out) p=0.795  fails (out) p=0.349
                                                      (b) concatenated         fails (out) p=0.149  fails (out) p=0.149  fails (out) p=0.154  fails (out) p=0.154
                                Pre all 10 (table S3) (a') per series, median  fails (out) p=0.078  fails (out) p=0.078        holds p=0.002        holds p=0.002
                                                      (a) per series, mean     fails (out) p=0.409   fails (out) p=0.71  fails (out) p=0.624  fails (out) p=0.262
                                                      (b) concatenated         fails (out) p=0.176  fails (out) p=0.176  fails (out) p=0.183  fails (out) p=0.183
                       unit     Pre 5-10 (code)       (a') per series, median        holds p=0.012        holds p=0.012        holds p=0.013        holds p=0.013
                                                      (a) per series, mean     fails (out) p=0.322  fails (out) p=0.322   fails (out) p=0.21   fails (out) p=0.21
                                                      (b) concatenated               holds p=0.044        holds p=0.044        holds p=0.046        holds p=0.046
                                Pre all 10 (table S3) (a') per series, median        holds p=0.002        holds p=0.002        holds p=0.003        holds p=0.003
                                                      (a) per series, mean     fails (out) p=0.304  fails (out) p=0.304  fails (out) p=0.165  fails (out) p=0.165
                                                      (b) concatenated               holds p=0.031        holds p=0.031        holds p=0.033        holds p=0.033
without sub-5 code     pooled   Pre 5-10 (code)       (a') per series, median        holds p=0.008        holds p=0.008        holds p=0.011        holds p=0.011
                                                      (a) per series, mean     holds (out) p=0.011  holds (out) p=0.011  holds (out) p=0.008  holds (out) p=0.008
                                                      (b) concatenated               holds p=0.004        holds p=0.004        holds p=0.004        holds p=0.004
                                Pre all 10 (table S3) (a') per series, median        holds p=0.001        holds p=0.001        holds p=0.003        holds p=0.003
                                                      (a) per series, mean     holds (out) p=0.018  holds (out) p=0.037  holds (out) p=0.012  holds (out) p=0.005
                                                      (b) concatenated               holds p=0.001        holds p=0.001        holds p=0.001        holds p=0.001
                       unit     Pre 5-10 (code)       (a') per series, median        holds p=0.003        holds p=0.003        holds p=0.006        holds p=0.006
                                                      (a) per series, mean     holds (out) p=0.021  holds (out) p=0.021        holds p=0.006        holds p=0.006
                                                      (b) concatenated               holds p=0.004        holds p=0.004        holds p=0.004        holds p=0.004
                                Pre all 10 (table S3) (a') per series, median        holds p=0.001        holds p=0.001        holds p=0.002        holds p=0.002
                                                      (a) per series, mean     holds (out) p=0.031  holds (out) p=0.031        holds p=0.003        holds p=0.003
                                                      (b) concatenated               holds p=0.002        holds p=0.002        holds p=0.002        holds p=0.002
              infusion pooled   Pre 5-10 (code)       (a') per series, median        holds p=0.016        holds p=0.016        holds p=0.012        holds p=0.012
                                                      (a) per series, mean     holds (out) p=0.006  holds (out) p=0.006  holds (out) p=0.003  holds (out) p=0.003
                                                      (b) concatenated               holds p=0.009        holds p=0.009        holds p=0.008        holds p=0.008
                                Pre all 10 (table S3) (a') per series, median        holds p=0.002        holds p=0.002        holds p=0.002        holds p=0.002
                                                      (a) per series, mean     holds (out) p=0.008  holds (out) p=0.018  holds (out) p=0.005  holds (out) p=0.002
                                                      (b) concatenated               holds p=0.001        holds p=0.001        holds p=0.001        holds p=0.001
                       unit     Pre 5-10 (code)       (a') per series, median        holds p=0.011        holds p=0.011        holds p=0.012        holds p=0.012
                                                      (a) per series, mean     holds (out) p=0.009  holds (out) p=0.009        holds p=0.007        holds p=0.007
                                                      (b) concatenated               holds p=0.008        holds p=0.008        holds p=0.008        holds p=0.008
                                Pre all 10 (table S3) (a') per series, median        holds p=0.002        holds p=0.002        holds p=0.002        holds p=0.002
                                                      (a) per series, mean     holds (out) p=0.008  holds (out) p=0.008        holds p=0.001        holds p=0.001
                                                      (b) concatenated               holds p=0.002        holds p=0.002        holds p=0.002        holds p=0.002

rows = []
for rd, ks0 in READINGS.items():
    ks = [k for k in ks0 if k in TONE_SESSIONS]
    for split in ("infusion", "code"):
        for var in ("unit", "pooled"):
            for lab, sname, meas, cp, ck in (("(a) Methods", "methods", "per_series_mean_r", "pre", "ket"),
                                             ("(b) code", "code", "concat_median_abs_r", "pre", "ket"),
                                             ("(b') code, puff windows", "code", "concat_median_abs_r", "pre_puffwin", "ket_puffwin")):
                f = (lambda d: np.nanmean(d[meas])) if meas == "per_series_mean_r" else (lambda d: d[meas])
                pv = np.array([f(res[k]["tone"][split]["coupling"][(var, sname)][cp]) for k in ks])
                kv = np.array([f(res[k]["tone"][split]["coupling"][(var, sname)][ck]) for k in ks])
                ok = ~(np.isnan(pv) | np.isnan(kv))
                t = scipy.stats.ttest_rel(pv[ok], kv[ok])
                ratio = kv[ok] / pv[ok]
                rows.append({"reading": rd, "split": split, "variance": var, "measure": lab, "sessions": int(ok.sum()),
                             "Ket/Pre mean": ratio.mean(), "t": t.statistic, "p": t.pvalue,
                             "holds": t.pvalue < 0.05 and kv[ok].mean() < pv[ok].mean(), "within": abs(ratio.mean() - 0.46) <= 0.16})
print(pd.DataFrame(rows).round(4).to_string(index=False))

      reading    split variance                 measure  sessions  Ket/Pre mean       t      p  holds  within
without sub-1 infusion     unit             (a) Methods         8        1.0055 -0.0895 0.9312  False   False
without sub-1 infusion     unit                (b) code         7        0.8315  1.4079 0.2088  False   False
without sub-1 infusion     unit (b') code, puff windows         7        0.6486  2.4483 0.0499   True   False
without sub-1 infusion   pooled             (a) Methods         8        0.7711  1.1668 0.2815  False   False
without sub-1 infusion   pooled                (b) code         8        0.6156  2.4817 0.0421   True    True
without sub-1 infusion   pooled (b') code, puff windows         8        0.4650  3.8310 0.0064   True    True
without sub-1     code     unit             (a) Methods         8        1.1107 -0.0226 0.9826  False   False
without sub-1     code     unit                (b) code         8        0.9332  0.7171 0.4965  False   False
without sub-1     code     unit (b') code, puff windows         8        0.6180  2.6320 0.0338   True    True
without sub-1     code   pooled             (a) Methods         8        0.8467  0.9236 0.3864  False   False
without sub-1     code   pooled                (b) code         8        0.8728  1.4690 0.1853  False   False
without sub-1     code   pooled (b') code, puff windows         8        0.4748  3.7229 0.0074   True    True
without sub-3 infusion     unit             (a) Methods         8        1.0455 -0.2623 0.8006  False   False
without sub-3 infusion     unit                (b) code         7        0.8762  1.1419 0.2970  False   False
without sub-3 infusion     unit (b') code, puff windows         7        0.7068  2.1008 0.0804  False   False
without sub-3 infusion   pooled             (a) Methods         8        0.7935  0.9408 0.3781  False   False
without sub-3 infusion   pooled                (b) code         8        0.6519  2.1845 0.0652  False   False
without sub-3 infusion   pooled (b') code, puff windows         8        0.5134  3.2941 0.0132   True    True
without sub-3     code     unit             (a) Methods         8        1.1659 -0.2286 0.8257  False   False
without sub-3     code     unit                (b) code         8        1.0031  0.2479 0.8113  False   False
without sub-3     code     unit (b') code, puff windows         8        0.6639  2.3376 0.0520  False   False
without sub-3     code   pooled             (a) Methods         8        0.8673  0.6720 0.5232  False   False
without sub-3     code   pooled                (b) code         8        0.9283  1.1010 0.3073  False   False
without sub-3     code   pooled (b') code, puff windows         8        0.5162  3.2648 0.0138   True    True
without sub-5 infusion     unit             (a) Methods         7        1.2282  0.1313 0.8998  False   False
without sub-5 infusion     unit                (b) code         7        0.8205  1.5155 0.1804  False   False
without sub-5 infusion     unit (b') code, puff windows         7        0.6699  2.4305 0.0511  False   False
without sub-5 infusion   pooled             (a) Methods         7        0.8733  0.9031 0.4013  False   False
without sub-5 infusion   pooled                (b) code         7        0.6605  2.3073 0.0605  False   False
without sub-5 infusion   pooled (b') code, puff windows         7        0.5335  3.1818 0.0190   True    True
without sub-5     code     unit             (a) Methods         7        1.2434 -0.0500 0.9617  False   False
without sub-5     code     unit                (b) code         7        1.0654  0.4303 0.6820  False   False
without sub-5     code     unit (b') code, puff windows         7        0.6470  2.6412 0.0385   True   False
without sub-5     code   pooled             (a) Methods         7        0.9114  0.7999 0.4543  False   False
without sub-5     code   pooled                (b) code         7        0.9976  0.9146 0.3957  False   False
without sub-5     code   pooled (b') code, puff windows         7        0.5362  3.1693 0.0193   True    True

def corr_table(pre_key, ket_key):
    rows = []
    for k in names:
        reg = res[k]["regions_all"]
        ok = reg != "NA"
        d = (res[k][ket_key] - res[k][pre_key])[np.ix_(ok, ok)].astype(float)
        r = reg[ok]
        same = r[:, None] == r[None, :]
        iu = np.triu(np.ones_like(same), 1).astype(bool)
        fin = np.isfinite(d)
        row = {"session": k, "units": int(ok.sum()), "within change": np.nanmean(d[same & iu & fin]), "between change": np.nanmean(d[~same & iu & fin]),
               "pairs dropped (silent unit)": int((iu & ~fin).sum())}
        for fam in ("RSP", "ORB"):
            m = np.array([K.under(x, fam, paths) for x in r])
            row[f"{fam} units"] = int(m.sum())
            row[f"{fam} within change"] = np.nanmean(d[np.outer(m, m) & same & iu & fin]) if m.sum() > 1 else np.nan
            row[f"{fam} family change"] = np.nanmean(d[np.outer(m, m) & iu & fin]) if m.sum() > 1 else np.nan
        rows.append(row)
    return pd.DataFrame(rows).set_index("session")

def corr_summary(c30):
    tw = scipy.stats.ttest_1samp(c30["within change"], 0); tb = scipy.stats.ttest_1samp(c30["between change"], 0)
    tp = scipy.stats.ttest_rel(c30["within change"], c30["between change"])
    print(f"within-region change: mean {c30['within change'].mean():.4f}, t = {tw.statistic:.2f}, p = {tw.pvalue:.4f}")
    print(f"between-region change: mean {c30['between change'].mean():.4f}, t = {tb.statistic:.2f}, p = {tb.pvalue:.4f}")
    print(f"within against between (paired): t = {tp.statistic:.2f}, p = {tp.pvalue:.4f}")
    for fam in ("RSP", "ORB"):
        for col in ("within change", "family change"):
            v = c30[f"{fam} {col}"].dropna()
            print(f"{fam} {col}: positive in {(v > 0).sum()} of {len(v)} sessions recording it")

# criterion: Ketamine is the 10 min from the code's ketamine onset (30 s after the last puff of series 20)
c30 = corr_table("corr_pre", "corr_ket")
print("Ketamine window from the code's onset (criterion):")
print(c30.round(4).to_string())
corr_summary(c30)
# verification variant: the 10 min from the file's infusion time (Methods: "a 10-minute period following the
# onset of ketamine infusion"), Pre everything before the infusion; the code's onset is 94 to 420 s earlier
c30i = corr_table("corr_pre_inf", "corr_ket_inf")
print("\nvariant, Ketamine window from the file's infusion time:")
print(c30i[["within change", "between change", "RSP within change", "ORB within change"]].round(4).to_string())
corr_summary(c30i)

Ketamine window from the code's onset (criterion):
                    units  within change  between change  pairs dropped (silent unit)  RSP units  RSP within change  RSP family change  ORB units  ORB within change  ORB family change
session                                                                                                                                                                                
sub-0_ses-20230901    818         0.0019         -0.0042                         4075         23             0.0422             0.0418         83            -0.0041            -0.0031
sub-1_ses-20230218    495        -0.0014         -0.0038                         3924         16             0.0680             0.0658          0                NaN                NaN
sub-1_ses-20230226    984         0.0050          0.0028                         3926         56             0.0433             0.0469          0                NaN                NaN
sub-2_ses-20230902    677         0.0213          0.0013                         6715        115             0.0647             0.0415        128             0.0165             0.0196
sub-3_ses-20230218    388         0.0101          0.0008                          387          8             0.0284             0.0204        108             0.0156             0.0121
sub-3_ses-20230225    864        -0.0002         -0.0006                            0          3            -0.0446             0.0355         89            -0.0060            -0.0043
sub-4_ses-20230907    869         0.0312          0.0048                          868        109             0.1392             0.1176        177             0.0135             0.0092
sub-5_ses-20230706    139        -0.0089         -0.0078                         2570          0                NaN                NaN          0                NaN                NaN
sub-5_ses-20230822    499         0.0048         -0.0000                         1986          0                NaN                NaN          0                NaN                NaN
sub-6_ses-20230823    480        -0.0019         -0.0019                         1434          0                NaN                NaN          0                NaN                NaN
sub-7_ses-20220831    858         0.0019          0.0028                         3422          2             0.0018             0.0018        128             0.0098             0.0113
sub-8_ses-20220228    762         0.0061         -0.0014                         2280         69             0.0663             0.0103          0                NaN                NaN
sub-9_ses-20230908    854         0.0120          0.0054                         1705        105             0.0171             0.0159        123             0.0037             0.0038
within-region change: mean 0.0063, t = 2.16, p = 0.0522
between-region change: mean -0.0001, t = -0.14, p = 0.8934
within against between (paired): t = 2.82, p = 0.0155
RSP within change: positive in 9 of 10 sessions recording it
RSP family change: positive in 10 of 10 sessions recording it
ORB within change: positive in 5 of 7 sessions recording it
ORB family change: positive in 5 of 7 sessions recording it

variant, Ketamine window from the file's infusion time:
                    within change  between change  RSP within change  ORB within change
session                                                                                
sub-0_ses-20230901         0.0002         -0.0054             0.0510            -0.0071
sub-1_ses-20230218         0.0020         -0.0036             0.0855                NaN
sub-1_ses-20230226         0.0121          0.0055             0.1331                NaN
sub-2_ses-20230902         0.0282          0.0030             0.0820             0.0228
sub-3_ses-20230218         0.0091          0.0008             0.0429             0.0169
sub-3_ses-20230225         0.0002          0.0009            -0.1472            -0.0065
sub-4_ses-20230907         0.0353          0.0052             0.1603             0.0159
sub-5_ses-20230706        -0.0025         -0.0109                NaN                NaN
sub-5_ses-20230822         0.0043          0.0003                NaN                NaN
sub-6_ses-20230823        -0.0003         -0.0006                NaN                NaN
sub-7_ses-20220831         0.0112          0.0047            -0.0198             0.0004
sub-8_ses-20220228         0.0097          0.0008             0.0813                NaN
sub-9_ses-20230908         0.0165          0.0020             0.0321             0.0110
within-region change: mean 0.0097, t = 3.05, p = 0.0100
between-region change: mean 0.0002, t = 0.17, p = 0.8656
within against between (paired): t = 3.80, p = 0.0025
RSP within change: positive in 8 of 10 sessions recording it
RSP family change: positive in 8 of 10 sessions recording it
ORB within change: positive in 5 of 7 sessions recording it
ORB family change: positive in 4 of 7 sessions recording it

rows = []
for k in names:
    reg = np.array([K.reduce_region(x) for x in res[k]["regions_all"]])
    for after in ("after5", "after10"):
        ch = res[k]["spec_peak"][after] - res[k]["spec_peak"]["before"]
        for r in np.unique(reg):
            m = reg == r
            rows.append({"session": k, "after": after, "region": r, "rsp": K.under(r, "RSP", paths), "units": int(m.sum()), "change": np.nanmean(ch[m]),
                         "values": ch[m]})
sp = pd.DataFrame(rows)
sp = sp[sp.region != "NA"]
for after in ("after5", "after10"):
    s = sp[sp.after == after]
    rsp_sess = s[s.rsp].groupby("session").apply(lambda g: np.nanmean(np.concatenate(g["values"].to_list())))
    pooled = s.groupby("region").apply(lambda g: pd.Series({"units": g["units"].sum(), "change": np.nanmean(np.concatenate(g["values"].to_list()))}))
    pooled = pooled[pooled.units >= 30].sort_values("change", ascending=False)
    print(f"\n{after}: RSP units' mean change in z-scored low-frequency peak per session recording RSP:")
    print(rsp_sess.round(3).to_string())
    undefined = list(rsp_sess.index[rsp_sess.isna()])
    print(f"rises in {(rsp_sess > 0).sum()} of {rsp_sess.notna().sum()} sessions with a defined RSP value "
          f"({len(rsp_sess)} record RSP; no value in {undefined})")
    print("pooled regions with >= 30 units, largest rises first:")
    print(pooled.head(10).round(3).to_string())
    print(f"rank of the first RSP region: {[i for i, r in enumerate(pooled.index) if K.under(r, 'RSP', paths)][:1]} of {len(pooled)}")
print("frequency resolution (Hz):", res[names[0]]["spec_freqs"][1])

after5: RSP units' mean change in z-scored low-frequency peak per session recording RSP:
session
sub-0_ses-20230901    0.613
sub-1_ses-20230218   -3.178
sub-1_ses-20230226    4.090
sub-2_ses-20230902    1.645
sub-3_ses-20230218    0.276
sub-3_ses-20230225   -2.949
sub-4_ses-20230907    4.503
sub-7_ses-20220831      NaN
sub-8_ses-20220228    1.541
sub-9_ses-20230908    0.095
rises in 7 of 9 sessions with a defined RSP value (10 record RSP; no value in ['sub-7_ses-20220831'])
pooled regions with >= 30 units, largest rises first:
        units  change
region               
RSPv    430.0   2.128
TH       62.0   1.714
RSPagl   37.0   1.237
RSPd     39.0   0.826
NLOT2    32.0   0.688
SF       40.0   0.680
PO      345.0   0.569
SUB      67.0   0.505
DP       37.0   0.141
PoT      95.0   0.099
rank of the first RSP region: [0] of 52

after10: RSP units' mean change in z-scored low-frequency peak per session recording RSP:
session
sub-0_ses-20230901    2.557
sub-1_ses-20230218   -3.232
sub-1_ses-20230226    3.928
sub-2_ses-20230902    1.295
sub-3_ses-20230218    2.102
sub-3_ses-20230225   -7.901
sub-4_ses-20230907   -0.239
sub-7_ses-20220831      NaN
sub-8_ses-20220228   -0.321
sub-9_ses-20230908    0.471
rises in 5 of 9 sessions with a defined RSP value (10 record RSP; no value in ['sub-7_ses-20220831'])
pooled regions with >= 30 units, largest rises first:
        units  change
region               
RSPd     39.0   2.133
RSPv    430.0   0.740
TH       62.0   0.611
STR      59.0   0.516
SUB      67.0   0.453
AON      89.0   0.280
DP       37.0   0.157
BST      35.0   0.031
ORBvl   168.0   0.007
SF       40.0  -0.086
rank of the first RSP region: [0] of 52
frequency resolution (Hz): 0.048828125

Claim	In this notebook
`fig1e-human-ketamine-late-closure`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs1b-human-early-closure-unchanged`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig1f-human-reports`	not tested (see the attempt record)
`fig2d-cadss-ketamine`	not tested (see the attempt record)
`fig2g-six-factors`	tested in reproduce_kauvar_2025_ieeg.ipynb (DANDI:001347)
`fig2i-factors-1-2-ketamine`	tested in reproduce_kauvar_2025_ieeg.ipynb (DANDI:001347)
`fig3b-cluster-test-networks`	tested in reproduce_kauvar_2025_ieeg.ipynb (DANDI:001347)
`figs12-erp-not-blink`	not tested (see the attempt record)
`figs14g-cadss-correlation`	not tested (see the attempt record)
`fig6c-human-persistent`	not tested (see the attempt record)
`fig6e-human-fast-unchanged`	not tested (see the attempt record)
`fig6g-human-timescale`	not tested (see the attempt record)
`fig6i-human-phase-locking`	not tested (see the attempt record)
`fig1k-mouse-ketamine-late-closure`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig1l-dissociatives-preinfusion`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs2bc-mouse-early-late`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig1m-accumulation-persistence`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs2g-scaling-puff-number`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs2h-dose-response`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig1p-anesthesia-abolishes-both`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs5-other-drugs`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs6-vpm-muscimol`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs4-reward-licking`	not tested (see the attempt record)
`fig5hi-behavior-model-doses`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig4-recordings-scale`	tested here
`fig4b-recorded-eye-closure-ketamine`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig5b-recorded-accumulation`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`fig4de-response-clusters`	tested here
`figs31b-cluster-overrepresentation`	tested here
`fig4g-ketamine-alters-slow-clusters`	tested here
`fig4i-rise-times-unchanged`	tested here
`fig4h-decay-shortened`	tested here
`fig4h-decay-anatomical-progression`	tested here
`fig5d-emotion-dimension-accumulates`	tested here
`fig5e-late-reduced`	tested here
`fig5e-rapid-unchanged`	tested here
`fig5g-tau-persist-model`	tested here
`figs33-alternative-models`	tested here
`figs32c-regional-emotion-dimension`	tested here
`fig6b-persistent-dimension-distributed`	tested here
`fig6d-mouse-persistent`	tested here
`fig6f-mouse-fast-unchanged`	tested here
`fig6h-mouse-timescale`	tested here
`fig6j-mouse-coupling`	tested here
`figs34c-tone-no-eye-closure`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`figs34d-tone-fast-unchanged`	tested here
`figs34e-tone-persistent`	tested here
`figs34f-tone-shorter-than-puff`	tested here
`figs34g-tone-timescale`	tested here
`figs34h-tone-coupling`	tested here
`figs30-correlation-structure`	tested here
`figs29-rsp-oscillations`	tested here
`link-001347-condition-labels`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)
`link-001347-recorded-sessions`	tested in reproduce_kauvar_2025_behavior.ipynb (DANDI:001347)

Claim	Outcome
fig4-recordings-scale	reproduced
fig4de-response-clusters	not reproduced
figs31b-cluster-overrepresentation	partially reproduced
fig4g-ketamine-alters-slow-clusters	not reproduced
fig4i-rise-times-unchanged	reproduced
fig4h-decay-shortened	reproduced
fig4h-decay-anatomical-progression	reproduced
fig5d-emotion-dimension-accumulates	partially reproduced
fig5e-late-reduced	partially reproduced
fig5e-rapid-unchanged	reproduced
figs32c-regional-emotion-dimension	partially reproduced
fig5g-tau-persist-model	reproduced
figs33-alternative-models	reproduced (assumes the Methods target window and per-unit variance where the readings disagree)
fig6d-mouse-persistent	reproduced (assumes sub-5 is the excluded animal)
fig6f-mouse-fast-unchanged	reproduced
fig6b-persistent-dimension-distributed	reproduced
fig6h-mouse-timescale	reproduced (assumes sub-5 is the excluded animal)
fig6j-mouse-coupling	partially reproduced (assumes sub-5 is the excluded animal)
figs34d-tone-fast-unchanged	reproduced (assumes sub-5 is the excluded animal)
figs34e-tone-persistent	reproduced
figs34f-tone-shorter-than-puff	reproduced
figs34g-tone-timescale	partially reproduced (assumes sub-5 is the excluded animal)
figs34h-tone-coupling	not reproduced
figs30-correlation-structure	partially reproduced
figs29-rsp-oscillations	reproduced

Kauvar, Richman, Liu et al. (2025) on dandiset 001326¶

Loading¶

Data check¶

fig4-recordings-scale¶

fig4de-response-clusters¶

figs31b-cluster-overrepresentation¶

fig4g-ketamine-alters-slow-clusters¶

Stability of the clustering (fig4de, figs31b, fig4g)¶

fig4i-rise-times-unchanged¶

fig4h-decay-shortened¶

fig4h-decay-anatomical-progression¶

Session readings for Figures 5 and 6¶

fig6d-mouse-persistent¶

fig6f-mouse-fast-unchanged¶

fig6b-persistent-dimension-distributed¶

fig6h-mouse-timescale¶

fig6j-mouse-coupling¶

fig5d-emotion-dimension-accumulates¶

fig5e-late-reduced¶

fig5e-rapid-unchanged¶

figs32c-regional-emotion-dimension¶

fig5g-tau-persist-model¶

figs33-alternative-models¶

figs34d-tone-fast-unchanged¶

figs34e-tone-persistent¶

figs34f-tone-shorter-than-puff¶

figs34g-tone-timescale¶

figs34h-tone-coupling¶

figs30-correlation-structure¶

figs29-rsp-oscillations¶

Summary¶