import hashlib
import os
import time

os.environ.setdefault("OMP_NUM_THREADS", "1")
import fsspec
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from scipy import stats

import pfc2

T0 = time.time()
raw, TAR = pfc2.read_raw()
with fsspec.open(pfc2.BASE.replace("data/", "checksums.md5")) as f:
    md5 = dict(reversed(line.split()) for line in f.read().decode().splitlines() if line.strip())
with fsspec.open(pfc2.BASE + "behaviors.zip", block_size=2**22) as f:
    h = hashlib.md5()
    while chunk := f.read(2**24):
        h.update(chunk)
print("data/behaviors.zip", md5["data/behaviors.zip"], "read", h.hexdigest())
for k in ["EE188_Behavior.mat", "EE188_Behavior_PFC.mat"]:
    print("data/EE0708/" + k, md5["data/EE0708/" + k], "read", hashlib.md5(raw[k]).hexdigest())
print("data/EE0708/EE.188.tar.gz", md5["data/EE0708/EE.188.tar.gz"], "(member EE188_example.mat read in the stream)")
S = pfc2.load_sessions(raw)
EX = pfc2._mat(raw["EE188_example.mat"])
print("loaded in", round(time.time() - T0), "s")

data/behaviors.zip 3a510fe0ef8a24ed4248d37bed3fd624 read 3a510fe0ef8a24ed4248d37bed3fd624
data/EE0708/EE188_Behavior.mat b7407ba84d282c33cf937f090f2ebb0c read b7407ba84d282c33cf937f090f2ebb0c
data/EE0708/EE188_Behavior_PFC.mat a3ca6359c45e2a8deeada4b8eb1d61c6 read a3ca6359c45e2a8deeada4b8eb1d61c6
data/EE0708/EE.188.tar.gz 115236aeb0b9890b4dd0a26fbabe0ed3 (member EE188_example.mat read in the stream)

loaded in 87 s

L23 = [k for k, s in S.items() if s.layer == "L2/3"]
L5 = [k for k, s in S.items() if s.layer == "L5"]
print(f"{'session':8s} {'rat':3s} {'layer':5s} {'cells':>5s} {'pyr/int/unk (file)':>18s} {'trials':>6s} {'L':>3s} {'R':>3s} "
      f"{'spikes':>8s} {'dur s':>6s} {'shanks':>6s}")
for k, s in S.items():
    ty = np.bincount(s.celltype_file, minlength=4)[1:]
    print(f"{k:8s} {s.rat:3s} {s.layer:5s} {s.n:5d} {str(ty):>18s} {s.ntrials:6d} {s.left.sum():3d} {(~s.left).sum():3d} "
          f"{len(s.t):8d} {s.duration:6.0f} {len(np.unique(s.shank)):6d}")
print("cells per rat:", {r: sum(s.n for s in S.values() if s.rat == r) for r in ("EE", "FF", "GG")})
print("sessions per layer:", {"L2/3": L23, "L5": L5})
allx = np.concatenate([s.x[s.x >= 0] for s in S.values()])
print(f"spike positions (spikeX): {allx.min():.4f} to {allx.max():.4f} (normalised length); "
      f"spikes outside trials marked -1: {np.mean(np.concatenate([s.x < 0 for s in S.values()])):.2f} of all")
# trial labels of spikes against the trial table
for k, s in S.items():
    assert s.max_ltr == s.left.sum() and s.max_rtr == (~s.left).sum()
    tt = s.trial >= 0
    st, en = s.NP[s.trial[tt], 0], s.NP[s.trial[tt], 2]
    inside = np.mean((s.t[tt] >= st - 0.05) & (s.t[tt] <= en + 0.05))
    d = s.whl[np.minimum((s.t[tt] * pfc2.WHL_RATE).astype(int), len(s.whl) - 1), 6]
    ok = d > 0
    r = np.corrcoef(s.x[tt][ok], d[ok] - 1)[0, 1]
    print(f"{k}: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): {inside:.4f}; "
          f"spikeX against the position file at the spike time r = {r:.3f}; trial durations "
          f"{np.min(s.NP[:, 2] - s.NP[:, 1]):.1f} to {np.max(s.NP[:, 2] - s.NP[:, 1]):.1f} s; "
          f"nose-poke {np.min(s.NP[:, 1] - s.NP[:, 0]):.2f} to {np.max(s.NP[:, 1] - s.NP[:, 0]):.1f} s")
w = S["EE188"].whl
print("position file EE188: x1", w[w[:, 0] > 0, 0].min(), "to", w[:, 0].max(), "y1", w[w[:, 1] > 0, 1].min(), "to",
      w[:, 1].max(), "(pixels); d", w[w[:, 6] > 0, 6].min(), "to", w[:, 6].max())

session  rat layer cells pyr/int/unk (file) trials   L   R   spikes  dur s shanks
EE049    EE  L2/3     44         [40  4  0]     19   9  10   193187   1789      8
EE088    EE  L2/3     28         [22  6  0]     31  11  20   114526   2165      7
EE152    EE  L2/3     61         [55  6  0]     36  15  21   389009   2552      8
EE165    EE  L2/3     70         [60 10  0]     42  19  23   469949   2944      8
EE188    EE  L2/3    117         [33 15 69]     42  18  24   763501   2812      8
EE196    EE  L2/3     66         [54 12  0]     48  21  27   760669   2844      8
EE210    EE  L2/3     60         [51  9  0]     53  25  28   278941   2412      8
FF103    FF  L5       42         [36  6  0]     39  19  20   323905   2867      6
FF130    FF  L5       64         [57  7  0]     46  23  23   434304   2030      7
FF167    FF  L5       40         [36  4  0]     43  21  22   242698   1978      7
GG069    GG  L2/3     14         [11  3  0]     40  20  20    49406   2376      5
GG097    GG  L2/3     12         [12  0  0]     38  20  18    19332   1750      4
GG149    GG  L2/3      9            [9 0 0]     31  16  15    29735   2846      6
cells per rat: {'EE': 446, 'FF': 146, 'GG': 35}
sessions per layer: {'L2/3': ['EE049', 'EE088', 'EE152', 'EE165', 'EE188', 'EE196', 'EE210', 'GG069', 'GG097', 'GG149'], 'L5': ['FF103', 'FF130', 'FF167']}
spike positions (spikeX): 0.0000 to 1.0000 (normalised length); spikes outside trials marked -1: 0.83 of all
EE049: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.999; trial durations 8.3 to 35.4 s; nose-poke 0.43 to 24.1 s
EE088: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.991; trial durations 4.0 to 11.1 s; nose-poke 0.49 to 267.3 s
EE152: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.998; trial durations 5.6 to 22.2 s; nose-poke 0.97 to 18.4 s
EE165: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.999; trial durations 5.5 to 12.6 s; nose-poke 0.77 to 12.4 s
EE188: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.999; trial durations 6.2 to 19.5 s; nose-poke 1.88 to 29.4 s
EE196: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.997; trial durations 3.6 to 12.7 s; nose-poke 1.49 to 16.9 s
EE210: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.997; trial durations 3.8 to 9.5 s; nose-poke 1.09 to 17.3 s
FF103: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.999; trial durations 5.2 to 9.5 s; nose-poke 0.03 to 315.6 s
FF130: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.998; trial durations 4.7 to 10.2 s; nose-poke 0.05 to 7.8 s
FF167: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.998; trial durations 4.9 to 9.9 s; nose-poke 0.15 to 15.9 s
GG069: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.999; trial durations 6.7 to 12.8 s; nose-poke 0.40 to 50.8 s
GG097: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 0.998; trial durations 4.2 to 8.4 s; nose-poke 0.19 to 61.8 s
GG149: labelled spikes inside their trial's [nose-poke start, end] (+-50 ms): 1.0000; spikeX against the position file at the spike time r = 1.000; trial durations 4.2 to 9.9 s; nose-poke 0.67 to 16.4 s
position file EE188: x1 142.0 to 571.8602 y1 69.33014 to 440.0 (pixels); d 1.0062149470873414 to 2.0

s = S["EE188"]
p105 = pfc2._mat(raw["EE188_Behavior_PFC.mat"])
c105 = p105["CellsInfo"][:, 0].astype(int)
print("EE188_example.mat cells:", s.n, "types", np.bincount(s.celltype_file, minlength=4)[1:],
      "| EE188_Behavior_PFC.mat cells:", len(c105), "types", np.bincount(p105["CellsInfo"][:, 2].astype(int), minlength=4)[1:])
print("cells of the example not in Behavior_PFC:", sorted(set(s.cells) - set(c105)),
      "| of Behavior_PFC not in the example:", sorted(set(c105) - set(s.cells)))
same = np.isin(s.cells, c105)
t1 = np.sort(np.asarray(p105["spiket"]).ravel())
print("spikes: example", len(s.t), "Behavior_PFC", len(t1), "; example spikes of the shared cells",
      int(np.isin(s.u, np.flatnonzero(same)).sum()))
# pairs of sessions with identical spike counts per cell (a sign of a recording deposited twice)
sig = {k: tuple(sorted(np.bincount(v.u))) for k, v in S.items()}
dups = [(a, b) for a in S for b in S if a < b and sig[a] == sig[b]]
print("sessions with identical spike count lists:", dups if dups else "none")
print("spikeph (all clusters of the day's sorting on the 8 PFC shanks) against the cells kept:")
for k in S:
    if k == "EE188":
        n_ph = p105["spikeph"].shape[0]
    else:
        n_ph = pfc2._mat(raw[k + "_Behavior_PFC.mat"])["spikeph"].shape[0]
    print(f"  {k}: {n_ph} clusters, {S[k].n} cells")

