import io
import time
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
import scipy.io as sio
import statsmodels.api as sm
from dandi.dandiapi import DandiAPIClient
from scipy import stats
from scipy.ndimage import gaussian_filter1d
from statsmodels.formula.api import ols

import gpivla as g

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

NB_START = started = time.time()
all_files = g.asset_urls(pre_only=False)
sessions = g.load_all()
print(f"{len(all_files)} files in dandiset {g.DANDISET_ID} version {g.DANDISET_VERSION}; "
      f"loaded the {len(sessions)} pre-MPTP files in {time.time() - started:.0f} s")

378 files in dandiset 000947 version 0.240510.2211; loaded the 254 pre-MPTP files in 298 s

with DandiAPIClient() as client:
    meta = client.get_dandiset(g.DANDISET_ID, g.DANDISET_VERSION).get_raw_metadata()
summary = meta["assetsSummary"]
print("Archive summary: subjects", summary["numberOfSubjects"], "| files", summary["numberOfFiles"],
      "| species", [s["name"] for s in summary["species"]])
print("Related resources:", [r.get("identifier") or r.get("url") for r in meta.get("relatedResource", [])])
paths = pd.DataFrame({"path": [p for p, _ in all_files]})
paths["subject"] = paths.path.str.extract(r"sub-(\w+)/")[0]
paths["condition"] = np.where(paths.path.str.contains("ses-pre-MPTP"), "pre-MPTP", "post-MPTP")
print("\nFiles by subject and condition:")
print(pd.crosstab(paths.subject, paths.condition).to_string())

info = pd.DataFrame([{
    "session": s["session_id"], "key": s["key"], "subject": s["subject"], "mptp_attr": s["mptp"],
    "start": s["session_start"][:10], "has_trials": s["trials"] is not None,
    "n_trials": 0 if s["trials"] is None else len(s["trials"]), "n_units": len(s["units"]),
    "publications": ";".join(s["related_publications"]), "n_raw": s["n_raw"], "groups": str(sorted(s["groups"].values()))}
    for s in sessions])
print("\nMPTP status attribute of the pre-MPTP files:", info.mptp_attr.value_counts().to_dict())
print("related_publications in the files:", info.publications.value_counts().to_dict())
print("Files without a trials table (session, curated units, raw samples):",
      [(r.session, r.n_units, r.n_raw) for r in info[~info.has_trials].itertuples()])
print("Files without curated units:", int((info.n_units == 0).sum()), "of", len(info))
pu_rows = []
for s in sessions:
    if s["units"]:
        continue
    f = g._open(s["url"])
    pu = f["processing/ecephys/units"] if "units" in f.get("processing/ecephys", {}) else None
    pu_rows.append({"session": s["session_id"], "sorted_clusters": 0 if pu is None else len(pu["id"]),
                    "colnames": "(no table)" if pu is None else ", ".join(g._s(c) for c in pu.attrs["colnames"])})
    f.close()
pu_rows = pd.DataFrame(pu_rows)
print("  their sorted clusters in processing/ecephys/units: files with clusters", int((pu_rows.sorted_clusters > 0).sum()),
      "; clusters", int(pu_rows.sorted_clusters.sum()), "; columns:", pu_rows.colnames.value_counts().to_dict())
print(pu_rows[pu_rows.session.str.contains("G-140813-2")].to_string(index=False))
print("Trials per file: min", info.n_trials[info.has_trials].min(), "median",
      info.n_trials[info.has_trials].median(), "max", info.n_trials.max())
day = pd.to_datetime(info.session.str.extract(r"-(\d{6})-")[0], format="%y%m%d")
offset = (pd.to_datetime(info.start) - day).dt.days
print("session_start_time minus the date in session_id (days):", offset.value_counts().to_dict())
print("Electrode group locations (files):")
print(info.groups.value_counts().head(12).to_string())

Archive summary: subjects 2 | files 378 | species ['Macaca mulatta - Rhesus monkey']
Related resources: ['https://doi.org/10.1371/journal.pbio.3000829', 'https://github.com/dandi/example-notebooks/blob/master/000947/TurnerLab/public_demo/000947_demo.ipynb', 'https://doi.org/10.1101/2023.09.08.556120']

Files by subject and condition:
condition  post-MPTP  pre-MPTP
subject                       
Gaia             101       119
Isis              23       135

MPTP status attribute of the pre-MPTP files: {'pre_MPTP': 254}
related_publications in the files: {'https://doi.org/10.1101/546598': 254}
Files without a trials table (session, curated units, raw samples): [('pre-MPTP-I-160616-4', 5, 10330112)]
Files without curated units: 21 of 254

  their sorted clusters in processing/ecephys/units: files with clusters 13 ; clusters 124 ; columns: {'spike_times, unit_name, location': 13, '(no table)': 8}
            session  sorted_clusters                         colnames
pre-MPTP-G-140813-2               16 spike_times, unit_name, location
Trials per file: min 11 median 103.0 max 215
session_start_time minus the date in session_id (days): {2: 254}
Electrode group locations (files):
groups
['GPi', 'VLa']             79
['GPi', 'VLp']             65
['VLa']                    21
['M1-G', 'M1-S']           18
['M1-S', 'SMA']            12
['GPi', 'Thalamus']        10
['VLp']                     7
['GP', 'VL']                5
['SMA']                     5
['M1-G', 'M1-S', 'SMA']     5
['M1-S']                    4
['GPi', 'VMb']              3

rows = []
for s in sessions:
    span = None
    for u in s["units"]:
        sp = u["spikes"]
        rows.append({"session": s["session_id"], "key": s["key"], "subject": s["subject"],
                     "unit_name": u.get("unit_name"), "label": u.get("location"),
                     "area": g.area(u.get("location")), "channel": pd.to_numeric(u.get("channel_ids"), errors="coerce"),
                     "sort": str(u.get("sort_label")), "quality": u.get("unit_quality_post_sorting"),
                     "n_spikes": len(sp), "first": sp[0] if len(sp) else np.nan, "last": sp[-1] if len(sp) else np.nan,
                     "sorted": u["sorted_in_file"], "rec_s": s["n_raw"] / 24414.0625})
U = pd.DataFrame(rows)
print("Curated units by subject and location label as written in the files:")
print(pd.crosstab(U.label.fillna("(none)").replace("", "(empty)"), U.subject).to_string())
print("\nUnits with no channel id:", int(U.channel.isna().sum()), "; with an empty location:", int((U.label.fillna("") == "").sum()))
print("Sort labels:", U["sort"].value_counts().to_dict())
print("Quality labels (GPi and VLa):", U[U.area.isin(["GPi", "VLa"])].quality.fillna("(column absent)").value_counts().to_dict())
print("\nSpike times stored in increasing order in every unit:", bool(U.sorted.all()))
print("Spike times within the raw recording (0 to its length):",
      bool(((U["first"] >= 0) & (U["last"] <= U.rec_s + 1)).all()))
# Recordings deposited twice: identical spike trains within or across files.
digest = {}
for s in sessions:
    for u in s["units"]:
        if len(u["spikes"]):
            digest.setdefault((len(u["spikes"]), round(float(u["spikes"][0]), 6), round(float(u["spikes"][-1]), 6)), []).append(
                (s["session_id"], u.get("unit_name")))
dups = {k: v for k, v in digest.items() if len(v) > 1}
print("Spike trains with the same count, first and last spike in more than one unit:", len(dups))
for k, v in list(dups.items())[:10]:
    us = [u for s in sessions if s["session_id"] == v[0][0] for u in s["units"] if u.get("unit_name") in [x[1] for x in v]]
    print("  ", v, "labels", [u.get("location") for u in us], "channels", [u.get("channel_ids") for u in us],
          "identical:", all(np.array_equal(us[0]["spikes"], x["spikes"]) for x in us[1:]))
same_name = U.groupby(["session", "unit_name"]).size()
print("Unit names used twice within a file:", int((same_name > 1).sum()))
grid = []
for s in sessions:
    sp = np.concatenate([u["spikes"] for u in s["units"]]) if s["units"] else np.array([])
    if len(sp):
        grid.append({"sample": np.mean(np.abs(sp * 24414.0625 - np.round(sp * 24414.0625)) < 1e-3),
                     "whole_ms": np.mean(np.abs(sp * 1000 - np.round(sp * 1000)) < 1e-6)})
grid = pd.DataFrame(grid)
print(f"Fraction of spike times on the 24414.0625 Hz sample grid: lowest over {len(grid)} files {grid['sample'].min():.4f}; "
      f"on whole ms: highest over files {grid.whole_ms.max():.4f}")
print("Files where a channel id is given to more than one electrode:",
      sum(bool(s["electrodes"].channel_ids.duplicated().any()) for s in sessions))
print("Electrode table columns read:", list(sessions[0]["electrodes"].columns))

Curated units by subject and location label as written in the files:
subject   Gaia  Isis
label               
(empty)      0    33
GPE          1     0
GPI         52     0
GPi         99   212
M1-G         0    98
M1-S         0   145
R-V          0    12
RTN          0    14
S-R          0     2
SMA          0    77
STR          0     2
Thalamus    54     0
VLA         66     0
VLP          8     0
VLa         19   219
VLp          6   118

Units with no channel id: 33 ; with an empty location: 33
Sort labels: {'a': 710, '1.0': 295, 'b': 130, '2.0': 48, '': 33, 'c': 11, '3.0': 7, '4.0': 3}
Quality labels (GPi and VLa): {'excellent': 353, '(column absent)': 179, 'F': 35, 'good': 25, 'changed to excellent from good based on post-sorting quality': 23, 'B -> C': 16, 'C': 10, 'B -> B': 10, 'no quality': 6, 'B -> F': 4, 'changed to good from excellent based on post-sorting quality': 3, '?': 1, 'C -> A': 1, 'M': 1}

Spike times stored in increasing order in every unit: True
Spike times within the raw recording (0 to its length): True
Spike trains with the same count, first and last spike in more than one unit: 1
   [('pre-MPTP-I-170518-2', '6'), ('pre-MPTP-I-170518-2', '7')] labels ['M1-S', 'M1-S'] channels [np.float64(40.0), np.float64(43.0)] identical: True
Unit names used twice within a file: 33
Fraction of spike times on the 24414.0625 Hz sample grid: lowest over 233 files 1.0000; on whole ms: highest over files 0.0009
Files where a channel id is given to more than one electrode: 0
Electrode table columns read: ['channel_ids', 'location', 'group_name', 'acx_x', 'acx_y', 'acx_z']

for s in sessions:
    s["vt"] = None if s["trials"] is None else g.valid_trials(s["trials"])
tr = pd.concat([s["vt"].assign(session=s["session_id"], subject=s["subject"]) for s in sessions if s["vt"] is not None])
print("Trial table columns (files):")
print(pd.Series([tuple(s["trials"].columns) for s in sessions if s["trials"] is not None]).value_counts().to_string())
print("\nTrials by subject: total, errors, outliers, valid")
print(tr.groupby("subject").agg(total=("valid", "size"), errors=("error", "sum"),
                                outliers=("outlier", "sum"), valid=("valid", "sum")).to_string())
print("Targets:", tr.target.value_counts().to_dict())
v = tr[tr.valid]
print("Reaction time of valid trials (s): min %.3f median %.3f max %.3f" % (v.rt.min(), v.rt.median(), v.rt.max()))
print("Movement duration of valid trials (s): min %.3f median %.3f max %.3f; trials without one: %d of %d"
      % (v.mt.min(), v.mt.median(), v.mt.max(), v.mt.isna().sum(), len(v)))
order_ok = (tr.cue_onset_time <= tr.movement_start_time) | ~tr.valid
print("Valid trials with cue onset before movement onset:", bool(order_ok.all()))
print("Files with movement ends, by subject:",
      {k: int(x) for k, x in pd.Series({s["session_id"]: ("movement_stop_time" in s["trials"]) for s in sessions if s["trials"] is not None}).groupby(lambda k: k.split("-")[2]).sum().items()})

Trial table columns (files):
(start_time, stop_time, error_onset_time, reward_start_time, reward_stop_time, movement_start_time, movement_stop_time, return_start_time, return_stop_time, cue_onset_time, target)    84
(start_time, stop_time, reward_start_time, reward_stop_time, movement_start_time, movement_stop_time, return_start_time, return_stop_time, cue_onset_time, target)                      82
(start_time, stop_time, reward_start_time, reward_stop_time, movement_start_time, cue_onset_time, target)                                                                               44
(start_time, stop_time, error_onset_time, reward_start_time, reward_stop_time, movement_start_time, cue_onset_time, target)                                                             43

Trials by subject: total, errors, outliers, valid
         total  errors  outliers  valid
subject                                
Gaia      9676     273       114   9289
Isis     15470     540       591  14061
Targets: {'Left': 12584, 'Right': 12562}
Reaction time of valid trials (s): min 0.104 median 0.239 max 0.529
Movement duration of valid trials (s): min 0.136 median 0.196 max 0.422; trials without one: 6131 of 23350
Valid trials with cue onset before movement onset: True
Files with movement ends, by subject: {'G': 32, 'I': 134}

def window_counts(spikes, align, lo, hi):
    a = np.searchsorted(spikes, align + lo)
    b = np.searchsorted(spikes, align + hi)
    return b - a

rng = np.random.default_rng(0)
rows = []
for s in sessions:
    if s["vt"] is None:
        continue
    t = s["vt"][s["vt"].valid]
    left = (t.target == "Left").to_numpy()
    permuted = rng.permutation(left)
    al = t.movement_start_time.to_numpy()
    for u in s["units"]:
        if g.area(u.get("location")) not in ("GPi", "VLa") or left.sum() < 10 or (~left).sum() < 10:
            continue
        r = {"session": s["session_id"], "subject": s["subject"], "area": g.area(u.get("location"))}
        for name, (lo, hi) in {"movement": (-0.2, 0.3), "pre_cue": (-1.5, -1.0)}.items():
            for lab, mask, shift in (("true", left, 0.0), ("permuted", permuted, 0.0), ("offset_1s", left, 1.0)):
                c = window_counts(u["spikes"], al + shift, lo, hi)
                r[f"{name}_{lab}"] = stats.mannwhitneyu(c[mask], c[~mask]).pvalue if c.std() > 0 else 1.0
        rows.append(r)
