import time
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import curve_fit
from statsmodels.stats.multitest import multipletests

import shield

pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore")
started = time.time()

# Stream the session files: tables first, then one pass over each session's spikes.
assets = shield.list_assets()
session_assets = assets[~assets["path"].str.contains("_probe-")].reset_index(drop=True)
sessions, units = shield.scan_all(session_assets["asset_id"].tolist())
sessions["path"] = session_assets["path"].values
print(f"tables of {len(sessions)} sessions read in {time.time() - started:.0f} s")
SP = shield.spike_pass_all(sessions, units)
print(f"spike pass done after {time.time() - started:.0f} s")

tables of 99 sessions read in 63 s

spike pass done after 715 s

# The paper's code: the 87 sessions passing ephys QC, its unit and session tables (pinned commit)
PASS = set(shield.passing_session_ids())
sessions["passing"] = sessions["session_id"].isin(PASS)
sessions["day"] = sessions["session_type"].str.replace("EPHYS_", "").astype(int)
units = units.join(sessions.set_index("session_id")[["subject_id", "session_type", "passing", "genotype"]], on="session_id")
units["qc"] = shield.qc_pass(units)
units["good"] = units["quality"] == "good"
units["in_area"] = ~units["structure_acronym"].isin(shield.NOT_IN_AREA)
code_units = shield.read_code_csv("metadata_tables/units_table.csv")
vbn_units = shield.read_code_csv("metadata_tables/VBN_units_table.csv")
vbn_sessions = shield.read_code_csv("metadata_tables/VBN_sessions.csv")
code_sessions = shield.read_code_csv("metadata_tables/dynamic_gating_session_metadata.csv")
print("code units table:", code_units.shape, "; VBN units:", vbn_units.shape, "; VBN sessions:", vbn_sessions.shape)

code units table: (201520, 47) ; VBN units: (319013, 27) ; VBN sessions: (153, 22)

meta = shield.version_metadata()
print("Dandiset", shield.DANDISET_ID, "version", shield.DANDISET_VERSION, ":", meta["name"])
print("assets:", len(assets), "; session files:", len(session_assets), "; LFP probe files:", len(assets) - len(session_assets),
      f"; total {assets['size'].sum() / 1e12:.2f} TB")
print("subjects in the archive summary:", meta["assetsSummary"].get("numberOfSubjects"))
print("\nsessions per subject and genotype:")
print(sessions.groupby(["genotype"]).agg(mice=("subject_id", "nunique"), sessions=("session_id", "size")).to_string())
print("\nsex by mouse:", sessions.groupby("subject_id")["sex"].first().value_counts().to_dict())
print("session types:", sessions["session_type"].value_counts().sort_index().to_dict())
print("duplicate session identifiers:", int(sessions["session_id"].duplicated().sum()))
print("sessions in the code's 87-session list found in the dandiset:", int(sessions["passing"].sum()), "of", len(PASS))
print("\nsessions outside the 87:")
print(sessions.loc[~sessions["passing"], ["path", "session_id", "session_type", "subject_id", "n_units"]].to_string(index=False))

Dandiset 001051 version 0.260922.1827 : Large-scale Neuropixels recordings through SHIELD implant during visual change detection task with dynamic gating of engagement
assets: 666 ; session files: 99 ; LFP probe files: 567 ; total 1.81 TB
subjects in the archive summary: 27

sessions per subject and genotype:
                                               mice  sessions
genotype                                                     
Ai32(RCL-ChR2(H134R)_EYFP)/wt                     1         3
Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt    16        58
Vip-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt     4        15
wt/wt                                             6        23

sex by mouse: {'M': 18, 'F': 9}
session types: {'EPHYS_1': 24, 'EPHYS_2': 26, 'EPHYS_3': 25, 'EPHYS_4': 24}
duplicate session identifiers: 0
sessions in the code's 87-session list found in the dandiset: 87 of 87

sessions outside the 87:
                                         path  session_id session_type subject_id  n_units
sub-604914/sub-604914_ses-20220427T041046.nwb  1173189336      EPHYS_2     604914     1758
sub-604914/sub-604914_ses-20220429T060155.nwb  1173741216      EPHYS_4     604914     2262
sub-614608/sub-614608_ses-20220629T220553.nwb  1187940755      EPHYS_3     614608     2356
sub-623322/sub-623322_ses-20220825T212029.nwb  1200878248      EPHYS_4     623322     1854
sub-623786/sub-623786_ses-20220915T202517.nwb  1211232400      EPHYS_4     623786     2351
sub-633232/sub-633232_ses-20221108T000611.nwb  1224385827      EPHYS_1     633232     1088
sub-633232/sub-633232_ses-20221108T230110.nwb  1224635971      EPHYS_2     633232     1933
sub-633232/sub-633232_ses-20221109T225948.nwb  1224928532      EPHYS_3     633232     1451
sub-633232/sub-633232_ses-20221110T225606.nwb  1225178442      EPHYS_4     633232     1461
sub-636740/sub-636740_ses-20221117T230548.nwb  1226780656      EPHYS_4     636740     2402
sub-637483/sub-637483_ses-20221123T172329.nwb  1228018071      EPHYS_3     637483     1817
sub-637488/sub-637488_ses-20221209T060640.nwb  1231476077      EPHYS_4     637488     2500

print("units:", len(units), "; quality values:", units["quality"].value_counts().to_dict())
print("isi_label values:", units["isi_label"].value_counts().to_dict())
cols = ["presence_ratio", "amplitude_cutoff", "isi_violations", "snr", "isolation_distance", "max_drift",
        "cumulative_drift", "amplitude", "firing_rate", "waveform_duration", "x", "y", "z"]
print(units[cols].describe(percentiles=[0.01, 0.5, 0.99]).T.round(3).to_string())
print("\nCCF coordinates (x anterior-posterior, y dorsal-ventral, z left-right, micrometres); midline z = 5700.")
inb = units[units["in_area"]]
print("annotated units with z > 5700:", int((inb["z"] > 5700).sum()), "; largest z:", inb["z"].max(),
      "; structures:", inb.loc[inb["z"] > 5700, "structure_acronym"].value_counts().head(8).to_dict())
print("units with coordinate -25 (no position):", int((units["x"] == -25).sum()),
      "; of those annotated in an area:", int(((units["x"] == -25) & units["in_area"]).sum()))

units: 201192 ; quality values: {'good': 163145, 'noise': 38047}
isi_label values: {'No label': 180372, 'VISp': 7834, 'VIS': 4715, 'VISlm': 2806, 'VISrl': 2139, 'VISam': 854, 'VISpm': 803, 'VISal': 731, 'VISl': 337, 'VISrll': 327, 'VISmma': 274}

                       count          mean           std     min      1%       50%           99%           max
presence_ratio      201192.0  8.800000e-01  2.220000e-01   0.010   0.080     0.990  9.900000e-01  9.900000e-01
amplitude_cutoff    201192.0  1.790000e-01  2.000000e-01   0.000   0.000     0.069  5.000000e-01  5.000000e-01
isi_violations      201192.0  1.228800e+01  6.511460e+02   0.000   0.000     0.072  6.280700e+01  1.473621e+05
snr                 201192.0  2.684000e+00  1.328000e+00   0.000   0.571     2.426  7.137000e+00  3.534900e+01
isolation_distance  201192.0  8.862258e+20  3.532465e+23   0.000   0.000    55.163  2.745429e+13  1.570126e+26
max_drift           201192.0  2.926000e+01  2.454000e+01   0.000   0.000    23.555  1.162310e+02  2.743800e+02
cumulative_drift    201192.0  1.408630e+02  1.421400e+02   0.000   0.000   108.800  6.986770e+02  2.316120e+03
amplitude           201192.0  1.518130e+02  1.091950e+02   0.847  29.944   122.620  5.773860e+02  2.551869e+03
firing_rate         201192.0  5.320000e+00  7.510000e+00   0.000   0.003     2.795  3.536300e+01  1.190590e+02
waveform_duration   201192.0  5.720000e-01  2.960000e-01   0.014   0.151     0.563  2.019000e+00  2.692000e+00
x                   201192.0  5.713993e+03  2.786315e+03 -25.000  -1.000  6400.000  9.300000e+03  9.800000e+03
y                   201192.0  2.124297e+03  1.272906e+03 -25.000  -1.000  2125.000  4.575000e+03  5.400000e+03
z                   201192.0  3.438140e+03  1.612368e+03 -25.000  -1.000  3825.000  5.525000e+03  6.250000e+03

CCF coordinates (x anterior-posterior, y dorsal-ventral, z left-right, micrometres); midline z = 5700.
annotated units with z > 5700: 464 ; largest z: 6250.0 ; structures: {'LSr': 142, 'ACB': 70, 'STR': 65, 'TTd': 64, 'SFO': 51, 'ILA': 34, 'SPA': 18, 'IMD': 12}
units with coordinate -25 (no position): 145 ; of those annotated in an area: 0

spt = pd.DataFrame([{"session_id": k, "first_stim": v["first_stim"], "behavior_start": v["beh_start"],
                     "behavior_end": v["beh_end"], "last_stim_end": v["last_stim_end"],
                     "images": ",".join(v["image_names"]), "im115_presentations": v["im115_n"]} for k, v in SP.items()])
sessions = sessions.merge(spt, on="session_id")
sessions["behavior_min"] = (sessions["behavior_end"] - sessions["behavior_start"]) / 60
print(sessions[["behavior_min", "first_stim", "last_stim_end", "im115_presentations"]].describe().round(1).to_string())
print("\nimage sets:", sessions["images"].value_counts().to_dict())
print("sessions without eye tracking:", sessions.loc[~sessions["has_eye"], "session_id"].tolist())
print("rewards per session (all / auto / earned):")
print(sessions[["n_rewards", "n_rewards_auto", "n_rewards_earned", "n_hit", "n_trials"]].describe().round(1).to_string())
print("earned rewards equal hit trials in", int((sessions["n_rewards_earned"] == sessions["n_hit"]).sum()), "of", len(sessions), "sessions")

       behavior_min  first_stim  last_stim_end  im115_presentations
count          99.0        99.0           99.0                 99.0
mean           65.1       101.0         5547.5               1219.6
std             0.0       101.6          103.2                 97.7
min            65.0        25.5         5459.5               1030.0
25%            65.0        35.1         5479.8               1153.5
50%            65.0        47.4         5499.2               1219.0
75%            65.1       137.0         5587.3               1279.0
max            65.2       635.7         6077.1               1506.0

image sets: {'im047_r-1.0,im114_r-1.0,im115_r-0.7,im115_r-1.0,omitted': 26, 'im005_r-1.0,im047_r-1.0,im115_r-0.7,im115_r-1.0,omitted': 25, 'im047_r-1.0,im087_r-1.0,im115_r-0.7,im115_r-1.0,omitted': 24, 'im047_r-1.0,im104_r-1.0,im115_r-0.7,im115_r-1.0,omitted': 24}
sessions without eye tracking: [1187475832, 1202438738]
rewards per session (all / auto / earned):
       n_rewards  n_rewards_auto  n_rewards_earned  n_hit  n_trials
count       99.0            99.0              99.0   99.0      99.0
mean        62.7             5.6              57.1   83.4     668.0
std         32.4             0.7              32.1   35.0     144.6
min          6.0             3.0               0.0    0.0     439.0
25%         38.5             6.0              32.5   63.5     566.0
50%         63.0             6.0              57.0   88.0     638.0
75%         82.5             6.0              76.5  106.0     730.0
max        138.0             6.0             132.0  162.0    1182.0
earned rewards equal hit trials in 3 of 99 sessions

tt_run = np.arange(500) * 0.01 - 2; tt_pup = np.arange(200) * 0.01 - 0.5
rows = []
for k, v in SP.items():
    c = v["clock"]
    h, m = c[("running", "hit")], c[("running", "miss")]
    d = (h - np.nanmean(h[(tt_run >= -1) & (tt_run < 0)])) - (m - np.nanmean(m[(tt_run >= -1) & (tt_run < 0)]))
    w, b = c[("pupil", "white")], c[("pupil", "black")]
    dp = (w - np.nanmean(w[tt_pup < 0])) - (b - np.nanmean(b[tt_pup < 0]))
    row = {"session_id": k, "hits": c[("n", "hit")], "misses": c[("n", "miss")], "response_latency_s": c["response_latency"]}
    if np.isfinite(d).any():
        post = (tt_run >= 0) & (tt_run < 1.5); pre = (tt_run >= -1) & (tt_run < 0)
        thr = 3 * np.nanmax(np.abs(d[pre]))
        cross = np.where(post & (d < -thr))[0]
        row.update(run_pre_max=np.nanmax(np.abs(d[pre])), run_post_min=np.nanmin(d[post]),
                   run_min_time_s=tt_run[post][np.nanargmin(d[post])], run_onset_s=tt_run[cross[0]] if len(cross) else np.nan)
    if np.isfinite(dp).any():
        post = tt_pup >= 0; pre = tt_pup < 0
        thr = 3 * np.nanmax(np.abs(dp[pre]))
        cross = np.where(post & (dp < -thr))[0]
        row.update(pupil_pre_max=np.nanmax(np.abs(dp[pre])), pupil_post_min=np.nanmin(dp[post]),
                   pupil_min_time_s=tt_pup[post][np.nanargmin(dp[post])], pupil_onset_s=tt_pup[cross[0]] if len(cross) else np.nan)
    rows.append(row)
ck = pd.DataFrame(rows).set_index("session_id")
print("running (cm/s; hit minus miss, each minus its -1 to 0 s mean) and pupil (area, pixels^2; white minus black flashes):")
print(ck.describe().round(3).to_string())
print("sessions without a running contrast:", ck.index[ck["run_post_min"].isna()].tolist(),
      "; without a pupil contrast:", ck.index[ck["pupil_post_min"].isna()].tolist())
print("sessions where the contrast after 0 is not more than 3 times the largest before 0 -- running:",
      ck.index[ck["run_post_min"].notna() & ck["run_onset_s"].isna()].tolist(),
      "; pupil:", ck.index[ck["pupil_post_min"].notna() & ck["pupil_onset_s"].isna()].tolist())
not_shown = sorted(set(ck.index[ck["run_post_min"].isna() | ck["run_onset_s"].isna()])
                   | set(ck.index[ck["pupil_post_min"].isna() | ck["pupil_onset_s"].isna()]))
print("sessions where alignment is not shown for at least one stream (no contrast, no eye tracking or a weak contrast):",
      len(not_shown), "of", len(ck), not_shown)
print("sessions where the minimum falls in the first 0.1 s -- running:", ck.index[ck["run_min_time_s"] < 0.1].tolist(),
      "; pupil:", ck.index[ck["pupil_min_time_s"] < 0.1].tolist())
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
R_ = np.array([v["clock"][("running", "hit")] - v["clock"][("running", "miss")] for v in SP.values()])
P_ = np.array([v["clock"][("pupil", "white")] - v["clock"][("pupil", "black")] for v in SP.values()])
ax[0].plot(tt_run, np.nanmean(R_, 0)); ax[0].axvline(0, color="k", lw=0.5); ax[0].set_xlabel("s from change"); ax[0].set_ylabel("hit minus miss running (cm/s)")
ax[1].plot(tt_pup, np.nanmean(P_, 0)); ax[1].axvline(0, color="k", lw=0.5); ax[1].set_xlabel("s from flash"); ax[1].set_ylabel("white minus black pupil area")
plt.show()
print(f"plotted ranges: running {np.nanmin(np.nanmean(R_, 0)):.2f} to {np.nanmax(np.nanmean(R_, 0)):.2f}, pupil {np.nanmin(np.nanmean(P_, 0)):.1f} to {np.nanmax(np.nanmean(P_, 0)):.1f}")

running (cm/s; hit minus miss, each minus its -1 to 0 s mean) and pupil (area, pixels^2; white minus black flashes):
          hits   misses  response_latency_s  pupil_pre_max  pupil_post_min  pupil_min_time_s  pupil_onset_s  run_pre_max  run_post_min  run_min_time_s  run_onset_s
count   99.000   99.000              98.000         97.000          97.000            97.000         94.000       97.000        97.000          97.000       92.000
mean    83.444  122.667               0.456         66.518        -827.689             0.602          0.366        2.751       -30.818           0.686        0.246
std     35.046   54.738               0.080        113.154         318.538             0.146          0.121        1.848        16.388           0.307        0.086
min      0.000    0.000               0.178         14.200       -1974.375             0.000          0.050        0.196       -59.045           0.370        0.130
25%     63.500   84.000               0.380         30.640        -963.825             0.570          0.320        1.550       -43.276           0.510        0.200
50%     88.000  113.000               0.479         38.703        -759.445             0.590          0.350        2.488       -35.702           0.570        0.225
75%    106.000  152.000               0.519         60.252        -624.976             0.610          0.388        3.644       -20.067           0.740        0.260
max    162.000  268.000               0.580       1062.558           7.917             1.490          1.330        9.176        -0.527           1.490        0.640
sessions without a running contrast: [1173189336, 1230732154] ; without a pupil contrast: [1187475832, 1202438738]
sessions where the contrast after 0 is not more than 3 times the largest before 0 -- running: [1182427414, 1182628226, 1178693650, 1191375883, 1226010097] ; pupil: [1173741216, 1231476077, 1244070169]
sessions where alignment is not shown for at least one stream (no contrast, no eye tracking or a weak contrast): 12 of 99 [1173189336, 1173741216, 1178693650, 1182427414, 1182628226, 1187475832, 1191375883, 1202438738, 1226010097, 1230732154, 1231476077, 1244070169]
sessions where the minimum falls in the first 0.1 s -- running: [] ; pupil: [1231476077]

plotted ranges: running -23.29 to 7.38, pupil -755.2 to 7.0

p_units = units[units["passing"]]
cu = code_units[code_units["session"].isin(PASS)]
print("code sessions in the list found as dandiset identifiers:", len(PASS & set(sessions["session_id"])), "of", len(PASS))
cs = code_sessions.set_index("ecephys_session_id")
mm = sessions.set_index("session_id").join(cs[["session_type", "mouse_id", "genotype"]], rsuffix="_code")
print("session type agrees:", int((mm["session_type"] == mm["session_type_code"]).sum()), "of", len(mm),
      "; mouse agrees:", int((mm["subject_id"].astype(int) == mm["mouse_id"]).sum()),
      "; genotype agrees:", int((mm["genotype"] == mm["genotype_code"]).sum()))
n_nwb = p_units.groupby("session_id").size()
n_code = cu.groupby("session").size()
cmp = pd.concat([n_nwb, n_code], axis=1, keys=["dandiset", "code"]).fillna(0).astype(int)
print("\nunits in the 87 sessions: dandiset", int(n_nwb.sum()), "; code", int(n_code.sum()))
print("sessions with the same unit count:", int((cmp["dandiset"] == cmp["code"]).sum()), "of", len(cmp))
print(cmp[cmp["dandiset"] != cmp["code"]].to_string())
ids_missing = cu[~cu["unit_id"].isin(units.index)]
print("code units not in the dandiset:", len(ids_missing), "; by structure:", ids_missing["structure_acronym"].value_counts().head(6).to_dict())
# which probes the code units without a dandiset unit sit on, and whether the session file has units on that probe
if "probe_id" in units.columns:
    d_probes = units.groupby("session_id")["probe_id"].apply(lambda x: set(x.dropna().astype(int)))
    pm = ids_missing.groupby(["session", "probe_id"]).size().rename("code_units").reset_index()
    pm["dandiset_units_on_probe"] = [int(((units["session_id"] == s_) & (units["probe_id"] == p_)).sum()) for s_, p_ in zip(pm["session"], pm["probe_id"])]
    pm["code_units_on_probe"] = [int(((cu["session"] == s_) & (cu["probe_id"] == p_)).sum()) for s_, p_ in zip(pm["session"], pm["probe_id"])]
    print("code units not in the dandiset, by session and the code's probe_id:")
    print(pm.to_string(index=False))
    print("dandiset probe_ids in those sessions:", {int(s_): sorted(d_probes.get(s_, set())) for s_ in pm["session"].unique()})
else:
    print("the session files' electrodes table has no probe_id column")
print("their peak channels among the dandiset's unit peak channels:", int(ids_missing["peak_channel_id"].isin(units["peak_channel_id"]).sum()), "of", len(ids_missing))
print("dandiset units (87 sessions) not in the code table:", int((~p_units.index.isin(cu["unit_id"])).sum()))
j = p_units.join(cu.set_index("unit_id")[["max_drift", "firing_rate", "quality", "presence_ratio", "snr"]], rsuffix="_code", how="inner")
rows = []
for sid, g in j.groupby("session_id"):
    rows.append({"session_id": sid, "units": len(g),
                 "max_drift_differs": float((np.abs(g["max_drift"] - g["max_drift_code"]) > 0.01).mean()),
                 "firing_rate_differs": float((np.abs(g["firing_rate"] - g["firing_rate_code"]) > 0.01).mean()),
                 "quality_differs": float((g["quality"] != g["quality_code"]).mean())})
agree = pd.DataFrame(rows)
bad = agree[(agree[["max_drift_differs", "firing_rate_differs", "quality_differs"]] > 0.01).any(axis=1)]
print("\nsessions where more than 1% of matched units differ in a metric:", len(bad), "of", len(agree))
print(bad.round(3).to_string(index=False))
print("in the other", len(agree) - len(bad), "sessions the largest fraction differing is",
      round(float(agree.drop(bad.index)[["max_drift_differs", "firing_rate_differs", "quality_differs"]].max().max()), 4))
# implied recording duration n_spikes / firing_rate per unit: constant within a session when the rate
# describes the same spikes; with the code's firing_rate it scatters if the code's units are others
jj = j.assign(dur_dandiset=j["n_spikes"] / j["firing_rate"], dur_code=j["n_spikes"] / j["firing_rate_code"])
dd = jj.groupby("session_id").agg(units=("dur_dandiset", "size"), dur_dandiset_median=("dur_dandiset", "median"),
                                  dandiset_cv=("dur_dandiset", lambda x: x.std() / x.mean()),
                                  code_cv=("dur_code", lambda x: x.std() / x.mean()),
                                  code_ratio_iqr=("dur_code", lambda x: float(np.subtract(*np.percentile(x / x.median(), [75, 25])))))
dd["differs"] = dd.index.isin(bad["session_id"])
print("\nspike counts against firing rates (implied duration n_spikes / firing_rate, s):")
print(dd.groupby("differs")[["dandiset_cv", "code_cv", "code_ratio_iqr"]].median().round(4).to_string())
print(dd[dd["differs"]].round(4).to_string())
print("largest across sessions:", dd.groupby("differs")[["dandiset_cv", "code_cv"]].max().round(4).to_dict("index"))

code sessions in the list found as dandiset identifiers: 87 of 87

session type agrees: 99 of 99 ; mouse agrees: 99 ; genotype agrees: 99

units in the 87 sessions: dandiset 177959 ; code 178287
sessions with the same unit count: 82 of 87
            dandiset  code
1176989662      1927  1988
1178693650      1278  1392
1179911454      1296  1353
1180266229      1867  1925
1181096406      2257  2295
code units not in the dandiset: 328 ; by structure: {'CP': 133, 'Track not annotated': 119, 'ILA': 16, 'VL': 12, 'MOp': 10, 'PL': 10}

code units not in the dandiset, by session and the code's probe_id:
   session  probe_id  code_units  dandiset_units_on_probe  code_units_on_probe