EE188_example.mat cells: 117 types [33 15 69] | EE188_Behavior_PFC.mat cells: 105 types [87 18  0]
cells of the example not in Behavior_PFC: [np.int64(33), np.int64(37), np.int64(54), np.int64(92), np.int64(134), np.int64(152), np.int64(183), np.int64(218), np.int64(250), np.int64(283), np.int64(285), np.int64(293)] | of Behavior_PFC not in the example: []
spikes: example 763501 Behavior_PFC 750910 ; example spikes of the shared cells 750910
sessions with identical spike count lists: none
spikeph (all clusters of the day's sorting on the 8 PFC shanks) against the cells kept:
  EE049: 78 clusters, 44 cells
  EE088: 185 clusters, 28 cells
  EE152: 241 clusters, 61 cells

  EE165: 234 clusters, 70 cells
  EE188: 297 clusters, 117 cells

  EE196: 133 clusters, 66 cells
  EE210: 109 clusters, 60 cells
  FF103: 51 clusters, 42 cells
  FF130: 68 clusters, 64 cells
  FF167: 43 clusters, 40 cells

  GG069: 19 clusters, 14 cells
  GG097: 19 clusters, 12 cells
  GG149: 18 clusters, 9 cells

def loo_centroid(X, y):
    X = (X - X.mean(0)) / (X.std(0) + 1e-9)
    hit = 0
    for i in range(len(y)):
        m = np.ones(len(y), bool)
        m[i] = False
        c1, c0 = X[m & y].mean(0), X[m & ~y].mean(0)
        hit += (np.sum((X[i] - c1) ** 2) < np.sum((X[i] - c0) ** 2)) == y[i]
    return hit / len(y)


def clock_check(s, shift):
    st = s.NP[:, 1] + shift
    en = s.NP[:, 2] + shift
    X = np.zeros((s.ntrials, s.n))
    for k, t0 in enumerate(en):
        a, b = np.searchsorted(s.t, [t0 - 3.0, t0])
        X[k] = np.bincount(s.u[a:b], minlength=s.n)
    acc_sp = loo_centroid(X, s.left)
    fr = ((st[:, None] + np.linspace(2, 4, 9)) * pfc2.WHL_RATE).astype(int)
    fr = np.clip(fr, 0, len(s.whl) - 1)
    acc_x = loo_centroid(s.whl[fr, 0], s.left)
    return acc_sp, acc_x


print(f"{'session':8s} {'spikes true':>11s} {'spikes +30s':>11s} {'x true':>7s} {'x +30s':>7s}")
for k, s in S.items():
    a0, x0 = clock_check(s, 0.0)
    a1, x1 = clock_check(s, 30.0)
    print(f"{k:8s} {a0:11.2f} {a1:11.2f} {x0:7.2f} {x1:7.2f}")

session  spikes true spikes +30s  x true  x +30s
EE049           0.79        0.58    0.53    0.42
EE088           0.94        0.65    0.52    0.42
EE152           0.94        0.67    0.56    0.50
EE165           0.90        0.57    0.57    0.55
EE188           0.95        0.40    0.74    0.29
EE196           1.00        0.65    0.69    0.58

EE210           1.00        0.38    0.77    0.47
FF103           0.95        0.62    0.62    0.59
FF130           0.96        0.80    0.63    0.30
FF167           0.98        0.58    0.67    0.60
GG069           0.82        0.50    0.88    0.47
GG097           0.61        0.29    0.68    0.45
GG149           0.94        0.58    0.71    0.55

n_cells = sum(s.n for s in S.values())
n_alt = n_cells - S["EE188"].n + len(c105)
print(f"cells over the 13 working memory sessions: {n_cells} (with EE188 as 105 cells: {n_alt}); paper 633, "
      f"within 633 +- 63.3: {abs(n_cells - 633) <= 63.3}")
print("per rat:", {r: sum(s.n for s in S.values() if s.rat == r) for r in ("EE", "FF", "GG")})
n_ph = sum((p105 if k == "EE188" else pfc2._mat(raw[k + "_Behavior_PFC.mat"]))["spikeph"].shape[0] for k in S)
print(f"all clusters listed in spikeph: {n_ph}")

# Which selection are the files' cells? The data description defines Cells as "the neuron # whose
# max firing rates are more than 1 Hz in the trials" (page 2); the files hold no variable of that
# name. Each cell's largest rate over the tenths of the maze on left and right trials (spikes over
# time spent, from the position file) is computed for every cell kept, and for the 12 cells of
# EE188_example.mat that EE188_Behavior_PFC.mat leaves out.
def max_trial_rate(s, cells):
    w = s.whl
    side_f = w[:, 4].astype(int)
    d = w[:, 6]
    inf = d > 0
    tenth_f = np.minimum(((d[inf] - 1) * 10).astype(int), 9)
    occ = np.zeros((3, 10))
    np.add.at(occ, (side_f[inf], tenth_f), 1 / pfc2.WHL_RATE)
    m = s.intrial
    side_s = np.where(s.left[np.maximum(s.trial, 0)], 2, 1)[m]
    tenth_s = np.minimum((s.x[m] * 10).astype(int), 9)
    cnt = np.zeros((s.n, 3, 10))
    np.add.at(cnt, (s.u[m], side_s, tenth_s), 1)
    with np.errstate(invalid="ignore", divide="ignore"):
        rate = np.where(occ[None, 1:] > 0, cnt[:, 1:] / occ[None, 1:], 0)
    return rate.reshape(s.n, -1).max(1)[cells]


for k, s in S.items():
    r = max_trial_rate(s, np.arange(s.n))
    extra = ""
    if k == "EE188":
        out = ~np.isin(s.cells, c105)
        extra = (f"; the 12 cells not in EE188_Behavior_PFC.mat: {np.sort(r[out]).round(2).tolist()} Hz, of the other 105 "
                 f"{np.sum(r[~out] > 1)} above 1 Hz")
    print(f"{k}: cells with a largest trial rate above 1 Hz {np.sum(r > 1)} of {s.n} (lowest {r.min():.2f} Hz){extra}")

cells over the 13 working memory sessions: 627 (with EE188 as 105 cells: 615); paper 633, within 633 +- 63.3: True
per rat: {'EE': 446, 'FF': 146, 'GG': 35}

all clusters listed in spikeph: 1495
EE049: cells with a largest trial rate above 1 Hz 44 of 44 (lowest 1.53 Hz)
EE088: cells with a largest trial rate above 1 Hz 28 of 28 (lowest 1.31 Hz)
EE152: cells with a largest trial rate above 1 Hz 59 of 61 (lowest 0.76 Hz)
EE165: cells with a largest trial rate above 1 Hz 66 of 70 (lowest 0.53 Hz)
EE188: cells with a largest trial rate above 1 Hz 107 of 117 (lowest 0.23 Hz); the 12 cells not in EE188_Behavior_PFC.mat: [0.23, 0.25, 0.39, 0.54, 0.72, 1.42, 1.71, 1.93, 1.96, 2.53, 2.63, 3.21] Hz, of the other 105 100 above 1 Hz
EE196: cells with a largest trial rate above 1 Hz 63 of 66 (lowest 0.64 Hz)
EE210: cells with a largest trial rate above 1 Hz 57 of 60 (lowest 0.69 Hz)
FF103: cells with a largest trial rate above 1 Hz 39 of 42 (lowest 0.81 Hz)
FF130: cells with a largest trial rate above 1 Hz 63 of 64 (lowest 0.96 Hz)
FF167: cells with a largest trial rate above 1 Hz 37 of 40 (lowest 0.48 Hz)
GG069: cells with a largest trial rate above 1 Hz 14 of 14 (lowest 1.04 Hz)
GG097: cells with a largest trial rate above 1 Hz 12 of 12 (lowest 1.56 Hz)
GG149: cells with a largest trial rate above 1 Hz 9 of 9 (lowest 1.29 Hz)

print(f"EE.188.tar.gz: {len(TAR)} members, for example "
      + ", ".join(f"{n.split('/')[-1]} ({sz / 1e6:.0f} MB)" for n, sz in TAR if n.endswith((".eeg", ".whl", ".whlrl", ".xml", ".res", "_example.mat"))))
print("member types:", sorted({n.split("/")[-1].split(".")[-2] if n.split("/")[-1].split(".")[-1].isdigit() else n.split(".")[-1]
                               for n, sz in TAR if sz > 0}))
s = S["EE188"]
sens = raw["EE188 nose-poke sensors"]
PK = {1: pfc2.poke_intervals(sens[:, 0]), 2: pfc2.poke_intervals(sens[:, 1])}
w = s.whl
t_f = np.arange(len(w)) / pfc2.WHL_RATE
ok_f = w[:, 0] > 0
two = w[:, 2] > 0
px = np.where(two, (w[:, 0] + w[:, 2]) / 2, w[:, 0])
py = np.where(two, (w[:, 1] + w[:, 3]) / 2, w[:, 1])
ends = {}
for side in (1, 2):
    f = np.minimum((s.NP[s.NP[:, 3] == side, 2] * pfc2.WHL_RATE).astype(int), len(w) - 1)
    ends[side] = (np.median(px[f]), np.median(py[f]))
print(f"arm ends (median position at the end of the table's correct trials): right {np.round(ends[1], 0)}, left "
      f"{np.round(ends[2], 0)} pixels; sensor pokes: right channel {len(PK[1])}, left channel {len(PK[2])}")


def all_trials(radius):
    ev = []
    for sd in (1, 2):
        near = ok_f & (np.hypot(px - ends[sd][0], py - ends[sd][1]) < radius)
        ev += [(t_f[i], "arm", sd, t_f[i]) for i in np.flatnonzero(np.diff(near.astype(int), prepend=0) == 1)]
        ev += [(a, "poke", sd, b) for a, b in PK[sd]]
    ev.sort()
    out, pending = [], []
    for t, kind, sd, t_end in ev:
        if kind == "poke":
            pending.append((t, sd, t_end))
        elif pending:
            # first poke, channels poked, arm reached, arrival, end of the last poke, number of pokes
            out.append((pending[0][0], {c for _, c, _ in pending}, sd, t, pending[-1][2], len(pending)))
            pending = []
    return out


PERF = {}
for radius in (60, 40, 80):
    tr = all_trials(radius)
    mixed = sum(len(x[1]) > 1 for x in tr)
    corr = [x for x in tr if len(x[1]) == 1 and min(x[1]) == x[2]]
    err = [x for x in tr if len(x[1]) == 1 and min(x[1]) != x[2]]
    matched = sum(any(abs(x[0] - r[0]) < 0.05 and x in corr and min(x[1]) == r[3] for x in tr) for r in s.NP)
    PERF[radius] = 100 * len(corr) / (len(corr) + len(err))
    print(f"radius {radius} px: {len(tr)} trials, {len(corr)} correct, {len(err)} error, {mixed} with pokes on both channels; "
          f"table rows matched to a detected correct trial of the same side by the first poke (within 50 ms): {matched} of "
          f"{len(s.NP)}; performance {PERF[radius]:.1f}% (within 91.9 +- 9.19: {abs(PERF[radius] - 91.9) <= 9.19})")
    if radius == 60:
        print(f"    table trials, first poke to trial end: {(s.NP[:, 2] - s.NP[:, 0]).min():.1f} to {(s.NP[:, 2] - s.NP[:, 0]).max():.1f} s")
        for x in tr:
            if not any(abs(x[0] - r[0]) < 0.05 for r in s.NP):
                print(f"    not in the table: first poke {x[0]:.1f} s, channel {sorted(x[1])}, {x[5]} pokes to {x[4]:.1f} s, arm {x[2]}, "
                      f"arrival {x[3]:.1f} s ({x[3] - x[0]:.1f} s after the first poke, {x[3] - x[4]:.1f} s after the last; "
                      f"{'correct' if x in corr else 'error'})")
        lmax = (s.NP[:, 2] - s.NP[:, 0]).max()
        nc = sum(x[3] - x[0] <= lmax for x in corr)
        ne = sum(x[3] - x[0] <= lmax for x in err)
        print(f"    without the detected trials whose arm arrival comes later after the first poke than the table's longest "
              f"trial ({lmax:.1f} s): {nc} correct, {ne} error, {100 * nc / (nc + ne):.1f}%")
        for r in s.NP:
            if not any(abs(x[0] - r[0]) < 0.05 and x in corr and min(x[1]) == r[3] for x in tr):
                inn = [x for x in tr if x[0] <= r[0] <= x[3]]
                print(f"    table row not matched: first poke {r[0]:.1f} s, side {int(r[3])}, trial end {r[2]:.1f} s; "
                      + "; ".join(f"detected trial from {x[0]:.1f} s ({r[0] - x[0]:.1f} s earlier), channels {sorted(x[1])}, "
                                  f"{x[5]} pokes, arm {x[2]} at {x[3]:.1f} s ({'correct' if x in corr else 'error' if x in err else 'both channels'})"
                                  for x in inn))
print(f"with the table's correct trials and the detected errors (radius 60 px): "
      f"{100 * len(s.NP) / (len(s.NP) + sum(1 for x in all_trials(60) if len(x[1]) == 1 and min(x[1]) != x[2])):.1f}%")

EE.188.tar.gz: 75 members, for example EE.188.xml (0 MB), EE.188.whlrl (3 MB), EE.188.eeg (731 MB), EE.188.whl (3 MB), EE.188.res (12 MB), EE188_example.mat (3 MB)
member types: ['RUN', 'clu', 'eeg', 'fet', 'led', 'mat', 'nrs', 'res', 'spk', 'threshold', 'whl', 'whlrl', 'xml']
arm ends (median position at the end of the table's correct trials): right [526. 134.], left [525. 385.] pixels; sensor pokes: right channel 78, left channel 111
radius 60 px: 46 trials, 45 correct, 1 error, 0 with pokes on both channels; table rows matched to a detected correct trial of the same side by the first poke (within 50 ms): 41 of 42; performance 97.8% (within 91.9 +- 9.19: True)
    table trials, first poke to trial end: 9.0 to 36.1 s
    not in the table: first poke 822.9 s, channel [2], 6 pokes to 829.3 s, arm 2, arrival 1022.7 s (199.8 s after the first poke, 193.4 s after the last; correct)
    not in the table: first poke 1130.0 s, channel [2], 5 pokes to 1135.6 s, arm 2, arrival 1141.2 s (11.3 s after the first poke, 5.7 s after the last; correct)
    not in the table: first poke 2245.1 s, channel [2], 4 pokes to 2251.4 s, arm 2, arrival 2270.9 s (25.9 s after the first poke, 19.6 s after the last; correct)
    not in the table: first poke 2370.7 s, channel [2], 8 pokes to 2453.7 s, arm 2, arrival 2630.6 s (259.9 s after the first poke, 176.9 s after the last; correct)
    not in the table: first poke 2687.2 s, channel [2], 5 pokes to 2693.1 s, arm 1, arrival 2699.2 s (12.0 s after the first poke, 6.1 s after the last; error)
    without the detected trials whose arm arrival comes later after the first poke than the table's longest trial (36.1 s): 43 correct, 1 error, 97.7%
    table row not matched: first poke 1133.5 s, side 2, trial end 1145.1 s; detected trial from 1130.0 s (3.5 s earlier), channels [2], 5 pokes, arm 2 at 1141.2 s (correct)
radius 40 px: 46 trials, 45 correct, 1 error, 0 with pokes on both channels; table rows matched to a detected correct trial of the same side by the first poke (within 50 ms): 41 of 42; performance 97.8% (within 91.9 +- 9.19: True)
radius 80 px: 48 trials, 46 correct, 2 error, 0 with pokes on both channels; table rows matched to a detected correct trial of the same side by the first poke (within 50 ms): 40 of 42; performance 95.8% (within 91.9 +- 9.19: True)
with the table's correct trials and the detected errors (radius 60 px): 97.7%

def session_psths(s):
    return np.array([pfc2.smooth(s.trial_hist(i)) for i in range(s.n)])


H = {k: session_psths(s) for k, s in S.items()}
LR = {}
rng = np.random.default_rng(101)
t0 = time.time()
for k, s in S.items():
    LR[k] = [pfc2.lr_test(H[k][i], s.left, rng) for i in range(s.n)]
print("left-right permutation tests:", sum(len(v) for v in LR.values()), "cells in", round(time.time() - t0), "s")
G = pfc2.GRID


def selective(res, lo, hi, reading):
    """(a) any part of a significant segment in [lo, hi]; (b) the global band broken in [lo, hi]."""
    D0, point, glob, gp = res
    inr = (G >= lo) & (G <= hi)
    if reading == "a":
        return any(inr[a:b + 1].any() for a, b in pfc2.segments(point, glob))
    return bool((glob & inr).any())


# null check: the same test with the trial labels shuffled once per cell
rng0 = np.random.default_rng(102)
null = [pfc2.lr_test(H[k][i], rng0.permutation(s.left), rng0)[3] <= 0.05 for k, s in S.items() for i in range(s.n)]
print(f"false positive rate with shuffled labels (global, anywhere): {np.mean(null):.3f} of {len(null)} cells (level 0.05)")

left-right permutation tests: 627 cells in 22 s

false positive rate with shuffled labels (global, anywhere): 0.051 of 627 cells (level 0.05)

def sync_windows(s, width=0.1):
    fr, pos = [], []
    for t0, t1 in s.NP[:, 1:3]:
        edges = np.arange(t0, t1 - width + 1e-9, width)
        for a in edges:
            i, j = np.searchsorted(s.t, [a, a + width])
            fr.append(len(np.unique(s.u[i:j])) / s.n)
            f = min(int((a + width / 2) * pfc2.WHL_RATE), len(s.whl) - 1)
            pos.append(s.whl[f, 6] - 1 if s.whl[f, 6] > 0 else np.nan)
    return np.array(fr), np.array(pos)


SYNC = {k: sync_windows(s) for k, s in S.items()}


def sync_summary(keys):
    fr = np.concatenate([SYNC[k][0] for k in keys])
    pos = np.concatenate([SYNC[k][1] for k in keys])
    ok = ~np.isnan(pos)
    tenth = [fr[ok & (pos >= a / 10) & (pos < (a + 1) / 10 + (a == 9))].mean() for a in range(10)]
    return 100 * fr.mean(), 100 * np.array(tenth), len(fr), np.mean(~ok)


for k in S:
    m, tn, nw, miss = sync_summary([k])
    print(f"{k} ({S[k].layer}, {S[k].n} cells): {m:5.1f}% over {nw} windows; per tenth {np.round(tn, 1)}")
print()
for label, keys in [("(a) EE188, layer 2/3 Figure 2a session", ["EE188"]), ("(b) all layer 2/3 sessions", L23),
                    ("(b) all layer 5 sessions", L5)]:
    m, tn, nw, miss = sync_summary(keys)
    within = np.all(np.abs(tn - m) <= 0.5 * m)
    print(f"{label}: mean {m:.1f}% of cells active per 100 ms window ({nw} windows, {100 * miss:.1f}% without a position); "
          f"tenths {tn.min():.1f} to {tn.max():.1f}%, all within 50% of the mean: {within}")

EE049 (L2/3, 44 cells):  15.1% over 2334 windows; per tenth [14.8 13.1 13.8 14.5 15.7 17.  17.5 16.  16.2 13.9]
EE088 (L2/3, 28 cells):  16.2% over 1669 windows; per tenth [17.4 15.8 15.4 14.5 16.5 17.  18.  17.1 14.1 16.4]
EE152 (L2/3, 61 cells):  11.9% over 2847 windows; per tenth [12.6 10.4 10.7 10.9 11.3 11.3 13.2 12.3 12.2 11.4]
EE165 (L2/3, 70 cells):  13.0% over 3393 windows; per tenth [16.3 13.8 10.4 11.5 13.3 11.8 13.2 12.7 13.  12.8]
EE188 (L2/3, 117 cells):  11.8% over 3575 windows; per tenth [12.8 11.6 10.9 12.2 14.1 12.6 12.5 12.4 12.8 10.2]
EE196 (L2/3, 66 cells):  19.1% over 2293 windows; per tenth [22.3 20.5 13.6 15.4 17.2 19.6 19.7 17.9 19.8 20.7]
EE210 (L2/3, 60 cells):  13.7% over 2368 windows; per tenth [16.9 13.6 12.6 11.4 13.5 13.4 13.4 13.5 14.4 13.8]
FF103 (L5, 42 cells):  17.9% over 2636 windows; per tenth [18.3 18.  17.6 16.7 16.8 17.5 18.6 18.6 18.1 17.5]
FF130 (L5, 64 cells):  21.9% over 2928 windows; per tenth [24.4 21.4 21.6 22.4 22.  21.1 20.3 20.1 19.8 20.3]
FF167 (L5, 40 cells):  21.0% over 2630 windows; per tenth [19.1 21.2 21.3 21.4 22.7 23.2 21.  22.  21.6 19.2]
GG069 (L2/3, 14 cells):   9.7% over 3761 windows; per tenth [ 8.5  7.6  8.8 10.  11.8 11.8 12.  11.  12.6  9.7]
GG097 (L2/3, 12 cells):   7.6% over 2284 windows; per tenth [6.6 9.  6.2 7.2 8.1 8.8 9.8 8.6 7.2 8.6]
GG149 (L2/3, 9 cells):  10.5% over 1837 windows; per tenth [ 9.7  6.8  7.4  6.6  8.2 12.6 15.1 16.2 14.2 12.2]

(a) EE188, layer 2/3 Figure 2a session: mean 11.8% of cells active per 100 ms window (3575 windows, 4.8% without a position); tenths 10.2 to 14.1%, all within 50% of the mean: True
(b) all layer 2/3 sessions: mean 12.6% of cells active per 100 ms window (26361 windows, 4.7% without a position); tenths 11.2 to 14.5%, all within 50% of the mean: True
(b) all layer 5 sessions: mean 20.3% of cells active per 100 ms window (8194 windows, 1.0% without a position); tenths 19.0 to 20.8%, all within 50% of the mean: True

def uniformity(k):
    s = S[k]
    PL = H[k][:, s.left].mean(axis=1)
    PR = H[k][:, ~s.left].mean(axis=1)
    peak = np.where(PL.max(1) >= PR.max(1), G[PL.argmax(1)], G[PR.argmax(1)])
    counts = np.histogram(peak, bins=np.linspace(0, 1, 11))[0]
    share = counts / (s.n / 10)
    chi = stats.chisquare(counts)
    pop = H[k].mean(axis=(0, 1))
    tenth = np.array([pop[(G >= a / 10) & (G < (a + 1) / 10)].mean() for a in range(10)])
    rel = tenth / pop.mean()
    ok = (share.max() <= 2) & (share.min() >= 0.5) & np.all(np.abs(rel - 1) <= 0.5)
    print(f"{k} ({s.layer}, {s.n} cells): peaks per tenth {counts.tolist()}, share of expected "
          f"{share.min():.2f} to {share.max():.2f}; chi-square against uniform {chi.statistic:.1f}, p = {chi.pvalue:.3f}; "
          f"population mean rate per tenth / overall mean {rel.min():.2f} to {rel.max():.2f}; criterion met: {ok}")
    return peak, PL, PR


peak188, PL188, PR188 = uniformity("EE188")
for k in L5:
    uniformity(k)
# time spent and speed per tenth of the maze in EE188 (LED midpoint, LED 1 alone where LED 2 is
# not tracked), as the PSTHs count spikes per unit of normalised length per trial
s = S["EE188"]
w = s.whl
occ, spd = np.zeros(10), [[] for _ in range(10)]
for t0_, t1_ in s.NP[:, 1:3]:
    f = np.arange(int(t0_ * pfc2.WHL_RATE), min(int(t1_ * pfc2.WHL_RATE) + 1, len(w)))
    f = f[(w[f, 6] > 0) & (w[f, 0] > 0)]
    tb = np.minimum(((w[f, 6] - 1) * 10).astype(int), 9)
    occ += np.bincount(tb, minlength=10) / pfc2.WHL_RATE
    two = w[f, 2] > 0
    mx = np.where(two, (w[f, 0] + w[f, 2]) / 2, w[f, 0])
    my = np.where(two, (w[f, 1] + w[f, 3]) / 2, w[f, 1])
    v = np.hypot(np.diff(mx), np.diff(my)) * pfc2.WHL_RATE
    for a in range(10):
        spd[a] += list(v[(tb[1:] == a) & (np.diff(f) == 1)])
print("EE188 time per trial in each tenth of the maze (s): " + ", ".join(f"{x:.2f}" for x in occ / len(s.NP)))
print("EE188 median speed in each tenth (pixels/s): " + ", ".join(f"{np.median(x):.0f}" for x in spd))
fig, ax = plt.subplots(1, 2, figsize=(9, 4))
o = np.argsort(peak188)
nrm = np.maximum(PL188.max(1), PR188.max(1))[:, None]
ax[0].imshow(PL188[o] / nrm, aspect="auto", extent=[0, 1, len(o), 0], vmin=0, vmax=1)
ax[1].imshow(PR188[o] / nrm, aspect="auto", extent=[0, 1, len(o), 0], vmin=0, vmax=1)
ax[0].set_title("EE188 left, rate / cell's peak (0 to 1)")
ax[1].set_title("EE188 right")
for a in ax:
    a.set_xlabel("normalised position")
plt.show()

EE188 (L2/3, 117 cells): peaks per tenth [17, 0, 2, 2, 16, 7, 7, 5, 0, 61], share of expected 0.00 to 5.21; chi-square against uniform 258.8, p = 0.000; population mean rate per tenth / overall mean 0.57 to 2.33; criterion met: False
FF103 (L5, 42 cells): peaks per tenth [31, 0, 1, 0, 2, 2, 3, 0, 1, 2], share of expected 0.00 to 7.38; chi-square against uniform 192.3, p = 0.000; population mean rate per tenth / overall mean 0.58 to 2.11; criterion met: False
FF130 (L5, 64 cells): peaks per tenth [44, 0, 0, 1, 7, 7, 1, 0, 2, 2], share of expected 0.00 to 6.88; chi-square against uniform 255.4, p = 0.000; population mean rate per tenth / overall mean 0.51 to 2.33; criterion met: False
FF167 (L5, 40 cells): peaks per tenth [19, 0, 1, 1, 8, 7, 1, 0, 1, 2], share of expected 0.00 to 4.75; chi-square against uniform 80.5, p = 0.000; population mean rate per tenth / overall mean 0.51 to 1.78; criterion met: False
EE188 time per trial in each tenth of the maze (s): 0.86, 0.49, 0.47, 0.50, 0.71, 0.67, 0.57, 0.51, 0.47, 2.89
EE188 median speed in each tenth (pixels/s): 81, 174, 194, 181, 138, 155, 192, 198, 188, 22

def widths(P, how):
    out = []
    for p in P:
        above = p >= 0.5 * p.max()
        if how == "ii":
            out.append(above.mean())
        else:
            k = p.argmax()
            a = b = k
            while a > 0 and above[a - 1]:
                a -= 1
            while b < len(p) - 1 and above[b + 1]:
                b += 1
            out.append((b - a + 1) / len(p))
    return np.array(out)


def psth_sets(keys, which):
    P = []
    for k in keys:
        s = S[k]
        PL = H[k][:, s.left].mean(1)
        PR = H[k][:, ~s.left].mean(1)
        if which == "a":
            P.append(H[k].mean(1))
        elif which == "b":
            P += [PL, PR]
        else:
            P.append(np.where((PL.max(1) >= PR.max(1))[:, None], PL, PR))
    P = np.concatenate(P)
    return P[P.max(1) > 0]


WIDTH = {}
for cells, keys in [("x", ["EE188"]), ("y", list(S))]:
    for which in "abc":
        for how in ("i", "ii"):
            w = widths(psth_sets(keys, which), how)
            WIDTH[(which, how, cells)] = (w.mean(), w.std(ddof=1))
            ok = abs(w.mean() - 0.27) <= 0.17 and abs(w.std(ddof=1) - 0.17) <= 0.017
            print(f"PSTH {which}, width {how:2s}, cells {cells}: n = {len(w):4d}, mean {w.mean():.3f}, s.d. "
                  f"{w.std(ddof=1):.3f}, range {w.min():.3f} to {w.max():.3f}; both within tolerance: {ok}")
w = widths(psth_sets(list(S), "c"), "i")
for r in ("EE", "FF", "GG"):
    wr = widths(psth_sets([k for k in S if S[k].rat == r], "c"), "i")
    print(f"rat {r}, reading c/i: mean {wr.mean():.3f}, s.d. {wr.std(ddof=1):.3f}, n = {len(wr)}")

PSTH a, width i , cells x: n =  117, mean 0.169, s.d. 0.125, range 0.080 to 0.885; both within tolerance: False
PSTH a, width ii, cells x: n =  117, mean 0.235, s.d. 0.158, range 0.085 to 0.885; both within tolerance: True
PSTH b, width i , cells x: n =  234, mean 0.166, s.d. 0.108, range 0.070 to 0.805; both within tolerance: False
PSTH b, width ii, cells x: n =  234, mean 0.228, s.d. 0.142, range 0.075 to 0.805; both within tolerance: False
PSTH c, width i , cells x: n =  117, mean 0.165, s.d. 0.117, range 0.075 to 0.805; both within tolerance: False
PSTH c, width ii, cells x: n =  117, mean 0.209, s.d. 0.138, range 0.075 to 0.805; both within tolerance: False
PSTH a, width i , cells y: n =  627, mean 0.209, s.d. 0.154, range 0.080 to 0.990; both within tolerance: True
PSTH a, width ii, cells y: n =  627, mean 0.286, s.d. 0.188, range 0.080 to 0.990; both within tolerance: False
PSTH b, width i , cells y: n = 1254, mean 0.199, s.d. 0.138, range 0.070 to 0.995; both within tolerance: False
PSTH b, width ii, cells y: n = 1254, mean 0.267, s.d. 0.171, range 0.075 to 0.995; both within tolerance: True
PSTH c, width i , cells y: n =  627, mean 0.193, s.d. 0.129, range 0.075 to 0.950; both within tolerance: False
PSTH c, width ii, cells y: n =  627, mean 0.239, s.d. 0.156, range 0.075 to 0.950; both within tolerance: True
rat EE, reading c/i: mean 0.195, s.d. 0.134, n = 446
rat FF, reading c/i: mean 0.198, s.d. 0.121, n = 146
rat GG, reading c/i: mean 0.143, s.d. 0.082, n = 35

NB = 100
GB = (np.arange(NB) + 0.5) / NB


def trial_paths(s):
    w = s.whl
    out = []
    for r, (t0, t1) in enumerate(s.NP[:, 1:3]):
        f = np.arange(int(t0 * pfc2.WHL_RATE), min(int(t1 * pfc2.WHL_RATE) + 1, len(w)))
        f = f[(w[f, 6] > 0) & (w[f, 0] > 0)]
        d = w[f, 6] - 1
        b = np.minimum((d * NB).astype(int), NB - 1)
        two = w[f, 2] > 0
        # LED midpoint where both LEDs are tracked, LED 1 alone otherwise
        mx = np.where(two, (w[f, 0] + w[f, 2]) / 2, w[f, 0])
        my = np.where(two, (w[f, 1] + w[f, 3]) / 2, w[f, 1])
        hd = np.where(two, np.arctan2(w[f, 3] - w[f, 1], w[f, 2] - w[f, 0]), np.nan)
        vals = []
        for v in (mx, my, np.cos(hd), np.sin(hd)):
            okv = ~np.isnan(v)
            sm = np.bincount(b[okv], v[okv], NB)
            n = np.bincount(b[okv], minlength=NB)
            m = np.where(n > 0, sm / np.maximum(n, 1), np.nan)
            ok = ~np.isnan(m)
            m = np.interp(np.arange(NB), np.flatnonzero(ok), m[ok]) if ok.sum() > 1 else np.full(NB, np.nan)
            vals.append(m)
        out.append(vals)
    return np.array(out)  # trials x 4 x NB


def perm_stat(X, labels, rng, M, stat, alpha):
    labels = np.asarray(labels, bool)
    D = np.empty((M + 1, X.shape[-1]))
    for j in range(M + 1):
        lab = labels if j == 0 else rng.permutation(labels)
        D[j] = stat(X, lab)
    point, glob, gp = pfc2.global_test(D, alpha, "two")
    return D[0], point, glob, gp


def stat_x(X, lab):
    return X[lab, 0].mean(0) - X[~lab, 0].mean(0)


def stat_y(X, lab):
    return X[lab, 1].mean(0) - X[~lab, 1].mean(0)


def stat_dist(X, lab):
    mx, my = X[:, 0].mean(0), X[:, 1].mean(0)
    dist = np.hypot(X[:, 0] - mx, X[:, 1] - my)
    return dist[lab].mean(0) - dist[~lab].mean(0)


def stat_hd(X, lab):
    aL = np.arctan2(X[lab, 3].mean(0), X[lab, 2].mean(0))
    aR = np.arctan2(X[~lab, 3].mean(0), X[~lab, 2].mean(0))
    return np.angle(np.exp(1j * (aL - aR)))


def first_sig(res):
    D0, point, glob, gp = res
    if not glob.any():
        return np.nan, np.nan
    g = np.flatnonzero(glob)[0]
    a = g
    while a > 0 and point[a - 1]:
        a -= 1
    return GB[g], GB[a]


def smooth100(X):
    from scipy.ndimage import gaussian_filter1d
    return gaussian_filter1d(X, 0.05 * NB, axis=-1, mode="nearest")


rng = np.random.default_rng(103)
TRAJ = {}
for k in ["EE188"] + [k for k in S if k != "EE188"]:
    s = S[k]
    X = trial_paths(s)
    good = ~np.isnan(X[:, :2]).any(axis=(1, 2))
    X, lab = X[good], s.left[good]
    hd_ok = not np.isnan(X[:, 2:]).any()
    row = []
    for name, stat in [("path a (lateral, y)", stat_y), ("path b (distance to mean path)", stat_dist),
                       ("path, x (along the stem)", stat_x), ("head direction", stat_hd)]:
        if name == "head direction" and not hd_ok:
            row.append("head direction: LED 2 not tracked")
            continue
        for sm in (False, True):
            r = perm_stat(smooth100(X) if sm else X, lab, rng, 1000, stat, 0.01)
            g, a = first_sig(r)
            TRAJ[(k, name, sm)] = (g, a)
            row.append(f"{name}{' smoothed' if sm else ''}: first global {g:.2f} (segment from {a:.2f})")
    print(f"{k} ({good.sum()} of {len(good)} trials with positions): " + "; ".join(row))
    if k == "EE188":
        print()
for name in ("path a (lateral, y)", "path b (distance to mean path)", "path, x (along the stem)"):
    g = TRAJ[("EE188", name, False)][0]
    h = TRAJ[("EE188", "head direction", False)][0]
    print(f"EE188 {name}: path {g:.2f} (0.30 +- 0.03: {abs(g - 0.30) <= 0.03}), head direction {h:.2f} "
          f"(0.24 +- 0.024: {abs(h - 0.24) <= 0.024})")
# Maze geometry in EE188: the mean path of left and right trials (LED midpoint, pixels), to show
# the orientation of the maze and whether both trial types start at the same place.
X = trial_paths(S["EE188"])
lab = S["EE188"].left
good = ~np.isnan(X[:, :2]).any(axis=(1, 2))
fig, ax = plt.subplots(figsize=(4, 4))
for m, name in ((lab & good, "left"), (~lab & good, "right")):
    mx, my = np.nanmean(X[m, 0], 0), np.nanmean(X[m, 1], 0)
    ax.plot(mx, my, label=name)
    ax.plot(mx[0], my[0], "o", color="k")
    pts = {p: (mx[int(p * NB) - (1 if p == 1 else 0)], my[int(p * NB) - (1 if p == 1 else 0)]) for p in (0, 0.3, 0.5, 1)}
    print(f"EE188 {name} trials, mean path (x, y) at positions 0, 0.3, 0.5, 1: "
          + "; ".join(f"{p}: ({v[0]:.0f}, {v[1]:.0f})" for p, v in pts.items())
          + f"; start position s.d. over trials x {np.nanstd(X[m, 0, 0]):.0f}, y {np.nanstd(X[m, 1, 0]):.0f}")
ax.set_xlabel("x (pixels)")
ax.set_ylabel("y (pixels)")
ax.legend()
ax.set_title("EE188 mean paths, dots at position 0")
plt.show()
st = np.hypot(np.nanmean(X[lab & good, 0, 0]) - np.nanmean(X[~lab & good, 0, 0]),
              np.nanmean(X[lab & good, 1, 0]) - np.nanmean(X[~lab & good, 1, 0]))
# head direction near the start: mean angle (LED 1 to LED 2, degrees) of left and right trials,
# their circular difference, and the spread over trials (circular s.d.)
for b in (0, 1, 2, 5, 10, 20, 30):
    row = []
    for m in (lab & good, ~lab & good):
        c, sn = np.nanmean(X[m, 2, b]), np.nanmean(X[m, 3, b])
        row.append((np.degrees(np.arctan2(sn, c)), np.degrees(np.sqrt(-2 * np.log(np.hypot(c, sn))))))
    dif = (row[0][0] - row[1][0] + 180) % 360 - 180
    print(f"EE188 head direction at position {GB[b]:.3f}: left {row[0][0]:.0f} (s.d. {row[0][1]:.0f}), right {row[1][0]:.0f} "
          f"(s.d. {row[1][1]:.0f}), difference {dif:.0f} degrees")
print(f"distance between the mean start positions of left and right trials: {st:.0f} pixels; plotted ranges x "
      f"{np.nanmin(X[good, 0]):.0f} to {np.nanmax(X[good, 0]):.0f}, y {np.nanmin(X[good, 1]):.0f} to {np.nanmax(X[good, 1]):.0f}")

EE188 (42 of 42 trials with positions): path a (lateral, y): first global 0.28 (segment from 0.27); path a (lateral, y) smoothed: first global 0.28 (segment from 0.27); path b (distance to mean path): first global 0.43 (segment from 0.42); path b (distance to mean path) smoothed: first global 0.40 (segment from 0.40); path, x (along the stem): first global 0.48 (segment from 0.47); path, x (along the stem) smoothed: first global 0.41 (segment from 0.38); head direction: first global 0.01 (segment from 0.01); head direction smoothed: first global 0.01 (segment from 0.01)

EE049 (19 of 19 trials with positions): path a (lateral, y): first global 0.43 (segment from 0.42); path a (lateral, y) smoothed: first global 0.35 (segment from 0.32); path b (distance to mean path): first global 0.52 (segment from 0.51); path b (distance to mean path) smoothed: first global 0.49 (segment from 0.48); path, x (along the stem): first global 0.10 (segment from 0.09); path, x (along the stem) smoothed: first global 0.10 (segment from 0.09); head direction: LED 2 not tracked

EE088 (31 of 31 trials with positions): path a (lateral, y): first global 0.35 (segment from 0.33); path a (lateral, y) smoothed: first global 0.30 (segment from 0.29); path b (distance to mean path): first global 0.43 (segment from 0.40); path b (distance to mean path) smoothed: first global 0.39 (segment from 0.34); path, x (along the stem): first global 0.41 (segment from 0.41); path, x (along the stem) smoothed: first global 0.47 (segment from 0.45); head direction: LED 2 not tracked

EE152 (36 of 36 trials with positions): path a (lateral, y): first global 0.17 (segment from 0.12); path a (lateral, y) smoothed: first global 0.12 (segment from 0.01); path b (distance to mean path): first global 0.46 (segment from 0.45); path b (distance to mean path) smoothed: first global 0.42 (segment from 0.41); path, x (along the stem): first global 0.04 (segment from 0.04); path, x (along the stem) smoothed: first global 0.07 (segment from 0.01); head direction: first global 0.17 (segment from 0.15); head direction smoothed: first global 0.17 (segment from 0.15)

EE165 (42 of 42 trials with positions): path a (lateral, y): first global 0.26 (segment from 0.24); path a (lateral, y) smoothed: first global 0.10 (segment from 0.07); path b (distance to mean path): first global 0.43 (segment from 0.42); path b (distance to mean path) smoothed: first global 0.41 (segment from 0.41); path, x (along the stem): first global 0.52 (segment from 0.51); path, x (along the stem) smoothed: first global 0.48 (segment from 0.47); head direction: first global 0.26 (segment from 0.23); head direction smoothed: first global 0.27 (segment from 0.24)

EE196 (48 of 48 trials with positions): path a (lateral, y): first global 0.29 (segment from 0.21); path a (lateral, y) smoothed: first global 0.24 (segment from 0.21); path b (distance to mean path): first global 0.48 (segment from 0.47); path b (distance to mean path) smoothed: first global 0.45 (segment from 0.45); path, x (along the stem): first global 0.01 (segment from 0.01); path, x (along the stem) smoothed: first global 0.01 (segment from 0.01); head direction: LED 2 not tracked

EE210 (53 of 53 trials with positions): path a (lateral, y): first global 0.24 (segment from 0.23); path a (lateral, y) smoothed: first global 0.24 (segment from 0.21); path b (distance to mean path): first global 0.12 (segment from 0.12); path b (distance to mean path) smoothed: first global 0.15 (segment from 0.14); path, x (along the stem): first global 0.01 (segment from 0.01); path, x (along the stem) smoothed: first global 0.01 (segment from 0.01); head direction: LED 2 not tracked

FF103 (39 of 39 trials with positions): path a (lateral, y): first global 0.38 (segment from 0.35); path a (lateral, y) smoothed: first global 0.34 (segment from 0.33); path b (distance to mean path): first global 0.49 (segment from 0.49); path b (distance to mean path) smoothed: first global 0.47 (segment from 0.46); path, x (along the stem): first global 0.41 (segment from 0.41); path, x (along the stem) smoothed: first global 0.35 (segment from 0.34); head direction: LED 2 not tracked

FF130 (46 of 46 trials with positions): path a (lateral, y): first global 0.38 (segment from 0.36); path a (lateral, y) smoothed: first global 0.34 (segment from 0.34); path b (distance to mean path): first global 0.04 (segment from 0.03); path b (distance to mean path) smoothed: first global nan (segment from nan); path, x (along the stem): first global 0.01 (segment from 0.01); path, x (along the stem) smoothed: first global 0.01 (segment from 0.01); head direction: LED 2 not tracked

FF167 (43 of 43 trials with positions): path a (lateral, y): first global 0.29 (segment from 0.28); path a (lateral, y) smoothed: first global 0.27 (segment from 0.24); path b (distance to mean path): first global 0.51 (segment from 0.47); path b (distance to mean path) smoothed: first global 0.47 (segment from 0.46); path, x (along the stem): first global 0.42 (segment from 0.41); path, x (along the stem) smoothed: first global 0.42 (segment from 0.41); head direction: LED 2 not tracked

def layer_fraction(keys, lo, hi, reading):
    v = [selective(LR[k][i], lo, hi, reading) for k in keys for i in range(S[k].n)]
    return 100 * np.mean(v), int(np.sum(v)), len(v)


for reading in "ab":
    f5 = layer_fraction(L5, 0.5, 1.0, reading)
    f23 = layer_fraction(L23, 0.5, 1.0, reading)
    print(f"reading {reading}: layer 5 {f5[0]:.1f}% ({f5[1]}/{f5[2]}), within 40 +- 4: {abs(f5[0] - 40) <= 4}; "
          f"layer 2/3 {f23[0]:.1f}% ({f23[1]}/{f23[2]}), within 22 +- 2.2: {abs(f23[0] - 22) <= 2.2}")
print("per session, reading a:")
for k in S:
    f = layer_fraction([k], 0.5, 1.0, "a")
    print(f"  {k} ({S[k].layer}): {f[0]:.1f}% ({f[1]}/{f[2]})")
print("per rat, reading a:", {r: round(layer_fraction([k for k in S if S[k].rat == r], 0.5, 1.0, 'a')[0], 1) for r in ("EE", "FF", "GG")})

reading a: layer 5 62.3% (91/146), within 40 +- 4: False; layer 2/3 54.3% (261/481), within 22 +- 2.2: False
reading b: layer 5 61.6% (90/146), within 40 +- 4: False; layer 2/3 51.8% (249/481), within 22 +- 2.2: False
per session, reading a:
  EE049 (L2/3): 34.1% (15/44)
  EE088 (L2/3): 46.4% (13/28)
  EE152 (L2/3): 63.9% (39/61)
  EE165 (L2/3): 44.3% (31/70)
  EE188 (L2/3): 50.4% (59/117)
  EE196 (L2/3): 74.2% (49/66)
  EE210 (L2/3): 66.7% (40/60)
  FF103 (L5): 31.0% (13/42)
  FF130 (L5): 78.1% (50/64)
  FF167 (L5): 70.0% (28/40)
  GG069 (L2/3): 28.6% (4/14)
  GG097 (L2/3): 41.7% (5/12)
  GG149 (L2/3): 66.7% (6/9)
per rat, reading a: {'EE': np.float64(55.2), 'FF': np.float64(62.3), 'GG': np.float64(42.9)}

for reading in "ab":
    f23 = layer_fraction(L23, 0.0, 0.3, reading)
    f5 = layer_fraction(L5, 0.0, 0.3, reading)
    print(f"reading {reading}: layer 2/3 {f23[0]:.1f}% ({f23[1]}/{f23[2]}), within 16 +- 1.6: {abs(f23[0] - 16) <= 1.6}; "
          f"layer 5 {f5[0]:.1f}% ({f5[1]}/{f5[2]}), within 18 +- 1.8: {abs(f5[0] - 18) <= 1.8}")
for k in S:
    f = layer_fraction([k], 0.0, 0.3, "a")
    print(f"  {k} ({S[k].layer}): {f[0]:.1f}% ({f[1]}/{f[2]})")
print("per rat, reading a:", {r: round(layer_fraction([k for k in S if S[k].rat == r], 0.0, 0.3, 'a')[0], 1) for r in ("EE", "FF", "GG")})

reading a: layer 2/3 31.4% (151/481), within 16 +- 1.6: False; layer 5 26.7% (39/146), within 18 +- 1.8: False
reading b: layer 2/3 28.7% (138/481), within 16 +- 1.6: False; layer 5 22.6% (33/146), within 18 +- 1.8: False
  EE049 (L2/3): 15.9% (7/44)
  EE088 (L2/3): 35.7% (10/28)
  EE152 (L2/3): 16.4% (10/61)
  EE165 (L2/3): 18.6% (13/70)
  EE188 (L2/3): 45.3% (53/117)
  EE196 (L2/3): 43.9% (29/66)
  EE210 (L2/3): 31.7% (19/60)
  FF103 (L5): 9.5% (4/42)
  FF130 (L5): 26.6% (17/64)
  FF167 (L5): 45.0% (18/40)
  GG069 (L2/3): 57.1% (8/14)
  GG097 (L2/3): 8.3% (1/12)
  GG149 (L2/3): 11.1% (1/9)
per rat, reading a: {'EE': np.float64(31.6), 'FF': np.float64(26.7), 'GG': np.float64(28.6)}

def nosepoke_rates(s):
    R = np.zeros((s.ntrials, s.n))
    for r, (a, b) in enumerate(s.NP[:, :2]):
        i, j = np.searchsorted(s.t, [a, b])
        R[r] = np.bincount(s.u[i:j], minlength=s.n) / (b - a)
    return R


ODOR = {}
for k, s in S.items():
    R = nosepoke_rates(s)
    for reading, eq in (("a", True), ("b", False)):
        with np.errstate(invalid="ignore", divide="ignore"):
            p = stats.ttest_ind(R[s.left], R[~s.left], equal_var=eq).pvalue
        ODOR[(k, reading)] = (p < 0.05, R[s.left].mean(0) - R[~s.left].mean(0), p)
    if k == "EE188":
        for reading in "ab":
            sig = ODOR[(k, reading)][0]
            print(f"EE188 reading {reading}: {100 * sig.mean():.1f}% ({sig.sum()}/{s.n}) odor selective; "
                  f"within 26 +- 2.6: {abs(100 * sig.mean() - 26) <= 2.6}")
        print(f"nose-poke durations EE188: median {np.median(s.NP[:, 1] - s.NP[:, 0]):.1f} s, range "
              f"{np.min(s.NP[:, 1] - s.NP[:, 0]):.1f} to {np.max(s.NP[:, 1] - s.NP[:, 0]):.1f} s")
print("other sessions, reading a:", {k: f"{100 * ODOR[(k, 'a')][0].mean():.0f}%" for k in S if k != "EE188"})
allsig = np.concatenate([ODOR[(k, "a")][0] for k in S])
print(f"all sessions, reading a: {100 * allsig.mean():.1f}% of {len(allsig)} cells")

# The sensor-detected nose-poke (Supplementary Figure 6A legend: "duration of nose-poking, detected
# by a photoelectric sensor"). The trial table gives only the first poke's onset and the last
# poke's end, so the window above includes the pauses between pokes. The session's .eeg file
# holds the two sensor channels (97, right trials; 98, left trials; data description, page 1),
# which give every poke. Rate: spikes during the pokes of the trial's channel within the trial
# table's window, over their summed duration.
sens = raw["EE188 nose-poke sensors"]
PK = {1: pfc2.poke_intervals(sens[:, 0]), 2: pfc2.poke_intervals(sens[:, 1])}
s = S["EE188"]
RS = np.zeros((s.ntrials, s.n))
npk, nother, mis, dur = [], [], [], []
for r, (a, b, _, side) in enumerate(s.NP):
    iv = PK[int(side)][(PK[int(side)][:, 1] > a - 0.01) & (PK[int(side)][:, 0] < b + 0.01)]
    other = PK[3 - int(side)]
    nother.append(int(((other[:, 1] > a) & (other[:, 0] < b)).sum()))
    npk.append(len(iv))
    mis.append(max(abs(iv[0, 0] - a), abs(iv[-1, 1] - b)))
    dur.append((iv[:, 1] - iv[:, 0]).sum())
    for x0, x1 in iv:
        i, j = np.searchsorted(s.t, [x0, x1])
        RS[r] += np.bincount(s.u[i:j], minlength=s.n)
    RS[r] /= dur[-1]
dur = np.array(dur)
print(f"EE188 sensors: pokes detected, right channel {len(PK[1])}, left channel {len(PK[2])}; per trial {min(npk)} to "
      f"{max(npk)} pokes on the trial's channel, {sum(nother)} on the other channel within the trial windows; first onset "
      f"and last end match the trial table to within {1000 * max(mis):.1f} ms")
print(f"summed poke duration per trial: median {np.median(dur):.2f} s, range {dur.min():.2f} to {dur.max():.2f} s "
      f"(left {np.median(dur[s.left]):.2f} s, right {np.median(dur[~s.left]):.2f} s); first onset to last end median "
      f"{np.median(s.NP[:, 1] - s.NP[:, 0]):.2f} s")
for reading, eq in (("a", True), ("b", False)):
    with np.errstate(invalid="ignore", divide="ignore"):
        p = stats.ttest_ind(RS[s.left], RS[~s.left], equal_var=eq).pvalue
    sig = p < 0.05
    ODOR[("EE188", reading + "-sensor")] = (sig, RS[s.left].mean(0) - RS[~s.left].mean(0), p)
    print(f"EE188 reading {reading}, sensor-detected pokes: {100 * sig.mean():.1f}% ({sig.sum()}/{s.n}) odor selective; "
          f"within 26 +- 2.6: {abs(100 * sig.mean() - 26) <= 2.6}")

EE188 reading a: 45.3% (53/117) odor selective; within 26 +- 2.6: False
EE188 reading b: 42.7% (50/117) odor selective; within 26 +- 2.6: False
nose-poke durations EE188: median 4.6 s, range 1.9 to 29.4 s
other sessions, reading a: {'EE049': '25%', 'EE088': '21%', 'EE152': '30%', 'EE165': '31%', 'EE196': '39%', 'EE210': '40%', 'FF103': '12%', 'FF130': '36%', 'FF167': '25%', 'GG069': '29%', 'GG097': '0%', 'GG149': '11%'}
all sessions, reading a: 32.4% of 627 cells

EE188 sensors: pokes detected, right channel 78, left channel 111; per trial 1 to 13 pokes on the trial's channel, 0 on the other channel within the trial windows; first onset and last end match the trial table to within 1.6 ms
summed poke duration per trial: median 2.43 s, range 1.01 to 10.00 s (left 1.52 s, right 3.05 s); first onset to last end median 4.59 s
EE188 reading a, sensor-detected pokes: 30.8% (36/117) odor selective; within 26 +- 2.6: False
EE188 reading b, sensor-detected pokes: 29.1% (34/117) odor selective; within 26 +- 2.6: False

s = S["EE188"]
for window, key in (("sensor-detected pokes", "a-sensor"), ("first poke to end of last (variant)", "a")):
    od_sig, od_diff, _ = ODOR[("EE188", key)]
    assoc = False
    print(f"odor preference from the {window}:")
    for name, lo, hi in (("stem 0-0.3", 0.0, 0.3), ("side arms 0.5-1", 0.5, 1.0)):
        reg = np.array([selective(LR["EE188"][i], lo, hi, "a") for i in range(s.n)])
        tab = np.array([[np.sum(od_sig & reg), np.sum(od_sig & ~reg)], [np.sum(~od_sig & reg), np.sum(~od_sig & ~reg)]])
        odds, pf = stats.fisher_exact(tab)
        inr = (G >= lo) & (G <= hi)
        dreg = np.array([LR["EE188"][i][0][inr].mean() for i in range(s.n)])
        rs = stats.spearmanr(od_diff, dreg)
        rp = stats.pearsonr(od_diff, dreg)
        same_sign = (odds > 1 and pf < 0.05) or (rs.statistic > 0 and rs.pvalue < 0.05)
        assoc |= same_sign
        print(f"  {name}: table [[odor & region, odor only], [region only, neither]] = {tab.tolist()}, Fisher odds "
              f"{odds:.2f}, p = {pf:.3f}; Spearman r = {rs.statistic:.3f}, p = {rs.pvalue:.3f}; Pearson r = "
              f"{rp.statistic:.3f}, p = {rp.pvalue:.2g}")
    print("  an association with the same sign of preference at p < 0.05 in either region:", assoc)

odor preference from the sensor-detected pokes:
  stem 0-0.3: table [[odor & region, odor only], [region only, neither]] = [[18, 18], [35, 46]], Fisher odds 1.31, p = 0.549; Spearman r = -0.155, p = 0.095; Pearson r = -0.278, p = 0.0024
  side arms 0.5-1: table [[odor & region, odor only], [region only, neither]] = [[22, 14], [37, 44]], Fisher odds 1.87, p = 0.161; Spearman r = 0.237, p = 0.010; Pearson r = -0.007, p = 0.94
  an association with the same sign of preference at p < 0.05 in either region: True
odor preference from the first poke to end of last (variant):
  stem 0-0.3: table [[odor & region, odor only], [region only, neither]] = [[28, 25], [25, 39]], Fisher odds 1.75, p = 0.191; Spearman r = -0.184, p = 0.047; Pearson r = -0.411, p = 4.2e-06
  side arms 0.5-1: table [[odor & region, odor only], [region only, neither]] = [[30, 23], [29, 35]], Fisher odds 1.57, p = 0.267; Spearman r = 0.134, p = 0.150; Pearson r = -0.176, p = 0.057
  an association with the same sign of preference at p < 0.05 in either region: False

rng = np.random.default_rng(104)
RETRO = {}
for k, s in S.items():
    prev = s.left[:-1]
    cur = s.left[1:]
    RETRO[k] = [(pfc2.lr_test(H[k][i][1:], prev, rng)[3] <= 0.05, pfc2.lr_test(H[k][i][1:], cur, rng)[3] <= 0.05)
                for i in range(s.n)]
    same = np.mean(prev == cur)
    if k == "EE188":
        print(f"EE188: previous side equals current side in {100 * same:.0f}% of trials")
for label, keys in (("(a) EE188", ["EE188"]), ("(b) layer 2/3", L23), ("(b) layer 5", L5)):
    v = np.array([x for k in keys for x in RETRO[k]])
    fp, fc = 100 * v[:, 0].mean(), 100 * v[:, 1].mean()
    print(f"{label}: previous-trial labels {fp:.1f}% ({v[:, 0].sum()}/{len(v)}), current-trial labels {fc:.1f}% "
          f"({v[:, 1].sum()}); previous below half of current: {fp < fc / 2}")
for k in S:
    v = np.array(RETRO[k])
    print(f"  {k}: previous {v[:, 0].sum()}/{len(v)}, current {v[:, 1].sum()}/{len(v)}, trials {S[k].ntrials}, "
          f"previous = current in {100 * np.mean(S[k].left[:-1] == S[k].left[1:]):.0f}%")

EE188: previous side equals current side in 44% of trials

(a) EE188: previous-trial labels 4.3% (5/117), current-trial labels 77.8% (91); previous below half of current: True
(b) layer 2/3: previous-trial labels 7.3% (35/481), current-trial labels 66.9% (322); previous below half of current: True
(b) layer 5: previous-trial labels 8.9% (13/146), current-trial labels 65.8% (96); previous below half of current: True
  EE049: previous 2/44, current 19/44, trials 19, previous = current in 44%
  EE088: previous 2/28, current 15/28, trials 31, previous = current in 47%
  EE152: previous 2/61, current 35/61, trials 36, previous = current in 46%
  EE165: previous 11/70, current 38/70, trials 42, previous = current in 44%
  EE188: previous 5/117, current 91/117, trials 42, previous = current in 44%
  EE196: previous 9/66, current 57/66, trials 48, previous = current in 47%
  EE210: previous 2/60, current 47/60, trials 53, previous = current in 44%
  FF103: previous 1/42, current 12/42, trials 39, previous = current in 45%
  FF130: previous 9/64, current 51/64, trials 46, previous = current in 44%
  FF167: previous 3/40, current 33/40, trials 43, previous = current in 45%
  GG069: previous 0/14, current 10/14, trials 40, previous = current in 41%
  GG097: previous 1/12, current 5/12, trials 38, previous = current in 49%
  GG149: previous 1/9, current 5/9, trials 31, previous = current in 50%

def speed_corr(s):
    rows, cnt = [], []
    w = s.whl
    for t0, t1 in s.NP[:, 1:3]:
        for a in np.arange(t0, t1 - 0.5 + 1e-9, 0.5):
            f = np.arange(int(a * pfc2.WHL_RATE), int((a + 0.5) * pfc2.WHL_RATE) + 1)
            f = f[f < len(w)]
            ok = w[f, 0] > 0
            if ok.sum() < 3:
                continue
            two = w[f, 2] > 0
            x = np.where(two, (w[f, 0] + w[f, 2]) / 2, w[f, 0])[ok]
            y = np.where(two, (w[f, 1] + w[f, 3]) / 2, w[f, 1])[ok]
            rows.append(np.sum(np.hypot(np.diff(x), np.diff(y))) / (len(x) - 1) * pfc2.WHL_RATE)
            i, j = np.searchsorted(s.t, [a, a + 0.5])
            cnt.append(np.bincount(s.u[i:j], minlength=s.n))
    sp = np.array(rows)
    C = np.array(cnt, float)
    r = np.array([np.corrcoef(sp, C[:, i])[0, 1] if C[:, i].std() > 0 else np.nan for i in range(s.n)])
    return r, len(sp)


SPEED = {}
for k, s in S.items():
    SPEED[k], nb = speed_corr(s)
r_all = np.concatenate([SPEED[k] for k in S])
r_all = r_all[~np.isnan(r_all)]
w_all = stats.wilcoxon(r_all)
print(f"all cells: n = {len(r_all)}, median r = {np.median(r_all):.3f}, mean {r_all.mean():.3f}, "
      f"fraction positive {np.mean(r_all > 0):.2f}, Wilcoxon p = {w_all.pvalue:.2g}; "
      f"median significantly above 0: {np.median(r_all) > 0 and w_all.pvalue < 0.05}")
print(f"range of r: {r_all.min():.2f} to {r_all.max():.2f}")
for k in S:
    r = SPEED[k][~np.isnan(SPEED[k])]
    print(f"  {k}: median r = {np.median(r):.3f}, n = {len(r)}")
for r_ in ("EE", "FF", "GG"):
    r = np.concatenate([SPEED[k] for k in S if S[k].rat == r_])
    r = r[~np.isnan(r)]
    print(f"  rat {r_}: median r = {np.median(r):.3f}, n = {len(r)}, Wilcoxon p = {stats.wilcoxon(r).pvalue:.2g}")
ty = np.concatenate([S[k].celltype_file for k in S])
rr = np.concatenate([SPEED[k] for k in S])
for t, name in ((1, "pyramidal"), (2, "interneuron")):
    v = rr[(ty == t) & ~np.isnan(rr)]
    print(f"  CellsInfo {name}: median r = {np.median(v):.3f}, n = {len(v)}, Wilcoxon p = {stats.wilcoxon(v).pvalue:.2g}")

all cells: n = 627, median r = -0.018, mean -0.026, fraction positive 0.45, Wilcoxon p = 0.00036; median significantly above 0: False
range of r: -0.52 to 0.64
  EE049: median r = 0.002, n = 44
  EE088: median r = -0.066, n = 28
  EE152: median r = -0.009, n = 61
  EE165: median r = -0.024, n = 70
  EE188: median r = -0.015, n = 117
  EE196: median r = -0.044, n = 66
  EE210: median r = -0.046, n = 60
  FF103: median r = 0.016, n = 42
  FF130: median r = -0.033, n = 64
  FF167: median r = -0.013, n = 40
  GG069: median r = 0.006, n = 14
  GG097: median r = 0.029, n = 12
  GG149: median r = 0.045, n = 9
  rat EE: median r = -0.023, n = 446, Wilcoxon p = 4.3e-05
  rat FF: median r = -0.007, n = 146, Wilcoxon p = 0.15
  rat GG: median r = 0.031, n = 35, Wilcoxon p = 0.13
  CellsInfo pyramidal: median r = -0.007, n = 476, Wilcoxon p = 0.065
  CellsInfo interneuron: median r = -0.078, n = 82, Wilcoxon p = 3.2e-05

def pv_width(keys, which, norm):
    P = []
    for k in keys:
        s = S[k]
        if which == "a":
            P.append(np.concatenate([H[k][:, s.left].mean(1), H[k][:, ~s.left].mean(1)], axis=0))
        else:
            P.append(H[k].mean(1))
    P = np.concatenate(P)
    P = P[P.std(1) > 0]
    if norm == "peak":
        P = P / P.max(1, keepdims=True)
    elif norm == "z":
        P = (P - P.mean(1, keepdims=True)) / P.std(1, keepdims=True)
    C = np.corrcoef(P.T)
    wd = []
    for x in range(len(C)):
        above = C[x] > 0.5
        a = b = x
        while a > 0 and above[a - 1]:
            a -= 1
        while b < len(C) - 1 and above[b + 1]:
            b += 1
        wd.append((b - a + 1) / len(C))
    return np.array(wd), C


PV = {}
for norm in ("raw", "peak", "z"):
  for label, keys in (("EE188 (L2/3)", ["EE188"]), ("all L2/3", L23), ("all L5", L5)):
    for which in "ab":
        wd, C = pv_width(keys, which, norm)
        PV[(label, which, norm)] = (wd.mean(), wd.std(ddof=1))
        print(f"{norm:4s} {label}, reading {which}: width {wd.mean():.3f} +- {wd.std(ddof=1):.3f} (range {wd.min():.3f} to "
              f"{wd.max():.3f}); correlations {C.min():.2f} to {C.max():.2f}")
for norm in ("raw", "peak", "z"):
  for which in "ab":
    l23 = PV[("EE188 (L2/3)", which, norm)][0]
    l23p = PV[("all L2/3", which, norm)][0]
    l5 = PV[("all L5", which, norm)][0]
    print(f"{norm:4s} reading {which}: L2/3 EE188 {l23:.3f} (0.18 +- 0.087: {abs(l23 - 0.18) <= 0.087}), L2/3 pooled {l23p:.3f} "
          f"({abs(l23p - 0.18) <= 0.087}), L5 pooled {l5:.3f} (0.27 +- 0.13: {abs(l5 - 0.27) <= 0.13})")

raw  EE188 (L2/3), reading a: width 1.000 +- 0.000 (range 1.000 to 1.000); correlations 0.66 to 1.00
raw  EE188 (L2/3), reading b: width 1.000 +- 0.000 (range 1.000 to 1.000); correlations 0.73 to 1.00
raw  all L2/3, reading a: width 0.988 +- 0.048 (range 0.730 to 1.000); correlations 0.47 to 1.00
raw  all L2/3, reading b: width 0.996 +- 0.030 (range 0.745 to 1.000); correlations 0.49 to 1.00
raw  all L5, reading a: width 0.959 +- 0.126 (range 0.410 to 1.000); correlations 0.35 to 1.00
raw  all L5, reading b: width 1.000 +- 0.000 (range 1.000 to 1.000); correlations 0.59 to 1.00
peak EE188 (L2/3), reading a: width 0.239 +- 0.060 (range 0.105 to 0.335); correlations -0.38 to 1.00
peak EE188 (L2/3), reading b: width 0.265 +- 0.079 (range 0.120 to 0.385); correlations -0.31 to 1.00
peak all L2/3, reading a: width 0.243 +- 0.059 (range 0.105 to 0.315); correlations -0.36 to 1.00
peak all L2/3, reading b: width 0.268 +- 0.072 (range 0.110 to 0.355); correlations -0.39 to 1.00
peak all L5, reading a: width 0.300 +- 0.078 (range 0.100 to 0.385); correlations -0.42 to 1.00
peak all L5, reading b: width 0.338 +- 0.086 (range 0.105 to 0.455); correlations -0.44 to 1.00
z    EE188 (L2/3), reading a: width 0.182 +- 0.028 (range 0.110 to 0.230); correlations -0.45 to 1.00
z    EE188 (L2/3), reading b: width 0.188 +- 0.033 (range 0.120 to 0.245); correlations -0.52 to 1.00
z    all L2/3, reading a: width 0.179 +- 0.026 (range 0.105 to 0.210); correlations -0.41 to 1.00
z    all L2/3, reading b: width 0.184 +- 0.028 (range 0.115 to 0.220); correlations -0.46 to 1.00
z    all L5, reading a: width 0.204 +- 0.036 (range 0.095 to 0.250); correlations -0.55 to 1.00
z    all L5, reading b: width 0.216 +- 0.039 (range 0.105 to 0.275); correlations -0.61 to 1.00
raw  reading a: L2/3 EE188 1.000 (0.18 +- 0.087: False), L2/3 pooled 0.988 (False), L5 pooled 0.959 (0.27 +- 0.13: False)
raw  reading b: L2/3 EE188 1.000 (0.18 +- 0.087: False), L2/3 pooled 0.996 (False), L5 pooled 1.000 (0.27 +- 0.13: False)
peak reading a: L2/3 EE188 0.239 (0.18 +- 0.087: True), L2/3 pooled 0.243 (True), L5 pooled 0.300 (0.27 +- 0.13: True)
peak reading b: L2/3 EE188 0.265 (0.18 +- 0.087: True), L2/3 pooled 0.268 (False), L5 pooled 0.338 (0.27 +- 0.13: True)
z    reading a: L2/3 EE188 0.182 (0.18 +- 0.087: True), L2/3 pooled 0.179 (True), L5 pooled 0.204 (0.27 +- 0.13: True)
z    reading b: L2/3 EE188 0.188 (0.18 +- 0.087: True), L2/3 pooled 0.184 (True), L5 pooled 0.216 (0.27 +- 0.13: True)

CONN = {}
JOINT = {"G-b": [], "G-a": []}
GS = ("G-b", "G-a")
t0 = time.time()
for bw, tag in ((0.001, "B-a"), (0.0005, "B-b")):
    for k, s in S.items():
        rng = np.random.default_rng(200 + list(S).index(k) + (0 if tag == "B-a" else 50))
        obs, sur = pfc2.ccg_jitter(s.t, s.u, s.n, rng, bw, nsur=1000)
        for gs, gside in (("G-b", "step"), ("G-a", 0.01)):
            c = pfc2.detect_connections(obs, sur, bw, gside=gside)
            CONN[(k, tag, gs)] = {"exc": c["exc"], "inh": c["inh"], "width": c["width"]}
            JOINT[gs].append(c["joint"][~np.eye(s.n, dtype=bool)])
        del obs, sur, c
print("jitter CCGs of", sum(s.n * (s.n - 1) for s in S.values()), "ordered pairs, two bin widths, in",
      round(time.time() - t0), "s")
for gs in GS:
    v = np.concatenate(JOINT[gs])
    print(f"{gs}: fraction of surrogates lying entirely within both global bands, over pairs and bin widths: median "
          f"{np.median(v):.4f}, 5th to 95th percentile {np.percentile(v, 5):.4f} to {np.percentile(v, 95):.4f} (target 0.99)")


def conns(k, B, W, G="G-b"):
    c = CONN[(k, B, G)]
    E, I = c["exc"].copy(), c["inh"].copy()
    if W == "W-b":
        E &= c["width"] <= 2
        I &= c["width"] <= 2
    return E, I


COMBOS = [(b, w) for b in ("B-a", "B-b") for w in ("W-a", "W-b")]

# check against the authors' list for EE188
s = S["EE188"]
idx = {c: i for i, c in enumerate(s.cells)}
PA = np.zeros((s.n, s.n), bool)
for a, b in EX["pairs_all"].astype(int):
    PA[idx[a], idx[b]] = True
pre_t = s.celltype_file[PA.any(1)]
post_t = s.celltype_file[PA.any(0)]
print(f"EE188 pairs_all: {PA.sum()} pairs; presynaptic cells {PA.any(1).sum()} (types {np.bincount(pre_t, minlength=4)[1:]}), "
      f"targets {PA.any(0).sum()} (types {np.bincount(post_t, minlength=4)[1:]}) [types: pyramidal, interneuron, unknown]")
for gs in GS:
    for B, W in COMBOS:
        E, I = conns("EE188", B, W, gs)
        print(f"  {gs} {B} {W}: excitatory {E.sum()}, inhibitory {I.sum()}; of the {PA.sum()} pairs_all found as "
              f"excitatory {(PA & E).sum()}, as inhibitory {(PA & I).sum()}")

jitter CCGs of 40240 ordered pairs, two bin widths, in 268 s
G-b: fraction of surrogates lying entirely within both global bands, over pairs and bin widths: median 0.9940, 5th to 95th percentile 0.9900 to 0.9990 (target 0.99)
G-a: fraction of surrogates lying entirely within both global bands, over pairs and bin widths: median 0.9940, 5th to 95th percentile 0.9830 to 0.9990 (target 0.99)
EE188 pairs_all: 43 pairs; presynaptic cells 34 (types [33  1  0]), targets 13 (types [4 7 2]) [types: pyramidal, interneuron, unknown]
  G-b B-a W-a: excitatory 95, inhibitory 35; of the 43 pairs_all found as excitatory 41, as inhibitory 0
  G-b B-a W-b: excitatory 91, inhibitory 30; of the 43 pairs_all found as excitatory 40, as inhibitory 0
  G-b B-b W-a: excitatory 81, inhibitory 9; of the 43 pairs_all found as excitatory 34, as inhibitory 0
  G-b B-b W-b: excitatory 80, inhibitory 9; of the 43 pairs_all found as excitatory 34, as inhibitory 0
  G-a B-a W-a: excitatory 105, inhibitory 43; of the 43 pairs_all found as excitatory 43, as inhibitory 0
  G-a B-a W-b: excitatory 101, inhibitory 38; of the 43 pairs_all found as excitatory 42, as inhibitory 0
  G-a B-b W-a: excitatory 88, inhibitory 9; of the 43 pairs_all found as excitatory 34, as inhibitory 0
  G-a B-b W-b: excitatory 87, inhibitory 9; of the 43 pairs_all found as excitatory 34, as inhibitory 0

s = S["EE188"]
rng = np.random.default_rng(250)
T = max(s.t.max(), s.duration) + 1
NULLSHIFT = {}
npairs = s.n * (s.n - 1)
for shift in ("circular", "rank x 25 ms"):
    off = rng.uniform(10, T - 10, s.n) if shift == "circular" else 0.025 * rng.permutation(s.n)
    tsh = (s.t + off[s.u]) % T
    o = np.argsort(tsh, kind="stable")
    obs, sur = pfc2.ccg_jitter(tsh[o], s.u[o], s.n, rng, 0.001, nsur=1000)
    for gs, gside in (("G-b", "step"), ("G-a", 0.01)):
        c = pfc2.detect_connections(obs, sur, 0.001, gside=gside)
        NULLSHIFT[(shift, gs)] = (c["exc"].sum() / npairs, c["inh"].sum() / npairs)
        E0, I0 = conns("EE188", "B-a", "W-a", gs)
        print(f"EE188 shifted ({shift}), {gs} B-a W-a: excitatory {c['exc'].sum()} ({100 * NULLSHIFT[(shift, gs)][0]:.2f}% "
              f"of {npairs} directed pairs), inhibitory {c['inh'].sum()} ({100 * NULLSHIFT[(shift, gs)][1]:.2f}%); "
              f"unshifted: excitatory {E0.sum()} ({100 * E0.sum() / npairs:.2f}%), inhibitory {I0.sum()} ({100 * I0.sum() / npairs:.2f}%)")
    del obs, sur, c

EE188 shifted (circular), G-b B-a W-a: excitatory 20 (0.15% of 13572 directed pairs), inhibitory 9 (0.07%); unshifted: excitatory 95 (0.70%), inhibitory 35 (0.26%)

EE188 shifted (circular), G-a B-a W-a: excitatory 25 (0.18% of 13572 directed pairs), inhibitory 12 (0.09%); unshifted: excitatory 105 (0.77%), inhibitory 43 (0.32%)

EE188 shifted (rank x 25 ms), G-b B-a W-a: excitatory 16 (0.12% of 13572 directed pairs), inhibitory 6 (0.04%); unshifted: excitatory 95 (0.70%), inhibitory 35 (0.26%)

EE188 shifted (rank x 25 ms), G-a B-a W-a: excitatory 18 (0.13% of 13572 directed pairs), inhibitory 10 (0.07%); unshifted: excitatory 105 (0.77%), inhibitory 43 (0.32%)

NPAIRS = sum(s.n * (s.n - 1) for s in S.values())
NCELLS = sum(s.n for s in S.values())
print(f"directed pairs of simultaneously recorded cells: {NPAIRS} (paper 62,408; within 10%: {abs(NPAIRS - 62408) <= 6240.8})")
for gs, B, W in [(g, b, w) for g in GS for b, w in COMBOS]:
    ne = ni = nany = nrec = 0
    for k in S:
        E, I = conns(k, B, W, gs)
        A = E | I
        ne += E.sum()
        ni += I.sum()
        nany += A.sum()
        nrec += np.triu(A & A.T, 1).sum()
    pa, pe, pi = 100 * nany / NPAIRS, 100 * ne / NPAIRS, 100 * ni / NPAIRS
    ra, rb = 100 * nrec / NPAIRS, 100 * nrec / (NPAIRS / 2)
    ok3 = abs(pa - 0.79) <= 0.079 and abs(pe - 0.55) <= 0.055 and abs(pi - 0.24) <= 0.024
    print(f"{gs} {B} {W}: any {nany} ({pa:.2f}%), excitatory {ne} ({pe:.2f}%), inhibitory {ni} ({pi:.2f}%), "
          f"reciprocal pairs {nrec} = {ra:.3f}% of directed (a), {rb:.3f}% of unordered (b); "
          f"three directed values within tolerance: {ok3}; reciprocal within 0.17 +- 0.017: (a) {abs(ra - 0.17) <= 0.017}, (b) {abs(rb - 0.17) <= 0.017}")
print("per rat (G-b B-a W-a):")
for r in ("EE", "FF", "GG"):
    ks = [k for k in S if S[k].rat == r]
    npr = sum(S[k].n * (S[k].n - 1) for k in ks)
    ne = sum(conns(k, "B-a", "W-a")[0].sum() for k in ks)
    ni = sum(conns(k, "B-a", "W-a")[1].sum() for k in ks)
    print(f"  {r}: {npr} pairs, excitatory {100 * ne / npr:.2f}%, inhibitory {100 * ni / npr:.2f}%")
for (shift, gs), (fe, fi) in NULLSHIFT.items():
    print(f"false positives on EE188 with shifted trains ({shift}, {gs} B-a W-a): excitatory {100 * fe:.2f}%, inhibitory "
          f"{100 * fi:.2f}% of directed pairs")

directed pairs of simultaneously recorded cells: 40240 (paper 62,408; within 10%: False)
G-b B-a W-a: any 530 (1.32%), excitatory 389 (0.97%), inhibitory 158 (0.39%), reciprocal pairs 80 = 0.199% of directed (a), 0.398% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
G-b B-a W-b: any 471 (1.17%), excitatory 364 (0.90%), inhibitory 118 (0.29%), reciprocal pairs 40 = 0.099% of directed (a), 0.199% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
G-b B-b W-a: any 398 (0.99%), excitatory 349 (0.87%), inhibitory 53 (0.13%), reciprocal pairs 27 = 0.067% of directed (a), 0.134% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
G-b B-b W-b: any 380 (0.94%), excitatory 339 (0.84%), inhibitory 44 (0.11%), reciprocal pairs 18 = 0.045% of directed (a), 0.089% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
G-a B-a W-a: any 605 (1.50%), excitatory 439 (1.09%), inhibitory 191 (0.47%), reciprocal pairs 91 = 0.226% of directed (a), 0.452% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
G-a B-a W-b: any 545 (1.35%), excitatory 413 (1.03%), inhibitory 151 (0.38%), reciprocal pairs 50 = 0.124% of directed (a), 0.249% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
G-a B-b W-a: any 447 (1.11%), excitatory 387 (0.96%), inhibitory 67 (0.17%), reciprocal pairs 36 = 0.089% of directed (a), 0.179% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) True
G-a B-b W-b: any 429 (1.07%), excitatory 377 (0.94%), inhibitory 57 (0.14%), reciprocal pairs 25 = 0.062% of directed (a), 0.124% of unordered (b); three directed values within tolerance: False; reciprocal within 0.17 +- 0.017: (a) False, (b) False
per rat (G-b B-a W-a):
  EE: 32540 pairs, excitatory 0.94%, inhibitory 0.38%
  FF: 7314 pairs, excitatory 1.07%, inhibitory 0.46%
  GG: 386 pairs, excitatory 1.30%, inhibitory 0.00%
false positives on EE188 with shifted trains (circular, G-b B-a W-a): excitatory 0.15%, inhibitory 0.07% of directed pairs
false positives on EE188 with shifted trains (circular, G-a B-a W-a): excitatory 0.18%, inhibitory 0.09% of directed pairs
false positives on EE188 with shifted trains (rank x 25 ms, G-b B-a W-a): excitatory 0.12%, inhibitory 0.04% of directed pairs
false positives on EE188 with shifted trains (rank x 25 ms, G-a B-a W-a): excitatory 0.13%, inhibitory 0.07% of directed pairs

def celltypes(k, B, W, C, G="G-b"):
    E, I = conns(k, B, W, G)
    pyr, inn = E.any(1), I.any(1)
    both = pyr & inn
    return pyr & ~both, inn & ~both, both


TYPES = {}
for gs in GS:
    for B, W in COMBOS:
        for C in ("C-a", "C-b"):
            npyr = nint = nboth = ntot = 0
            tgt_int = tgt = 0
            per_s, per_r = [], {}
            for k in S:
                pyr, inn, both = celltypes(k, B, W, C, gs)
                TYPES[(k, B, W, C, gs)] = (pyr, inn)
                E, I = conns(k, B, W, gs)
                npyr += pyr.sum()
                nint += inn.sum()
                nboth += both.sum()
                ntot += S[k].n - (both.sum() if C == "C-b" else 0)
                targets = E.any(0)
                tgt += targets.sum()
                tgt_int += (targets & I.any(1)).sum()
                if inn.sum() > 0:
                    per_s.append(pyr.sum() / inn.sum())
                r = S[k].rat
                per_r.setdefault(r, [0, 0])
                per_r[r][0] += pyr.sum()
                per_r[r][1] += inn.sum()
            rs = np.array(per_s)
            rr = np.array([a / b for a, b in per_r.values() if b > 0])
            fp, fi, ft = 100 * npyr / ntot, 100 * nint / ntot, 100 * tgt_int / tgt
            ok = (abs(fp - 32.5) <= 3.25, abs(fi - 12.5) <= 1.25, abs(ft - 50.7) <= 5.07,
                  abs(rs.mean() - 2.82) <= 0.51 if len(rs) else False, abs(rr.mean() - 2.82) <= 0.51 if len(rr) else False)
            print(f"{gs} {B} {W} {C}: pyramidal {fp:.1f}% ({npyr}/{ntot}), interneurons {fi:.1f}% ({nint}), both {nboth}; "
                  f"targets of excitation that are interneurons {ft:.1f}% ({tgt_int}/{tgt}); ratio per session "
                  f"{rs.mean() if len(rs) else np.nan:.2f} +- {rs.std(ddof=1) if len(rs) > 1 else np.nan:.2f} (n = {len(rs)} of {len(S)}), per rat "
                  f"{rr.mean() if len(rr) else np.nan:.2f} +- {rr.std(ddof=1) if len(rr) > 1 else np.nan:.2f} (n = {len(rr)}); within tolerance "
                  f"[pyr, int, targets, ratio a, ratio b]: {ok}")
ci = S["EE188"].celltype_file
print(f"EE188 CellsInfo of the example file: pyramidal {np.mean(ci == 1) * 100:.1f}%, interneuron "
      f"{np.mean(ci == 2) * 100:.1f}%, unknown {np.mean(ci == 3) * 100:.1f}% of {len(ci)}")
for gs in GS:
    for B, W in COMBOS:
        pyr, inn = TYPES[("EE188", B, W, "C-a", gs)]
        print(f"  {gs} {B} {W}: agreement with CellsInfo; our pyramidal cells typed 1/2/3 {np.bincount(ci[pyr], minlength=4)[1:]}, "
              f"our interneurons {np.bincount(ci[inn], minlength=4)[1:]}")
# speed-not-correlated, variant: the rate-speed correlations by the connection-based types (G-b B-a W-a C-a)
for j, name in ((0, "pyramidal"), (1, "interneuron")):
    v = np.concatenate([SPEED[k][TYPES[(k, "B-a", "W-a", "C-a", "G-b")][j]] for k in S])
    v = v[~np.isnan(v)]
    print(f"speed-not-correlated, connection-based {name} cells: median r = {np.median(v):.3f}, n = {len(v)}, "
          f"Wilcoxon p = {stats.wilcoxon(v).pvalue:.2g}")

G-b B-a W-a C-a: pyramidal 36.5% (229/627), interneurons 5.3% (33), both 37; targets of excitation that are interneurons 29.2% (52/178); ratio per session 6.90 +- 1.80 (n = 9 of 13), per rat 7.06 +- 0.62 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-b B-a W-a C-b: pyramidal 38.8% (229/590), interneurons 5.6% (33), both 37; targets of excitation that are interneurons 29.2% (52/178); ratio per session 6.90 +- 1.80 (n = 9 of 13), per rat 7.06 +- 0.62 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-b B-a W-b C-a: pyramidal 36.2% (227/627), interneurons 6.1% (38), both 28; targets of excitation that are interneurons 27.8% (49/176); ratio per session 6.28 +- 2.23 (n = 9 of 13), per rat 5.41 +- 1.29 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-b B-a W-b C-b: pyramidal 37.9% (227/599), interneurons 6.3% (38), both 28; targets of excitation that are interneurons 27.8% (49/176); ratio per session 6.28 +- 2.23 (n = 9 of 13), per rat 5.41 +- 1.29 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-b B-b W-a C-a: pyramidal 35.6% (223/627), interneurons 2.2% (14), both 21; targets of excitation that are interneurons 18.2% (28/154); ratio per session 16.48 +- 11.15 (n = 8 of 13), per rat 14.24 +- 7.98 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.True_, np.False_, np.False_, np.False_, np.False_)
G-b B-b W-a C-b: pyramidal 36.8% (223/606), interneurons 2.3% (14), both 21; targets of excitation that are interneurons 18.2% (28/154); ratio per session 16.48 +- 11.15 (n = 8 of 13), per rat 14.24 +- 7.98 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-b B-b W-b C-a: pyramidal 35.2% (221/627), interneurons 2.2% (14), both 19; targets of excitation that are interneurons 16.9% (26/154); ratio per session 16.57 +- 11.60 (n = 7 of 13), per rat 14.15 +- 5.16 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.True_, np.False_, np.False_, np.False_, np.False_)
G-b B-b W-b C-b: pyramidal 36.3% (221/608), interneurons 2.3% (14), both 19; targets of excitation that are interneurons 16.9% (26/154); ratio per session 16.57 +- 11.60 (n = 7 of 13), per rat 14.15 +- 5.16 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-a W-a C-a: pyramidal 38.4% (241/627), interneurons 5.4% (34), both 51; targets of excitation that are interneurons 31.8% (61/192); ratio per session 7.73 +- 3.53 (n = 9 of 13), per rat 7.00 +- 0.21 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-a W-a C-b: pyramidal 41.8% (241/576), interneurons 5.9% (34), both 51; targets of excitation that are interneurons 31.8% (61/192); ratio per session 7.73 +- 3.53 (n = 9 of 13), per rat 7.00 +- 0.21 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-a W-b C-a: pyramidal 38.0% (238/627), interneurons 5.9% (37), both 45; targets of excitation that are interneurons 30.4% (58/191); ratio per session 6.53 +- 2.27 (n = 9 of 13), per rat 6.22 +- 0.13 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-a W-b C-b: pyramidal 40.9% (238/582), interneurons 6.4% (37), both 45; targets of excitation that are interneurons 30.4% (58/191); ratio per session 6.53 +- 2.27 (n = 9 of 13), per rat 6.22 +- 0.13 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-b W-a C-a: pyramidal 38.1% (239/627), interneurons 2.2% (14), both 26; targets of excitation that are interneurons 19.5% (31/159); ratio per session 18.94 +- 10.32 (n = 9 of 13), per rat 15.65 +- 4.45 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-b W-a C-b: pyramidal 39.8% (239/601), interneurons 2.3% (14), both 26; targets of excitation that are interneurons 19.5% (31/159); ratio per session 18.94 +- 10.32 (n = 9 of 13), per rat 15.65 +- 4.45 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-b W-b C-a: pyramidal 37.8% (237/627), interneurons 2.2% (14), both 24; targets of excitation that are interneurons 18.2% (29/159); ratio per session 18.54 +- 9.74 (n = 8 of 13), per rat 16.67 +- 0.47 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
G-a B-b W-b C-b: pyramidal 39.3% (237/603), interneurons 2.3% (14), both 24; targets of excitation that are interneurons 18.2% (29/159); ratio per session 18.54 +- 9.74 (n = 8 of 13), per rat 16.67 +- 0.47 (n = 2); within tolerance [pyr, int, targets, ratio a, ratio b]: (np.False_, np.False_, np.False_, np.False_, np.False_)
EE188 CellsInfo of the example file: pyramidal 28.2%, interneuron 12.8%, unknown 59.0% of 117
  G-b B-a W-a: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [33  3 20], our interneurons [0 5 3]
  G-b B-a W-b: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [32  3 20], our interneurons [0 5 4]
  G-b B-b W-a: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [28  5 19], our interneurons [0 3 0]
  G-b B-b W-b: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [28  5 19], our interneurons [0 4 0]
  G-a B-a W-a: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [33  2 21], our interneurons [0 6 4]
  G-a B-a W-b: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [32  2 21], our interneurons [0 6 5]
  G-a B-b W-a: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [28  5 22], our interneurons [0 2 0]
  G-a B-b W-b: agreement with CellsInfo; our pyramidal cells typed 1/2/3 [28  5 22], our interneurons [0 3 0]