clock = pd.DataFrame(rows)
cols = [f"{w}_{l}" for w in ("movement", "pre_cue") for l in ("true", "permuted", "offset_1s")]
print("Fraction of units whose spike counts differ between left and right trials (p < 0.05):")
print(clock.groupby("area")[cols].agg(lambda p: (p < 0.05).mean()).round(3).T.to_string())
per_session = clock.groupby("session").movement_true.apply(lambda p: (p < 0.05).any())
print(f"\nFiles with at least one direction-tuned unit in the movement window: {per_session.sum()} of {len(per_session)}")
print("Files where none is tuned (the clock cannot be checked this way):", sorted(per_session[~per_session].index))

Fraction of units whose spike counts differ between left and right trials (p < 0.05):
area                  GPi    VLa
movement_true       0.485  0.386
movement_permuted   0.067  0.034
movement_offset_1s  0.372  0.319
pre_cue_true        0.087  0.078
pre_cue_permuted    0.038  0.041
pre_cue_offset_1s   0.253  0.214

Files with at least one direction-tuned unit in the movement window: 129 of 163
Files where none is tuned (the clock cannot be checked this way): ['pre-MPTP-G-140511-3', 'pre-MPTP-G-140511-6', 'pre-MPTP-G-140511-8', 'pre-MPTP-G-140511-9', 'pre-MPTP-G-140512-1', 'pre-MPTP-G-140512-8', 'pre-MPTP-G-140512-9', 'pre-MPTP-G-140514-10', 'pre-MPTP-G-140514-8', 'pre-MPTP-G-140521-3', 'pre-MPTP-G-140524-1', 'pre-MPTP-G-140524-3', 'pre-MPTP-G-140524-5', 'pre-MPTP-G-140524-6', 'pre-MPTP-G-140525-6', 'pre-MPTP-G-140611-1', 'pre-MPTP-G-140611-3', 'pre-MPTP-G-140612-5', 'pre-MPTP-G-140612-9', 'pre-MPTP-G-140617-3', 'pre-MPTP-G-140617-4', 'pre-MPTP-G-140627-2', 'pre-MPTP-G-140701-1', 'pre-MPTP-G-140723-1', 'pre-MPTP-G-140723-2', 'pre-MPTP-G-140902-3', 'pre-MPTP-G-140902-4', 'pre-MPTP-I-160624-1', 'pre-MPTP-I-160714-1', 'pre-MPTP-I-160714-4', 'pre-MPTP-I-160718-1', 'pre-MPTP-I-170615-2', 'pre-MPTP-I-170616-2', 'pre-MPTP-I-170623-2']

started = time.time()
def run_units(records, alpha=g.ALPHA):
    """records: list of dicts with unit id, area, subject and segments [(spikes, valid trials)]."""
    out = []
    for r in records:
        for target in ("Left", "Right"):
            d = g.unit_direction(r["segments"], target, alpha)
            if d is None:
                continue
            out.append({"uid": r["uid"], "area": r["area"], "subject": r["subject"], "target": target,
                        **{k: d[k] for k in ("n_trials", "median_rt_ms", "rest_rate", "integrated", "inc_onset",
                                             "dec_onset", "slope_per_s2", "slope_sig")},
                        "sdf": d["sdf"], "counts": d["trial_counts"]})
    df = pd.DataFrame(out)
    df["form"] = [g.response_form(i, d) for i, d in zip(df.inc_onset, df.dec_onset)]
    return df

dand_records, dropped, rate_rule = [], {"rate": 0, "trials": 0}, []
for s in sessions:
    if s["vt"] is None:
        continue
    for k, u in enumerate(s["units"]):
        a = g.area(u.get("location"))
        if a not in ("GPi", "VLa") or len(u["spikes"]) < 2:
            continue
        rate = len(u["spikes"]) / (u["spikes"][-1] - u["spikes"][0])
        rec_s = s["n_raw"] / 24414.0625
        rate_rule.append({"area": a, "span": rate >= g.MIN_RATE[a],
                          "recording": "no raw data" if rec_s == 0 else len(u["spikes"]) / rec_s >= g.MIN_RATE[a]})
        if rate < g.MIN_RATE[a]:
            dropped["rate"] += 1
            continue
        dand_records.append({"uid": f"{s['session_id']}/{k}", "area": a, "subject": s["subject"], "session": s["session_id"],
                             "key": s["key"], "channel": pd.to_numeric(u.get("channel_ids"), errors="coerce"),
                             "sort": str(u.get("sort_label")), "rate": rate, "segments": [(u["spikes"], s["vt"])],
                             "unit": u, "sess": s})
ud = run_units(dand_records)
kept = set(ud.uid)
dropped["trials"] = sum(r["uid"] not in kept for r in dand_records)
units = pd.DataFrame([{k: r[k] for k in ("uid", "area", "subject", "session", "key", "channel", "sort", "rate")}
                      for r in dand_records if r["uid"] in kept])
print(f"Dandiset rule: {len(units)} units, {len(ud)} unit-directions ({time.time() - started:.0f} s)")
print("Left out: below the rate criterion", dropped["rate"], "; fewer than 10 valid trials in both directions", dropped["trials"])
print(units.groupby(["area", "subject"]).size().unstack().to_string())
print("Variant, the rate over the length of the raw recording instead of the span of the spike train; units passing "
      "(rows: span rule, columns: recording rule):")
print(pd.crosstab([pd.DataFrame(rate_rule).area, pd.DataFrame(rate_rule).span], pd.DataFrame(rate_rule).recording.astype(str)).to_string())

Dandiset rule: 547 units, 1083 unit-directions (7 s)
Left out: below the rate criterion 107 ; fewer than 10 valid trials in both directions 8
subject  Gaia  Isis
area               
GPi       137   131
VLa        78   201
Variant, the rate over the length of the raw recording instead of the span of the spike train; units passing (rows: span rule, columns: recording rule):
recording   False  True
area span              
GPi  False     84     0
     True       2   272
VLa  False     23     0
     True       0   281

started = time.time()
ZENODO = "https://zenodo.org/records/4000719/files/{}?download=1"
analy = {}
for monk in ("Gaia", "Isis"):
    r = requests.get(ZENODO.format(f"{monk}_Analy.mat"), timeout=600)
    r.raise_for_status()
    m = sio.loadmat(io.BytesIO(r.content), squeeze_me=True, struct_as_record=False)
    analy[monk] = (m["RespSDFIncr"], m["RespISIIncr"])

def sort_letter(x):
    return {"1.0": "a", "2.0": "b", "3.0": "c", "4.0": "d"}.get(str(x), str(x))

by_key = {s["key"]: s for s in sessions}
paper_rows, paper_records = [], []
for monk, (A, B) in analy.items():
    for i in range(len(A.fname)):
        utype = str(A.unit_type[i])
        files = str(A.fname[i]).split("/")
        keep = utype != "VLp" and not np.isnan(A.evt_sdf[i, :, 0]).any()
        found = []
        for f in files:
            s = by_key.get(f)
            hits = [] if s is None else [u for u in s["units"]
                    if pd.to_numeric(u.get("channel_ids"), errors="coerce") == float(A.Uchan[i])
                    and sort_letter(u.get("sort_label")) == A.Usort[i]]
            found.append((s, hits[0]) if len(hits) == 1 else None)
        complete = all(x is not None for x in found)
        paper_rows.append({"monk": monk, "i": i, "fname": str(A.fname[i]), "area": utype, "keep": keep,
                           "n_files": len(files), "n_found": sum(x is not None for x in found), "complete": complete,
                           "labels": sorted({g.area(x[1].get("location")) or "(empty)" for x in found if x})})
        if keep and complete:
            paper_records.append({"uid": f"paper/{monk}/{i}", "area": utype, "subject": monk, "monk": monk, "i": i,
                                  "segments": [(u["spikes"], s["vt"]) for s, u in found],
                                  "sessions": [s["session_id"] for s, u in found]})
paper = pd.DataFrame(paper_rows)
ud_paper = run_units(paper_records)
print(f"Paper's records: {len(paper)}; kept by its code: {int(paper.keep.sum())}; found in every file: "
      f"{int((paper.keep & paper.complete).sum())}; unit-directions analysed: {len(ud_paper)} ({time.time() - started:.0f} s)")

Paper's records: 406; kept by its code: 394; found in every file: 364; unit-directions analysed: 724 (68 s)

def same_sort(x):
    t = str(x)
    try:
        return "abcd"[int(float(t)) - 1]
    except (ValueError, IndexError):
        return t

groups = {}
for r in dand_records:
    if r["uid"] not in kept:
        continue
    day = r["session"].split("-")[-2]
    key = (r["subject"], day, r["area"], r["channel"], same_sort(r["sort"])) if np.isfinite(r["channel"]) else (r["uid"],)
    groups.setdefault(key, []).append(r)
merged_records, uid2neuron = [], {}
for n, (key, rs) in enumerate(groups.items()):
    rs = sorted(rs, key=lambda r: r["session"])
    for r in rs:
        uid2neuron[r["uid"]] = f"neuron/{n}"
    merged_records.append({"uid": f"neuron/{n}", "area": rs[0]["area"], "subject": rs[0]["subject"], "key": key,
                           "segments": [seg for r in rs for seg in r["segments"]], "files": tuple(sorted(r["key"] for r in rs))})
ud_m = run_units(merged_records)
ud["neuron"] = ud.uid.map(uid2neuron)
print(f"Presumed neurons: {len(merged_records)} {ud_m.groupby('area').uid.nunique().to_dict()}; unit-directions analysed: {len(ud_m)}")
print("Files per presumed neuron:", pd.Series([len(r["files"]) for r in merged_records]).value_counts().sort_index().to_dict())
paper_sets = {tuple(sorted(str(x).split("/"))) for x in paper[paper.keep].fname}
multi = [r for r in merged_records if len(r["files"]) > 1]
print(f"Distinct file sets of the paper's records spanning several files: {sum(len(x) > 1 for x in paper_sets)}; "
      f"presumed neurons spanning several files: {len(multi)}, with exactly the file set of a paper record: "
      f"{sum(r['files'] in paper_sets for r in multi)}")
prec = [(set(str(p.fname).split("/")), float(analy[p.monk][0].Uchan[p.i]), str(analy[p.monk][0].Usort[p.i]))
        for p in paper[paper.keep].itertuples()]
n_rec = [sum(ch == r["key"][3] and srt == r["key"][4] and fs <= set(r["files"]) for fs, ch, srt in prec) for r in multi]
print("Presumed neurons spanning several files that hold more than one of the paper's records (merged too far):",
      sum(n > 1 for n in n_rec), "of", len(multi))

Presumed neurons: 318 {'GPi': 132, 'VLa': 186}; unit-directions analysed: 632
Files per presumed neuron: {1: 201, 2: 64, 3: 27, 4: 11, 5: 4, 6: 7, 7: 2, 8: 1, 9: 1}
Distinct file sets of the paper's records spanning several files: 53; presumed neurons spanning several files: 117, with exactly the file set of a paper record: 53
Presumed neurons spanning several files that hold more than one of the paper's records (merged too far): 47 of 117

started = time.time()
rng = np.random.default_rng(11)

def random_alignment(vt, spikes):
    t = vt.copy()
    t["movement_start_time"] = rng.uniform(spikes[0] + 1.6, spikes[-1] - 1.1, len(t))
    return t

nul = pd.concat([run_units([{**r, "segments": [(sp, random_alignment(vt, sp)) for sp, vt in r["segments"]]}
                            for r in dand_records if r["uid"] in kept]).assign(rep=rep) for rep in range(3)], ignore_index=True)

def detection_summary(df, key):
    per_dir = df.groupby("area").agg(increase=("inc_onset", lambda x: x.notna().mean()),
                                     decrease=("dec_onset", lambda x: x.notna().mean()))
    per_unit = (df.assign(resp=df.form != "none").groupby(key).agg(responsive=("resp", "any"), slope_sig=("slope_sig", "any"))
                .groupby("area").mean())
    return pd.concat([per_dir, per_unit], axis=1).round(3)

print(f"Random alignment, 3 draws ({time.time() - started:.0f} s). Columns: unit-directions with an increase, with a "
      "decrease; units responsive, with a significant control slope, in at least one direction.")
print("Data, dandiset rule:"); print(detection_summary(ud, ["uid", "area"]).to_string())
print("Chance, mean of 3 draws:"); print(detection_summary(nul, ["rep", "uid", "area"]).to_string())
print("Responsive units by draw:", nul.assign(resp=nul.form != "none").groupby(["rep", "uid", "area"]).resp.any()
      .groupby(["rep", "area"]).mean().unstack().round(3).to_dict("list"))
for k in ("inc", "dec"):
    print(f"Chance {k}rease onsets (ms): median and count by area",
          nul.groupby("area")[f"{k}_onset"].agg(["median", "count"]).to_dict("index"))
print("Chance forms, percent of unit-directions:")
print((100 * pd.crosstab(nul.area, nul.form, normalize="index")).round(1).to_string())

Random alignment, 3 draws (22 s). Columns: unit-directions with an increase, with a decrease; units responsive, with a significant control slope, in at least one direction.
Data, dandiset rule:
      increase  decrease  responsive  slope_sig
area                                           
GPi      0.710     0.497       0.985      0.985
VLa      0.668     0.442       0.975      0.986
Chance, mean of 3 draws:
      increase  decrease  responsive  slope_sig
area                                           
GPi      0.148     0.252       0.623      0.984
VLa      0.162     0.337       0.708      0.990
Responsive units by draw: {'GPi': [0.597, 0.646, 0.627], 'VLa': [0.706, 0.713, 0.706]}
Chance increase onsets (ms): median and count by area {'GPi': {'median': 83.0, 'count': 236}, 'VLa': {'median': 70.0, 'count': 268}}
Chance decrease onsets (ms): median and count by area {'GPi': {'median': 52.0, 'count': 402}, 'VLa': {'median': -8.0, 'count': 558}}
Chance forms, percent of unit-directions:
form  dec_first  decrease  inc_first  increase  none
area                                                
GPi         0.5      24.1        0.6      13.7  61.1
VLa         1.1      31.6        1.0      14.1  52.2

