import time
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from dandi.dandiapi import DandiAPIClient
from scipy import signal, stats
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

import dg

pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)

notebook_started = started = time.time()
S = dg.load_all(workers=12)
print(f"loaded {len(S)} files from dandiset {dg.DANDISET_ID} version {dg.DANDISET_VERSION} in {time.time() - started:.0f} s")

loaded 101 files from dandiset 000003 version 0.260218.2052 in 225 s

# Versions: compare the assets of every published version with the pinned one.
with DandiAPIClient() as client:
    versions = [v.identifier for v in client.get_dandiset(dg.DANDISET_ID).get_versions() if v.identifier != "draft"]
    listing = {}
    for v in versions:
        ds = client.get_dandiset(dg.DANDISET_ID, v)
        listing[v] = {a.path: (a.size, a.get_raw_metadata().get("digest", {}).get("dandi:dandi-etag")) for a in ds.get_assets()}
ref = listing[dg.DANDISET_VERSION]
for v in versions:
    same = listing[v] == ref
    print(f"version {v}: {len(listing[v])} assets; same paths, sizes and etags as {dg.DANDISET_VERSION}: {same}")

version 0.210812.1448: 101 assets; same paths, sizes and etags as 0.260218.2052: True
version 0.230629.1955: 101 assets; same paths, sizes and etags as 0.260218.2052: True
version 0.250624.0409: 101 assets; same paths, sizes and etags as 0.260218.2052: True
version 0.260218.2052: 101 assets; same paths, sizes and etags as 0.260218.2052: True

info = pd.DataFrame([{"session": s["session"], "subject": s["subject"], "genotype": s["genotype"], "age": s["age"][:9],
                      "format": s["format"], "units": len(s["units_all"]), "kept": len(s["units"]),
                      "electrodes": s["n_electrodes"], "raw": s["has_raw"],
                      "span_h": (s["spike_time_range"][1] - s["spike_time_range"][0]) / 3600} for s in S])
print("Files:", len(info), "; subjects:", info.subject.nunique())
print("\nSubjects by genotype (paper: 7 SST-Cre, 5 DRD2-Cre::Ai35, 3 POMC-Cre::Ai35, 1 wild type):")
print(info.groupby("genotype").subject.nunique().to_string())
print("\nFiles per subject:")
print(info.groupby(["subject", "genotype"]).size().to_string())
print("\nFile formats (units table columns):", info.format.value_counts().to_dict())
print("Files with the raw wide-band signal (used for unit sites and waveforms):", int(info.raw.sum()))
print("Recording days (subject and date):", info.assign(day=info.session.str[-6:]).groupby(["subject", "day"]).ngroups)
print("Recording span per file (h): median %.1f, range %.1f to %.1f" % (info.span_h.median(), info.span_h.min(), info.span_h.max()))

Files: 101 ; subjects: 16

Subjects by genotype (paper: 7 SST-Cre, 5 DRD2-Cre::Ai35, 3 POMC-Cre::Ai35, 1 wild type):
genotype
Drd2-Cre::Ai35    5
POMC-Cre::Ai35    1
POMC-Cre::Arch    4
SST-Cre           6

Files per subject:
subject      genotype      
YutaMouse20  SST-Cre            5
YutaMouse23  SST-Cre            3
YutaMouse33  POMC-Cre::Arch     9
YutaMouse37  SST-Cre            5
YutaMouse38  SST-Cre            1
YutaMouse39  SST-Cre            7
YutaMouse40  SST-Cre            3
YutaMouse41  POMC-Cre::Arch    10
YutaMouse42  POMC-Cre::Arch    15
YutaMouse44  POMC-Cre::Arch     3
YutaMouse45  Drd2-Cre::Ai35     5
YutaMouse51  Drd2-Cre::Ai35     4
YutaMouse54  Drd2-Cre::Ai35     5
YutaMouse55  Drd2-Cre::Ai35     9
YutaMouse56  Drd2-Cre::Ai35     9
YutaMouse57  POMC-Cre::Ai35     8

File formats (units table columns): {'standard': 83, 'electrode_group': 18}
Files with the raw wide-band signal (used for unit sites and waveforms): 17
Recording days (subject and date): 94
Recording span per file (h): median 9.7, range 4.7 to 17.0

U = pd.concat([s["units_all"].assign(session=s["session"], subject=s["subject"]) for s in S], ignore_index=True)
print("Units by cell_type label (all files):")
print(U.cell_type.value_counts().to_string())
print("\nglobal_id: %d units without one (label 'missing' has %d); duplicated values: %d" % (
    U.global_id.isna().sum(), (U.cell_type == "missing").sum(), U.global_id.dropna().duplicated().sum()))
print("\nKept units by type and subject:")
print(pd.crosstab(U[U.type != "other"].subject, U[U.type != "other"].type, margins=True).to_string())

Units by cell_type label (all files):
cell_type
unknown                                                        1104
pyramidal cell                                                  649
granule cell                                                    609
mossy cell                                                      568
narrow waveform cell                                            537
missing                                                         237
wide waveform cell (narrower, exclude opto tagged SST cell)     191
wide waveform cell (wider)                                      184
positive waveform unit (bursty)                                  60
positive waveform unit (non-bursty)                              31
optogenetically tagged SST cell                                  21
positive negative waveform unit                                  18

global_id: 237 units without one (label 'missing' has 237); duplicated values: 0

Kept units by type and subject:
type         CA3   GC   MC  SST   nI   wI   All
subject                                        
YutaMouse20    0    3   26    6   16    7    58
YutaMouse23    0   27   12    2   12    6    59
YutaMouse33    0   35   74    0   61   66   236
YutaMouse37    0   22    9    0    5    6    42
YutaMouse38    0    5    5    2    3    3    18
YutaMouse39   46   40   30   11   23   12   162
YutaMouse40    0   23    0    0    4    3    30
YutaMouse41  104  130   86    0   58   72   450
YutaMouse42  182   65  146    0  126   51   570
YutaMouse44    0    8    1    0    6    3    18
YutaMouse45    0   24    7    0    4    2    37
YutaMouse51    0   32    3    0    7    9    51
YutaMouse54   68   45   55    0   41   25   234
YutaMouse55  102   37    6    0   44   21   210
YutaMouse56   66   87   43    0   47   34   277
YutaMouse57   81   26   65    0   80   55   307
All          649  609  568   21  537  375  2759

rep = pd.DataFrame([{"session": s["session"], **r} for s in S for r in s["repairs"]])
print("Files where the first unit's spike train continues at the start of the second unit's:", len(rep))
print(rep.to_string(index=False))
print("\nUnits with decreasing spike times as stored:", sum(len(s["unsorted_units_stored"]) for s in S),
      "in", sum(bool(s["unsorted_units_stored"]) for s in S), "files; after the repair:",
      sum(len(s["unsorted_units_repaired"]) for s in S))
lab = {s["session"]: s["units_all"].set_index("unit").cell_type for s in S}
print("\nLabels of the units affected (unit 0 gains, unit 1 loses the moved spikes):")
print(pd.Series([f"{lab[r.session][0]} / {lab[r.session][1]}" for r in rep.itertuples()]).value_counts().to_string())
print("\nLargest distance of a spike time from the 20 kHz sample grid (samples):", max(s["spike_grid_error"] for s in S))

Files where the first unit's spike train continues at the start of the second unit's: 24
           session  unit  moved_spikes   gap_s
YutaMouse20-140325     0        196608 0.00230
YutaMouse20-140327     0        327680 0.02665
YutaMouse20-140328     0        131072 0.00370
YutaMouse23-140805     0        262144 0.20970
YutaMouse33-150223     0        327680 0.02715
YutaMouse37-150609     0         65536 0.46670
YutaMouse37-150610     0        131072 0.11735
YutaMouse37-150611     0        131072 0.20955
YutaMouse39-150729     0         65536 0.29530
YutaMouse40-150812     0        196608 0.08245
YutaMouse41-150831     0        131072 0.26830
YutaMouse41-150901     0        458752 0.00975
YutaMouse41-150903     0         65536 0.01210
YutaMouse45-160302     0        393216 0.07310
YutaMouse51-160516     0         65536 0.02290
YutaMouse54-160701     0         65536 0.50395
YutaMouse55-160910     0         65536 0.19400
YutaMouse55b160907     0         65536 0.02895
YutaMouse56-160914     0         65536 0.02810
YutaMouse56-160916     0         65536 0.30635
YutaMouse56b160915     0         65536 0.07195
YutaMouse57-161010     0        196608 0.01940
YutaMouse57-161011     0         65536 0.52640
YutaMouse57-161014     0         65536 0.00790

Units with decreasing spike times as stored: 24 in 24 files; after the repair: 0

Labels of the units affected (unit 0 gains, unit 1 loses the moved spikes):
unknown / unknown                                                                     6
unknown / granule cell                                                                3
narrow waveform cell / pyramidal cell                                                 3
narrow waveform cell / unknown                                                        2
wide waveform cell (narrower, exclude opto tagged SST cell) / narrow waveform cell    2
pyramidal cell / pyramidal cell                                                       2
positive waveform unit (non-bursty) / positive waveform unit (non-bursty)             1
missing / missing                                                                     1
wide waveform cell (wider) / mossy cell                                               1
unknown / pyramidal cell                                                              1
positive waveform unit (bursty) / narrow waveform cell                                1
narrow waveform cell / granule cell                                                   1

Largest distance of a spike time from the 20 kHz sample grid (samples): 1.1920928955078125e-07

st_all = pd.concat([s["states"].assign(session=s["session"]) for s in S], ignore_index=True)
st_all["dur"] = st_all.stop_time - st_all.start_time
print("State periods by label: count, total hours, median minutes")
print(st_all.groupby("label").dur.agg(["size", lambda d: d.sum() / 3600, lambda d: d.median() / 60]).round(2).to_string())
print("Files without a states table:", sum(len(s["states"]) == 0 for s in S))
def overlaps(g):
    g = g.sort_values("start_time")
    return int((g.start_time.values[1:] < g.stop_time.values[:-1] - 1e-3).sum())
print("Overlapping periods within a file:", sum(overlaps(g) for _, g in st_all.groupby("session")))

# Clock check: rates of GC and CA3 units in NREM against awake periods, with the true labels and with
# the labels shifted by 30 min. If the states share the spike clock, the contrast should shrink.
def nrem_awake_contrast(shift):
    out = []
    for s in S:
        st = s["states"].copy()
        st[["start_time", "stop_time"]] += shift
        a, n = st[st.label == "awake"], st[st.label == "nrem"]
        da, dn = (a.stop_time - a.start_time).sum(), (n.stop_time - n.start_time).sum()
        if da < 300 or dn < 300:
            continue
        for u in s["units"][s["units"].type.isin(["GC", "CA3"])].unit:
            t = dg.times(s, u)
            ra = dg.interval_mask_count(t, a.start_time, a.stop_time) / da
            rn = dg.interval_mask_count(t, n.start_time, n.stop_time) / dn
            if ra > 0 and rn > 0:
                out.append(np.log2(rn / ra))
    return np.median(out), len(out)
for shift in (0, 1800, -1800):
    m, n = nrem_awake_contrast(shift)
    print(f"labels shifted by {shift:+5d} s: median log2(NREM rate / awake rate) of GC and CA3 units = {m:.2f} (n = {n})")

State periods by label: count, total hours, median minutes

         size  <lambda_0>  <lambda_1>
label                                
awake     566      278.34        8.34
nrem     1100      193.51        8.89
quiet     413       20.38        1.83
rem       642       16.29        1.38
transit   643        2.73        0.22
Files without a states table: 0

Overlapping periods within a file: 0

labels shifted by    +0 s: median log2(NREM rate / awake rate) of GC and CA3 units = 1.27 (n = 1229)

labels shifted by +1800 s: median log2(NREM rate / awake rate) of GC and CA3 units = 0.30 (n = 1229)

labels shifted by -1800 s: median log2(NREM rate / awake rate) of GC and CA3 units = 0.33 (n = 1229)

# Positions. Open fields: normalised 0..1, scaled by the stored conversion (0.46 m = the 46 cm arena).
rows = []
for s in S:
    for k, p in s["positions"].items():
        xy = p["xy"]
        ok = np.all(np.isfinite(xy), axis=1)
        rows.append({"session": s["session"], "series": k.replace("_twhl", ""), "n": len(xy), "nan_frac": 1 - ok.mean(),
                     "dt_ms": 1000 * np.median(np.diff(p["t"])), "c0_min": np.nanmin(xy[:, 0]), "c0_max": np.nanmax(xy[:, 0]),
                     "c1_min": np.nanmin(xy[:, 1]), "c1_max": np.nanmax(xy[:, 1]), "conversion": p["conversion"], "unit": p["unit"]})
pos = pd.DataFrame(rows)
print(pos.groupby("series").agg(files=("session", "size"), dt_ms=("dt_ms", "median"), nan_frac=("nan_frac", "max"),
                                c0_min=("c0_min", "min"), c0_max=("c0_max", "max"), c1_min=("c1_min", "min"),
                                c1_max=("c1_max", "max"), conversion=("conversion", "first"), unit=("unit", "first")).round(3).to_string())
print("\nSessions with open fields A (Old), B (New) and A' (Oldlast):",
      sum(all(f"{dg.OPEN_FIELDS[k]}_twhl_norm" in s["positions"] for k in ("A", "B", "A2")) for s in S))
print("Sessions with T-maze trials:", sum(s["trials"] is not None and "EightMazePosition_twhl_linearized" in s["positions"] for s in S))

                                    files  dt_ms  nan_frac  c0_min  c0_max  c1_min  c1_max  conversion    unit
series                                                                                                        
EightMazePosition_linearized           45   25.6     0.587   0.000   2.000  -3.142   3.142       0.325  meters
EightMazePosition_norm                 45   25.6     0.587  -2.637   1.587  -1.602   1.329       0.325  meters
OpenFieldPosition_ExtraLarge_norm      11   25.6     0.000   0.000   1.000   0.000   1.000         NaN  meters
OpenFieldPosition_New_Curtain_norm     21   25.6     0.000   0.000   1.000   0.000   1.000       0.460  meters
OpenFieldPosition_New_norm             43   25.6     0.000   0.000   1.000   0.000   1.000       0.460  meters
OpenFieldPosition_Old_Curtain_norm     21   25.6     0.000   0.000   1.000   0.000   1.000       0.460  meters
OpenFieldPosition_Old_norm             43   25.6     0.000   0.000   1.000   0.000   1.000       0.460  meters
OpenFieldPosition_Oldlast_norm         35   25.6     0.000   0.000   1.000   0.000   1.000       0.460  meters

Sessions with open fields A (Old), B (New) and A' (Oldlast): 35
Sessions with T-maze trials: 42

# T-maze: what the two linearised columns are. Column 0 is 0..2 on the central arm and 0 elsewhere;
# column 1 is an angle (-pi..pi) on the side arms and 0 on the central arm.
TM = [s for s in S if s["trials"] is not None and "EightMazePosition_twhl_linearized" in s["positions"]]
both = zero = n = 0
for s in TM:
    a = s["positions"]["EightMazePosition_twhl_linearized"]["xy"]
    a = a[np.all(np.isfinite(a), axis=1)]
    both += int(((a[:, 0] > 0) & (a[:, 1] != 0)).sum()); zero += int(((a[:, 0] == 0) & (a[:, 1] == 0)).sum()); n += len(a)
print(f"T-maze samples: {n}; with both columns non-zero: {both}; with both zero: {zero}")
print("Central arm length (2 x 0.325 m): %.0f cm; side arm at |angle| = pi: %.0f cm (Methods: 65 and 102 cm)" % (2 * dg.RADIUS_CM, np.pi * dg.RADIUS_CM))
runs_all = []
for s in TM:
    for a, b, side, cond in dg.tmaze_runs(s):
        runs_all.append({"session": s["session"], "subject": s["subject"], "dur": b - a, "side": side, "cond": cond})
runs_all = pd.DataFrame(runs_all)
trials_all = sum(len(s["trials"]) for s in TM)
nonstim = sum(int((s["trials"].stim_run == 0).sum()) for s in TM)
print(f"Trials: {trials_all}; without light: {nonstim}; outbound runs found: {len(runs_all)}")
print("Side arm reached against the trial's condition:")
print(pd.crosstab(runs_all.cond, runs_all.side).to_string())
print("Run duration (s): median %.1f, 5th-95th percentile %.1f to %.1f" % (runs_all.dur.median(), *np.percentile(runs_all.dur, [5, 95])))
print("Runs per session: median %d, range %d to %d" % (runs_all.groupby("session").size().median(), runs_all.groupby("session").size().min(), runs_all.groupby("session").size().max()))

T-maze samples: 7182639; with both columns non-zero: 0; with both zero: 0
Central arm length (2 x 0.325 m): 65 cm; side arm at |angle| = pi: 102 cm (Methods: 65 and 102 cm)

Trials: 2266; without light: 1841; outbound runs found: 1701
Side arm reached against the trial's condition:
side         1    2
cond               
run_left   877    7
run_right   28  789
Run duration (s): median 8.0, 5th-95th percentile 4.5 to 71.7
Runs per session: median 40, range 2 to 84

# Clock check for positions: spatial information of MC and CA3 units on the T-maze with true positions,
# and with the spike train shifted by 20 s (circularly within the session).
def info_with_shift(shift):
    vals = []
    for s in TM:
        units = s["units"][s["units"].type.isin(["MC", "CA3"])].unit.tolist()
        if not units:
            continue
        saved = {u: s["spikes"][u] for u in units}
        span = int((s["spike_time_range"][1] - s["spike_time_range"][0]) * dg.SPIKE_HZ)
        for u in units:
            s["spikes"][u] = np.sort((saved[u] + int(shift * dg.SPIKE_HZ)) % span).astype(np.int32)
        occ, maps, _ = dg.tmaze_maps(s, units)
        for u in units:
            vals.append(dg.spatial_info(np.concatenate(maps[u]), np.concatenate(occ)))
        s["spikes"].update(saved)
    return np.nanmedian(vals), len(vals)
for shift in (0, 20):
    m, n = info_with_shift(shift)
    print(f"spikes shifted by {shift:2d} s: median spatial information of MC and CA3 units = {m:.3f} bits/s (n = {n})")

spikes shifted by  0 s: median spatial information of MC and CA3 units = 0.440 bits/s (n = 644)

spikes shifted by 20 s: median spatial information of MC and CA3 units = 0.089 bits/s (n = 644)

# DS2 events and light pulses.
ds2n = pd.DataFrame([{"session": s["session"], "series": ",".join(s["ds2_series"]), "n_ds2": len(s["ds2"]),
                      "per_min_nrem": len(s["ds2"]) / max(1e-9, (s["states"][s["states"].label == "nrem"].eval("stop_time - start_time").sum() / 60))}
                     for s in S])
print("DS2 peaks per file: median %d, range %d to %d; files without DS2 events: %d" % (
    ds2n.n_ds2.median(), ds2n.n_ds2.min(), ds2n.n_ds2.max(), (ds2n.n_ds2 == 0).sum()))