speed-not-correlated, connection-based pyramidal cells: median r = -0.015, n = 229, Wilcoxon p = 0.041
speed-not-correlated, connection-based interneuron cells: median r = -0.088, n = 33, Wilcoxon p = 0.065

for gs in GS:
    for B, W in COMBOS:
        for C in ("C-a", "C-b"):
            out = []
            for lay, keys in (("L2/3", L23), ("L5", L5)):
                npp = nall = npyrpairs = 0
                for k in keys:
                    E, I = conns(k, B, W, gs)
                    pyr, inn = TYPES[(k, B, W, C, gs)]
                    npp += (E & pyr[:, None] & pyr[None, :]).sum()
                    nall += S[k].n * (S[k].n - 1)
                    npyrpairs += pyr.sum() * (pyr.sum() - 1)
                out.append((lay, npp, 100 * npp / nall, 100 * npp / max(npyrpairs, 1)))
            ok_a = abs(out[0][2] - 0.12) <= 0.012 and abs(out[1][2] - 0.27) <= 0.027
            ok_b = abs(out[0][3] - 0.12) <= 0.012 and abs(out[1][3] - 0.27) <= 0.027
            print(f"{gs} {B} {W} {C}: " + "; ".join(f"{l} {n} connections, {a:.3f}% of all pairs, {b:.2f}% of pyramidal pairs" for l, n, a, b in out)
                  + f"; within tolerance (a) {ok_a}, (b) {ok_b}")