1176989662        88          61                        0                  491
1178693650        46         114                        0                  131
1179911454       121          57                        0                  413
1180266229       133          58                        0                  490
1181096406       139          38                        0                  649
dandiset probe_ids in those sessions: {1176989662: [83, 84, 85, 86, 87], 1178693650: [41, 42, 43, 44], 1179911454: [116, 117, 118, 119, 120], 1180266229: [128, 129, 130, 131, 132], 1181096406: [134, 135, 136, 137, 138]}
their peak channels among the dandiset's unit peak channels: 0 of 328
dandiset units (87 sessions) not in the code table: 0

sessions where more than 1% of matched units differ in a metric: 6 of 87
 session_id  units  max_drift_differs  firing_rate_differs  quality_differs
 1176989662   1927              0.993                0.994            0.280
 1178460518   1969              0.993                0.993            0.349
 1178693650   1278              0.996                0.995            0.350
 1179911454   1296              0.990                0.985            0.382
 1180266229   1867              0.998                0.994            0.284
 1181096406   2257              0.553                0.553            0.144
in the other 81 sessions the largest fraction differing is 0.0

spike counts against firing rates (implied duration n_spikes / firing_rate, s):
         dandiset_cv  code_cv  code_ratio_iqr
differs                                      
False            0.0   0.0000          0.0000
True             0.0   9.1413          4.9251
            units  dur_dandiset_median  dandiset_cv  code_cv  code_ratio_iqr  differs
session_id                                                                           
1176989662   1927            5247.0036       0.0759   9.3586          4.8571     True
1178460518   1969            6165.1432       0.0000   7.4913          4.9931     True
1178693650   1278            6297.1002       0.0000   5.3501          4.6123     True
1179911454   1296            6196.4103       0.0000  11.8224          6.5466     True
1180266229   1867            6318.4278       0.0000   8.9241          5.7766     True
1181096406   2257            6128.5493       0.0000  10.2842          0.6353     True
largest across sessions: {False: {'dandiset_cv': 0.0003, 'code_cv': 0.0003}, True: {'dandiset_cv': 0.0759, 'code_cv': 11.8224}}

# receptive-field p values computed here (chi-square bootstrap) against the code's p_value_rf
rf = pd.concat([v["rf"].assign(session_id=k) for k, v in SP.items()])
rfj = rf.join(code_units.set_index("unit_id")[["p_value_rf"]], how="inner").dropna()
ok_s = ~rfj["session_id"].isin(bad["session_id"])
for lab, g in (("sessions whose metrics agree", rfj[ok_s]), ("the sessions that differ", rfj[~ok_s])):
    print(f"{lab}: {len(g)} units; Pearson r of p values {np.corrcoef(g['p_rf'], g['p_value_rf'])[0, 1]:.4f}; "
          f"same significance at 0.01 for {np.mean((g['p_rf'] < 0.01) == (g['p_value_rf'] < 0.01)):.3f}; "
          f"significant here {np.mean(g['p_rf'] < 0.01):.3f}, in the code {np.mean(g['p_value_rf'] < 0.01):.3f}")
# per session, in the sessions whose metrics agree: fraction of units whose significance at 0.01 differs
ps_ = rfj[ok_s].assign(diff=(rfj["p_rf"] < 0.01) != (rfj["p_value_rf"] < 0.01)).groupby("session_id")["diff"].agg(["size", "mean"])
print(f"per session (sessions whose metrics agree, {len(ps_)} sessions with RF units): fraction differing median {ps_['mean'].median():.3f}, "
      f"max {ps_['mean'].max():.3f}; sessions with more than 1% differing: {int((ps_['mean'] > 0.01).sum())}")
print("those sessions (units, fraction differing):", ps_[ps_["mean"] > 0.01].round(3).to_dict("index"))
near = rfj[ok_s]
dif = near[(near["p_rf"] < 0.01) != (near["p_value_rf"] < 0.01)]
print(f"units whose significance differs: {len(dif)}; code p of these from {dif['p_value_rf'].min():.3f} to {dif['p_value_rf'].max():.3f}, "
      f"p here from {dif['p_rf'].min():.3f} to {dif['p_rf'].max():.3f}; with both p values within 0.005 of 0.01: "
      f"{int(((dif['p_rf'] - 0.01).abs() <= 0.005).sum() if len(dif) else 0)} of those with p here near, "
      f"{int((((dif['p_rf'] - 0.01).abs() <= 0.005) & ((dif['p_value_rf'] - 0.01).abs() <= 0.005)).sum())} both")

sessions whose metrics agree: 7095 units; Pearson r of p values 0.9990; same significance at 0.01 for 0.988; significant here 0.551, in the code 0.551
the sessions that differ: 251 units; Pearson r of p values 0.0343; same significance at 0.01 for 0.450; significant here 0.701, in the code 0.271
per session (sessions whose metrics agree, 79 sessions with RF units): fraction differing median 0.008, max 0.038; sessions with more than 1% differing: 33
those sessions (units, fraction differing): {1174553025: {'size': 49, 'mean': 0.02}, 1175253205: {'size': 57, 'mean': 0.035}, 1181731440: {'size': 60, 'mean': 0.033}, 1182871514: {'size': 84, 'mean': 0.012}, 1187668018: {'size': 137, 'mean': 0.015}, 1187940755: {'size': 126, 'mean': 0.016}, 1190892228: {'size': 28, 'mean': 0.036}, 1194850630: {'size': 84, 'mean': 0.012}, 1197967151: {'size': 173, 'mean': 0.017}, 1198130743: {'size': 95, 'mean': 0.021}, 1198995158: {'size': 37, 'mean': 0.027}, 1199249150: {'size': 81, 'mean': 0.012}, 1200471383: {'size': 72, 'mean': 0.014}, 1202438738: {'size': 79, 'mean': 0.038}, 1202644967: {'size': 144, 'mean': 0.028}, 1203574877: {'size': 42, 'mean': 0.024}, 1208053773: {'size': 29, 'mean': 0.034}, 1210343162: {'size': 132, 'mean': 0.03}, 1210863396: {'size': 117, 'mean': 0.026}, 1218577463: {'size': 136, 'mean': 0.015}, 1218795542: {'size': 176, 'mean': 0.011}, 1219127153: {'size': 39, 'mean': 0.026}, 1219401111: {'size': 135, 'mean': 0.015}, 1224635971: {'size': 99, 'mean': 0.03}, 1224928532: {'size': 131, 'mean': 0.031}, 1226251663: {'size': 198, 'mean': 0.025}, 1227855488: {'size': 225, 'mean': 0.013}, 1229429400: {'size': 111, 'mean': 0.018}, 1232430325: {'size': 190, 'mean': 0.011}, 1232954220: {'size': 176, 'mean': 0.011}, 1234203203: {'size': 141, 'mean': 0.014}, 1243523806: {'size': 104, 'mean': 0.029}, 1243799530: {'size': 104, 'mean': 0.019}}
units whose significance differs: 83; code p of these from 0.003 to 0.019, p here from 0.003 to 0.020; with both p values within 0.005 of 0.01: 78 of those with p here near, 66 both

ins = units.groupby(["session_id", "probe"]).size().rename("n").reset_index()
ins = ins.merge(sessions[["session_id", "subject_id", "passing", "session_type"]], on="session_id")
code_probes = shield.read_code_csv("metadata_tables/probes_table.csv")
cpk = code_units[["session", "probe_id"]].drop_duplicates().merge(code_probes[["id", "name"]], left_on="probe_id", right_on="id")
code_keys = set(zip(cpk["session"], cpk["name"]))
ins["in_code"] = [k in code_keys for k in zip(ins["session_id"], ins["probe"])]
for lab, g in (("87 QC-passing sessions", ins[ins["passing"]]), ("all 99 sessions", ins)):
    print(f"{lab}: mice {g['subject_id'].nunique()}, experiments {g['session_id'].nunique()}, insertions {len(g)}; "
          f"of these also in the code's units table {int(g['in_code'].sum())}")
pa_ = ins[ins["passing"] & ~ins["in_code"]]
print(f"insertions of the 87 sessions not in the code's table: {len(pa_)}, with {int(pa_['n'].sum())} units, labelled",
      units.set_index(["session_id", "probe"]).loc[list(zip(pa_["session_id"], pa_["probe"])), "structure_acronym"].value_counts().to_dict())
print("probes in the code's probes table for the 87 sessions:",
      int(shield.read_code_csv("metadata_tables/probes_table.csv").query("session in @PASS").shape[0]))

87 QC-passing sessions: mice 25, experiments 87, insertions 507; of these also in the code's units table 467
all 99 sessions: mice 27, experiments 99, insertions 578; of these also in the code's units table 533
insertions of the 87 sessions not in the code's table: 40, with 2004 units, labelled {'Track not annotated': 1873, 'CP': 92, 'VL': 22, 'V3': 14, 'No Area': 2, 'MOp': 1}

probes in the code's probes table for the 87 sessions: 467

for lab, g in (("(a) 87 QC-passing sessions", ins[ins["passing"]]), ("(b) all 99 sessions", ins),
               ("(a), insertions in the code's table only", ins[ins["passing"] & ins["in_code"]])):
    per = g.groupby("session_id").size()
    print(f"{lab}: {len(g)} insertions / {per.size} sessions = {len(g) / per.size:.3f} probes per session; "
          f"sessions with 6, 5, 4 or fewer probes: {(per == 6).sum()}, {(per == 5).sum()}, {(per <= 4).sum()}")
for lab, q in (("87 sessions", code_probes[code_probes["session"].isin(PASS)]), ("99 sessions", code_probes)):
    print(f"code's probes table, {lab}: {len(q)} probes in {q['session'].nunique()} sessions, inserted {q['inserted'].value_counts().to_dict()}, "
          f"{len(q) / q['session'].nunique():.3f} per session")

(a) 87 QC-passing sessions: 507 insertions / 87 sessions = 5.828 probes per session; sessions with 6, 5, 4 or fewer probes: 77, 8, 2
(b) all 99 sessions: 578 insertions / 99 sessions = 5.838 probes per session; sessions with 6, 5, 4 or fewer probes: 88, 9, 2
(a), insertions in the code's table only: 467 insertions / 87 sessions = 5.368 probes per session; sessions with 6, 5, 4 or fewer probes: 54, 24, 9
code's probes table, 87 sessions: 467 probes in 87 sessions, inserted {True: 467}, 5.368 per session
code's probes table, 99 sessions: 533 probes in 99 sessions, inserted {True: 533}, 5.384 per session

ia = units[units["passing"] & units["in_area"]]
nprobe = ia.groupby("session_id")["probe"].nunique()
six = nprobe[nprobe == 6].index
d6 = ia[ia["session_id"].isin(six)]
areas = d6.groupby("session_id")["structure_acronym"].nunique()
nunits = d6[d6["good"]].groupby("session_id").size()
print(f"6-probe sessions: {len(six)} of {nprobe.size}")
print(f"areas per session: mean {areas.mean():.2f}, SD {areas.std():.2f} (range {areas.min()} to {areas.max()})")
print(f"non-noise units per session: mean {nunits.mean():.1f}, SD {nunits.std():.1f} (range {nunits.min()} to {nunits.max()})")
by_mouse = d6[d6["good"]].groupby(["subject_id", "session_id"]).size().groupby("subject_id").mean()
print("mean non-noise units per 6-probe session, by mouse:", by_mouse.round(0).astype(int).to_dict())
d6c = d6[~d6["session_id"].isin(bad["session_id"])]
print(f"leaving out the 6 sessions whose units differ from the code's: {d6c['session_id'].nunique()} sessions, units mean "
      f"{d6c[d6c['good']].groupby('session_id').size().mean():.1f}, areas mean {d6c.groupby('session_id')['structure_acronym'].nunique().mean():.2f}")

6-probe sessions: 44 of 77
areas per session: mean 28.11, SD 5.12 (range 19 to 41)
non-noise units per session: mean 1821.0, SD 341.9 (range 1100 to 2553)
mean non-noise units per 6-probe session, by mouse: {'607660': 2007, '608672': 1844, '612090': 2553, '614547': 2036, '614608': 2298, '615047': 1456, '615048': 2072, '615563': 1856, '615564': 1640, '623322': 1598, '623784': 1463, '623786': 1693, '626279': 1607, '632295': 1935, '636740': 1723, '637483': 1916, '640887': 1734, '640890': 2064}
leaving out the 6 sessions whose units differ from the code's: 44 sessions, units mean 1821.0, areas mean 28.11

for lab, g in (("(a) 87 QC-passing sessions", ins[ins["passing"]]), ("(b) all sessions", ins),
               ("(a), insertions in the code's table only", ins[ins["passing"] & ins["in_code"]])):
    per = g.groupby("subject_id").size()
    print(f"{lab}: {per.size} mice, insertions per mouse mean {per.mean():.2f}, SD {per.std():.2f}, "
          f"range {per.min()} to {per.max()}")
print("sessions per mouse (all sessions):", sessions.groupby("subject_id").size().value_counts().sort_index().to_dict())

(a) 87 QC-passing sessions: 25 mice, insertions per mouse mean 20.28, SD 3.34, range 12 to 24
(b) all sessions: 27 mice, insertions per mouse mean 21.41, SD 3.33, range 12 to 24
(a), insertions in the code's table only: 25 mice, insertions per mouse mean 18.68, SD 4.18, range 8 to 24
sessions per mouse (all sessions): {2: 1, 3: 7, 4: 19}

dg_good = units[units["passing"] & units["good"]]
vbn = vbn_units.merge(vbn_sessions[["ecephys_session_id", "session_number"]], on="ecephys_session_id")
vbn_good = vbn[vbn["quality"] == "good"]
print("non-noise units: SHIELD", len(dg_good), "(paper 143,828); VBN", len(vbn_good), "(paper 262,177)")
def boot_sem_median(x, rng, n=1000):
    return np.std([np.median(rng.choice(x, len(x))) for _ in range(n)])
rng = np.random.default_rng(4)
x, y = dg_good["max_drift"].dropna().values, vbn_good["max_drift"].dropna().values
print(f"median max_drift: SHIELD {np.median(x):.2f} (bootstrap SE {boot_sem_median(x, rng, 200):.3f}), "
      f"VBN {np.median(y):.2f} (SE {boot_sem_median(y, rng, 200):.3f}); difference {np.median(x) - np.median(y):.2f} um")
print(f"two-sided Mann-Whitney p = {stats.mannwhitneyu(x, y).pvalue:.2e}")
xc = dg_good.loc[~dg_good["session_id"].isin(bad["session_id"]), "max_drift"]
in_code_keys = set(zip(ins.loc[ins["in_code"], "session_id"], ins.loc[ins["in_code"], "probe"]))
dg_in = np.array([k in in_code_keys for k in zip(dg_good["session_id"], dg_good["probe"])])
print("non-noise units in the code's units table for the 87 sessions:", int((cu["quality"] == "good").sum()),
      "; SHIELD non-noise units on insertions in the code's table:", int(dg_in.sum()), ", on the 40 others:", int((~dg_in).sum()))
print("non-noise units outside the 6 sessions that differ: dandiset", int((~dg_good["session_id"].isin(bad["session_id"])).sum()),
      ", code", int(((cu["quality"] == "good") & ~cu["session"].isin(bad["session_id"])).sum()))
print(f"SHIELD median without the 6 sessions that differ from the code: {xc.median():.2f}")
print("SHIELD median by mouse:", dg_good.groupby("subject_id")["max_drift"].median().round(1).to_dict())

non-noise units: SHIELD 144030 (paper 143,828); VBN 262177 (paper 262,177)

median max_drift: SHIELD 25.05 (bootstrap SE 0.063), VBN 22.82 (SE 0.043); difference 2.23 um
two-sided Mann-Whitney p = 3.92e-292
non-noise units in the code's units table for the 87 sessions: 144068 ; SHIELD non-noise units on insertions in the code's table: 143828 , on the 40 others: 202
non-noise units outside the 6 sessions that differ: dandiset 135334 , code 135334
SHIELD median without the 6 sessions that differ from the code: 25.22
SHIELD median by mouse: {'607660': 27.4, '608671': 23.7, '608672': 29.2, '611160': 29.7, '611166': 26.4, '612090': 27.4, '614547': 21.8, '614608': 20.2, '615047': 25.8, '615048': 29.2, '615563': 26.6, '615564': 24.5, '623319': 30.6, '623322': 24.4, '623784': 27.7, '623786': 32.3, '626279': 25.7, '632295': 21.8, '636740': 21.9, '637483': 20.5, '637488': 24.4, '638387': 27.2, '640887': 20.8, '640890': 19.6, '642504': 20.5}

metrics = ["presence_ratio", "amplitude_cutoff", "isi_violations", "snr", "isolation_distance", "max_drift", "amplitude"]
rows = []
for m in metrics:
    a, b = dg_good[m].dropna().values, vbn_good[m].dropna().values
    rows.append({"metric": m, "SHIELD_median": np.median(a), "VBN_median": np.median(b),
                 "ratio": np.median(a) / np.median(b), "within_10pct": abs(np.median(a) / np.median(b) - 1) <= 0.10,
                 "mannwhitney_p": stats.mannwhitneyu(a, b).pvalue,
                 "rank_biserial": 1 - 2 * stats.mannwhitneyu(a, b).statistic / (len(a) * len(b))})
qm = pd.DataFrame(rows)
print(qm.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("(a) metrics not significant at 0.05:", int((qm["mannwhitney_p"] >= 0.05).sum()), "of 7;",
      " (b) metrics with medians within 10%:", int(qm["within_10pct"].sum()), "of 7")
fig, ax = plt.subplots(1, 7, figsize=(16, 3))
for i, m in enumerate(metrics):
    ax[i].boxplot([dg_good[m].dropna(), vbn_good[m].dropna()], showfliers=False)
    ax[i].set_xticklabels(["SHIELD", "VBN"]); ax[i].set_title(m, fontsize=8)
plt.tight_layout(); plt.show()
print("Boxes without outliers; whiskers at 1.5 IQR.")

            metric  SHIELD_median  VBN_median  ratio  within_10pct  mannwhitney_p  rank_biserial
    presence_ratio           0.99        0.99      1          True      2.09e-157        0.04061
  amplitude_cutoff        0.04408     0.04947  0.891         False      3.233e-16        0.01543
    isi_violations        0.06981     0.08875 0.7866         False     4.134e-232        0.06143
               snr          2.471       2.497 0.9898          True      3.861e-06       0.008746
isolation_distance          56.11       58.22 0.9638          True     2.991e-111         0.0426
         max_drift          25.05       22.82  1.098          True     3.918e-292       -0.06948
         amplitude          125.2       130.2 0.9616          True     3.739e-161        0.05122
(a) metrics not significant at 0.05: 0 of 7;  (b) metrics with medians within 10%: 5 of 7

Boxes without outliers; whiskers at 1.5 IQR.

vq = (vbn_good["isi_violations"] < 0.5) & (vbn_good["amplitude_cutoff"] < 0.1) & (vbn_good["presence_ratio"] > 0.95)
print(f"SHIELD: {dg_good['qc'].sum()} of {len(dg_good)} non-noise units pass = {100 * dg_good['qc'].mean():.2f}%")
print(f"VBN: {vq.sum()} of {len(vbn_good)} = {100 * vq.mean():.2f}%")
print("SHIELD per mouse (%):", (100 * dg_good.groupby("subject_id")["qc"].mean()).round(1).to_dict())
print(f"SHIELD without the 6 differing sessions: {100 * dg_good.loc[~dg_good['session_id'].isin(bad['session_id']), 'qc'].mean():.2f}%")

SHIELD: 59188 of 144030 non-noise units pass = 41.09%
VBN: 106680 of 262177 = 40.69%
SHIELD per mouse (%): {'607660': 38.3, '608671': 39.8, '608672': 39.1, '611160': 42.2, '611166': 41.8, '612090': 39.8, '614547': 44.9, '614608': 43.8, '615047': 40.9, '615048': 40.7, '615563': 41.8, '615564': 38.6, '623319': 37.9, '623322': 40.0, '623784': 41.3, '623786': 37.6, '626279': 42.6, '632295': 42.1, '636740': 42.4, '637483': 43.7, '637488': 39.4, '638387': 39.7, '640887': 41.7, '640890': 46.5, '642504': 43.7}
SHIELD without the 6 differing sessions: 40.93%

yc = dg_good.groupby(["subject_id", "probe", "session_type"]).size().rename("units").reset_index()
days = ["EPHYS_1", "EPHYS_2", "EPHYS_3", "EPHYS_4"]
groups = [yc.loc[yc["session_type"] == d, "units"].values for d in days]
for d, g in zip(days, groups):
    print(f"{d}: {len(g)} insertions, median {np.median(g):.1f} non-noise units (IQR {np.percentile(g, 25):.0f} to {np.percentile(g, 75):.0f})")
kw = stats.kruskal(*groups)
print(f"Kruskal-Wallis across the four days: H = {kw.statistic:.3f}, p = {kw.pvalue:.3f}")
dg_code = dg_good[[k in code_keys for k in zip(dg_good["session_id"], dg_good["probe"])]]
yca = dg_code.groupby(["subject_id", "probe", "session_type"]).size().rename("units").reset_index()
ga = [yca.loc[yca["session_type"] == d, "units"].values for d in days]
print("insertions in the code's table only: insertions", [len(g) for g in ga], ", medians", [float(np.median(g)) for g in ga],
      f", Kruskal-Wallis p = {stats.kruskal(*ga).pvalue:.3f}")
print("(paper: 129, 126, 121, 91 insertions; p = 0.48; the code printed medians 314, 305.5, 295, 303)")
ycc = dg_good[~dg_good["session_id"].isin(bad["session_id"])].groupby(["subject_id", "probe", "session_type"]).size().rename("units").reset_index()
kw2 = stats.kruskal(*[ycc.loc[ycc["session_type"] == d, "units"].values for d in days])
print(f"without the 6 differing sessions: p = {kw2.pvalue:.3f}, medians",
      [float(np.median(ycc.loc[ycc['session_type'] == d, 'units'])) for d in days])

EPHYS_1: 136 insertions, median 308.0 non-noise units (IQR 224 to 378)
EPHYS_2: 135 insertions, median 302.0 non-noise units (IQR 213 to 374)
EPHYS_3: 128 insertions, median 287.5 non-noise units (IQR 222 to 358)
EPHYS_4: 101 insertions, median 277.0 non-noise units (IQR 193 to 365)
Kruskal-Wallis across the four days: H = 3.409, p = 0.333
insertions in the code's table only: insertions [129, 126, 121, 91] , medians [314.0, 305.5, 295.0, 303.0] , Kruskal-Wallis p = 0.477
(paper: 129, 126, 121, 91 insertions; p = 0.48; the code printed medians 314, 305.5, 295, 303)
without the 6 differing sessions: p = 0.282, medians [304.0, 304.5, 287.5, 275.0]

vy = vbn_good.groupby(["session_number", "ecephys_session_id", "ecephys_probe_id"]).size().rename("units").reset_index()
v1, v2 = [vy.loc[vy["session_number"] == k, "units"].values for k in (1, 2)]
print(f"VBN day 1: {len(v1)} insertions, median {np.median(v1):.1f}; day 2: {len(v2)} insertions, median {np.median(v2):.1f}")
ka = stats.kruskal(v1, v2)
kb = stats.kruskal(v1, v2, *groups)
print(f"(a) Kruskal-Wallis VBN day 1 against day 2: H = {ka.statistic:.2f}, p = {ka.pvalue:.2e}")
print(f"(b) Kruskal-Wallis over the code's six groups (two VBN days, four SHIELD days): H = {kb.statistic:.2f}, p = {kb.pvalue:.2e}")

VBN day 1: 459 insertions, median 321.0; day 2: 446 insertions, median 261.0
(a) Kruskal-Wallis VBN day 1 against day 2: H = 105.07, p = 1.18e-24
(b) Kruskal-Wallis over the code's six groups (two VBN days, four SHIELD days): H = 98.30, p = 1.20e-19

fr = dg_good.groupby(["session_type", "subject_id", "session_id"])["firing_rate"].median().rename("fr").reset_index()
fg = [fr.loc[fr["session_type"] == d, "fr"].values for d in days]
for d, g in zip(days, fg):
    print(f"{d}: {len(g)} sessions, median of session medians {np.median(g):.3f} Hz (range {g.min():.2f} to {g.max():.2f})")
k = stats.kruskal(*fg)
print(f"Kruskal-Wallis: H = {k.statistic:.3f}, p = {k.pvalue:.3f} (text 0.21; code and figure 0.26)")

EPHYS_1: 23 sessions, median of session medians 3.154 Hz (range 2.11 to 4.39)
EPHYS_2: 24 sessions, median of session medians 3.176 Hz (range 2.45 to 4.36)
EPHYS_3: 22 sessions, median of session medians 3.468 Hz (range 2.40 to 6.32)
EPHYS_4: 18 sessions, median of session medians 3.107 Hz (range 2.34 to 4.41)
Kruskal-Wallis: H = 4.385, p = 0.223 (text 0.21; code and figure 0.26)

v1u = units[units["passing"] & units["qc"] & (units["isi_label"] == "VISp")].join(rf[["p_rf"]], how="left")
print("QC VISp units:", len(v1u), "; with an RF p value:", int(v1u["p_rf"].notna().sum()))
v1u = v1u.join(code_units.set_index("unit_id")[["p_value_rf"]], how="left")
res = {}
for thr in (0.01, 0.05):
    for src in ("p_rf", "p_value_rf"):
        frac = v1u.assign(sig=v1u[src] < thr).groupby(["session_type", "subject_id", "session_id"])["sig"].mean().reset_index()
        g = [frac.loc[frac["session_type"] == d, "sig"].values for d in days]
        k = stats.kruskal(*g)
        res[(thr, src)] = k.pvalue
        print(f"threshold {thr}, p values {'computed here' if src == 'p_rf' else 'of the code table'}: sessions per day "
              f"{[len(x) for x in g]}, median fraction {[round(float(np.median(x)), 3) for x in g]}, Kruskal-Wallis p = {k.pvalue:.3f}")

QC VISp units: 2469 ; with an RF p value: 2469
threshold 0.01, p values computed here: sessions per day [12, 11, 12, 11], median fraction [0.9, 0.877, 0.913, 0.9], Kruskal-Wallis p = 0.671
threshold 0.01, p values of the code table: sessions per day [12, 11, 12, 11], median fraction [0.9, 0.877, 0.904, 0.865], Kruskal-Wallis p = 0.864
threshold 0.05, p values computed here: sessions per day [12, 11, 12, 11], median fraction [0.907, 0.896, 0.931, 0.964], Kruskal-Wallis p = 0.368
threshold 0.05, p values of the code table: sessions per day [12, 11, 12, 11], median fraction [0.915, 0.896, 0.923, 0.944], Kruskal-Wallis p = 0.813

ps = pd.concat([v["im115_psth"] for v in SP.values()])
sel = units[units["passing"] & units["qc"] & (units["isi_label"] == "VISp")].join(rf[["p_rf"]], how="inner")
sel = sel[sel["p_rf"] < 0.01]
P = ps.loc[sel.index].values
t_ms = np.arange(2500) - 1000
base = P[:, 990:1000].mean(1, keepdims=True)
resp = (P[:, 1000:1250] - base).mean(1)
peak = np.argmax(P[:, 1000:1250], axis=1)
summ = []
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5))
for d in days:
    m = (sel["session_type"] == d).values
    mean = P[m].mean(0); mb = mean - mean[990:1000].mean(); mn = mb / mb.max()
    ax[0].plot(t_ms, mb, label=d); ax[1].plot(t_ms, mn, label=d)
    summ.append({"day": d, "units": int(m.sum()), "sessions": sel.loc[m, "session_id"].nunique(), "mice": sel.loc[m, "subject_id"].nunique(),
                 "median_response_0_250ms_Hz": np.median(resp[m]), "median_peak_ms": np.median(peak[m]),
                 "population_peak_ms": int(np.argmax(mn[1000:1250]))})