started = time.time()
rng = np.random.default_rng(12)
STRICT = 0.0001
ud_strict = run_units([r for r in dand_records if r["uid"] in kept], alpha=STRICT)
nul_strict = run_units([{**r, "segments": [(sp, random_alignment(vt, sp)) for sp, vt in r["segments"]]}
                        for r in dand_records if r["uid"] in kept], alpha=STRICT)
print(f"Detection at alpha {STRICT} ({time.time() - started:.0f} s). Columns as above.")
print("Data:"); print(detection_summary(ud_strict, ["uid", "area"]).to_string())
print("Chance, one draw:"); print(detection_summary(nul_strict, ["uid", "area"]).to_string())
for k in ("inc", "dec"):
    print(f"Chance {k}rease onsets at alpha {STRICT} (ms): median and count by area",
          nul_strict.groupby("area")[f"{k}_onset"].agg(["median", "count"]).to_dict("index"))

Detection at alpha 0.0001 (15 s). Columns as above.
Data:
      increase  decrease  responsive  slope_sig
area                                           
GPi      0.678     0.475       0.981      0.985
VLa      0.634     0.420       0.953      0.986
Chance, one draw:
      increase  decrease  responsive  slope_sig
area                                           
GPi      0.115     0.196       0.519      0.985
VLa      0.147     0.308       0.667      0.993
Chance increase onsets at alpha 0.0001 (ms): median and count by area {'GPi': {'median': 126.0, 'count': 61}, 'VLa': {'median': 142.0, 'count': 81}}
Chance decrease onsets at alpha 0.0001 (ms): median and count by area {'GPi': {'median': 98.0, 'count': 104}, 'VLa': {'median': 8.5, 'count': 170}}

k = paper[paper.keep]
print("Records kept by the paper's code, by animal and area:")
print(k.groupby(["monk", "area"]).size().unstack().to_string())
print("\nFound in every file they span:")
print(k.groupby(["monk", "area"]).complete.agg(["sum", "size", "mean"]).round(3).to_string())
print(f"Overall: {int(k.complete.sum())} of {len(k)} ({k.complete.mean():.1%}); records spanning 2 to 5 files: {int((k.n_files > 1).sum())}")
missing = k[~k.complete]
print("\nRecords not found, by file:")
print(missing.groupby("fname").size().to_string())
print("\nDandiset labels of the matched records, by the paper's area:")
print(pd.crosstab(k[k.complete].area, k[k.complete].labels.map(lambda l: "/".join(l))).to_string())

rows = []
for _, r in ud_paper.iterrows():
    monk, i = r.uid.split("/")[1], int(r.uid.split("/")[2])
    A, B = analy[monk]
    t = 0 if r.target == "Left" else 1
    ref = A.evt_sdf[i, t, :]
    rows.append({"area": r.area, "r": np.corrcoef(r.sdf, ref)[0, 1], "mean_ratio": r.sdf.mean() / ref.mean(),
                 "n": r.n_trials, "n_ref": np.shape(A.evt_spkt[i, t])[0],
                 "inc": r.inc_onset, "inc_ref": A.on_ind[i, t] - 1500, "dec": r.dec_onset, "dec_ref": B.on_ind[i, t] - 1500})
cmp_ = pd.DataFrame(rows)
print(f"\nUnit-directions compared: {len(cmp_)}")
print(f"SDF correlation with the stored SDF: median {cmp_.r.median():.4f}; r > 0.9 for {(cmp_.r > 0.9).mean():.1%}; "
      f"r > 0.99 for {(cmp_.r > 0.99).mean():.1%}; mean rate ratio median {cmp_.mean_ratio.median():.4f}")
print(f"Valid trial counts equal to the paper's: {(cmp_.n == cmp_.n_ref).mean():.1%}; "
      f"median absolute difference {np.median(np.abs(cmp_.n - cmp_.n_ref)):.0f}, max {np.max(np.abs(cmp_.n - cmp_.n_ref)):.0f}")
for kind in ("inc", "dec"):
    a, b = cmp_[kind].notna(), cmp_[kind + "_ref"].notna()
    both = a & b
    print(f"{kind}reases: detected here {a.sum()}, by the paper {b.sum()}, agreement on presence {(a == b).mean():.1%}; "
          f"onset difference where both, median absolute {np.median(np.abs(cmp_[kind][both] - cmp_[kind + '_ref'][both])):.0f} ms")

Records kept by the paper's code, by animal and area:
area  GPi  VLa
monk          
Gaia  105   66
Isis  104  119

Found in every file they span:
           sum  size   mean
monk area                  
Gaia GPi    95   105  0.905
     VLa    60    66  0.909
Isis GPi   104   104  1.000
     VLa   105   119  0.882
Overall: 364 of 394 (92.4%); records spanning 2 to 5 files: 97

Records not found, by file:
fname
G_140524_3                          1
G_140524_8                          1
G_140525_1                          1
G_140731_1/G_140731_2               1
G_140731_2                          2
G_140813_1                          1
G_140813_2                          3
G_140813_3                          2
G_140815_1                          2
G_140815_1/G_140815_2               1
G_140815_3                          1
I_160620_3/I_160620_4               1
I_160621_1                          1
I_160621_3                          1
I_160622_1                          1
I_160622_3                          2
I_160624_1/I_160624_2               1
I_160629_4                          1
I_160630_2/I_160630_3               1
I_160701_1                          1
I_160707_1/I_160707_2/I_160707_3    1
I_160817_1                          2
I_160818_2                          1

Dandiset labels of the matched records, by the paper's area:
labels  GPi  Thalamus  VLa
area                      
GPi     199         0    0
VLa       0        25  140

Unit-directions compared: 724
SDF correlation with the stored SDF: median 0.9992; r > 0.9 for 94.6%; r > 0.99 for 89.5%; mean rate ratio median 1.0003
Valid trial counts equal to the paper's: 62.6%; median absolute difference 0, max 64
increases: detected here 538, by the paper 524, agreement on presence 90.9%; onset difference where both, median absolute 7 ms
decreases: detected here 330, by the paper 286, agreement on presence 80.7%; onset difference where both, median absolute 11 ms

counts = units.groupby(["area", "subject"]).size().unstack()
counts["total"] = counts.sum(axis=1)
print("Dandiset rule (one unit per file):")
print(counts.to_string())
for a, ref in (("GPi", [209]), ("VLa", [185, 218])):
    n = int(counts.loc[a, "total"])
    print(f"{a}: {n}; paper {ref}; within 10%: {[abs(n - x) <= 0.1 * x for x in ref]}")
day = units.session.str.extract(r"-(\d{6})-")[0]
distinct = units.assign(day=day).dropna(subset=["channel"]).drop_duplicates(["subject", "day", "area", "channel", "sort"])
print("\nVariant, units with the same channel and sort letter on the same day taken as one neuron (an assumption):")
print(distinct.groupby(["area", "subject"]).size().unstack().assign(total=lambda d: d.sum(axis=1)).to_string())
print("\nVariant, the paper's records found in the dandiset:")
print(paper[paper.keep & paper.complete].groupby(["area", "monk"]).size().unstack().assign(total=lambda d: d.sum(axis=1)).to_string())
in_paper_file = units.key.isin({f for x in paper.fname for f in x.split("/")})
print(f"\nDandiset-rule units in files the paper used: {int(in_paper_file.sum())} of {len(units)}")

Dandiset rule (one unit per file):
subject  Gaia  Isis  total
area                      
GPi       137   131    268
VLa        78   201    279
GPi: 268; paper [209]; within 10%: [False]
VLa: 279; paper [185, 218]; within 10%: [False, False]

Variant, units with the same channel and sort letter on the same day taken as one neuron (an assumption):
subject  Gaia  Isis  total
area                      
GPi        46    86    132
VLa        42   144    186

Variant, the paper's records found in the dandiset:
monk  Gaia  Isis  total
area                   
GPi     95   104    199
VLa     60   105    165

Dandiset-rule units in files the paper used: 518 of 547

def rest_table(df):
    t = df.groupby(["area", "subject"]).rest_rate.agg(["mean", "std", "size"])
    t2 = df.groupby("area").rest_rate.agg(["mean", "std", "size"]).assign(subject="all").set_index("subject", append=True)
    return pd.concat([t, t2]).sort_index()

print("Rest rate (spikes/s), mean SDF from -1500 to -301 ms, over unit-directions. Dandiset rule:")
print(rest_table(ud).round(2).to_string())
gp, vl = ud[ud.area == "GPi"].rest_rate, ud[ud.area == "VLa"].rest_rate
res = stats.ranksums(gp, vl)
print(f"GPi {gp.mean():.2f} against VLa {vl.mean():.2f}; rank sum z = {res.statistic:.2f}, p = {res.pvalue:.2g}")
for a in ("GPi", "VLa"):
    x = ud[ud.area == a]
    p = stats.ranksums(x[x.subject == "Gaia"].rest_rate, x[x.subject == "Isis"].rest_rate).pvalue
    print(f"{a}: Gaia against Isis, rank sum p = {p:.3g}")
print("\nRange of rest rates (spikes/s):", ud.groupby("area").rest_rate.agg(["min", "max"]).round(1).to_dict("index"))
print("\nVariant, the paper's records:")
print(rest_table(ud_paper).round(2).to_string())
print(f"rank sum p = {stats.ranksums(ud_paper[ud_paper.area == 'GPi'].rest_rate, ud_paper[ud_paper.area == 'VLa'].rest_rate).pvalue:.2g}")
print("\nVariant, one record per presumed neuron:")
print(rest_table(ud_m).round(2).to_string())
print(f"rank sum p = {stats.ranksums(ud_m[ud_m.area == 'GPi'].rest_rate, ud_m[ud_m.area == 'VLa'].rest_rate).pvalue:.2g}; "
      f"Gaia against Isis in GPi p = {stats.ranksums(ud_m[(ud_m.area == 'GPi') & (ud_m.subject == 'Gaia')].rest_rate, ud_m[(ud_m.area == 'GPi') & (ud_m.subject == 'Isis')].rest_rate).pvalue:.3g}")

Rest rate (spikes/s), mean SDF from -1500 to -301 ms, over unit-directions. Dandiset rule:
               mean    std  size
area subject                    
GPi  Gaia     73.17  21.55   272
     Isis     78.25  29.12   259
     all      75.65  25.63   531
VLa  Gaia     12.29  11.85   154
     Isis     12.97  11.54   398
     all      12.78  11.62   552
GPi 75.65 against VLa 12.78; rank sum z = 27.93, p = 1.2e-171
GPi: Gaia against Isis, rank sum p = 0.0409
VLa: Gaia against Isis, rank sum p = 0.531

Range of rest rates (spikes/s): {'GPi': {'min': 13.2, 'max': 152.3}, 'VLa': {'min': 0.3, 'max': 80.1}}

Variant, the paper's records:
               mean    std  size
area subject                    
GPi  Gaia     69.12  25.02   190
     Isis     59.75  33.70   205
     all      64.26  30.17   395
VLa  Gaia     15.54  11.12   120
     Isis     12.38   9.89   209
     all      13.53  10.45   329
rank sum p = 4.3e-94

Variant, one record per presumed neuron:
               mean    std  size
area subject                    
GPi  Gaia     74.23  21.63    92
     Isis     74.21  28.25   171
     all      74.22  26.08   263
VLa  Gaia     11.22   8.25    83
     Isis     13.22  11.49   286
     all      12.77  10.86   369
rank sum p = 3.9e-98; Gaia against Isis in GPi p = 0.691

from concurrent.futures import ThreadPoolExecutor
started = time.time()
jobs = {}
for r in dand_records:
    if r["uid"] not in kept or np.isnan(r["channel"]):
        continue
    s = r["sess"]
    rows_ = np.flatnonzero(pd.to_numeric(s["electrodes"].channel_ids, errors="coerce").to_numpy() == r["channel"])
    for name, p in s["processed"].items():
        cols = np.flatnonzero(np.isin(p["electrodes"], rows_))
        if len(cols) == 1:
            jobs.setdefault((s["url"], name), (p, []))[1].append((r["uid"], int(cols[0]), r["segments"][0][0]))
            break

def width_job(item):
    (url, name), (p, us) = item
    return g.ap_widths(url, name, p["rate"], p["start"], p["chunks"][0], us)

with ThreadPoolExecutor(16) as ex:
    widths = pd.DataFrame([w for ws in ex.map(width_job, jobs.items()) for w in ws])
widths = widths.merge(units[["uid", "area", "subject"]], on="uid")
print(f"Units with a channel in a filtered series: {len(widths)} of {len(units)}, in {len(jobs)} file-series; widths measured: "
      f"{int(widths.width_ms.notna().sum())} ({time.time() - started:.0f} s)")
print("Spikes per unit used: median", int(widths.n_spikes.median()), ", min", int(widths.n_spikes.min()))
print(widths.groupby(["area", "subject"]).width_ms.agg(["mean", "std", "median", "size"]).round(3).to_string())
print(widths.groupby("area").width_ms.agg(["mean", "std", "median", "min", "max", "size"]).round(3).to_string())
w_gp, w_vl = widths[widths.area == "GPi"].width_ms.dropna(), widths[widths.area == "VLa"].width_ms.dropna()
res = stats.ranksums(w_gp, w_vl)
print(f"GPi {w_gp.mean():.3f} ms against VLa {w_vl.mean():.3f} ms; rank sum z = {res.statistic:.2f}, p = {res.pvalue:.2g}")
print("Sampling interval of the trace: %.4f ms" % (1000 / list(sessions[0]["processed"].values())[0]["rate"]))
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist([w_gp, w_vl], bins=np.arange(0, 1.2, 0.041), label=["GPi", "VLa"])
ax.set_xlabel("trough to peak (ms)"); ax.set_ylabel("units"); ax.legend()
print("Units outside the axes (above 1.2 ms):", int((widths.width_ms > 1.2).sum()))
print(f"Widths of 0 (peak at the trough sample): {int((widths.width_ms == 0).sum())}; peak at the last sample of the 1.2 ms "
      f"window (width capped): {int(widths.at_end.sum())}; not measured (fewer than 10 spikes in the chunk): {int(widths.width_ms.isna().sum())}")