G-b B-a W-a C-a: L2/3 58 connections, 0.176% of all pairs, 0.98% of pyramidal pairs; L5 11 connections, 0.150% of all pairs, 1.19% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-a W-a C-b: L2/3 58 connections, 0.176% of all pairs, 0.98% of pyramidal pairs; L5 11 connections, 0.150% of all pairs, 1.19% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-a W-b C-a: L2/3 57 connections, 0.173% of all pairs, 0.98% of pyramidal pairs; L5 16 connections, 0.219% of all pairs, 1.89% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-a W-b C-b: L2/3 57 connections, 0.173% of all pairs, 0.98% of pyramidal pairs; L5 16 connections, 0.219% of all pairs, 1.89% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-b W-a C-a: L2/3 48 connections, 0.146% of all pairs, 0.81% of pyramidal pairs; L5 9 connections, 0.123% of all pairs, 1.21% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-b W-a C-b: L2/3 48 connections, 0.146% of all pairs, 0.81% of pyramidal pairs; L5 9 connections, 0.123% of all pairs, 1.21% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-b W-b C-a: L2/3 50 connections, 0.152% of all pairs, 0.86% of pyramidal pairs; L5 8 connections, 0.109% of all pairs, 1.15% of pyramidal pairs; within tolerance (a) False, (b) False
G-b B-b W-b C-b: L2/3 50 connections, 0.152% of all pairs, 0.86% of pyramidal pairs; L5 8 connections, 0.109% of all pairs, 1.15% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-a W-a C-a: L2/3 55 connections, 0.167% of all pairs, 0.88% of pyramidal pairs; L5 14 connections, 0.191% of all pairs, 1.33% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-a W-a C-b: L2/3 55 connections, 0.167% of all pairs, 0.88% of pyramidal pairs; L5 14 connections, 0.191% of all pairs, 1.33% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-a W-b C-a: L2/3 54 connections, 0.164% of all pairs, 0.88% of pyramidal pairs; L5 16 connections, 0.219% of all pairs, 1.67% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-a W-b C-b: L2/3 54 connections, 0.164% of all pairs, 0.88% of pyramidal pairs; L5 16 connections, 0.219% of all pairs, 1.67% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-b W-a C-a: L2/3 54 connections, 0.164% of all pairs, 0.82% of pyramidal pairs; L5 18 connections, 0.246% of all pairs, 1.78% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-b W-a C-b: L2/3 54 connections, 0.164% of all pairs, 0.82% of pyramidal pairs; L5 18 connections, 0.246% of all pairs, 1.78% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-b W-b C-a: L2/3 58 connections, 0.176% of all pairs, 0.90% of pyramidal pairs; L5 17 connections, 0.232% of all pairs, 1.78% of pyramidal pairs; within tolerance (a) False, (b) False
G-a B-b W-b C-b: L2/3 58 connections, 0.176% of all pairs, 0.90% of pyramidal pairs; L5 17 connections, 0.232% of all pairs, 1.78% of pyramidal pairs; within tolerance (a) False, (b) False