for a in ax: a.set_xlim(-100, 400); a.legend(fontsize=7); a.set_xlabel("ms from im115_r onset")
ax[0].set_ylabel("baseline-subtracted rate (Hz)"); ax[1].set_ylabel("normalised")
plt.show()
summ = pd.DataFrame(summ); print(summ.round(2).to_string(index=False))
ka = stats.kruskal(*[resp[(sel["session_type"] == d).values] for d in days])
kb = stats.kruskal(*[peak[(sel["session_type"] == d).values] for d in days])
print(f"(a) mean baseline-subtracted response 0-250 ms: Kruskal-Wallis p = {ka.pvalue:.3g}")
print(f"(b) time of peak 0-250 ms: Kruskal-Wallis p = {kb.pvalue:.3g}")
day_means = np.array([P[(sel["session_type"] == d).values].mean(0)[1000:1250] for d in days])
print("Pearson r between day population responses (0-250 ms):", np.round(np.corrcoef(day_means), 3).tolist())
pre, post = P[:, 900:1000].mean(1), P[:, 1020:1120].mean(1)
print(f"clock check: VISp rate -100 to 0 ms {np.median(pre):.2f} Hz, 20 to 120 ms {np.median(post):.2f} Hz "
      f"(Wilcoxon p = {stats.wilcoxon(post, pre).pvalue:.1e}); sessions where the post-onset rate is higher:",
      int((sel.assign(d=post - pre).groupby('session_id')['d'].median() > 0).sum()), "of", sel["session_id"].nunique())

    day  units  sessions  mice  median_response_0_250ms_Hz  median_peak_ms  population_peak_ms
EPHYS_1    595        12    12                        0.82            53.0                  42
EPHYS_2    482        11    11                        1.51            48.0                  44
EPHYS_3    534        12    12                        1.19            48.0                  44
EPHYS_4    446        11    11                        0.83            49.0                  42
(a) mean baseline-subtracted response 0-250 ms: Kruskal-Wallis p = 0.343
(b) time of peak 0-250 ms: Kruskal-Wallis p = 0.7
Pearson r between day population responses (0-250 ms): [[1.0, 0.991, 0.983, 0.98], [0.991, 1.0, 0.982, 0.977], [0.983, 0.982, 1.0, 0.983], [0.98, 0.977, 0.983, 1.0]]
clock check: VISp rate -100 to 0 ms 4.05 Hz, 20 to 120 ms 5.66 Hz (Wilcoxon p = 1.2e-41); sessions where the post-onset rate is higher: 39 of 46

for lab, sset in (("87 QC-passing sessions", sessions[sessions["passing"]]), ("all sessions (variant)", sessions)):
    for col, name in (("n_rewards", "(a) all rewards"), ("n_rewards_earned", "(b) earned rewards")):
        g = [sset.loc[sset["session_type"] == d, col].values for d in days]
        k = stats.kruskal(*g)
        print(f"{lab}, {name}: sessions {[len(x) for x in g]}, medians {[float(np.median(x)) for x in g]}, Kruskal-Wallis p = {k.pvalue:.3f}")

87 QC-passing sessions, (a) all rewards: sessions [23, 24, 22, 18], medians [73.0, 65.0, 61.0, 53.5], Kruskal-Wallis p = 0.731
87 QC-passing sessions, (b) earned rewards: sessions [23, 24, 22, 18], medians [67.0, 59.0, 55.0, 48.0], Kruskal-Wallis p = 0.736
all sessions (variant), (a) all rewards: sessions [24, 26, 25, 24], medians [68.0, 63.5, 60.0, 56.5], Kruskal-Wallis p = 0.817
all sessions (variant), (b) earned rewards: sessions [24, 26, 25, 24], medians [62.5, 57.5, 54.0, 50.5], Kruskal-Wallis p = 0.816

def sat(t, a, b, tau):
    return a - b * np.exp(-t / tau)