wz = widths[(widths.width_ms > 0) & ~widths.at_end]
print(f"Variant, without widths of 0 or capped: GPi {wz[wz.area == 'GPi'].width_ms.mean():.3f} ms (n {int((wz.area == 'GPi').sum())}), "
      f"VLa {wz[wz.area == 'VLa'].width_ms.mean():.3f} ms (n {int((wz.area == 'VLa').sum())}); rank sum p = "
      f"{stats.ranksums(wz[wz.area == 'GPi'].width_ms, wz[wz.area == 'VLa'].width_ms).pvalue:.2g}")

Units with a channel in a filtered series: 547 of 547, in 205 file-series; widths measured: 532 (263 s)
Spikes per unit used: median 273 , min 0
               mean    std  median  size
area subject                            
GPi  Gaia     0.239  0.064   0.246   137
     Isis     0.256  0.133   0.205   131
VLa  Gaia     0.489  0.223   0.532    78
     Isis     0.518  0.218   0.573   201
       mean    std  median    min    max  size
area                                          
GPi   0.247  0.104   0.246  0.000  0.983   268
VLa   0.509  0.219   0.553  0.082  1.065   279
GPi 0.247 ms against VLa 0.509 ms; rank sum z = -12.76, p = 2.9e-37
Sampling interval of the trace: 0.0410 ms
Units outside the axes (above 1.2 ms): 0
Widths of 0 (peak at the trough sample): 1; peak at the last sample of the 1.2 ms window (width capped): 9; not measured (fewer than 10 spikes in the chunk): 15
Variant, without widths of 0 or capped: GPi 0.249 ms (n 266), VLa 0.499 ms (n 257); rank sum p = 1.9e-35

def ramps(df):
    per_unit = df.groupby(["uid", "area"]).slope_sig.any().groupby("area").agg(["mean", "sum", "size"])
    return per_unit

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper)):
    t = ramps(df)
    print(name, "- units with a significant control-period slope in at least one direction:")
    print(t.assign(percent=lambda d: (100 * d["mean"]).round(1)).to_string())
    gs, vs = df[df.area == "GPi"].slope_per_s2, df[df.area == "VLa"].slope_per_s2
    print(f"  mean slope (spikes/s2) GPi {gs.mean():.3f}, VLa {vs.mean():.3f}; medians {gs.median():.3f}, {vs.median():.3f}; "
          f"rank sum p = {stats.ranksums(gs, vs).pvalue:.3g}")
    sig = df[df.slope_sig]
    print(f"  significant unit-direction slopes positive: {(sig.slope_per_s2 > 0).mean():.1%} of {len(sig)}")
print("By animal (dandiset rule):")
print(ud.groupby(["uid", "area", "subject"]).slope_sig.any().groupby(["area", "subject"]).mean().round(3).to_string())
print("Variant, significance at alpha 0.01 (99% CI) instead of 0.05:")
def sig01(df):
    out = []
    for _, r in df.iterrows():
        cue = int(round(1500 - r.median_rt_ms))
        x = np.arange(cue - 700, cue)
        lr = stats.linregress(x, r.sdf[x])
        out.append(lr.pvalue < 0.01)
    return pd.Series(out, index=df.index)
ud["slope_sig01"] = sig01(ud)
print(ud.groupby(["uid", "area"]).slope_sig01.any().groupby("area").mean().round(3).to_string())
print("Chance (random alignment, 3 draws), units with a significant slope:",
      nul.groupby(["rep", "uid", "area"]).slope_sig.any().groupby("area").mean().round(3).to_dict())
pn = ud.groupby(["neuron", "area"]).slope_per_s2.mean().reset_index()
print(f"One mean slope per presumed neuron (dandiset rule): GPi against VLa rank sum p = "
      f"{stats.ranksums(pn[pn.area == 'GPi'].slope_per_s2, pn[pn.area == 'VLa'].slope_per_s2).pvalue:.3g}")

Dandiset rule - units with a significant control-period slope in at least one direction:
          mean  sum  size  percent
area                              
GPi   0.985075  264   268     98.5
VLa   0.985663  275   279     98.6
  mean slope (spikes/s2) GPi -0.598, VLa 0.233; medians -0.531, 0.176; rank sum p = 0.0415
  significant unit-direction slopes positive: 49.9% of 969
Variant, one record per presumed neuron - units with a significant control-period slope in at least one direction:
          mean  sum  size  percent
area                              
GPi   0.992424  131   132     99.2
VLa   0.983871  183   186     98.4
  mean slope (spikes/s2) GPi -0.581, VLa 0.269; medians -0.739, 0.199; rank sum p = 0.0128
  significant unit-direction slopes positive: 50.5% of 562
Variant, the paper's records - units with a significant control-period slope in at least one direction:
          mean  sum  size  percent
area                              
GPi   0.994975  198   199     99.5
VLa   0.987879  163   165     98.8
  mean slope (spikes/s2) GPi -0.383, VLa -0.011; medians -0.353, 0.199; rank sum p = 0.187
  significant unit-direction slopes positive: 49.6% of 639
By animal (dandiset rule):
area  subject
GPi   Gaia       0.993
      Isis       0.977
VLa   Gaia       0.974
      Isis       0.990
Variant, significance at alpha 0.01 (99% CI) instead of 0.05:

area
GPi    0.970
VLa    0.975
Chance (random alignment, 3 draws), units with a significant slope: {'GPi': 0.984, 'VLa': 0.99}
One mean slope per presumed neuron (dandiset rule): GPi against VLa rank sum p = 0.0326

def responsive(df):
    return df.assign(resp=df.form != "none").groupby(["uid", "area", "subject"]).resp.any()

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper)):
    r = responsive(df)
    t = r.groupby("area").agg(["sum", "size", "mean"])
    print(name); print(t.assign(percent=lambda d: (100 * d["mean"]).round(1)).to_string())
    print(r.groupby(["area", "subject"]).mean().round(3).to_string())
print("\nUnit-directions with a response (dandiset rule):", (ud.form != "none").groupby(ud.area).agg(["sum", "size"]).to_dict("index"))
print("Variant, alpha 0.0001 instead of 0.001 (stricter):")
strict = []
for _, r in ud.iterrows():
    f_ = g.detect(r.sdf, r.median_rt_ms, alpha=0.0001)
    strict.append(np.isfinite(f_["onset"]))
print("  units with an increase in at least one direction:", ud.assign(s=strict).groupby(["uid", "area"]).s.any().groupby("area").mean().round(3).to_dict())
print("Chance (random alignment, 3 draws), responsive units:",
      nul.assign(resp=nul.form != "none").groupby(["rep", "uid", "area"]).resp.any().groupby("area").mean().round(3).to_dict())

Dandiset rule
      sum  size      mean  percent
area                              
GPi   264   268  0.985075     98.5
VLa   272   279  0.974910     97.5
area  subject
GPi   Gaia       0.971
      Isis       1.000
VLa   Gaia       0.974
      Isis       0.975
Variant, one record per presumed neuron
      sum  size      mean  percent
area                              
GPi   132   132  1.000000    100.0
VLa   183   186  0.983871     98.4
area  subject
GPi   Gaia       1.000
      Isis       1.000
VLa   Gaia       0.976
      Isis       0.986
Variant, the paper's records
      sum  size      mean  percent
area                              
GPi   199   199  1.000000    100.0
VLa   164   165  0.993939     99.4

area  subject
GPi   Gaia       1.000
      Isis       1.000
VLa   Gaia       0.983
      Isis       1.000

Unit-directions with a response (dandiset rule): {'GPi': {'sum': 497, 'size': 531}, 'VLa': {'sum': 497, 'size': 552}}
Variant, alpha 0.0001 instead of 0.001 (stricter):

  units with an increase in at least one direction: {'GPi': 0.81, 'VLa': 0.774}
Chance (random alignment, 3 draws), responsive units: {'GPi': 0.623, 'VLa': 0.708}

FORMS = ["increase", "decrease", "inc_first", "dec_first"]
def forms_table(df):
    r = df[df.form != "none"]
    t = pd.crosstab(r.area, r.form)[FORMS]
    pct = (100 * t.div(t.sum(axis=1), axis=0)).round(1)
    inc = t.increase + t.inc_first + t.dec_first
    changes = t.increase + t.decrease + 2 * (t.inc_first + t.dec_first)
    return t, pct, inc, changes

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper),
                 ("Variant, alpha 0.0001", ud_strict)):
    t, pct, inc, changes = forms_table(df)
    print(name, "- responses by form (unit-directions):"); print(t.to_string()); print("percent:"); print(pct.to_string())
    for a in t.index:
        print(f"  {a}: increases {inc[a]} of {changes[a]} changes ({100 * inc[a] / changes[a]:.1f}%)")
    chi = stats.chi2_contingency(np.array([[inc[a], changes[a] - inc[a]] for a in t.index]))
    allinc, allchg = inc.sum(), changes.sum()
    print(f"  increases against decreases over both areas: {allinc} of {allchg}, binomial p = {stats.binomtest(int(allinc), int(allchg)).pvalue:.2g}; "
          f"GPi against VLa chi-squared p = {chi.pvalue:.3g}")
print("By animal (dandiset rule), percent monophasic increases:")
r = ud[ud.form != "none"]
print((100 * pd.crosstab([r.area, r.subject], r.form, normalize="index")).round(1)[FORMS].to_string())

Dandiset rule - responses by form (unit-directions):
form  increase  decrease  inc_first  dec_first
area                                          
GPi        233       120         86         58
VLa        253       128         86         30
percent:
form  increase  decrease  inc_first  dec_first
area                                          
GPi       46.9      24.1       17.3       11.7
VLa       50.9      25.8       17.3        6.0
  GPi: increases 377 of 641 changes (58.8%)
  VLa: increases 369 of 613 changes (60.2%)
  increases against decreases over both areas: 746 of 1254, binomial p = 1.9e-11; GPi against VLa chi-squared p = 0.66
Variant, one record per presumed neuron - responses by form (unit-directions):
form  increase  decrease  inc_first  dec_first
area                                          
GPi        106        53         56         39
VLa        175        85         67         22
percent:
form  increase  decrease  inc_first  dec_first
area                                          
GPi       41.7      20.9       22.0       15.4
VLa       50.1      24.4       19.2        6.3
  GPi: increases 201 of 349 changes (57.6%)
  VLa: increases 264 of 438 changes (60.3%)
  increases against decreases over both areas: 465 of 787, binomial p = 3.9e-07; GPi against VLa chi-squared p = 0.492
Variant, the paper's records - responses by form (unit-directions):
form  increase  decrease  inc_first  dec_first
area                                          
GPi        192        69         69         48
VLa        170        85         43         16
percent:
form  increase  decrease  inc_first  dec_first
area                                          
GPi       50.8      18.3       18.3       12.7
VLa       54.1      27.1       13.7        5.1
  GPi: increases 309 of 495 changes (62.4%)
  VLa: increases 229 of 373 changes (61.4%)
  increases against decreases over both areas: 538 of 868, binomial p = 1.7e-12; GPi against VLa chi-squared p = 0.811
Variant, alpha 0.0001 - responses by form (unit-directions):
form  increase  decrease  inc_first  dec_first
area                                          
GPi        235       127         76         49
VLa        250       132         77         23
percent:
form  increase  decrease  inc_first  dec_first
area                                          
GPi       48.3      26.1       15.6       10.1
VLa       51.9      27.4       16.0        4.8
  GPi: increases 360 of 612 changes (58.8%)
  VLa: increases 350 of 582 changes (60.1%)
  increases against decreases over both areas: 710 of 1194, binomial p = 6.6e-11; GPi against VLa chi-squared p = 0.687
By animal (dandiset rule), percent monophasic increases:
form          increase  decrease  inc_first  dec_first
area subject                                          
GPi  Gaia         53.6      23.0       15.3        8.1
     Isis         40.2      25.3       19.3       15.3
VLa  Gaia         48.9      19.5       22.6        9.0
     Isis         51.6      28.0       15.4        4.9

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper)):
    r = df[df.form != "none"]
    mono = r.form.isin(["increase", "decrease"])
    t = pd.crosstab(r.area, mono)
    chi = stats.chi2_contingency(t.to_numpy())
    print(name); print((100 * mono.groupby(r.area).mean()).round(1).to_string())
    print(f"  counts (polyphasic, monophasic): {t.to_dict('index')}; chi-squared p = {chi.pvalue:.2g}")
r = ud[ud.form != "none"]
print("By animal (dandiset rule):", (100 * r.form.isin(["increase", "decrease"]).groupby([r.area, r.subject]).mean()).round(1).to_dict())

Dandiset rule
area
GPi    71.0
VLa    76.7
  counts (polyphasic, monophasic): {'GPi': {False: 144, True: 353}, 'VLa': {False: 116, True: 381}}; chi-squared p = 0.051
Variant, one record per presumed neuron
area
GPi    62.6
VLa    74.5
  counts (polyphasic, monophasic): {'GPi': {False: 95, True: 159}, 'VLa': {False: 89, True: 260}}; chi-squared p = 0.0023
Variant, the paper's records
area
GPi    69.0
VLa    81.2
  counts (polyphasic, monophasic): {'GPi': {False: 117, True: 261}, 'VLa': {False: 59, True: 255}}; chi-squared p = 0.00036
By animal (dandiset rule): {('GPi', 'Gaia'): 76.6, ('GPi', 'Isis'): 65.5, ('VLa', 'Gaia'): 68.4, ('VLa', 'Isis'): 79.7}

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper)):
    print(name)
    for a in ("GPi", "VLa"):
        x = df[df.area == a].integrated
        w = stats.wilcoxon(x)
        print(f"  {a}: mean {x.mean():.2f} spikes (SD {x.std():.2f}), median {x.median():.2f}, n {len(x)}; signed rank p = {w.pvalue:.2g}; "
              f"range {x.min():.1f} to {x.max():.1f}")
    print("  by animal:", df.groupby(["area", "subject"]).integrated.mean().round(2).to_dict())