s = S["EE188"]
print(f"EE188: {s.n} cells, {s.n * (s.n - 1)} ordered pairs (paper: 1172 pairs, 13,572 matrix points)")
for gs in GS:
    for B, W in COMBOS:
        E, I = conns("EE188", B, W, gs)
        pyr, inn = TYPES[("EE188", B, W, "C-a", gs)]
        A = E | I
        a_p, a_i = A[pyr].sum(), A[inn].sum()
        b_p, b_i = pyr.sum(), inn.sum()
        print(f"{gs} {B} {W}: (a) connections from pyramidal references {a_p}, from interneurons {a_i}; (b) reference cells "
              f"{b_p} and {b_i}; within 48 +- 4.8 and 30 +- 3: (a) {abs(a_p - 48) <= 4.8 and abs(a_i - 30) <= 3}, "
              f"(b) {abs(b_p - 48) <= 4.8 and abs(b_i - 30) <= 3}")
ci = s.celltype_file
print(f"authors' pairs_all with CellsInfo types: {PA.sum()} connections, from pyramidal references {PA[ci == 1].sum()}, "
      f"from interneurons {PA[ci == 2].sum()}, from unknown {PA[ci == 3].sum()}; reference cells "
      f"{PA.any(1)[ci == 1].sum()}, {PA.any(1)[ci == 2].sum()} and {PA.any(1)[ci == 3].sum()}")

EE188: 117 cells, 13572 ordered pairs (paper: 1172 pairs, 13,572 matrix points)
G-b B-a W-a: (a) connections from pyramidal references 81, from interneurons 18; (b) reference cells 56 and 8; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-b B-a W-b: (a) connections from pyramidal references 79, from interneurons 17; (b) reference cells 55 and 9; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-b B-b W-a: (a) connections from pyramidal references 78, from interneurons 6; (b) reference cells 52 and 3; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-b B-b W-b: (a) connections from pyramidal references 78, from interneurons 7; (b) reference cells 52 and 4; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-a B-a W-a: (a) connections from pyramidal references 86, from interneurons 21; (b) reference cells 56 and 10; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-a B-a W-b: (a) connections from pyramidal references 84, from interneurons 20; (b) reference cells 55 and 11; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-a B-b W-a: (a) connections from pyramidal references 84, from interneurons 5; (b) reference cells 55 and 2; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
G-a B-b W-b: (a) connections from pyramidal references 84, from interneurons 6; (b) reference cells 55 and 3; within 48 +- 4.8 and 30 +- 3: (a) False, (b) False
authors' pairs_all with CellsInfo types: 43 connections, from pyramidal references 40, from interneurons 3, from unknown 0; reference cells 33, 1 and 0

def largest_component(A):
    g = nx.Graph()
    g.add_edges_from(zip(*np.nonzero(A)))
    return max((len(c) for c in nx.connected_components(g)), default=0)