print("Files with more than one DS2 series (pooled here):", ds2n[ds2n.series.str.contains(",")][["session", "series"]].values.tolist())
pul = pd.DataFrame([{"session": s["session"], "genotype": s["genotype"], "n": len(s["pulse20"]),
                     "ipi_ms_median": 1000 * np.median(np.diff(s["pulse20"])), "ipi_ms_min": 1000 * np.min(np.diff(s["pulse20"]))}
                    for s in S if len(s["pulse20"]) > 1])
print("\nFiles with 20 ms light pulses (PulseStim_5V_20ms_LD12), by genotype:", pul.genotype.value_counts().to_dict())
print("Pulses per file: %d to %d (Methods: 60,000 to 100,000 repetitions)" % (pul.n.min(), pul.n.max()))
print("Interval between pulses after sorting (ms): median %.1f, smallest %.1f" % (pul.ipi_ms_median.median(), pul.ipi_ms_min.min()))
print("Files whose stored pulse times decrease or repeat (session, decreases, duplicates):",
      [(s["session"], s["pulse20_decreasing"], s["pulse20_duplicates"]) for s in S if s["pulse20_decreasing"] or s["pulse20_duplicates"]])

DS2 peaks per file: median 4690, range 44 to 47574; files without DS2 events: 0
Files with more than one DS2 series (pooled here): [['YutaMouse20-140325', 'DS2.ch11,DS2.ch33'], ['YutaMouse55-160906', 'DS2.ch48,DS2.ch60']]

Files with 20 ms light pulses (PulseStim_5V_20ms_LD12), by genotype: {'Drd2-Cre::Ai35': 26, 'POMC-Cre::Ai35': 8}
Pulses per file: 12478 to 97191 (Methods: 60,000 to 100,000 repetitions)
Interval between pulses after sorting (ms): median 200.0, smallest 199.2
Files whose stored pulse times decrease or repeat (session, decreases, duplicates): [('YutaMouse57-161006', 1, 24286)]

# Recordings deposited twice: sessions of one mouse on the same date (the 'b' files) against each other.
from itertools import combinations
def coincidence(a, b, tol=0.0005):
    j = np.clip(np.searchsorted(b, a), 1, len(b) - 1)
    return float((np.minimum(np.abs(b[j] - a), np.abs(b[j - 1] - a)) < tol).mean())
by_day = {}
for s in S:
    by_day.setdefault((s["subject"], s["session"][-6:]), []).append(s)
rows = []
for (subj, day), group in by_day.items():
    for s1, s2 in combinations(group, 2):
        best = 0.0
        for u in s1["spikes"]:
            a = dg.times(s1, u)[:20000]
            for v in s2["spikes"]:
                b = dg.times(s2, v)
                if len(a) > 100 and len(b) > 100:
                    best = max(best, coincidence(a, b))
        n1, n2 = len(s1["positions"]), len(s2["positions"])
        rows.append({"pair": f"{s1['session']} / {s2['session']}", "units": f"{len(s1['spikes'])} / {len(s2['spikes'])}",
                     "max_coincident_spike_fraction": round(best, 3),
                     "same_states": s1["states"][["start_time", "stop_time"]].round(0).equals(s2["states"][["start_time", "stop_time"]].round(0)),
                     "same_unit_count": len(s1["units_all"]) == len(s2["units_all"])})
dups = pd.DataFrame(rows)
print(dups.to_string(index=False))
# Across days: identical spike trains between any two files of a mouse (first 50 spikes of each unit).
sig = {}
for s in S:
    for u in s["spikes"]:
        key = (s["subject"], tuple(s["spikes"][u][:50]))
        sig.setdefault(key, []).append(s["session"])
print("Units whose first 50 spike times appear in more than one file:", sum(len(v) > 1 for v in sig.values()))

                                   pair   units  max_coincident_spike_fraction  same_states  same_unit_count
YutaMouse51-160516 / YutaMouse51b160516 17 / 13                          0.014        False            False
YutaMouse54-160701 / YutaMouse54b160701 40 / 52                          0.032        False            False
YutaMouse54-160702 / YutaMouse54b160702 67 / 69                          0.033        False            False
YutaMouse55-160907 / YutaMouse55b160907 19 / 13                          0.030        False            False
YutaMouse56-160913 / YutaMouse56b160913 30 / 12                          0.043        False            False
YutaMouse56-160914 / YutaMouse56b160914  3 / 36                          0.035        False            False
YutaMouse56-160915 / YutaMouse56b160915 29 / 26                          0.025        False            False
Units whose first 50 spike times appear in more than one file: 0

# Theta-reference LFP channel and the phase data used for theta-modulation and gamma-modulation.
ph = pd.DataFrame([{"session": s["session"], "awake_s": s["phases"]["seconds"] if s["phases"] else 0,
                    "theta_above": s["phases"]["theta_fraction_above"] if s["phases"] else np.nan,
                    "gamma_above": s["phases"]["gamma_fraction_above"] if s["phases"] else np.nan} for s in S])
print("Files with one theta-reference LFP channel and awake periods used:", int((ph.awake_s > 0).sum()), "of", len(S))
print("Awake LFP used per file (s): median %.0f, range %.0f to %.0f" % (ph[ph.awake_s > 0].awake_s.median(), ph[ph.awake_s > 0].awake_s.min(), ph[ph.awake_s > 0].awake_s.max()))
print("Fraction of that time above the amplitude criterion: theta %.2f, gamma %.2f" % (ph.theta_above.median(), ph.gamma_above.median()))
print("Files without:", ph[ph.awake_s == 0].session.tolist())

Files with one theta-reference LFP channel and awake periods used: 98 of 101
Awake LFP used per file (s): median 3600, range 297 to 3600
Fraction of that time above the amplitude criterion: theta 0.31, gamma 0.18
Files without: ['YutaMouse33-150218', 'YutaMouse55-160901', 'YutaMouse55-160902']

# Files that give each unit's shank: DS2 amplitude per LFP channel along each shank (the LFP over
# the 2 LFP chunks with the most DS2 peaks), and where the units sit.
SITE = [s for s in S if s["unit_site"] is not None]
print("Files with unit sites:", len(SITE), "; mice:", sorted({s["subject"] for s in SITE}))
print("DS2 events used per file: median %d, range %d to %d; LFP window %.0f s" % (
    np.median([s["unit_site"]["n_events"] for s in SITE]), min(s["unit_site"]["n_events"] for s in SITE),
    max(s["unit_site"]["n_events"] for s in SITE), np.median([s["unit_site"]["window_s"] for s in SITE])))
rev = []
for s in SITE:
    us = s["unit_site"]
    g = us["groups"][us["region"]]
    for sh in np.unique(g):
        a = us["channel_amp"][g == sh]
        rev.append({"session": s["session"], "shank": sh, "min_uV": 1e6 * a.min(), "max_uV": 1e6 * a.max(), "reverses": bool(a.min() < 0 < a.max())})
rev = pd.DataFrame(rev)
print("Shanks: %d; with a DS2 polarity reversal along the shank: %d" % (len(rev), rev.reverses.sum()))
print("DS2 amplitude range across channels (uV): %.0f to %.0f" % (rev.min_uV.min(), rev.max_uV.max()))
wch = pd.DataFrame([{"session": s["session"], **v} for s in SITE for v in s["unit_site"]["units"].values()])
print("Units with a site:", len(wch), "; waveform channels equal to the shank's electrodes for all units:", bool((wch.waveform_channels == wch.n_channels).all()))
print("\nElectrodes table of", SITE[0]["session"], "(first 10 rows; x, y, z are NaN):")
print(SITE[0]["electrodes"].head(10).to_string())
print("Columns of the LFP and of the wide-band signal: electrode rows in table order in every site file:",
      all(np.array_equal(s["unit_site"]["region"], np.arange(len(s["electrodes"]))) for s in SITE))
# DS2 amplitude against the site's number along the shank (shank_electrode_number); if the numbers
# followed depth, each profile would change smoothly and cross zero once.
prof = []
for s in SITE:
    us, el = s["unit_site"], s["electrodes"]
    for sh, g in el.groupby("group_name"):
        a = 1e6 * us["channel_amp"][g.index.values]
        order = np.argsort(g.shank_electrode_number.values)
        a = a[order]
        prof.append({"session": s["session"], "shank": sh, "rho": stats.spearmanr(np.arange(len(a)), a)[0],
                     "sign_changes": int(np.sum(np.diff(np.sign(a)) != 0)), "amp": a})
prof = pd.DataFrame(prof)
print("\nShanks: %d; Spearman rho of DS2 amplitude against site number: median %.2f, |rho| >= 0.8 on %d; "
      "shanks whose DS2 changes sign more than once along the numbering: %d of %d with any change" % (
      len(prof), prof.rho.median(), int((prof.rho.abs() >= 0.8).sum()), int((prof.sign_changes > 1).sum()), int((prof.sign_changes > 0).sum())))
fig, axs = plt.subplots(1, 8, figsize=(14, 2.5), sharey=True)
for ax, r in zip(axs, prof[prof.session == SITE[4]["session"]].itertuples()):
    ax.plot(r.amp, np.arange(len(r.amp)), "o-"); ax.axvline(0, c="k", lw=0.5); ax.set_title(r.shank, fontsize=8)
axs[0].set_ylabel("site number"); axs[0].invert_yaxis(); fig.suptitle(f"DS2 amplitude (uV) along each shank, {SITE[4]['session']}", fontsize=9)
print("Plotted file: DS2 amplitude range %.0f to %.0f uV" % (np.min(np.concatenate(prof[prof.session == SITE[4]['session']].amp.values)),
      np.max(np.concatenate(prof[prof.session == SITE[4]['session']].amp.values))))

# The same units in the wide-band signal: the column where each unit's spike-triggered mean is largest,
# against the column the electrodes table gives for the largest channel of its stored waveform.
mp = []
for s in SITE:
    ru, us, el = s["raw_units"], s["unit_site"]["units"], s["electrodes"]
    for u, v in ru["units"].items():
        if u not in us or v["n"] < 20:
            continue
        w = s["waveforms"][u]
        mp.append({"session": s["session"], "unit": u, "n": v["n"], "raw_col": v["site_col"], "table_col": us[u]["col"],
                   "same_col": v["site_col"] == us[u]["col"], "raw_on_table_shank": el.group_name[v["site_col"]] == us[u]["shank"],
                   "raw_uV": 1e6 * v["profile"][v["site_col"]], "stored_uV": 1e6 * np.ptp(w, axis=1).max(),
                   "second_to_first": np.sort(v["profile"])[-2] / v["profile"].max()})
MP = pd.DataFrame(mp)
print(f"\nWide-band signal: {SITE[0]['raw_units']['seconds']:.0f} s per file; units with at least 20 spikes there: {len(MP)}")
lm = [s["raw_units"]["lfp_match"] for s in SITE]
print("Wide-band columns whose decimated signal correlates best with the LFP column of the same index: %d of %d (median r %.3f)" % (
      sum(m[0] for m in lm), sum(s["raw_units"]["columns"] for s in SITE), np.median([m[1] for m in lm])))
print("Largest wide-band column equal to the column the table gives: %.2f; on the shank the table gives: %.2f" % (MP.same_col.mean(), MP.raw_on_table_shank.mean()))
print("By file (fraction on the table's shank):", MP.groupby("session").raw_on_table_shank.mean().round(2).to_dict())
print("Peak-to-peak amplitude, median: stored waveforms %.0f uV, wide-band (300 Hz high-pass) %.0f uV; ratio %.1f" % (
      MP.stored_uV.median(), MP.raw_uV.median(), (MP.stored_uV / MP.raw_uV).median()))
print("Wide-band: second-largest column over the largest, median %.2f" % MP.second_to_first.median())
th = [{"session": s["session"], "theta_ref_col": int(np.flatnonzero(s["electrodes"].theta_reference.values)[0]),
       "table_shank": s["electrodes"].group_name[int(np.flatnonzero(s["electrodes"].theta_reference.values)[0])],
       "ds2_uV_there": 1e6 * s["unit_site"]["channel_amp"][int(np.flatnonzero(s["electrodes"].theta_reference.values)[0])]}
      for s in SITE if s["electrodes"].theta_reference.sum() == 1]
print("\nTheta-reference channel in the site files (column, table shank, DS2 amplitude there in uV):")
print(pd.DataFrame(th).round(0).to_string(index=False))

Files with unit sites: 17 ; mice: ['YutaMouse33', 'YutaMouse41', 'YutaMouse42', 'YutaMouse54']
DS2 events used per file: median 121, range 8 to 365; LFP window 220 s
Shanks: 144; with a DS2 polarity reversal along the shank: 64
DS2 amplitude range across channels (uV): -632 to 786
Units with a site: 860 ; waveform channels equal to the shank's electrodes for all units: True

Electrodes table of YutaMouse33-150222 (first 10 rows; x, y, z are NaN):
  group_name  shank_electrode_number  theta_reference filtering location
0     shank1                       0             True      none  unknown
1     shank1                       1            False      none  unknown
2     shank1                       2            False      none  unknown
3     shank1                       3            False      none  unknown
4     shank1                       4            False      none  unknown
5     shank1                       5            False      none  unknown
6     shank1                       6            False      none  unknown
7     shank1                       7            False      none  unknown
8     shank2                       0            False      none  unknown
9     shank2                       1            False      none  unknown
Columns of the LFP and of the wide-band signal: electrode rows in table order in every site file: True

Shanks: 144; Spearman rho of DS2 amplitude against site number: median 0.13, |rho| >= 0.8 on 15; shanks whose DS2 changes sign more than once along the numbering: 43 of 64 with any change

Plotted file: DS2 amplitude range -632 to 617 uV