Dandiset rule
  GPi: mean 3.07 spikes (SD 11.20), median 3.21, n 531; signed rank p = 1.1e-10; range -35.4 to 41.7
  VLa: mean 1.50 spikes (SD 4.34), median 0.78, n 552; signed rank p = 3.7e-18; range -10.9 to 38.2
  by animal: {('GPi', 'Gaia'): 4.22, ('GPi', 'Isis'): 1.87, ('VLa', 'Gaia'): 2.48, ('VLa', 'Isis'): 1.12}
Variant, one record per presumed neuron
  GPi: mean 2.13 spikes (SD 10.87), median 2.53, n 263; signed rank p = 0.00039; range -35.4 to 40.4
  VLa: mean 1.34 spikes (SD 3.68), median 0.92, n 369; signed rank p = 1.9e-12; range -9.3 to 20.1
  by animal: {('GPi', 'Gaia'): 1.59, ('GPi', 'Isis'): 2.42, ('VLa', 'Gaia'): 2.47, ('VLa', 'Isis'): 1.02}
Variant, the paper's records
  GPi: mean 3.42 spikes (SD 10.10), median 2.64, n 395; signed rank p = 4.9e-12; range -35.4 to 41.7
  VLa: mean 1.97 spikes (SD 5.45), median 1.22, n 329; signed rank p = 2.6e-11; range -9.3 to 38.2
  by animal: {('GPi', 'Gaia'): 4.41, ('GPi', 'Isis'): 2.5, ('VLa', 'Gaia'): 3.95, ('VLa', 'Isis'): 0.83}

def onsets(df, area, kind):
    col = {"inc": "inc_onset", "dec": "dec_onset"}
    if kind == "all":
        return pd.concat([df[df.area == area].inc_onset, df[df.area == area].dec_onset]).dropna()
    return df[df.area == area][col[kind]].dropna()

def compare(df, a1, k1, a2, k2, label):
    x, y = onsets(df, a1, k1), onsets(df, a2, k2)
    res = stats.ranksums(x, y)
    print(f"  {label}: {a1} {k1} median {x.median():.1f} ms (n {len(x)}), {a2} {k2} median {y.median():.1f} ms (n {len(y)}); "
          f"difference {x.median() - y.median():.1f} ms; rank sum z = {res.statistic:.2f}, p = {res.pvalue:.2g}")
    return x, y

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper),
                 ("Variant, alpha 0.0001", ud_strict)):
    print(name)
    x, y = compare(df, "GPi", "all", "VLa", "all", "all changes")
for subj in ("Gaia", "Isis"):
    print("Dandiset rule,", subj)
    compare(ud[ud.subject == subj], "GPi", "all", "VLa", "all", "all changes")

def compare_per_neuron(df, a1, k1, a2, k2, label):
    """One onset per presumed neuron: the median of its onsets of the kind over directions and files."""
    def per(a, k):
        d = df[df.area == a]
        cols = ["inc_onset", "dec_onset"] if k == "all" else [f"{k}_onset"]
        v = pd.concat([d[["neuron", c]].rename(columns={c: "o"}) for c in cols]).dropna()
        return v.groupby("neuron").o.median()
    x, y = per(a1, k1), per(a2, k2)
    print(f"  {label}, one value per presumed neuron: {a1} {k1} median {x.median():.1f} ms (n {len(x)}), "
          f"{a2} {k2} median {y.median():.1f} ms (n {len(y)}); rank sum p = {stats.ranksums(x, y).pvalue:.2g}")

print("Dandiset rule, variants")
compare_per_neuron(ud, "GPi", "all", "VLa", "all", "all changes")
x, y = onsets(ud, "GPi", "all"), onsets(ud, "VLa", "all")
print(f"  all changes, one-sided rank sum (GPi later): p = {stats.ranksums(x, y, alternative='greater').pvalue:.2g}")
x, y = onsets(ud, "GPi", "all"), onsets(ud, "VLa", "all")
fig, ax = plt.subplots(figsize=(5, 3))
for v, lab in ((x, "GPi"), (y, "VLa")):
    ax.step(np.sort(v), np.arange(1, len(v) + 1) / len(v), label=lab)
ax.set_xlabel("onset from movement onset (ms)"); ax.set_ylabel("cumulative fraction"); ax.legend()
print("Onset range (ms):", int(min(x.min(), y.min())), "to", int(max(x.max(), y.max())), "; all inside the axes")

Dandiset rule
  all changes: GPi all median -48.0 ms (n 641), VLa all median -67.0 ms (n 613); difference 19.0 ms; rank sum z = 1.75, p = 0.079
Variant, one record per presumed neuron
  all changes: GPi all median -73.0 ms (n 349), VLa all median -72.0 ms (n 438); difference -1.0 ms; rank sum z = 0.06, p = 0.95
Variant, the paper's records
  all changes: GPi all median -49.0 ms (n 495), VLa all median -85.0 ms (n 373); difference 36.0 ms; rank sum z = 2.88, p = 0.004
Variant, alpha 0.0001
  all changes: GPi all median -42.5 ms (n 612), VLa all median -63.0 ms (n 582); difference 20.5 ms; rank sum z = 1.94, p = 0.052
Dandiset rule, Gaia
  all changes: GPi all median -48.5 ms (n 306), VLa all median -48.0 ms (n 175); difference -0.5 ms; rank sum z = 0.43, p = 0.67
Dandiset rule, Isis
  all changes: GPi all median -45.0 ms (n 335), VLa all median -77.0 ms (n 438); difference 32.0 ms; rank sum z = 1.39, p = 0.16
Dandiset rule, variants
  all changes, one value per presumed neuron: GPi all median -47.0 ms (n 132), VLa all median -53.0 ms (n 184); rank sum p = 0.44
  all changes, one-sided rank sum (GPi later): p = 0.04
Onset range (ms): -272 to 439 ; all inside the axes

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper),
                 ("Variant, alpha 0.0001", ud_strict)):
    print(name)
    compare(df, "GPi", "dec", "VLa", "inc", "GPi decreases against VLa increases")
for subj in ("Gaia", "Isis"):
    print("Dandiset rule,", subj)
    compare(ud[ud.subject == subj], "GPi", "dec", "VLa", "inc", "GPi decreases against VLa increases")
print("Dandiset rule, variant")
compare_per_neuron(ud, "GPi", "dec", "VLa", "inc", "GPi decreases against VLa increases")

Dandiset rule
  GPi decreases against VLa increases: GPi dec median -29.5 ms (n 264), VLa inc median -84.0 ms (n 369); difference 54.5 ms; rank sum z = 3.70, p = 0.00021
Variant, one record per presumed neuron
  GPi decreases against VLa increases: GPi dec median -46.5 ms (n 148), VLa inc median -92.0 ms (n 264); difference 45.5 ms; rank sum z = 2.56, p = 0.011
Variant, the paper's records
  GPi decreases against VLa increases: GPi dec median -36.0 ms (n 186), VLa inc median -85.0 ms (n 229); difference 49.0 ms; rank sum z = 3.39, p = 0.00069
Variant, alpha 0.0001
  GPi decreases against VLa increases: GPi dec median -30.0 ms (n 252), VLa inc median -81.5 ms (n 350); difference 51.5 ms; rank sum z = 3.46, p = 0.00054
Dandiset rule, Gaia
  GPi decreases against VLa increases: GPi dec median 0.0 ms (n 115), VLa inc median -95.0 ms (n 107); difference 95.0 ms; rank sum z = 3.20, p = 0.0014
Dandiset rule, Isis
  GPi decreases against VLa increases: GPi dec median -44.0 ms (n 149), VLa inc median -81.0 ms (n 262); difference 37.0 ms; rank sum z = 1.92, p = 0.055
Dandiset rule, variant
  GPi decreases against VLa increases, one value per presumed neuron: GPi dec median -1.5 ms (n 102), VLa inc median -71.2 ms (n 154); rank sum p = 0.00074

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper),
                 ("Variant, alpha 0.0001", ud_strict)):
    print(name)
    compare(df, "GPi", "inc", "VLa", "dec", "GPi increases against VLa decreases")
for subj in ("Gaia", "Isis"):
    print("Dandiset rule,", subj)
    compare(ud[ud.subject == subj], "GPi", "inc", "VLa", "dec", "GPi increases against VLa decreases")
print("Dandiset rule, variant")
compare_per_neuron(ud, "GPi", "inc", "VLa", "dec", "GPi increases against VLa decreases")

Dandiset rule
  GPi increases against VLa decreases: GPi inc median -69.0 ms (n 377), VLa dec median -19.0 ms (n 244); difference -50.0 ms; rank sum z = -1.38, p = 0.17
Variant, one record per presumed neuron
  GPi increases against VLa decreases: GPi inc median -92.0 ms (n 201), VLa dec median -8.0 ms (n 174); difference -84.0 ms; rank sum z = -2.64, p = 0.0082
Variant, the paper's records
  GPi increases against VLa decreases: GPi inc median -73.0 ms (n 309), VLa dec median -86.0 ms (n 144); difference 13.0 ms; rank sum z = 1.28, p = 0.2
Variant, alpha 0.0001
  GPi increases against VLa decreases: GPi inc median -57.5 ms (n 360), VLa dec median -21.5 ms (n 232); difference -36.0 ms; rank sum z = -1.03, p = 0.3
Dandiset rule, Gaia
  GPi increases against VLa decreases: GPi inc median -75.0 ms (n 191), VLa dec median 12.5 ms (n 68); difference -87.5 ms; rank sum z = -2.99, p = 0.0028
Dandiset rule, Isis
  GPi increases against VLa decreases: GPi inc median -46.0 ms (n 186), VLa dec median -68.0 ms (n 176); difference 22.0 ms; rank sum z = -0.08, p = 0.94
Dandiset rule, variant
  GPi increases against VLa decreases, one value per presumed neuron: GPi inc median -45.0 ms (n 115), VLa dec median 20.5 ms (n 124); rank sum p = 0.021

def gating(df, lo=0, hi=50):
    def changes(a):
        d = df[df.area == a]
        return np.r_[d.inc_onset.dropna(), d.dec_onset.dropna()], np.r_[np.ones(d.inc_onset.notna().sum()), -np.ones(d.dec_onset.notna().sum())]
    gt, gs = changes("GPi")
    vt, vs = changes("VLa")
    lead = vt[None, :] - gt[:, None]
    fits = (lead >= lo) & (lead <= hi) & (gs[:, None] != vs[None, :])
    return fits.size, int(fits.sum())

for name, df in (("Dandiset rule", ud), ("Variant, one record per presumed neuron", ud_m), ("Variant, the paper's records", ud_paper)):
    n, fit = gating(df)
    print(f"{name}: {n} comparisons, {fit} fit gating; conflicting {100 * (1 - fit / n):.1f}%")
for lo, hi in ((0, 20), (0, 100), (0, 200)):
    n, fit = gating(ud, lo, hi)
    print(f"  variant, GPi leading by {lo} to {hi} ms: conflicting {100 * (1 - fit / n):.1f}%")
n, fit = gating(ud)
rng = np.random.default_rng(1)
d = ud.copy()
null = []
for _ in range(200):
    d["inc_onset"] = rng.permutation(ud.inc_onset.to_numpy())
    d["dec_onset"] = rng.permutation(ud.dec_onset.to_numpy())
    nn, ff = gating(d)
    null.append(100 * (1 - ff / nn))
print(f"Reference: onsets shuffled across unit-directions, conflicting {np.mean(null):.1f}% (2.5 to 97.5 percentile {np.percentile(null, 2.5):.1f} to {np.percentile(null, 97.5):.1f})")

Dandiset rule: 392933 comparisons, 17921 fit gating; conflicting 95.4%
Variant, one record per presumed neuron: 152862 comparisons, 7493 fit gating; conflicting 95.1%
Variant, the paper's records: 184635 comparisons, 8098 fit gating; conflicting 95.6%

  variant, GPi leading by 0 to 20 ms: conflicting 98.0%
  variant, GPi leading by 0 to 100 ms: conflicting 91.9%
  variant, GPi leading by 0 to 200 ms: conflicting 86.6%

Reference: onsets shuffled across unit-directions, conflicting 95.1% (2.5 to 97.5 percentile 95.0 to 95.2)

def directional(df, records):
    seg = {r["uid"]: r["segments"] for r in records}
    out = []
    for uid, d in df[df.form != "none"].groupby("uid"):
        first = np.nanmin(np.r_[d.inc_onset, d.dec_onset]) / 1000.0
        c = {}
        for target in ("Left", "Right"):
            vals = []
            for spikes, t in seg[uid]:
                tt = t[t.valid & (t.target == target)]
                vals.append(window_counts(spikes, tt.movement_start_time.to_numpy(), first, first + 0.3))
            c[target] = np.concatenate(vals)
        p = stats.ranksums(c["Left"], c["Right"]).pvalue if min(len(c["Left"]), len(c["Right"])) else np.nan
        out.append({"uid": uid, "area": d.area.iloc[0], "subject": d.subject.iloc[0], "p": p})
    return pd.DataFrame(out)

for name, df, recs in (("Dandiset rule", ud, dand_records), ("Variant, one record per presumed neuron", ud_m, merged_records),
                       ("Variant, the paper's records", ud_paper, paper_records)):
    dd = directional(df, recs)
    t = dd.assign(dir=dd.p < 0.005).groupby("area").dir.agg(["sum", "size", "mean"])
    print(name); print(t.assign(percent=lambda x: (100 * x["mean"]).round(1)).to_string())
    chi = stats.chi2_contingency(pd.crosstab(dd.area, dd.p < 0.005).to_numpy())
    print(f"  GPi against VLa chi-squared p = {chi.pvalue:.2g}")
    if name == "Dandiset rule":
        print("  by animal:", dd.assign(dir=dd.p < 0.005).groupby(["area", "subject"]).dir.mean().round(3).to_dict())
        print("  variant p < 0.05:", dd.assign(dir=dd.p < 0.05).groupby("area").dir.mean().round(3).to_dict())

Dandiset rule
      sum  size      mean  percent