for gs in GS:
    for B, W in COMBOS:
        E, I = conns("EE188", B, W, gs)
        n = largest_component(E | I)
        print(f"{gs} {B} {W}: largest connected circuit {n} of {S['EE188'].n} cells; within 39 +- 3.9: {abs(n - 39) <= 3.9}")
print(f"authors' pairs_all graph: largest component {largest_component(PA)} cells")

G-b B-a W-a: largest connected circuit 83 of 117 cells; within 39 +- 3.9: False
G-b B-a W-b: largest connected circuit 82 of 117 cells; within 39 +- 3.9: False
G-b B-b W-a: largest connected circuit 59 of 117 cells; within 39 +- 3.9: False
G-b B-b W-b: largest connected circuit 59 of 117 cells; within 39 +- 3.9: False
G-a B-a W-a: largest connected circuit 86 of 117 cells; within 39 +- 3.9: False
G-a B-a W-b: largest connected circuit 85 of 117 cells; within 39 +- 3.9: False
G-a B-b W-a: largest connected circuit 64 of 117 cells; within 39 +- 3.9: False
G-a B-b W-b: largest connected circuit 64 of 117 cells; within 39 +- 3.9: False
authors' pairs_all graph: largest component 39 cells

def cp_interval(k, n, level=0.99):
    a = 1 - level
    lo = stats.beta.ppf(a / 2, k, n - k + 1) if k > 0 else 0.0
    hi = stats.beta.ppf(1 - a / 2, k + 1, n - k) if k < n else 1.0
    return lo, hi


for gs, (B, W) in [(g, c) for g in GS for c in COMBOS]:
    cnt = {}
    for k, s in S.items():
        E, I = conns(k, B, W, gs)
        sep = np.abs(s.shank[:, None] - s.shank[None, :]) * 200
        off = ~np.eye(s.n, dtype=bool)
        for d in np.unique(sep[off]):
            m = (sep == d) & off
            c = cnt.setdefault(d, [0, 0, 0])
            c[0] += m.sum()
            c[1] += (E & m).sum()
            c[2] += (I & m).sum()
    ds = sorted(cnt)
    same = cnt[0][1] + cnt[0][2]
    other = sum(cnt[d][1] + cnt[d][2] for d in ds if d > 0)
    p0 = (cnt[0][1] + cnt[0][2]) / cnt[0][0]
    lo0, hi0 = cp_interval(cnt[0][1] + cnt[0][2], cnt[0][0])
    lo2, hi2 = cp_interval(cnt[200][1] + cnt[200][2], cnt[200][0])
    p2 = (cnt[200][1] + cnt[200][2]) / cnt[200][0]
    far = max(d for d in ds if cnt[d][1] + cnt[d][2] > 0)
    print(f"{gs} {B} {W}: connected pairs on the same shank {same}, on other shanks {other}; same-shank probability "
          f"{100 * p0:.2f}% [{100 * lo0:.2f}, {100 * hi0:.2f}] against 200 um {100 * p2:.2f}% [{100 * lo2:.2f}, {100 * hi2:.2f}]; "
          f"largest separation with a connection {far} um; criterion met: {same > other and lo0 > hi2 and far == 1200}")
    if (B, W) == ("B-a", "W-a"):
        print(f"    {gs} B-a W-a by separation:")
        for d in ds:
            e_lo, e_hi = cp_interval(cnt[d][1], cnt[d][0])
            i_lo, i_hi = cp_interval(cnt[d][2], cnt[d][0])
            print(f"    {d:5d} um: {cnt[d][0]:6d} pairs, excitatory {100 * cnt[d][1] / cnt[d][0]:.2f}% [{100 * e_lo:.2f}, "
                  f"{100 * e_hi:.2f}], inhibitory {100 * cnt[d][2] / cnt[d][0]:.2f}% [{100 * i_lo:.2f}, {100 * i_hi:.2f}]")

G-b B-a W-a: connected pairs on the same shank 352, on other shanks 195; same-shank probability 6.20% [5.41, 7.07] against 200 um 1.24% [0.97, 1.56]; largest separation with a connection 1400 um; criterion met: False
    G-b B-a W-a by separation:
        0 um:   5676 pairs, excitatory 4.40% [3.73, 5.15], inhibitory 1.80% [1.37, 2.30]
      200 um:   9678 pairs, excitatory 0.88% [0.65, 1.15], inhibitory 0.36% [0.22, 0.55]
      400 um:   7656 pairs, excitatory 0.31% [0.17, 0.52], inhibitory 0.14% [0.06, 0.30]
      600 um:   6682 pairs, excitatory 0.21% [0.09, 0.40], inhibitory 0.01% [0.00, 0.11]
      800 um:   4960 pairs, excitatory 0.20% [0.07, 0.43], inhibitory 0.08% [0.01, 0.25]
     1000 um:   3318 pairs, excitatory 0.06% [0.00, 0.28], inhibitory 0.09% [0.01, 0.33]
     1200 um:   1766 pairs, excitatory 0.11% [0.01, 0.52], inhibitory 0.00% [0.00, 0.30]
     1400 um:    504 pairs, excitatory 0.40% [0.02, 1.83], inhibitory 0.40% [0.02, 1.83]
G-b B-a W-b: connected pairs on the same shank 295, on other shanks 187; same-shank probability 5.20% [4.47, 6.00] against 200 um 1.16% [0.90, 1.47]; largest separation with a connection 1400 um; criterion met: False
G-b B-b W-a: connected pairs on the same shank 264, on other shanks 138; same-shank probability 4.65% [3.96, 5.42] against 200 um 0.83% [0.61, 1.09]; largest separation with a connection 1400 um; criterion met: False
G-b B-b W-b: connected pairs on the same shank 248, on other shanks 135; same-shank probability 4.37% [3.70, 5.12] against 200 um 0.80% [0.58, 1.06]; largest separation with a connection 1400 um; criterion met: False
G-a B-a W-a: connected pairs on the same shank 396, on other shanks 234; same-shank probability 6.98% [6.13, 7.89] against 200 um 1.44% [1.14, 1.78]; largest separation with a connection 1400 um; criterion met: False
    G-a B-a W-a by separation:
        0 um:   5676 pairs, excitatory 4.90% [4.19, 5.68], inhibitory 2.08% [1.62, 2.62]
      200 um:   9678 pairs, excitatory 1.02% [0.78, 1.32], inhibitory 0.41% [0.26, 0.61]
      400 um:   7656 pairs, excitatory 0.37% [0.21, 0.58], inhibitory 0.22% [0.11, 0.40]
      600 um:   6682 pairs, excitatory 0.24% [0.11, 0.44], inhibitory 0.03% [0.00, 0.14]
      800 um:   4960 pairs, excitatory 0.24% [0.10, 0.49], inhibitory 0.14% [0.04, 0.35]
     1000 um:   3318 pairs, excitatory 0.06% [0.00, 0.28], inhibitory 0.12% [0.02, 0.38]
     1200 um:   1766 pairs, excitatory 0.11% [0.01, 0.52], inhibitory 0.06% [0.00, 0.42]
     1400 um:    504 pairs, excitatory 0.40% [0.02, 1.83], inhibitory 0.40% [0.02, 1.83]
G-a B-a W-b: connected pairs on the same shank 339, on other shanks 225; same-shank probability 5.97% [5.19, 6.83] against 200 um 1.34% [1.06, 1.67]; largest separation with a connection 1400 um; criterion met: False
G-a B-b W-a: connected pairs on the same shank 288, on other shanks 166; same-shank probability 5.07% [4.35, 5.87] against 200 um 0.94% [0.71, 1.22]; largest separation with a connection 1400 um; criterion met: False
G-a B-b W-b: connected pairs on the same shank 272, on other shanks 162; same-shank probability 4.79% [4.09, 5.57] against 200 um 0.90% [0.67, 1.18]; largest separation with a connection 1400 um; criterion met: False

for gs in GS:
    for B, W in COMBOS:
        for C in ("C-a", "C-b"):
            na = ka = kb = 0
            for k, s in S.items():
                E, I = conns(k, B, W, gs)
                pyr, inn = TYPES[(k, B, W, C, gs)]
                m = (s.shank[:, None] == s.shank[None, :]) & pyr[:, None] & inn[None, :]
                na += m.sum()
                ka += (E & m).sum()
                kb += ((E | I | E.T | I.T) & m).sum()
            pa_, pb_ = 100 * ka / max(na, 1), 100 * kb / max(na, 1)
            print(f"{gs} {B} {W} {C}: {na} same-shank pyramidal-interneuron pairs; (a) {ka} = {pa_:.1f}%, (b) {kb} = {pb_:.1f}%; "
                  f"in 0.9 to 2.2%: (a) {0.9 <= pa_ <= 2.2}, (b) {0.9 <= pb_ <= 2.2}")
# variant: CellsInfo types and the authors' pairs_all, EE188 only
s = S["EE188"]
ci = s.celltype_file
m = (s.shank[:, None] == s.shank[None, :]) & (ci == 1)[:, None] & (ci == 2)[None, :]
print(f"EE188, CellsInfo types and pairs_all: {m.sum()} same-shank pyramidal-interneuron pairs; (a) {(PA & m).sum()} = "
      f"{100 * (PA & m).sum() / m.sum():.1f}%, (b) {((PA | PA.T) & m).sum()} = {100 * ((PA | PA.T) & m).sum() / m.sum():.1f}%")
for B, W in COMBOS:
    E, I = conns("EE188", B, W)
    print(f"EE188, CellsInfo types and the connections of G-b {B} {W}: (a) {(E & m).sum()} = {100 * (E & m).sum() / m.sum():.1f}%, "
          f"(b) {((E | I | E.T | I.T) & m).sum()} = {100 * ((E | I | E.T | I.T) & m).sum() / m.sum():.1f}%")

G-b B-a W-a C-a: 134 same-shank pyramidal-interneuron pairs; (a) 52 = 38.8%, (b) 59 = 44.0%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-a W-a C-b: 134 same-shank pyramidal-interneuron pairs; (a) 52 = 38.8%, (b) 59 = 44.0%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-a W-b C-a: 163 same-shank pyramidal-interneuron pairs; (a) 62 = 38.0%, (b) 70 = 42.9%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-a W-b C-b: 163 same-shank pyramidal-interneuron pairs; (a) 62 = 38.0%, (b) 70 = 42.9%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-b W-a C-a: 76 same-shank pyramidal-interneuron pairs; (a) 47 = 61.8%, (b) 48 = 63.2%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-b W-a C-b: 76 same-shank pyramidal-interneuron pairs; (a) 47 = 61.8%, (b) 48 = 63.2%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-b W-b C-a: 79 same-shank pyramidal-interneuron pairs; (a) 50 = 63.3%, (b) 51 = 64.6%; in 0.9 to 2.2%: (a) False, (b) False
G-b B-b W-b C-b: 79 same-shank pyramidal-interneuron pairs; (a) 50 = 63.3%, (b) 51 = 64.6%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-a W-a C-a: 139 same-shank pyramidal-interneuron pairs; (a) 47 = 33.8%, (b) 54 = 38.8%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-a W-a C-b: 139 same-shank pyramidal-interneuron pairs; (a) 47 = 33.8%, (b) 54 = 38.8%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-a W-b C-a: 157 same-shank pyramidal-interneuron pairs; (a) 53 = 33.8%, (b) 61 = 38.9%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-a W-b C-b: 157 same-shank pyramidal-interneuron pairs; (a) 53 = 33.8%, (b) 61 = 38.9%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-b W-a C-a: 64 same-shank pyramidal-interneuron pairs; (a) 35 = 54.7%, (b) 36 = 56.2%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-b W-a C-b: 64 same-shank pyramidal-interneuron pairs; (a) 35 = 54.7%, (b) 36 = 56.2%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-b W-b C-a: 69 same-shank pyramidal-interneuron pairs; (a) 41 = 59.4%, (b) 42 = 60.9%; in 0.9 to 2.2%: (a) False, (b) False
G-a B-b W-b C-b: 69 same-shank pyramidal-interneuron pairs; (a) 41 = 59.4%, (b) 42 = 60.9%; in 0.9 to 2.2%: (a) False, (b) False
EE188, CellsInfo types and pairs_all: 54 same-shank pyramidal-interneuron pairs; (a) 23 = 42.6%, (b) 23 = 42.6%
EE188, CellsInfo types and the connections of G-b B-a W-a: (a) 25 = 46.3%, (b) 28 = 51.9%
EE188, CellsInfo types and the connections of G-b B-a W-b: (a) 24 = 44.4%, (b) 27 = 50.0%
EE188, CellsInfo types and the connections of G-b B-b W-a: (a) 25 = 46.3%, (b) 25 = 46.3%
EE188, CellsInfo types and the connections of G-b B-b W-b: (a) 25 = 46.3%, (b) 25 = 46.3%

def traj_frac(keys, B, W, C, gs):
    out = {}
    for name, j in (("pyr", 0), ("int", 1)):
        v = np.concatenate([[LR[k][i][3] <= 0.05 for i in np.flatnonzero(TYPES[(k, B, W, C, gs)][j])] for k in keys])
        out[name] = (100 * v.mean() if len(v) else np.nan, int(v.sum()), len(v))
    return out


for gs, (B, W) in [(g, c) for g in GS for c in COMBOS]:
    for C in ("C-a",):
        a23 = traj_frac(["EE188"], B, W, C, gs)
        b23 = traj_frac(L23, B, W, C, gs)
        b5 = traj_frac(L5, B, W, C, gs)
        ok_b = (abs(b23["pyr"][0] - 72.2) <= 7.22 and abs(b5["pyr"][0] - 84.7) <= 8.47 and
                abs(b23["int"][0] - 58.3) <= 5.83 and abs(b5["int"][0] - 83) <= 8.3)
        print(f"{gs} {B} {W} {C}: (a) EE188 pyramidal {a23['pyr'][0]:.1f}% ({a23['pyr'][1]}/{a23['pyr'][2]}), interneurons "
              f"{a23['int'][0]:.1f}% ({a23['int'][1]}/{a23['int'][2]}); (b) L2/3 pyramidal {b23['pyr'][0]:.1f}% "
              f"({b23['pyr'][1]}/{b23['pyr'][2]}), interneurons {b23['int'][0]:.1f}% ({b23['int'][1]}/{b23['int'][2]}); "
              f"L5 pyramidal {b5['pyr'][0]:.1f}% ({b5['pyr'][1]}/{b5['pyr'][2]}), interneurons {b5['int'][0]:.1f}% "
              f"({b5['int'][1]}/{b5['int'][2]}); all four within tolerance under (b): {ok_b}")
allv = np.concatenate([[LR[k][i][3] <= 0.05 for i in range(S[k].n)] for k in S])
print(f"all cells, any type: {100 * allv.mean():.1f}% trajectory specific")

G-b B-a W-a C-a: (a) EE188 pyramidal 76.8% (43/56), interneurons 50.0% (4/8); (b) L2/3 pyramidal 69.0% (127/184), interneurons 59.3% (16/27); L5 pyramidal 73.3% (33/45), interneurons 83.3% (5/6); all four within tolerance under (b): False
G-b B-a W-b C-a: (a) EE188 pyramidal 76.4% (42/55), interneurons 55.6% (5/9); (b) L2/3 pyramidal 69.8% (127/182), interneurons 60.7% (17/28); L5 pyramidal 75.6% (34/45), interneurons 90.0% (9/10); all four within tolerance under (b): False
G-b B-b W-a C-a: (a) EE188 pyramidal 69.2% (36/52), interneurons 100.0% (3/3); (b) L2/3 pyramidal 64.4% (116/180), interneurons 88.9% (8/9); L5 pyramidal 67.4% (29/43), interneurons 100.0% (5/5); all four within tolerance under (b): False
G-b B-b W-b C-a: (a) EE188 pyramidal 69.2% (36/52), interneurons 100.0% (4/4); (b) L2/3 pyramidal 65.4% (117/179), interneurons 90.0% (9/10); L5 pyramidal 66.7% (28/42), interneurons 100.0% (4/4); all four within tolerance under (b): False
G-a B-a W-a C-a: (a) EE188 pyramidal 75.0% (42/56), interneurons 50.0% (5/10); (b) L2/3 pyramidal 70.2% (134/191), interneurons 55.6% (15/27); L5 pyramidal 72.0% (36/50), interneurons 85.7% (6/7); all four within tolerance under (b): False
G-a B-a W-b C-a: (a) EE188 pyramidal 74.5% (41/55), interneurons 54.5% (6/11); (b) L2/3 pyramidal 70.9% (134/189), interneurons 55.2% (16/29); L5 pyramidal 73.5% (36/49), interneurons 87.5% (7/8); all four within tolerance under (b): False
G-a B-b W-a C-a: (a) EE188 pyramidal 70.9% (39/55), interneurons 100.0% (2/2); (b) L2/3 pyramidal 65.6% (124/189), interneurons 70.0% (7/10); L5 pyramidal 70.0% (35/50), interneurons 100.0% (4/4); all four within tolerance under (b): False
G-a B-b W-b C-a: (a) EE188 pyramidal 70.9% (39/55), interneurons 100.0% (3/3); (b) L2/3 pyramidal 66.5% (125/188), interneurons 72.7% (8/11); L5 pyramidal 69.4% (34/49), interneurons 100.0% (3/3); all four within tolerance under (b): False
all cells, any type: 66.5% trajectory specific

def lr_class(res):
    D0, point, glob, gp = res
    segs = pfc2.segments(point, glob)
    left = any(D0[a:b + 1].mean() > 0 for a, b in segs)
    right = any(D0[a:b + 1].mean() < 0 for a, b in segs)
    return "both" if left and right else "left" if left else "right" if right else "none"


def cluster_test(k):
    s = S[k]
    cls = np.array([lr_class(r) for r in LR[k]])
    labs = ["left", "right", "both", "none"]
    sh = np.unique(s.shank)
    tab = np.array([[np.sum((s.shank == h) & (cls == c)) for c in labs] for h in sh])
    tab = tab[:, tab.sum(0) > 0]
    tab = tab[tab.sum(1) > 0]
    chi2, p, dof, _ = stats.chi2_contingency(tab)
    v = np.sqrt(chi2 / (tab.sum() * (min(tab.shape) - 1)))
    return chi2, p, v, tab, cls


chi2, p, v, tab, cls = cluster_test("EE188")
print(f"EE188: classes {dict(zip(*np.unique(cls, return_counts=True)))}; shanks x classes table {tab.tolist()}; "
      f"chi-square {chi2:.1f}, p = {p:.3f}, Cramer's V = {v:.2f}; association at p < 0.05: {p < 0.05}")
for k in S:
    if k != "EE188":
        chi2, p, v, tab, cls = cluster_test(k)
        print(f"  {k}: chi-square {chi2:.1f}, p = {p:.3f}, V = {v:.2f}, {S[k].n} cells")

EE188: classes {np.str_('both'): np.int64(15), np.str_('left'): np.int64(37), np.str_('none'): np.int64(36), np.str_('right'): np.int64(29)}; shanks x classes table [[5, 5, 0, 4], [4, 5, 5, 5], [4, 5, 1, 6], [8, 1, 2, 4], [4, 5, 2, 7], [3, 4, 2, 1], [7, 2, 2, 4], [2, 2, 1, 5]]; chi-square 19.1, p = 0.577, Cramer's V = 0.23; association at p < 0.05: False
  EE049: chi-square 16.6, p = 0.279, V = 0.43, 44 cells
  EE088: chi-square 15.2, p = 0.652, V = 0.42, 28 cells
  EE152: chi-square 11.9, p = 0.942, V = 0.26, 61 cells
  EE165: chi-square 21.0, p = 0.460, V = 0.32, 70 cells
  EE196: chi-square 14.6, p = 0.843, V = 0.27, 66 cells
  EE210: chi-square 25.4, p = 0.230, V = 0.38, 60 cells
  FF103: chi-square 16.5, p = 0.351, V = 0.36, 42 cells
  FF130: chi-square 29.1, p = 0.047, V = 0.39, 64 cells
  FF167: chi-square 25.8, p = 0.105, V = 0.46, 40 cells
  GG069: chi-square 7.0, p = 0.858, V = 0.41, 14 cells
  GG097: chi-square 8.9, p = 0.449, V = 0.50, 12 cells
  GG149: chi-square 7.8, p = 0.653, V = 0.66, 9 cells

POSDEP = {}
t0 = time.time()
rng = np.random.default_rng(300)
pairs_needed = set()
# the G-b connections are a subset of the G-a ones, so the tests are run on G-a and selected below
for B in ("B-a", "B-b"):
    for k in S:
        assert not (conns(k, B, "W-a")[0] & ~conns(k, B, "W-a", "G-a")[0]).any()
        E, _ = conns(k, B, "W-a", "G-a")
        pairs_needed |= {(k, i, j) for i, j in zip(*np.nonzero(E))}
for k, i, j in sorted(pairs_needed):
    s = S[k]
    res = {}
    for name, m in (("L", s.left), ("R", ~s.left), ("all", np.ones(s.ntrials, bool))):
        S0, mu, point, glob, gp = pfc2.coincidence_test(s, i, j, rng, trials_mask=m)
        res[name] = (gp <= 0.01, bool(glob.any() and (~point).any()), point, glob, S0, mu)
    POSDEP[(k, i, j)] = res