Wide-band signal: 149 s per file; units with at least 20 spikes there: 786
Wide-band columns whose decimated signal correlates best with the LFP column of the same index: 1216 of 1216 (median r 0.999)
Largest wide-band column equal to the column the table gives: 0.09; on the shank the table gives: 0.45
By file (fraction on the table's shank): {'YutaMouse33-150222': 0.28, 'YutaMouse33-150223': 0.24, 'YutaMouse33-150224': 0.39, 'YutaMouse33-150225': 0.49, 'YutaMouse41-150902': 0.38, 'YutaMouse41-150903': 0.39, 'YutaMouse41-150908': 0.41, 'YutaMouse41-150909': 0.37, 'YutaMouse41-150910': 0.27, 'YutaMouse41-150911': 0.11, 'YutaMouse42-151105': 0.76, 'YutaMouse42-151108': 0.77, 'YutaMouse42-151110': 0.69, 'YutaMouse42-151111': 0.63, 'YutaMouse42-151112': 0.79, 'YutaMouse54-160701': 0.46, 'YutaMouse54-160702': 0.29}
Peak-to-peak amplitude, median: stored waveforms 421 uV, wide-band (300 Hz high-pass) 81 uV; ratio 5.3
Wide-band: second-largest column over the largest, median 0.57

Theta-reference channel in the site files (column, table shank, DS2 amplitude there in uV):
           session  theta_ref_col table_shank  ds2_uV_there
YutaMouse33-150222              0      shank1         671.0
YutaMouse33-150223              2      shank1         631.0
YutaMouse33-150224             52      shank7         571.0
YutaMouse33-150225             48      shank7         556.0
YutaMouse41-150902              0      shank1         617.0
YutaMouse41-150903              0      shank1         610.0
YutaMouse41-150908             13      shank2         693.0
YutaMouse41-150909             13      shank2         743.0
YutaMouse41-150910             13      shank2         704.0
YutaMouse41-150911             26      shank4         526.0
YutaMouse42-151105             60      shank8         786.0
YutaMouse42-151108             10      shank2         584.0
YutaMouse42-151110             56      shank8         750.0
YutaMouse42-151111             56      shank8         526.0
YutaMouse42-151112             56      shank8         683.0
YutaMouse54-160701             31      shank4         459.0
YutaMouse54-160702             21      shank3         680.0

# Which side of the mismatch is wrong? Each unit's shank (`electrode_group`) and the channel of its
# stored waveform where the spike is largest are paired with the wide-band column where its spike is
# largest. If the shanks and stored channels are right, each (shank, channel) should point to one
# column, the same in every file of the 4 mice (one probe layout). The map is built from the other
# mice and tested on each mouse in turn.
VOT = pd.DataFrame([{"session": s["session"], "subject": s["subject"], "unit": u, "shank": s["unit_site"]["units"][u]["shank"],
                     "channel": s["unit_site"]["units"][u]["channel"], "col": v["site_col"], "table_col": s["unit_site"]["units"][u]["col"]}
                    for s in SITE for u, v in s["raw_units"]["units"].items() if u in s["unit_site"]["units"] and v["n"] >= 20])
def vote_map(d):
    return d.groupby(["shank", "channel"]).col.agg(lambda c: c.value_counts().index[0])
CMAP = vote_map(VOT)
rng_map = np.random.default_rng(20170301)
held = []
for subj in sorted(VOT.subject.unique()):
    train, test = VOT[VOT.subject != subj], VOT[VOT.subject == subj]
    m = vote_map(train)
    perm = train.copy()
    for sess, idx in perm.groupby("session").groups.items():
        perm.loc[idx, ["shank", "channel"]] = perm.loc[idx, ["shank", "channel"]].values[rng_map.permutation(len(idx))]
    mp_ = vote_map(perm)
    key = list(zip(test.shank, test.channel))
    held.append({"subject": subj, "units": len(test),
                 "map_from_other_mice": np.mean([m.get(k, -1) == c for k, c in zip(key, test.col)]),
                 "same_with_units_shuffled": np.mean([mp_.get(k, -1) == c for k, c in zip(key, test.col)]),
                 "electrodes_table": float((test.table_col == test.col).mean())})
held = pd.DataFrame(held)
print(f"Units with at least 20 wide-band spikes: {len(VOT)}; (shank, stored channel) pairs seen: {len(CMAP)}, "
      f"pointing to {CMAP.nunique()} distinct columns")
print(f"Units whose wide-band column is the one their (shank, channel) points to, map from all mice: {(VOT.col.values == CMAP.loc[list(zip(VOT.shank, VOT.channel))].values).mean():.2f}")
print("Map from the other mice, tested on each mouse, against the same map with shanks and channels shuffled among units within each file, and the electrodes table:")
print(held.round(2).to_string(index=False))
el0 = SITE[0]["electrodes"]
print(f"\nRecovered map (stored channel 0 to 7 -> wide-band column; the shank of that column in the table of {SITE[0]['session']}):")
for sh in sorted(CMAP.index.get_level_values(0).unique()):
    cols = [int(CMAP.get((sh, ch), -1)) for ch in range(8)]
    print(f"  {sh}: {cols}  table: {[el0.group_name[c] if c >= 0 else '-' for c in cols]}")

# DS2 amplitude along the stored channel order of the recovered map, against the table's site numbers
# (the profiles of the data check above).
prof2 = []
for s in SITE:
    amp = 1e6 * s["unit_site"]["channel_amp"]
    for sh in sorted(CMAP.index.get_level_values(0).unique()):
        cols = [int(CMAP[(sh, ch)]) for ch in range(8) if (sh, ch) in CMAP.index]
        a = amp[cols]
        prof2.append({"session": s["session"], "shank": sh, "sites": len(a), "rho": stats.spearmanr(np.arange(len(a)), a)[0],
                      "sign_changes": int(np.sum(np.diff(np.sign(a)) != 0))})
prof2 = pd.DataFrame(prof2)
for name, d in (("recovered stored-channel order", prof2), ("table's site numbers", prof)):
    print(f"DS2 along the {name}: {len(d)} shanks, |Spearman rho| median {d.rho.abs().median():.2f}, >= 0.8 on {int((d.rho.abs() >= 0.8).sum())}; "
          f"crossing zero once on {int((d.sign_changes == 1).sum())} and more than once on {int((d.sign_changes > 1).sum())} of {int((d.sign_changes > 0).sum())} shanks that cross")

# The stored waveforms against the wide-band spike at the unit's site: which filtering of the
# wide-band mean they match, and their scale. Units with at least 100 wide-band spikes; the stored
# 32 samples are aligned to the 72-sample wide-band mean at the offset of best correlation.
wf = []
for s in SITE:
    conv = s["raw_units"]["conversion"]
    for u, v in s["raw_units"]["units"].items():
        if u not in s["waveforms"] or v["n"] < 100:
            continue
        w = s["waveforms"][u]
        st = w[int(np.argmax(np.ptp(w, axis=1)))].astype(float)
        versions = {"unfiltered": v["wave"].astype(float), "300 Hz": v["wave_hp"].astype(float)}
        for hz in (500, 800):
            versions[f"{hz} Hz"] = signal.filtfilt(*signal.butter(3, hz / (dg.SPIKE_HZ / 2), btype="high"), v["wave"].astype(float), padtype="even")
        r = {"session": s["session"], "unit": u, "conv_uV": 1e6 * conv}
        for name, x in versions.items():
            fits = []
            for o in range(len(x) - len(st) + 1):
                seg = x[o:o + len(st)] - x[o:o + len(st)].mean()
                fits.append((np.corrcoef(seg, st)[0, 1], o, np.dot(seg, st - st.mean()) / np.dot(seg, seg)))
            c, o, g = max(fits)
            r[f"r {name}"], r[f"offset {name}"], r[f"gain {name}"] = c, o, g
        wf.append(r)
WF = pd.DataFrame(wf)
print(f"\nStored against wide-band waveforms: {len(WF)} units with at least 100 wide-band spikes")
for name in ("unfiltered", "300 Hz", "500 Hz", "800 Hz"):
    print(f"  wide-band {name:>10}: correlation median {WF[f'r {name}'].median():.3f}; offset (samples) median {WF[f'offset {name}'].median():.0f}; "
          f"gain median {WF[f'gain {name}'].median():.2f} (IQR {WF[f'gain {name}'].quantile(0.25):.2f} to {WF[f'gain {name}'].quantile(0.75):.2f})")
print("Best-matching filtering per unit:", WF[[f"r {n}" for n in ("unfiltered", "300 Hz", "500 Hz", "800 Hz")]].idxmax(axis=1).str[2:].value_counts().to_dict())
print(f"Wide-band conversion factor: {WF.conv_uV.unique()} uV per count; its inverse {1 / WF.conv_uV.median():.2f}")

Units with at least 20 wide-band spikes: 786; (shank, stored channel) pairs seen: 63, pointing to 63 distinct columns
Units whose wide-band column is the one their (shank, channel) points to, map from all mice: 0.88
Map from the other mice, tested on each mouse, against the same map with shanks and channels shuffled among units within each file, and the electrodes table:
    subject  units  map_from_other_mice  same_with_units_shuffled  electrodes_table
YutaMouse33    149                 0.89                      0.08              0.08
YutaMouse41    311                 0.83                      0.00              0.05
YutaMouse42    232                 0.78                      0.09              0.13
YutaMouse54     94                 0.78                      0.05              0.09

Recovered map (stored channel 0 to 7 -> wide-band column; the shank of that column in the table of YutaMouse33-150222):
  shank1: [33, -1, 32, 39, 35, 36, 34, 37]  table: ['shank5', '-', 'shank5', 'shank5', 'shank5', 'shank5', 'shank5', 'shank5']
  shank2: [41, 46, 40, 47, 43, 44, 42, 45]  table: ['shank6', 'shank6', 'shank6', 'shank6', 'shank6', 'shank6', 'shank6', 'shank6']
  shank3: [16, 25, 17, 24, 18, 22, 20, 21]  table: ['shank3', 'shank4', 'shank3', 'shank4', 'shank3', 'shank3', 'shank3', 'shank3']
  shank4: [26, 19, 28, 23, 29, 27, 30, 31]  table: ['shank4', 'shank3', 'shank4', 'shank3', 'shank4', 'shank4', 'shank4', 'shank4']
  shank5: [13, 4, 9, 2, 5, 3, 1, 0]  table: ['shank2', 'shank1', 'shank2', 'shank1', 'shank1', 'shank1', 'shank1', 'shank1']
  shank6: [7, 14, 6, 15, 8, 12, 11, 10]  table: ['shank1', 'shank2', 'shank1', 'shank2', 'shank2', 'shank2', 'shank2', 'shank2']
  shank7: [48, 55, 49, 54, 50, 53, 51, 52]  table: ['shank7', 'shank7', 'shank7', 'shank7', 'shank7', 'shank7', 'shank7', 'shank7']
  shank8: [56, 63, 57, 62, 58, 61, 59, 60]  table: ['shank8', 'shank8', 'shank8', 'shank8', 'shank8', 'shank8', 'shank8', 'shank8']
DS2 along the recovered stored-channel order: 136 shanks, |Spearman rho| median 0.88, >= 0.8 on 83; crossing zero once on 43 and more than once on 8 of 51 shanks that cross
DS2 along the table's site numbers: 144 shanks, |Spearman rho| median 0.35, >= 0.8 on 15; crossing zero once on 21 and more than once on 43 of 64 shanks that cross

Stored against wide-band waveforms: 534 units with at least 100 wide-band spikes
  wide-band unfiltered: correlation median 0.981; offset (samples) median 8; gain median 4.82 (IQR 4.43 to 5.15)
  wide-band     300 Hz: correlation median 0.987; offset (samples) median 8; gain median 5.08 (IQR 4.72 to 5.45)
  wide-band     500 Hz: correlation median 0.972; offset (samples) median 8; gain median 5.63 (IQR 5.17 to 6.08)
  wide-band     800 Hz: correlation median 0.867; offset (samples) median 8; gain median 6.27 (IQR 5.67 to 6.97)
Best-matching filtering per unit: {'300 Hz': 289, 'unfiltered': 190, '500 Hz': 53, '800 Hz': 2}
Wide-band conversion factor: [0.195] uV per count; its inverse 5.13

started = time.time()
rows = []
for s in S:
    if not s["spikes"]:
        continue
    r = dg.state_rates(s)
    ut = s["units"].set_index("unit")
    r["type"] = [ut.type[u] for u in r.unit]
    r["bi"] = [dg.burst_index(dg.times(s, u)) for u in r.unit]
    tp = [dg.trough_to_peak_ms(s["waveforms"][u]) if u in s["waveforms"] else (np.nan, False, None) for u in r.unit]
    r["tp"] = [x[0] for x in tp]
    r["tp_at_end"] = [x[1] for x in tp]
    r["trough_ms"] = [np.argmin(x[2]) / 100.0 if x[2] is not None and np.isfinite(x[0]) else np.nan for x in tp]
    ru = s["raw_units"]["units"] if s["raw_units"] else {}
    raw = [(dg.trough_to_peak_ms(ru[u]["wave"][None])[:2] + dg.trough_to_peak_ms(ru[u]["wave_hp"][None])[:2]
            + (ru[u]["n"], ru[u]["site_col"], 1e6 * s["unit_site"]["channel_amp"][ru[u]["site_col"]]))
           if u in ru else (np.nan, False, np.nan, False, 0, -1, np.nan) for u in r.unit]
    for k, name in enumerate(("tp_raw", "tp_raw_at_end", "tp_rawhp", "tp_rawhp_at_end", "raw_n", "raw_col", "ds2_rawsite_uV")):
        r[name] = [x[k] for x in raw]
    for band in ("theta", "gamma"):
        st = [dg.phase_stats(s["phases"]["units"][u][band]) if s["phases"] else dg.phase_stats(np.array([])) for u in r.unit]
        for k in ("n", "mrl", "phase", "p", "kappa"):
            r[f"{band}_{k}"] = [x[k] for x in st]
    if len(s["ds2"]):
        d = [dg.ds2_response(dg.times(s, u), s["ds2"]) for u in r.unit]
        for k in ("peak_ratio", "peak_ms", "base_hz"):
            r[f"ds2_{k}"] = [x[k] for x in d]
    r["session"], r["subject"], r["genotype"] = s["session"], s["subject"], s["genotype"]
    rows.append(r)
R = pd.concat(rows, ignore_index=True)
for band in ("theta", "gamma"):
    R[f"{band}_sig"] = (R[f"{band}_p"] < 0.05) & (R[f"{band}_kappa"] > 0.1)
site = pd.DataFrame([{"session": s["session"], "unit": u, "ds2_site_uV": 1e6 * v["ds2_amp"], "shank": v["shank"], "site_channel": v["channel"], "table_col": v["col"]}
                     for s in S if s["unit_site"] for u, v in s["unit_site"]["units"].items()])
site["dentate_shank"] = [not (U[(U.session == a) & (U.electrode_group == b)].type == "CA3").any()
                         and (U[(U.session == a) & (U.electrode_group == b)].type.isin(["GC", "MC"])).any()
                         for a, b in zip(site.session, site.shank)]
R = R.merge(site, on=["session", "unit"], how="left")
print(f"{len(R)} units measured in {time.time() - started:.0f} s")
print(R.groupby("type")[["rate_overall", "rate_awake", "rate_quiet", "rate_nrem", "rate_rem", "bi", "tp"]].median().round(3).to_string())
print("\nDS2 measures missing (no events, or no spikes in the baseline):", R.ds2_peak_ratio.isna().groupby(R.type).sum().to_dict())
print("\nAwake time in the home cage (outside behaviour epochs and position series): %.0f of %.0f h" % (
    R.drop_duplicates("session").dur_awake_home.sum() / 3600, R.drop_duplicates("session").dur_awake.sum() / 3600))
print("Units with a wide-band waveform (at least 20 spikes):", (R.raw_n >= 20).groupby(R.type).sum().to_dict())

2759 units measured in 48 s
      rate_overall  rate_awake  rate_quiet  rate_nrem  rate_rem      bi    tp
type                                                                         
CA3          0.441       0.441       0.265      0.459     0.287   9.818  0.74
GC           0.234       0.077       0.139      0.391     0.198  11.927  0.71
MC           0.803       0.736       0.612      0.913     0.784   9.582  0.74
SST          4.653       5.544       4.641      1.532     4.465   0.422  0.54
nI           5.963       8.748       4.779      4.584     7.499   0.858  0.33
wI           2.574       3.083       1.755      1.732     2.319   0.193  0.60

DS2 measures missing (no events, or no spikes in the baseline): {'CA3': 27, 'GC': 10, 'MC': 14, 'SST': 0, 'nI': 5, 'wI': 3}

Awake time in the home cage (outside behaviour epochs and position series): 174 of 272 h
Units with a wide-band waveform (at least 20 spikes): {'CA3': 161, 'GC': 139, 'MC': 226, 'SST': 0, 'nI': 139, 'wI': 121}

cnt = R.type.value_counts()
print("GC: %d (paper 586, tolerance 527.4 to 644.6); MC: %d (paper 500, tolerance 450 to 550)" % (cnt["GC"], cnt["MC"]))
print("GC + MC: %d (paper 1086 excitatory dentate units)" % (cnt["GC"] + cnt["MC"]))
print("Subjects with GC or MC:", R[R.type.isin(["GC", "MC"])].subject.nunique())
print("\nBy subject:")
print(pd.crosstab(R.subject, R.type)[["GC", "MC"]].to_string())

GC: 609 (paper 586, tolerance 527.4 to 644.6); MC: 568 (paper 500, tolerance 450 to 550)
GC + MC: 1177 (paper 1086 excitatory dentate units)
Subjects with GC or MC: 16

By subject:
type          GC   MC
subject              
YutaMouse20    3   26
YutaMouse23   27   12
YutaMouse33   35   74
YutaMouse37   22    9
YutaMouse38    5    5
YutaMouse39   40   30
YutaMouse40   23    0
YutaMouse41  130   86
YutaMouse42   65  146
YutaMouse44    8    1
YutaMouse45   24    7
YutaMouse51   32    3
YutaMouse54   45   55
YutaMouse55   37    6
YutaMouse56   87   43
YutaMouse57   26   65

rows = []
for thr in dg.BURST_THRESHOLDS:
    n = R[R.type == "nI"]; e = R[R.type.isin(["GC", "MC"])]; w = R[R.type == "wI"]
    rows.append({"burst threshold": thr,
                 "nI below 0.425 ms": (n.tp < dg.TP_NARROW_MS).mean(),
                 "GC+MC wide and bursty": ((e.tp >= dg.TP_NARROW_MS) & (e.bi > thr)).mean(),
                 "wI wide and not bursty": ((w.tp >= dg.TP_NARROW_MS) & (w.bi <= thr)).mean(),
                 "CA3 wide and bursty": ((R[R.type == "CA3"].tp >= dg.TP_NARROW_MS) & (R[R.type == "CA3"].bi > thr)).mean()})
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("\nUnits without a latency (no waveforms, or a positive spike):", R.tp.isna().groupby(R.type).sum().to_dict())
print("Burst index medians:", R.groupby("type").bi.median().round(2).to_dict())
print("Latency medians (ms):", R.groupby("type").tp.median().round(3).to_dict())
print("nI below 0.425 ms by subject:", R[R.type == "nI"].groupby("subject").tp.apply(lambda x: round((x < 0.425).mean(), 2)).to_dict())
fig, ax = plt.subplots(figsize=(6, 4))
for ty, c in (("nI", "tab:blue"), ("wI", "tab:cyan"), ("GC", "gold"), ("MC", "m")):
    x = R[R.type == ty]
    ax.scatter(x.tp, np.log10(x.bi.clip(lower=1e-2)), s=4, c=c, label=ty)
ax.axvline(0.425, c="k", lw=0.5); ax.axhline(np.log10(1.2), c="k", lw=0.5); ax.axhline(np.log10(1.8), c="k", lw=0.5, ls="--")
ax.set_xlabel("trough-to-peak (ms)"); ax.set_ylabel("log10 burst index (clipped at 0.01)"); ax.legend()
print("Plotted latency range: %.2f to %.2f ms; burst index below 0.01 drawn at 0.01: %d units" % (R.tp.min(), R.tp.max(), int((R.bi < 0.01).sum())))
print("\nnI below 0.525 ms (the paper's CA3 threshold; the labels pool dentate and CA3 interneurons): %.3f" % (R[R.type == "nI"].tp < 0.525).mean())
print("Trough of the stored waveform (ms from the start of the 1.6 ms): median %.2f, 5th to 95th percentile %.2f to %.2f;"
      " longest latency that can be measured (1.55 ms minus the trough): median %.2f ms" % (
      R.trough_ms.median(), *np.nanpercentile(R.trough_ms, [5, 95]), (1.55 - R.trough_ms).median()))
print("Peak at the end of the stored waveform:", R.groupby("type").tp_at_end.mean().round(2).to_dict())
X = R[R.raw_n >= 20]
print("\nVariant, 17 files, latency on the wide-band waveform (-1.2 to +2.4 ms) at the unit's wide-band site:")
for col, name in (("tp_raw", "unfiltered"), ("tp_rawhp", "300 Hz high-pass")):
    print(f"  {name}: nI below 0.425 ms {(X[X.type == 'nI'][col] < 0.425).mean():.3f} (n = {(X.type == 'nI').sum()}); "
          f"GC+MC at or above {(X[X.type.isin(['GC', 'MC'])][col] >= 0.425).mean():.3f}; wI at or above {(X[X.type == 'wI'][col] >= 0.425).mean():.3f}; "
          f"peak at the window's end {X[col + '_at_end'].mean():.2f}")
D = R[R.dentate_shank.notna()].astype({"dentate_shank": bool})
d = D[D.dentate_shank & (D.type == "nI")]
print("Variant, 17 files, nI on shanks with GC or MC and no CA3 unit: below 0.425 ms %.3f (n = %d)" % ((d.tp < 0.425).mean(), len(d)))

 burst threshold  nI below 0.425 ms  GC+MC wide and bursty  wI wide and not bursty  CA3 wide and bursty
             1.2              0.825                  0.982                   0.917                0.983
             1.8              0.825                  0.982                   0.933                0.963

Units without a latency (no waveforms, or a positive spike): {'CA3': 5, 'GC': 1, 'MC': 5, 'SST': 0, 'nI': 14, 'wI': 5}
Burst index medians: {'CA3': 9.82, 'GC': 11.93, 'MC': 9.58, 'SST': 0.42, 'nI': 0.86, 'wI': 0.19}
Latency medians (ms): {'CA3': 0.74, 'GC': 0.71, 'MC': 0.74, 'SST': 0.54, 'nI': 0.33, 'wI': 0.6}
nI below 0.425 ms by subject: {'YutaMouse20': 0.81, 'YutaMouse23': 0.83, 'YutaMouse33': 0.89, 'YutaMouse37': 0.4, 'YutaMouse38': 0.67, 'YutaMouse39': 0.91, 'YutaMouse40': 1.0, 'YutaMouse41': 0.84, 'YutaMouse42': 0.88, 'YutaMouse44': 0.83, 'YutaMouse45': 1.0, 'YutaMouse51': 0.71, 'YutaMouse54': 0.71, 'YutaMouse55': 0.66, 'YutaMouse56': 0.85, 'YutaMouse57': 0.81}
Plotted latency range: 0.13 to 0.75 ms; burst index below 0.01 drawn at 0.01: 16 units

nI below 0.525 ms (the paper's CA3 threshold; the labels pool dentate and CA3 interneurons): 0.927
Trough of the stored waveform (ms from the start of the 1.6 ms): median 0.80, 5th to 95th percentile 0.80 to 0.81; longest latency that can be measured (1.55 ms minus the trough): median 0.75 ms
Peak at the end of the stored waveform: {'CA3': 0.78, 'GC': 0.43, 'MC': 0.67, 'SST': 0.05, 'nI': 0.02, 'wI': 0.16}

Variant, 17 files, latency on the wide-band waveform (-1.2 to +2.4 ms) at the unit's wide-band site:
  unfiltered: nI below 0.425 ms 0.647 (n = 139); GC+MC at or above 0.989; wI at or above 0.950; peak at the window's end 0.07
  300 Hz high-pass: nI below 0.425 ms 0.885 (n = 139); GC+MC at or above 0.995; wI at or above 0.909; peak at the window's end 0.00
Variant, 17 files, nI on shanks with GC or MC and no CA3 unit: below 0.425 ms 0.944 (n = 71)

for col, name in (("rate_overall", "(a) spikes over the recording span"), ("rate_awake", "(b) awake periods")):
    n, w = R[R.type == "nI"][col].dropna(), R[R.type == "wI"][col].dropna()
    print(f"{name}: nI median {n.median():.2f} Hz (n = {len(n)}), wI median {w.median():.2f} Hz (n = {len(w)}), "
          f"Kruskal-Wallis P = {stats.kruskal(n, w).pvalue:.2g}")
print("Paper: 7.97 (tolerance 7.17 to 8.77) and 2.11 Hz (1.90 to 2.32).")
print("\nBy subject, reading (a):")
print(R[R.type.isin(["nI", "wI"])].groupby(["subject", "type"]).rate_overall.median().unstack().round(2).to_string())
D = R[R.dentate_shank.notna()].astype({"dentate_shank": bool})
for flag, name in ((True, "shanks with GC or MC and no CA3 unit"), (False, "other shanks")):
    n, w = D[(D.type == "nI") & (D.dentate_shank == flag)].rate_overall, D[(D.type == "wI") & (D.dentate_shank == flag)].rate_overall
    print(f"Variant, 17 files, {name}: nI median {n.median():.2f} Hz (n = {len(n)}), wI median {w.median():.2f} Hz (n = {len(w)}), "
          f"Kruskal-Wallis P = {stats.kruskal(n, w).pvalue:.2g}")

(a) spikes over the recording span: nI median 5.96 Hz (n = 537), wI median 2.57 Hz (n = 375), Kruskal-Wallis P = 5.3e-17
(b) awake periods: nI median 8.75 Hz (n = 537), wI median 3.08 Hz (n = 375), Kruskal-Wallis P = 1.4e-15
Paper: 7.97 (tolerance 7.17 to 8.77) and 2.11 Hz (1.90 to 2.32).

By subject, reading (a):
type            nI    wI
subject                 
YutaMouse20   5.84  3.54
YutaMouse23  12.36  3.21
YutaMouse33   3.54  2.13
YutaMouse37   0.72  3.35
YutaMouse38   5.01  4.42
YutaMouse39  13.94  4.76
YutaMouse40  23.51  2.80
YutaMouse41  14.56  2.41
YutaMouse42   5.27  1.11
YutaMouse44   7.34  0.87
YutaMouse45   1.27  0.56
YutaMouse51   3.26  1.77
YutaMouse54   8.34  4.04
YutaMouse55   1.38  1.77
YutaMouse56   9.11  1.20
YutaMouse57   7.25  3.99
Variant, 17 files, shanks with GC or MC and no CA3 unit: nI median 8.56 Hz (n = 71), wI median 2.29 Hz (n = 101), Kruskal-Wallis P = 3.7e-08
Variant, 17 files, other shanks: nI median 6.11 Hz (n = 73), wI median 1.93 Hz (n = 33), Kruskal-Wallis P = 0.0005

g, m = R[R.type == "GC"].tp.dropna(), R[R.type == "MC"].tp.dropna()
print(f"Trough-to-peak latency: GC median {g.median():.3f} ms (n = {len(g)}), MC median {m.median():.3f} ms (n = {len(m)}), "
      f"Mann-Whitney P = {stats.mannwhitneyu(g, m).pvalue:.2g}")
print("Units whose peak falls at the end of the stored 1.6 ms (latency a lower bound): GC %.2f, MC %.2f" % (
    R[R.type == "GC"].tp_at_end.mean(), R[R.type == "MC"].tp_at_end.mean()))
gg, mm = R[(R.type == "GC") & ~R.tp_at_end].tp.dropna(), R[(R.type == "MC") & ~R.tp_at_end].tp.dropna()
print(f"Variant, units with the peak inside the window: GC {gg.median():.3f} ms (n = {len(gg)}), MC {mm.median():.3f} ms (n = {len(mm)}), P = {stats.mannwhitneyu(gg, mm).pvalue:.2g}")
print("\nBy subject (median ms):")
print(R[R.type.isin(["GC", "MC"])].groupby(["subject", "type"]).tp.median().unstack().round(3).to_string())
X = R[(R.raw_n >= 20) & R.type.isin(["GC", "MC"])]
print("\nVariant, the 17 files with the wide-band signal (window -1.2 to +2.4 ms):")
for col, name in (("tp", "stored waveforms, same units"), ("tp_raw", "wide-band, unfiltered (as the Methods)"), ("tp_rawhp", "wide-band, 300 Hz high-pass")):
    g, m = X[X.type == "GC"][col].dropna(), X[X.type == "MC"][col].dropna()
    print(f"  {name}: GC median {g.median():.3f} ms (n = {len(g)}), MC {m.median():.3f} ms (n = {len(m)}), Mann-Whitney P = {stats.mannwhitneyu(g, m).pvalue:.2g}; "
          f"peak at the window's end GC {X[X.type == 'GC'][col + '_at_end'].mean():.2f}, MC {X[X.type == 'MC'][col + '_at_end'].mean():.2f}")

Trough-to-peak latency: GC median 0.710 ms (n = 608), MC median 0.740 ms (n = 563), Mann-Whitney P = 2.3e-14
Units whose peak falls at the end of the stored 1.6 ms (latency a lower bound): GC 0.43, MC 0.67
Variant, units with the peak inside the window: GC 0.660 ms (n = 344), MC 0.645 ms (n = 182), P = 0.057

By subject (median ms):
type            GC     MC
subject                  
YutaMouse20  0.680  0.740
YutaMouse23  0.650  0.740
YutaMouse33  0.710  0.730
YutaMouse37  0.730  0.740
YutaMouse38  0.570  0.730
YutaMouse39  0.705  0.730
YutaMouse40  0.720    NaN
YutaMouse41  0.690  0.730
YutaMouse42  0.740  0.740
YutaMouse44  0.710  0.740
YutaMouse45  0.715  0.720
YutaMouse51  0.740  0.740
YutaMouse54  0.630  0.740
YutaMouse55  0.670  0.535
YutaMouse56  0.720  0.730
YutaMouse57  0.730  0.740

Variant, the 17 files with the wide-band signal (window -1.2 to +2.4 ms):
  stored waveforms, same units: GC median 0.710 ms (n = 138), MC 0.740 ms (n = 222), Mann-Whitney P = 5.8e-06; peak at the window's end GC 0.44, MC 0.68
  wide-band, unfiltered (as the Methods): GC median 0.755 ms (n = 138), MC 0.840 ms (n = 226), Mann-Whitney P = 0.017; peak at the window's end GC 0.04, MC 0.09
  wide-band, 300 Hz high-pass: GC median 0.710 ms (n = 139), MC 0.720 ms (n = 226), Mann-Whitney P = 0.4; peak at the window's end GC 0.00, MC 0.00

rows = []
for s in S:
    if "Drd2" not in s["genotype"] or len(s["pulse20"]) == 0:
        continue
    ut = s["units"].set_index("unit")
    for u in s["spikes"]:
        t = dg.times(s, u)
        a = dg.opto_response(t, s["pulse20"])
        b = dg.opto_response(t, s["pulse20"], base=(0.050, 0.200))
        rows.append({"session": s["session"], "subject": s["subject"], "type": ut.type[u],
                     "change_pre": a["change"], "p_pre": a["p_decrease"], "change_cycle": b["change"], "base_hz": a["base_hz"]})
O = pd.DataFrame(rows)
O["supp_pre"] = O.change_pre <= -dg.OPTO_DROP
O["supp_cycle"] = O.change_cycle <= -dg.OPTO_DROP
tab = O.groupby("type").agg(n=("type", "size"), pre=("supp_pre", "sum"), frac_pre=("supp_pre", "mean"),
                            cycle=("supp_cycle", "sum"), frac_cycle=("supp_cycle", "mean"),
                            median_change_pre=("change_pre", "median")).round(3)
print("Sessions:", O.session.nunique(), "; mice:", O.subject.nunique())
print(tab.to_string())
print("\nPaper: MC 29 of 107 (0.271, tolerance 0.244 to 0.298), GC 0 of 202, wI 8 of 68, nI 0 of 69.")
for thr in (0.1, 0.2, 0.3, 0.5):
    print(f"threshold {thr:.0%} decrease, baseline before the pulse: MC {np.mean(O[O.type == 'MC'].change_pre <= -thr):.3f}, "
          f"GC {np.mean(O[O.type == 'GC'].change_pre <= -thr):.3f}")
print("With the 20% rule and a one-sided Poisson P < 0.05 (baseline before the pulse):",
      O[O.supp_pre & (O.p_pre < 0.05)].groupby("type").size().to_dict())
print("\nMC suppressed by subject (baseline before the pulse):")
print(O[O.type == "MC"].groupby("subject").agg(n=("type", "size"), suppressed=("supp_pre", "sum")).to_string())

Sessions: 26 ; mice: 4
        n  pre  frac_pre  cycle  frac_cycle  median_change_pre
type                                                          
CA3   219    6     0.027      9       0.041              0.041
GC    198    6     0.030     13       0.066              0.052
MC    106   68     0.642     71       0.670             -0.273
nI    132   14     0.106     24       0.182              0.017
wI     88   20     0.227     23       0.261             -0.029

Paper: MC 29 of 107 (0.271, tolerance 0.244 to 0.298), GC 0 of 202, wI 8 of 68, nI 0 of 69.
threshold 10% decrease, baseline before the pulse: MC 0.802, GC 0.116
threshold 20% decrease, baseline before the pulse: MC 0.642, GC 0.030
threshold 30% decrease, baseline before the pulse: MC 0.472, GC 0.020
threshold 50% decrease, baseline before the pulse: MC 0.066, GC 0.010
With the 20% rule and a one-sided Poisson P < 0.05 (baseline before the pulse): {'CA3': 6, 'GC': 5, 'MC': 68, 'nI': 14, 'wI': 20}

MC suppressed by subject (baseline before the pulse):
              n  suppressed
subject                    
YutaMouse51   3           0
YutaMouse54  55          44
YutaMouse55   6           3
YutaMouse56  42          21

for col, name in (("rate_awake", "wake (a): awake"), ("rate_wake_all", "wake (b): awake + quiet"), ("rate_awake_home", "variant: home-cage awake"), ("rate_nrem", "NREM"), ("rate_rem", "REM")):
    g, m = R[R.type == "GC"][col].dropna(), R[R.type == "MC"][col].dropna()
    print(f"{name:24s} GC median {g.median():.3f} Hz (n = {len(g)}), MC median {m.median():.3f} Hz (n = {len(m)}), Kruskal-Wallis P = {stats.kruskal(g, m).pvalue:.2g}")
print("\nBy subject (median Hz, awake and NREM):")
print(R[R.type.isin(["GC", "MC"])].groupby(["subject", "type"])[["rate_awake", "rate_nrem"]].median().unstack().round(3).to_string())

wake (a): awake          GC median 0.077 Hz (n = 609), MC median 0.736 Hz (n = 568), Kruskal-Wallis P = 1.8e-114
wake (b): awake + quiet  GC median 0.083 Hz (n = 609), MC median 0.727 Hz (n = 568), Kruskal-Wallis P = 8.7e-113
variant: home-cage awake GC median 0.077 Hz (n = 609), MC median 0.692 Hz (n = 568), Kruskal-Wallis P = 8.1e-102
NREM                     GC median 0.391 Hz (n = 609), MC median 0.913 Hz (n = 568), Kruskal-Wallis P = 2e-29
REM                      GC median 0.198 Hz (n = 609), MC median 0.784 Hz (n = 568), Kruskal-Wallis P = 4.4e-60

By subject (median Hz, awake and NREM):
            rate_awake        rate_nrem       
type                GC     MC        GC     MC
subject                                       
YutaMouse20      0.548  0.903     1.981  1.238
YutaMouse23      0.180  1.078     0.713  0.984
YutaMouse33      0.035  0.845     0.111  0.610
YutaMouse37      0.058  0.611     0.412  0.956
YutaMouse38      0.126  0.482     0.368  0.326
YutaMouse39      0.116  0.636     0.416  0.795
YutaMouse40      0.052    NaN     0.335    NaN
YutaMouse41      0.152  0.799     0.571  0.927
YutaMouse42      0.048  0.837     0.259  1.091
YutaMouse44      0.065  0.037     0.278  0.102
YutaMouse45      0.193  0.384     0.559  0.231
YutaMouse51      0.040  0.233     0.510  0.106
YutaMouse54      0.052  0.416     0.232  0.596
YutaMouse55      0.087  0.491     0.339  0.328
YutaMouse56      0.048  0.521     0.432  0.792
YutaMouse57      0.059  0.843     0.233  1.352

def three_states(ty, wake_col):
    x = R[R.type == ty]
    groups = {"wake": x[wake_col].dropna(), "nrem": x.rate_nrem.dropna(), "rem": x.rate_rem.dropna()}
    kw = stats.kruskal(*groups.values()).pvalue
    post = {f"{a}-{b}": min(1.0, 3 * stats.mannwhitneyu(groups[a], groups[b]).pvalue) for a, b in (("nrem", "wake"), ("nrem", "rem"), ("wake", "rem"))}
    med = {k: round(v.median(), 3) for k, v in groups.items()}
    return med, kw, post
for wake_col, name in (("rate_awake", "(a) awake"), ("rate_wake_all", "(b) awake + quiet"), ("rate_awake_home", "variant, home-cage awake")):
    med, kw, post = three_states("GC", wake_col)
    print(f"GC, wake {name}: medians {med}; Kruskal-Wallis P = {kw:.2g}; Bonferroni-adjusted Mann-Whitney P: " + ", ".join(f"{k} {v:.2g}" for k, v in post.items()))
x = R[R.type == "GC"]
print("\nPaired view: GC with NREM above awake: %.2f; with NREM above REM: %.2f (n = %d)" % (
    (x.rate_nrem > x.rate_awake).mean(), (x.rate_nrem > x.rate_rem).mean(), len(x)))
print("By subject (GC medians, Hz):")
print(x.groupby("subject")[["rate_awake", "rate_quiet", "rate_nrem", "rate_rem"]].median().round(3).to_string())

GC, wake (a) awake: medians {'wake': np.float64(0.077), 'nrem': np.float64(0.391), 'rem': np.float64(0.198)}; Kruskal-Wallis P = 7.5e-93; Bonferroni-adjusted Mann-Whitney P: nrem-wake 1.3e-86, nrem-rem 8.7e-31, wake-rem 1.4e-25
GC, wake (b) awake + quiet: medians {'wake': np.float64(0.083), 'nrem': np.float64(0.391), 'rem': np.float64(0.198)}; Kruskal-Wallis P = 3.4e-89; Bonferroni-adjusted Mann-Whitney P: nrem-wake 1.1e-83, nrem-rem 8.7e-31, wake-rem 6.3e-23
GC, wake variant, home-cage awake: medians {'wake': np.float64(0.077), 'nrem': np.float64(0.391), 'rem': np.float64(0.198)}; Kruskal-Wallis P = 3e-89; Bonferroni-adjusted Mann-Whitney P: nrem-wake 1.2e-82, nrem-rem 8.7e-31, wake-rem 3.3e-24

Paired view: GC with NREM above awake: 0.97; with NREM above REM: 0.86 (n = 609)
By subject (GC medians, Hz):
             rate_awake  rate_quiet  rate_nrem  rate_rem
subject                                                 
YutaMouse20       0.548       0.774      1.981     1.154
YutaMouse23       0.180       0.153      0.713     0.403
YutaMouse33       0.035       0.036      0.111     0.062
YutaMouse37       0.058       0.088      0.412     0.280
YutaMouse38       0.126       0.264      0.368     0.143
YutaMouse39       0.116       0.157      0.416     0.316
YutaMouse40       0.052       0.063      0.335     0.152
YutaMouse41       0.152       0.306      0.571     0.272
YutaMouse42       0.048       0.112      0.259     0.133
YutaMouse44       0.065       0.204      0.278     0.157
YutaMouse45       0.193       0.243      0.559     0.304
YutaMouse51       0.040       0.161      0.510     0.111
YutaMouse54       0.052       0.106      0.232     0.194
YutaMouse55       0.087       0.121      0.339     0.104
YutaMouse56       0.048       0.082      0.432     0.267
YutaMouse57       0.059       0.063      0.233     0.126

for wake_col, name in (("rate_awake", "(a) awake"), ("rate_wake_all", "(b) awake + quiet"), ("rate_awake_home", "variant, home-cage awake")):
    med, kw, post = three_states("CA3", wake_col)
    print(f"CA3, wake {name}: medians {med}; Kruskal-Wallis P = {kw:.2g}; adjusted P: " + ", ".join(f"{k} {v:.2g}" for k, v in post.items()))
print("\nBy subject (CA3 medians, Hz):")
print(R[R.type == "CA3"].groupby("subject")[["rate_awake", "rate_quiet", "rate_nrem", "rate_rem"]].median().round(3).to_string())
x = R[R.type == "CA3"]
print("CA3 awake rate in the home cage and during behaviour (units with both): medians %.3f and %.3f Hz (n = %d)" % (
    x.dropna(subset=["rate_awake_home", "rate_awake_maze"]).rate_awake_home.median(), x.dropna(subset=["rate_awake_home", "rate_awake_maze"]).rate_awake_maze.median(),
    len(x.dropna(subset=["rate_awake_home", "rate_awake_maze"]))))

CA3, wake (a) awake: medians {'wake': np.float64(0.441), 'nrem': np.float64(0.459), 'rem': np.float64(0.287)}; Kruskal-Wallis P = 9.9e-16; adjusted P: nrem-wake 0.77, nrem-rem 1.1e-15, wake-rem 5.4e-09
CA3, wake (b) awake + quiet: medians {'wake': np.float64(0.433), 'nrem': np.float64(0.459), 'rem': np.float64(0.287)}; Kruskal-Wallis P = 1.2e-15; adjusted P: nrem-wake 0.59, nrem-rem 1.1e-15, wake-rem 8.5e-09
CA3, wake variant, home-cage awake: medians {'wake': np.float64(0.412), 'nrem': np.float64(0.459), 'rem': np.float64(0.287)}; Kruskal-Wallis P = 3.2e-14; adjusted P: nrem-wake 0.14, nrem-rem 1.1e-15, wake-rem 4.1e-06

By subject (CA3 medians, Hz):
             rate_awake  rate_quiet  rate_nrem  rate_rem
subject                                                 
YutaMouse39       0.425       0.331      0.629     0.350
YutaMouse41       0.621       0.411      0.568     0.407
YutaMouse42       0.541       0.362      0.462     0.314
YutaMouse54       0.109       0.093      0.187     0.123
YutaMouse55       0.806       0.405      0.641     0.308
YutaMouse56       0.100       0.144      0.356     0.305
YutaMouse57       0.391       0.139      0.357     0.180
CA3 awake rate in the home cage and during behaviour (units with both): medians 0.424 and 0.415 Hz (n = 476)

for wake_col, name in (("rate_awake", "(a) awake"), ("rate_wake_all", "(b) awake + quiet"), ("rate_awake_home", "variant, home-cage awake")):
    ok = R[wake_col].fillna(0) > 0
    ratio = (R.rate_nrem / R[wake_col])[ok]
    g, m = ratio[R.type[ok] == "GC"].dropna(), ratio[R.type[ok] == "MC"].dropna()
    mc = R[ok & (R.type == "MC")]
    print(f"wake {name}: NREM / wake ratio GC median {g.median():.2f} (n = {len(g)}), MC median {m.median():.2f} (n = {len(m)}); "
          f"Mann-Whitney P = {stats.mannwhitneyu(g, m).pvalue:.2g}; MC NREM against wake, Wilcoxon P = {stats.wilcoxon(mc.rate_nrem, mc[wake_col]).pvalue:.2g}")
ok = R.rate_awake > 0
print("\nBy subject, reading (a), median ratio:")
print((R.rate_nrem / R.rate_awake)[ok].groupby([R.subject[ok], R.type[ok]]).median().unstack()[["GC", "MC"]].round(2).to_string())

wake (a) awake: NREM / wake ratio GC median 4.63 (n = 609), MC median 1.11 (n = 568); Mann-Whitney P = 3.9e-141; MC NREM against wake, Wilcoxon P = 1.4e-06
wake (b) awake + quiet: NREM / wake ratio GC median 4.28 (n = 609), MC median 1.11 (n = 568); Mann-Whitney P = 4e-140; MC NREM against wake, Wilcoxon P = 1e-08
wake variant, home-cage awake: NREM / wake ratio GC median 4.56 (n = 609), MC median 1.10 (n = 568); Mann-Whitney P = 2e-120; MC NREM against wake, Wilcoxon P = 3.2e-05

By subject, reading (a), median ratio:
type           GC    MC
subject                
YutaMouse20  3.62  1.27
YutaMouse23  3.43  0.84
YutaMouse33  3.71  0.75
YutaMouse37  5.41  1.11
YutaMouse38  2.76  0.41
YutaMouse39  4.06  1.15
YutaMouse40  4.94   NaN
YutaMouse41  3.70  1.15
YutaMouse42  5.66  1.14
YutaMouse44  5.53  2.74
YutaMouse45  3.43  0.74
YutaMouse51  8.56  1.04
YutaMouse54  4.28  1.06
YutaMouse55  3.17  0.68
YutaMouse56  8.72  1.42
YutaMouse57  4.42  1.31

def circ_mean_deg(v):
    return float(np.degrees(np.angle(np.exp(1j * v).mean())))

def phase_claim(band):
    x = R[R.type.isin(["GC", "MC"])]
    print(f"Units with a {band} estimate: " + ", ".join(f"{t} {int((x[x.type == t][f'{band}_n'] >= 10).sum())}" for t in ("GC", "MC")))
    print(f"Modulated (Rayleigh P < 0.05, kappa > 0.1): " + ", ".join(f"{t} {int(x[x.type == t][f'{band}_sig'].sum())}" for t in ("GC", "MC")))
    sig = x[x[f"{band}_sig"]]
    g, m = sig[sig.type == "GC"], sig[sig.type == "MC"]
    print(f"MRL of modulated units: GC median {g[f'{band}_mrl'].median():.3f}, MC median {m[f'{band}_mrl'].median():.3f}, "
          f"Kruskal-Wallis P = {stats.kruskal(g[f'{band}_mrl'], m[f'{band}_mrl']).pvalue:.2g}")
    allg, allm = x[(x.type == 'GC') & x[f'{band}_mrl'].notna()], x[(x.type == 'MC') & x[f'{band}_mrl'].notna()]
    print(f"Variant, all units with an estimate: GC median MRL {allg[f'{band}_mrl'].median():.3f}, MC {allm[f'{band}_mrl'].median():.3f}, "
          f"P = {stats.kruskal(allg[f'{band}_mrl'], allm[f'{band}_mrl']).pvalue:.2g}")
    print(f"Spikes per modulated unit: GC median {g[f'{band}_n'].median():.0f}, MC median {m[f'{band}_n'].median():.0f} "
          "(the MRL of a uniform phase with n spikes is about 0.89 / sqrt(n))")
    print(f"Variant, MRL above the uniform-phase expectation 0.89/sqrt(n): GC median {(g[f'{band}_mrl'] - 0.89 / np.sqrt(g[f'{band}_n'])).median():.3f}, "
          f"MC {(m[f'{band}_mrl'] - 0.89 / np.sqrt(m[f'{band}_n'])).median():.3f}, P = "
          f"{stats.kruskal(g[f'{band}_mrl'] - 0.89 / np.sqrt(g[f'{band}_n']), m[f'{band}_mrl'] - 0.89 / np.sqrt(m[f'{band}_n'])).pvalue:.2g}")
    gp, mp = g[f"{band}_phase"], m[f"{band}_phase"]
    print(f"Circular mean of preferred phases (0 = peak, degrees): GC {circ_mean_deg(gp):.1f}, MC {circ_mean_deg(mp):.1f}; "
          f"GC minus MC {np.degrees(np.angle(np.exp(1j * (np.angle(np.exp(1j * gp).mean()) - np.angle(np.exp(1j * mp).mean()))))):.1f}")
    print(f"Kruskal-Wallis on phases: (a) on (-pi, pi] P = {stats.kruskal(gp, mp).pvalue:.2g}; "
          f"(b) on [0, 2 pi) P = {stats.kruskal(np.mod(gp, 2 * np.pi), np.mod(mp, 2 * np.pi)).pvalue:.2g}; "
          f"medians (a) GC {np.degrees(gp.median()):.1f}, MC {np.degrees(mp.median()):.1f}; "
          f"(b) GC {np.degrees(np.mod(gp, 2 * np.pi).median()):.1f}, MC {np.degrees(np.mod(mp, 2 * np.pi).median()):.1f}")
    print("By subject, median MRL of modulated units:")
    print(sig.groupby(["subject", "type"])[f"{band}_mrl"].median().unstack().round(3).to_string())
phase_claim("theta")

Units with a theta estimate: GC 557, MC 563
Modulated (Rayleigh P < 0.05, kappa > 0.1): GC 483, MC 532
MRL of modulated units: GC median 0.422, MC median 0.339, Kruskal-Wallis P = 1.5e-20
Variant, all units with an estimate: GC median MRL 0.387, MC 0.334, P = 6.1e-12
Spikes per modulated unit: GC median 119, MC median 747 (the MRL of a uniform phase with n spikes is about 0.89 / sqrt(n))
Variant, MRL above the uniform-phase expectation 0.89/sqrt(n): GC median 0.324, MC 0.295, P = 0.00056
Circular mean of preferred phases (0 = peak, degrees): GC 30.8, MC 12.0; GC minus MC 18.8
Kruskal-Wallis on phases: (a) on (-pi, pi] P = 1e-17; (b) on [0, 2 pi) P = 0.82; medians (a) GC 29.4, MC 10.1; (b) GC 46.0, MC 42.0
By subject, median MRL of modulated units:
type            GC     MC
subject                  
YutaMouse20  0.629  0.436
YutaMouse23  0.307  0.213
YutaMouse33  0.606  0.373
YutaMouse37  0.426  0.403
YutaMouse38  0.464  0.429
YutaMouse39  0.484  0.276
YutaMouse40  0.325    NaN
YutaMouse41  0.402  0.342
YutaMouse42  0.464  0.358
YutaMouse44  0.296    NaN
YutaMouse45  0.301  0.389
YutaMouse51  0.536  0.480
YutaMouse54  0.392  0.289
YutaMouse55  0.428  0.345
YutaMouse56  0.391  0.301
YutaMouse57  0.372  0.245

phase_claim("gamma")

Units with a gamma estimate: GC 540, MC 559
Modulated (Rayleigh P < 0.05, kappa > 0.1): GC 345, MC 475
MRL of modulated units: GC median 0.319, MC median 0.212, Kruskal-Wallis P = 4.8e-43
Variant, all units with an estimate: GC median MRL 0.271, MC 0.203, P = 1.7e-21
Spikes per modulated unit: GC median 103, MC median 591 (the MRL of a uniform phase with n spikes is about 0.89 / sqrt(n))
Variant, MRL above the uniform-phase expectation 0.89/sqrt(n): GC median 0.219, MC 0.166, P = 4.1e-18
Circular mean of preferred phases (0 = peak, degrees): GC -20.2, MC 1.7; GC minus MC -21.9
Kruskal-Wallis on phases: (a) on (-pi, pi] P = 1.2e-23; (b) on [0, 2 pi) P = 0.18; medians (a) GC -22.3, MC -4.2; (b) GC 322.8, MC 323.6
By subject, median MRL of modulated units:

type            GC     MC
subject                  
YutaMouse20  0.297  0.229
YutaMouse23  0.359  0.253
YutaMouse33  0.424  0.263
YutaMouse37  0.329  0.167
YutaMouse38  0.266  0.155
YutaMouse39  0.262  0.284
YutaMouse40  0.267    NaN
YutaMouse41  0.274  0.200
YutaMouse42  0.352  0.207
YutaMouse44  0.389    NaN
YutaMouse45  0.265  0.348
YutaMouse51  0.376  0.442
YutaMouse54  0.328  0.182
YutaMouse55  0.354  0.255
YutaMouse56  0.295  0.226
YutaMouse57  0.337  0.154

for reading, sel in (("(a) all units with an MRL", R.gamma_mrl.notna()), ("(b) modulated units", R.gamma_sig)):
    for ty in ("GC", "MC"):
        x = R[sel & (R.type == ty)]
        r, p = stats.pearsonr(x.rate_overall, x.gamma_mrl)
        print(f"{reading}, {ty}: r = {r:.3f}, P = {p:.2g}, n = {len(x)}")
print("Paper: GC -0.176 (tolerance -0.194 to -0.158), MC -0.184 (-0.202 to -0.166).")
for ty in ("GC", "MC"):
    x = R[R.gamma_mrl.notna() & (R.type == ty)]
    print(f"Variant, log10 rate, {ty}: r = {stats.pearsonr(np.log10(x.rate_overall), x.gamma_mrl)[0]:.3f}; "
          f"Spearman rho = {stats.spearmanr(x.rate_overall, x.gamma_mrl)[0]:.3f}")

(a) all units with an MRL, GC: r = -0.177, P = 3.6e-05, n = 540
(a) all units with an MRL, MC: r = -0.161, P = 0.00013, n = 559
(b) modulated units, GC: r = -0.397, P = 1.9e-14, n = 345
(b) modulated units, MC: r = -0.253, P = 2.3e-08, n = 475
Paper: GC -0.176 (tolerance -0.194 to -0.158), MC -0.184 (-0.202 to -0.166).

Variant, log10 rate, GC: r = -0.211; Spearman rho = -0.188
Variant, log10 rate, MC: r = -0.283; Spearman rho = -0.177

x = R[R.type.isin(["GC", "MC"]) & R.ds2_peak_ratio.notna() & np.isfinite(R.ds2_peak_ratio)]
g, m = x[x.type == "GC"], x[x.type == "MC"]
print(f"Peak ratio: GC median {g.ds2_peak_ratio.median():.2f} (n = {len(g)}), MC median {m.ds2_peak_ratio.median():.2f} (n = {len(m)}), "
      f"Kruskal-Wallis P = {stats.kruskal(g.ds2_peak_ratio, m.ds2_peak_ratio).pvalue:.2g}")
print(f"Peak time (ms from the DS2 peak): GC median {g.ds2_peak_ms.median():.1f}, MC median {m.ds2_peak_ms.median():.1f}, "
      f"Kruskal-Wallis P = {stats.kruskal(g.ds2_peak_ms, m.ds2_peak_ms).pvalue:.2g}")
print("Dropped (no DS2 events or a zero baseline): GC %d, MC %d" % ((R.type == "GC").sum() - len(g), (R.type == "MC").sum() - len(m)))
print("Units peaking within +-4 ms of the DS2 peak: GC %.2f, MC %.2f" % ((g.ds2_peak_ms.abs() <= 4).mean(), (m.ds2_peak_ms.abs() <= 4).mean()))
print("\nBy subject (median peak ratio, median peak ms):")
print(x.groupby(["subject", "type"])[["ds2_peak_ratio", "ds2_peak_ms"]].median().unstack().round(2).to_string())

Peak ratio: GC median 6.88 (n = 599), MC median 4.82 (n = 554), Kruskal-Wallis P = 1.1e-10
Peak time (ms from the DS2 peak): GC median -3.0, MC median 1.0, Kruskal-Wallis P = 2.2e-09
Dropped (no DS2 events or a zero baseline): GC 10, MC 14
Units peaking within +-4 ms of the DS2 peak: GC 0.54, MC 0.72

By subject (median peak ratio, median peak ms):
            ds2_peak_ratio        ds2_peak_ms      
type                    GC     MC          GC    MC
subject                                            
YutaMouse20           7.30   7.98        11.0   2.0
YutaMouse23           9.17   5.49        -3.0  -1.0
YutaMouse33          28.95  15.21        -1.0   3.0
YutaMouse37           7.82   8.58        -3.0   1.0
YutaMouse38          13.00   4.47         1.0  -1.0
YutaMouse39           8.09   5.98        -3.0  -1.0
YutaMouse40          14.29    NaN        -3.0   NaN
YutaMouse41           4.53   3.11        -3.0  -1.0
YutaMouse42           6.00   4.45        -1.0  -1.0
YutaMouse44           6.86   0.00        -3.0 -19.0
YutaMouse45           3.29   5.56        -4.0  -3.0
YutaMouse51          11.16  42.22         1.0  -1.0
YutaMouse54          13.82  13.73        -7.0  -1.0
YutaMouse55          23.94  41.25         1.0   3.0
YutaMouse56           4.11   3.95        -3.0   1.0
YutaMouse57           9.42   3.37        -3.0   1.0

started = time.time()
rows = []
for s in S:
    ut = s["units"].set_index("unit")
    T = {u: dg.times(s, u) for u in s["spikes"]}
    gcs = [u for u in T if ut.type[u] == "GC"]
    mcs = [u for u in T if ut.type[u] == "MC"]
    for g in gcs:
        for m in mcs:
            a = dg.ccg_test(dg.ccg(T[g], T[m]))
            b = dg.ccg_test(dg.ccg(T[m], T[g]))
            rows.append({"session": s["session"], "subject": s["subject"], "gc": g, "mc": m, "gc_to_mc": a["mono"], "mc_to_gc": b["mono"],
                         "peak_ms": a["mono_peak_ms"], "tp": a["mono_excess"] / len(T[g]), "n_gc": len(T[g]), "n_mc": len(T[m]),
                         "gc_mc_shank_id_same": ut.shank_id[g] == ut.shank_id[m]})
CC = pd.DataFrame(rows)
n_gc = int((R.type == "GC").sum())
pairs = CC[CC.gc_to_mc]
print(f"GC-MC pairs tested: {len(CC)} in {CC.session.nunique()} sessions ({time.time() - started:.0f} s)")
print(f"GC to MC pairs: {len(pairs)}; MC to GC pairs: {int(CC.mc_to_gc.sum())}")
print(f"(a) pairs / all GC = {len(pairs)} / {n_gc} = {len(pairs) / n_gc:.4f}")
print(f"(b) GC with a target / all GC = {pairs.groupby(['session', 'gc']).ngroups} / {n_gc} = {pairs.groupby(['session', 'gc']).ngroups / n_gc:.4f}")
print(f"(c) pairs / tested pairs = {len(pairs)} / {len(CC)} = {len(pairs) / len(CC):.4f}")
print("Paper: 0.055 (tolerance 0.0495 to 0.0605), 32 pairs, no MC to GC pair.")
print("Peak lag of detected GC to MC pairs (ms):", pairs.peak_ms.value_counts().sort_index().to_dict())
print("\nVariant thresholds (Poisson percentile): ", end="")
for pct in (0.9999, 0.99999, 0.9999999):
    k1 = k2 = 0
    for s in S:
        for r in CC[CC.session == s["session"]].itertuples():
            k1 += dg.ccg_test(dg.ccg(dg.times(s, r.gc), dg.times(s, r.mc)), pctl=pct)["mono"]
            k2 += dg.ccg_test(dg.ccg(dg.times(s, r.mc), dg.times(s, r.gc)), pctl=pct)["mono"]
    print(f"{pct}: GC to MC {k1}, MC to GC {k2}; ", end="")
print("\n\nBy subject:")
print(CC.groupby("subject").agg(pairs_tested=("gc", "size"), gc_to_mc=("gc_to_mc", "sum"), mc_to_gc=("mc_to_gc", "sum")).to_string())
eg = U.dropna(subset=["electrode_group"]).set_index(["session", "unit"]).electrode_group if "electrode_group" in U else pd.Series(dtype=str)
CC["in_site_files"] = CC.session.isin([s["session"] for s in SITE])
CC["same_shank"] = [eg.get((a, g)) == eg.get((a, m)) if c else np.nan for a, g, m, c in zip(CC.session, CC.gc, CC.mc, CC.in_site_files)]
mc_gc = CC[CC.mc_to_gc]
x = CC[CC.in_site_files]
print(f"\nIn the 17 files with shanks: {len(x)} pairs tested, GC to MC {int(x.gc_to_mc.sum())} ({int(x[x.gc_to_mc].same_shank.sum())} on one shank), MC to GC {int(x.mc_to_gc.sum())}")
first_bin = pairs.peak_ms < 2.0
print(f"GC to MC pairs whose largest bin is the first of the window (1.5 to 2 ms): {int(first_bin.sum())}; without them: {int((~first_bin).sum())} pairs, "
      f"(a) {(~first_bin).sum() / n_gc:.4f}")
smap = {s["session"]: s for s in S}
print("MC to GC pairs (session, GC, MC, also GC to MC, synchrony at -1 to +1 ms):")
for r in mc_gc.itertuples():
    s = smap[r.session]
    sy = dg.ccg_test(dg.ccg(dg.times(s, r.gc), dg.times(s, r.mc)))["sync"]
    print(f"  {r.session} GC {r.gc} MC {r.mc}: GC to MC {r.gc_to_mc}, synchrony {sy}, GC spikes {r.n_gc}, MC spikes {r.n_mc}")
show = [("MC to GC", r.session, r.mc, r.gc) for r in mc_gc.itertuples()] + [("GC to MC, first bin", r.session, r.gc, r.mc) for r in pairs[first_bin].itertuples()]
fig, axs = plt.subplots(int(np.ceil(len(show) / 6)), 6, figsize=(15, 1.9 * np.ceil(len(show) / 6)), squeeze=False)
peak_ratio = []
for ax, (kind, sess, a, b) in zip(axs.ravel(), show):
    s = smap[sess]
    c = dg.ccg(dg.times(s, a), dg.times(s, b))
    t = dg.ccg_test(c)
    lag = (t["left"] + dg.CCG_BIN / 2) * 1000
    m = np.abs(lag) <= 20
    ax.bar(lag[m], c[m], width=0.5, color="0.4"); ax.plot(lag[m], t["hi"][m], "r", lw=0.7)
    ax.set_title(f"{kind}\n{sess[-13:]} {a}->{b}", fontsize=6); ax.tick_params(labelsize=6)
    w = (lag > 1.5) & (lag < 5)
    peak_ratio.append(c[w].max() / t["base"][w].mean())
for ax in axs.ravel()[len(show):]:
    ax.axis("off")
print(f"Plotted: CCGs at -20 to +20 ms of the {len(mc_gc)} MC to GC pairs and the {int(first_bin.sum())} first-bin GC to MC pairs; "
      f"largest count over baseline at +1.5 to +5 ms: median {np.median(peak_ratio):.1f}, range {np.min(peak_ratio):.1f} to {np.max(peak_ratio):.1f}")

GC-MC pairs tested: 4713 in 66 sessions (20 s)
GC to MC pairs: 65; MC to GC pairs: 8
(a) pairs / all GC = 65 / 609 = 0.1067
(b) GC with a target / all GC = 49 / 609 = 0.0805
(c) pairs / tested pairs = 65 / 4713 = 0.0138
Paper: 0.055 (tolerance 0.0495 to 0.0605), 32 pairs, no MC to GC pair.
Peak lag of detected GC to MC pairs (ms): {1.750000000000043: 22, 2.2500000000000435: 33, 2.750000000000044: 3, 3.2500000000000444: 3, 3.750000000000045: 1, 4.250000000000045: 2, 4.750000000000046: 1}

Variant thresholds (Poisson percentile):

0.9999: GC to MC 123, MC to GC 12;

0.99999: GC to MC 78, MC to GC 11;

0.9999999: GC to MC 56, MC to GC 6; 

By subject:
             pairs_tested  gc_to_mc  mc_to_gc
subject                                      
YutaMouse20            16         2         0
YutaMouse23           138         0         1
YutaMouse33           308        15         6
YutaMouse37            34         1         0
YutaMouse38            25         0         0
YutaMouse39           246         0         0
YutaMouse41          1480        24         0
YutaMouse42           817         8         0
YutaMouse44             3         0         0
YutaMouse45            23         0         0
YutaMouse51            17         1         1
YutaMouse54           515         5         0
YutaMouse55            21         0         0
YutaMouse56           782         7         0
YutaMouse57           288         2         0

In the 17 files with shanks: 2147 pairs tested, GC to MC 39 (17 on one shank), MC to GC 6
GC to MC pairs whose largest bin is the first of the window (1.5 to 2 ms): 22; without them: 43 pairs, (a) 0.0706
MC to GC pairs (session, GC, MC, also GC to MC, synchrony at -1 to +1 ms):
  YutaMouse23-140804 GC 5 MC 4: GC to MC False, synchrony False, GC spikes 22605, MC spikes 41544
  YutaMouse33-150223 GC 45 MC 14: GC to MC False, synchrony False, GC spikes 27478, MC spikes 442812
  YutaMouse33-150223 GC 45 MC 55: GC to MC True, synchrony False, GC spikes 27478, MC spikes 159282
  YutaMouse33-150224 GC 3 MC 16: GC to MC True, synchrony True, GC spikes 47327, MC spikes 144372
  YutaMouse33-150224 GC 3 MC 19: GC to MC True, synchrony True, GC spikes 47327, MC spikes 105739
  YutaMouse33-150224 GC 3 MC 26: GC to MC False, synchrony True, GC spikes 47327, MC spikes 135017
  YutaMouse33-150224 GC 3 MC 29: GC to MC True, synchrony True, GC spikes 47327, MC spikes 149571
  YutaMouse51b160516 GC 17 MC 8: GC to MC True, synchrony False, GC spikes 26562, MC spikes 12566

Plotted: CCGs at -20 to +20 ms of the 8 MC to GC pairs and the 22 first-bin GC to MC pairs; largest count over baseline at +1.5 to +5 ms: median 2.4, range 1.4 to 8.0

rows = []
smap = {s["session"]: s for s in S}
for r in pairs.itertuples():
    s = smap[r.session]
    tg, tm = dg.times(s, r.gc), dg.times(s, r.mc)
    pos = dg.burst_positions(tg)
    rows.append({"single": dg.transmission(tg[pos == 0], tm), "first": dg.transmission(tg[pos == 1], tm),
                 "second": dg.transmission(tg[pos == 2], tm), "third": dg.transmission(tg[pos == 3], tm),
                 "n_first": int((pos == 1).sum()), "n_second": int((pos == 2).sum())})
B = pd.DataFrame(rows)
print(f"Pairs: {len(B)}; median GC spikes per pair that are first / second in a burst: {B.n_first.median():.0f} / {B.n_second.median():.0f}")
print("Median transmission probability:", B[["single", "first", "second", "third"]].median().round(4).to_dict())
ok = B[["first", "second"]].notna().all(axis=1)
w = stats.wilcoxon(B.second[ok], B.first[ok])
print(f"Second against first: {int((B.second[ok] > B.first[ok]).sum())} of {int(ok.sum())} pairs higher, paired Wilcoxon P = {w.pvalue:.2g}")
for isi in (0.004, 0.010):
    vals = []
    for r in pairs.itertuples():
        s = smap[r.session]; tg, tm = dg.times(s, r.gc), dg.times(s, r.mc); pos = dg.burst_positions(tg, isi)
        vals.append((dg.transmission(tg[pos == 1], tm), dg.transmission(tg[pos == 2], tm)))
    v = np.array(vals); ok2 = np.all(np.isfinite(v), axis=1)
    print(f"Variant, bursts at most {isi * 1000:.0f} ms apart: median first {np.median(v[ok2, 0]):.4f}, second {np.median(v[ok2, 1]):.4f}, P = {stats.wilcoxon(v[ok2, 1], v[ok2, 0]).pvalue:.2g}")
vals = []
for r in pairs.itertuples():
    s = smap[r.session]; tg, tm = dg.times(s, r.gc), dg.times(s, r.mc); pos = dg.burst_positions(tg)
    nxt, prv = np.r_[np.diff(tg), np.inf], np.r_[np.inf, np.diff(tg)]
    first = tg[(pos == 1) & (nxt >= 0.0035)]
    second = tg[(pos == 2) & (prv >= 0.0035)]
    vals.append((dg.transmission(first, tm), dg.transmission(second, tm), dg.transmission(tg[pos == 1], tm, hollow=0.6),
                 dg.transmission(tg[pos == 2], tm, hollow=0.6), dg.transmission(tg[pos == 0], tm, hollow=0.6)))
v = np.array(vals)
ok2 = np.all(np.isfinite(v[:, :2]), axis=1)
print(f"Variant, burst intervals of 3.5 to 6 ms only (no overlap of the two +1.5 to +5 ms windows): median first {np.median(v[ok2, 0]):.4f}, "
      f"second {np.median(v[ok2, 1]):.4f}, second higher in {int((v[ok2, 1] > v[ok2, 0]).sum())} of {int(ok2.sum())}, P = {stats.wilcoxon(v[ok2, 1], v[ok2, 0]).pvalue:.2g}")
ok3 = np.all(np.isfinite(v[:, 2:4]), axis=1)
print(f"Variant, baseline hollowed by 60%: median single {np.nanmedian(v[:, 4]):.4f}, first {np.median(v[ok3, 2]):.4f}, second {np.median(v[ok3, 3]):.4f}, "
      f"P = {stats.wilcoxon(v[ok3, 3], v[ok3, 2]).pvalue:.2g}")
nf = ~first_bin.values
okn = ok & nf
print(f"Variant, without the {int(first_bin.sum())} first-bin pairs: second higher in {int((B.second[okn] > B.first[okn]).sum())} of {int(okn.sum())}, "
      f"P = {stats.wilcoxon(B.second[okn], B.first[okn]).pvalue:.2g}")

Pairs: 65; median GC spikes per pair that are first / second in a burst: 1267 / 1267
Median transmission probability: {'single': 0.0147, 'first': 0.0066, 'second': 0.0334, 'third': 0.0427}
Second against first: 58 of 65 pairs higher, paired Wilcoxon P = 1.2e-09

Variant, bursts at most 4 ms apart: median first 0.0086, second 0.0310, P = 4.4e-08

Variant, bursts at most 10 ms apart: median first 0.0050, second 0.0335, P = 5.5e-10

Variant, burst intervals of 3.5 to 6 ms only (no overlap of the two +1.5 to +5 ms windows): median first 0.0062, second 0.0344, second higher in 58 of 65, P = 7.6e-10
Variant, baseline hollowed by 60%: median single 0.0149, first 0.0067, second 0.0338, P = 1.2e-09
Variant, without the 22 first-bin pairs: second higher in 41 of 43, P = 3.8e-11

rows = []
for r in pairs.itertuples():
    s = smap[r.session]
    tg, tm = dg.times(s, r.gc), dg.times(s, r.mc)
    out = {}
    for lab in ("awake", "nrem"):
        st = s["states"][s["states"].label == lab]
        keep = np.zeros(len(tg), bool)
        for a, b in zip(st.start_time, st.stop_time):
            keep[np.searchsorted(tg, a):np.searchsorted(tg, b)] = True
        out[lab] = dg.transmission(tg[keep], tm) if keep.sum() >= 50 else np.nan
        out[f"n_{lab}"] = int(keep.sum())
    rows.append(out)
W = pd.DataFrame(rows)
ok = W[["awake", "nrem"]].notna().all(axis=1)
print(f"Pairs with at least 50 GC spikes in both states: {int(ok.sum())} of {len(W)}")
print(f"Median transmission: awake {W.awake[ok].median():.4f}, NREM {W.nrem[ok].median():.4f}; NREM lower in {int((W.nrem[ok] < W.awake[ok]).sum())} pairs; "
      f"paired Wilcoxon P = {stats.wilcoxon(W.nrem[ok], W.awake[ok]).pvalue:.2g}")
okn = ok & ~first_bin.values
print(f"Variant, without the first-bin pairs: median awake {W.awake[okn].median():.4f}, NREM {W.nrem[okn].median():.4f}, NREM lower in "
      f"{int((W.nrem[okn] < W.awake[okn]).sum())} of {int(okn.sum())}, P = {stats.wilcoxon(W.nrem[okn], W.awake[okn]).pvalue:.2g}")
xs = CC[CC.gc_to_mc].in_site_files.values
oks = ok & xs
if oks.sum() >= 5:
    print(f"Variant, pairs in the 17 files with shanks: NREM lower in {int((W.nrem[oks] < W.awake[oks]).sum())} of {int(oks.sum())}, P = {stats.wilcoxon(W.nrem[oks], W.awake[oks]).pvalue:.2g}")

Pairs with at least 50 GC spikes in both states: 65 of 65
Median transmission: awake 0.0125, NREM 0.0108; NREM lower in 34 pairs; paired Wilcoxon P = 0.39
Variant, without the first-bin pairs: median awake 0.0168, NREM 0.0108, NREM lower in 24 of 43, P = 0.27
Variant, pairs in the 17 files with shanks: NREM lower in 20 of 39, P = 0.86

started = time.time()
rows = []
for s in S:
    ut = s["units"].set_index("unit")
    for ty in ("GC", "MC", "CA3"):
        us = [u for u in s["spikes"] if ut.type[u] == ty]
        T = {u: dg.times(s, u) for u in us}
        for i, a in enumerate(us):
            for b in us[i + 1:]:
                rows.append({"session": s["session"], "type": ty, "a": a, "b": b,
                             "sync": dg.ccg_test(dg.ccg(T[a], T[b]))["sync"], "same_shank_id": ut.shank_id[a] == ut.shank_id[b],
                             "same_group": (eg.get((s["session"], a)) == eg.get((s["session"], b))) if (s["session"], a) in eg.index else np.nan})
Z = pd.DataFrame(rows)
print(f"Pairs tested in {time.time() - started:.0f} s")
for ty in ("GC", "MC", "CA3"):
    z = Z[Z.type == ty]
    cells = pd.concat([z[["session", "a"]].rename(columns={"a": "u"}), z[["session", "b"]].rename(columns={"b": "u"})]).drop_duplicates()
    synced = pd.concat([z[z.sync][["session", "a"]].rename(columns={"a": "u"}), z[z.sync][["session", "b"]].rename(columns={"b": "u"})]).drop_duplicates()
    print(f"{ty}: (a) pairs {int(z.sync.sum())} of {len(z)} = {z.sync.mean():.4f}; (b) cells with a synchronous partner {len(synced)} of {len(cells)} = {len(synced) / max(1, len(cells)):.4f}")
print("Paper: MC 0.23 (0.207 to 0.253), CA3 0.09 (0.081 to 0.099), GC pairs below 0.002 (up to 0.0022).")
def sync_fracs(z):
    cells = pd.concat([z[["session", "a"]].rename(columns={"a": "u"}), z[["session", "b"]].rename(columns={"b": "u"})]).drop_duplicates()
    synced = pd.concat([z[z.sync][["session", "a"]].rename(columns={"a": "u"}), z[z.sync][["session", "b"]].rename(columns={"b": "u"})]).drop_duplicates()
    return z.sync.mean(), len(synced) / max(1, len(cells)), len(z), len(cells)
print("\nVariant, the 17 files that give each unit's shank:")
Zs = Z[Z.same_group.notna()]
for ty in ("GC", "MC", "CA3"):
    z = Zs[Zs.type == ty]
    a1 = sync_fracs(z); a2 = sync_fracs(z[z.same_group == False])
    print(f"  {ty}: all pairs {a1[0]:.4f} of {a1[2]} pairs, cells {a1[1]:.3f} of {a1[3]}; pairs on different shanks only {a2[0]:.4f} of {a2[2]}, cells {a2[1]:.3f}; "
          f"synchronous pairs on one shank {int(z[z.same_group == True].sync.sum())} of {int(z.sync.sum())}")
def shank_gap(r):
    return abs(int(eg[(r.session, r.a)].replace("shank", "")) - int(eg[(r.session, r.b)].replace("shank", "")))
Zd = Zs[Zs.same_group == False].copy()
Zd["gap"] = [shank_gap(r) for r in Zd.itertuples()]
print("Synchronous pairs by distance between shanks (shank numbers apart; 200 um per step), 17 files:")
print(pd.crosstab(Zd.type, Zd.gap, values=Zd.sync, aggfunc="mean").round(3).to_string())
print(pd.crosstab(Zd[Zd.sync].type, Zd[Zd.sync].gap).to_string())
# One neuron sorted twice would share most of its spikes with its copy; synchrony from common input
# shares a few percent. Coincidence: fraction of the sparser unit's spikes with a partner spike within 0.5 ms.
smap = {s["session"]: s for s in S}
def coincident(s, a, b, tol=0.0005):
    x, y = dg.times(s, a), dg.times(s, b)
    if len(x) > len(y):
        x, y = y, x
    j = np.clip(np.searchsorted(y, x), 1, len(y) - 1)
    return float((np.minimum(np.abs(y[j] - x), np.abs(y[j - 1] - x)) < tol).mean())
Zy = Z[Z.sync].copy()
Zy["coinc"] = [coincident(smap[r.session], r.a, r.b) for r in Zy.itertuples()]
print("\nSynchronous pairs: fraction of the sparser unit's spikes within 0.5 ms of the other's:")
print(Zy.groupby("type").coinc.describe(percentiles=[0.5, 0.9]).round(3).to_string())
print("Pairs above 0.2:", Zy[Zy.coinc > 0.2].groupby("type").size().to_dict())
Zc = Z[~Z.index.isin(Zy[Zy.coinc > 0.2].index)]
for ty in ("MC", "CA3"):
    f = sync_fracs(Zc[Zc.type == ty])
    print(f"Variant without pairs above 0.2, {ty}: pairs {f[0]:.4f}, cells {f[1]:.3f}")

Pairs tested in 29 s
GC: (a) pairs 4 of 3155 = 0.0013; (b) cells with a synchronous partner 7 of 605 = 0.0116
MC: (a) pairs 281 of 3667 = 0.0766; (b) cells with a synchronous partner 240 of 554 = 0.4332
CA3: (a) pairs 65 of 5395 = 0.0120; (b) cells with a synchronous partner 66 of 645 = 0.1023
Paper: MC 0.23 (0.207 to 0.253), CA3 0.09 (0.081 to 0.099), GC pairs below 0.002 (up to 0.0022).

Variant, the 17 files that give each unit's shank:
  GC: all pairs 0.0010 of 1005 pairs, cells 0.013 of 152; pairs on different shanks only 0.0013 of 773, cells 0.013; synchronous pairs on one shank 0 of 1
  MC: all pairs 0.0852 of 1877 pairs, cells 0.494 of 243; pairs on different shanks only 0.1385 of 1155, cells 0.494; synchronous pairs on one shank 0 of 160
  CA3: all pairs 0.0281 of 1815 pairs, cells 0.249 of 185; pairs on different shanks only 0.0521 of 978, cells 0.267; synchronous pairs on one shank 0 of 51
Synchronous pairs by distance between shanks (shank numbers apart; 200 um per step), 17 files:
gap       1      2    3    4    5
type                             
CA3   0.042  0.097  NaN  NaN  NaN
GC    0.003  0.000  0.0  0.0  0.0
MC    0.155  0.078  0.0  NaN  NaN
gap     1   2
type         
CA3    33  18
GC      1   0
MC    143  17

Synchronous pairs: fraction of the sparser unit's spikes within 0.5 ms of the other's:
      count   mean    std    min    50%    90%    max
type                                                 
CA3    65.0  0.016  0.019  0.003  0.010  0.028  0.146
GC      4.0  0.015  0.006  0.007  0.016  0.020  0.020
MC    281.0  0.010  0.008  0.002  0.008  0.018  0.058
Pairs above 0.2: {}
Variant without pairs above 0.2, MC: pairs 0.0766, cells 0.433
Variant without pairs above 0.2, CA3: pairs 0.0120, cells 0.102

started = time.time()
TMAPS = {}
rows = []
for s in TM:
    ut = s["units"].set_index("unit")
    units = [u for u in s["spikes"] if ut.type[u] in ("GC", "MC", "CA3")]
    if not units:
        continue
    occ, maps, runs = dg.tmaze_maps(s, units)
    TMAPS[s["session"]] = (occ, maps, runs)
    for u in units:
        fields = [(g, f) for g in range(3) for f in dg.linear_fields(maps[u][g])]
        flat = np.concatenate(maps[u])
        rows.append({"session": s["session"], "subject": s["subject"], "unit": u, "type": ut.type[u], "n_fields": len(fields),
                     "field_peak": max([f[2] for _, f in fields], default=np.nan),
                     "info_bits": dg.spatial_info(flat, np.concatenate(occ)), "fields": fields})
F = pd.DataFrame(rows)
print(f"T-maze maps of {len(F)} units in {len(TMAPS)} sessions, {F.subject.nunique()} mice ({time.time() - started:.0f} s)")
print("Units per type:", F.type.value_counts().to_dict(), "(paper: GC 213, MC 225, CA3 249 in 9 mice)")
for ty in ("GC", "MC", "CA3"):
    stack = [TMAPS[r.session][1][r.unit] for r in F[F.type == ty].itertuples()]
    frac = [np.nanmean(np.array([m[g] for m in stack]) > 1.0, axis=0) for g in range(3)]
    central, side = frac[0], np.concatenate(frac[1:])
    print(f"{ty}: fraction active per bin, central median {np.median(central):.3f} ({len(central)} bins), "
          f"side median {np.median(side):.3f} ({len(side)} bins), Kruskal-Wallis P = {stats.kruskal(central, side).pvalue:.2g}")
print("\nVariant, per unit: fraction of central-arm bins above 1 Hz against fraction of side-arm bins above 1 Hz (Wilcoxon over units):")
for ty in ("GC", "MC", "CA3"):
    cen, sid = [], []
    for r in F[F.type == ty].itertuples():
        m = TMAPS[r.session][1][r.unit]
        cen.append(np.nanmean(m[0] > 1.0)); sid.append(np.nanmean(np.concatenate(m[1:]) > 1.0))
    cen, sid = np.array(cen), np.array(sid)
    d = cen - sid
    print(f"  {ty}: mean central {cen.mean():.3f}, side {sid.mean():.3f}; central higher in {int((d > 0).sum())}, lower in {int((d < 0).sum())} of {len(d)}; "
          f"P = {stats.wilcoxon(cen, sid).pvalue if np.any(d != 0) else np.nan:.2g}")

T-maze maps of 897 units in 40 sessions, 7 mice (3 s)
Units per type: {'CA3': 341, 'MC': 303, 'GC': 253} (paper: GC 213, MC 225, CA3 249 in 9 mice)
GC: fraction active per bin, central median 0.020 (35 bins), side median 0.018 (92 bins), Kruskal-Wallis P = 0.043
MC: fraction active per bin, central median 0.314 (35 bins), side median 0.272 (92 bins), Kruskal-Wallis P = 2.2e-06
CA3: fraction active per bin, central median 0.214 (35 bins), side median 0.189 (92 bins), Kruskal-Wallis P = 0.048

Variant, per unit: fraction of central-arm bins above 1 Hz against fraction of side-arm bins above 1 Hz (Wilcoxon over units):
  GC: mean central 0.024, side 0.018; central higher in 17, lower in 31 of 253; P = 0.83
  MC: mean central 0.312, side 0.247; central higher in 146, lower in 104 of 303; P = 2.8e-05
  CA3: mean central 0.199, side 0.174; central higher in 102, lower in 123 of 341; P = 0.21

for ty in ("GC", "MC", "CA3"):
    x = F[F.type == ty].n_fields
    print(f"{ty}: n = {len(x)}; field counts {x.value_counts().sort_index().to_dict()}; (a) two or more of all {np.mean(x >= 2):.3f}; "
          f"(b) of units with a field {np.sum(x >= 2) / max(1, np.sum(x >= 1)):.3f}")
print("Paper: GC 0.066 (0.059 to 0.073), MC 0.578 (0.520 to 0.636), CA3 0.361 (0.325 to 0.397).")
print("Kolmogorov-Smirnov on field counts: GC against MC P = %.2g; GC against CA3 P = %.2g; MC against CA3 P = %.2g" % (
    stats.ks_2samp(F[F.type == "GC"].n_fields, F[F.type == "MC"].n_fields).pvalue,
    stats.ks_2samp(F[F.type == "GC"].n_fields, F[F.type == "CA3"].n_fields).pvalue,
    stats.ks_2samp(F[F.type == "MC"].n_fields, F[F.type == "CA3"].n_fields).pvalue))
print("\nVariant field thresholds, reading (a):")
for frac, nb in ((0.5, 5), (0.5, 3), (0.3, 5), (0.7, 5)):
    out = {}
    for ty in ("GC", "MC", "CA3"):
        cnt = [sum(len(dg.linear_fields(TMAPS[r.session][1][r.unit][g], min_bins=nb, frac=frac)) for g in range(3)) for r in F[F.type == ty].itertuples()]
        out[ty] = round(float(np.mean(np.array(cnt) >= 2)), 3)
    print(f"  above {frac:.0%} of the local peak, at least {nb} bins: {out}")
print("\nBy subject, reading (a):")
print(F.groupby(["subject", "type"]).n_fields.apply(lambda x: round(float(np.mean(x >= 2)), 2)).unstack().to_string())

GC: n = 253; field counts {0: 208, 1: 34, 2: 7, 3: 2, 4: 1, 7: 1}; (a) two or more of all 0.043; (b) of units with a field 0.244
MC: n = 303; field counts {0: 56, 1: 54, 2: 50, 3: 46, 4: 27, 5: 23, 6: 22, 7: 9, 8: 6, 9: 8, 10: 2}; (a) two or more of all 0.637; (b) of units with a field 0.781
CA3: n = 341; field counts {0: 121, 1: 73, 2: 39, 3: 29, 4: 22, 5: 18, 6: 21, 7: 6, 8: 9, 10: 3}; (a) two or more of all 0.431; (b) of units with a field 0.668
Paper: GC 0.066 (0.059 to 0.073), MC 0.578 (0.520 to 0.636), CA3 0.361 (0.325 to 0.397).
Kolmogorov-Smirnov on field counts: GC against MC P = 3.2e-53; GC against CA3 P = 3.3e-29; MC against CA3 P = 1.9e-06

Variant field thresholds, reading (a):
  above 50% of the local peak, at least 5 bins: {'GC': 0.043, 'MC': 0.637, 'CA3': 0.431}

  above 50% of the local peak, at least 3 bins: {'GC': 0.095, 'MC': 0.713, 'CA3': 0.537}

  above 30% of the local peak, at least 5 bins: {'GC': 0.051, 'MC': 0.673, 'CA3': 0.472}
  above 70% of the local peak, at least 5 bins: {'GC': 0.032, 'MC': 0.597, 'CA3': 0.399}

By subject, reading (a):
type          CA3    GC    MC
subject                      
YutaMouse33   NaN  0.00  0.71
YutaMouse37   NaN   NaN  0.00
YutaMouse39  0.51  0.00  0.75
YutaMouse41  0.52  0.06  0.56
YutaMouse42  0.40  0.00  0.68
YutaMouse45   NaN  0.08  0.14
YutaMouse57  0.32  0.04  0.68

pk = {ty: F[(F.type == ty) & (F.n_fields > 0)].field_peak for ty in ("GC", "MC", "CA3")}
inf = {ty: F[F.type == ty].info_bits.dropna() for ty in ("GC", "MC", "CA3")}
print("Within-field peak (Hz), units with a field: " + ", ".join(f"{t} median {v.median():.2f} (n = {len(v)})" for t, v in pk.items()))
print("  KS GC against MC P = %.2g; MC against CA3 P = %.2g" % (stats.ks_2samp(pk["GC"], pk["MC"]).pvalue, stats.ks_2samp(pk["MC"], pk["CA3"]).pvalue))
print("Spatial information (bits/s), all units: " + ", ".join(f"{t} median {v.median():.3f} (n = {len(v)})" for t, v in inf.items()))
print("  KS GC against MC P = %.2g; GC against CA3 P = %.2g; MC against CA3 P = %.2g" % (
    stats.ks_2samp(inf["GC"], inf["MC"]).pvalue, stats.ks_2samp(inf["GC"], inf["CA3"]).pvalue, stats.ks_2samp(inf["MC"], inf["CA3"]).pvalue))
fig, axs = plt.subplots(1, 2, figsize=(9, 3))
for ty, c in (("GC", "gold"), ("MC", "m"), ("CA3", "r")):
    axs[0].hist(pk[ty], bins=np.linspace(0, 40, 41), histtype="step", density=True, color=c, label=ty)
    axs[1].hist(inf[ty], bins=np.linspace(0, 3, 31), histtype="step", density=True, color=c, label=ty)
axs[0].set_xlabel("within-field peak (Hz)"); axs[1].set_xlabel("information (bits/s)"); axs[0].legend()
print("Outside the plotted ranges: peaks above 40 Hz %d, information above 3 bits/s %d" % (
    sum(int((v > 40).sum()) for v in pk.values()), sum(int((v > 3).sum()) for v in inf.values())))

Within-field peak (Hz), units with a field: GC median 2.63 (n = 45), MC median 5.68 (n = 247), CA3 median 4.36 (n = 220)

  KS GC against MC P = 4.5e-05; MC against CA3 P = 0.024
Spatial information (bits/s), all units: GC median 0.035 (n = 232), MC median 0.603 (n = 301), CA3 median 0.293 (n = 329)
  KS GC against MC P = 1.6e-49; GC against CA3 P = 1.5e-29; MC against CA3 P = 3.2e-05
Outside the plotted ranges: peaks above 40 Hz 6, information above 3 bits/s 58

rows = []
for r in F.itertuples():
    central = [f for g, f in r.fields if g == 0]
    if not central:
        continue
    s = smap[r.session]
    runs = TMAPS[r.session][2]
    split = False
    for a, b, _ in central:
        rates = dg.tmaze_run_rates(s, r.unit, runs, 0, a, b)
        side = np.array([x[2] for x in runs])
        df = pd.DataFrame([{"rate": rates[i, j], "pos": j, "choice": side[i]} for i in range(len(runs)) for j in range(rates.shape[1]) if np.isfinite(rates[i, j])])
        if df.choice.nunique() < 2 or len(df) < 10:
            continue
        p = anova_lm(ols("rate ~ C(pos) + C(choice)", df).fit(), typ=2).loc["C(choice)", "PR(>F)"]
        split |= bool(p < 0.001)
    rows.append({"type": r.type, "splitter": split})
SP = pd.DataFrame(rows)
tab = pd.crosstab(SP.type, SP.splitter)
print("Units with a central-arm field, and splitters among them:")
print(tab.to_string())
print("Fractions:", SP.groupby("type").splitter.mean().round(3).to_dict())
chi = stats.chi2_contingency(tab)
print(f"Chi-square across the three types: chi2 = {chi[0]:.2f}, P = {chi[1]:.2g}")
for a, b in (("GC", "MC"), ("GC", "CA3"), ("MC", "CA3")):
    t2 = tab.loc[[a, b]].reindex(columns=[False, True], fill_value=0)
    print(f"Fisher's exact test {a} against {b}: P = {stats.fisher_exact(t2.values)[1]:.2g}")
print("Expected GC splitters under equal fractions: %.2f" % chi[3][list(tab.index).index("GC")][-1])

Units with a central-arm field, and splitters among them:
splitter  False  True 
type                  
CA3         137      4
GC           17      0
MC          175      4
Fractions: {'CA3': 0.028, 'GC': 0.0, 'MC': 0.022}
Chi-square across the three types: chi2 = 0.56, P = 0.76
Fisher's exact test GC against MC: P = 1
Fisher's exact test GC against CA3: P = 1
Fisher's exact test MC against CA3: P = 0.73
Expected GC splitters under equal fractions: 0.40

started = time.time()
OFS = [s for s in S if all(f"{dg.OPEN_FIELDS[k]}_twhl_norm" in s["positions"] for k in ("A", "B", "A2"))]
OMAPS = {}
rows = []
for s in OFS:
    ut = s["units"].set_index("unit")
    units = [u for u in s["spikes"] if ut.type[u] in ("GC", "MC", "CA3")]
    if not units:
        continue
    M = {k: dg.open_field_maps(s, dg.OPEN_FIELDS[k], units, split=True) for k in ("A", "B", "A2")}
    OMAPS[s["session"]] = M
    for u in units:
        r = {"session": s["session"], "subject": s["subject"], "unit": u, "type": ut.type[u]}
        place = False
        for k in ("A", "B", "A2"):
            full, h1, h2 = M[k][u]
            stab = dg.map_corr(h1, h2)
            nf = len(dg.open_field_fields(full))
            r[f"stab_{k}"], r[f"nf_{k}"] = stab, nf
            place |= bool(stab >= 0.5 and nf > 0)
        r["place"] = place
        r["AB"], r["BA2"], r["AA2"] = dg.map_corr(M["A"][u][0], M["B"][u][0]), dg.map_corr(M["B"][u][0], M["A2"][u][0]), dg.map_corr(M["A"][u][0], M["A2"][u][0])
        r["remap"] = r["AA2"] - np.nanmean([r["AB"], r["BA2"]])
        rows.append(r)
OFD = pd.DataFrame(rows)
print(f"Open-field maps of {len(OFD)} units in {len(OMAPS)} sessions, {OFD.subject.nunique()} mice ({time.time() - started:.0f} s)")
print("Units per type:", OFD.type.value_counts().to_dict(), "(paper: GC 247, MC 181, CA3 317 in 5 mice)")
print("Place cells:", OFD[OFD.place].type.value_counts().to_dict(), "(paper: GC 30, MC 111, CA3 144)")
print("Units without a remapping index (no spikes at speed in a session):", OFD.remap.isna().groupby(OFD.type).sum().to_dict())
P = OFD[OFD.place & OFD.remap.notna()]
for ty in ("GC", "MC", "CA3"):
    x = P[P.type == ty].remap
    print(f"{ty}: place cells {len(x)}, median remapping index {x.median():.3f}, Wilcoxon P = {stats.wilcoxon(x).pvalue:.2g}")
print("\nAll units (variant): " + ", ".join(f"{t} median {OFD[OFD.type == t].remap.median():.3f}" for t in ("GC", "MC", "CA3")))
print("By subject, place cells (n, median index):")
print(P.groupby(["subject", "type"]).remap.agg(["size", "median"]).round(3).unstack().to_string())

Open-field maps of 745 units in 35 sessions, 5 mice (4 s)
Units per type: {'CA3': 317, 'GC': 246, 'MC': 182} (paper: GC 247, MC 181, CA3 317 in 5 mice)
Place cells: {'CA3': 145, 'MC': 116, 'GC': 27} (paper: GC 30, MC 111, CA3 144)
Units without a remapping index (no spikes at speed in a session): {'CA3': 17, 'GC': 30, 'MC': 1}
GC: place cells 27, median remapping index 0.128, Wilcoxon P = 0.0013
MC: place cells 116, median remapping index 0.293, Wilcoxon P = 1.5e-13
CA3: place cells 143, median remapping index 0.207, Wilcoxon P = 8.6e-14

All units (variant): GC median 0.010, MC median 0.209, CA3 median 0.121
By subject, place cells (n, median index):
             size             median              
type          CA3    GC    MC    CA3     GC     MC
subject                                           
YutaMouse41  50.0  11.0  45.0  0.150 -0.011  0.308
YutaMouse42  21.0   1.0  36.0  0.079  0.383  0.208
YutaMouse45   NaN   5.0   3.0    NaN  0.253  0.076
YutaMouse55  54.0   3.0   5.0  0.449  0.125  0.057
YutaMouse56  18.0   7.0  27.0  0.124  0.310  0.502

g, m, c = (P[P.type == t].remap for t in ("GC", "MC", "CA3"))
print(f"Medians: GC {g.median():.3f}, MC {m.median():.3f}, CA3 {c.median():.3f}")
print(f"Kruskal-Wallis across the three types: P = {stats.kruskal(g, m, c).pvalue:.3g}")
print(f"Mann-Whitney GC against MC: P = {stats.mannwhitneyu(g, m).pvalue:.3g} (Bonferroni x3: {min(1, 3 * stats.mannwhitneyu(g, m).pvalue):.3g}); "
      f"GC against CA3: P = {stats.mannwhitneyu(g, c).pvalue:.3g}")
print("\nVariant place-cell rule, stability above 0.5 in A or A' only:")
alt = OFD[((OFD.stab_A >= 0.5) & (OFD.nf_A > 0)) | ((OFD.stab_A2 >= 0.5) & (OFD.nf_A2 > 0))].dropna(subset=["remap"])
ga, ma, ca = (alt[alt.type == t].remap for t in ("GC", "MC", "CA3"))
print(f"  n {len(ga)}, {len(ma)}, {len(ca)}; medians {ga.median():.3f}, {ma.median():.3f}, {ca.median():.3f}; KW P = {stats.kruskal(ga, ma, ca).pvalue:.3g}; GC against MC P = {stats.mannwhitneyu(ga, ma).pvalue:.3g}")

Medians: GC 0.128, MC 0.293, CA3 0.207
Kruskal-Wallis across the three types: P = 0.306
Mann-Whitney GC against MC: P = 0.114 (Bonferroni x3: 0.342); GC against CA3: P = 0.259

Variant place-cell rule, stability above 0.5 in A or A' only:
  n 24, 96, 119; medians 0.126, 0.355, 0.271; KW P = 0.194; GC against MC P = 0.0648

for ty, paper in (("GC", "0.367 (0.330 to 0.404)"), ("MC", "0.613 (0.552 to 0.674)"), ("CA3", "0.556 (0.500 to 0.612)")):
    x = P[P.type == ty].remap
    print(f"{ty}: {int((x >= 0.2).sum())} of {len(x)} = {np.mean(x >= 0.2):.3f}; paper {paper}")
for thr in (0.1, 0.15, 0.25, 0.3):
    print(f"Variant threshold {thr}: " + ", ".join(f"{t} {np.mean(P[P.type == t].remap >= thr):.3f}" for t in ("GC", "MC", "CA3")))

GC: 11 of 27 = 0.407; paper 0.367 (0.330 to 0.404)
MC: 65 of 116 = 0.560; paper 0.613 (0.552 to 0.674)
CA3: 73 of 143 = 0.510; paper 0.556 (0.500 to 0.612)
Variant threshold 0.1: GC 0.556, MC 0.664, CA3 0.664
Variant threshold 0.15: GC 0.444, MC 0.612, CA3 0.552
Variant threshold 0.25: GC 0.407, MC 0.517, CA3 0.469
Variant threshold 0.3: GC 0.296, MC 0.491, CA3 0.434

rng = np.random.default_rng(20170201)
pmaps = {ty: [tuple(OMAPS[r.session][k][r.unit][0] for k in ("A", "B", "A2")) for r in P[P.type == ty].itertuples()] for ty in ("GC", "MC", "CA3")}
print("Place cells pooled per type:", {k: len(v) for k, v in pmaps.items()})
for k in (5, 10, 15, 20):
    res = {ty: dg.pv_remap(pmaps[ty], k, rng) for ty in ("GC", "MC", "CA3")}
    kw = stats.kruskal(*[v[np.isfinite(v)] for v in res.values()]).pvalue
    print(f"k = {k:2d}: median PV remapping index GC {np.nanmedian(res['GC']):.3f}, MC {np.nanmedian(res['MC']):.3f}, CA3 {np.nanmedian(res['CA3']):.3f}; "
          f"Kruskal-Wallis P = {kw:.2g}; GC below MC in {np.mean(res['GC'] < np.median(res['MC'])):.2f} of GC draws")
allpv = {ty: dg.pv_remap(pmaps[ty], len(pmaps[ty]), rng, n=1)[0] for ty in ("GC", "MC", "CA3")}
print("All place cells together:", {k: round(v, 3) for k, v in allpv.items()})
prng = np.random.default_rng(20170202)
def pv_index(maps):
    return dg.pv_remap(maps, len(maps), prng, n=1)[0]
pool = pmaps["GC"] + pmaps["MC"]
obs = pv_index(pmaps["MC"]) - pv_index(pmaps["GC"])
null = []
for _ in range(200):
    j = prng.permutation(len(pool))
    null.append(pv_index([pool[i] for i in j[len(pmaps["GC"]):]]) - pv_index([pool[i] for i in j[:len(pmaps["GC"])]]))
null = np.array(null)
print(f"All place cells, MC minus GC PV remapping index {obs:.3f}; label permutation (200): one-sided P = {(1 + np.sum(null >= obs)) / (1 + len(null)):.3f}, "
      f"null 95th percentile {np.percentile(null, 95):.3f}")

Place cells pooled per type: {'GC': 27, 'MC': 116, 'CA3': 143}

k =  5: median PV remapping index GC 0.162, MC 0.253, CA3 0.251; Kruskal-Wallis P = 5.8e-35; GC below MC in 0.71 of GC draws

k = 10: median PV remapping index GC 0.177, MC 0.302, CA3 0.311; Kruskal-Wallis P = 4e-99; GC below MC in 0.84 of GC draws

k = 15: median PV remapping index GC 0.188, MC 0.296, CA3 0.338; Kruskal-Wallis P = 2.7e-149; GC below MC in 0.88 of GC draws

k = 20: median PV remapping index GC 0.210, MC 0.282, CA3 0.352; Kruskal-Wallis P = 8.1e-190; GC below MC in 0.89 of GC draws
All place cells together: {'GC': np.float64(0.22), 'MC': np.float64(0.25), 'CA3': np.float64(0.428)}

All place cells, MC minus GC PV remapping index 0.030; label permutation (200): one-sided P = 0.323, null 95th percentile 0.197

Fi = F.set_index(["session", "unit"])
tm_pairs = pairs[pairs.session.isin(TMAPS.keys())]
rows = []
for r in tm_pairs.itertuples():
    gf = Fi.loc[(r.session, r.gc), "fields"]
    mf = Fi.loc[(r.session, r.mc), "fields"]
    for g, (a, b, _) in gf:
        overlap = any(g == g2 and not (b2 < a or a2 > b) for g2, (a2, b2, _) in mf)
        rows.append({"session": r.session, "gc": r.gc, "mc": r.mc, "overlap": overlap})
IH = pd.DataFrame(rows, columns=["session", "gc", "mc", "overlap"])
print(f"Detected GC to MC pairs in T-maze sessions: {len(tm_pairs)} ({tm_pairs.groupby(['session', 'gc']).ngroups} GC, {tm_pairs.groupby(['session', 'mc']).ngroups} MC)")
print("GC partners without a T-maze field:", sum(len(Fi.loc[(r.session, r.gc), 'fields']) == 0 for r in tm_pairs.drop_duplicates(['session', 'gc']).itertuples()),
      "of", len(tm_pairs.drop_duplicates(['session', 'gc'])))
print("MC fields of the partners:", sum(len(Fi.loc[(r.session, r.mc), 'fields']) for r in tm_pairs.drop_duplicates(['session', 'mc']).itertuples()))
print(f"GC fields (counted per pair): {len(IH)}; overlapping a partner MC field: {int(IH.overlap.sum())} ({IH.overlap.mean() if len(IH) else np.nan:.2f})")
srng = np.random.default_rng(20170203)
mcs_by_session = {k: [u for u in g.unit] for k, g in F[F.type == "MC"].groupby("session")}
null = []
for _ in range(1000):
    k = n = 0
    for r in tm_pairs.itertuples():
        others = [m for m in mcs_by_session.get(r.session, []) if m != r.mc]
        if not others:
            continue
        m2 = others[srng.integers(len(others))]
        mf = Fi.loc[(r.session, m2), "fields"]
        for g, (a, b, _) in Fi.loc[(r.session, r.gc), "fields"]:
            n += 1
            k += any(g == g2 and not (b2 < a or a2 > b) for g2, (a2, b2, _) in mf)
    null.append(k / n if n else np.nan)
null = np.array(null)
print(f"Chance level, partner replaced by a random other MC of the same session (1000 draws): median {np.nanmedian(null):.2f}, "
      f"5th to 95th percentile {np.nanpercentile(null, 5):.2f} to {np.nanpercentile(null, 95):.2f}; observed {IH.overlap.mean():.2f}, "
      f"fraction of draws at or below it {np.nanmean(null <= IH.overlap.mean()):.2f}")

Detected GC to MC pairs in T-maze sessions: 30 (23 GC, 28 MC)
GC partners without a T-maze field: 15 of 23
MC fields of the partners: 112
GC fields (counted per pair): 16; overlapping a partner MC field: 5 (0.31)

Chance level, partner replaced by a random other MC of the same session (1000 draws): median 0.38, 5th to 95th percentile 0.19 to 0.62; observed 0.31, fraction of draws at or below it 0.33

rows = []
for r in pairs.itertuples():
    s = smap[r.session]
    for k, name in dg.OPEN_FIELDS.items():
        if f"{name}_twhl_norm" not in s["positions"]:
            continue
        maps = dg.open_field_maps(s, name, [r.gc, r.mc])
        fields = dg.open_field_fields(maps[r.gc])
        if not fields:
            continue
        mask = np.any(fields, axis=0)
        a, b = maps[r.gc][mask], maps[r.mc][mask]
        ok = np.isfinite(a) & np.isfinite(b)
        if ok.sum() >= 10 and np.std(a[ok]) > 0 and np.std(b[ok]) > 0:
            rows.append({"session": r.session, "gc": r.gc, "mc": r.mc, "field": k, "r": float(np.corrcoef(a[ok], b[ok])[0, 1])})
IO = pd.DataFrame(rows, columns=["session", "gc", "mc", "field", "r"])
print(f"Comparisons (pair x open-field session with a GC field): {len(IO)} from {IO.groupby(['session', 'gc', 'mc']).ngroups if len(IO) else 0} pairs, {IO.session.nunique() if len(IO) else 0} sessions")
if len(IO) >= 5:
    print(f"Median r = {IO.r.median():.3f}; above 0.5: {int((IO.r > 0.5).sum())}; Wilcoxon against 0 P = {stats.wilcoxon(IO.r).pvalue:.2g}")
    print(IO.groupby("field").r.agg(["size", "median"]).round(3).to_string())

Comparisons (pair x open-field session with a GC field): 18 from 6 pairs, 4 sessions
Median r = -0.038; above 0.5: 5; Wilcoxon against 0 P = 0.67
       size  median
field              
A         5  -0.188
A2        4  -0.134
Ac        2   0.360
B         4   0.041
Bc        3  -0.055

# Main: the site read from the wide-band signal (the column where the unit's spike is largest), as the
# electrodes table's channel map disagrees with the wide-band signal (data check). Variant: the table's site.
x = R[R.type.isin(["GC", "MC"]) & (R.raw_n >= 20) & R.ds2_rawsite_uV.notna()]
g, m = x[x.type == "GC"].ds2_rawsite_uV, x[x.type == "MC"].ds2_rawsite_uV
print(f"Units with a wide-band site: GC {len(g)}, MC {len(m)}, in {x.session.nunique()} files of {x.subject.nunique()} mice")
print(f"DS2 amplitude at the unit's wide-band site (uV): GC median {g.median():.0f}, MC median {m.median():.0f}; Mann-Whitney P = {stats.mannwhitneyu(g, m).pvalue:.2g}")
print("Fraction at a site where DS2 is positive (below the reversal): GC %.2f, MC %.2f" % ((g > 0).mean(), (m > 0).mean()))
print("By subject (median uV):")
print(x.groupby(["subject", "type"]).ds2_rawsite_uV.agg(["size", "median"]).round(0).unstack().to_string())
v = R[R.type.isin(["GC", "MC"]) & R.ds2_site_uV.notna()]
gv, mv = v[v.type == "GC"].ds2_site_uV, v[v.type == "MC"].ds2_site_uV
print(f"\nVariant, site from the electrodes table: GC median {gv.median():.0f}, MC median {mv.median():.0f} uV (n = {len(gv)}, {len(mv)}); P = {stats.mannwhitneyu(gv, mv).pvalue:.2g}")
both = x.dropna(subset=["ds2_site_uV"])
print("Units with both sites: %d; same LFP column %.2f; same DS2 sign %.2f" % (
    len(both), np.mean(both.raw_col == both.table_col),
    np.mean(np.sign(both.ds2_site_uV) == np.sign(both.ds2_rawsite_uV))))
print("DS2 events per file used for the amplitudes: median %d, range %d to %d" % (
    np.median([s["unit_site"]["n_events"] for s in SITE]), min(s["unit_site"]["n_events"] for s in SITE), max(s["unit_site"]["n_events"] for s in SITE)))
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist([g, m], bins=np.linspace(-1500, 1500, 41), color=["gold", "m"], label=["GC", "MC"])
ax.set_xlabel("DS2 amplitude at the unit's wide-band site (uV)"); ax.legend()
print("Outside the plotted range: %d units" % int(((x.ds2_rawsite_uV < -1500) | (x.ds2_rawsite_uV > 1500)).sum()))

Units with a wide-band site: GC 139, MC 226, in 17 files of 4 mice
DS2 amplitude at the unit's wide-band site (uV): GC median 232, MC median 580; Mann-Whitney P = 2e-40
Fraction at a site where DS2 is positive (below the reversal): GC 0.80, MC 1.00
By subject (median uV):
            size     median       
type          GC  MC     GC     MC
subject                           
YutaMouse33   10  59   49.0  571.0
YutaMouse41   84  66  214.0  603.0
YutaMouse42   32  71  301.0  602.0
YutaMouse54   13  30  200.0  543.0

Variant, site from the electrodes table: GC median -38, MC median 271 uV (n = 153, 243); P = 5.7e-09
Units with both sites: 365; same LFP column 0.09; same DS2 sign 0.63
DS2 events per file used for the amplitudes: median 121, range 8 to 365
Outside the plotted range: 0 units

from sklearn.cluster import KMeans
from sklearn.decomposition import PCA

def d2_features(wave):
    """Second derivative of the waveform normalised to its trough, upsampled to 100 kHz, from the trough
    to 0.8 ms after it; None when the window runs past the waveform's end."""
    tp, at_end, y = dg.trough_to_peak_ms(wave)
    if y is None or not np.isfinite(tp):
        return None, 0
    i = int(np.argmin(y))
    seg = y[i:i + 80]
    short = 80 - len(seg)
    return np.gradient(np.gradient(np.pad(seg, (0, short), mode="edge"))), short

rows = []
for s in SITE:
    ru = s["raw_units"]["units"]
    ty = s["units"].set_index("unit").type
    for u in s["spikes"]:
        if ty[u] not in ("GC", "MC"):
            continue
        r = {"session": s["session"], "unit": u}
        if u in ru and ru[u]["n"] >= 20:
            r["d2_raw"], r["pad_raw"] = d2_features(ru[u]["wave"][None])
        if u in s["waveforms"]:
            r["d2_stored"], r["pad_stored"] = d2_features(s["waveforms"][u])
        rows.append(r)
WF = pd.DataFrame(rows)
K = WF.merge(R[["session", "unit", "type", "subject", "ds2_site_uV", "ds2_rawsite_uV", "rate_nrem", "rate_awake", "rate_awake_home"]], on=["session", "unit"])
n_all = len(K)

def kmeans_agreement(K, feats):
    z = (K[feats].values - K[feats].values.mean(0)) / K[feats].values.std(0)
    lab = KMeans(n_clusters=2, n_init=50, random_state=0).fit_predict(z)
    truth = (K.type == "MC").values.astype(int)
    flip = np.mean(lab == truth) < 0.5
    return max(np.mean(lab == truth), np.mean(lab != truth)), (1 - lab if flip else lab)

def prepare(K, d2col, ds2col, wake_col):
    k = K[K[d2col].notna() & K[ds2col].notna() & (K[wake_col] > 0)].copy()
    pcs = PCA(n_components=2).fit_transform(np.stack(k[d2col].values))
    k["w_pc1"], k["w_pc2"] = pcs[:, 0], pcs[:, 1]
    k["ratio"] = k.rate_nrem / k[wake_col]
    k["log_ratio"] = np.log10((k.rate_nrem + 1e-3) / k[wake_col])
    k["ds2"] = k[ds2col]
    return k

# Main: wide-band site and unfiltered wide-band waveform, the Methods' NREM / wake ratio with wake as `awake`.
KM = prepare(K, "d2_raw", "ds2_rawsite_uV", "rate_awake")
agree, lab = kmeans_agreement(KM, ["ds2", "ratio", "w_pc1", "w_pc2"])
print(f"GC and MC in the 17 files: {n_all}; used: {len(KM)} (GC {int((KM.type == 'GC').sum())}, MC {int((KM.type == 'MC').sum())}); "
      f"dropped for fewer than 20 spikes in the wide-band window, no DS2 site or a zero awake rate: {n_all - len(KM)}")
print(f"Agreement of the two clusters with the labels: {agree:.3f}")
print(pd.crosstab(KM.type, lab, rownames=["label"], colnames=["cluster (1 = MC-like)"]).to_string())
print("Feature medians by label:")
print(KM.groupby("type")[["ds2", "ratio", "w_pc1", "w_pc2"]].median().round(3).to_string())
print("NREM / wake ratio range: %.2f to %.2f" % (KM.ratio.min(), KM.ratio.max()))
print("By subject, agreement:", KM.assign(ok=(lab == (KM.type == "MC").values.astype(int))).groupby("subject").ok.mean().round(2).to_dict())
print("\nVariants:")
for drop in ("ds2", "ratio", "w_pc1"):
    print(f"  without {drop}: {kmeans_agreement(KM, [c for c in ('ds2', 'ratio', 'w_pc1', 'w_pc2') if c != drop])[0]:.3f}")
print(f"  log10 ratio in place of the ratio: {kmeans_agreement(KM, ['ds2', 'log_ratio', 'w_pc1', 'w_pc2'])[0]:.3f}")
KH = prepare(K, "d2_raw", "ds2_rawsite_uV", "rate_awake_home")
print(f"  wake as home-cage awake (n = {len(KH)}): ratio {kmeans_agreement(KH, ['ds2', 'ratio', 'w_pc1', 'w_pc2'])[0]:.3f}, log ratio {kmeans_agreement(KH, ['ds2', 'log_ratio', 'w_pc1', 'w_pc2'])[0]:.3f}")
KS = prepare(K, "d2_stored", "ds2_site_uV", "rate_awake")
print(f"  stored waveforms and the table's site (n = {len(KS)}; second-derivative window padded past the waveform's end for {int((KS.pad_stored > 0).sum())} units, "
      f"median {KS.pad_stored[KS.pad_stored > 0].median():.0f} of 80 samples): ratio {kmeans_agreement(KS, ['ds2', 'ratio', 'w_pc1', 'w_pc2'])[0]:.3f}, "
      f"log ratio {kmeans_agreement(KS, ['ds2', 'log_ratio', 'w_pc1', 'w_pc2'])[0]:.3f}")
print("  wide-band window padded for %d units" % int((KM.pad_raw > 0).sum()))

GC and MC in the 17 files: 396; used: 364 (GC 138, MC 226); dropped for fewer than 20 spikes in the wide-band window, no DS2 site or a zero awake rate: 32
Agreement of the two clusters with the labels: 0.953
cluster (1 = MC-like)    0    1
label                          
GC                     130    8
MC                       9  217
Feature medians by label:
          ds2  ratio  w_pc1  w_pc2
type                              
GC    232.004  4.203  0.008  0.001
MC    579.593  1.073 -0.006  0.000
NREM / wake ratio range: 0.09 to 36.92
By subject, agreement: {'YutaMouse33': 0.96, 'YutaMouse41': 0.95, 'YutaMouse42': 0.97, 'YutaMouse54': 0.93}

Variants:
  without ds2: 0.890
  without ratio: 0.942
  without w_pc1: 0.885
  log10 ratio in place of the ratio: 0.956

  wake as home-cage awake (n = 364): ratio 0.948, log ratio 0.951
  stored waveforms and the table's site (n = 391; second-derivative window padded past the waveform's end for 391 units, median 5 of 80 samples): ratio 0.928, log ratio 0.944
  wide-band window padded for 1 units

print("Finite x, y or z values in the electrodes tables of all files:", sum(s["electrode_xyz_finite"] for s in S))
print("Files whose units name their shank:", sum(s["format"] == "electrode_group" for s in S), "; of which with waveforms:", len(SITE))
print("shank_id range over all units:", int(U.shank_id.min()), "to", int(U.shank_id.max()))
print("Wide-band check (data check): largest column on the shank the table gives for %.2f of units; DS2 profiles crossing zero more than once along the site numbers: %d of %d shanks with a crossing" % (
    MP.raw_on_table_shank.mean(), int((prof.sign_changes > 1).sum()), int((prof.sign_changes > 0).sum())))
print("Recovered stored-channel order (data check): DS2 profiles crossing zero once %d, more than once %d, of %d shanks that cross" % (
    int((prof2.sign_changes == 1).sum()), int((prof2.sign_changes > 1).sum()), int((prof2.sign_changes > 0).sum())))
print(f"Total run time: {(time.time() - notebook_started) / 60:.1f} min")

Finite x, y or z values in the electrodes tables of all files: 0
Files whose units name their shank: 18 ; of which with waveforms: 17
shank_id range over all units: 0 to 38
Wide-band check (data check): largest column on the shank the table gives for 0.45 of units; DS2 profiles crossing zero more than once along the site numbers: 43 of 64 shanks with a crossing
Recovered stored-channel order (data check): DS2 profiles crossing zero once 43, more than once 8, of 51 shanks that cross
Total run time: 8.0 min

Claim	Status
unit-counts	tested
first-stage-labels	tested
interneuron-rates	tested
gc-narrower-waveform	tested
optotag-drd2	tested
mc-more-active	tested
gc-nrem-highest	tested
ca3-wake-highest	tested
nrem-wake-ratio	tested
theta-modulation	tested
gamma-modulation	tested
gamma-rate-correlation	tested
ds2-recruitment	tested
monosyn-gc-to-mc	tested
burst-facilitation	tested (depends on monosyn-gc-to-mc)
state-transmission	tested (depends on monosyn-gc-to-mc)
zero-lag-synchrony	tested
tmaze-central-active	tested
tmaze-multiple-fields	tested
tmaze-peak-info	tested
tmaze-splitters	tested
remap-positive	tested
remap-gc-lt-mc	tested
remap-fractions	tested
pv-remap	tested
inheritance-tmaze	tested (depends on monosyn-gc-to-mc, tmaze-multiple-fields)
inheritance-openfield	tested (depends on monosyn-gc-to-mc)
mc-deeper	tested on the 17 files that give each unit's shank
gc-mc-classification	tested on the 17 files that give each unit's shank
gc-mc-distance	blocked: no recording site per unit
gc-backpropagation	blocked: no recording site per unit

Claim	Outcome
unit-counts	partially reproduced
first-stage-labels	partially reproduced
interneuron-rates	partially reproduced
gc-narrower-waveform	reproduced
optotag-drd2	partially reproduced
mc-more-active	reproduced
gc-nrem-highest	partially reproduced
ca3-wake-highest	not reproduced
nrem-wake-ratio	reproduced
theta-modulation	reproduced
gamma-modulation	reproduced
gamma-rate-correlation	partially reproduced
ds2-recruitment	reproduced
monosyn-gc-to-mc	partially reproduced
burst-facilitation	reproduced
state-transmission	not reproduced
zero-lag-synchrony	partially reproduced
tmaze-central-active	partially reproduced
tmaze-multiple-fields	partially reproduced
tmaze-peak-info	reproduced
tmaze-splitters	reproduced
remap-positive	reproduced
remap-gc-lt-mc	not reproduced
remap-fractions	partially reproduced
pv-remap	reproduced
inheritance-tmaze	reproduced
inheritance-openfield	reproduced
mc-deeper	reproduced
gc-mc-classification	reproduced
gc-mc-distance	blocked: no site coordinates
gc-backpropagation	blocked: no site coordinates

Reproducing Senzai & Buzsaki (2017) from dandiset 000003¶

Data check¶

Unit measures¶

unit-counts¶

first-stage-labels¶

interneuron-rates¶

gc-narrower-waveform¶

optotag-drd2¶

mc-more-active¶

gc-nrem-highest¶

ca3-wake-highest¶

nrem-wake-ratio¶

theta-modulation¶

gamma-modulation¶

gamma-rate-correlation¶

ds2-recruitment¶

monosyn-gc-to-mc¶

burst-facilitation¶

state-transmission¶

zero-lag-synchrony¶

tmaze-central-active¶

tmaze-multiple-fields¶

tmaze-peak-info¶

tmaze-splitters¶

remap-positive¶

remap-gc-lt-mc¶

remap-fractions¶

pv-remap¶

inheritance-tmaze¶

inheritance-openfield¶

mc-deeper¶

gc-mc-classification¶

Blocked claims¶

Summary¶