area                              
GPi   122   264  0.462121     46.2
VLa    66   272  0.242647     24.3
  GPi against VLa chi-squared p = 1.7e-07
  by animal: {('GPi', 'Gaia'): 0.391, ('GPi', 'Isis'): 0.534, ('VLa', 'Gaia'): 0.237, ('VLa', 'Isis'): 0.245}
  variant p < 0.05: {'GPi': 0.564, 'VLa': 0.419}

Variant, one record per presumed neuron
      sum  size      mean  percent
area                              
GPi    71   132  0.537879     53.8
VLa    40   183  0.218579     21.9
  GPi against VLa chi-squared p = 9.8e-09

Variant, the paper's records
      sum  size      mean  percent
area                              
GPi    96   199  0.482412     48.2
VLa    44   164  0.268293     26.8
  GPi against VLa chi-squared p = 4.9e-05

pairs = []
for s in sessions:
    us = units[units.session == s["session_id"]]
    gps, vls = us[us.area == "GPi"], us[us.area == "VLa"]
    for _, a in gps.iterrows():
        for _, b in vls.iterrows():
            pairs.append({"session": s["session_id"], "subject": s["subject"], "gpi": a.uid, "vla": b.uid})
pairs = pd.DataFrame(pairs)
spk = {r["uid"]: r["segments"][0][0] for r in dand_records}
sess_of = {s["session_id"]: s for s in sessions}
print(f"GPi-VLa pairs in the same file: {len(pairs)} in {pairs.session.nunique()} files")
print(pairs.groupby("subject").size().to_string())

# Variant, the paper's records: a GPi and a VLa record that share a file make a pair in that file.
analysed = set(ud_paper.uid)
pairs_paper, pair_neuron = [], dict(uid2neuron)
for s in sessions:
    sid = s["session_id"]
    if s["vt"] is None:
        continue
    here = {a: [(r, r["sessions"].index(sid)) for r in paper_records
                if r["uid"] in analysed and r["area"] == a and sid in r["sessions"]] for a in ("GPi", "VLa")}
    for ra, ka in here["GPi"]:
        for rb, kb in here["VLa"]:
            ga, vb = f"{ra['uid']}@{sid}", f"{rb['uid']}@{sid}"
            spk[ga], spk[vb] = ra["segments"][ka][0], rb["segments"][kb][0]
            pair_neuron[ga], pair_neuron[vb] = ra["uid"], rb["uid"]
            pairs_paper.append({"session": sid, "subject": s["subject"], "gpi": ga, "vla": vb})
pairs_paper = pd.DataFrame(pairs_paper)
n_rec_pairs = (pairs_paper.gpi.str.split("@").str[0] + "|" + pairs_paper.vla.str.split("@").str[0]).nunique()
print(f"Variant, the paper's records: GPi-VLa pairs in the same file {len(pairs_paper)} in {pairs_paper.session.nunique()} files "
      f"({n_rec_pairs} distinct pairs of records)")
print(pairs_paper.groupby("subject").size().to_string())

GPi-VLa pairs in the same file: 338 in 76 files
subject
Gaia     88
Isis    250
Variant, the paper's records: GPi-VLa pairs in the same file 329 in 73 files (313 distinct pairs of records)
subject
Gaia    104
Isis    225

started = time.time()
N_SURR = 400

def ccf_run(pairs_df, rng, rng_c0):
    ccf_rows, ctrl_peaks, same_time = [], {"whole": [], "rest": [], "movement": []}, []
    for _, p in pairs_df.iterrows():
        s = sess_of[p.session]
        a, b = spk[p.gpi], spk[p.vla]
        span = max(a[-1], b[-1]) - min(a[0], b[0])
        same_time.append({"same": len(np.intersect1d(np.round(a * 24414.0625), np.round(b * 24414.0625))),
                          "expected": len(a) * len(b) / (span * 24414.0625)})
        t = s["vt"][s["vt"].valid]
        periods = {"whole": None,
                   "rest": np.c_[t.movement_start_time + g.REST_WIN[0], t.movement_start_time + g.REST_WIN[1]],
                   "movement": np.c_[t.movement_start_time + g.MVT_WIN[0], t.movement_start_time + g.MVT_WIN[1]]}
        for name, win in periods.items():
            if win is None:
                x, y = a, b
                dur = span
            else:
                x, y = a[g._in_windows(a, win)], b[g._in_windows(b, win)]
                dur = float(np.sum(win[:, 1] - win[:, 0]))
            if min(len(x), len(y)) < 100:
                continue
            r = g.ccf_test(x, y, dur, rng, N_SURR)
            ccf_rows.append({"period": name, "subject": p.subject, "neuron_pair": pair_neuron[p.gpi] + "|" + pair_neuron[p.vla],
                             "peak": r["peak"], "sign": r["sign"], "lag_ms": r["lag_ms"], "thr95": r["thr95"],
                             "significant": r["significant"], "c0": g.min_detectable(r["expected"], r["thr95"], rng_c0)})
            ctrl_peaks[name].append(r["ctrl_peaks"])
    return pd.DataFrame(ccf_rows), ctrl_peaks, pd.DataFrame(same_time)

cc, ctrl_peaks, st = ccf_run(pairs, np.random.default_rng(4), np.random.default_rng(14))
print(f"CCFs computed in {time.time() - started:.0f} s; {N_SURR} jittered controls per pair and period")
for name in ("whole", "rest", "movement"):
    d = cc[cc.period == name]
    ctrl = np.vstack(ctrl_peaks[name])
    k, n = int(d.significant.sum()), len(d)
    ks_p = g.ks_permutation(d.peak.to_numpy(), ctrl)
    print(f"{name}: {n} pairs; significant {k} ({100 * k / n:.1f}%), binomial p against 5% (greater) = "
          f"{stats.binomtest(k, n, 0.05, alternative='greater').pvalue:.3g}; signs of significant peaks "
          f"+{int((d[d.significant].sign > 0).sum())} / -{int((d[d.significant].sign < 0).sum())}; "
          f"KS permutation p = {ks_p:.3f}; median peak {d.peak.median():.3f}, median control peak {np.median(ctrl):.3f}, "
          f"median control 95th percentile {d.thr95.median():.3f}")
    print(f"   lags of the significant peaks (ms): {d[d.significant].lag_ms.value_counts().sort_index().to_dict()}")
    keep_c0 = (d.c0 <= 0.3).to_numpy()
    k0 = int(d[keep_c0].significant.sum())
    print(f"   variant, pairs with c0 <= 0.3 (approximate c0, median {d.c0.median():.2f}): {keep_c0.sum()} pairs, significant {k0} "
          f"({100 * k0 / keep_c0.sum():.1f}%), KS permutation p = {g.ks_permutation(d[keep_c0].peak.to_numpy(), ctrl[keep_c0]):.3f}")
    first = (~d.neuron_pair.duplicated()).to_numpy()
    print(f"   variant, one file per pair of presumed neurons: {first.sum()} pairs, KS permutation p = "
          f"{g.ks_permutation(d[first].peak.to_numpy(), ctrl[first]):.3f}")
print("By animal, percent significant:", (100 * cc.groupby(["period", "subject"]).significant.mean()).round(1).to_dict())
print(f"GPi and VLa spikes at the same sample, summed over {len(st)} pairs: {int(st.same.sum())}, expected for independent trains "
      f"{st.expected.sum():.0f}; largest ratio in a pair {(st.same / st.expected).max():.2f}")

started = time.time()
cc_p, ctrl_p, _ = ccf_run(pairs_paper, np.random.default_rng(24), np.random.default_rng(34))
print(f"\nVariant, the paper's records ({time.time() - started:.0f} s):")
for name in ("whole", "rest", "movement"):
    d = cc_p[cc_p.period == name]
    ctrl = np.vstack(ctrl_p[name])
    k, n = int(d.significant.sum()), len(d)
    keep_c0 = (d.c0 <= 0.3).to_numpy()
    print(f"  {name}: {n} pairs; significant {k} ({100 * k / n:.1f}%), signs +{int((d[d.significant].sign > 0).sum())} / "
          f"-{int((d[d.significant].sign < 0).sum())}; KS permutation p = {g.ks_permutation(d.peak.to_numpy(), ctrl):.3f}; "
          f"with c0 <= 0.3: {keep_c0.sum()} pairs, KS permutation p = {g.ks_permutation(d[keep_c0].peak.to_numpy(), ctrl[keep_c0]):.3f}")

CCFs computed in 63 s; 400 jittered controls per pair and period

whole: 338 pairs; significant 19 (5.6%), binomial p against 5% (greater) = 0.333; signs of significant peaks +15 / -4; KS permutation p = 0.003; median peak 0.042, median control peak 0.038, median control 95th percentile 0.059
   lags of the significant peaks (ms): {0: 1, 1: 4, 2: 1, 3: 3, 4: 1, 5: 1, 7: 3, 8: 1, 9: 2, 10: 2}

   variant, pairs with c0 <= 0.3 (approximate c0, median 0.09): 332 pairs, significant 19 (5.7%), KS permutation p = 0.003

   variant, one file per pair of presumed neurons: 282 pairs, KS permutation p = 0.003

rest: 298 pairs; significant 10 (3.4%), binomial p against 5% (greater) = 0.932; signs of significant peaks +8 / -2; KS permutation p = 0.705; median peak 0.134, median control peak 0.134, median control 95th percentile 0.215
   lags of the significant peaks (ms): {0: 2, 1: 2, 2: 1, 4: 1, 5: 1, 9: 3}

   variant, pairs with c0 <= 0.3 (approximate c0, median 0.34): 121 pairs, significant 4 (3.3%), KS permutation p = 0.333

   variant, one file per pair of presumed neurons: 257 pairs, KS permutation p = 0.627

movement: 306 pairs; significant 20 (6.5%), binomial p against 5% (greater) = 0.136; signs of significant peaks +15 / -5; KS permutation p = 0.045; median peak 0.116, median control peak 0.113, median control 95th percentile 0.179
   lags of the significant peaks (ms): {0: 2, 1: 1, 2: 1, 3: 2, 4: 1, 5: 4, 6: 2, 7: 2, 8: 1, 9: 1, 10: 3}

   variant, pairs with c0 <= 0.3 (approximate c0, median 0.29): 172 pairs, significant 12 (7.0%), KS permutation p = 0.025

   variant, one file per pair of presumed neurons: 262 pairs, KS permutation p = 0.037
By animal, percent significant: {('movement', 'Gaia'): 6.8, ('movement', 'Isis'): 6.5, ('rest', 'Gaia'): 1.4, ('rest', 'Isis'): 4.0, ('whole', 'Gaia'): 2.3, ('whole', 'Isis'): 6.8}
GPi and VLa spikes at the same sample, summed over 338 pairs: 11376, expected for independent trains 11435; largest ratio in a pair 3.02

def detrend(c, x):
    return c - np.polyval(np.polyfit(x, c, 1), x)

def nc_run(pairs_df, rng, rng_d):
    nc_rows, nc_ctrl, nc_ctrl_d = [], {"rest": [], "movement": []}, {"rest": [], "movement": []}
    for _, p in pairs_df.iterrows():
        s = sess_of[p.session]
        for target in ("Left", "Right"):
            t = s["vt"][s["vt"].valid & (s["vt"].target == target)]
            al = t.movement_start_time.to_numpy()
            for name, (lo, hi) in (("rest", g.REST_WIN), ("movement", g.MVT_WIN)):
                ca = window_counts(spk[p.gpi], al, lo, hi).astype(float)
                cb = window_counts(spk[p.vla], al, lo, hi).astype(float)
                if ca.sum() < 10 or cb.sum() < 10 or len(al) < 10 or ca.std() == 0 or cb.std() == 0:
                    continue
                za, zb = stats.zscore(ca), stats.zscore(cb)
                ok = (np.abs(za) <= 3) & (np.abs(zb) <= 3)
                ca, cb, x = ca[ok], cb[ok], al[ok]
                if ca.std() == 0 or cb.std() == 0:
                    continue
                r = np.corrcoef(ca, cb)[0, 1]
                sh = np.array([np.corrcoef(ca, rng.permutation(cb))[0, 1] for _ in range(400)])
                da, db = detrend(ca, x), detrend(cb, x)
                shd = np.array([np.corrcoef(da, rng_d.permutation(db))[0, 1] for _ in range(400)])
                own = np.percentile(sh, [2.5, 97.5])
                nc_rows.append({"period": name, "pair": f"{p.gpi}|{p.vla}", "neuron_pair": pair_neuron[p.gpi] + "|" + pair_neuron[p.vla],
                                "subject": p.subject, "target": target, "r": r, "r_detrended": np.corrcoef(da, db)[0, 1],
                                "drift_gpi": stats.spearmanr(x, ca).pvalue < 0.05, "drift_vla": stats.spearmanr(x, cb).pvalue < 0.05,
                                "sig_own": bool(r < own[0] or r > own[1])})
                nc_ctrl[name].append(sh)
                nc_ctrl_d[name].append(shd)
    return pd.DataFrame(nc_rows), nc_ctrl, nc_ctrl_d