curves = {k: v["s4i_counts"] for k, v in SP.items() if k in PASS and len(v["s4i_counts"]) > 600}
L = min(len(c) for c in curves.values())
C = np.array([c[:L] for c in curves.values()], float)
tmin = np.arange(L) / 60
print(f"sessions: {len(C)}; common length {L / 60:.1f} min from the first stimulus")
fits = {}
fig, ax = plt.subplots(figsize=(6, 3))
for lab, norm in (("(a) divided by the session mean", C.mean(1, keepdims=True)),
                  ("(b) divided by the mean of the last 10 min", C[:, -600:].mean(1, keepdims=True))):
    y = (C / norm).mean(0)
    ym = y.reshape(-1, 60)[: L // 60].mean(1) if L % 60 == 0 else y[: (L // 60) * 60].reshape(-1, 60).mean(1)
    tm = np.arange(len(ym)) + 0.5
    popt, pcov = curve_fit(sat, tm, ym, p0=[1, 0.3, 20], maxfev=20000)
    fits[lab] = popt
    ax.plot(tm, ym, label=lab[:3])
    print(f"{lab}: tau = {popt[2]:.1f} min (SE {np.sqrt(pcov[2, 2]):.1f}), a = {popt[0]:.3f}, b = {popt[1]:.3f}; "
          f"rate at minute 1 {ym[0]:.3f}, at minute 30 {ym[29]:.3f}, last minute {ym[-1]:.3f}")
ax.set_xlabel("minutes from first stimulus"); ax.set_ylabel("normalised rate"); ax.legend(); plt.show()
# the same per session, to show the spread
taus = []
for c in C:
    ym = (c / c.mean())[: (L // 60) * 60].reshape(-1, 60).mean(1)
    try:
        taus.append(curve_fit(sat, np.arange(len(ym)) + 0.5, ym, p0=[1, 0.3, 20], maxfev=20000)[0][2])
    except Exception:
        taus.append(np.nan)
taus = np.array(taus)
print(f"per-session fits (a): median tau {np.nanmedian(taus):.1f} min, IQR {np.nanpercentile(taus, 25):.1f} to {np.nanpercentile(taus, 75):.1f}, failed {int(np.isnan(taus).sum())}")
first = sessions.set_index("session_id").loc[list(curves), "first_stim"]
print(f"first stimulus at {first.min():.0f} to {first.max():.0f} s of the recording (median {first.median():.0f})")

# the fit window: the fits above use the common length (77.3 min, which runs past the active block
# into the Gabor and flash epochs); the same fit over shorter windows from the first stimulus,
# including the shortest active block, and each window normalised by its own mean
be = np.array([(SP[k]["beh_end"] - SP[k]["first_stim"]) / 60 for k in curves])
bs = np.array([(SP[k]["beh_start"] - SP[k]["first_stim"]) / 60 for k in curves])
print(f"active block from {bs.min():.1f} to {bs.max():.1f} min after the first stimulus (median {np.median(bs):.1f}), "
      f"ending at {be.min():.1f} to {be.max():.1f} min (median {np.median(be):.1f})")
win_rows = []
for W in (30, 40, 45, 50, 60, int(be.min()), L // 60):
    Cw = C[:, : W * 60]
    for lab, normw in (("a", Cw.mean(1, keepdims=True)), ("b", Cw[:, -600:].mean(1, keepdims=True))):
        ym = (Cw / normw).mean(0).reshape(-1, 60).mean(1)
        tm = np.arange(len(ym)) + 0.5
        try:
            popt, pcov = curve_fit(sat, tm, ym, p0=[1, 0.3, 20], maxfev=20000)
            tau, se = popt[2], np.sqrt(pcov[2, 2])
        except Exception:
            tau, se = np.nan, np.nan
        win_rows.append({"window_min": W, "reading": lab, "tau_min": tau, "SE": se, "within_22.5_33": 22.5 <= tau <= 33})
print(pd.DataFrame(win_rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
# per-session fits over the active block of each session (from its first stimulus to its end, at most
# the common length)
pts = []
for k, c in zip(curves, C):
    W = min(int((SP[k]["beh_end"] - SP[k]["first_stim"]) // 60), L // 60)
    ym = (c[: W * 60] / c[: W * 60].mean()).reshape(-1, 60).mean(1)
    try:
        pts.append(curve_fit(sat, np.arange(len(ym)) + 0.5, ym, p0=[1, 0.3, 20], maxfev=20000)[0][2])
    except Exception:
        pts.append(np.nan)
pts = np.array(pts)
print(f"per-session fits over each session's active block (a): median tau {np.nanmedian(pts):.1f} min, "
      f"IQR {np.nanpercentile(pts, 25):.1f} to {np.nanpercentile(pts, 75):.1f}, sessions within 22.5 to 33: "
      f"{int(((pts >= 22.5) & (pts <= 33)).sum())} of {len(pts)}, failed {int(np.isnan(pts).sum())}")

sessions: 87; common length 77.3 min from the first stimulus
(a) divided by the session mean: tau = 36.0 min (SE 2.3), a = 1.147, b = 0.358; rate at minute 1 0.826, at minute 30 0.999, last minute 1.077
(b) divided by the mean of the last 10 min: tau = 35.4 min (SE 2.3), a = 1.053, b = 0.322; rate at minute 1 0.763, at minute 30 0.921, last minute 0.990

per-session fits (a): median tau 34.8 min, IQR 23.1 to 52.5, failed 0
first stimulus at 25 to 636 s of the recording (median 64)
active block from 0.0 to 25.6 min after the first stimulus (median 0.0), ending at 65.0 to 90.7 min (median 65.1)
 window_min reading  tau_min       SE  within_22.5_33
         30       a 2.32e+05  4.3e+07           False
         30       b 9.23e+04 2.01e+07           False
         40       a     79.1       27           False
         40       b     66.9     19.8           False
         45       a       84     22.7           False
         45       b     76.2     19.4           False
         50       a     69.6     12.2           False
         50       b     65.2       11           False
         60       a     47.1      4.3           False
         60       b     45.9      4.2           False
         64       a     42.2     3.34           False
         64       b     41.2     3.26           False
         77       a       36      2.3           False
         77       b     35.4     2.27           False
per-session fits over each session's active block (a): median tau 39.8 min, IQR 25.7 to 66.8, sessions within 22.5 to 33: 26 of 87, failed 0

COMBOS = [(L, E, R, S) for L in ("L1", "L2") for E in ("E1", "E2") for R in ("R1", "R2") for S in ("a", "b")]
SETS = {"a": PASS, "b": set(sessions["session_id"])}

def events_table(ref, L, E, S):
    rows = []
    for sid, v in SP.items():
        if sid not in SETS[S]:
            continue
        e = v["alpha"].get((L, ref, E))
        if e is None or not len(e["events"]):
            continue
        ev = e["events"].copy()
        ev["session_id"] = sid
        rows.append(ev)
    return pd.concat(rows, ignore_index=True) if rows else pd.DataFrame()

det = []
for L in ("L1", "L2"):
    for E in ("E1", "E2"):
        for S in ("a", "b"):
            for ref in shield.REFERENCE_AREAS:
                ev = events_table(ref, L, E, S)
                elig = sum(1 for sid, v in SP.items() if sid in SETS[S] and ref in v["areas"][L])
                q = ev[ev["qual"]] if len(ev) else ev
                det.append({"L": L, "E": E, "set": S, "area": ref, "eligible_sessions": elig,
                            "sessions_with_events": ev["session_id"].nunique() if len(ev) else 0,
                            "events": len(ev), "qualifying": len(q),
                            "qualifying_smoothed_acg": int((ev["auto_power_smooth"] > shield.AUTO_THRESHOLD).sum()) if len(ev) else 0,
                            "sessions_over_5_qualifying": int((q.groupby("session_id").size() > 5).sum()) if len(q) else 0,
                            "qualifying_in_those": int(q.groupby("session_id").size().loc[lambda s: s > 5].sum()) if len(q) else 0})
det = pd.DataFrame(det)
print(det.to_string(index=False))
print("(the code printed 50 sessions with 1653 qualifying VISp events and 14 sessions with 113 MOp events; not values of the paper)")

# area labels under L1 and L2 (87 sessions): sessions with at least the minimum of QC units, median units
from collections import Counter
for L in ("L1", "L2"):
    cnt, nu = Counter(), {}
    for sid, v in SP.items():
        if sid in PASS:
            for a, n in v["areas"][L].items():
                cnt[a] += 1; nu.setdefault(a, []).append(n)
    print(f"\n{L} areas (sessions, median units):", {a: (cnt[a], int(np.median(nu[a]))) for a in sorted(cnt)})

 L  E set area  eligible_sessions  sessions_with_events  events  qualifying  qualifying_smoothed_acg  sessions_over_5_qualifying  qualifying_in_those
L1 E1   a VISp                 47                    38    2965        1518                     2840                          11                 1483
L1 E1   a  MOp                 38                     8     106          80                       98                           4                   68
L1 E1   b VISp                 51                    42    3011        1522                     2872                          11                 1483
L1 E1   b  MOp                 44                    10     119          82                      110                           4                   68
L1 E2   a VISp                 47                    41    3581        1694                     3364                          11                 1651
L1 E2   a  MOp                 38                    12     162         119                      152                           5                  108
L1 E2   b VISp                 51                    45    3665        1704                     3426                          11                 1651
L1 E2   b  MOp                 44                    14     182         122                      170                           5                  108
L2 E1   a VISp                 51                    46    2673        1469                     2554                          11                 1428
L2 E1   a  MOp                 38                     8     106          80                       98                           4                   68
L2 E1   b VISp                 56                    51    2737        1478                     2603                          11                 1428
L2 E1   b  MOp                 44                    10     119          82                      110                           4                   68
L2 E2   a VISp                 51                    47    3244        1643                     3035                          13                 1604
L2 E2   a  MOp                 38                    12     162         119                      152                           5                  108
L2 E2   b VISp                 56                    52    3355        1660                     3120                          14                 1611
L2 E2   b  MOp                 44                    14     182         122                      170                           5                  108
(the code printed 50 sessions with 1653 qualifying VISp events and 14 sessions with 113 MOp events; not values of the paper)

L1 areas (sessions, median units): {'ACAd': (25, 27), 'ACAv': (19, 55), 'ACB': (3, 32), 'APN': (17, 34), 'AV': (1, 20), 'BST': (1, 28), 'CA1': (56, 35), 'CA3': (24, 31), 'CL': (5, 32), 'CP': (42, 52), 'DG': (31, 33), 'DP': (8, 31), 'ENTm': (1, 38), 'Eth': (1, 21), 'FRP': (1, 25), 'ILA': (8, 38), 'LD': (14, 39), 'LGd': (20, 38), 'LGv': (5, 27), 'LP': (26, 52), 'LSc': (17, 29), 'LSr': (30, 56), 'LSv': (4, 35), 'MB': (14, 31), 'MD': (5, 28), 'MGd': (6, 35), 'MGm': (6, 35), 'MGv': (14, 42), 'MOp': (38, 64), 'MOs': (46, 65), 'MRN': (5, 32), 'MS': (1, 28), 'NOT': (1, 21), 'OLF': (4, 31), 'ORBm': (2, 49), 'ORBvl': (3, 119), 'PAR': (1, 38), 'PF': (2, 24), 'PIL': (1, 36), 'PL': (12, 45), 'PO': (17, 51), 'POL': (3, 44), 'POST': (5, 30), 'PPT': (3, 22), 'PVT': (1, 48), 'ProS': (8, 37), 'RSPagl': (26, 42), 'RSPd': (28, 37), 'RSPv': (24, 53), 'RT': (10, 24), 'SCig': (7, 29), 'SCiw': (6, 28), 'SCsg': (1, 28), 'SF': (4, 26), 'SFO': (1, 25), 'SGN': (8, 24), 'SSp': (61, 62), 'SSs': (2, 99), 'STR': (7, 27), 'SUB': (20, 35), 'TH': (4, 32), 'TRS': (3, 30), 'TTd': (3, 41), 'Track not annotated': (21, 267), 'VAL': (2, 22), 'VIS': (27, 50), 'VISa': (3, 29), 'VISal': (5, 54), 'VISam': (6, 44), 'VISl': (2, 43), 'VISlm': (15, 52), 'VISmma': (2, 41), 'VISp': (47, 50), 'VISpm': (5, 55), 'VISrl': (11, 48), 'VISrll': (1, 60), 'VL': (1, 26), 'VPL': (6, 27), 'VPM': (11, 41), 'ZI': (1, 29)}

L2 areas (sessions, median units): {'ACAd': (21, 28), 'ACAv': (17, 57), 'ACB': (3, 32), 'APN': (17, 34), 'BST': (1, 28), 'CA1': (56, 35), 'CA3': (22, 31), 'CL': (4, 36), 'CP': (41, 55), 'DG': (30, 33), 'DP': (7, 31), 'ENTm': (1, 38), 'Eth': (1, 21), 'FRP': (1, 25), 'ILA': (8, 38), 'LD': (13, 39), 'LGd': (20, 38), 'LGv': (5, 27), 'LP': (24, 56), 'LSc': (15, 32), 'LSr': (29, 58), 'LSv': (4, 35), 'MB': (14, 31), 'MD': (4, 39), 'MGd': (6, 35), 'MGm': (5, 36), 'MGv': (14, 42), 'MOp': (38, 64), 'MOs': (46, 65), 'MRN': (4, 48), 'MS': (1, 28), 'NOT': (1, 21), 'OLF': (4, 31), 'ORBm': (2, 49), 'ORBvl': (3, 119), 'PAR': (1, 38), 'PF': (2, 24), 'PIL': (1, 36), 'PL': (11, 47), 'PO': (17, 63), 'POL': (3, 44), 'POST': (5, 30), 'PPT': (3, 22), 'PVT': (1, 48), 'ProS': (8, 37), 'RSPagl': (24, 46), 'RSPd': (28, 39), 'RSPv': (26, 52), 'RT': (6, 34), 'SCig': (7, 29), 'SCiw': (6, 28), 'SCsg': (1, 28), 'SF': (4, 26), 'SFO': (1, 25), 'SGN': (8, 24), 'SSp': (62, 64), 'SSs': (2, 99), 'STR': (6, 30), 'SUB': (20, 35), 'TEa': (1, 69), 'TH': (4, 32), 'TRS': (3, 30), 'TTd': (3, 41), 'Track not annotated': (21, 267), 'VAL': (2, 22), 'VISa': (5, 37), 'VISal': (14, 48), 'VISam': (8, 40), 'VISl': (18, 46), 'VISli': (4, 41), 'VISp': (51, 55), 'VISpm': (9, 51), 'VISpor': (1, 28), 'VISrl': (11, 47), 'VL': (1, 26), 'VPL': (6, 27), 'VPM': (11, 41), 'ZI': (1, 29)}

NOT_FAMILY = ["CP-RF", "CP-noRF", "VIS", *shield.NOT_IN_AREA]

def pair_table(ref, L, E, R, S, min_events=25):
    """The Figure 7F pair test: per pair the per-event band coherence minus random, pooled over
    sessions with more than 5 qualifying events; pairs with at least 25 events; Wilcoxon against
    zero; Benjamini-Hochberg across the reference area's regions. The family holds brain areas
    only: CP-RF and CP-noRF (Figure 8F groups), the unannotated labels and the ISI label "VIS"
    (visual cortex without an area) are tested but left out of it; `p_bh_all` is the correction
    over every pair, as in the first version of this notebook."""
    vals, ccg_sum, ccg_n = {}, {}, {}
    for sid, v in SP.items():
        if sid not in SETS[S]:
            continue
        e = v["alpha"].get((L, ref, E))
        if e is None or "coh" not in e or len(e["coh"]["event"]) <= 5:
            continue
        d = e["coh"]["event"] - e["coh"][R]
        for j, a in enumerate(e["coh_names"]):
            if a == ref:
                continue
            vals.setdefault(a, []).append(d[:, j])
            ccg_sum[a] = ccg_sum.get(a, 0) + e["ccg_mean"][a] * e["ccg_n"]
            ccg_n[a] = ccg_n.get(a, 0) + e["ccg_n"]
    rows = []
    lag_ms = (np.arange(200) - 100) * 10
    for a, lst in vals.items():
        x = np.concatenate(lst)
        x = x[np.isfinite(x)]
        if len(x) < min_events:
            continue
        c = ccg_sum[a] / ccg_n[a]
        rows.append({"area": a, "events": len(x), "sessions": len(lst), "mean_diff": x.mean(),
                     "p": stats.wilcoxon(x).pvalue, "ccg_peak": c.max(), "ccg_peak_lag_ms": lag_ms[np.argmax(c)]})
    t = pd.DataFrame(rows)
    if len(t):
        fam = ~t["area"].isin(NOT_FAMILY)
        t["p_bh_all"] = multipletests(t["p"], method="fdr_bh")[1]
        t["p_bh"] = np.nan
        t.loc[fam, "p_bh"] = multipletests(t.loc[fam, "p"], method="fdr_bh")[1]
        t.loc[~fam, "p_bh"] = t.loc[~fam, "p_bh_all"]
        t["sig_pos"] = (t["p_bh"] < 0.05) & (t["mean_diff"] > 0)
        t = t.sort_values("area").reset_index(drop=True)
    return t

PAIRS = {(ref, L, E, R, S): pair_table(ref, L, E, R, S) for ref in shield.REFERENCE_AREAS for (L, E, R, S) in COMBOS}
print("pairs tested per combination (VISp, MOp):",
      {f"{L}{E}{R}{S}": (len(PAIRS[('VISp', L, E, R, S)]), len(PAIRS[('MOp', L, E, R, S)])) for (L, E, R, S) in COMBOS})
changed = [(k, x.area) for k, t in PAIRS.items() if len(t) for x in t.itertuples()
           if x.area not in NOT_FAMILY and ((x.p_bh < 0.05) != (x.p_bh_all < 0.05))]
print("pairs whose significance at q < 0.05 changes with the family:", len(changed), "of",
      sum(int((~t["area"].isin(NOT_FAMILY)).sum()) for t in PAIRS.values() if len(t)), changed)

pairs tested per combination (VISp, MOp): {'L1E1R1a': (37, 12), 'L1E1R1b': (37, 12), 'L1E1R2a': (37, 12), 'L1E1R2b': (37, 12), 'L1E2R1a': (38, 22), 'L1E2R1b': (38, 22), 'L1E2R2a': (38, 22), 'L1E2R2b': (38, 22), 'L2E1R1a': (33, 11), 'L2E1R1b': (33, 11), 'L2E1R2a': (33, 11), 'L2E1R2b': (33, 11), 'L2E2R1a': (33, 21), 'L2E2R1b': (34, 21), 'L2E2R2a': (33, 21), 'L2E2R2b': (34, 21)}
pairs whose significance at q < 0.05 changes with the family: 4 of 734 [(('VISp', 'L2', 'E2', 'R2', 'a'), 'PO'), (('VISp', 'L2', 'E2', 'R2', 'b'), 'PO'), (('MOp', 'L1', 'E2', 'R2', 'a'), 'CA1'), (('MOp', 'L1', 'E2', 'R2', 'b'), 'CA1')]

rows = []
for L in ("L1", "L2"):
    for E in ("E1", "E2"):
        for S in ("a", "b"):
            ev_v = events_table("VISp", L, E, S); ev_m = events_table("MOp", L, E, S)
            ev_v, ev_m = ev_v[ev_v["qual"]], ev_m[ev_m["qual"]]
            for fcol, lab in (("freq_peaks", "a"), ("freq_auto", "b"), ("freq_auto_pad", "b, zero-padded")):
                xv, xm = ev_v[fcol].dropna(), ev_m[fcol].dropna()
                p = stats.ranksums(xm, xv, alternative="greater").pvalue if len(xm) and len(xv) else np.nan
                rows.append({"L": L, "E": E, "set": S, "reading": lab, "VISp_n": len(xv), "VISp_median_Hz": xv.median(),
                             "MOp_n": len(xm), "MOp_median_Hz": xm.median(), "p_MOp_higher": p,
                             "holds": (2.7 <= xv.median() <= 5.5) and (4.5 <= xm.median() <= 7.7) and p < 0.05})
fq = pd.DataFrame(rows)
print(fq.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
for lab, g in fq.groupby("reading"):
    print(f"reading ({lab}): holds in {int(g['holds'].sum())} of {len(g)}; MOp median {g['MOp_median_Hz'].min():.3g} to {g['MOp_median_Hz'].max():.3g} Hz, "
          f"VISp median {g['VISp_median_Hz'].min():.3g} to {g['VISp_median_Hz'].max():.3g} Hz, p {g['p_MOp_higher'].min():.2g} to {g['p_MOp_higher'].max():.2g}")
ev_m = events_table("MOp", "L2", "E2", "a"); ev_m = ev_m[ev_m["qual"]]
print("L2 E2 a, MOp events: share of freq_auto values at exactly 4.5 Hz", round(float((ev_m["freq_auto"] == 4.5).mean()), 3),
      "; values on the 0.5 Hz grid:", ev_m["freq_auto"].value_counts().sort_index().to_dict())
# the detector's band: synthetic sinusoidal rates of 2 to 8 Hz, 20 s each, from 0 to 20 spikes/s (the
# minima must reach the height rule, -1 on the inverted rate)
syn = []
for fz in np.arange(2, 8.01, 0.5):
    tt_ = np.arange(2000) * shield.BIN
    for E in ("E1", "E2"):
        syn.append({"Hz": fz, "rule": E, "events": len(shield.find_events((10 * (1 + np.sin(2 * np.pi * fz * tt_))).astype(np.float32), E)[0])})
syn = pd.DataFrame(syn).pivot(index="Hz", columns="rule", values="events")
print("events found in 20 s of a pure sinusoidal rate, by frequency:", syn.T.to_dict("index"))

 L  E set        reading  VISp_n  VISp_median_Hz  MOp_n  MOp_median_Hz  p_MOp_higher  holds
L1 E1   a              a    1518            4.17     80           4.55        0.0784  False
L1 E1   a              b    1518               4     80            4.5      0.000298   True
L1 E1   a b, zero-padded    1518            4.05     80           4.64      0.000668   True
L1 E1   b              a    1522            4.17     82           4.55        0.0972  False
L1 E1   b              b    1522               4     82            4.5      0.000235   True
L1 E1   b b, zero-padded    1522            4.05     82           4.64      0.000377   True
L1 E2   a              a    1694            4.17    119           4.35         0.568  False
L1 E2   a              b    1694               4    119            4.5      9.62e-07   True
L1 E2   a b, zero-padded    1694            4.05    119           4.74      1.69e-06   True
L1 E2   b              a    1704            4.17    122           4.35          0.67  False
L1 E2   b              b    1704               4    122            4.5      5.28e-07   True
L1 E2   b b, zero-padded    1704            4.05    122           4.74      7.06e-07   True
L2 E1   a              a    1469            4.17     80           4.55         0.116  False
L2 E1   a              b    1469               4     80            4.5       0.00158   True
L2 E1   a b, zero-padded    1469            4.05     80           4.64       0.00385   True
L2 E1   b              a    1478            4.17     82           4.55         0.138  False
L2 E1   b              b    1478               4     82            4.5        0.0012   True
L2 E1   b b, zero-padded    1478            4.05     82           4.64       0.00222   True
L2 E2   a              a    1643            4.17    119           4.35         0.684  False
L2 E2   a              b    1643               4    119            4.5      1.44e-05   True
L2 E2   a b, zero-padded    1643             4.1    119           4.74      4.13e-05   True
L2 E2   b              a    1660            4.21    122           4.35         0.765  False
L2 E2   b              b    1660               4    122            4.5      7.27e-06   True
L2 E2   b b, zero-padded    1660             4.1    122           4.74      1.67e-05   True
reading (a): holds in 0 of 8; MOp median 4.35 to 4.55 Hz, VISp median 4.17 to 4.21 Hz, p 0.078 to 0.76
reading (b): holds in 8 of 8; MOp median 4.5 to 4.5 Hz, VISp median 4 to 4 Hz, p 5.3e-07 to 0.0016
reading (b, zero-padded): holds in 8 of 8; MOp median 4.64 to 4.74 Hz, VISp median 4.05 to 4.1 Hz, p 7.1e-07 to 0.0038
L2 E2 a, MOp events: share of freq_auto values at exactly 4.5 Hz 0.227 ; values on the 0.5 Hz grid: {2.0: 2, 2.5: 1, 3.0: 3, 3.5: 8, 4.0: 21, 4.5: 27, 5.0: 27, 5.5: 14, 6.0: 9, 6.5: 5, 7.0: 2}

events found in 20 s of a pure sinusoidal rate, by frequency: {'E1': {2.0: 0, 2.5: 0, 3.0: 7, 3.5: 8, 4.0: 8, 4.5: 9, 5.0: 9, 5.5: 0, 6.0: 0, 6.5: 0, 7.0: 0, 7.5: 0, 8.0: 0}, 'E2': {2.0: 0, 2.5: 0, 3.0: 7, 3.5: 8, 4.0: 8, 4.5: 9, 5.0: 9, 5.5: 0, 6.0: 0, 6.5: 0, 7.0: 0, 7.5: 0, 8.0: 0}}

rows = []
for L in ("L1", "L2"):
    for E in ("E1", "E2"):
        for S in ("a", "b"):
            per = {}
            for ref in shield.REFERENCE_AREAS:
                ev = events_table(ref, L, E, S)
                ev = ev[ev["qual"]].copy()
                ref_mean = ev.groupby("session_id")["pupil_win_mean_rand"].transform("mean")
                ev["chg_area"] = (ev["pupil"] - ref_mean) / ref_mean
                ev["chg_diam"] = (np.sqrt(ev["pupil"]) - np.sqrt(ref_mean)) / np.sqrt(ref_mean)
                per[ref] = ev
                run_e, run_r = ev["running"].dropna(), ev["running_rand"].dropna()
                for meas in ("area", "diam"):
                    c = ev["chg_" + meas].dropna()
                    rows.append({"L": L, "E": E, "set": S, "events": ref, "measure": meas, "n": len(c),
                                 "mean_pct": 100 * c.mean(), "median_pct": 100 * c.median(),
                                 "p_pupil": stats.wilcoxon(c).pvalue if len(c) > 5 else np.nan,
                                 "running_event": run_e.mean(), "running_random": run_r.mean(),
                                 "p_running_lower": stats.ranksums(run_e, run_r, alternative="less").pvalue if len(run_e) > 5 else np.nan})
            both = pd.concat(per.values())
            for meas in ("area", "diam"):
                c = both["chg_" + meas].dropna()
                rows.append({"L": L, "E": E, "set": S, "events": "pooled", "measure": meas, "n": len(c),
                             "mean_pct": 100 * c.mean(), "median_pct": 100 * c.median(),
                             "p_pupil": stats.wilcoxon(c).pvalue if len(c) > 5 else np.nan})
ar = pd.DataFrame(rows)
ar["within_tol_mean"] = (ar["mean_pct"] >= -27.5) & (ar["mean_pct"] <= -22.5)
ar["within_tol_median"] = (ar["median_pct"] >= -27.5) & (ar["median_pct"] <= -22.5)
print(ar.to_string(index=False, float_format=lambda v: f"{v:.3g}"))

 L  E set events measure    n  mean_pct  median_pct   p_pupil  running_event  running_random  p_running_lower  within_tol_mean  within_tol_median
L1 E1   a   VISp    area 1516     -12.9       -14.8  6.67e-97          0.235            4.19         8.45e-76            False              False
L1 E1   a   VISp    diam 1516     -7.48       -7.69  2.2e-104          0.235            4.19         8.45e-76            False              False
L1 E1   a    MOp    area   80     -36.5       -38.9  8.15e-15         0.0323            7.11         7.77e-10            False              False
L1 E1   a    MOp    diam   80     -20.8       -21.9  8.15e-15         0.0323            7.11         7.77e-10            False              False
L1 E1   a pooled    area 1596     -14.1       -16.1 9.62e-112            NaN             NaN              NaN            False              False
L1 E1   a pooled    diam 1596     -8.15       -8.38 1.16e-119            NaN             NaN              NaN            False              False
L1 E1   b   VISp    area 1520     -12.9       -14.8   3.9e-97          0.235            4.27         1.14e-76            False              False
L1 E1   b   VISp    diam 1520     -7.46       -7.68 1.33e-104          0.235            4.27         1.14e-76            False              False
L1 E1   b    MOp    area   82     -35.7       -37.9   3.8e-15           1.16            7.66         2.67e-09            False              False
L1 E1   b    MOp    diam   82     -20.4       -21.2   3.8e-15           1.16            7.66         2.67e-09            False              False
L1 E1   b pooled    area 1602       -14       -15.9 4.12e-112            NaN             NaN              NaN            False              False
L1 E1   b pooled    diam 1602     -8.12       -8.28 5.28e-120            NaN             NaN              NaN            False              False
L1 E2   a   VISp    area 1692     -13.2       -15.1 8.23e-115          0.271            4.44         4.61e-79            False              False
L1 E2   a   VISp    diam 1692      -7.6       -7.85 4.29e-123          0.271            4.44         4.61e-79            False              False
L1 E2   a    MOp    area  119     -35.6       -36.9  3.19e-21           0.04            5.67         5.35e-11            False              False
L1 E2   a    MOp    diam  119     -20.5       -20.6  3.19e-21           0.04            5.67         5.35e-11            False              False
L1 E2   a pooled    area 1811     -14.7       -16.3 1.83e-136            NaN             NaN              NaN            False              False
L1 E2   a pooled    diam 1811     -8.44       -8.51  2.6e-145            NaN             NaN              NaN            False              False
L1 E2   b   VISp    area 1702     -13.1       -14.9 3.18e-115          0.281            4.65         5.46e-81            False              False
L1 E2   b   VISp    diam 1702     -7.56       -7.78 1.69e-123          0.281            4.65         5.46e-81            False              False
L1 E2   b    MOp    area  122     -34.7       -35.8  1.21e-21           1.11            6.58          2.3e-10            False              False
L1 E2   b    MOp    diam  122       -20       -19.9  1.21e-21           1.11            6.58          2.3e-10            False              False
L1 E2   b pooled    area 1824     -14.6       -16.2 5.84e-137            NaN             NaN              NaN            False              False
L1 E2   b pooled    diam 1824     -8.39       -8.45 8.84e-146            NaN             NaN              NaN            False              False
L2 E1   a   VISp    area 1467     -15.4       -17.5 9.26e-126          0.212            4.15         2.15e-80            False              False
L2 E1   a   VISp    diam 1467     -8.81       -9.19 1.64e-132          0.212            4.15         2.15e-80            False              False
L2 E1   a    MOp    area   80     -31.7       -31.6  9.13e-15         0.0323            5.16         3.96e-06            False              False
L2 E1   a    MOp    diam   80     -17.8       -17.3  9.13e-15         0.0323            5.16         3.96e-06            False              False
L2 E1   a pooled    area 1547     -16.3       -18.5 5.97e-140            NaN             NaN              NaN            False              False
L2 E1   a pooled    diam 1547     -9.27       -9.73 5.06e-147            NaN             NaN              NaN            False              False
L2 E1   b   VISp    area 1476     -15.4       -17.5 1.52e-126          0.248            4.38         6.48e-81            False              False
L2 E1   b   VISp    diam 1476     -8.78       -9.17 2.71e-133          0.248            4.38         6.48e-81            False              False
L2 E1   b    MOp    area   82     -30.8       -31.2  6.84e-15           1.16            6.23         7.98e-06            False              False
L2 E1   b    MOp    diam   82     -17.3         -17   6.6e-15           1.16            6.23         7.98e-06            False              False
L2 E1   b pooled    area 1558     -16.2       -18.3 1.12e-140            NaN             NaN              NaN            False              False
L2 E1   b pooled    diam 1558     -9.23        -9.6 9.87e-148            NaN             NaN              NaN            False              False
L2 E2   a   VISp    area 1641     -14.5       -16.1 7.89e-134          0.297            4.18         1.16e-82            False              False
L2 E2   a   VISp    diam 1641     -8.26       -8.39 4.75e-141          0.297            4.18         1.16e-82            False              False
L2 E2   a    MOp    area  119     -36.4       -38.6  3.19e-21           0.04            5.25         3.32e-12            False              False
L2 E2   a    MOp    diam  119     -20.9       -21.7  3.19e-21           0.04            5.25         3.32e-12            False              False
L2 E2   a pooled    area 1760       -16       -17.3 3.38e-156            NaN             NaN              NaN            False              False
L2 E2   a pooled    diam 1760     -9.11       -9.07 8.01e-164            NaN             NaN              NaN            False              False
L2 E2   b   VISp    area 1658     -14.4       -15.9 7.85e-135          0.359            4.55         4.99e-84            False              False
L2 E2   b   VISp    diam 1658      -8.2       -8.29 5.01e-142          0.359            4.55         4.99e-84            False              False
L2 E2   b    MOp    area  122     -35.5       -38.2  1.44e-21           1.11            5.83         1.73e-11            False              False
L2 E2   b    MOp    diam  122     -20.4       -21.4  1.44e-21           1.11            5.83         1.73e-11            False              False
L2 E2   b pooled    area 1780     -15.8       -17.1 3.07e-157            NaN             NaN              NaN            False              False
L2 E2   b pooled    diam 1780     -9.03       -8.96 8.05e-165            NaN             NaN              NaN            False              False

rows = []
for L in ("L1", "L2"):
    for E in ("E1", "E2"):
        for S in ("a", "b"):
            for ref in shield.REFERENCE_AREAS:
                ev = events_table(ref, L, E, S)
                x = ev.loc[ev["qual"], "since_offset"].dropna()
                x = x[(x >= 0) & (x < 0.5)].values
                h = np.histogram(x, np.arange(0, 0.5001, 0.05))[0]
                rows.append({"L": L, "E": E, "set": S, "area": ref, "n_in_0_500ms": len(x), "bins": h.tolist(),
                             "largest_bin_ms": int(np.argmax(h) * 50) if len(x) else -1,
                             "p_chi2": stats.chisquare(h).pvalue if len(x) else np.nan,
                             "p_ks": stats.kstest(x, stats.uniform(0, 0.5).cdf).pvalue if len(x) else np.nan})
off = pd.DataFrame(rows)
print(off.to_string(index=False, float_format=lambda v: f"{v:.3g}"))

 L  E set area  n_in_0_500ms                                       bins  largest_bin_ms  p_chi2      p_ks
L1 E1   a VISp          1034   [8, 41, 634, 243, 23, 6, 22, 22, 15, 20]             100       0 1.44e-248
L1 E1   a  MOp            53            [10, 5, 5, 6, 2, 3, 7, 6, 4, 5]               0   0.502     0.375
L1 E1   b VISp          1035   [8, 41, 635, 243, 23, 6, 22, 22, 15, 20]             100       0 6.09e-249
L1 E1   b  MOp            55            [10, 6, 5, 6, 2, 3, 7, 6, 4, 6]               0   0.525     0.283
L1 E2   a VISp          1079   [9, 45, 644, 244, 24, 8, 25, 27, 23, 30]             100       0 9.15e-233
L1 E2   a  MOp            73           [11, 7, 6, 9, 3, 6, 9, 8, 4, 10]               0   0.511     0.623
L1 E2   b VISp          1084   [9, 45, 649, 244, 24, 8, 25, 27, 23, 30]             100       0  1.4e-234
L1 E2   b  MOp            76          [11, 8, 6, 10, 3, 6, 9, 8, 4, 11]               0   0.413     0.549
L2 E1   a VISp          1054  [8, 48, 627, 261, 14, 12, 24, 30, 14, 16]             100       0  8.1e-248
L2 E1   a  MOp            53            [10, 5, 5, 6, 2, 3, 7, 6, 4, 5]               0   0.502     0.375
L2 E1   b VISp          1059  [8, 48, 629, 264, 14, 12, 24, 30, 14, 16]             100       0  1.3e-249
L2 E1   b  MOp            55            [10, 6, 5, 6, 2, 3, 7, 6, 4, 6]               0   0.525     0.283
L2 E2   a VISp          1110 [10, 52, 651, 261, 15, 11, 27, 37, 21, 25]             100       0 1.71e-239
L2 E2   a  MOp            73           [11, 7, 6, 9, 3, 6, 9, 8, 4, 10]               0   0.511     0.623
L2 E2   b VISp          1119 [10, 52, 657, 264, 15, 11, 27, 37, 21, 25]             100       0 1.16e-242
L2 E2   b  MOp            76          [11, 8, 6, 10, 3, 6, 9, 8, 4, 11]               0   0.413     0.549

rows = []
for L in ("L1", "L2"):
    for E in ("E1", "E2"):
        for S in ("a", "b"):
            for ew in (0.25, 0.5):
                for vw, post in ((1.0, False), (1.25, False), (1.0, True), (1.25, True)):
                    xs, ys, ux, uy, nsess, nev = [], [], [], [], 0, 0
                    for sid, v in SP.items():
                        if sid not in SETS[S]:
                            continue
                        e = v["alpha"].get((L, "VISp", E))
                        if e is None or "replay" not in e:
                            continue
                        rp = e["replay"]
                        prev = np.array(rp["event_prev_image"], dtype=object)
                        so = rp["event_since_offset"]
                        keep = (so >= 0) & (so < 0.5) if post else np.ones(len(prev), bool)
                        imgs = rp["images"]
                        if not all(((prev == im) & keep).any() for im in imgs):
                            continue
                        evk = rp["evoked"][ew]
                        evr = np.stack([rp["event_rates"][vw][:, (prev == im) & keep].mean(1) for im in imgs], 1)
                        with np.errstate(all="ignore"):
                            a = evk / evk.max(1, keepdims=True)
                            b = evr / evr.max(1, keepdims=True)
                        ok = np.isfinite(a).all(1) & np.isfinite(b).all(1)
                        if ok.sum() == 0:
                            continue
                        xs += list(a[ok].mean(0)); ys += list(b[ok].mean(0)); ux += list(a[ok].ravel()); uy += list(b[ok].ravel())
                        nsess += 1; nev += int(sum(((prev == im) & keep).sum() for im in imgs))
                    r, p = stats.pearsonr(xs, ys) if len(xs) > 3 else (np.nan, np.nan)
                    ru, pu = stats.pearsonr(ux, uy) if len(ux) > 3 else (np.nan, np.nan)
                    rows.append({"L": L, "E": E, "set": S, "evoked_s": ew, "event_s": vw, "post_offset_only": post, "sessions": nsess, "events": nev, "points": len(xs),
                                 "r": r, "p": p, "per_unit_r": ru, "per_unit_points": len(ux), "per_unit_p": pu})
rep = pd.DataFrame(rows)
print(rep.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
rep["holds"] = (rep["r"] >= 0.693) & (rep["r"] <= 0.847) & (rep["p"] < 0.05)
for post, g in rep.groupby("post_offset_only"):
    print(f"{'events 0 to 0.5 s after the offset only' if post else 'all events'}: r {g['r'].min():.3f} to {g['r'].max():.3f}; "
          f"p < 0.05 in {int((g['p'] < 0.05).sum())} of {len(g)}; within 0.77 +- 0.077 with p < 0.05 in {int(g['holds'].sum())}")
so_all = np.concatenate([v["alpha"][("L2", "VISp", "E2")]["replay"]["event_since_offset"] for v in SP.values()
                         if "replay" in v["alpha"].get(("L2", "VISp", "E2"), {})])
print("L2 E2 qualifying VISp events by time since the last active presentation's offset (s):",
      dict(zip(["< -0.25", "-0.25 to 0 (during the image)", "0 to 0.5", "0.5 to 5", "over 5 (after the active block)"],
               np.histogram(so_all, [-np.inf, -0.25, 0, 0.5, 5, np.inf])[0].tolist())))

 L  E set  evoked_s  event_s  post_offset_only  sessions  events  points     r        p  per_unit_r  per_unit_points  per_unit_p
L1 E1   a      0.25        1             False         7    1429      28  0.39   0.0402       0.546             1464   2.05e-114
L1 E1   a      0.25     1.25             False         7    1429      28 0.633 0.000302         0.6             1464   4.22e-144
L1 E1   a      0.25        1              True         7     990      28 0.335   0.0812       0.547             1452   3.18e-114
L1 E1   a      0.25     1.25              True         7     990      28 0.527  0.00396       0.608             1464   1.56e-148
L1 E1   a       0.5        1             False         7    1429      28 0.366   0.0552       0.627             1464   8.74e-161
L1 E1   a       0.5     1.25             False         7    1429      28 0.599 0.000762       0.674             1464   1.92e-194
L1 E1   a       0.5        1              True         7     990      28 0.326   0.0905       0.611             1452   4.91e-149
L1 E1   a       0.5     1.25              True         7     990      28 0.498  0.00695       0.655             1464   6.46e-180
L1 E1   b      0.25        1             False         7    1429      28  0.39   0.0402       0.546             1464   2.05e-114
L1 E1   b      0.25     1.25             False         7    1429      28 0.633 0.000302         0.6             1464   4.22e-144
L1 E1   b      0.25        1              True         7     990      28 0.335   0.0812       0.547             1452   3.18e-114
L1 E1   b      0.25     1.25              True         7     990      28 0.527  0.00396       0.608             1464   1.56e-148
L1 E1   b       0.5        1             False         7    1429      28 0.366   0.0552       0.627             1464   8.74e-161
L1 E1   b       0.5     1.25             False         7    1429      28 0.599 0.000762       0.674             1464   1.92e-194
L1 E1   b       0.5        1              True         7     990      28 0.326   0.0905       0.611             1452   4.91e-149
L1 E1   b       0.5     1.25              True         7     990      28 0.498  0.00695       0.655             1464   6.46e-180
L1 E2   a      0.25        1             False         7    1578      28 0.274    0.159       0.543             1464   3.81e-113
L1 E2   a      0.25     1.25             False         7    1578      28 0.526  0.00403        0.59             1464   2.76e-138
L1 E2   a      0.25        1              True         7    1027      28 0.254    0.193       0.537             1460    5.6e-110
L1 E2   a      0.25     1.25              True         7    1027      28 0.449   0.0165       0.595             1464   5.28e-141
L1 E2   a       0.5        1             False         7    1578      28 0.275    0.156       0.624             1464    1.6e-158
L1 E2   a       0.5     1.25             False         7    1578      28  0.51  0.00555       0.667             1464   7.41e-189
L1 E2   a       0.5        1              True         7    1027      28 0.274    0.159         0.6             1460   1.32e-143
L1 E2   a       0.5     1.25              True         7    1027      28 0.442   0.0185       0.643             1464   1.46e-171
L1 E2   b      0.25        1             False         7    1578      28 0.274    0.159       0.543             1464   3.81e-113
L1 E2   b      0.25     1.25             False         7    1578      28 0.526  0.00403        0.59             1464   2.76e-138
L1 E2   b      0.25        1              True         7    1027      28 0.254    0.193       0.537             1460    5.6e-110
L1 E2   b      0.25     1.25              True         7    1027      28 0.449   0.0165       0.595             1464   5.28e-141
L1 E2   b       0.5        1             False         7    1578      28 0.275    0.156       0.624             1464    1.6e-158
L1 E2   b       0.5     1.25             False         7    1578      28  0.51  0.00555       0.667             1464   7.41e-189
L1 E2   b       0.5        1              True         7    1027      28 0.274    0.159         0.6             1460   1.32e-143
L1 E2   b       0.5     1.25              True         7    1027      28 0.442   0.0185       0.643             1464   1.46e-171
L2 E1   a      0.25        1             False         8    1396      32 0.308   0.0861       0.463             1960   1.72e-104
L2 E1   a      0.25     1.25             False         8    1396      32 0.473  0.00626       0.508             1964   2.65e-129
L2 E1   a      0.25        1              True         7    1023      28 0.292    0.131       0.533             1716   2.75e-126
L2 E1   a      0.25     1.25              True         7    1023      28 0.461   0.0137       0.588             1720   2.26e-160
L2 E1   a       0.5        1             False         8    1396      32 0.275    0.128       0.537             1960   1.29e-146
L2 E1   a       0.5     1.25             False         8    1396      32 0.431   0.0139       0.578             1964   8.98e-176
L2 E1   a       0.5        1              True         7    1023      28 0.302    0.118       0.593             1716   3.44e-163
L2 E1   a       0.5     1.25              True         7    1023      28 0.465   0.0126       0.634             1720   7.08e-194
L2 E1   b      0.25        1             False         8    1396      32 0.308   0.0861       0.463             1960   1.72e-104
L2 E1   b      0.25     1.25             False         8    1396      32 0.473  0.00626       0.508             1964   2.65e-129
L2 E1   b      0.25        1              True         7    1023      28 0.292    0.131       0.533             1716   2.75e-126
L2 E1   b      0.25     1.25              True         7    1023      28 0.461   0.0137       0.588             1720   2.26e-160
L2 E1   b       0.5        1             False         8    1396      32 0.275    0.128       0.537             1960   1.29e-146
L2 E1   b       0.5     1.25             False         8    1396      32 0.431   0.0139       0.578             1964   8.98e-176
L2 E1   b       0.5        1              True         7    1023      28 0.302    0.118       0.593             1716   3.44e-163
L2 E1   b       0.5     1.25              True         7    1023      28 0.465   0.0126       0.634             1720   7.08e-194
L2 E2   a      0.25        1             False         8    1549      32 0.233      0.2       0.493             1964   7.49e-121
L2 E2   a      0.25     1.25             False         8    1549      32 0.389   0.0277       0.532             1964   4.82e-144
L2 E2   a      0.25        1              True         7    1077      28 0.193    0.326       0.524             1716   1.26e-121
L2 E2   a      0.25     1.25              True         7    1077      28 0.368   0.0539       0.582             1720   1.54e-156
L2 E2   a       0.5        1             False         8    1549      32 0.208    0.253       0.566             1964   8.42e-167
L2 E2   a       0.5     1.25             False         8    1549      32 0.358   0.0442       0.604             1964   1.78e-195
L2 E2   a       0.5        1              True         7    1077      28 0.215    0.272       0.582             1716   7.89e-156
L2 E2   a       0.5     1.25              True         7    1077      28 0.384   0.0434       0.625             1720   9.79e-187
L2 E2   b      0.25        1             False         8    1549      32 0.233      0.2       0.493             1964   7.49e-121
L2 E2   b      0.25     1.25             False         8    1549      32 0.389   0.0277       0.532             1964   4.82e-144
L2 E2   b      0.25        1              True         7    1077      28 0.193    0.326       0.524             1716   1.26e-121
L2 E2   b      0.25     1.25              True         7    1077      28 0.368   0.0539       0.582             1720   1.54e-156
L2 E2   b       0.5        1             False         8    1549      32 0.208    0.253       0.566             1964   8.42e-167
L2 E2   b       0.5     1.25             False         8    1549      32 0.358   0.0442       0.604             1964   1.78e-195
L2 E2   b       0.5        1              True         7    1077      28 0.215    0.272       0.582             1716   7.89e-156
L2 E2   b       0.5     1.25              True         7    1077      28 0.384   0.0434       0.625             1720   9.79e-187

all events: r 0.208 to 0.633; p < 0.05 in 18 of 32; within 0.77 +- 0.077 with p < 0.05 in 0
events 0 to 0.5 s after the offset only: r 0.193 to 0.527; p < 0.05 in 14 of 32; within 0.77 +- 0.077 with p < 0.05 in 0
L2 E2 qualifying VISp events by time since the last active presentation's offset (s): {'< -0.25': 0, '-0.25 to 0 (during the image)': 473, '0 to 0.5': 1119, '0.5 to 5': 0, 'over 5 (after the active block)': 64}

def lookup(t, names):
    """Rows of a pair table for an area given as one or more labels (subdivisions)."""
    return t[t["area"].isin(names)] if len(t) else t

def area_state(t, names, mode="any"):
    r = lookup(t, names)
    if not len(r):
        return "absent"
    s = r["sig_pos"].all() if mode == "all" else r["sig_pos"].any()
    return "sig" if s else "ns"

def show_pairs(ref, areas):
    out = []
    for (L, E, R, S) in COMBOS:
        t = PAIRS[(ref, L, E, R, S)]
        row = {"combo": f"{L} {E} {R} {S}"}
        for a, names in areas.items():
            r = lookup(t, names)
            row[a] = ("; ".join(f"{x.area} {x.mean_diff:+.3f} q={x.p_bh:.2g}{'*' if x.sig_pos else ''}" for x in r.itertuples())
                      if len(r) else "absent")
        out.append(row)
    return pd.DataFrame(out)

A = {"LP": ["LP"], "VISl": ["VISl", "VISlm"], "PO": ["PO"], "SSp": ["SSp"]}
print(show_pairs("VISp", A).to_string(index=False))
rows = []
for (L, E, R, S) in COMBOS:
    t = PAIRS[("VISp", L, E, R, S)]
    st = {a: area_state(t, n) for a, n in A.items()}
    coh_ok = st["LP"] == "sig" and st["VISl"] == "sig" and st["PO"] != "sig" and st["SSp"] != "sig"
    untested = ", ".join(a for a in A if st[a] == "absent") or "none"
    c = {a: lookup(t, n) for a, n in A.items()}
    pk = {a: (c[a]["ccg_peak"].max() if len(c[a]) else np.nan) for a in A}
    lg = {a: (c[a].sort_values("ccg_peak")["ccg_peak_lag_ms"].iloc[-1] if len(c[a]) else np.nan) for a in A}
    for w in (20, 50):
        near = all(abs(lg[a]) <= w for a in ("LP", "VISl") if np.isfinite(lg[a]))
        higher = min(pk["LP"], pk["VISl"]) > np.nanmax([pk["PO"], pk["SSp"]]) if np.isfinite([pk["PO"], pk["SSp"]]).any() else np.nan
        rows.append({"combo": f"{L} {E} {R} {S}", "ccg_window_ms": w, "coherence_part": coh_ok, "untested_parts": untested,
                     **{f"peak_{a}": round(pk[a], 3) for a in A}, **{f"lag_{a}": lg[a] for a in A},
                     "ccg_near_zero": near, "ccg_higher": higher})
v7 = pd.DataFrame(rows)
print(v7.to_string(index=False))
print("combinations where the coherence part holds:", v7.loc[v7["coherence_part"], "combo"].unique().tolist())

     combo                    LP                                            VISl                PO                   SSp
L1 E1 R1 a    LP +0.068 q=1e-20* VISl +0.158 q=3.6e-10*; VISlm +0.130 q=1.8e-08*  PO +0.005 q=0.33    SSp +0.005 q=0.097
L1 E1 R1 b    LP +0.068 q=1e-20* VISl +0.158 q=3.6e-10*; VISlm +0.130 q=1.8e-08*  PO +0.005 q=0.33    SSp +0.005 q=0.097
L1 E1 R2 a LP +0.266 q=3.6e-122*   VISl +0.375 q=4e-17*; VISlm +0.306 q=9.4e-17* PO +0.007 q=0.098 SSp +0.011 q=2.5e-05*
L1 E1 R2 b LP +0.266 q=3.6e-122*   VISl +0.375 q=4e-17*; VISlm +0.306 q=9.4e-17* PO +0.007 q=0.098 SSp +0.011 q=2.5e-05*
L1 E2 R1 a  LP +0.064 q=5.8e-19* VISl +0.151 q=4.5e-11*; VISlm +0.123 q=1.6e-08*  PO +0.005 q=0.29     SSp +0.004 q=0.23
L1 E2 R1 b  LP +0.064 q=5.8e-19* VISl +0.151 q=4.5e-11*; VISlm +0.123 q=1.6e-08*  PO +0.005 q=0.29     SSp +0.004 q=0.23
L1 E2 R2 a LP +0.264 q=2.1e-134* VISl +0.346 q=9.6e-19*; VISlm +0.294 q=1.1e-20*  PO +0.000 q=0.96  SSp +0.008 q=0.0032*
L1 E2 R2 b LP +0.264 q=2.1e-134* VISl +0.346 q=9.6e-19*; VISlm +0.294 q=1.1e-20*  PO +0.000 q=0.96  SSp +0.008 q=0.0032*
L2 E1 R1 a  LP +0.091 q=3.6e-27*                            VISl +0.124 q=2e-10*  PO -0.004 q=0.68     SSp +0.003 q=0.68
L2 E1 R1 b  LP +0.091 q=3.6e-27*                            VISl +0.124 q=2e-10*  PO -0.004 q=0.68     SSp +0.003 q=0.68
L2 E1 R2 a LP +0.305 q=1.1e-110*                          VISl +0.343 q=2.1e-24*  PO +0.005 q=0.27    SSp +0.006 q=0.075
L2 E1 R2 b LP +0.305 q=1.1e-110*                          VISl +0.343 q=2.1e-24*  PO +0.005 q=0.27    SSp +0.006 q=0.075
L2 E2 R1 a  LP +0.086 q=5.2e-29*                          VISl +0.137 q=3.7e-13*  PO +0.003 q=0.47      SSp +0.004 q=0.3
L2 E2 R1 b  LP +0.086 q=5.4e-29*                          VISl +0.137 q=3.9e-13*  PO +0.003 q=0.49     SSp +0.004 q=0.31
L2 E2 R2 a LP +0.302 q=5.5e-121*                          VISl +0.347 q=6.5e-29* PO +0.011 q=0.051   SSp +0.007 q=0.048*
L2 E2 R2 b LP +0.302 q=5.6e-121*                          VISl +0.347 q=6.7e-29* PO +0.011 q=0.051   SSp +0.007 q=0.036*

     combo  ccg_window_ms  coherence_part untested_parts  peak_LP  peak_VISl  peak_PO  peak_SSp  lag_LP  lag_VISl  lag_PO  lag_SSp  ccg_near_zero  ccg_higher
L1 E1 R1 a             20            True           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R1 a             50            True           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R1 b             20            True           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R1 b             50            True           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R2 a             20           False           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R2 a             50           False           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R2 b             20           False           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E1 R2 b             50           False           none    0.340      0.484    0.058     0.055       0         0      30        0           True        True
L1 E2 R1 a             20            True           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R1 a             50            True           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R1 b             20            True           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R1 b             50            True           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R2 a             20           False           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R2 a             50           False           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R2 b             20           False           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L1 E2 R2 b             50           False           none    0.335      0.477    0.056     0.055       0         0      20       20           True        True
L2 E1 R1 a             20            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R1 a             50            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R1 b             20            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R1 b             50            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R2 a             20            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R2 a             50            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R2 b             20            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E1 R2 b             50            True           none    0.384      0.508    0.092     0.062       0         0      20        0           True        True
L2 E2 R1 a             20            True           none    0.378      0.498    0.092     0.061       0         0      20        0           True        True
L2 E2 R1 a             50            True           none    0.378      0.498    0.092     0.061       0         0      20        0           True        True
L2 E2 R1 b             20            True           none    0.378      0.498    0.092     0.060       0         0      20        0           True        True
L2 E2 R1 b             50            True           none    0.378      0.498    0.092     0.060       0         0      20        0           True        True
L2 E2 R2 a             20           False           none    0.378      0.498    0.092     0.061       0         0      20        0           True        True
L2 E2 R2 a             50           False           none    0.378      0.498    0.092     0.061       0         0      20        0           True        True
L2 E2 R2 b             20           False           none    0.378      0.498    0.092     0.060       0         0      20        0           True        True
L2 E2 R2 b             50           False           none    0.378      0.498    0.092     0.060       0         0      20        0           True        True
combinations where the coherence part holds: ['L1 E1 R1 a', 'L1 E1 R1 b', 'L1 E2 R1 a', 'L1 E2 R1 b', 'L2 E1 R1 a', 'L2 E1 R1 b', 'L2 E1 R2 a', 'L2 E1 R2 b', 'L2 E2 R1 a', 'L2 E2 R1 b']

A2 = {"PO": ["PO"], "SSp": ["SSp"], "LP": ["LP"], "VISl": ["VISl", "VISlm"]}
print(show_pairs("MOp", A2).to_string(index=False))
rows = []
for (L, E, R, S) in COMBOS:
    t = PAIRS[("MOp", L, E, R, S)]
    st = {a: area_state(t, n) for a, n in A2.items()}
    c = {a: lookup(t, n) for a, n in A2.items()}
    pk = {a: (c[a]["ccg_peak"].max() if len(c[a]) else np.nan) for a in A2}
    lg = {a: (c[a].sort_values("ccg_peak")["ccg_peak_lag_ms"].iloc[-1] if len(c[a]) else np.nan) for a in A2}
    rows.append({"combo": f"{L} {E} {R} {S}", **st,
                 "coherence_part": st["PO"] == "sig" and st["SSp"] == "sig" and st["LP"] != "sig" and st["VISl"] != "sig",
                 "untested_parts": ", ".join(a for a in A2 if st[a] == "absent") or "none",
                 "contradicted": st["PO"] == "ns" or st["SSp"] == "ns" or st["LP"] == "sig" or st["VISl"] == "sig",
                 **{f"peak_{a}": round(pk[a], 3) for a in A2}, **{f"lag_{a}": lg[a] for a in A2}})
print(pd.DataFrame(rows).to_string(index=False))
m7 = pd.DataFrame(rows)
print("combinations with every part tested and holding:", m7.loc[m7["coherence_part"] & (m7["untested_parts"] == "none"), "combo"].tolist())
print("combinations where a tested part contradicts the claim:", m7.loc[m7["contradicted"], "combo"].tolist())

     combo                   PO                   SSp               LP               VISl
L1 E1 R1 a               absent SSp +0.181 q=8.8e-07* LP +0.000 q=0.79             absent
L1 E1 R1 b               absent SSp +0.181 q=8.8e-07* LP +0.000 q=0.79             absent
L1 E1 R2 a               absent   SSp +0.259 q=8e-14* LP +0.025 q=0.14             absent
L1 E1 R2 b               absent   SSp +0.259 q=8e-14* LP +0.025 q=0.14             absent
L1 E2 R1 a PO +0.174 q=1.1e-05* SSp +0.190 q=4.3e-09* LP -0.002 q=0.86             absent
L1 E2 R1 b PO +0.174 q=1.1e-05* SSp +0.190 q=4.3e-09* LP -0.002 q=0.86             absent
L1 E2 R2 a PO +0.195 q=5.3e-06*   SSp +0.250 q=1e-11* LP -0.006 q=0.77             absent
L1 E2 R2 b PO +0.195 q=5.3e-06*   SSp +0.250 q=1e-11* LP -0.006 q=0.77             absent
L2 E1 R1 a               absent SSp +0.173 q=2.7e-10*  LP +0.026 q=0.2             absent
L2 E1 R1 b               absent SSp +0.173 q=2.7e-10*  LP +0.026 q=0.2             absent
L2 E1 R2 a               absent SSp +0.246 q=8.6e-12* LP +0.018 q=0.49             absent
L2 E1 R2 b               absent SSp +0.246 q=8.6e-12* LP +0.018 q=0.49             absent
L2 E2 R1 a  PO +0.147 q=0.0015* SSp +0.187 q=5.4e-11* LP +0.001 q=0.89 VISl +0.005 q=0.89
L2 E2 R1 b  PO +0.147 q=0.0015* SSp +0.187 q=5.4e-11* LP +0.001 q=0.89 VISl +0.005 q=0.89
L2 E2 R2 a PO +0.169 q=1.6e-05* SSp +0.247 q=2.2e-11* LP -0.009 q=0.56 VISl +0.027 q=0.34
L2 E2 R2 b PO +0.169 q=1.6e-05* SSp +0.247 q=2.2e-11* LP -0.009 q=0.56 VISl +0.027 q=0.34
     combo     PO SSp LP   VISl  coherence_part untested_parts  contradicted  peak_PO  peak_SSp  peak_LP  peak_VISl  lag_PO  lag_SSp  lag_LP  lag_VISl
L1 E1 R1 a absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.035        NaN     NaN        0     180       NaN
L1 E1 R1 b absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.035        NaN     NaN        0     180       NaN
L1 E1 R2 a absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.035        NaN     NaN        0     180       NaN
L1 E1 R2 b absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.035        NaN     NaN        0     180       NaN
L1 E2 R1 a    sig sig ns absent            True           VISl         False    0.250     0.329    0.029        NaN    20.0        0     210       NaN
L1 E2 R1 b    sig sig ns absent            True           VISl         False    0.250     0.329    0.029        NaN    20.0        0     210       NaN
L1 E2 R2 a    sig sig ns absent            True           VISl         False    0.250     0.329    0.029        NaN    20.0        0     210       NaN
L1 E2 R2 b    sig sig ns absent            True           VISl         False    0.250     0.329    0.029        NaN    20.0        0     210       NaN
L2 E1 R1 a absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.033        NaN     NaN        0    -400       NaN
L2 E1 R1 b absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.033        NaN     NaN        0    -400       NaN
L2 E1 R2 a absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.033        NaN     NaN        0    -400       NaN
L2 E1 R2 b absent sig ns absent           False       PO, VISl         False      NaN     0.325    0.033        NaN     NaN        0    -400       NaN
L2 E2 R1 a    sig sig ns     ns            True           none         False    0.199     0.329    0.037      0.062    20.0        0    -400     320.0
L2 E2 R1 b    sig sig ns     ns            True           none         False    0.199     0.329    0.037      0.062    20.0        0    -400     320.0
L2 E2 R2 a    sig sig ns     ns            True           none         False    0.199     0.329    0.037      0.062    20.0        0    -400     320.0
L2 E2 R2 b    sig sig ns     ns            True           none         False    0.199     0.329    0.037      0.062    20.0        0    -400     320.0
combinations with every part tested and holding: ['L2 E2 R1 a', 'L2 E2 R1 b', 'L2 E2 R2 a', 'L2 E2 R2 b']
combinations where a tested part contradicts the claim: []

rows = []
for (L, E, R, S) in COMBOS:
    for ref in shield.REFERENCE_AREAS:
        t = PAIRS[(ref, L, E, R, S)]
        lagv = {}
        for sid, v in SP.items():
            if sid not in SETS[S]:
                continue
            e = v["alpha"].get((L, ref, E))
            if e is None or "lags" not in e or len(e["coh"]["event"]) <= 5:
                continue
            for a in e["lags"].columns:
                lagv.setdefault(a, []).append(e["lags"][a].values)
        res = []
        for a in shield.LAG_PAIRS[ref]:
            if a not in lagv or not len(lookup(t, [a])) or not lookup(t, [a])["sig_pos"].any():
                continue
            x = np.concatenate(lagv[a]); x = x[np.isfinite(x)]
            res.append({"pair": f"{ref}-{a}", "n": len(x), "median_lag_ms": 1000 * np.median(x),
                        "fraction_positive": np.mean(x > 0), "p": stats.wilcoxon(x).pvalue})
        if res:
            r = pd.DataFrame(res); r["p_bh"] = multipletests(r["p"], method="fdr_bh")[1]
            for x in r.itertuples():
                rows.append({"combo": f"{L} {E} {R} {S}", **x._asdict()})
lagt = pd.DataFrame(rows).drop(columns="Index")
lagt["not_significant"] = lagt["p_bh"] >= 0.05
print(lagt.to_string(index=False, float_format=lambda v: f"{v:.3g}"))
print("pairs not significant (straddle zero under the test):", int(lagt["not_significant"].sum()), "of", len(lagt))

     combo       pair    n  median_lag_ms  fraction_positive        p     p_bh  not_significant
L1 E1 R1 a  VISp-VISl   97          -2.12              0.402    0.102    0.102             True
L1 E1 R1 a VISp-VISal  159           13.6              0.962 1.77e-22 4.43e-22            False
L1 E1 R1 a   VISp-LGd 1054          -5.13              0.267 1.52e-65 7.58e-65            False
L1 E1 R1 a    VISp-LP  760         -0.687               0.48 4.22e-06 5.28e-06            False
L1 E1 R1 a    VISp-CP  438           10.5              0.742  2.8e-21 4.66e-21            False
L1 E1 R1 a    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L1 E1 R1 b  VISp-VISl   97          -2.12              0.402    0.102    0.102             True
L1 E1 R1 b VISp-VISal  159           13.6              0.962 1.77e-22 4.43e-22            False
L1 E1 R1 b   VISp-LGd 1054          -5.13              0.267 1.52e-65 7.58e-65            False
L1 E1 R1 b    VISp-LP  760         -0.687               0.48 4.22e-06 5.28e-06            False
L1 E1 R1 b    VISp-CP  438           10.5              0.742  2.8e-21 4.66e-21            False
L1 E1 R1 b    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L1 E1 R2 a  VISp-VISl   97          -2.12              0.402    0.102    0.102             True
L1 E1 R2 a VISp-VISal  159           13.6              0.962 1.77e-22 4.43e-22            False
L1 E1 R2 a   VISp-LGd 1054          -5.13              0.267 1.52e-65 7.58e-65            False
L1 E1 R2 a    VISp-LP  760         -0.687               0.48 4.22e-06 5.28e-06            False
L1 E1 R2 a    VISp-CP  438           10.5              0.742  2.8e-21 4.66e-21            False
L1 E1 R2 a    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L1 E1 R2 b  VISp-VISl   97          -2.12              0.402    0.102    0.102             True
L1 E1 R2 b VISp-VISal  159           13.6              0.962 1.77e-22 4.43e-22            False
L1 E1 R2 b   VISp-LGd 1054          -5.13              0.267 1.52e-65 7.58e-65            False
L1 E1 R2 b    VISp-LP  760         -0.687               0.48 4.22e-06 5.28e-06            False
L1 E1 R2 b    VISp-CP  438           10.5              0.742  2.8e-21 4.66e-21            False
L1 E1 R2 b    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L1 E2 R1 a  VISp-VISl  107          -2.12              0.411   0.0618   0.0618             True
L1 E2 R1 a VISp-VISal  199           13.5              0.965 5.49e-29  1.1e-28            False
L1 E2 R1 a   VISp-LGd 1143          -5.16               0.27 2.16e-69 8.63e-69            False
L1 E2 R1 a    VISp-LP  834          -0.66              0.483 2.83e-06 3.78e-06            False
L1 E2 R1 a    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L1 E2 R1 a     MOp-PO   25          -23.3               0.12 0.000103 0.000206            False
L1 E2 R1 b  VISp-VISl  107          -2.12              0.411   0.0618   0.0618             True
L1 E2 R1 b VISp-VISal  199           13.5              0.965 5.49e-29  1.1e-28            False
L1 E2 R1 b   VISp-LGd 1143          -5.16               0.27 2.16e-69 8.63e-69            False
L1 E2 R1 b    VISp-LP  834          -0.66              0.483 2.83e-06 3.78e-06            False
L1 E2 R1 b    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L1 E2 R1 b     MOp-PO   25          -23.3               0.12 0.000103 0.000206            False
L1 E2 R2 a  VISp-VISl  107          -2.12              0.411   0.0618   0.0618             True
L1 E2 R2 a VISp-VISal  199           13.5              0.965 5.49e-29 1.37e-28            False
L1 E2 R2 a   VISp-LGd 1143          -5.16               0.27 2.16e-69 1.08e-68            False
L1 E2 R2 a    VISp-LP  834          -0.66              0.483 2.83e-06 3.54e-06            False
L1 E2 R2 a    VISp-CP  493           10.3              0.736 1.52e-23 2.53e-23            False
L1 E2 R2 a    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L1 E2 R2 a     MOp-PO   25          -23.3               0.12 0.000103 0.000206            False
L1 E2 R2 b  VISp-VISl  107          -2.12              0.411   0.0618   0.0618             True
L1 E2 R2 b VISp-VISal  199           13.5              0.965 5.49e-29 1.37e-28            False
L1 E2 R2 b   VISp-LGd 1143          -5.16               0.27 2.16e-69 1.08e-68            False
L1 E2 R2 b    VISp-LP  834          -0.66              0.483 2.83e-06 3.54e-06            False
L1 E2 R2 b    VISp-CP  493           10.3              0.736 1.52e-23 2.53e-23            False
L1 E2 R2 b    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L1 E2 R2 b     MOp-PO   25          -23.3               0.12 0.000103 0.000206            False
L2 E1 R1 a  VISp-VISl  141           1.79              0.596 0.000218 0.000272            False
L2 E1 R1 a VISp-VISal  388           2.11              0.595 7.99e-11 1.33e-10            False
L2 E1 R1 a   VISp-LGd 1025           -4.2              0.314 3.47e-48 1.73e-47            False
L2 E1 R1 a    VISp-LP  678         -0.865              0.485  0.00107  0.00107            False
L2 E1 R1 a    VISp-CP  387           10.6              0.762  5.8e-22 1.45e-21            False
L2 E1 R1 a    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L2 E1 R1 b  VISp-VISl  141           1.79              0.596 0.000218 0.000272            False
L2 E1 R1 b VISp-VISal  388           2.11              0.595 7.99e-11 1.33e-10            False
L2 E1 R1 b   VISp-LGd 1025           -4.2              0.314 3.47e-48 1.73e-47            False
L2 E1 R1 b    VISp-LP  678         -0.865              0.485  0.00107  0.00107            False
L2 E1 R1 b    VISp-CP  387           10.6              0.762  5.8e-22 1.45e-21            False
L2 E1 R1 b    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L2 E1 R2 a  VISp-VISl  141           1.79              0.596 0.000218 0.000272            False
L2 E1 R2 a VISp-VISal  388           2.11              0.595 7.99e-11 1.33e-10            False
L2 E1 R2 a   VISp-LGd 1025           -4.2              0.314 3.47e-48 1.73e-47            False
L2 E1 R2 a    VISp-LP  678         -0.865              0.485  0.00107  0.00107            False
L2 E1 R2 a    VISp-CP  387           10.6              0.762  5.8e-22 1.45e-21            False
L2 E1 R2 a    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L2 E1 R2 b  VISp-VISl  141           1.79              0.596 0.000218 0.000272            False
L2 E1 R2 b VISp-VISal  388           2.11              0.595 7.99e-11 1.33e-10            False
L2 E1 R2 b   VISp-LGd 1025           -4.2              0.314 3.47e-48 1.73e-47            False
L2 E1 R2 b    VISp-LP  678         -0.865              0.485  0.00107  0.00107            False
L2 E1 R2 b    VISp-CP  387           10.6              0.762  5.8e-22 1.45e-21            False
L2 E1 R2 b    MOp-SSp   50          -6.23               0.34  0.00833  0.00833            False
L2 E2 R1 a  VISp-VISl  169            1.8              0.592 0.000128  0.00016            False
L2 E2 R1 a VISp-VISal  457           2.92              0.621 7.28e-16 1.21e-15            False
L2 E2 R1 a   VISp-LGd 1114          -4.08               0.32 2.42e-49 1.21e-48            False
L2 E2 R1 a    VISp-LP  740         -0.168              0.497   0.0103   0.0103            False
L2 E2 R1 a    VISp-CP  439           10.5              0.754 3.01e-24 7.52e-24            False
L2 E2 R1 a    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L2 E2 R1 a     MOp-PO   25          -19.1               0.12 0.000217 0.000434            False
L2 E2 R1 b  VISp-VISl  169            1.8              0.592 0.000128  0.00016            False
L2 E2 R1 b VISp-VISal  464           2.98              0.625 1.59e-16 2.65e-16            False
L2 E2 R1 b   VISp-LGd 1114          -4.08               0.32 2.42e-49 1.21e-48            False
L2 E2 R1 b    VISp-LP  740         -0.168              0.497   0.0103   0.0103            False
L2 E2 R1 b    VISp-CP  439           10.5              0.754 3.01e-24 7.52e-24            False
L2 E2 R1 b    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L2 E2 R1 b     MOp-PO   25          -19.1               0.12 0.000217 0.000434            False
L2 E2 R2 a  VISp-VISl  169            1.8              0.592 0.000128  0.00016            False
L2 E2 R2 a VISp-VISal  457           2.92              0.621 7.28e-16 1.21e-15            False
L2 E2 R2 a   VISp-LGd 1114          -4.08               0.32 2.42e-49 1.21e-48            False
L2 E2 R2 a    VISp-LP  740         -0.168              0.497   0.0103   0.0103            False
L2 E2 R2 a    VISp-CP  439           10.5              0.754 3.01e-24 7.52e-24            False
L2 E2 R2 a    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L2 E2 R2 a     MOp-PO   25          -19.1               0.12 0.000217 0.000434            False
L2 E2 R2 b  VISp-VISl  169            1.8              0.592 0.000128  0.00016            False
L2 E2 R2 b VISp-VISal  464           2.98              0.625 1.59e-16 2.65e-16            False
L2 E2 R2 b   VISp-LGd 1114          -4.08               0.32 2.42e-49 1.21e-48            False
L2 E2 R2 b    VISp-LP  740         -0.168              0.497   0.0103   0.0103            False
L2 E2 R2 b    VISp-CP  439           10.5              0.754 3.01e-24 7.52e-24            False
L2 E2 R2 b    MOp-SSp   69          -5.86              0.319  0.00236  0.00236            False
L2 E2 R2 b     MOp-PO   25          -19.1               0.12 0.000217 0.000434            False
pairs not significant (straddle zero under the test): 8 of 102

A3 = {"ACA": ["ACA", "ACAd", "ACAv"], "RSP": ["RSP", "RSPd", "RSPv", "RSPagl"], "CP": ["CP"], "LGd": ["LGd"],
      "LP": ["LP"], "VISl": ["VISl", "VISlm"]}
print(show_pairs("VISp", A3).to_string(index=False))
rows = []
for (L, E, R, S) in COMBOS:
    t = PAIRS[("VISp", L, E, R, S)]
    for mode in ("any", "all"):
        st = {a: area_state(t, n, mode) for a, n in A3.items()}
        rows.append({"combo": f"{L} {E} {R} {S}", "subdivisions": mode, **st, "holds": all(s == "sig" for s in st.values())})
f7v = pd.DataFrame(rows)
print(f7v.to_string(index=False))
print("combinations where all six are significant and positive:", int(f7v["holds"].sum()), "of", len(f7v))
print("\nall VISp pairs, L1 E2 R2 a:")
print(PAIRS[("VISp", "L1", "E2", "R2", "a")].round(4).to_string(index=False))

     combo                                            ACA                                                                      RSP                   CP                    LGd                    LP                                            VISl
L1 E1 R1 a    ACAd +0.028 q=0.00032*; ACAv +0.013 q=0.091                           RSPagl +0.040 q=0.002*; RSPv +0.062 q=8.8e-10*  CP +0.021 q=0.0097*  LGd +0.104 q=2.2e-55*    LP +0.068 q=1e-20* VISl +0.158 q=3.6e-10*; VISlm +0.130 q=1.8e-08*
L1 E1 R1 b    ACAd +0.028 q=0.00032*; ACAv +0.013 q=0.091                           RSPagl +0.040 q=0.002*; RSPv +0.062 q=8.8e-10*  CP +0.021 q=0.0097*  LGd +0.104 q=2.2e-55*    LP +0.068 q=1e-20* VISl +0.158 q=3.6e-10*; VISlm +0.130 q=1.8e-08*
L1 E1 R2 a ACAd +0.071 q=4.8e-22*; ACAv +0.036 q=4.7e-12*                         RSPagl +0.044 q=0.00023*; RSPv +0.137 q=1.5e-32* CP +0.107 q=4.5e-40* LGd +0.347 q=5.5e-172* LP +0.266 q=3.6e-122*   VISl +0.375 q=4e-17*; VISlm +0.306 q=9.4e-17*
L1 E1 R2 b ACAd +0.071 q=4.8e-22*; ACAv +0.036 q=4.7e-12*                         RSPagl +0.044 q=0.00023*; RSPv +0.137 q=1.5e-32* CP +0.107 q=4.5e-40* LGd +0.347 q=5.5e-172* LP +0.266 q=3.6e-122*   VISl +0.375 q=4e-17*; VISlm +0.306 q=9.4e-17*
L1 E2 R1 a   ACAd +0.021 q=0.0027*; ACAv +0.017 q=0.0013*     RSPagl +0.045 q=0.0004*; RSPd +0.074 q=0.057; RSPv +0.065 q=2.3e-11*      CP +0.014 q=0.1  LGd +0.104 q=3.4e-63*  LP +0.064 q=5.8e-19* VISl +0.151 q=4.5e-11*; VISlm +0.123 q=1.6e-08*
L1 E2 R1 b   ACAd +0.021 q=0.0027*; ACAv +0.017 q=0.0013*     RSPagl +0.045 q=0.0004*; RSPd +0.074 q=0.057; RSPv +0.065 q=2.3e-11*      CP +0.014 q=0.1  LGd +0.104 q=3.4e-63*  LP +0.064 q=5.8e-19* VISl +0.151 q=4.5e-11*; VISlm +0.123 q=1.6e-08*
L1 E2 R2 a ACAd +0.054 q=6.5e-16*; ACAv +0.034 q=6.8e-13* RSPagl +0.061 q=2.6e-06*; RSPd +0.113 q=0.00071*; RSPv +0.133 q=3.2e-34* CP +0.103 q=2.3e-42* LGd +0.343 q=1.2e-186* LP +0.264 q=2.1e-134* VISl +0.346 q=9.6e-19*; VISlm +0.294 q=1.1e-20*
L1 E2 R2 b ACAd +0.054 q=6.5e-16*; ACAv +0.034 q=6.8e-13* RSPagl +0.061 q=2.6e-06*; RSPd +0.113 q=0.00071*; RSPv +0.133 q=3.2e-34* CP +0.103 q=2.3e-42* LGd +0.343 q=1.2e-186* LP +0.264 q=2.1e-134* VISl +0.346 q=9.6e-19*; VISlm +0.294 q=1.1e-20*
L2 E1 R1 a  ACAd +0.047 q=2.7e-06*; ACAv +0.017 q=0.0049*                             RSPagl +0.025 q=0.12; RSPv +0.079 q=4.8e-13*  CP +0.028 q=0.0049*  LGd +0.109 q=3.8e-58*  LP +0.091 q=3.6e-27*                            VISl +0.124 q=2e-10*
L2 E1 R1 b  ACAd +0.047 q=2.7e-06*; ACAv +0.017 q=0.0049*                             RSPagl +0.025 q=0.12; RSPv +0.079 q=4.8e-13*  CP +0.028 q=0.0049*  LGd +0.109 q=3.8e-58*  LP +0.091 q=3.6e-27*                            VISl +0.124 q=2e-10*
L2 E1 R2 a ACAd +0.108 q=1.6e-24*; ACAv +0.037 q=1.2e-13*                         RSPagl +0.057 q=1.2e-05*; RSPv +0.150 q=1.9e-34* CP +0.109 q=3.1e-37* LGd +0.353 q=1.7e-167* LP +0.305 q=1.1e-110*                          VISl +0.343 q=2.1e-24*
L2 E1 R2 b ACAd +0.108 q=1.6e-24*; ACAv +0.037 q=1.2e-13*                         RSPagl +0.057 q=1.2e-05*; RSPv +0.150 q=1.9e-34* CP +0.109 q=3.1e-37* LGd +0.353 q=1.7e-167* LP +0.305 q=1.1e-110*                          VISl +0.343 q=2.1e-24*
L2 E2 R1 a ACAd +0.045 q=9.5e-07*; ACAv +0.018 q=2.4e-05*                            RSPagl +0.025 q=0.093; RSPv +0.073 q=1.1e-13* CP +0.030 q=0.00066*  LGd +0.107 q=2.6e-66*  LP +0.086 q=5.2e-29*                          VISl +0.137 q=3.7e-13*
L2 E2 R1 b ACAd +0.045 q=9.9e-07*; ACAv +0.018 q=2.4e-05*       RSPagl +0.032 q=0.022*; RSPd +0.038 q=0.34; RSPv +0.073 q=1.1e-13* CP +0.030 q=0.00068*  LGd +0.107 q=2.6e-66*  LP +0.086 q=5.4e-29*                          VISl +0.137 q=3.9e-13*
L2 E2 R2 a ACAd +0.107 q=2.4e-25*; ACAv +0.039 q=3.3e-16*                           RSPagl +0.060 q=2e-07*; RSPv +0.143 q=9.6e-39* CP +0.107 q=1.2e-37* LGd +0.352 q=3.6e-182* LP +0.302 q=5.5e-121*                          VISl +0.347 q=6.5e-29*
L2 E2 R2 b ACAd +0.107 q=2.4e-25*; ACAv +0.039 q=3.4e-16* RSPagl +0.063 q=5.9e-08*; RSPd +0.092 q=0.00051*; RSPv +0.143 q=9.9e-39* CP +0.107 q=1.3e-37* LGd +0.352 q=3.7e-182* LP +0.302 q=5.6e-121*                          VISl +0.347 q=6.7e-29*

     combo subdivisions ACA RSP  CP LGd  LP VISl  holds
L1 E1 R1 a          any sig sig sig sig sig  sig   True
L1 E1 R1 a          all  ns sig sig sig sig  sig  False
L1 E1 R1 b          any sig sig sig sig sig  sig   True
L1 E1 R1 b          all  ns sig sig sig sig  sig  False
L1 E1 R2 a          any sig sig sig sig sig  sig   True
L1 E1 R2 a          all sig sig sig sig sig  sig   True
L1 E1 R2 b          any sig sig sig sig sig  sig   True
L1 E1 R2 b          all sig sig sig sig sig  sig   True
L1 E2 R1 a          any sig sig  ns sig sig  sig  False
L1 E2 R1 a          all sig  ns  ns sig sig  sig  False
L1 E2 R1 b          any sig sig  ns sig sig  sig  False
L1 E2 R1 b          all sig  ns  ns sig sig  sig  False
L1 E2 R2 a          any sig sig sig sig sig  sig   True
L1 E2 R2 a          all sig sig sig sig sig  sig   True
L1 E2 R2 b          any sig sig sig sig sig  sig   True
L1 E2 R2 b          all sig sig sig sig sig  sig   True
L2 E1 R1 a          any sig sig sig sig sig  sig   True
L2 E1 R1 a          all sig  ns sig sig sig  sig  False
L2 E1 R1 b          any sig sig sig sig sig  sig   True
L2 E1 R1 b          all sig  ns sig sig sig  sig  False
L2 E1 R2 a          any sig sig sig sig sig  sig   True
L2 E1 R2 a          all sig sig sig sig sig  sig   True
L2 E1 R2 b          any sig sig sig sig sig  sig   True
L2 E1 R2 b          all sig sig sig sig sig  sig   True
L2 E2 R1 a          any sig sig sig sig sig  sig   True
L2 E2 R1 a          all sig  ns sig sig sig  sig  False
L2 E2 R1 b          any sig sig sig sig sig  sig   True
L2 E2 R1 b          all sig  ns sig sig sig  sig  False
L2 E2 R2 a          any sig sig sig sig sig  sig   True
L2 E2 R2 a          all sig sig sig sig sig  sig   True
L2 E2 R2 b          any sig sig sig sig sig  sig   True
L2 E2 R2 b          all sig sig sig sig sig  sig   True
combinations where all six are significant and positive: 22 of 32

all VISp pairs, L1 E2 R2 a:
               area  events  sessions  mean_diff      p  ccg_peak  ccg_peak_lag_ms  p_bh_all   p_bh  sig_pos
               ACAd     450         3     0.0536 0.0000    0.1398                0    0.0000 0.0000     True
               ACAv     652         4     0.0337 0.0000    0.0699               10    0.0000 0.0000     True
                CA1    1156         7     0.0053 0.0868    0.0230              -60    0.1100 0.1136    False
                CA3     610         5     0.0065 0.0640    0.0139              -50    0.0838 0.0870    False
                 CP     493         4     0.1030 0.0000    0.1960              -10    0.0000 0.0000     True
              CP-RF    1021         5     0.1270 0.0000    0.2188              -20    0.0000 0.0000     True
            CP-noRF    1067         6     0.1066 0.0000    0.1623              -20    0.0000 0.0000     True
                 DG     690         6     0.0285 0.0000    0.0452              -40    0.0000 0.0000     True
                 DP     476         4     0.0099 0.0316    0.0651              -20    0.0444 0.0467     True
                LGd    1143         5     0.3433 0.0000    0.4256                0    0.0000 0.0000     True
                 LP     834         3     0.2642 0.0000    0.3354                0    0.0000 0.0000     True
                LSc     258         1     0.0078 0.1170    0.0199              190    0.1435 0.1474    False
                LSr     696         6     0.0134 0.0017    0.0629              -40    0.0027 0.0030     True
                 MB     283         3     0.0026 0.7827    0.0443              -10    0.8262 0.8316    False
                MOp     751         4     0.0023 0.3437    0.0274               80    0.3958 0.4030    False
                MOs     696         6     0.0008 0.8110    0.0364               10    0.8329 0.8356    False
                 MS     202         1     0.0060 0.4586    0.0376              -50    0.4979 0.5030    False
                OLF     265         2     0.0116 0.0543    0.0450              -40    0.0737 0.0769    False
               ORBm      26         1    -0.0133 0.3666    0.0832                0    0.4097 0.4155    False
                 PL     670         5     0.0104 0.0078    0.0517              -20    0.0113 0.0120     True
                 PO     413         3     0.0002 0.9584    0.0564               20    0.9584 0.9584    False
               POST     125         2     0.0259 0.0068    0.1087                0    0.0103 0.0110     True
                PPT     258         1     0.2460 0.0000    0.2964                0    0.0000 0.0000     True
               ProS     265         2     0.0234 0.0002    0.0433              -50    0.0003 0.0003     True
             RSPagl     102         1     0.0609 0.0000    0.0863               20    0.0000 0.0000     True
               RSPd      25         2     0.1133 0.0004    0.2271                0    0.0007 0.0007     True
               RSPv     309         4     0.1333 0.0000    0.2538              -10    0.0000 0.0000     True
               SCig      26         1     0.0186 0.1499    0.0522             -130    0.1780 0.1820    False
                SGN     834         3     0.2347 0.0000    0.2709                0    0.0000 0.0000     True
                SSp    1651        11     0.0080 0.0019    0.0552               20    0.0030 0.0032     True
                STR     461         3     0.0198 0.0000    0.0259              -10    0.0000 0.0000     True
                SUB     746         3     0.0147 0.0000    0.0904              -90    0.0000 0.0000     True
Track not annotated     209         2     0.4065 0.0000    0.5114              -10    0.0000 0.0000     True
                VIS     360         2     0.4085 0.0000    0.5532                0    0.0000 0.0000     True
              VISal     199         2     0.3941 0.0000    0.4645              -10    0.0000 0.0000     True
              VISam     406         2     0.3451 0.0000    0.4347                0    0.0000 0.0000     True
               VISl     107         1     0.3460 0.0000    0.4770                0    0.0000 0.0000     True
              VISlm     120         2     0.2935 0.0000    0.4546                0    0.0000 0.0000     True

A4 = {"SSp": ["SSp"], "PO": ["PO"], "MRN": ["MRN"], "SCiw": ["SCiw"]}
print(show_pairs("MOp", A4).to_string(index=False))
rows = []
for (L, E, R, S) in COMBOS:
    t = PAIRS[("MOp", L, E, R, S)]
    st = {a: area_state(t, n) for a, n in A4.items()}
    rows.append({"combo": f"{L} {E} {R} {S}", **st, "holds": all(s == "sig" for s in st.values())})
f7m = pd.DataFrame(rows)
print(f7m.to_string(index=False))
print("\nall MOp pairs, L1 E2 R2 a:")
print(PAIRS[("MOp", "L1", "E2", "R2", "a")].round(4).to_string(index=False))

     combo                   SSp                   PO                   MRN                   SCiw
L1 E1 R1 a SSp +0.181 q=8.8e-07*               absent                absent                 absent
L1 E1 R1 b SSp +0.181 q=8.8e-07*               absent                absent                 absent
L1 E1 R2 a   SSp +0.259 q=8e-14*               absent                absent                 absent
L1 E1 R2 b   SSp +0.259 q=8e-14*               absent                absent                 absent
L1 E2 R1 a SSp +0.190 q=4.3e-09* PO +0.174 q=1.1e-05*   MRN +0.112 q=0.035*  SCiw +0.114 q=0.0038*
L1 E2 R1 b SSp +0.190 q=4.3e-09* PO +0.174 q=1.1e-05*   MRN +0.112 q=0.035*  SCiw +0.114 q=0.0038*
L1 E2 R2 a   SSp +0.250 q=1e-11* PO +0.195 q=5.3e-06* MRN +0.171 q=2.5e-05* SCiw +0.179 q=2.4e-07*
L1 E2 R2 b   SSp +0.250 q=1e-11* PO +0.195 q=5.3e-06* MRN +0.171 q=2.5e-05* SCiw +0.179 q=2.4e-07*
L2 E1 R1 a SSp +0.173 q=2.7e-10*               absent                absent                 absent
L2 E1 R1 b SSp +0.173 q=2.7e-10*               absent                absent                 absent
L2 E1 R2 a SSp +0.246 q=8.6e-12*               absent                absent                 absent
L2 E1 R2 b SSp +0.246 q=8.6e-12*               absent                absent                 absent
L2 E2 R1 a SSp +0.187 q=5.4e-11*  PO +0.147 q=0.0015*     MRN +0.082 q=0.11   SCiw +0.121 q=0.003*
L2 E2 R1 b SSp +0.187 q=5.4e-11*  PO +0.147 q=0.0015*     MRN +0.082 q=0.11   SCiw +0.121 q=0.003*
L2 E2 R2 a SSp +0.247 q=2.2e-11* PO +0.169 q=1.6e-05* MRN +0.190 q=1.7e-05* SCiw +0.201 q=4.4e-09*
L2 E2 R2 b SSp +0.247 q=2.2e-11* PO +0.169 q=1.6e-05* MRN +0.190 q=1.7e-05* SCiw +0.201 q=4.4e-09*
     combo SSp     PO    MRN   SCiw  holds
L1 E1 R1 a sig absent absent absent  False
L1 E1 R1 b sig absent absent absent  False
L1 E1 R2 a sig absent absent absent  False
L1 E1 R2 b sig absent absent absent  False
L1 E2 R1 a sig    sig    sig    sig   True
L1 E2 R1 b sig    sig    sig    sig   True
L1 E2 R2 a sig    sig    sig    sig   True
L1 E2 R2 b sig    sig    sig    sig   True
L2 E1 R1 a sig absent absent absent  False
L2 E1 R1 b sig absent absent absent  False
L2 E1 R2 a sig absent absent absent  False
L2 E1 R2 b sig absent absent absent  False
L2 E2 R1 a sig    sig     ns    sig  False
L2 E2 R1 b sig    sig     ns    sig  False
L2 E2 R2 a sig    sig    sig    sig   True
L2 E2 R2 b sig    sig    sig    sig   True

all MOp pairs, L1 E2 R2 a:
   area  events  sessions  mean_diff      p  ccg_peak  ccg_peak_lag_ms  p_bh_all   p_bh  sig_pos
    CA1     108         5     0.0250 0.0116    0.0300              170    0.0513 0.0443     True
    CA3      25         2     0.0408 0.0551    0.0427              990    0.1514 0.1472    False
     CP      25         2     0.0098 0.9578    0.0546              -20    0.9578 0.9578    False
  CP-RF      69         3     0.0082 0.2700    0.0307               10    0.4569 0.4569    False
CP-noRF      76         4     0.0040 0.8684    0.0349             -230    0.9417 0.9417    False
     DG      64         4     0.0236 0.0501    0.0195              940    0.1514 0.1472    False
    LGd      76         4     0.0054 0.7719    0.0229             -160    0.8937 0.8627    False
     LP      60         2    -0.0058 0.6219    0.0293              210    0.7958 0.7732    False
    LSc      39         2     0.0155 0.4185    0.0332             -180    0.6138 0.6117    False
    LSr      32         1     0.0025 0.5608    0.0492              -10    0.7710 0.7610    False
    MRN      32         1     0.1713 0.0000    0.2314                0    0.0000 0.0000     True
     PO      25         2     0.1952 0.0000    0.2503               20    0.0000 0.0000     True
   POST      41         2     0.0188 0.0779    0.0353             -650    0.1558 0.1480    False
 RSPagl      48         2     0.0140 0.8989    0.0429              170    0.9417 0.9489    False
   RSPd      39         2     0.0285 0.0620    0.0603                0    0.1515 0.1472    False
   RSPv      32         1     0.0111 0.6511    0.0732               10    0.7958 0.7732    False
   SCiw      32         1     0.1795 0.0000    0.1620                0    0.0000 0.0000     True
    SGN      60         2     0.0164 0.0760    0.0310             -170    0.1558 0.1480    False
    SSp      69         3     0.2503 0.0000    0.3293                0    0.0000 0.0000     True
    SUB      44         1     0.0200 0.1378    0.0577             -270    0.2527 0.2381    False
    VIS      55         3     0.0250 0.0435    0.0449              150    0.1514 0.1514    False
   VISp     108         5     0.0023 0.3147    0.0211              680    0.4946 0.4983    False

tree = shield.read_code_csv("analysis_code/ccf_structure_tree_2017.csv")
hpf_id = tree.loc[tree["acronym"] == "HPF", "id"].iloc[0]
hpf = set(tree.loc[tree["structure_id_path"].astype(str).str.contains(f"/{hpf_id}/"), "acronym"])
rows = []
for (L, E, R, S) in COMBOS:
    for ref in shield.REFERENCE_AREAS:
        t = PAIRS[(ref, L, E, R, S)]
        h = t[t["area"].isin(hpf)] if len(t) else t
        rows.append({"combo": f"{L} {E} {R} {S}", "ref": ref, "hippocampal_pairs": len(h),
                     "significant_positive": ", ".join(h.loc[h["sig_pos"], "area"]) if len(h) else "",
                     "tested": ", ".join(f"{x.area} {x.mean_diff:+.3f} q={x.p_bh:.2g}" for x in h.itertuples())})
hp = pd.DataFrame(rows)
print(hp.to_string(index=False))
print("combinations and references with a significant positive hippocampal pair:", int((hp["significant_positive"] != "").sum()), "of", len(hp))

     combo  ref  hippocampal_pairs     significant_positive                                                                                                                        tested
L1 E1 R1 a VISp                  6                       DG          CA1 +0.003 q=0.48, CA3 +0.004 q=0.49, DG +0.019 q=8.3e-05, POST +0.006 q=0.71, ProS +0.009 q=0.29, SUB +0.003 q=0.64
L1 E1 R1 a  MOp                  3                                                                                                 CA1 +0.015 q=0.52, DG +0.036 q=0.27, SUB -0.004 q=0.95
L1 E1 R1 b VISp                  6                       DG          CA1 +0.003 q=0.48, CA3 +0.004 q=0.49, DG +0.019 q=8.3e-05, POST +0.006 q=0.71, ProS +0.009 q=0.29, SUB +0.003 q=0.64
L1 E1 R1 b  MOp                  3                                                                                                 CA1 +0.015 q=0.52, DG +0.036 q=0.27, SUB -0.004 q=0.95
L1 E1 R2 a VISp                  6       CA1, DG, ProS, SUB   CA1 +0.007 q=0.016, CA3 +0.008 q=0.073, DG +0.031 q=6.2e-14, POST +0.023 q=0.063, ProS +0.022 q=8e-05, SUB +0.018 q=3.9e-07
L1 E1 R2 a  MOp                  3                                                                                                 CA1 +0.024 q=0.14, DG +0.019 q=0.62, SUB +0.017 q=0.15
L1 E1 R2 b VISp                  6       CA1, DG, ProS, SUB   CA1 +0.007 q=0.016, CA3 +0.008 q=0.073, DG +0.031 q=6.2e-14, POST +0.023 q=0.063, ProS +0.022 q=8e-05, SUB +0.018 q=3.9e-07
L1 E1 R2 b  MOp                  3                                                                                                 CA1 +0.024 q=0.14, DG +0.019 q=0.62, SUB +0.017 q=0.15
L1 E2 R1 a VISp                  6                 DG, ProS          CA1 -0.001 q=0.65, CA3 +0.006 q=0.15, DG +0.015 q=0.0013, POST +0.014 q=0.42, ProS +0.011 q=0.027, SUB +0.002 q=0.65
L1 E2 R1 a  MOp                  5                                                           CA1 +0.019 q=0.13, CA3 +0.030 q=0.23, DG +0.021 q=0.2, POST -0.003 q=0.86, SUB +0.007 q=0.77
L1 E2 R1 b VISp                  6                 DG, ProS          CA1 -0.001 q=0.65, CA3 +0.006 q=0.15, DG +0.015 q=0.0013, POST +0.014 q=0.42, ProS +0.011 q=0.027, SUB +0.002 q=0.65
L1 E2 R1 b  MOp                  5                                                           CA1 +0.019 q=0.13, CA3 +0.030 q=0.23, DG +0.021 q=0.2, POST -0.003 q=0.86, SUB +0.007 q=0.77
L1 E2 R2 a VISp                  6      DG, POST, ProS, SUB   CA1 +0.005 q=0.11, CA3 +0.006 q=0.087, DG +0.029 q=6.5e-13, POST +0.026 q=0.011, ProS +0.023 q=0.0003, SUB +0.015 q=4.9e-05
L1 E2 R2 a  MOp                  5                      CA1                                CA1 +0.025 q=0.044, CA3 +0.041 q=0.15, DG +0.024 q=0.15, POST +0.019 q=0.15, SUB +0.020 q=0.24
L1 E2 R2 b VISp                  6      DG, POST, ProS, SUB   CA1 +0.005 q=0.11, CA3 +0.006 q=0.087, DG +0.029 q=6.5e-13, POST +0.026 q=0.011, ProS +0.023 q=0.0003, SUB +0.015 q=4.9e-05
L1 E2 R2 b  MOp                  5                      CA1                                CA1 +0.025 q=0.044, CA3 +0.041 q=0.15, DG +0.024 q=0.15, POST +0.019 q=0.15, SUB +0.020 q=0.24
L2 E1 R1 a VISp                  6                       DG           CA1 -0.001 q=0.68, CA3 -0.001 q=0.81, DG +0.012 q=0.041, POST +0.005 q=0.81, ProS +0.011 q=0.089, SUB +0.008 q=0.16
L2 E1 R1 a  MOp                  3                                                                                                  CA1 +0.005 q=0.58, DG +0.037 q=0.2, SUB +0.007 q=0.74
L2 E1 R1 b VISp                  6                       DG           CA1 -0.001 q=0.68, CA3 -0.001 q=0.81, DG +0.012 q=0.041, POST +0.005 q=0.81, ProS +0.011 q=0.089, SUB +0.008 q=0.16
L2 E1 R1 b  MOp                  3                                                                                                  CA1 +0.005 q=0.58, DG +0.037 q=0.2, SUB +0.007 q=0.74
L2 E1 R2 a VISp                  6 CA3, DG, POST, ProS, SUB  CA1 +0.005 q=0.075, CA3 +0.011 q=0.0064, DG +0.032 q=1e-12, POST +0.046 q=0.008, ProS +0.023 q=0.00021, SUB +0.022 q=8.8e-08
L2 E1 R2 a  MOp                  3                  CA1, DG                                                                     CA1 +0.038 q=0.0097, DG +0.041 q=0.014, SUB +0.020 q=0.49
L2 E1 R2 b VISp                  6 CA3, DG, POST, ProS, SUB  CA1 +0.005 q=0.075, CA3 +0.011 q=0.0064, DG +0.032 q=1e-12, POST +0.046 q=0.008, ProS +0.023 q=0.00021, SUB +0.022 q=8.8e-08
L2 E1 R2 b  MOp                  3                  CA1, DG                                                                     CA1 +0.038 q=0.0097, DG +0.041 q=0.014, SUB +0.020 q=0.49
L2 E2 R1 a VISp                  6                  DG, SUB           CA1 +0.001 q=0.94, CA3 +0.002 q=0.75, DG +0.013 q=0.002, POST -0.000 q=0.94, ProS +0.011 q=0.07, SUB +0.008 q=0.026
L2 E2 R1 a  MOp                  5                                                         CA1 +0.012 q=0.25, CA3 +0.039 q=0.096, DG +0.028 q=0.11, POST -0.009 q=0.63, SUB +0.011 q=0.48
L2 E2 R1 b VISp                  6                  DG, SUB          CA1 +0.001 q=0.94, CA3 +0.002 q=0.7, DG +0.013 q=0.0021, POST -0.000 q=0.94, ProS +0.011 q=0.068, SUB +0.008 q=0.022
L2 E2 R1 b  MOp                  5                                                         CA1 +0.012 q=0.25, CA3 +0.039 q=0.096, DG +0.028 q=0.11, POST -0.009 q=0.63, SUB +0.011 q=0.48
L2 E2 R2 a VISp                  6  CA1, CA3, DG, ProS, SUB CA1 +0.009 q=0.0084, CA3 +0.010 q=0.014, DG +0.029 q=1.6e-13, POST +0.021 q=0.064, ProS +0.017 q=0.0041, SUB +0.022 q=5.5e-09
L2 E2 R2 a  MOp                  5                                                          CA1 +0.013 q=0.48, CA3 +0.026 q=0.48, DG +0.014 q=0.53, POST +0.012 q=0.77, SUB +0.011 q=0.61
L2 E2 R2 b VISp                  6  CA1, CA3, DG, ProS, SUB  CA1 +0.009 q=0.0082, CA3 +0.010 q=0.012, DG +0.029 q=1.7e-13, POST +0.021 q=0.063, ProS +0.017 q=0.004, SUB +0.022 q=5.4e-09
L2 E2 R2 b  MOp                  5                                                          CA1 +0.013 q=0.48, CA3 +0.026 q=0.48, DG +0.014 q=0.53, POST +0.012 q=0.77, SUB +0.011 q=0.61
combinations and references with a significant positive hippocampal pair: 20 of 32

cpu = units[units["passing"] & units["qc"] & (units["structure_acronym"] == "CP")].join(rf[["p_rf"]], how="inner")
cpu["rf_sig"] = cpu["p_rf"] < 0.01
r_, n_ = cpu[cpu["rf_sig"]], cpu[~cpu["rf_sig"]]
print("orientation check, median z (um): VISp units", units.loc[units["structure_acronym"] == "VISp", "z"].median(),
      "; SSp", units.loc[units["structure_acronym"] == "SSp", "z"].median(), "; midline structures (LSr, SFO)",
      units.loc[units["structure_acronym"].isin(["LSr", "SFO"]), "z"].median(), "(the midline is at 5700)")
print(f"QC CP units with an RF p value: {len(cpu)} ({cpu['session_id'].nunique()} sessions, {cpu['subject_id'].nunique()} mice); "
      f"CP-RF {len(r_)}, CP-noRF {len(n_)}")
for ax_, lab, alt in (("z", "medial-lateral (z)", "greater"), ("y", "dorsal-ventral (y)", "less")):
    p = stats.ranksums(r_[ax_], n_[ax_], alternative=alt).pvalue
    print(f"{lab}: CP-RF median {r_[ax_].median():.0f}, CP-noRF median {n_[ax_].median():.0f} um; one-sided rank-sum p = {p:.2e}"
          f" (two-sided {stats.ranksums(r_[ax_], n_[ax_]).pvalue:.2e})")
print(f"anterior-posterior (x): CP-RF median {r_['x'].median():.0f} um; CP units sampled from {cpu['x'].min():.0f} to {cpu['x'].max():.0f} "
      f"(median {cpu['x'].median():.0f}); CP-RF fraction by AP tercile:",
      cpu.groupby(pd.qcut(cpu["x"], 3))["rf_sig"].mean().round(3).tolist())
print("CP-RF units by mouse:", r_.groupby("subject_id").size().to_dict())

orientation check, median z (um): VISp units 3050.0 ; SSp 3150.0 ; midline structures (LSr, SFO) 5300.0 (the midline is at 5700)
QC CP units with an RF p value: 2376 (52 sessions, 21 mice); CP-RF 688, CP-noRF 1688
medial-lateral (z): CP-RF median 4075, CP-noRF median 4350 um; one-sided rank-sum p = 1.00e+00 (two-sided 1.82e-33)
dorsal-ventral (y): CP-RF median 3150, CP-noRF median 3500 um; one-sided rank-sum p = 1.89e-58 (two-sided 3.78e-58)
anterior-posterior (x): CP-RF median 5388 um; CP units sampled from 3700 to 6500 (median 5038); CP-RF fraction by AP tercile: [0.102, 0.325, 0.447]
CP-RF units by mouse: {'607660': 57, '608671': 6, '608672': 19, '612090': 28, '614608': 57, '615047': 59, '615048': 49, '615563': 25, '615564': 55, '623319': 8, '623784': 63, '623786': 49, '626279': 8, '632295': 33, '636740': 23, '637483': 20, '638387': 14, '640887': 58, '640890': 8, '642504': 49}

f8 = pd.concat([v["fig8d"] for v in SP.values() if "fig8d" in v])
u8 = units.loc[f8.index].join(rf[["p_rf"]], how="left")
# every CP unit of the presence-0.9 set, noise ones included, has an RF p value (shield.py), so the
# variant without a noise filter splits noise CP units by RF as the code does
print("CP units in the five sessions without an RF p value:", int((u8["structure_acronym"].eq("CP") & u8["p_rf"].isna()).sum()))
res8 = {}
for filt in ("qc", "presence 0.9 (code variant)"):
    keep = u8["qc"] if filt == "qc" else pd.Series(True, index=u8.index)
    cp_rf = u8.index[keep & (u8["structure_acronym"] == "CP") & (u8["p_rf"] < 0.01)]
    cp_no = u8.index[keep & (u8["structure_acronym"] == "CP") & ~(u8["p_rf"] < 0.01)]
    v1i = u8.index[keep & (u8["isi_label"] == "VISp")]
    def prep(ix):
        X = f8.loc[ix].values.astype(float)
        X = X - X[:, 200:250].mean()
        return shield.exponential_convolve(X, tau=2, axis=1)
    Xr, Xn, Xv = prep(cp_rf), prep(cp_no), prep(v1i)
    p = np.array([stats.ranksums(a, b).pvalue for a, b in zip(Xr.T, Xn.T)])
    sig = multipletests(p, method="fdr_bh")[0]
    t = np.arange(750) - 250
    win = (t >= 0) & (t < 250)
    higher = Xr.mean(0) > Xn.mean(0)
    res8[filt] = dict(Xr=Xr, Xn=Xn, Xv=Xv, sig=sig, t=t)
    print(f"{filt}: CP-RF {len(cp_rf)}, CP-noRF {len(cp_no)}, VISp {len(v1i)} units; significant bins 0-250 ms: {int((sig & win).sum())}, "
          f"of which CP-RF higher: {int((sig & win & higher).sum())}; significant bins before onset: {int((sig & (t < 0)).sum())}; "
          f"first significant bin {t[sig][0] if sig.any() else None} ms")
    print("   units by session (CP-RF, CP-noRF, VISp):",
          {int(s): (int((u8.loc[cp_rf, 'session_id'] == s).sum()), int((u8.loc[cp_no, 'session_id'] == s).sum()), int((u8.loc[v1i, 'session_id'] == s).sum())) for s in shield.CP_RF_SESSIONS})
R8 = res8["qc"]
fig, ax = plt.subplots(figsize=(6, 3)); ax2 = ax.twinx()
ax.plot(R8["t"], R8["Xr"].mean(0), "r", label="CP-RF"); ax.plot(R8["t"], R8["Xn"].mean(0), "b", label="CP-noRF")
ax2.plot(R8["t"], R8["Xv"].mean(0), "k", label="VISp"); ax.set_xlim(-50, 400); ax.legend(loc="upper left"); ax2.legend(loc="upper right")
ax.set_xlabel("ms from image onset"); ax.set_ylabel("CP rate (Hz)"); ax2.set_ylabel("VISp rate (Hz)"); plt.show()

CP units in the five sessions without an RF p value: 0
qc: CP-RF 70, CP-noRF 162, VISp 218 units; significant bins 0-250 ms: 53, of which CP-RF higher: 53; significant bins before onset: 0; first significant bin 55 ms
   units by session (CP-RF, CP-noRF, VISp): {1182628226: (16, 48, 57), 1182871514: (22, 13, 42), 1183071525: (19, 23, 56), 1194850630: (12, 23, 45), 1195036753: (1, 55, 18)}

presence 0.9 (code variant): CP-RF 79, CP-noRF 206, VISp 239 units; significant bins 0-250 ms: 78, of which CP-RF higher: 78; significant bins before onset: 0; first significant bin 51 ms
   units by session (CP-RF, CP-noRF, VISp): {1182628226: (18, 57, 67), 1182871514: (23, 17, 46), 1183071525: (23, 36, 60), 1194850630: (14, 24, 48), 1195036753: (1, 72, 18)}

def latencies(X, t):
    m = X.mean(0); w = (t >= 0) & (t < 250)
    mw, tw = m[w], t[w]
    pk = tw[np.argmax(mw)]
    half = tw[np.argmax(mw >= mw.max() / 2)]
    return pk, half
for filt, R8_ in res8.items():
    lv, lr, ln = (latencies(R8_[k], R8_["t"]) for k in ("Xv", "Xr", "Xn"))
    print(f"{filt}: (a) peak  VISp {lv[0]} ms, CP-RF {lr[0]} ms, CP-noRF {ln[0]} ms -> CP-RF minus VISp {lr[0] - lv[0]}, CP-noRF minus CP-RF {ln[0] - lr[0]}")
    print(f"{' ' * len(filt)}  (b) half  VISp {lv[1]} ms, CP-RF {lr[1]} ms, CP-noRF {ln[1]} ms -> CP-RF minus VISp {lr[1] - lv[1]}, CP-noRF minus CP-RF {ln[1] - lr[1]}")

qc: (a) peak  VISp 48 ms, CP-RF 82 ms, CP-noRF 91 ms -> CP-RF minus VISp 34, CP-noRF minus CP-RF 9
    (b) half  VISp 38 ms, CP-RF 54 ms, CP-noRF 74 ms -> CP-RF minus VISp 16, CP-noRF minus CP-RF 20
presence 0.9 (code variant): (a) peak  VISp 48 ms, CP-RF 82 ms, CP-noRF 91 ms -> CP-RF minus VISp 34, CP-noRF minus CP-RF 9
                             (b) half  VISp 38 ms, CP-RF 54 ms, CP-noRF 74 ms -> CP-RF minus VISp 16, CP-noRF minus CP-RF 20

from scipy.signal import argrelextrema
rows = []
for L in ("L1", "L2"):
    for E in ("E1", "E2"):
        for S in ("a", "b"):
            tr = []
            for sid, v in SP.items():
                e = v["alpha"].get((L, "VISp", E))
                if sid not in SETS[S] or e is None or "trig_cp" not in e or len(e["coh"]["event"]) <= 5 or "CP" not in v["areas"][L]:
                    continue
                if len(e["trig_cp"]):
                    tr.append(e["trig_cp"])
            if not tr:
                continue
            T = shield.exponential_convolve(np.concatenate(tr), tau=1, axis=2)
            tt = np.arange(250) * 0.01 - 0.5
            vis, crf, cno = T[:, 0], T[:, 1], T[:, 2]
            p = np.array([stats.ranksums(a, b).pvalue for a, b in zip(crf.T, cno.T)])
            sig = multipletests(p, method="fdr_bh")[0]
            pw = np.array([stats.wilcoxon(a, b).pvalue for a, b in zip(crf.T, cno.T)])
            sigw = multipletests(pw, method="fdr_bh")[0]
            vm = vis.mean(0)
            w01 = (tt >= 0) & (tt <= 1)
            troughs = [i for i in argrelextrema(vm, np.less, order=5)[0] if w01[i]]
            near_trough = np.zeros(250, bool)
            for i in troughs:
                near_trough[max(i - 2, 0): i + 3] = True
            base_r, base_n = crf[:, tt < 0].mean(), cno[:, tt < 0].mean()
            rng_r = np.ptp(crf.mean(0)[w01]); rng_n = np.ptp(cno.mean(0)[w01])
            rows.append({"L": L, "E": E, "set": S, "events": len(T), "sessions": len(tr), "VISp_troughs_s": [round(tt[i], 2) for i in troughs],
                         "sig_bins_ranksum": int(sig.sum()), "sig_bins_at_troughs": int((sig & near_trough).sum()),
                         "sig_bins_signedrank": int(sigw.sum()), "sig_at_troughs_signedrank": int((sigw & near_trough).sum()),
                         "range_CP_RF": rng_r, "range_CP_noRF": rng_n, "range_norm_CP_RF": rng_r / base_r, "range_norm_CP_noRF": rng_n / base_n})
f8e = pd.DataFrame(rows)
f8e["holds_a"] = (f8e["sig_bins_at_troughs"] > 0) & (f8e["range_CP_RF"] > f8e["range_CP_noRF"])
f8e["holds_b"] = (f8e["sig_bins_at_troughs"] > 0) & (f8e["range_norm_CP_RF"] > f8e["range_norm_CP_noRF"])
print(f8e.to_string(index=False, float_format=lambda v: f"{v:.3g}"))

 L  E set  events  sessions                VISp_troughs_s  sig_bins_ranksum  sig_bins_at_troughs  sig_bins_signedrank  sig_at_troughs_signedrank  range_CP_RF  range_CP_noRF  range_norm_CP_RF  range_norm_CP_noRF  holds_a  holds_b
L1 E1   a     438         4 [0.0, 0.23, 0.49, 0.74, 0.98]               112                   24                   98                         23         1.98           1.29               0.5               0.323     True     True
L1 E1   b     438         4 [0.0, 0.23, 0.49, 0.74, 0.98]               112                   24                   98                         23         1.98           1.29               0.5               0.323     True     True
L1 E2   a     493         4 [0.0, 0.23, 0.49, 0.73, 0.97]               117                   25                  102                         25         1.96           1.18               0.5               0.297     True     True
L1 E2   b     493         4 [0.0, 0.23, 0.49, 0.73, 0.97]               117                   25                  102                         25         1.96           1.18               0.5               0.297     True     True
L2 E1   a     387         5  [0.0, 0.23, 0.5, 0.74, 0.98]               210                   25                  198                         25         2.04           1.34              0.53               0.323     True     True
L2 E1   b     387         5  [0.0, 0.23, 0.5, 0.74, 0.98]               210                   25                  198                         25         2.04           1.34              0.53               0.323     True     True
L2 E2   a     439         5 [0.0, 0.23, 0.49, 0.74, 0.98]               222                   25                  216                         25         2.11           1.29             0.555               0.314     True     True
L2 E2   b     439         5 [0.0, 0.23, 0.49, 0.74, 0.98]               222                   25                  216                         25         2.11           1.29             0.555               0.314     True     True

rows = []
for (L, E, R, S) in COMBOS:
    dr, dn = [], []
    for sid, v in SP.items():
        e = v["alpha"].get((L, "VISp", E))
        if sid not in SETS[S] or e is None or "coh" not in e or len(e["coh"]["event"]) <= 5 or "CP" not in v["areas"][L]:
            continue
        nm = e["coh_names"]
        if "CP-RF" not in nm or "CP-noRF" not in nm:
            continue
        d = e["coh"]["event"] - e["coh"][R]
        dr.append(d[:, nm.index("CP-RF")]); dn.append(d[:, nm.index("CP-noRF")])
    a, b = np.concatenate(dr), np.concatenate(dn)
    ok = np.isfinite(a) & np.isfinite(b)
    w = stats.wilcoxon(a[ok], b[ok])
    rows.append({"combo": f"{L} {E} {R} {S}", "events": int(ok.sum()), "sessions": len(dr), "CP_RF_mean": a[ok].mean(),
                 "CP_noRF_mean": b[ok].mean(), "p": w.pvalue, "holds": (a[ok].mean() > b[ok].mean()) and w.pvalue < 0.05})
f8f = pd.DataFrame(rows)
print(f8f.to_string(index=False, float_format=lambda v: f"{v:.3g}"))

     combo  events  sessions  CP_RF_mean  CP_noRF_mean        p  holds
L1 E1 R1 a     438         4      0.0116        0.0177    0.441  False
L1 E1 R1 b     438         4      0.0116        0.0177    0.441  False
L1 E1 R2 a     438         4      0.0998        0.0507 6.47e-11   True
L1 E1 R2 b     438         4      0.0998        0.0507 6.47e-11   True
L1 E2 R1 a     493         4     0.00679       0.00214    0.364  False
L1 E2 R1 b     493         4     0.00679       0.00214    0.364  False
L1 E2 R2 a     493         4      0.0956        0.0416 6.32e-14   True
L1 E2 R2 b     493         4      0.0956        0.0416 6.32e-14   True
L2 E1 R1 a     387         5      0.0171        0.0189    0.693  False
L2 E1 R1 b     387         5      0.0171        0.0189    0.693  False
L2 E1 R2 a     387         5        0.11        0.0454 1.22e-14   True
L2 E1 R2 b     387         5        0.11        0.0454 1.22e-14   True
L2 E2 R1 a     439         5      0.0221        0.0175    0.754  False
L2 E2 R1 b     439         5      0.0221        0.0175    0.754  False
L2 E2 R2 a     439         5       0.101        0.0517 1.94e-10   True
L2 E2 R2 b     439         5       0.101        0.0517 1.94e-10   True

f = shield.open_nwb(sessions["asset_id"].iloc[0])
groups = []
f.visititems(lambda n, o: groups.append(n) if o.attrs.get("neurodata_type") and n.count("/") <= 2 else None)
print("groups with a neurodata type:", [g for g in groups if not g.startswith("specifications")])
print("probe device attributes:", {k: str(v) for k, v in f["general/devices/probeA"].attrs.items()})
print("electrode group attributes:", {k: str(v) for k, v in f["general/extracellular_ephys/probeA"].attrs.items()})
print("session metadata:", {k: str(v) for k, v in f["general/metadata"].attrs.items()})
f.close()
print("\ngenotypes:", sorted(sessions["genotype"].unique()))
print("code session table columns:", list(code_sessions.columns))
print("code probes table columns:", list(code_probes.columns))

groups with a neurodata type: ['acquisition/EyeTracking', 'acquisition/EyeTracking/corneal_reflection_tracking', 'acquisition/EyeTracking/eye_tracking', 'acquisition/EyeTracking/likely_blink', 'acquisition/EyeTracking/pupil_tracking', 'acquisition/v_in', 'acquisition/v_sig', 'general/devices/NP.1', 'general/devices/probeA', 'general/devices/probeB', 'general/devices/probeC', 'general/devices/probeD', 'general/devices/probeE', 'general/devices/probeF', 'general/extracellular_ephys/electrodes', 'general/extracellular_ephys/probeA', 'general/extracellular_ephys/probeB', 'general/extracellular_ephys/probeC', 'general/extracellular_ephys/probeD', 'general/extracellular_ephys/probeE', 'general/extracellular_ephys/probeF', 'general/metadata', 'general/subject', 'general/task_parameters', 'intervals/dynamic_routing_image_set_presentations', 'intervals/dynamic_routing_image_set_presentations/active', 'intervals/dynamic_routing_image_set_presentations/duration', 'intervals/dynamic_routing_image_set_presentations/end_frame', 'intervals/dynamic_routing_image_set_presentations/flashes_since_change', 'intervals/dynamic_routing_image_set_presentations/id', 'intervals/dynamic_routing_image_set_presentations/image_name', 'intervals/dynamic_routing_image_set_presentations/is_change', 'intervals/dynamic_routing_image_set_presentations/is_image_novel', 'intervals/dynamic_routing_image_set_presentations/omitted', 'intervals/dynamic_routing_image_set_presentations/rewarded', 'intervals/dynamic_routing_image_set_presentations/start_frame', 'intervals/dynamic_routing_image_set_presentations/start_time', 'intervals/dynamic_routing_image_set_presentations/stimulus_block', 'intervals/dynamic_routing_image_set_presentations/stimulus_name', 'intervals/dynamic_routing_image_set_presentations/stop_time', 'intervals/dynamic_routing_image_set_presentations/tags', 'intervals/dynamic_routing_image_set_presentations/tags_index', 'intervals/dynamic_routing_image_set_presentations/timeseries', 'intervals/dynamic_routing_image_set_presentations/timeseries_index', 'intervals/flash_250ms_presentations', 'intervals/flash_250ms_presentations/active', 'intervals/flash_250ms_presentations/color', 'intervals/flash_250ms_presentations/contrast', 'intervals/flash_250ms_presentations/duration', 'intervals/flash_250ms_presentations/end_frame', 'intervals/flash_250ms_presentations/id', 'intervals/flash_250ms_presentations/is_change', 'intervals/flash_250ms_presentations/is_image_novel', 'intervals/flash_250ms_presentations/omitted', 'intervals/flash_250ms_presentations/rewarded', 'intervals/flash_250ms_presentations/start_frame', 'intervals/flash_250ms_presentations/start_time', 'intervals/flash_250ms_presentations/stimulus_block', 'intervals/flash_250ms_presentations/stimulus_index', 'intervals/flash_250ms_presentations/stimulus_name', 'intervals/flash_250ms_presentations/stop_time', 'intervals/flash_250ms_presentations/tags', 'intervals/flash_250ms_presentations/tags_index', 'intervals/flash_250ms_presentations/timeseries', 'intervals/flash_250ms_presentations/timeseries_index', 'intervals/gabor_20_deg_250ms_presentations', 'intervals/gabor_20_deg_250ms_presentations/active', 'intervals/gabor_20_deg_250ms_presentations/contrast', 'intervals/gabor_20_deg_250ms_presentations/duration', 'intervals/gabor_20_deg_250ms_presentations/end_frame', 'intervals/gabor_20_deg_250ms_presentations/id', 'intervals/gabor_20_deg_250ms_presentations/is_change', 'intervals/gabor_20_deg_250ms_presentations/is_image_novel', 'intervals/gabor_20_deg_250ms_presentations/omitted', 'intervals/gabor_20_deg_250ms_presentations/orientation', 'intervals/gabor_20_deg_250ms_presentations/position_x', 'intervals/gabor_20_deg_250ms_presentations/position_y', 'intervals/gabor_20_deg_250ms_presentations/rewarded', 'intervals/gabor_20_deg_250ms_presentations/spatial_frequency', 'intervals/gabor_20_deg_250ms_presentations/start_frame', 'intervals/gabor_20_deg_250ms_presentations/start_time', 'intervals/gabor_20_deg_250ms_presentations/stimulus_block', 'intervals/gabor_20_deg_250ms_presentations/stimulus_index', 'intervals/gabor_20_deg_250ms_presentations/stimulus_name', 'intervals/gabor_20_deg_250ms_presentations/stop_time', 'intervals/gabor_20_deg_250ms_presentations/tags', 'intervals/gabor_20_deg_250ms_presentations/tags_index', 'intervals/gabor_20_deg_250ms_presentations/temporal_frequency', 'intervals/gabor_20_deg_250ms_presentations/timeseries', 'intervals/gabor_20_deg_250ms_presentations/timeseries_index', 'intervals/spontaneous_presentations', 'intervals/spontaneous_presentations/active', 'intervals/spontaneous_presentations/duration', 'intervals/spontaneous_presentations/end_frame', 'intervals/spontaneous_presentations/id', 'intervals/spontaneous_presentations/is_change', 'intervals/spontaneous_presentations/is_image_novel', 'intervals/spontaneous_presentations/omitted', 'intervals/spontaneous_presentations/rewarded', 'intervals/spontaneous_presentations/start_frame', 'intervals/spontaneous_presentations/start_time', 'intervals/spontaneous_presentations/stimulus_block', 'intervals/spontaneous_presentations/stimulus_name', 'intervals/spontaneous_presentations/stop_time', 'intervals/spontaneous_presentations/tags', 'intervals/spontaneous_presentations/tags_index', 'intervals/spontaneous_presentations/timeseries', 'intervals/spontaneous_presentations/timeseries_index', 'intervals/trials', 'intervals/trials/aborted', 'intervals/trials/auto_rewarded', 'intervals/trials/catch', 'intervals/trials/change_frame', 'intervals/trials/change_image_name', 'intervals/trials/change_time', 'intervals/trials/correct_reject', 'intervals/trials/false_alarm', 'intervals/trials/go', 'intervals/trials/hit', 'intervals/trials/id', 'intervals/trials/initial_image_name', 'intervals/trials/is_change', 'intervals/trials/is_sham_change', 'intervals/trials/lick_times', 'intervals/trials/lick_times_index', 'intervals/trials/miss', 'intervals/trials/no_reward_epoch', 'intervals/trials/omitted_reward', 'intervals/trials/response_latency', 'intervals/trials/response_time', 'intervals/trials/reward_time', 'intervals/trials/reward_volume', 'intervals/trials/start_time', 'intervals/trials/stop_time', 'intervals/trials/trial_length', 'processing/eye_tracking_rig_metadata', 'processing/eye_tracking_rig_metadata/eye_tracking_rig_metadata', 'processing/licking', 'processing/licking/licks', 'processing/optotagging', 'processing/optotagging/optogenetic_stimulation', 'processing/optotagging/optotagging', 'processing/rewards', 'processing/rewards/autorewarded', 'processing/rewards/volume', 'processing/running', 'processing/running/dx', 'processing/running/speed', 'processing/running/speed_unfiltered', 'processing/stimulus', 'processing/stimulus/timestamps', 'stimulus/templates/dynamic_routing_image_set', 'units', 'units/PT_ratio', 'units/amplitude', 'units/amplitude_cutoff', 'units/cluster_id', 'units/cumulative_drift', 'units/d_prime', 'units/firing_rate', 'units/id', 'units/isi_label', 'units/isi_violations', 'units/isolation_distance', 'units/l_ratio', 'units/local_index', 'units/max_drift', 'units/nn_hit_rate', 'units/nn_miss_rate', 'units/peak_channel_id', 'units/presence_ratio', 'units/quality', 'units/recovery_slope', 'units/repolarization_slope', 'units/silhouette_score', 'units/snr', 'units/spike_amplitudes', 'units/spike_amplitudes_index', 'units/spike_times', 'units/spike_times_index', 'units/spread', 'units/structure_layer', 'units/velocity_above', 'units/velocity_below', 'units/waveform_duration', 'units/waveform_mean', 'units/waveform_mean_index']
probe device attributes: {'description': 'Neuropixels 1.0 Probe', 'manufacturer': 'imec', 'namespace': 'ndx-aibs-ecephys', 'neurodata_type': 'EcephysProbe', 'object_id': 'ade98fd2-9111-470c-a0c1-2001a1f27d87', 'probe_id': '627', 'sampling_rate': '30000.1328473262'}
electrode group attributes: {'description': 'Ecephys Electrode Group', 'has_lfp_data': 'False', 'lfp_sampling_rate': 'nan', 'location': 'See electrode locations', 'namespace': 'ndx-aibs-ecephys', 'neurodata_type': 'EcephysElectrodeGroup', 'object_id': 'b8453190-8a0c-4959-a9b0-1d2f44a42516', 'probe_id': '627'}
session metadata: {'behavior_session_id': '1173237009', 'behavior_session_uuid': 'afc5511f-a84f-424b-a32b-0c76ffbad82f', 'equipment_name': 'NP.1', 'namespace': 'ndx-aibs-behavior-ophys', 'neurodata_type': 'BehaviorMetadata', 'object_id': '493f1382-563b-44c4-9011-8d1e7e2603f7', 'project_code': 'Not Available', 'session_type': 'EPHYS_2', 'stimulus_frame_rate': '60.0'}

genotypes: ['Ai32(RCL-ChR2(H134R)_EYFP)/wt', 'Sst-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt', 'Vip-IRES-Cre/wt;Ai32(RCL-ChR2(H134R)_EYFP)/wt', 'wt/wt']
code session table columns: ['ecephys_session_id', 'behavior_session_id', 'date_of_acquisition', 'equipment_name', 'session_type', 'mouse_id', 'genotype', 'sex', 'project_code', 'age_in_days', 'probes_annotated_in_ccf', 'unit_count', 'channel_count', 'structure_acronyms', 'session_number', 'novel_image_id', 'has_lfp']
code probes table columns: ['Unnamed: 0.1', 'Unnamed: 0', 'id', 'name', 'location', 'sampling_rate', 'lfp_sampling_rate', 'has_lfp_data', 'session', 'inserted']

Claim	Tested here
link-code-tables (added link)	tested
fig3b-surgical-outcomes, fig3c-isi-through-implant, fig3d-histology-intact, tracks-rare-bleeding	blocked, data missing
fig3fh-mri-deformation, fig3i-learning-similar	blocked, data missing
fig4b-units-areas-per-session, figs4ab-insertions-per-mouse, scale-mice-sessions-insertions, penetration-success	tested
fig4b-trajectory-consistency	blocked, metadata insufficient
fig4c-drift-discontinuities-rare	blocked, methods underspecified
fig4d-unit-drift, fig4d-qc-metrics-comparable, qc-pass-fraction-similar	tested
fig4e-unit-yield-stable, fig4e-vbn-yield-drops, fig4e-firing-rate-stable, fig4ef-rf-fraction-stable	tested
fig4f-rf-match-retinotopy	blocked, data missing
fig4g-visual-responses-similar, behavior-stable-across-days, figs4i-insertion-suppression	tested
figs4de-retraction-reduces-drift, figs4f-retraction-yield-similar, figs4gh-sharpened-probes	blocked, metadata insufficient
figs3d-duragel-yield-stable, fig5-bilateral-recordings, fig6bc-cortical-silencing, fig6ef-optotagging-mrn	blocked, data missing
fig7a to fig7f, fig7-hippocampus-not-engaged, figs6-replay (10 claims)	tested
fig8bc, fig8d (two), fig8e, fig8f (5 claims)	tested

Claim	Outcome
link-code-tables	partially reproduced (6 of 87 sessions hold other units than the code's table)
fig3b-surgical-outcomes	blocked, data missing
fig3c-isi-through-implant	blocked, data missing
fig3d-histology-intact	blocked, data missing
tracks-rare-bleeding	blocked, data missing
fig3fh-mri-deformation	blocked, data missing
fig3i-learning-similar	blocked, data missing
fig4b-units-areas-per-session	reproduced
figs4ab-insertions-per-mouse	partially reproduced (holes not testable)
scale-mice-sessions-insertions	reproduced
penetration-success	reproduced
fig4b-trajectory-consistency	blocked, metadata insufficient
fig4c-drift-discontinuities-rare	blocked, methods underspecified
fig4d-unit-drift	reproduced
fig4d-qc-metrics-comparable	not reproduced
qc-pass-fraction-similar	reproduced
fig4e-unit-yield-stable	reproduced
fig4e-vbn-yield-drops	reproduced
fig4e-firing-rate-stable	reproduced
fig4ef-rf-fraction-stable	reproduced
fig4f-rf-match-retinotopy	blocked, data missing
fig4g-visual-responses-similar	reproduced
behavior-stable-across-days	partially reproduced
figs4i-insertion-suppression	not reproduced
figs4de-retraction-reduces-drift	blocked, metadata insufficient
figs4f-retraction-yield-similar	blocked, metadata insufficient
figs4gh-sharpened-probes	blocked, metadata insufficient
figs3d-duragel-yield-stable	blocked, data missing
fig5-bilateral-recordings	blocked, data missing
fig6bc-cortical-silencing	blocked, data missing
fig6ef-optotagging-mrn	blocked, data missing
fig7a-mop-higher-frequency	reproduced
fig7b-low-arousal	partially reproduced
fig7b-visp-offset-locked	reproduced
figs6-replay	not reproduced
fig7cd-visp-coupling	blocked, methods underspecified
fig7cd-mop-coupling	blocked, methods underspecified
fig7e-lags-straddle-zero	not reproduced
fig7f-visp-network	blocked, methods underspecified
fig7f-mop-network	blocked, methods underspecified
fig7-hippocampus-not-engaged	not reproduced
fig8bc-cp-rf-location	partially reproduced
fig8d-cp-rf-larger-response	partially reproduced (effect holds; unit counts match only with the code's filter)
fig8d-response-latencies	blocked, methods underspecified
fig8e-cp-rf-more-modulated	reproduced
fig8f-cp-rf-higher-coherence	blocked, methods underspecified

Reproducing Bennett, Ouellette, Ramirez et al. (2024) from dandiset 001051¶

Data check¶

link-code-tables: the paper's session list and unit table against the dandiset¶

scale-mice-sessions-insertions¶

penetration-success¶

fig4b-units-areas-per-session¶

figs4ab-insertions-per-mouse¶

fig4d-unit-drift¶

fig4d-qc-metrics-comparable¶

qc-pass-fraction-similar¶

fig4e-unit-yield-stable¶

fig4e-vbn-yield-drops¶

fig4e-firing-rate-stable¶

fig4ef-rf-fraction-stable¶

fig4g-visual-responses-similar¶

behavior-stable-across-days¶

figs4i-insertion-suppression¶

Alpha events: detection and the shared method¶

fig7a-mop-higher-frequency¶

fig7b-low-arousal¶

fig7b-visp-offset-locked¶

figs6-replay¶

fig7cd-visp-coupling¶

fig7cd-mop-coupling¶

fig7e-lags-straddle-zero¶

fig7f-visp-network¶

fig7f-mop-network¶

fig7-hippocampus-not-engaged¶

fig8bc-cp-rf-location¶

fig8d-cp-rf-larger-response¶

fig8d-response-latencies¶

fig8e-cp-rf-more-modulated¶

fig8f-cp-rf-higher-coherence¶

Blocked claims: what the files hold¶

Summary¶