print(f"coincidence tests on {len(POSDEP)} excitatory connections x 3 trial sets in {round(time.time() - t0)} s")


def posdep(key, trials, dep):
    r = POSDEP[key]
    sets = ("L", "R") if trials == "a" else ("all",)
    return any(r[t][0] if dep == "i" else r[t][1] for t in sets)

coincidence tests on 573 excitatory connections x 3 trial sets in 157 s

for gs in GS:
    for B in ("B-a", "B-b"):
        keys = [(k, i, j) for k in S for i, j in zip(*np.nonzero(conns(k, B, "W-a", gs)[0]))]
        for trials in "ab":
            for dep in ("i", "ii"):
                n = sum(posdep(x, trials, dep) for x in keys)
                ok = (abs(len(keys) - 343) <= 34.3 and abs(n - 67) <= 6.7) or abs(100 * n / len(keys) - 19.5) <= 2.0
                print(f"{gs} {B} W-a, trials {trials}, dependence {dep:2s}: {n} of {len(keys)} excitatory connections "
                      f"({100 * n / len(keys):.1f}%); counts or fraction within tolerance: {ok}")
keys = [(k, i, j) for k in S for i, j in zip(*np.nonzero(conns(k, "B-a", "W-a")[0]))]
for r in ("EE", "FF", "GG"):
    kk = [x for x in keys if S[x[0]].rat == r]
    print(f"  rat {r}: {sum(posdep(x, 'a', 'i') for x in kk)} of {len(kk)} (G-b B-a, a, i)")
# variant: jittering the presynaptic cell only (Supplementary Figure 12B), G-b B-a W-a
rng1 = np.random.default_rng(302)
PRE = {}
for k, i, j in keys:
    s = S[k]
    PRE[(k, i, j)] = {name: pfc2.coincidence_test(s, i, j, rng1, trials_mask=m, jitter_post=False)[4] <= 0.01
                      for name, m in (("L", s.left), ("R", ~s.left), ("all", np.ones(s.ntrials, bool)))}
for trials in "ab":
    n = sum(any(PRE[x][t] for t in (("L", "R") if trials == "a" else ("all",))) for x in keys)
    print(f"presynaptic jitter only, G-b B-a W-a, trials {trials}, dependence i: {n} of {len(keys)} ({100 * n / len(keys):.1f}%)")
# null check: the same test with pre spikes shifted by 20 ms, which removes the monosynaptic lag
rng0 = np.random.default_rng(301)
sub = keys[::10]
fp = 0
for k, i, j in sub:
    s = S[k]
    ps = pfc2.pair_spikes(s, i, j)
    ps["tpre"] = ps["tpre"] + 0.02
    fp += pfc2.coincidence_test(s, i, j, rng0, ps=ps)[4] <= 0.01
print(f"null check, presynaptic spikes shifted by 20 ms: {fp} of {len(sub)} pairs significant (level 0.01)")

G-b B-a W-a, trials a, dependence i : 31 of 389 excitatory connections (8.0%); counts or fraction within tolerance: False
G-b B-a W-a, trials a, dependence ii: 31 of 389 excitatory connections (8.0%); counts or fraction within tolerance: False
G-b B-a W-a, trials b, dependence i : 34 of 389 excitatory connections (8.7%); counts or fraction within tolerance: False
G-b B-a W-a, trials b, dependence ii: 33 of 389 excitatory connections (8.5%); counts or fraction within tolerance: False
G-b B-b W-a, trials a, dependence i : 30 of 349 excitatory connections (8.6%); counts or fraction within tolerance: False
G-b B-b W-a, trials a, dependence ii: 30 of 349 excitatory connections (8.6%); counts or fraction within tolerance: False
G-b B-b W-a, trials b, dependence i : 36 of 349 excitatory connections (10.3%); counts or fraction within tolerance: False
G-b B-b W-a, trials b, dependence ii: 35 of 349 excitatory connections (10.0%); counts or fraction within tolerance: False
G-a B-a W-a, trials a, dependence i : 32 of 439 excitatory connections (7.3%); counts or fraction within tolerance: False
G-a B-a W-a, trials a, dependence ii: 32 of 439 excitatory connections (7.3%); counts or fraction within tolerance: False
G-a B-a W-a, trials b, dependence i : 35 of 439 excitatory connections (8.0%); counts or fraction within tolerance: False
G-a B-a W-a, trials b, dependence ii: 34 of 439 excitatory connections (7.7%); counts or fraction within tolerance: False
G-a B-b W-a, trials a, dependence i : 31 of 387 excitatory connections (8.0%); counts or fraction within tolerance: False
G-a B-b W-a, trials a, dependence ii: 31 of 387 excitatory connections (8.0%); counts or fraction within tolerance: False
G-a B-b W-a, trials b, dependence i : 36 of 387 excitatory connections (9.3%); counts or fraction within tolerance: False
G-a B-b W-a, trials b, dependence ii: 35 of 387 excitatory connections (9.0%); counts or fraction within tolerance: False
  rat EE: 27 of 306 (G-b B-a, a, i)
  rat FF: 4 of 78 (G-b B-a, a, i)
  rat GG: 0 of 5 (G-b B-a, a, i)

presynaptic jitter only, G-b B-a W-a, trials a, dependence i: 46 of 389 (11.8%)
presynaptic jitter only, G-b B-a W-a, trials b, dependence i: 53 of 389 (13.6%)

null check, presynaptic spikes shifted by 20 ms: 0 of 39 pairs significant (level 0.01)

for gs in GS:
    for B in ("B-a", "B-b"):
        kk = [("EE188", i, j) for i, j in zip(*np.nonzero(conns("EE188", B, "W-a", gs)[0]))]
        for trials in "ab":
            for dep in ("i", "ii"):
                n = sum(posdep(x, trials, dep) for x in kk)
                print(f"{gs} {B} W-a, trials {trials}, dependence {dep:2s}: {n} of {len(kk)} connections; within 15 +- 1.5 and "
                      f"36 +- 3.6: {abs(n - 15) <= 1.5 and abs(len(kk) - 36) <= 3.6}")
kk = [("EE188", i, j) for i, j in zip(*np.nonzero(PA))]
for trials in "ab":
    for dep in ("i", "ii"):
        n = sum(posdep(x, trials, dep) for x in kk if x in POSDEP)
        print(f"pairs_all ({len(kk)} pairs, {sum(x in POSDEP for x in kk)} tested), trials {trials}, dependence {dep}: {n}")

G-b B-a W-a, trials a, dependence i : 7 of 95 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-a W-a, trials a, dependence ii: 7 of 95 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-a W-a, trials b, dependence i : 6 of 95 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-a W-a, trials b, dependence ii: 6 of 95 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-b W-a, trials a, dependence i : 9 of 81 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-b W-a, trials a, dependence ii: 9 of 81 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-b W-a, trials b, dependence i : 7 of 81 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-b B-b W-a, trials b, dependence ii: 7 of 81 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-a W-a, trials a, dependence i : 7 of 105 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-a W-a, trials a, dependence ii: 7 of 105 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-a W-a, trials b, dependence i : 6 of 105 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-a W-a, trials b, dependence ii: 6 of 105 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-b W-a, trials a, dependence i : 10 of 88 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-b W-a, trials a, dependence ii: 10 of 88 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-b W-a, trials b, dependence i : 7 of 88 connections; within 15 +- 1.5 and 36 +- 3.6: False
G-a B-b W-a, trials b, dependence ii: 7 of 88 connections; within 15 +- 1.5 and 36 +- 3.6: False
pairs_all (43 pairs, 43 tested), trials a, dependence i: 6
pairs_all (43 pairs, 43 tested), trials a, dependence ii: 6
pairs_all (43 pairs, 43 tested), trials b, dependence i: 5
pairs_all (43 pairs, 43 tested), trials b, dependence ii: 5

def example_pair(a, b, trials, seed):
    s = S["EE188"]
    i, j = idx[a], idx[b]
    m = s.left if trials == "L" else ~s.left
    S0, mu, point, glob, gp = pfc2.coincidence_test(s, i, j, np.random.default_rng(seed), trials_mask=m)
    segs = [(round(G[x] - 0.5 / pfc2.NGRID, 3), round(G[y] + 0.5 / pfc2.NGRID, 3)) for x, y in pfc2.segments(point, glob)]
    gl = np.flatnonzero(glob)
    gext = (round(G[gl[0]], 3), round(G[gl[-1]], 3)) if len(gl) else None
    return segs, gext, gp, S0, mu


s = S["EE188"]
print("shanks: 156 on", s.shank[idx[156]], ", 136 on", s.shank[idx[136]], "; types (CellsInfo)",
      s.celltype_file[idx[156]], s.celltype_file[idx[136]])
E, I = conns("EE188", "B-a", "W-a")
print("connection 156 -> 136 (G-b B-a W-a): excitatory", E[idx[156], idx[136]], "; in pairs_all:", PA[idx[156], idx[136]])
segs, gext, gp, S0, mu = example_pair(156, 136, "L", 401)
print(f"left trials: global p = {gp:.4f}; significant regions (pointwise extent) {segs}; global band broken "
      f"from {gext[0] if gext else None} to {gext[1] if gext else None}")
PLp = H["EE188"][idx[156]][s.left].mean(0)
PLq = H["EE188"][idx[136]][s.left].mean(0)
r25 = (G >= 0.2) & (G <= 0.5)
for c, P in ((156, PLp), (136, PLq)):
    v = P[r25]
    print(f"cell {c} left PSTH 0.2-0.5: {v.min():.1f} to {v.max():.1f} (spikes per unit length per trial); "
          f"variation {(v.max() - v.min()) / v.max() * 100:.0f}% of the maximum (steady if < 50%)")
if segs:
    st, en = segs[0]
    print(f"only one region: {len(segs) == 1}; start {st} within 0.4 +- 0.04: {abs(st - 0.4) <= 0.04}; end {en} within "
          f"0.6 +- 0.06: {abs(en - 0.6) <= 0.06}")
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(G, S0, label="coincidences")
ax.plot(G, mu, label="jitter mean")
ax.set_xlabel("normalised position")
ax.set_ylabel("events per unit length per trial")
ax.legend()
ax.set_title("156-136, left trials")
plt.show()
print(f"plotted range: coincidences {S0.min():.2f} to {S0.max():.2f}, jitter mean {mu.min():.2f} to {mu.max():.2f}")

shanks: 156 on 5 , 136 on 4 ; types (CellsInfo) 1 3
connection 156 -> 136 (G-b B-a W-a): excitatory True ; in pairs_all: True
left trials: global p = 0.0010; significant regions (pointwise extent) [(np.float64(0.325), np.float64(0.645))]; global band broken from 0.372 to 0.622
cell 156 left PSTH 0.2-0.5: 66.6 to 129.4 (spikes per unit length per trial); variation 49% of the maximum (steady if < 50%)
cell 136 left PSTH 0.2-0.5: 24.7 to 60.8 (spikes per unit length per trial); variation 59% of the maximum (steady if < 50%)
only one region: True; start 0.325 within 0.4 +- 0.04: False; end 0.645 within 0.6 +- 0.06: True

plotted range: coincidences 0.00 to 13.56, jitter mean 0.00 to 5.42

res5b = {}
for (a, b), seed in (((197, 201), 402), ((49, 52), 403), ((186, 201), 404)):
    for tr in ("L", "R"):
        segs, gext, gp, S0, mu = example_pair(a, b, tr, seed)
        res5b[(a, b, tr)] = segs
        e0 = mu[(G >= 0) & (G <= 0.5)]
        print(f"{a}-{b} {tr}: global p = {gp:.4f}, regions {segs}, global band broken {gext}; jitter mean (expected "
              f"coincidences) 0-0.5 from {e0.min():.2f} to {e0.max():.2f}; excitatory connection (G-b B-a W-a) "
              f"{E[idx[a], idx[b]]}")
s1 = res5b[(197, 201, "L")]
c1 = len(s1) == 1 and abs(s1[0][0] - 0.4) <= 0.04 and abs(s1[0][1] - 0.5) <= 0.05
c2 = len(res5b[(49, 52, "L")]) > 0 and len(res5b[(49, 52, "R")]) == 0
c3 = len(res5b[(186, 201, "R")]) == 1 and res5b[(186, 201, "R")][0][0] > 0.75 and len(res5b[(186, 201, "L")]) == 0
c3b = len(res5b[(186, 201, "R")]) >= 1 and all(x[0] > 0.75 for x in res5b[(186, 201, "R")])
print(f"197-201 meets its expected values: {c1}; 49-52 left only: {c2}; 186-201 only region on right trials beyond 0.75: "
      f"{c3} (right-trial regions all beyond 0.75: {c3b}); pairs holding: {int(c1) + int(c2) + int(c3)}")

197-201 L: global p = 0.0070, regions [(np.float64(0.355), np.float64(0.525))], global band broken (np.float64(0.398), np.float64(0.432)); jitter mean (expected coincidences) 0-0.5 from 1.07 to 2.38; excitatory connection (G-b B-a W-a) True

197-201 R: global p = 0.1299, regions [], global band broken None; jitter mean (expected coincidences) 0-0.5 from 0.12 to 4.06; excitatory connection (G-b B-a W-a) True

49-52 L: global p = 0.5465, regions [], global band broken None; jitter mean (expected coincidences) 0-0.5 from 0.00 to 1.05; excitatory connection (G-b B-a W-a) True
49-52 R: global p = 0.7982, regions [], global band broken None; jitter mean (expected coincidences) 0-0.5 from 0.00 to 0.60; excitatory connection (G-b B-a W-a) True

186-201 L: global p = 0.0649, regions [], global band broken None; jitter mean (expected coincidences) 0-0.5 from 0.26 to 3.82; excitatory connection (G-b B-a W-a) True

186-201 R: global p = 0.0010, regions [(np.float64(0.715), np.float64(0.975))], global band broken (np.float64(0.752), np.float64(0.922)); jitter mean (expected coincidences) 0-0.5 from 0.31 to 5.15; excitatory connection (G-b B-a W-a) True
197-201 meets its expected values: False; 49-52 left only: False; 186-201 only region on right trials beyond 0.75: False (right-trial regions all beyond 0.75: False); pairs holding: 0

def thin_pair(key, tr, lo, hi, rng):
    k, i, j = key
    s = S[k]
    m = s.left if tr == "L" else ~s.left
    ps = pfc2.pair_spikes(s, i, j, m)
    kp, fp_ = pfc2.thin(ps["tpre"], ps["xpre"], ps["inpre"], lo, hi, rng)
    kq, fq_ = pfc2.thin(ps["tpost"], ps["xpost"], ps["inpost"], lo, hi, rng)
    a = kp | ~ps["inpre"]
    b = kq | ~ps["inpost"]
    p2 = {"tpre": ps["tpre"][a], "xpre": ps["xpre"][a], "inpre": kp[a], "tpost": ps["tpost"][b], "N": ps["N"]}
    S0, mu, point, glob, gp = pfc2.coincidence_test(s, i, j, rng, ps=p2)
    reg = (G >= lo) & (G <= hi)
    return bool((glob & reg).any() and (~point & reg).any()), min(fp_, fq_)


rng = np.random.default_rng(500)
cand = []
for key in keys:
    for tr in ("L", "R"):
        if POSDEP[key][tr][0]:
            cand.append((key, tr))
print(f"{len(cand)} position-dependent connection and trial-set combinations (G-b B-a W-a)")
THIN = {}
for key, tr in cand:
    point, glob = POSDEP[key][tr][2], POSDEP[key][tr][3]
    gl = np.flatnonzero(glob)
    a0, a1 = G[gl[0]], G[gl[-1]]
    span = a1 - a0 + 0.05
    lo_b, hi_b = (a0 - span, a1 + 0.025) if a0 - span >= 0 else (a0 - 0.025, min(a1 + span, 1.0))
    for reading, (lo, hi) in (("a", (0.2, 0.5)), ("b", (max(lo_b, 0.0), hi_b))):
        if reading == "a" and not (glob & (G >= 0.2) & (G <= 0.5)).any():
            continue
        THIN[(key, tr, reading)] = [thin_pair(key, tr, lo, hi, rng) for rep in range(10)]
for reading in "ab":
    v = {x: int(sum(ok and fl >= 0.9 for ok, fl in r)) for x, r in THIN.items() if x[2] == reading}
    v0 = {x: int(sum(ok for ok, fl in r)) for x, r in THIN.items() if x[2] == reading}
    held = sum(n >= 6 for n in v.values())
    print(f"reading {reading}: {len(v)} pairs tested, {held} keep a significant and a non-significant part, with both cells "
          f"flattened, in the majority of 10 thinnings; counts per pair {sorted(v.values(), reverse=True)}; at least two: "
          f"{held >= 2}. Without the flatness condition {sum(n >= 6 for n in v0.values())} pairs, counts "
          f"{sorted(v0.values(), reverse=True)}")
    fl = np.array([f for x, r in THIN.items() if x[2] == reading for ok, f in r])
    print(f"  min/max ratio of the less flat of the two thinned PSTHs: median {np.nanmedian(fl):.2f}, range {np.nanmin(fl):.2f} "
          f"to {np.nanmax(fl):.2f}; target reached in {np.mean(fl >= 0.9) * 100:.0f}% of {len(fl)} thinnings")
    for x, r in THIN.items():
        if x[2] == reading and v0[x] >= 6:
            k, i, j = x[0]
            print(f"  {k} {S[k].cells[i]}-{S[k].cells[j]} {x[1]}: held in {v0[x]} of 10, ratios "
                  f"{[round(float(f), 2) for ok, f in r]}")

36 position-dependent connection and trial-set combinations (G-b B-a W-a)