nc, nc_ctrl, nc_ctrl_d = nc_run(pairs, np.random.default_rng(5), np.random.default_rng(15))
for name in ("rest", "movement"):
    d = nc[nc.period == name].copy()
    ctrl, ctrl_d = np.vstack(nc_ctrl[name]), np.vstack(nc_ctrl_d[name])
    lo, hi = np.percentile(ctrl, [2.5, 97.5])
    d["sig"] = (d.r < lo) | (d.r > hi)
    pos, neg = int((d.r > hi).sum()), int((d.r < lo).sum())
    ks_p = g.ks_permutation(np.abs(d.r.to_numpy()), np.abs(ctrl))
    by_pair = d.groupby("pair").sig.any()
    print(f"{name}: (a) {len(d)} pair-targets, significant {int(d.sig.sum())} ({100 * d.sig.mean():.1f}%), +{pos} / -{neg} "
          f"(binomial p of the sign balance {stats.binomtest(pos, pos + neg).pvalue:.3g}); "
          f"(b) {len(by_pair)} pairs, significant in either target {int(by_pair.sum())} ({100 * by_pair.mean():.1f}%)")
    print(f"   KS permutation p (absolute correlations against shuffles) = {ks_p:.3f}; median r {d.r.median():.3f}; "
          f"shuffle band {lo:.3f} to {hi:.3f}")
    print(f"   variant, each pair-target against the 2.5 to 97.5 percentiles of its own shuffles: {100 * d.sig_own.mean():.1f}% significant")
    print(f"   drift: spike counts correlated with trial time (Spearman p < 0.05) for GPi in {100 * d.drift_gpi.mean():.1f}%, "
          f"VLa in {100 * d.drift_vla.mean():.1f}%, both in {100 * (d.drift_gpi & d.drift_vla).mean():.1f}% of pair-targets")
    lo_d, hi_d = np.percentile(ctrl_d, [2.5, 97.5])
    print(f"   variant, counts linearly detrended over trial time: median |r| {d.r_detrended.abs().median():.3f} (raw {d.r.abs().median():.3f}); "
          f"significant {100 * ((d.r_detrended < lo_d) | (d.r_detrended > hi_d)).mean():.1f}%; KS permutation p = "
          f"{g.ks_permutation(np.abs(d.r_detrended.to_numpy()), np.abs(ctrl_d)):.3f}")
    first = d.pair.isin(d.drop_duplicates("neuron_pair").pair)
    print(f"   variant, one file per pair of presumed neurons: {int(first.sum())} pair-targets, KS permutation p = "
          f"{g.ks_permutation(np.abs(d[first].r.to_numpy()), np.abs(ctrl[first.to_numpy()])):.3f}")
nc_p, nc_ctrl_p, _ = nc_run(pairs_paper, np.random.default_rng(25), np.random.default_rng(35))
print("Variant, the paper's records:")
for name in ("rest", "movement"):
    d = nc_p[nc_p.period == name]
    ctrl = np.vstack(nc_ctrl_p[name])
    lo, hi = np.percentile(ctrl, [2.5, 97.5])
    sig = (d.r < lo) | (d.r > hi)
    print(f"  {name}: {len(d)} pair-targets, significant {100 * sig.mean():.1f}%, +{int((d.r > hi).sum())} / -{int((d.r < lo).sum())}; "
          f"KS permutation p = {g.ks_permutation(np.abs(d.r.to_numpy()), np.abs(ctrl)):.3f}")

rest: (a) 665 pair-targets, significant 67 (10.1%), +38 / -29 (binomial p of the sign balance 0.328); (b) 338 pairs, significant in either target 59 (17.5%)
   KS permutation p (absolute correlations against shuffles) = 0.000; median r 0.016; shuffle band -0.291 to 0.291
   variant, each pair-target against the 2.5 to 97.5 percentiles of its own shuffles: 10.8% significant
   drift: spike counts correlated with trial time (Spearman p < 0.05) for GPi in 18.9%, VLa in 36.1%, both in 7.8% of pair-targets

   variant, counts linearly detrended over trial time: median |r| 0.114 (raw 0.108); significant 9.3%; KS permutation p = 0.000

   variant, one file per pair of presumed neurons: 555 pair-targets, KS permutation p = 0.000

movement: (a) 669 pair-targets, significant 72 (10.8%), +39 / -33 (binomial p of the sign balance 0.556); (b) 338 pairs, significant in either target 63 (18.6%)
   KS permutation p (absolute correlations against shuffles) = 0.000; median r 0.001; shuffle band -0.291 to 0.292
   variant, each pair-target against the 2.5 to 97.5 percentiles of its own shuffles: 14.2% significant
   drift: spike counts correlated with trial time (Spearman p < 0.05) for GPi in 29.0%, VLa in 38.3%, both in 12.1% of pair-targets

   variant, counts linearly detrended over trial time: median |r| 0.111 (raw 0.121); significant 9.7%; KS permutation p = 0.000

   variant, one file per pair of presumed neurons: 559 pair-targets, KS permutation p = 0.000

Variant, the paper's records:

  rest: 647 pair-targets, significant 9.1%, +28 / -31; KS permutation p = 0.000

  movement: 648 pair-targets, significant 7.4%, +24 / -24; KS permutation p = 0.000

started = time.time()
LAGS = np.arange(-800, 801)
win = (LAGS >= 0) & (LAGS <= 300)

def fig5_run(pairs_df, rng):
    ta = {"burst onset": [], "burst offset": [], "pause onset": []}
    ctrl_sum = {k: 0.0 for k in ta}
    ta_pct = {k: [] for k in ta}
    ctrl_pct_sum = {k: 0.0 for k in ta}
    burst_cache, burst_stats = {}, []
    pair_sig = {k: [] for k in ta}
    pair_sig_ctrl = {k: [] for k in ta}
    n_excluded, n_gap = {k: 0 for k in ta}, {k: 0 for k in ta}
    n_events = {k: [] for k in ta}
    BASE = (LAGS >= -800) & (LAGS < -50)
    POST = (LAGS >= 0) & (LAGS <= 100)
    for _, p in pairs_df.iterrows():
        s = sess_of[p.session]
        t = s["trials"]
        end = t["return_stop_time"] if "return_stop_time" in t else t["stop_time"]
        task = np.c_[t.movement_start_time + g.TASK_EXCL[0], end.fillna(t.stop_time) + g.TASK_EXCL[1]]
        task = task[np.isfinite(task).all(1)]
        a, b = spk[p.gpi], spk[p.vla]
        T = int(np.ceil(max(a[-1], b[-1]) * 1000)) + 2000
        in_task = np.zeros(T, bool)
        for lo_, hi_ in (task * 1000).astype(int):
            in_task[max(lo_, 0):min(hi_, T)] = True
        if p.gpi not in burst_cache:
            rest_a = a[~in_task[(a * 1000).astype(int)]]
            rest_dur = (~in_task[int(a[0] * 1000):int(a[-1] * 1000)]).sum() / 1000
            rate = len(rest_a) / rest_dur
            edges = np.flatnonzero(np.diff(np.r_[0, (~in_task).astype(int), 0]))
            segs = [a[(a >= lo_) & (a < hi_)] for lo_, hi_ in edges.reshape(-1, 2) / 1000]      # rest segments
            bl = [x for seg in segs for x in g.bursts(seg, rate)]
            pl = [x for seg in segs for x in g.pauses(seg, rate)]
            bl10 = [x for seg in segs for x in g.bursts(seg, rate, log=np.log10)]
            pl10 = [x for seg in segs for x in g.pauses(seg, rate, log=np.log10)]
            burst_cache[p.gpi] = (bl, pl)
            burst_stats.append({"gpi": p.gpi, "rest_s": rest_dur, "rate": rate, "n_bursts": len(bl), "n_pauses": len(pl),
                                "burst_ms": np.mean([1000 * (e - b0) for b0, e in bl]) if bl else np.nan,
                                "pause_ms": np.mean([1000 * (e - b0) for b0, e in pl]) if pl else np.nan,
                                "burst_ms_log10": np.mean([1000 * (e - b0) for b0, e in bl10]) if bl10 else np.nan,
                                "pause_ms_log10": np.mean([1000 * (e - b0) for b0, e in pl10]) if pl10 else np.nan})
        bl, pl = burst_cache[p.gpi]
        events = {"burst onset": np.array([e[0] for e in bl]), "burst offset": np.array([e[1] for e in bl]),
                  "pause onset": np.array([e[0] for e in pl])}
        counts = np.bincount((b * 1000).astype(int), minlength=T)[:T].astype(float)
        rest = (~in_task).astype(float)
        vrate = gaussian_filter1d(counts * 1000.0, 10.0) * rest       # VLa rate, rest samples only
        for k, ev in events.items():
            ev = ev[(ev > 1.0) & (ev < T / 1000 - 1.0)]
            if len(ev) < 5:
                continue
            idx = (ev * 1000).astype(int)
            X, W = g.event_sums(vrate, idx), g.event_sums(rest, idx)
            ratio = np.where(W > 0.5, X / np.maximum(W, 0.5), np.nan)     # lags where no event has a rest sample: NaN
            M = len(ratio)
            avg = ratio[LAGS % M]
            if np.isnan(avg).any():
                n_gap[k] += 1
                continue
            base = avg[BASE]
            lo_b, hi_b = np.percentile(base, [2.5, 97.5])
            if hi_b - lo_b > 0.6 * abs(base.mean()):
                n_excluded[k] += 1
                continue
            ta[k].append(avg)
            ta_pct[k].append(100 * avg / base.mean())
            n_events[k].append(len(idx))
            post = avg[POST]
            pair_sig[k].append("decrease" if (post < lo_b).any() else ("increase" if (post > hi_b).any() else "none"))
            shifts = rng.integers(2000, 20001, size=400) * rng.choice([-1, 1], size=400)
            c_ = np.array([ratio[(LAGS + sh) % M] for sh in shifts])
            c_ = np.where(np.isnan(c_), np.nanmean(avg), c_)
            ctrl_sum[k] = ctrl_sum[k] + c_
            ctrl_pct_sum[k] = ctrl_pct_sum[k] + 100 * c_ / c_[:, BASE].mean(1, keepdims=True)
            # the per-pair test applied to each shifted control: its chance level
            clo, chi_ = np.percentile(c_[:, BASE], [2.5, 97.5], axis=1)
            cpost = c_[:, POST]
            pair_sig_ctrl[k].append(np.mean((cpost < clo[:, None]).any(1) | (cpost > chi_[:, None]).any(1)))
    return dict(ta=ta, ctrl_sum=ctrl_sum, ta_pct=ta_pct, ctrl_pct_sum=ctrl_pct_sum, burst_stats=pd.DataFrame(burst_stats),
                pair_sig=pair_sig, pair_sig_ctrl=pair_sig_ctrl, n_excluded=n_excluded, n_gap=n_gap, n_events=n_events)

def population_band(ta_k, ctrl_sum_k):
    """Population average and the control band, widened so that 5% of controls cross it in 0 to 300 ms."""
    pop = np.mean(ta_k, 0)
    c = ctrl_sum_k / len(ta_k)                   # control population averages, 400 x lags
    lo_c, hi_c = np.percentile(c, [2.5, 97.5], axis=0)
    shift = min(np.percentile((c[:, win] - lo_c[win]).min(1), 5), 0)
    return pop, lo_c + shift, hi_c - shift

F = fig5_run(pairs, np.random.default_rng(6))
ta, ctrl_sum, ta_pct, ctrl_pct_sum = F["ta"], F["ctrl_sum"], F["ta_pct"], F["ctrl_pct_sum"]
pair_sig, pair_sig_ctrl, n_excluded, n_gap, n_events = F["pair_sig"], F["pair_sig_ctrl"], F["n_excluded"], F["n_gap"], F["n_events"]
bs = F["burst_stats"]
print(f"Bursts and pauses in {time.time() - started:.0f} s")
print(f"GPi units: {len(bs)}; rest per unit median {bs.rest_s.median():.0f} s; bursts per minute of rest median "
      f"{(60 * bs.n_bursts / bs.rest_s).median():.2f}; pauses per minute median {(60 * bs.n_pauses / bs.rest_s).median():.2f}")
print(f"Mean burst duration {np.nanmean(bs.burst_ms):.1f} ms (paper 24 ms); mean pause duration {np.nanmean(bs.pause_ms):.1f} ms (paper 195 ms)")
print(f"Variant, surprise with log10: mean burst duration {np.nanmean(bs.burst_ms_log10):.1f} ms; mean pause duration {np.nanmean(bs.pause_ms_log10):.1f} ms")
fig, axes = plt.subplots(1, 3, figsize=(12, 3), sharey=True)
for ax, k in zip(axes, ta):
    pop, lo_corr, hi_corr = population_band(ta[k], ctrl_sum[k])
    below = win & (pop < lo_corr)
    above = win & (pop > hi_corr)
    first = LAGS[below][0] if below.any() else None
    sig = pd.Series(pair_sig[k]).value_counts()
    print(f"{k}: {len(ta[k])} pairs of {len(ta[k]) + n_excluded[k] + n_gap[k]} with at least 5 events (excluded: noisy baseline "
          f"{n_excluded[k]}, lags without rest samples {n_gap[k]}), events per pair median {np.median(n_events[k]):.0f}; "
          f"population rate at baseline {pop[LAGS < -50].mean():.2f} spikes/s, minimum in 0 to 300 ms {pop[win].min():.2f} "
          f"({100 * (pop[win].min() / pop[LAGS < -50].mean() - 1):.1f}%); below the corrected band: {bool(below.any())}, "
          f"from {first} ms; above it: {bool(above.any())}")
    print(f"   pairs with a change in 0 to 100 ms: {int(sig.drop('none', errors='ignore').sum())} of {len(pair_sig[k])} "
          f"({100 * sig.drop('none', errors='ignore').sum() / len(pair_sig[k]):.1f}%); decreases {int(sig.get('decrease', 0))}, "
          f"increases {int(sig.get('increase', 0))}; the same test on the shifted controls: {100 * np.mean(pair_sig_ctrl[k]):.1f}%")
    # Variant, as Fig5_S14.m: each pair in percent of its baseline, band per lag from the controls, lower
    # edge shifted by the 2.5th percentile over controls of their lowest excursion from 0 to +800 ms.
    pop_p = np.mean(ta_pct[k], 0)
    cp = ctrl_pct_sum[k] / len(ta_pct[k])
    down = np.percentile(cp, 2.5, axis=0)
    post800 = LAGS >= 0
    down_mc = down + np.percentile((cp[:, post800] - down[post800]).min(1), 2.5)
    below_p = post800 & (pop_p < down_mc)
    print(f"   variant as Fig5_S14.m (percent of baseline): minimum in 0 to 300 ms {pop_p[win].min():.1f}%, "
          f"below the band in 0 to 300 ms: {bool((below_p & win).any())}; first lag below it in 0 to 800 ms: "
          f"{LAGS[below_p][0] if below_p.any() else None}")
    ax.plot(LAGS, pop, "k"); ax.fill_between(LAGS, lo_corr, hi_corr, color="0.8")
    ax.set_title(k); ax.set_xlabel("ms from event")