reading a: 16 pairs tested, 0 keep a significant and a non-significant part, with both cells flattened, in the majority of 10 thinnings; counts per pair [3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; at least two: False. Without the flatness condition 2 pairs, counts [10, 6, 5, 5, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  min/max ratio of the less flat of the two thinned PSTHs: median 0.64, range 0.00 to 0.91; target reached in 26% of 160 thinnings
  FF130 57-51 L: held in 6 of 10, ratios [0.04, 0.08, 0.04, 0.04, 0.08, 0.08, 0.04, 0.04, 0.04, 0.08]
  FF130 57-60 L: held in 10 of 10, ratios [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
reading b: 36 pairs tested, 1 keep a significant and a non-significant part, with both cells flattened, in the majority of 10 thinnings; counts per pair [10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; at least two: False. Without the flatness condition 1 pairs, counts [10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  min/max ratio of the less flat of the two thinned PSTHs: median 0.57, range 0.00 to 0.91; target reached in 23% of 360 thinnings
  EE196 8-7 R: held in 10 of 10, ratios [0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9]

# Synthetic check of the test: a presynaptic cell firing trains of 1 to 3 spikes 10 ms apart after
# silences of 0.3 to 1.5 s, a target firing at 5 Hz plus a spike 2 ms after a presynaptic spike with
# probability q_first (first spikes) or q_later (later spikes), 1,800 s.
def synth(q_first, q_later, rng):
    starts = np.cumsum(rng.uniform(0.3, 1.5, 1500))
    pre, lab = [], []
    for t in starts:
        n = rng.integers(1, 4)
        pre += [t + 0.01 * m for m in range(n)]
        lab += [m == 0 for m in range(n)]
    pre, lab = np.array(pre), np.array(lab)
    q = np.where(lab, q_first, q_later)
    post = np.sort(np.concatenate([rng.uniform(0, starts[-1] + 1, int(5 * (starts[-1] + 1))),
                                   pre[rng.random(len(pre)) < q] + 0.002]))
    return pre, post


rng = np.random.default_rng(602)
for qf, ql, name in ((0.1, 0.1, "equal efficacy"), (0.05, 0.2, "later spikes 4 times as effective"),
                     (0.2, 0.05, "first spikes 4 times as effective")):
    out = []
    for rep_ in range(20):
        pre, post = synth(qf, ql, rng)
        f, l = pfc2.train_labels(pre)
        r = pfc2.order_effect(pre, post, f, l, rng, "i", nsur=100, nperm=500)
        out.append((r[2], r[7], r[3] < 0.10, r[6]))
    d, m, sg, sd = (np.array(x) for x in zip(*out))
    print(f"synthetic, {name}, 20 runs: difference first - later mean {d.mean():.1f}, mean permuted difference {m.mean():.1f}; "
          f"significant at P < 0.10 as depression {np.sum(sg & (sd > 0))}, as facilitation {np.sum(sg & (sd < 0))}")

t0 = time.time()
rng = np.random.default_rng(600)
F6 = {}
for k, s in S.items():
    E, I = conns(k, "B-a", "W-a", "G-a")  # the G-b connections are a subset
    for i, j in zip(*np.nonzero(E | I)):
        tpre, tpost = s.spikes_of(i), s.spikes_of(j)
        for ch in ("x", "y"):
            f, l = pfc2.train_labels(tpre, chained=(ch == "x"))
            if f.sum() < 5 or l.sum() < 5:
                continue
            for how in ("i", "ii"):
                F6[(k, i, j, ch, how)] = pfc2.order_effect(tpre, tpost, f, l, rng, how, nsur=100, nperm=500) + (bool(E[i, j]),)
print(f"order effects for {len(F6)} pair-reading combinations in {round(time.time() - t0)} s")
for gs in GS:
    for pairs in "ab":
        for ch in ("x", "y"):
            for how in ("i", "ii"):
                CE = {k: conns(k, "B-a", "W-a", gs) for k in S}
                inc = lambda key: CE[key[0]][0][key[1], key[2]] or (pairs == "b" and CE[key[0]][1][key[1], key[2]])
                rows = [(key, v) for key, v in F6.items() if key[3] == ch and key[4] == how and inc(key)]
                ntot = sum(int(CE[k][0].sum() if pairs == "a" else (CE[k][0] | CE[k][1]).sum()) for k in S)
                dep = [key for key, v in rows if v[3] < 0.10 and v[6] > 0]
                fac = [key for key, v in rows if v[3] < 0.10 and v[6] < 0]
                sig = dep + fac
                diff = [key for key in sig if S[key[0]].shank[key[1]] != S[key[0]].shank[key[2]]]
                pd, pf_ = 100 * len(dep) / ntot, 100 * len(fac) / ntot
                psh = 100 * len(diff) / max(len(sig), 1)
                ok = abs(pd - 12.7) <= 1.27 and abs(pf_ - 10.7) <= 1.07 and abs(psh - 32.2) <= 3.22
                off = np.median([v[7] for key, v in rows if np.isfinite(v[7])])
                print(f"{gs} pairs {pairs}, later {ch}, height {how:2s}: {len(rows)} of {ntot} pairs testable; depression {len(dep)} "
                      f"({pd:.1f}%), facilitation {len(fac)} ({pf_:.1f}%); significant pairs on different shanks {len(diff)} "
                      f"({psh:.1f}%); significant either way {len(sig)} ({100 * len(sig) / ntot:.1f}%); all three within tolerance: {ok}; median of the mean permuted difference {off:.2f}")
# null rate: labels assigned at random to the same numbers of spikes
rng0 = np.random.default_rng(601)
nul = []
for key in [x for x in F6 if x[3] == "x" and x[4] == "i"][::5]:
    k, i, j = key[:3]
    tpre, tpost = S[k].spikes_of(i), S[k].spikes_of(j)
    f, l = pfc2.train_labels(tpre)
    sel = np.flatnonzero(f | l)
    perm = rng0.permutation(sel)
    f2 = np.zeros(len(tpre), bool)
    l2 = np.zeros(len(tpre), bool)
    f2[perm[:f.sum()]] = True
    l2[perm[f.sum():]] = True
    nul.append(pfc2.order_effect(tpre, tpost, f2, l2, rng0, "i", nsur=100, nperm=500)[3] < 0.10)
bt = stats.binomtest(int(np.sum(nul)), len(nul), 0.10)
ci = bt.proportion_ci(0.95)
print(f"null check, labels assigned at random: {np.sum(nul)} of {len(nul)} pairs significant at P < 0.10 "
      f"({100 * np.mean(nul):.1f}%, 95% interval {100 * ci.low:.1f} to {100 * ci.high:.1f}%; expected {0.1 * len(nul):.1f} "
      f"+- {np.sqrt(0.09 * len(nul)):.1f}; binomial test against 10%: p = {bt.pvalue:.2f})")

synthetic, equal efficacy, 20 runs: difference first - later mean 1.0, mean permuted difference -2.6; significant at P < 0.10 as depression 3, as facilitation 0

synthetic, later spikes 4 times as effective, 20 runs: difference first - later mean -45.5, mean permuted difference -9.9; significant at P < 0.10 as depression 0, as facilitation 20

synthetic, first spikes 4 times as effective, 20 runs: difference first - later mean 42.4, mean permuted difference 5.9; significant at P < 0.10 as depression 20, as facilitation 0

order effects for 2416 pair-reading combinations in 393 s
G-b pairs a, later x, height i : 388 of 389 pairs testable; depression 23 (5.9%), facilitation 20 (5.1%); significant pairs on different shanks 12 (27.9%); significant either way 43 (11.1%); all three within tolerance: False; median of the mean permuted difference 1.38
G-b pairs a, later x, height ii: 388 of 389 pairs testable; depression 20 (5.1%), facilitation 26 (6.7%); significant pairs on different shanks 15 (32.6%); significant either way 46 (11.8%); all three within tolerance: False; median of the mean permuted difference 0.83
G-b pairs a, later y, height i : 388 of 389 pairs testable; depression 21 (5.4%), facilitation 36 (9.3%); significant pairs on different shanks 19 (33.3%); significant either way 57 (14.7%); all three within tolerance: False; median of the mean permuted difference 0.38
G-b pairs a, later y, height ii: 388 of 389 pairs testable; depression 17 (4.4%), facilitation 41 (10.5%); significant pairs on different shanks 16 (27.6%); significant either way 58 (14.9%); all three within tolerance: False; median of the mean permuted difference 0.24
G-b pairs b, later x, height i : 529 of 530 pairs testable; depression 29 (5.5%), facilitation 29 (5.5%); significant pairs on different shanks 16 (27.6%); significant either way 58 (10.9%); all three within tolerance: False; median of the mean permuted difference 0.78
G-b pairs b, later x, height ii: 529 of 530 pairs testable; depression 24 (4.5%), facilitation 38 (7.2%); significant pairs on different shanks 19 (30.6%); significant either way 62 (11.7%); all three within tolerance: False; median of the mean permuted difference 0.54
G-b pairs b, later y, height i : 529 of 530 pairs testable; depression 26 (4.9%), facilitation 46 (8.7%); significant pairs on different shanks 23 (31.9%); significant either way 72 (13.6%); all three within tolerance: False; median of the mean permuted difference 0.19
G-b pairs b, later y, height ii: 529 of 530 pairs testable; depression 21 (4.0%), facilitation 56 (10.6%); significant pairs on different shanks 21 (27.3%); significant either way 77 (14.5%); all three within tolerance: False; median of the mean permuted difference 0.31
G-a pairs a, later x, height i : 438 of 439 pairs testable; depression 28 (6.4%), facilitation 24 (5.5%); significant pairs on different shanks 17 (32.7%); significant either way 52 (11.8%); all three within tolerance: False; median of the mean permuted difference 1.27
G-a pairs a, later x, height ii: 438 of 439 pairs testable; depression 26 (5.9%), facilitation 27 (6.2%); significant pairs on different shanks 19 (35.8%); significant either way 53 (12.1%); all three within tolerance: False; median of the mean permuted difference 0.77
G-a pairs a, later y, height i : 438 of 439 pairs testable; depression 26 (5.9%), facilitation 38 (8.7%); significant pairs on different shanks 23 (35.9%); significant either way 64 (14.6%); all three within tolerance: False; median of the mean permuted difference 0.34
G-a pairs a, later y, height ii: 438 of 439 pairs testable; depression 23 (5.2%), facilitation 43 (9.8%); significant pairs on different shanks 19 (28.8%); significant either way 66 (15.0%); all three within tolerance: False; median of the mean permuted difference 0.18
G-a pairs b, later x, height i : 604 of 605 pairs testable; depression 36 (6.0%), facilitation 34 (5.6%); significant pairs on different shanks 23 (32.9%); significant either way 70 (11.6%); all three within tolerance: False; median of the mean permuted difference 0.67
G-a pairs b, later x, height ii: 604 of 605 pairs testable; depression 31 (5.1%), facilitation 42 (6.9%); significant pairs on different shanks 25 (34.2%); significant either way 73 (12.1%); all three within tolerance: False; median of the mean permuted difference 0.49
G-a pairs b, later y, height i : 604 of 605 pairs testable; depression 31 (5.1%), facilitation 49 (8.1%); significant pairs on different shanks 27 (33.8%); significant either way 80 (13.2%); all three within tolerance: False; median of the mean permuted difference 0.14
G-a pairs b, later y, height ii: 604 of 605 pairs testable; depression 28 (4.6%), facilitation 60 (9.9%); significant pairs on different shanks 25 (28.4%); significant either way 88 (14.5%); all three within tolerance: False; median of the mean permuted difference 0.26

null check, labels assigned at random: 16 of 121 pairs significant at P < 0.10 (13.2%, 95% interval 7.8 to 20.6%; expected 12.1 +- 3.3; binomial test against 10%: p = 0.23)

s = S["EE188"]
ok197 = ok186 = False
for ch in ("x", "y"):
    for how in ("i", "ii"):
        r1 = F6.get(("EE188", idx[197], idx[201], ch, how))
        r2 = F6.get(("EE188", idx[186], idx[201], ch, how))
        for name, r in (("197-201", r1), ("186-201", r2)):
            if r is None:
                print(f"{name} later {ch} height {how}: not a connection under G-a B-a W-a")
                continue
            print(f"{name} later {ch} height {how}: first z = {r[0]:.2f} ({r[4]} spikes), later z = {r[1]:.2f} ({r[5]} spikes), "
                  f"difference {r[2]:.2f}, mean permuted difference {r[7]:.2f}, P = {r[3]:.3f}, "
                  f"{'depression' if r[6] > 0 else 'facilitation'} side")
        ok197 |= r1 is not None and r1[6] > 0 and r1[3] < 0.10
        ok186 |= r2 is not None and r2[6] < 0 and r2[3] < 0.10
E, I = conns("EE188", "B-a", "W-a")
print("connections under G-b B-a W-a: 197-201", E[idx[197], idx[201]], ", 186-201", E[idx[186], idx[201]])
print("CellsInfo types: 197", s.celltype_file[idx[197]], "186", s.celltype_file[idx[186]], "201", s.celltype_file[idx[201]])
print(f"197-201 depressing at P < 0.10 under one reading: {ok197}; 186-201 facilitating at P < 0.10 under one reading: {ok186}")

197-201 later x height i: first z = 10.23 (2311 spikes), later z = 2.10 (252 spikes), difference 8.13, mean permuted difference 7.69, P = 0.910, depression side
186-201 later x height i: first z = 1.42 (2134 spikes), later z = 5.76 (571 spikes), difference -4.35, mean permuted difference 2.51, P = 0.012, facilitation side
197-201 later x height ii: first z = 8.42 (2311 spikes), later z = 3.12 (252 spikes), difference 5.30, mean permuted difference 5.73, P = 0.707, facilitation side
186-201 later x height ii: first z = 3.55 (2134 spikes), later z = 3.45 (571 spikes), difference 0.10, mean permuted difference 1.85, P = 0.160, facilitation side
197-201 later y height i: first z = 8.44 (2311 spikes), later z = 3.73 (464 spikes), difference 4.71, mean permuted difference 4.88, P = 0.878, facilitation side
186-201 later y height i: first z = 1.05 (2134 spikes), later z = 9.26 (2389 spikes), difference -8.21, mean permuted difference 0.76, P = 0.004, facilitation side
197-201 later y height ii: first z = 9.10 (2311 spikes), later z = 4.03 (464 spikes), difference 5.07, mean permuted difference 4.99, P = 0.914, depression side
186-201 later y height ii: first z = 4.01 (2134 spikes), later z = 7.92 (2389 spikes), difference -3.91, mean permuted difference 0.70, P = 0.004, facilitation side
connections under G-b B-a W-a: 197-201 True , 186-201 True
CellsInfo types: 197 1 186 1 201 2
197-201 depressing at P < 0.10 under one reading: False; 186-201 facilitating at P < 0.10 under one reading: True

def satellite_curve(s, pres, tgt, w):
    tl = [s.spikes_of(p) for p in pres]
    t = np.concatenate(tl)
    c = np.concatenate([np.full(len(x), n) for n, x in enumerate(tl)])
    o = np.argsort(t)
    t, c = t[o], c[o]
    tt = s.spikes_of(tgt)
    nxt = np.append(np.diff(t), np.inf)
    a = np.searchsorted(t, t - w, side="right")
    k = np.array([len(np.unique(c[a[m]:m + 1])) for m in range(len(t))])
    i0 = np.searchsorted(tt, t + 0.001, side="left")
    i1 = np.searchsorted(tt, t + 0.004, side="right")
    hit = i1 > i0
    clean = nxt > 0.004
    prev_other = np.array([np.any(c[a[m]:m] != c[m]) for m in range(len(t))])
    p = {}
    for kk in (1, 2, 3):
        m = clean & (k == kk) & ((~prev_other) if kk == 1 else True)
        p[kk] = (hit[m].mean() if m.sum() >= 20 else np.nan, int(m.sum()))
    return p


def satellites(gs):
    out = []
    for k, s in S.items():
        E, I = conns(k, "B-a", "W-a", gs)
        pyr, inn = TYPES[(k, "B-a", "W-a", "C-a", gs)]
        for tgt in np.flatnonzero(inn):
            pres = np.flatnonzero(E[:, tgt] & pyr)
            if len(pres) >= 2:
                out.append((k, tgt, pres))
    return out


for gs in GS:
    sat = satellites(gs)
    print(f"{gs} B-a W-a: satellites {len(sat)} (paper 14); within 14 +- 1.4: {abs(len(sat) - 14) <= 1.4}; presynaptic cells "
          f"per satellite {[len(x[2]) for x in sat]}")
SAT = satellites("G-b")
print("supralinearity, G-b satellites:")
for w in (0.002, 0.005, 0.01, 0.02):
    rows = []
    for k, tgt, pres in SAT:
        p = satellite_curve(S[k], pres, tgt, w)
        rows.append([p[1][0], p[2][0], p[3][0], p[2][1], p[3][1]])
    R = np.array(rows, float)
    p1, p2, p3 = R[:, 0], R[:, 1], R[:, 2]
    for kk, pk in ((2, p2), (3, p3)):
        ok = ~np.isnan(pk) & ~np.isnan(p1)
        lin_a = kk * p1[ok]
        lin_b = 1 - (1 - p1[ok]) ** kk
        ex_a = pk[ok] - lin_a
        ex_b = pk[ok] - lin_b
        print(f"window {w * 1000:4.0f} ms, k = {kk}: {ok.sum()} satellites with values; mean p{kk} {np.mean(pk[ok]):.3f}, "
              f"mean p1 {np.mean(p1[ok]):.3f}; above k*p1 in {np.sum(ex_a > 0)}, mean excess {np.mean(ex_a):.3f} "
              f"(s.d. {np.std(ex_a, ddof=1) if ok.sum() > 1 else np.nan:.3f}); above 1-(1-p1)^k in {np.sum(ex_b > 0)}, "
              f"mean excess {np.mean(ex_b):.3f} (s.d. {np.std(ex_b, ddof=1) if ok.sum() > 1 else np.nan:.3f})")

G-b B-a W-a: satellites 16 (paper 14); within 14 +- 1.4: False; presynaptic cells per satellite [5, 5, 3, 7, 3, 3, 9, 2, 5, 4, 2, 8, 2, 6, 2, 5]
G-a B-a W-a: satellites 15 (paper 14); within 14 +- 1.4: True; presynaptic cells per satellite [7, 5, 3, 6, 2, 4, 3, 10, 3, 3, 9, 2, 6, 2, 5]
supralinearity, G-b satellites:

window    2 ms, k = 2: 15 satellites with values; mean p2 0.072, mean p1 0.055; above k*p1 in 3, mean excess -0.037 (s.d. 0.028); above 1-(1-p1)^k in 3, mean excess -0.033 (s.d. 0.027)
window    2 ms, k = 3: 0 satellites with values; mean p3 nan, mean p1 nan; above k*p1 in 0, mean excess nan (s.d. nan); above 1-(1-p1)^k in 0, mean excess nan (s.d. nan)

/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/fromnumeric.py:3862: RuntimeWarning: Mean of empty slice
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/numpy/_core/_methods.py:142: RuntimeWarning: invalid value encountered in scalar divide
  ret = ret.dtype.type(ret / rcount)

window    5 ms, k = 2: 16 satellites with values; mean p2 0.064, mean p1 0.054; above k*p1 in 0, mean excess -0.044 (s.d. 0.025); above 1-(1-p1)^k in 0, mean excess -0.040 (s.d. 0.023)
window    5 ms, k = 3: 2 satellites with values; mean p3 0.119, mean p1 0.090; above k*p1 in 0, mean excess -0.152 (s.d. 0.004); above 1-(1-p1)^k in 0, mean excess -0.128 (s.d. 0.005)

window   10 ms, k = 2: 16 satellites with values; mean p2 0.063, mean p1 0.054; above k*p1 in 0, mean excess -0.045 (s.d. 0.024); above 1-(1-p1)^k in 0, mean excess -0.041 (s.d. 0.021)
window   10 ms, k = 3: 5 satellites with values; mean p3 0.126, mean p1 0.083; above k*p1 in 0, mean excess -0.123 (s.d. 0.049); above 1-(1-p1)^k in 0, mean excess -0.101 (s.d. 0.036)

window   20 ms, k = 2: 16 satellites with values; mean p2 0.062, mean p1 0.053; above k*p1 in 0, mean excess -0.045 (s.d. 0.024); above 1-(1-p1)^k in 0, mean excess -0.041 (s.d. 0.021)
window   20 ms, k = 3: 8 satellites with values; mean p3 0.085, mean p1 0.066; above k*p1 in 0, mean excess -0.112 (s.d. 0.052); above 1-(1-p1)^k in 0, mean excess -0.097 (s.d. 0.041)

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

notebook run time 20.5 min

Claim	Tested here
unit-count	tested
task-performance	tested on EE188 only (sensor channels of its tar); the trial table holds correct trials only
fig2c-synchrony-by-layer	tested
fig2-population-uniform	tested (layer 2/3 session only; the layer 5 session is not in the data set)
lifetime-width	tested
trajectory-divergence	tested
side-arm-selectivity-by-layer	tested
central-arm-selectivity	tested
odor-selective-fraction	tested
odor-unrelated-to-stem	tested
retrospective-coding-rare	tested
connection-counts	tested
cell-type-fractions	tested
pyr-pyr-connections	tested
fig3b-session-connections	tested (layer 2/3 session only)
fig3d-hub	tested (layer 2/3 session only)
fig3f-distance	tested
nearby-pyr-int-probability	tested
connected-cells-trajectory-specific	tested
no-spatial-clustering	tested
position-dependent-transmission	tested
fig5a-session-selective-connections	tested
fig4-pair-156-136	tested
fig5b-example-pairs	tested
thinning-preserves-dependence	tested
fig6-depression-facilitation	tested
fig6b-convergent-examples	tested
fig7-supralinear-coincidence	tested
speed-not-correlated	tested
population-vector-width	tested

Claim	Outcome
unit-count	partially reproduced
task-performance	partially reproduced (EE188 only)
fig2c-synchrony-by-layer	partially reproduced
fig2-population-uniform	not reproduced
lifetime-width	reproduced
trajectory-divergence	partially reproduced
side-arm-selectivity-by-layer	not reproduced
central-arm-selectivity	not reproduced
odor-selective-fraction	not reproduced
odor-unrelated-to-stem	not reproduced (resting on an assumption: Spearman correlation)
retrospective-coding-rare	reproduced
connection-counts	not reproduced
cell-type-fractions	not reproduced
pyr-pyr-connections	not reproduced
fig3b-session-connections	not reproduced
fig3d-hub	not reproduced
fig3f-distance	partially reproduced
nearby-pyr-int-probability	not reproduced
connected-cells-trajectory-specific	partially reproduced
no-spatial-clustering	reproduced
position-dependent-transmission	not reproduced
fig5a-session-selective-connections	not reproduced
fig4-pair-156-136	partially reproduced
fig5b-example-pairs	not reproduced
thinning-preserves-dependence	not reproduced
fig6-depression-facilitation	not reproduced
fig6b-convergent-examples	partially reproduced
fig7-supralinear-coincidence	not reproduced
speed-not-correlated	reproduced
population-vector-width	reproduced (resting on an assumption: PSTHs normalised per cell)

Fujisawa, Amarasingham, Harrison and Buzsaki (2008) on CRCNS pfc-2¶

Data check¶

unit-count¶

task-performance¶

Firing rates along the maze¶

fig2c-synchrony-by-layer¶

fig2-population-uniform¶

lifetime-width¶

trajectory-divergence¶

side-arm-selectivity-by-layer¶

central-arm-selectivity¶

odor-selective-fraction¶

odor-unrelated-to-stem¶

retrospective-coding-rare¶

speed-not-correlated¶

population-vector-width¶

Monosynaptic connections¶

connection-counts¶

cell-type-fractions¶

pyr-pyr-connections¶

fig3b-session-connections¶

fig3d-hub¶

fig3f-distance¶

nearby-pyr-int-probability¶

connected-cells-trajectory-specific¶

no-spatial-clustering¶

Position dependence of spike transmission¶

position-dependent-transmission¶

fig5a-session-selective-connections¶

fig4-pair-156-136¶

fig5b-example-pairs¶

thinning-preserves-dependence¶

Spike order and convergence¶

fig6-depression-facilitation¶

fig6b-convergent-examples¶

fig7-supralinear-coincidence¶

Summary¶