axes[0].set_ylabel("VLa rate (spikes/s)")
print("All plotted values are inside the axes.")
started = time.time()
Fp = fig5_run(pairs_paper, np.random.default_rng(26))
print(f"\nVariant, the paper's records ({time.time() - started:.0f} s); mean burst duration {np.nanmean(Fp['burst_stats'].burst_ms):.1f} ms:")
for k in Fp["ta"]:
    pop, lo_corr, hi_corr = population_band(Fp["ta"][k], Fp["ctrl_sum"][k])
    below = win & (pop < lo_corr)
    sig = pd.Series(Fp["pair_sig"][k])
    print(f"  {k}: {len(Fp['ta'][k])} pairs (excluded: noisy {Fp['n_excluded'][k]}, lags without rest {Fp['n_gap'][k]}); "
          f"minimum in 0 to 300 ms {100 * (pop[win].min() / pop[LAGS < -50].mean() - 1):.1f}%; below the corrected band: "
          f"{bool(below.any())}; pairs with a change in 0 to 100 ms {100 * (sig != 'none').mean():.1f}% "
          f"(shifted controls {100 * np.mean(Fp['pair_sig_ctrl'][k]):.1f}%)")

Bursts and pauses in 45 s
GPi units: 144; rest per unit median 312 s; bursts per minute of rest median 17.75; pauses per minute median 11.66
Mean burst duration 67.1 ms (paper 24 ms); mean pause duration 145.3 ms (paper 195 ms)
Variant, surprise with log10: mean burst duration 254.0 ms; mean pause duration 237.0 ms
burst onset: 150 pairs of 307 with at least 5 events (excluded: noisy baseline 157, lags without rest samples 0), events per pair median 224; population rate at baseline 21.26 spikes/s, minimum in 0 to 300 ms 20.73 (-2.5%); below the corrected band: False, from None ms; above it: False
   pairs with a change in 0 to 100 ms: 76 of 150 (50.7%); decreases 43, increases 33; the same test on the shifted controls: 47.5%
   variant as Fig5_S14.m (percent of baseline): minimum in 0 to 300 ms 98.1%, below the band in 0 to 300 ms: False; first lag below it in 0 to 800 ms: None
burst offset: 150 pairs of 307 with at least 5 events (excluded: noisy baseline 157, lags without rest samples 0), events per pair median 241; population rate at baseline 21.20 spikes/s, minimum in 0 to 300 ms 20.80 (-1.9%); below the corrected band: False, from None ms; above it: False
   pairs with a change in 0 to 100 ms: 75 of 150 (50.0%); decreases 38, increases 37; the same test on the shifted controls: 47.8%
   variant as Fig5_S14.m (percent of baseline): minimum in 0 to 300 ms 98.4%, below the band in 0 to 300 ms: False; first lag below it in 0 to 800 ms: None
pause onset: 110 pairs of 253 with at least 5 events (excluded: noisy baseline 143, lags without rest samples 0), events per pair median 208; population rate at baseline 22.18 spikes/s, minimum in 0 to 300 ms 21.76 (-1.9%); below the corrected band: False, from None ms; above it: False
   pairs with a change in 0 to 100 ms: 49 of 110 (44.5%); decreases 28, increases 21; the same test on the shifted controls: 47.0%
   variant as Fig5_S14.m (percent of baseline): minimum in 0 to 300 ms 96.9%, below the band in 0 to 300 ms: False; first lag below it in 0 to 800 ms: None
All plotted values are inside the axes.

Variant, the paper's records (39 s); mean burst duration 106.5 ms:
  burst onset: 131 pairs (excluded: noisy 163, lags without rest 0); minimum in 0 to 300 ms -1.7%; below the corrected band: False; pairs with a change in 0 to 100 ms 51.9% (shifted controls 47.3%)
  burst offset: 132 pairs (excluded: noisy 162, lags without rest 0); minimum in 0 to 300 ms -2.5%; below the corrected band: False; pairs with a change in 0 to 100 ms 54.5% (shifted controls 47.8%)
  pause onset: 86 pairs (excluded: noisy 159, lags without rest 0); minimum in 0 to 300 ms -2.0%; below the corrected band: False; pairs with a change in 0 to 100 ms 57.0% (shifted controls 46.8%)

v = tr[tr.valid].copy()
v["key"] = v.session.map({s["session_id"]: s["key"] for s in sessions})
per_all = v.groupby(["subject", "session", "key", "target"]).agg(rt=("rt", "median"), mt=("mt", "median"), n=("rt", "size")).reset_index()

def behaviour(per, label):
    print(f"{label}: {per.session.nunique()} files ({per.groupby('subject').session.nunique().to_dict()})")
    print(per.groupby(["subject", "target"])[["rt", "mt"]].agg(["median", "count"]).round(3).to_string())
    for col in ("rt", "mt"):
        d = per.dropna(subset=[col])
        a = sm.stats.anova_lm(ols(f"{col} ~ C(subject) * C(target)", data=d).fit(), typ=2)
        print(f"  two-way ANOVA of {col} (type II), {len(d)} file-directions, residual df {a.loc['Residual', 'df']:.0f}: "
              + "; ".join(f"{k.replace('C(', '').replace(')', '')} F = {a.loc[k, 'F']:.2f}, p = {a.loc[k, 'PR(>F)']:.3g}" for k in a.index[:-1]))

behaviour(per_all[per_all.key.isin(set(units.key))], "Files with units of the dandiset rule")
print()
behaviour(per_all, "Variant, all pre-MPTP files with trials")
print()
behaviour(per_all[per_all.key.isin({f for x in paper.fname for f in str(x).split("/")})], "Variant, the files of the paper's records")
print("\nVariant, trials pooled over all files: RT animal p = %.3g, direction p = %.3g" % (
    stats.ranksums(v[v.subject == "Gaia"].rt, v[v.subject == "Isis"].rt).pvalue,
    stats.ranksums(v[v.target == "Left"].rt, v[v.target == "Right"].rt).pvalue))

Files with units of the dandiset rule: 165 files ({'Gaia': 91, 'Isis': 74})
                   rt           mt      
               median count median count
subject target                          
Gaia    Left    0.242    91  0.216    20
        Right   0.238    91  0.273    20
Isis    Left    0.239    74  0.174    74
        Right   0.245    74  0.195    74
  two-way ANOVA of rt (type II), 330 file-directions, residual df 326: subject F = 0.02, p = 0.89; target F = 0.69, p = 0.408; subject:target F = 2.48, p = 0.116
  two-way ANOVA of mt (type II), 188 file-directions, residual df 184: subject F = 526.30, p = 7.31e-56; target F = 270.38, p = 5.76e-38; subject:target F = 33.39, p = 3.16e-08

Variant, all pre-MPTP files with trials: 250 files ({'Gaia': 119, 'Isis': 131})
                   rt           mt      
               median count median count
subject target                          
Gaia    Left    0.243   119  0.218    32
        Right   0.242   119  0.273    32
Isis    Left    0.235   131  0.175   131
        Right   0.241   131  0.204   131

  two-way ANOVA of rt (type II), 500 file-directions, residual df 496: subject F = 15.19, p = 0.000111; target F = 3.23, p = 0.073; subject:target F = 4.24, p = 0.04
  two-way ANOVA of mt (type II), 326 file-directions, residual df 322: subject F = 691.86, p = 3.41e-82; target F = 483.80, p = 4.17e-66; subject:target F = 30.86, p = 5.81e-08

Variant, the files of the paper's records: 162 files ({'Gaia': 97, 'Isis': 65})
                   rt           mt      
               median count median count
subject target                          
Gaia    Left    0.243    97  0.217    28
        Right   0.242    97  0.273    28
Isis    Left    0.243    65  0.174    65
        Right   0.249    65  0.195    65
  two-way ANOVA of rt (type II), 324 file-directions, residual df 320: subject F = 3.79, p = 0.0525; target F = 0.45, p = 0.503; subject:target F = 2.77, p = 0.097
  two-way ANOVA of mt (type II), 186 file-directions, residual df 182: subject F = 1862.70, p = 1.55e-97; target F = 620.06, p = 1.62e-60; subject:target F = 136.49, p = 6.88e-24

Variant, trials pooled over all files: RT animal p = 2.47e-59, direction p = 2.23e-13

rows = []
for s in sessions:
    e = s["electrodes"]
    e = e[e.location.map(g.area).isin(["VLa", "VLp"])]
    rows += [{"subject": s["subject"], "area": g.area(r.location), "y": r.acx_y, "z": r.acx_z} for r in e.itertuples()]
loc = pd.DataFrame(rows).drop_duplicates()
f = g._open(sessions[0]["url"])
for c in ("acx_y", "acx_z"):
    print(f"Description of the electrode column {c}:", g._s(f["general/extracellular_ephys/electrodes"][c].attrs.get("description")))
f.close()
print("Electrodes labelled VLa or VLp without coordinates:", int(loc[["y", "z"]].isna().any(axis=1).sum()), "of", len(loc))
loc = loc.dropna()
print("Distinct electrode positions labelled VLa or VLp (mm from the anterior commissure; y posterior negative, z ventral positive):")
print(loc.groupby(["subject", "area"])[["y", "z"]].agg(["median", "min", "max", "count"]).round(2).to_string())
for subj, d in loc.groupby("subject"):
    va, vp = d[d.area == "VLa"], d[d.area == "VLp"]
    if len(va) and len(vp):
        print(f"{subj}: VLa anterior of VLp, one-sided rank sum p = {stats.ranksums(va.y, vp.y, alternative='greater').pvalue:.2g}; "
              f"VLa dorsal of VLp (smaller z), p = {stats.ranksums(va.z, vp.z, alternative='less').pvalue:.2g}; "
              f"VLa ventral of VLp (larger z), p = {stats.ranksums(va.z, vp.z, alternative='greater').pvalue:.2g}")

Description of the electrode column acx_y: Antero-posterior coordinate relative to the anterior commissure crossing, posterior negative, in millimeters.
Description of the electrode column acx_z: Dorso-ventral coordinate relative to the anterior commissure crossing, ventral positive, in millimeters.
Electrodes labelled VLa or VLp without coordinates: 1 of 617
Distinct electrode positions labelled VLa or VLp (mm from the anterior commissure; y posterior negative, z ventral positive):
                  y                        z                  
             median   min   max count median   min   max count
subject area                                                  
Gaia    VLa   -7.44 -9.15 -5.80   135   1.93 -1.62  4.16   135
        VLp   -8.28 -9.05 -6.98    42   0.50 -2.17  2.49    42
Isis    VLa   -6.41 -7.71 -4.30   259   4.03  1.71  6.91   259
        VLp   -7.82 -9.78 -6.28   180   2.84  0.31  4.26   180
Gaia: VLa anterior of VLp, one-sided rank sum p = 6.4e-10; VLa dorsal of VLp (smaller z), p = 1; VLa ventral of VLp (larger z), p = 1.8e-11
Isis: VLa anterior of VLp, one-sided rank sum p = 5.5e-58; VLa dorsal of VLp (smaller z), p = 1; VLa ventral of VLp (larger z), p = 2e-27

print(f"Notebook ran in {(time.time() - NB_START) / 60:.1f} min")

Notebook ran in 16.7 min

Claim	Data	Status
unit-link	with inputs	tested
table1-unit-counts	dandiset	tested
table1-rest-rate	dandiset	tested
table1-ap-width	dandiset	tested
hold-period-ramps	dandiset	tested
mvt-responsive-fraction	dandiset	tested
increases-most-common	dandiset	tested
monophasic-more-in-vla	dandiset	tested
integrated-change-positive	dandiset	tested
fig3a-vla-leads-gpi	dandiset	tested
fig3b-gpi-decr-lag-vla-incr	dandiset	tested
fig3c-vla-decr-lead-gpi-incr	dandiset	tested
gating-conflict-fraction	dandiset	tested
directional-units	dandiset	tested
fig4-ccf-rare	dandiset	tested
fig4-noise-corr-rare	dandiset	tested
fig5-vla-decrease-after-gpi-bursts	dandiset	tested
behavior-rt-mt	dandiset	tested
figs3-emg-onset	elsewhere	blocked: methods underspecified
fig1-vla-localization	dandiset	blocked: data missing
fig6-simulation	elsewhere	blocked: other

Claim	Data	Outcome
unit-link	with inputs	reproduced
table1-unit-counts	dandiset	not reproduced
table1-rest-rate	dandiset	reproduced
table1-ap-width	dandiset	reproduced
hold-period-ramps	dandiset	partially reproduced
mvt-responsive-fraction	dandiset	reproduced
increases-most-common	dandiset	partially reproduced
monophasic-more-in-vla	dandiset	not reproduced
integrated-change-positive	dandiset	partially reproduced
fig3a-vla-leads-gpi	dandiset	not reproduced
fig3b-gpi-decr-lag-vla-incr	dandiset	reproduced
fig3c-vla-decr-lead-gpi-incr	dandiset	not reproduced
gating-conflict-fraction	dandiset	reproduced
directional-units	dandiset	partially reproduced
fig4-ccf-rare	dandiset	not reproduced
fig4-noise-corr-rare	dandiset	not reproduced
fig5-vla-decrease-after-gpi-bursts	dandiset	not reproduced
behavior-rt-mt	dandiset	reproduced
figs3-emg-onset	elsewhere	blocked: methods underspecified
fig1-vla-localization	dandiset	blocked: data missing
fig6-simulation	elsewhere	blocked: other

Reproducing Schwab et al. (2020) from dandiset 000947¶

Data check¶

Unit selection and peri-movement responses¶

unit-link¶

table1-unit-counts¶

table1-rest-rate¶

table1-ap-width¶

hold-period-ramps¶

mvt-responsive-fraction¶

increases-most-common¶

monophasic-more-in-vla¶

integrated-change-positive¶

fig3a-vla-leads-gpi¶

fig3b-gpi-decr-lag-vla-incr¶

fig3c-vla-decr-lead-gpi-incr¶

gating-conflict-fraction¶

directional-units¶

Pairs of simultaneously recorded GPi and VLa units¶

fig4-ccf-rare¶

fig4-noise-corr-rare¶

fig5-vla-decrease-after-gpi-bursts¶

behavior-rt-mt¶

figs3-emg-onset¶

fig1-vla-localization¶

fig6-simulation¶

Summary¶