import time
import warnings
from collections import Counter

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import seqlick as sl

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

import os

started = time.time()
print(f"CPU cores: {os.cpu_count()}; workers for parallel steps: {sl.cpu_workers()}")
assets = sl.list_assets()
print(f"dandiset {sl.DANDISET_ID} version {sl.VERSION}: {len(assets)} files")
# Trial tables of every file, and unit quality where there are units (spike times are not kept).
tables = sl.load_many(assets, workers=sl.cpu_workers(), video=False, spikes=False)
print(f"read the trial tables of {len(tables)} files in {time.time() - started:.0f} s")
FIG3 = ["S1TJ", "M1TJ", "ALM", "S1BF"]
EPHYS57 = FIG3 + ["M1B", "S1L", "ZZ"]
for t in tables:
    kind, area = sl.group_label(t["description"])
    t["kind"], t["area"] = kind, area
full_paths = {t["path"] for t in tables if t["kind"] == "ephys" and t["area"] in FIG3 + ["ZZ"]}
t1 = time.time()
sessions = sl.load_many([a for a in assets if a["path"] in full_paths], workers=sl.cpu_workers())
for s in sessions:
    s["kind"], s["area"] = sl.group_label(s["description"])
print(f"read video and spikes of {len(sessions)} recording sessions in {time.time() - t1:.0f} s")

CPU cores: 4; workers for parallel steps: 16

dandiset 000239 version 0.230607.1752: 754 files

read the trial tables of 754 files in 31 s

read video and spikes of 41 recording sessions in 21 s

info = pd.DataFrame([{
    "path": t["path"], "subject": t["subject"], "group": t["group"], "kind": t["kind"], "area": t["area"],
    "trials": t["trials"]["n"], "units": t["n_units"], "size_mb": t["size"] / 1e6,
    "behavior_file": t["path"].endswith("_behavior.nwb")} for t in tables])
print("Files, subjects and units by session description (the name of the lab's data folder):")
print(info.groupby("group").agg(files=("path", "size"), subjects=("subject", "nunique"), trials=("trials", "sum"),
                                files_with_units=("units", lambda u: int((u > 0).sum())), units=("units", "sum"),
                                behavior_files=("behavior_file", "sum")).to_string())
print(f"\nAll files: {len(info)}; subjects: {info.subject.nunique()}; files with units: {(info.units > 0).sum()}; "
      f"units: {info.units.sum()}")
cols = Counter(tuple(t["trials"].keys()) for t in tables)
print("\nTrial table columns (sets) and how many files have them:")
for c, n in cols.most_common():
    print(n, [x for x in c if x != "n"])

Files, subjects and units by session description (the name of the lab's data folder):
                                    files  subjects  trials  files_with_units  units  behavior_files
group                                                                                               
Data behav air                          6         6    1454                 0      0               6
Data behav earplug                     10         5    1972                 0      0              10
Data behav numbing                     11         6    1643                 0      0              11
Data ephys ALM                         13         8    3457                13    450              13
Data ephys M1B                          4         2    1096                 4    117               4
Data ephys M1TJ                         9         4    2093                 9    328               9
Data ephys S1BF                         5         3    1592                 5     78               5
Data ephys S1L                         12         3    2508                12    453              12
Data ephys S1TJ                         8         4    2375                 8    265               8
Data ephys ZZ                           6         3    1549                 6    226               6
Data ephys others                      22         8    5696                22    742              22
Data learning                         574        15  141232                 0      0               0
Data opto VGAT-CRE Ai32 2s 2.5V        25         5    8222                 0      0              25
Data opto VGAT-CRE Ai32 2s 5V          35         7   10023                 0      0              35
Data opto VGAT-CRE Ai32 efficiency     14         4    3207                14    414              14

All files: 754; subjects: 33; files with units: 93; units: 3073

Trial table columns (sets) and how many files have them:
574 ['start_time', 'stop_time', 'cue', 'posIndex', 'water', 'opto', 'nolickITI', 'waterDelay', 'seqId', 'optoDur1', 'optoDur2', 'optoFreq1', 'optoFreq2', 'optoMod1', 'optoMod2', 'cue_times', 'posIndex_times', 'water_times', 'opto_times', 'lickOn_times', 'lickOff_times', 'cueOff_times', 'waterOff_times', 'waterTrig_times', 'endConsump_times']
163 ['start_time', 'stop_time', 'cue', 'posIndex', 'water', 'opto', 'cue_times', 'posIndex_times', 'water_times', 'opto_times', 'lickOn_times', 'lickOff_times']
17 ['start_time', 'stop_time', 'cue', 'posIndex', 'water', 'opto', 'nolickITI', 'cue_times', 'posIndex_times', 'water_times', 'opto_times', 'lickOn_times', 'lickOff_times']

# Recordings deposited more than once: compare the first 20 touch times of every pair of files.
keys = {}
for t in tables:
    touches = np.concatenate(t["trials"]["lickOn_times"]) if t["trials"]["n"] else np.array([])
    keys.setdefault(tuple(np.round(touches[:20], 4)), []).append(t)
dups = [v for v in keys.values() if len(v) > 1]
print(f"Groups of files with identical first 20 touch times: {len(dups)}")
for v in dups:
    for t in v:
        tt = np.concatenate(t["trials"]["lickOn_times"])
        print(f"  {t['path']}  {t['description']}  trials {t['trials']['n']}  touches {len(tt)}  "
              f"touch range {tt.min():.1f} to {tt.max():.1f} s")
    a, b = [np.concatenate(t["trials"]["lickOn_times"]) for t in v[:2]]
    n = min(len(a), len(b))
    print(f"  identical touch times over the first {n} touches: {np.mean(np.isclose(a[:n], b[:n])):.3f} of them; "
          f"the longer file's extra touches start at {max(a, b, key=len)[n]:.1f} s")
dup_paths = {t["path"] for v in dups for t in v}

Groups of files with identical first 20 touch times: 1
  sub-MX181002/sub-MX181002_ses-20181002T153057.nwb  Data learning - MX181002 2018-10-02 15:30:57  trials 259  touches 4028  touch range 2266762.2 to 2269692.5 s
  sub-MX181003/sub-MX181003_ses-20181002T153057.nwb  Data learning - MX181003 2018-10-02 15:30:57  trials 499  touches 8441  touch range 2266762.2 to 2274595.9 s
  identical touch times over the first 4028 touches: 1.000 of them; the longer file's extra touches start at 2270389.9 s

# Clocks of the trial table: start_time against cue_times, per kind of file.
rows = []
for t in tables:
    tr = t["trials"]
    if not tr["n"]:
        continue
    touch = np.concatenate(tr["lickOn_times"])
    rows.append({"kind": t["kind"], "behavior_file": t["path"].endswith("_behavior.nwb"),
                 "start_minus_cue_max": np.nanmax(np.abs(tr["start_time"] - tr["cue_times"])),
                 "first_start": tr["start_time"][0], "first_cue": tr["cue_times"][0],
                 "touch_in_trials": np.mean((touch >= tr["start_time"][0]) & (touch <= tr["stop_time"][-1]))})
clk = pd.DataFrame(rows)
print("Largest |start_time - cue_times| per file, and the fraction of touches inside the span of the trial table:")
print(clk.groupby(["kind", "behavior_file"]).agg(files=("kind", "size"),
      start_minus_cue_median=("start_minus_cue_max", "median"), start_minus_cue_max=("start_minus_cue_max", "max"),
      first_cue_median=("first_cue", "median"), touch_in_trials_median=("touch_in_trials", "median")).to_string())
learn = [t for t in tables if t["kind"] == "learning"]
print("\nseqId in learning files, first value:", repr(learn[0]["trials"]["seqId"][0])[:160])
print("Distinct seqId values over all learning files:", len({v for t in learn for v in t["trials"]["seqId"]}))

Largest |start_time - cue_times| per file, and the fraction of touches inside the span of the trial table:
                                  files  start_minus_cue_median  start_minus_cue_max  first_cue_median  touch_in_trials_median
kind               behavior_file                                                                                              
behav              True              27                   0.000                0.000      0.000000e+00                     1.0
ephys              True              79                   0.000                0.001      1.843993e+01                     1.0
learning           False            574             1023587.282          4274364.864      1.023587e+06                     0.0
opto VGAT-CRE Ai32 True              74                   0.001                0.001      0.000000e+00                     1.0

seqId in learning files, first value: np.str_("[(b'', b'MCOS', b'categorical', array([3707764736,          2,          1,          1,          1,\n                 1], dtype=uint32))                
Distinct seqId values over all learning files: 499

# Trial patterns (posIndex: the port positions visited after each port-moving touch).
pat = Counter()
for t in tables:
    for k in t["trials"]["posIndex"]:
        pat[(t["kind"], sl.trial_type(k)[0])] += 1
print("Trials by kind of file and trial type:")
print(pd.Series(pat).unstack(fill_value=0).to_string())
print("\nposIndex patterns not classified as standard, backtracking or zigzag (top 8):")
print(Counter(k for t in tables for k in t["trials"]["posIndex"] if sl.trial_type(k)[0] == "other").most_common(8))
# Opto and condition labels in the opto and control sessions.
for kind in ["opto VGAT-CRE Ai32", "behav"]:
    sub = [t for t in tables if t["kind"] == kind]
    lab = Counter()
    for t in sub:
        tr = t["trials"]
        lab["files"] += 1
        lab["opto column non-NaN"] += int(np.sum(~np.isnan(tr["opto"])))
        lab["opto_times non-NaN"] += int(np.sum(~np.isnan(tr["opto_times"])))
    print(f"\n{kind}: {dict(lab)}; areas in descriptions: {Counter(t['area'] for t in sub)}")
print("Columns that could name the inhibited region or the control condition:",
      sorted({c for t in tables if t["kind"] in ("opto VGAT-CRE Ai32", "behav") for c in t["trials"]} -
             {"n", "start_time", "stop_time", "cue", "posIndex", "water", "opto", "cue_times", "posIndex_times",
              "water_times", "opto_times", "lickOn_times", "lickOff_times", "nolickITI"}))

Trials by kind of file and trial type:
                    backtrack  other  standard  zigzag
behav                     819     10      4240       0
ephys                    3484    357     14976    1549
learning                10370   8644    115115    7103
opto VGAT-CRE Ai32        546    153     20753       0

posIndex patterns not classified as standard, backtracking or zigzag (top 8):
[(np.str_('[3 6 3 0 3 6 3 0 3 6]'), 2280), (np.str_('[3 0 3 6 3 0 3 6 3 0]'), 2272), (np.str_('[1 2 3 5 6]'), 1031), (np.str_('[5 4 3 1 0]'), 974), (np.str_('[1 2 3 2 1 2 3 2 1 2 3 4]'), 970), (np.str_('[3 2 1 2 3 2 1 2 3 2 1 0]'), 968), (np.str_('[3 6 3 0 3 6]'), 271), (np.str_('[3 0 3 6 3 0]'), 271)]

opto VGAT-CRE Ai32: {'files': 74, 'opto column non-NaN': 6221, 'opto_times non-NaN': 6583}; areas in descriptions: Counter({'2s 5V': 35, '2s 2.5V': 25, 'efficiency': 14})

behav: {'files': 27, 'opto column non-NaN': 0, 'opto_times non-NaN': 0}; areas in descriptions: Counter({'numbing': 11, 'earplug': 10, 'air': 6})
Columns that could name the inhibited region or the control condition: []

# Units: quality, sortedness and duplicates.
uq = pd.DataFrame([{"group": t["group"], "area": t["area"], "kind": t["kind"], "path": t["path"],
                    "rpv": q[0], "contam": q[1], "spikes": n}
                   for t in tables for q, n in zip(t["unit_quality"], t["n_spikes"])])
print("Units with unsorted spike times:", sum(t["unsorted"] for t in tables),
      "; duplicated spike times within a unit:", sum(t["duplicate_spikes"] for t in tables))
dd = np.concatenate([t["duplicate_boundary_distance"] for t in tables])
print(f"Distance of duplicated spike times to the nearest trial start or stop (s): median {np.median(dd):.2f}, "
      f"5th percentile {np.percentile(dd, 5):.3f}; within 10 ms of a boundary: {int(np.sum(dd < 0.01))} of {len(dd)}")
print("Files with duplicated spike times:", sum(t["duplicate_spikes"] > 0 for t in tables))
print("Units with fewer than 2 spikes:", int((uq.spikes < 2).sum()))
print("ISI violation rate (2.5 ms) and contamination, quantiles over all units:")
print(uq[["rpv", "contam", "spikes"]].quantile([0, 0.25, 0.5, 0.75, 0.95, 1]).round(4).to_string())
print("Unit table columns:", Counter(tuple(t["unit_columns"]) for t in tables if t["n_units"]))
print("Files with an electrodes table or extracellular ephys metadata:", sum(t["has_electrodes"] for t in tables))
print("\nOpto sessions per mouse and description (the inhibited region is not stated):")
print(pd.crosstab(info[info.kind == "opto VGAT-CRE Ai32"].subject, info[info.kind == "opto VGAT-CRE Ai32"].area).to_string())
print("\nControl sessions per mouse (which session or block is the manipulation is not stated):")
print(pd.crosstab(info[info.kind == "behav"].subject, info[info.kind == "behav"].area).to_string())

Units with unsorted spike times: 0 ; duplicated spike times within a unit: 1446
Distance of duplicated spike times to the nearest trial start or stop (s): median 2.90, 5th percentile 0.276; within 10 ms of a boundary: 2 of 1446
Files with duplicated spike times: 35
Units with fewer than 2 spikes: 0
ISI violation rate (2.5 ms) and contamination, quantiles over all units:
         rpv  contam    spikes
0.00  0.0000  0.0000      84.0
0.25  0.0003  0.0071    5968.0
0.50  0.0011  0.0285   14390.0
0.75  0.0033  0.0952   29235.0
0.95  0.0333  1.0000   65409.4
1.00  0.5968  1.0000  336559.0
Unit table columns: Counter({('spike_times', 'obs_intervals'): 93})
Files with an electrodes table or extracellular ephys metadata: 0

Opto sessions per mouse and description (the inhibited region is not stated):
area      2s 2.5V  2s 5V  efficiency
subject                             
MX180803        5      5           0
MX180804        5      5           0
MX181002        5      5           0
MX181003        0      5           0
MX190101        0      0           4
MX190103        0      5           0
MX190201        0      0           4
MX190401        0      0           2
MX190402        0      0           4
MX200101        5      5           0
MX200102        5      5           0

Control sessions per mouse (which session or block is the manipulation is not stated):
area      air  earplug  numbing
subject                        
MX190201    1        0        0
MX190602    1        2        0
MX200102    0        0        2
MX200103    0        0        2
MX210203    0        0        2
MX210502    0        0        1
MX210601    0        0        2
MX210602    0        0        2
VC010102    1        2        0
VC010103    1        2        0
WO010401    1        2        0
WO010402    1        2        0

# Everything in the opto and control files that could name the inhibited region or the condition:
# free-text fields, general, stimulus and acquisition, and the opto and force series' attributes.
cond_assets = [a for a in assets if sl.group_label(a["description"])[0] in ("opto VGAT-CRE Ai32", "behav")]
meta = sl.run_parallel(sl.file_metadata, cond_assets, workers=sl.cpu_workers())
mdf = pd.DataFrame([{"group": sl.group_label(a["description"])[1], "session_description": m.get("session_description"),
                     "other_text": {k: m[k] for k in ("experiment_description", "notes", "protocol", "stimulus_notes") if k in m},
                     "general": tuple(m["general"]), "stimulus": tuple(m["stimulus"]), "acquisition": tuple(m["acquisition"]),
                     "series": tuple(sorted((k, v["unit"], v["conversion"], v["offset"], v["description"], v["comments"])
                                            for k, v in m["series"].items()))} for a, m in zip(cond_assets, meta)])
print(f"Files read: {len(mdf)}")
print("Session descriptions as stored in the files (first of each group):")
print(mdf.groupby("group").session_description.first().to_string())
print("Descriptions ending in ' se' (cut short in the file itself, not only in the archive):",
      int(mdf.session_description.str.endswith(" se").sum()), "of", len(mdf))
print("Other free-text fields (experiment_description, notes, protocol, stimulus_notes) present in any file:",
      sorted({k for d in mdf.other_text for k in d}))
for col in ("general", "stimulus", "acquisition", "series"):
    print(f"\nDistinct contents of {col} (number of files):")
    for v, n in Counter(mdf[col]).most_common():
        print(f"  {n}: {v}")
print("Objects inside the groups of stimulus (presentation, templates), over all files:",
      sorted({n for v in mdf.stimulus for n in v if "/" in n}))
opto_vals, opto_nan = Counter(), 0
for t in tables:
    if t["kind"] == "opto VGAT-CRE Ai32":
        o = t["trials"]["opto"]
        opto_vals.update(np.unique(o[~np.isnan(o)]).tolist())
        opto_nan += int(np.isnan(o).any())
print(f"\nDistinct values of the trial column opto over the 74 opto files (files holding each): {dict(opto_vals)}; "
      f"files with NaN in it: {opto_nan}")
print("Values of opto1 and opto2 in the first light trial of each file (first 10), distinct sets:",
      len({tuple(m["series"]["opto1"].get("values_first_light_trial", [])) for m in meta if "opto1" in m["series"]}),
      "; largest value over files:",
      max(max(m["series"][k].get("values_first_light_trial", [0])) for m in meta for k in ("opto1", "opto2") if k in m["series"]))

Files read: 101
Session descriptions as stored in the files (first of each group):
group
2s 2.5V       Data opto VGAT-CRE Ai32 2s 2.5V - MX180803 201...
2s 5V         Data opto VGAT-CRE Ai32 2s 5V - MX180803 2018-...
air                     Data behav air - MX190201 2019-10-08 se
earplug             Data behav earplug - MX190602 2019-12-03 se
efficiency    Data opto VGAT-CRE Ai32 efficiency - MX190101 ...
numbing             Data behav numbing - MX200102 2021-01-17 se
Descriptions ending in ' se' (cut short in the file itself, not only in the archive): 101 of 101
Other free-text fields (experiment_description, notes, protocol, stimulus_notes) present in any file: []

Distinct contents of general (number of files):
  101: ('experimenter', 'institution', 'lab', 'subject', 'subject/date_of_birth', 'subject/description', 'subject/genotype', 'subject/sex', 'subject/species', 'subject/subject_id')

Distinct contents of stimulus (number of files):
  101: ('presentation', 'templates')

Distinct contents of acquisition (number of files):
  101: ()

Distinct contents of series (number of files):
  74: (('lickportH', 'volt', 1.0, 0.0, 'no description', 'no comments'), ('opto1', 'categorical', 1.0, 0.0, 'no description', 'no comments'), ('opto2', 'categorical', 1.0, 0.0, 'no description', 'no comments'), ('tubeH', 'volt', 1.0, 0.0, 'no description', 'no comments'), ('tubeV', 'volt', 1.0, 0.0, 'no description', 'no comments'))
  21: (('lickportH', 'volt', 1.0, 0.0, 'no description', 'no comments'), ('tubeH', 'volt', 1.0, 0.0, 'no description', 'no comments'), ('tubeV', 'volt', 1.0, 0.0, 'no description', 'no comments'))
  6: ()
Objects inside the groups of stimulus (presentation, templates), over all files: []

Distinct values of the trial column opto over the 74 opto files (files holding each): {0.0: 56, 1.0: 70, 2.0: 56}; files with NaN in it: 74
Values of opto1 and opto2 in the first light trial of each file (first 10), distinct sets: 28 ; largest value over files: 0.01600000075995922

# Video series: rates, lengths and the scale of L and theta.
vr = pd.DataFrame([{"area": s["area"], "rate": v["rate"], "dur": len(v["out"]) / v["rate"],
                    "start_minus_cue": v["start"] - s["trials"]["cue_times"][i]}
                   for s in sessions for i, v in enumerate(s["video"])])
print("Video series per trial: sampling rate (Hz), duration (s), start minus cue (s):")
print(vr.describe().round(3).to_string())
lk = pd.DataFrame(np.concatenate([np.concatenate([l for l in s["licks"] if len(l)]) for s in sessions]),
                  columns=["on", "off", "mid", "Lmax", "theta_shoot"])
print(f"\nVideo licks: {len(lk)}; duration (ms) quantiles 1/50/99%: "
      f"{np.percentile((lk.off - lk.on) * 1000, [1, 50, 99]).round(1)}")
print(f"Lmax (mm, at 32 um/pixel) quantiles 1/50/99%: {np.nanpercentile(lk.Lmax, [1, 50, 99]).round(2)}  "
      "(Figure 3a: L up to about 2.7 mm)")
print(f"theta_shoot (deg) quantiles 1/50/99%: {np.nanpercentile(lk.theta_shoot, [1, 50, 99]).round(1)}  "
      "(Figure 3a: about -30 to 30 degrees)")
print("Units of the behavioural series as stored: tongue landmarks 'pixels', lick port and force sensors 'volt'.")

Video series per trial: sampling rate (Hz), duration (s), start minus cue (s):
            rate        dur  start_minus_cue
count  11066.000  11066.000        11066.000
mean     400.143      4.564           -0.000
std        0.125      1.431            0.000
min      399.744      2.927           -0.001
25%      400.000      3.790            0.000
50%      400.176      4.064            0.000
75%      400.256      4.637            0.000
max      400.416     10.000            0.000

Video licks: 210903; duration (ms) quantiles 1/50/99%: [ 20.   65.  119.9]
Lmax (mm, at 32 um/pixel) quantiles 1/50/99%: [0.7  2.51 3.55]  (Figure 3a: L up to about 2.7 mm)
theta_shoot (deg) quantiles 1/50/99%: [-38.4   0.3  47.1]  (Figure 3a: about -30 to 30 degrees)
Units of the behavioural series as stored: tongue landmarks 'pixels', lick port and force sensors 'volt'.

# Do the video series share the clock of the touch times? For each session, the fraction of touch
# onsets that fall in a video frame classified tongue-out, with the touch times shifted by d.
shifts = np.round(np.arange(-0.2, 0.2001, 0.0025), 4)
vid_rows = []
for s in sessions:
    hits = np.zeros(len(shifts))
    count = 0
    for i, v in enumerate(s["video"]):
        touch = np.asarray(s["trials"]["lickOn_times"][i])
        t_end = v["start"] + len(v["out"]) / v["rate"]
        touch = touch[(touch > v["start"] + 0.25) & (touch < t_end - 0.25)]
        if not len(touch):
            continue
        k = np.round((touch[None, :] + shifts[:, None] - v["start"]) * v["rate"]).astype(int)
        hits += v["out"][np.clip(k, 0, len(v["out"]) - 1)].sum(1)
        count += len(touch)
    frac = hits / count
    vid_rows.append({"session": s["path"].split("/")[1][:32], "area": s["area"], "touches": count,
                     "at_0": frac[shifts == 0][0], "best_shift_ms": 1000 * shifts[np.argmax(frac)],
                     "best": frac.max(), "at_-77ms": frac[np.isclose(shifts, -0.0775)][0],
                     "at_+77ms": frac[np.isclose(shifts, 0.0775)][0]})
vid = pd.DataFrame(vid_rows)
print("Fraction of touches inside tongue-out frames, at zero shift, at the best shift, and half a lick cycle away:")
print(vid.round(3).to_string(index=False))
print(f"\nBest shift over sessions (ms): median {vid.best_shift_ms.median():.1f}, range {vid.best_shift_ms.min():.1f} "
      f"to {vid.best_shift_ms.max():.1f}; at zero shift median {vid.at_0.median():.3f}, lowest {vid.at_0.min():.3f}; "
      f"at +/-77 ms medians {vid['at_-77ms'].median():.3f} and {vid['at_+77ms'].median():.3f}")

Fraction of touches inside tongue-out frames, at zero shift, at the best shift, and half a lick cycle away:
                         session area  touches  at_0  best_shift_ms  best  at_-77ms  at_+77ms
sub-MX170903_ses-20180304T135947  ALM     3327 0.974          -15.0 0.992     0.037     0.103
sub-MX180202_ses-20180504T134624  ALM     1903 1.000           -2.5 1.000     0.011     0.111
sub-MX180203_ses-20180526T163118  ALM     3283 1.000           -2.5 1.000     0.133     0.002
sub-MX180202_ses-20180505T144012  ALM     3766 0.999            5.0 0.999     0.009     0.099
sub-MX180203_ses-20180527T135221  ALM     2977 1.000            0.0 1.000     0.166     0.013
sub-MX180203_ses-20180529T131703  ALM     4426 1.000            0.0 1.000     0.196     0.007
sub-MX180701_ses-20180916T121311  ALM     2522 0.999            2.5 1.000     0.066     0.103
sub-MX180803_ses-20181128T130139 S1BF     4922 1.000            0.0 1.000     0.036     0.002
sub-MX180804_ses-20181125T151623  ALM     3420 1.000            0.0 1.000     0.002     0.096
sub-MX180804_ses-20181128T151110 S1BF     3831 0.998           10.0 1.000     0.008     0.071
sub-MX180804_ses-20181202T124714 S1TJ     3641 1.000            2.5 1.000     0.010     0.052
sub-MX180804_ses-20181201T130007 S1TJ     3060 0.997           10.0 1.000     0.008     0.069
sub-MX180804_ses-20181129T150336 S1BF     4211 0.999            5.0 1.000     0.010     0.052
sub-MX181002_ses-20181224T154648  ALM     5940 1.000            0.0 1.000     0.027     0.008
sub-MX181002_ses-20181229T151032 S1TJ     4932 1.000           -2.5 1.000     0.045     0.007
sub-MX181002_ses-20181228T150236 S1TJ     5553 1.000            0.0 1.000     0.056     0.015
sub-MX181002_ses-20181230T141519 S1TJ     5159 1.000           -2.5 1.000     0.031     0.005
sub-MX181002_ses-20181219T132039 S1BF     4045 1.000            0.0 1.000     0.027     0.010
sub-MX181002_ses-20181218T140101 S1BF     4948 1.000            2.5 1.000     0.028     0.012
sub-MX181002_ses-20181223T145116  ALM     5628 1.000            2.5 1.000     0.034     0.023
sub-MX181101_ses-20190329T105021 S1TJ     2753 1.000            0.0 1.000     0.031     0.182
sub-MX181101_ses-20190401T114039 M1TJ     2572 1.000            0.0 1.000     0.042     0.181
sub-MX181101_ses-20190328T114227 S1TJ     3188 1.000            0.0 1.000     0.023     0.103
sub-MX181101_ses-20190331T103702  ALM     3341 1.000           -2.5 1.000     0.054     0.189
sub-MX181101_ses-20190330T112923  ALM     3424 1.000            2.5 1.000     0.042     0.191
sub-MX181101_ses-20190402T123529 M1TJ     2960 1.000            0.0 1.000     0.048     0.174
sub-MX181301_ses-20190522T105041 M1TJ     2144 1.000            0.0 1.000     0.079     0.074
sub-MX181301_ses-20190524T135431 M1TJ     2696 0.998            7.5 1.000     0.071     0.091
sub-MX181301_ses-20190521T122250 M1TJ     2037 0.997            5.0 1.000     0.080     0.043
sub-MX181302_ses-20190210T163258 M1TJ     3353 0.987            7.5 0.998     0.065     0.018
sub-MX181302_ses-20190212T154624 M1TJ     4182 0.991           10.0 0.998     0.046     0.027
sub-MX190103_ses-20190513T110307 M1TJ     3771 1.000            2.5 1.000     0.031     0.081
sub-MX190103_ses-20190512T114257 M1TJ     4255 0.999            5.0 1.000     0.045     0.081
sub-MX190103_ses-20190515T110643  ALM     2701 1.000            0.0 1.000     0.021     0.092
sub-MX190103_ses-20190514T112816 S1TJ     4108 1.000            0.0 1.000     0.024     0.077
sub-MX200101_ses-20201126T155559   ZZ     5530 0.998            0.0 0.998     0.035     0.001
sub-MX200101_ses-20201127T154735   ZZ     5699 0.998            0.0 0.998     0.035     0.001
sub-MX200105_ses-20201210T175011   ZZ     4737 1.000           -2.5 1.000     0.018     0.003
sub-MX200105_ses-20201212T152521   ZZ     4091 1.000            0.0 1.000     0.027     0.000
sub-MX210301_ses-20210501T160243   ZZ     5088 0.977           10.0 0.994     0.011     0.060
sub-MX210301_ses-20210506T151300   ZZ     5329 0.983            7.5 0.996     0.004     0.054

Best shift over sessions (ms): median 0.0, range -15.0 to 10.0; at zero shift median 1.000, lowest 0.974; at +/-77 ms medians 0.034 and 0.054

# Do the spike times share the clock of the video? As in Extended Data Fig. 5a, tongue length L is
# decoded with the spike times shifted by each lag (positive: spikes taken later than the
# behaviour). This cell also runs the Figure 3 fits used below (one process per session).
import functools
t0 = time.time()
fig3_sessions = [s for s in sessions if s["area"] in FIG3]
LAGS = (-0.1, -0.05, -0.025, 0.0, 0.025, 0.05, 0.1)
f3 = sl.run_parallel(functools.partial(sl.fig3_session, seed=3, shifts=LAGS), fig3_sessions)
print(f"Figure 3 fits for {len(f3)} sessions in {time.time() - t0:.0f} s")
sh = pd.DataFrame([{"area": o["area"], **{f"{1000 * k:+.0f} ms": v for k, v in o["shift_r2"].items()}} for o in f3])
print("Cross-validated R2 of L (5 folds over trials) with spike times shifted, mean over sessions:")
print(sh.groupby("area").mean().round(3).loc[FIG3].to_string())
best = sh.drop(columns="area").idxmax(axis=1)
print("\nLag of the best R2, sessions per lag, S1TJ and M1TJ:")
print(best[sh.area.isin(["S1TJ", "M1TJ"])].value_counts().to_string())

Figure 3 fits for 35 sessions in 321 s
Cross-validated R2 of L (5 folds over trials) with spike times shifted, mean over sessions:
      -100 ms  -50 ms  -25 ms  +0 ms  +25 ms  +50 ms  +100 ms
area                                                         
S1TJ    0.387   0.390   0.486  0.555   0.537   0.451    0.419
M1TJ    0.307   0.374   0.392  0.398   0.457   0.416    0.351
ALM     0.072   0.092   0.086  0.090   0.111   0.097    0.098
S1BF    0.006   0.026   0.020  0.019   0.040   0.029    0.038

Lag of the best R2, sessions per lag, S1TJ and M1TJ:
+25 ms    7
+0 ms     5
-25 ms    2
+50 ms    2
-50 ms    1

# Shared by the sections below: the other per-session analyses, each with its own seed.
t0 = time.time()
bt_sessions = [s for s in sessions if s["area"] in FIG3 and all(
    sum(1 for k in s["trials"]["posIndex"] if sl.BACKTRACK.get(k) == d) >= 10 for d in ("LR", "RL"))]
bt = sl.run_parallel(functools.partial(sl.backtrack_session, seed=5), bt_sessions)
zz = sl.run_parallel(functools.partial(sl.zigzag_session, seed=6), [s for s in sessions if s["area"] == "ZZ"])
iti = sl.run_parallel(functools.partial(sl.iti_session, seed=7), fig3_sessions)
pe = sl.run_parallel(functools.partial(sl.peth_session, seed=8), fig3_sessions)
print(f"SVM, ITI and PETH analyses in {time.time() - t0:.0f} s")
quality = pd.DataFrame([{"path": t["path"], "area": t["area"], "kind": t["kind"], "subject": t["subject"],
                         "rpv": q[0], "contam": q[1], "contam_censored": q[2], "contam_span": q[3]}
                        for t in tables for q in t["unit_quality"]])
quality["single"] = (quality.rpv <= sl.RPV_MAX) & (quality.contam <= sl.CONTAM_MAX)


def perm_paired(a, b):
    """Exact two-tailed sign-flip permutation p of the mean paired difference."""
    d = np.asarray(a) - np.asarray(b)
    signs = np.array(np.meshgrid(*[[1, -1]] * len(d))).reshape(len(d), -1)
    null = np.abs((d[:, None] * signs).mean(0))
    return float(np.mean(null >= abs(d.mean()) - 1e-12))

SVM, ITI and PETH analyses in 471 s

e57 = quality[(quality.kind == "ephys") & quality.area.isin(EPHYS57)]
n_sessions = sum(1 for t in tables if t["kind"] == "ephys" and t["area"] in EPHYS57)
print(f"Recording sessions (ephys ALM, M1TJ, S1TJ, S1BF, M1B, S1L, ZZ): {n_sessions}; units {len(e57)}")
print(f"Single units {int(e57.single.sum())}, multi-units {int((~e57.single).sum())} "
      "(paper: 1537 and 303, total 1840)")
print("\nBy area:")
print(e57.groupby("area").single.agg(units="size", single="sum").assign(multi=lambda d: d.units - d.single).to_string())
print("\nBy subject:")
print(e57.groupby("subject").single.agg(units="size", single="sum").to_string())
print("\nOther classifications of the same units:")
for name, m in [("RPV <= 1% only", e57.rpv <= 0.01), ("contamination <= 15% only", e57.contam <= 0.15),
                ("RPV <= 0.5% only (the stricter rule the Methods mention)", e57.rpv <= 0.005),
                ("RPV <= 1% and contamination <= 10%", (e57.rpv <= 0.01) & (e57.contam <= 0.10)),
                ("RPV <= 1% and contamination <= 20%", (e57.rpv <= 0.01) & (e57.contam <= 0.20))]:
    print(f"  {name}: single {int(m.sum())}, multi {int((~m).sum())}")
print("Other estimates of the contamination, with RPV <= 1% and contamination <= 15%:")
for name, col in [("0.5 ms censored period", "contam_censored"), ("spikes inside the trial span, T the whole span", "contam_span")]:
    m = (e57.rpv <= 0.01) & (e57[col] <= 0.15)
    print(f"  {name}: single {int(m.sum())}, multi {int((~m).sum())}")
other = quality[(quality.kind == "ephys") & ~quality.area.isin(EPHYS57)]
print(f"\nNot among the 57: 'ephys others' {len(other)} units in 22 sessions; opto efficiency "
      f"{int((quality.kind == 'opto VGAT-CRE Ai32').sum())} units in 14 sessions")

Recording sessions (ephys ALM, M1TJ, S1TJ, S1BF, M1B, S1L, ZZ): 57; units 1917
Single units 1578, multi-units 339 (paper: 1537 and 303, total 1840)

By area:
      units  single  multi
area                      
ALM     450     411     39
M1B     117     110      7
M1TJ    328     277     51
S1BF     78      70      8
S1L     453     333    120
S1TJ    265     165    100
ZZ      226     212     14

By subject:
          units  single
subject                
MX170903     58      57
MX180201     95      78
MX180202     62      48
MX180203    117     111
MX180601    317     225
MX180602     91      74
MX180701     57      56
MX180803     23      20
MX180804    159     119
MX181002    158     124
MX181101    234     182
MX181301    104      90
MX181302     45      44
MX190103    171     138
MX200101     87      78
MX200105     73      71
MX210301     66      63

Other classifications of the same units:
  RPV <= 1% only: single 1732, multi 185
  contamination <= 15% only: single 1586, multi 331
  RPV <= 0.5% only (the stricter rule the Methods mention): single 1608, multi 309
  RPV <= 1% and contamination <= 10%: single 1490, multi 427
  RPV <= 1% and contamination <= 20%: single 1622, multi 295
Other estimates of the contamination, with RPV <= 1% and contamination <= 15%:
  0.5 ms censored period: single 1576, multi 341
  spikes inside the trial span, T the whole span: single 1439, multi 478

Not among the 57: 'ephys others' 742 units in 22 sessions; opto efficiency 414 units in 14 sessions

# The split per session (Extended Data Fig. 4d shows units per session).
per = e57.groupby(["area", "path"]).single.agg(units="size", single="sum")
per["multi"] = per.units - per.single
per.index = [f"{a:5s} {p.split('/')[1][:32]}" for a, p in per.index]
print(per.to_string())
print(f"\nSessions: {len(per)}; multi-units per session median {per.multi.median():.0f}, range {per.multi.min()} to {per.multi.max()}")

                                        units  single  multi
ALM   sub-MX170903_ses-20180304T135947     58      57      1
ALM   sub-MX180202_ses-20180504T134624     37      25     12
ALM   sub-MX180202_ses-20180505T144012     25      23      2
ALM   sub-MX180203_ses-20180526T163118     14      12      2
ALM   sub-MX180203_ses-20180527T135221     27      24      3
ALM   sub-MX180203_ses-20180529T131703      9       9      0
ALM   sub-MX180701_ses-20180916T121311     57      56      1
ALM   sub-MX180804_ses-20181125T151623     33      30      3
ALM   sub-MX181002_ses-20181223T145116     16      15      1
ALM   sub-MX181002_ses-20181224T154648     44      42      2
ALM   sub-MX181101_ses-20190330T112923     48      48      0
ALM   sub-MX181101_ses-20190331T103702     39      29     10
ALM   sub-MX190103_ses-20190515T110643     43      41      2
M1B   sub-MX180201_ses-20180403T143829     50      44      6
M1B   sub-MX180203_ses-20180601T162617     27      26      1
M1B   sub-MX180203_ses-20180602T153218     30      30      0
M1B   sub-MX180203_ses-20180603T173845     10      10      0
M1TJ  sub-MX181101_ses-20190401T114039     43      32     11
M1TJ  sub-MX181101_ses-20190402T123529     39      32      7
M1TJ  sub-MX181301_ses-20190521T122250     50      41      9
M1TJ  sub-MX181301_ses-20190522T105041     30      26      4
M1TJ  sub-MX181301_ses-20190524T135431     24      23      1
M1TJ  sub-MX181302_ses-20190210T163258     32      32      0
M1TJ  sub-MX181302_ses-20190212T154624     13      12      1
M1TJ  sub-MX190103_ses-20190512T114257     49      41      8
M1TJ  sub-MX190103_ses-20190513T110307     48      38     10
S1BF  sub-MX180803_ses-20181128T130139     23      20      3
S1BF  sub-MX180804_ses-20181128T151110     17      14      3
S1BF  sub-MX180804_ses-20181129T150336     17      17      0
S1BF  sub-MX181002_ses-20181218T140101     12      11      1
S1BF  sub-MX181002_ses-20181219T132039      9       8      1
S1L   sub-MX180201_ses-20180405T154810     45      34     11
S1L   sub-MX180601_ses-20180720T123927     34      21     13
S1L   sub-MX180601_ses-20180721T130219     37      20     17
S1L   sub-MX180601_ses-20180722T130108     40      31      9
S1L   sub-MX180601_ses-20180723T124346     26      20      6
S1L   sub-MX180601_ses-20180725T123219     56      35     21
S1L   sub-MX180601_ses-20180726T131239     37      30      7
S1L   sub-MX180601_ses-20180727T130005     49      37     12
S1L   sub-MX180601_ses-20180728T140319     38      31      7
S1L   sub-MX180602_ses-20180707T103113     21      16      5
S1L   sub-MX180602_ses-20180713T123246     36      31      5
S1L   sub-MX180602_ses-20180714T100456     34      27      7
S1TJ  sub-MX180804_ses-20181201T130007     59      39     20
S1TJ  sub-MX180804_ses-20181202T124714     33      19     14
S1TJ  sub-MX181002_ses-20181228T150236     16       9      7
S1TJ  sub-MX181002_ses-20181229T151032     29      18     11
S1TJ  sub-MX181002_ses-20181230T141519     32      21     11
S1TJ  sub-MX181101_ses-20190328T114227     28      21      7
S1TJ  sub-MX181101_ses-20190329T105021     37      20     17
S1TJ  sub-MX190103_ses-20190514T112816     31      18     13
ZZ    sub-MX200101_ses-20201126T155559     40      35      5
ZZ    sub-MX200101_ses-20201127T154735     47      43      4
ZZ    sub-MX200105_ses-20201210T175011     33      33      0
ZZ    sub-MX200105_ses-20201212T152521     40      38      2
ZZ    sub-MX210301_ses-20210501T160243     27      25      2
ZZ    sub-MX210301_ses-20210506T151300     39      38      1

Sessions: 57; multi-units per session median 5, range 0 to 21

# The contamination estimate on synthetic units of known contamination (20 units per level,
# 10 spikes/s over 1 h, refractory 2.5 ms; contaminating spikes are Poisson).
rng = np.random.default_rng(21)
syn = []
for c in (0.0, 0.05, 0.10, 0.15, 0.20, 0.30):
    for _ in range(20):
        u = sl.simulate_unit(10.0, c, 3600.0, rng)
        q = sl.unit_quality(u, 0.0, 3600.0)
        syn.append({"true": c, "rpv": q[0], "estimate": q[1], "censored_0.5ms": q[2]})
syn = pd.DataFrame(syn)
print(syn.groupby("true")[["rpv", "estimate", "censored_0.5ms"]].agg(["mean", "std"]).round(3).to_string())
print("Fraction classified multi (RPV > 1% or contamination > 15%) by true contamination:")
print(syn.assign(multi=(syn.rpv > 0.01) | (syn.estimate > 0.15)).groupby("true").multi.mean().round(2).to_string())

        rpv        estimate        censored_0.5ms       
       mean    std     mean    std           mean    std
true                                                    
0.00  0.000  0.000    0.000  0.000          0.000  0.000
0.05  0.002  0.000    0.050  0.005          0.050  0.006
0.10  0.005  0.000    0.104  0.009          0.103  0.010
0.15  0.007  0.000    0.163  0.014          0.162  0.017
0.20  0.009  0.000    0.233  0.016          0.236  0.020
0.30  0.013  0.001    0.762  0.300          0.764  0.298
Fraction classified multi (RPV > 1% or contamination > 15%) by true contamination:
true
0.00    0.00
0.05    0.00
0.10    0.00
0.15    0.85
0.20    1.00
0.30    1.00

cnt = pd.DataFrame([{"area": o["area"], "subject": o["subject"], "units": o["n_units"],
                     "LR": o["n_trials"]["LR"], "RL": o["n_trials"]["RL"]} for o in f3])
print(cnt.groupby("area").agg(sessions=("area", "size"), mice=("subject", "nunique"),
                              min_LR=("LR", "min"), min_RL=("RL", "min"), min_units=("units", "min")).loc[FIG3].to_string())
print("\nEvery session has a units table and selected standard trials in both directions:",
      bool((cnt.units > 0).all() and (cnt.LR >= 10).all() and (cnt.RL >= 10).all()))

      sessions  mice  min_LR  min_RL  min_units
area                                           
S1TJ         8     4      23      23         16
M1TJ         9     4      16      18         13
ALM         13     8      21      13          9
S1BF         5     3      31      24          9

Every session has a units table and selected standard trials in both directions: True

u35 = cnt.units
u57 = pd.Series([t["n_units"] for t in tables if t["kind"] == "ephys" and t["area"] in EPHYS57])
print(f"(a) 35 Figure 3 sessions: mean {u35.mean():.1f}, SD {u35.std():.1f} units (range {u35.min()} to {u35.max()})")
print(f"(b) 57 recording sessions: mean {u57.mean():.1f}, SD {u57.std():.1f} units")
print("\nPer area, Figure 3 sessions:")
print(cnt.groupby("area").units.agg(["mean", "std", "min", "max"]).round(1).loc[FIG3].to_string())
sq = quality[quality.path.isin([o["path"] for o in f3])].groupby("path").single.sum()
print(f"\nSingle units only, 35 sessions: mean {sq.mean():.1f}, SD {sq.std():.1f}")

(a) 35 Figure 3 sessions: mean 32.0, SD 14.4 units (range 9 to 59)
(b) 57 recording sessions: mean 33.6, SD 13.0 units

Per area, Figure 3 sessions:
      mean   std  min  max
area                      
S1TJ  33.1  12.1   16   59
M1TJ  36.4  12.7   13   50
ALM   34.6  15.8    9   58
S1BF  15.6   5.4    9   23

Single units only, 35 sessions: mean 26.4, SD 13.2

nstd = sum(sum(o["n_trials"].values()) for o in f3)
nbt = sum(sum(o["n_bt"].values()) for o in f3)
print(f"Selected trials over the 35 sessions: standard {nstd}, backtracking {nbt}")
print(f"(a) standard and backtracking: {nstd + nbt}; (b) standard only: {nstd}  (paper: 2684)")
print("\nBy area (standard, backtracking):")
print(pd.DataFrame([{"area": o["area"], "standard": sum(o["n_trials"].values()), "backtracking": sum(o["n_bt"].values())}
                    for o in f3]).groupby("area").sum().loc[FIG3].to_string())

Selected trials over the 35 sessions: standard 2169, backtracking 717
(a) standard and backtracking: 2886; (b) standard only: 2169  (paper: 2684)

By area (standard, backtracking):
      standard  backtracking
area                        
S1TJ       518           205
M1TJ       443           183
ALM        860           164
S1BF       348           165

# Does the greedy search find the most similar subset? (1) For every Figure 3 session and direction,
# exhaustive search over the first 16 candidate trials for the best 6, against greedy on the same
# 16. (2) On the full candidate sets, greedy against a local search that swaps one chosen and one
# unchosen trial while the mean pairwise distance falls.
import itertools


def local_search(D, chosen):
    chosen = list(chosen)
    rest = [k for k in range(len(D)) if k not in chosen]
    best = sl.subset_cost(D, chosen)
    improved = True
    while improved:
        improved = False
        for a in range(len(chosen)):
            for b in range(len(rest)):
                cand = chosen[:a] + [rest[b]] + chosen[a + 1:]
                c = sl.subset_cost(D, cand)
                if c < best - 1e-12:
                    chosen[a], rest[b] = rest[b], chosen[a]
                    best, improved = c, True
    return sorted(chosen), best


gx = []
for s in fig3_sessions:
    for d in ("LR", "RL"):
        D = sl.selection_features(s, "standard", d)
        if len(D) >= 16:
            Ds = D[:16, :16]
            combos = np.array(list(itertools.combinations(range(16), 6)))
            costs = Ds[combos[:, :, None], combos[:, None, :]].sum((1, 2)) / 30
            g = sl.greedy_subset(Ds, 6)
            ex = (sl.subset_cost(Ds, g), costs.min(), float(np.mean(costs <= sl.subset_cost(Ds, g) + 1e-12)))
        else:
            ex = (np.nan, np.nan, np.nan)
        n = max(10, int(round(len(D) / 3)))
        if len(D) <= 10:
            continue
        g = sl.greedy_subset(D, n)
        ls, lc = local_search(D, g)
        gx.append({"area": s["area"], "dir": d, "candidates": len(D), "chosen": n, "greedy16": ex[0], "exhaustive16": ex[1],
                   "rank16": ex[2], "greedy": sl.subset_cost(D, g), "local": lc, "shared": len(set(g) & set(ls)) / n})
gx = pd.DataFrame(gx)
e16 = gx.dropna(subset=["greedy16"])
print(f"(1) Exhaustive, 16 candidates choose 6: {len(e16)} session-directions; greedy finds the best subset in "
      f"{int(np.isclose(e16.greedy16, e16.exhaustive16).sum())}; greedy cost over best: median "
      f"{(e16.greedy16 / e16.exhaustive16).median():.3f}, max {(e16.greedy16 / e16.exhaustive16).max():.3f}; "
      f"fraction of the 8008 subsets at least as good as greedy: median {e16.rank16.median():.5f}, max {e16.rank16.max():.5f}")
print(f"(2) Full sets: {len(gx)} session-directions; local search improves on greedy in {int((gx.local < gx.greedy - 1e-12).sum())}; "
      f"greedy cost over local optimum median {(gx.greedy / gx.local).median():.3f}, max {(gx.greedy / gx.local).max():.3f}; "
      f"trials shared median {gx.shared.median():.2f}, min {gx.shared.min():.2f}")

(1) Exhaustive, 16 candidates choose 6: 70 session-directions; greedy finds the best subset in 43; greedy cost over best: median 1.000, max 1.084; fraction of the 8008 subsets at least as good as greedy: median 0.00012, max 0.00400
(2) Full sets: 70 session-directions; local search improves on greedy in 37; greedy cost over local optimum median 1.000, max 1.049; trials shared median 0.97, min 0.53

r2 = pd.DataFrame([{"area": o["area"], "subject": o["subject"], "session": o["path"].split("/")[1][:32],
                    **{v: o["r2"][v] for v in sl.VARS}, **{v + "_trialfolds": o["r2_trialfolds"][v] for v in sl.VARS}}
                   for o in f3])
PAPER_R2 = {"L": (0.56, 0.42, 0.11, 0.04), "Ldot": (0.35, 0.32, 0.07, 0.01), "theta": (0.59, 0.57, 0.54, 0.12),
            "I": (0.23, 0.42, 0.67, 0.14), "tau": (0.17, 0.45, 0.60, 0.08)}
PAPER_SD = {"L": (0.10, 0.19, 0.09, 0.03), "Ldot": (0.12, 0.10, 0.06, 0.02), "theta": (0.12, 0.11, 0.12, 0.10),
            "I": (0.13, 0.15, 0.16, 0.05), "tau": (0.10, 0.15, 0.10, 0.03)}


def show_r2(v):
    g = r2.groupby("area")
    tab = pd.DataFrame({"mean": g[v].mean(), "sd": g[v].std(), "min": g[v].min(), "max": g[v].max(),
                        "trialfolds_mean": g[v + "_trialfolds"].mean(), "sessions": g[v].size()}).loc[FIG3]
    tab["paper"] = PAPER_R2[v]
    tab["paper_sd"] = PAPER_SD[v]
    tab["within"] = (tab["mean"] - tab["paper"]).abs() <= tab["paper_sd"]
    print(f"Cross-validated R2 of {v} by area (paper: mean and SD read from Figure 3b):")
    print(tab.round(3).to_string())
    print("\nBy subject (mean over sessions):")
    print(r2.groupby(["area", "subject"])[v].agg(["mean", "size"]).round(3).to_string())
    fig, ax = plt.subplots(figsize=(4, 2.5))
    for k, a in enumerate(FIG3):
        y = r2[r2.area == a][v]
        ax.scatter(np.full(len(y), k), y, s=12)
        ax.bar(k, y.mean(), alpha=0.3)
    ax.set_xticks(range(4), FIG3)
    ax.set_ylabel(f"R2 of {v}")
    ax.set_ylim(-0.05, 1)
    plt.show()
    print(f"Points outside the axes (below -0.05 or above 1): {int(((r2[v] < -0.05) | (r2[v] > 1)).sum())}")
    return tab

tab = show_r2("L")
print("\nOrder (min of S1TJ, M1TJ > max of ALM, S1BF):", min(tab.loc[["S1TJ", "M1TJ"], "mean"]) > max(tab.loc[["ALM", "S1BF"], "mean"]),
      "; trial-fold variant:", min(tab.loc[["S1TJ", "M1TJ"], "trialfolds_mean"]) > max(tab.loc[["ALM", "S1BF"], "trialfolds_mean"]))

Cross-validated R2 of L by area (paper: mean and SD read from Figure 3b):
       mean     sd    min    max  trialfolds_mean  sessions  paper  paper_sd  within
area                                                                                
S1TJ  0.568  0.106  0.475  0.743            0.551         8   0.56      0.10    True
M1TJ  0.423  0.175  0.125  0.683            0.392         9   0.42      0.19    True
ALM   0.117  0.089  0.005  0.288            0.087        13   0.11      0.09    True
S1BF  0.032  0.020  0.006  0.061            0.016         5   0.04      0.03    True

By subject (mean over sessions):
                mean  size
area subject              
ALM  MX170903  0.065     1
     MX180202  0.157     2
     MX180203  0.049     3
     MX180701  0.288     1
     MX180804  0.133     1
     MX181002  0.149     2
     MX181101  0.056     2
     MX190103  0.164     1
M1TJ MX181101  0.585     2
     MX181301  0.328     3
     MX181302  0.351     2
     MX190103  0.475     2
S1BF MX180803  0.028     1
     MX180804  0.032     2
     MX181002  0.034     2
S1TJ MX180804  0.724     2
     MX181002  0.493     3
     MX181101  0.501     2
     MX190103  0.618     1

Points outside the axes (below -0.05 or above 1): 0

Order (min of S1TJ, M1TJ > max of ALM, S1BF): True ; trial-fold variant: True

tab = show_r2("Ldot")
print("\nOrder (min of S1TJ, M1TJ > max of ALM, S1BF):", min(tab.loc[["S1TJ", "M1TJ"], "mean"]) > max(tab.loc[["ALM", "S1BF"], "mean"]),
      "; trial-fold variant:", min(tab.loc[["S1TJ", "M1TJ"], "trialfolds_mean"]) > max(tab.loc[["ALM", "S1BF"], "trialfolds_mean"]))

Cross-validated R2 of Ldot by area (paper: mean and SD read from Figure 3b):
       mean     sd    min    max  trialfolds_mean  sessions  paper  paper_sd  within
area                                                                                
S1TJ  0.361  0.111  0.245  0.537            0.344         8   0.35      0.12    True
M1TJ  0.326  0.102  0.112  0.454            0.305         9   0.32      0.10    True
ALM   0.081  0.060 -0.000  0.203            0.064        13   0.07      0.06    True
S1BF  0.019  0.018  0.005  0.048            0.014         5   0.01      0.02    True

By subject (mean over sessions):
                mean  size
area subject              
ALM  MX170903  0.026     1
     MX180202  0.123     2
     MX180203  0.041     3
     MX180701  0.203     1
     MX180804  0.062     1
     MX181002  0.111     2
     MX181101  0.036     2
     MX190103  0.106     1
M1TJ MX181101  0.408     2
     MX181301  0.238     3
     MX181302  0.331     2
     MX190103  0.371     2
S1BF MX180803  0.048     1
     MX180804  0.015     2
     MX181002  0.007     2
S1TJ MX180804  0.399     2
     MX181002  0.303     3
     MX181101  0.383     2
     MX190103  0.419     1

Points outside the axes (below -0.05 or above 1): 0

Order (min of S1TJ, M1TJ > max of ALM, S1BF): True ; trial-fold variant: True

tab = show_r2("theta")
three = tab.loc[["S1TJ", "M1TJ", "ALM"], "mean"]
print(f"\nS1BF below each of the others: {bool((tab.loc['S1BF', 'mean'] < three).all())}; "
      f"spread of the other three (max - min): {three.max() - three.min():.3f}")
t3 = tab.loc[["S1TJ", "M1TJ", "ALM"], "trialfolds_mean"]
print(f"Trial-fold variant: spread {t3.max() - t3.min():.3f}, S1BF {tab.loc['S1BF', 'trialfolds_mean']:.3f}")

Cross-validated R2 of theta by area (paper: mean and SD read from Figure 3b):
       mean     sd    min    max  trialfolds_mean  sessions  paper  paper_sd  within
area                                                                                
S1TJ  0.569  0.074  0.415  0.640            0.528         8   0.59      0.12    True
M1TJ  0.534  0.107  0.345  0.630            0.476         9   0.57      0.11    True
ALM   0.520  0.130  0.289  0.709            0.469        13   0.54      0.12    True
S1BF  0.111  0.107  0.015  0.235            0.071         5   0.12      0.10    True

By subject (mean over sessions):
                mean  size
area subject              
ALM  MX170903  0.533     1
     MX180202  0.477     2
     MX180203  0.485     3
     MX180701  0.699     1
     MX180804  0.659     1
     MX181002  0.426     2
     MX181101  0.544     2
     MX190103  0.523     1
M1TJ MX181101  0.498     2
     MX181301  0.506     3
     MX181302  0.564     2
     MX190103  0.580     2
S1BF MX180803  0.219     1
     MX180804  0.147     2
     MX181002  0.021     2
S1TJ MX180804  0.517     2
     MX181002  0.575     3
     MX181101  0.577     2
     MX190103  0.640     1

Points outside the axes (below -0.05 or above 1): 0

S1BF below each of the others: True; spread of the other three (max - min): 0.049
Trial-fold variant: spread 0.059, S1BF 0.071

tab = show_r2("I")
m = tab["mean"]
print("\nOrder ALM > M1TJ > S1TJ > S1BF:", bool(m["ALM"] > m["M1TJ"] > m["S1TJ"] > m["S1BF"]))
mt = tab["trialfolds_mean"]
print("Trial-fold variant order:", bool(mt["ALM"] > mt["M1TJ"] > mt["S1TJ"] > mt["S1BF"]))

Cross-validated R2 of I by area (paper: mean and SD read from Figure 3b):
       mean     sd    min    max  trialfolds_mean  sessions  paper  paper_sd  within
area                                                                                
S1TJ  0.270  0.141  0.126  0.601            0.188         8   0.23      0.13    True
M1TJ  0.412  0.167  0.126  0.606            0.323         9   0.42      0.15    True
ALM   0.693  0.142  0.376  0.892            0.657        13   0.67      0.16    True
S1BF  0.113  0.047  0.069  0.180            0.049         5   0.14      0.05    True

By subject (mean over sessions):
                mean  size
area subject              
ALM  MX170903  0.775     1
     MX180202  0.583     2
     MX180203  0.635     3
     MX180701  0.850     1
     MX180804  0.764     1
     MX181002  0.649     2
     MX181101  0.706     2
     MX190103  0.834     1
M1TJ MX181101  0.525     2
     MX181301  0.421     3
     MX181302  0.266     2
     MX190103  0.431     2
S1BF MX180803  0.180     1
     MX180804  0.103     2
     MX181002  0.089     2
S1TJ MX180804  0.420     2
     MX181002  0.254     3
     MX181101  0.166     2
     MX190103  0.225     1

Points outside the axes (below -0.05 or above 1): 0

Order ALM > M1TJ > S1TJ > S1BF: True
Trial-fold variant order: True

tab = show_r2("tau")
m = tab["mean"]
print("\nOrder ALM > M1TJ > S1TJ > S1BF:", bool(m["ALM"] > m["M1TJ"] > m["S1TJ"] > m["S1BF"]))
mt = tab["trialfolds_mean"]
print("Trial-fold variant order:", bool(mt["ALM"] > mt["M1TJ"] > mt["S1TJ"] > mt["S1BF"]))

Cross-validated R2 of tau by area (paper: mean and SD read from Figure 3b):
       mean     sd    min    max  trialfolds_mean  sessions  paper  paper_sd  within
area                                                                                
S1TJ  0.198  0.120  0.068  0.463            0.145         8   0.17      0.10    True
M1TJ  0.440  0.140  0.176  0.650            0.381         9   0.45      0.15    True
ALM   0.637  0.083  0.529  0.807            0.601        13   0.60      0.10    True
S1BF  0.087  0.016  0.068  0.105            0.063         5   0.08      0.03    True

By subject (mean over sessions):
                mean  size
area subject              
ALM  MX170903  0.544     1
     MX180202  0.612     2
     MX180203  0.630     3
     MX180701  0.709     1
     MX180804  0.574     1
     MX181002  0.677     2
     MX181101  0.676     2
     MX190103  0.635     1
M1TJ MX181101  0.553     2
     MX181301  0.422     3
     MX181302  0.363     2
     MX190103  0.430     2
S1BF MX180803  0.088     1
     MX180804  0.103     2
     MX181002  0.071     2
S1TJ MX180804  0.325     2
     MX181002  0.125     3
     MX181101  0.203     2
     MX190103  0.155     1

Points outside the axes (below -0.05 or above 1): 0

Order ALM > M1TJ > S1TJ > S1BF: True
Trial-fold variant order: True

cr = pd.DataFrame([o["corr"] for o in f3])
summ = pd.DataFrame({"mean_abs_r": cr.mean(), "max_over_sessions": cr.max()})
summ.index = [f"{a}-{b}" for a, b in summ.index]
print("Absolute Pearson r per pair, over 35 sessions:")
print(summ.round(3).to_string())
print(f"\nPairs with mean |r| >= 0.3: {int((summ.mean_abs_r >= 0.3).sum())} of {len(summ)}")

Absolute Pearson r per pair, over 35 sessions:
            mean_abs_r  max_over_sessions
L-Ldot           0.074              0.259
L-theta          0.103              0.368
L-I              0.042              0.117
L-tau            0.062              0.173
Ldot-theta       0.078              0.229
Ldot-I           0.013              0.033
Ldot-tau         0.033              0.095
theta-I          0.118              0.257
theta-tau        0.083              0.348
I-tau            0.000              0.000

Pairs with mean |r| >= 0.3: 0 of 10

ve = pd.DataFrame([{"area": o["area"], "axes_sum": sum(o["ve_axes"].values()), "pc5": o["ve_pc5"],
                    **{"ve_" + v: o["ve_axes"][v] for v in sl.VARS}} for o in f3])
ve["ratio"] = ve.axes_sum / ve.pc5
print(ve.groupby("area")[["axes_sum", "pc5", "ratio"]].median().round(3).loc[FIG3].to_string())
print(f"\nAll 35 sessions: median ratio {ve.ratio.median():.3f} (IQR {ve.ratio.quantile(0.25):.3f} to "
      f"{ve.ratio.quantile(0.75):.3f}); mean {ve.ratio.mean():.3f}")
print(f"Coding axes with all weights zero: {sum(1 for o in f3 for v in sl.VARS if not np.any(o['weights'][v]))} of {5 * len(f3)}")

      axes_sum    pc5  ratio
area                        
S1TJ     0.202  0.459  0.442
M1TJ     0.189  0.452  0.476
ALM      0.187  0.462  0.454
S1BF     0.228  0.565  0.433

All 35 sessions: median ratio 0.454 (IQR 0.387 to 0.544); mean 0.490
Coding axes with all weights zero: 1 of 175

cs = pd.DataFrame([o["cos"] for o in f3])
summ = pd.DataFrame({"mean_abs_cos": cs.mean(), "sessions": cs.count()})
summ.index = [f"{a}-{b}" for a, b in summ.index]
print(summ.round(3).to_string())
print(f"\nPairs with mean |cos| >= 0.3: {int((summ.mean_abs_cos >= 0.3).sum())} of {len(summ)}; "
      f"random directions in the same sessions: mean |cos| {np.mean([o['cos_random'] for o in f3]):.3f}")

            mean_abs_cos  sessions
L-Ldot             0.225        34
L-theta            0.151        35
L-I                0.151        35
L-tau              0.182        35
Ldot-theta         0.302        34
Ldot-I             0.114        34
Ldot-tau           0.223        34
theta-I            0.215        35
theta-tau          0.212        35
I-tau              0.209        35

Pairs with mean |cos| >= 0.3: 1 of 10; random directions in the same sessions: mean |cos| 0.154

t0 = time.time()
SEEDS = (11, 12, 13)
f3s = {3: f3}
for sd in SEEDS:
    f3s[sd] = sl.run_parallel(functools.partial(sl.fig3_session, seed=sd, trialfolds=False), fig3_sessions)
print(f"Fits with {len(SEEDS)} more seeds in {time.time() - t0:.0f} s")
srows = []
for sd, res in f3s.items():
    ratio = pd.Series([sum(o["ve_axes"].values()) / o["ve_pc5"] for o in res])
    cosm = pd.DataFrame([o["cos"] for o in res]).mean()
    r2m = pd.DataFrame([{"area": o["area"], **o["r2"]} for o in res]).groupby("area").mean()
    srows.append({"seed": sd, "ed5e_median_ratio": ratio.median(), "ed5f_pairs_ge_0.3": int((cosm >= 0.3).sum()),
                  "ed5f_largest_pair": f"{'-'.join(cosm.idxmax())} {cosm.max():.3f}",
                  "Ldot_theta_cos": cosm[("Ldot", "theta")],
                  **{f"{v}_{a}": r2m.loc[a, v] for v in sl.VARS for a in FIG3}})
sv = pd.DataFrame(srows).set_index("seed")
print(sv[["ed5e_median_ratio", "ed5f_pairs_ge_0.3", "ed5f_largest_pair", "Ldot_theta_cos"]].round(3).to_string())
print(f"\ned5e median ratio over seeds: {sv.ed5e_median_ratio.min():.3f} to {sv.ed5e_median_ratio.max():.3f}; "
      f"seeds inside 0.45 to 0.55: {int(sv.ed5e_median_ratio.between(0.45, 0.55).sum())} of {len(sv)}")
print(f"ed5f: seeds with no pair >= 0.3: {int((sv['ed5f_pairs_ge_0.3'] == 0).sum())}, with one pair: "
      f"{int((sv['ed5f_pairs_ge_0.3'] == 1).sum())}, with two or more: {int((sv['ed5f_pairs_ge_0.3'] >= 2).sum())} of {len(sv)}")
print("\nMean R2 by area and variable, range over the 4 seeds:")
print(pd.DataFrame({v: [f"{sv[f'{v}_{a}'].min():.3f} to {sv[f'{v}_{a}'].max():.3f}" for a in FIG3] for v in sl.VARS},
                   index=FIG3).to_string())
print("\nUnits per session do not depend on the seed: SD over the 35 sessions",
      round(float(np.std([o["n_units"] for o in f3s[11]], ddof=1)), 1), "with seed 11, as with seed 3.")

Fits with 3 more seeds in 367 s
      ed5e_median_ratio  ed5f_pairs_ge_0.3 ed5f_largest_pair  Ldot_theta_cos
seed                                                                        
3                 0.454                  1  Ldot-theta 0.302           0.302
11                0.455                  1  Ldot-theta 0.320           0.320
12                0.455                  1  Ldot-theta 0.306           0.306
13                0.453                  1  Ldot-theta 0.317           0.317

ed5e median ratio over seeds: 0.453 to 0.455; seeds inside 0.45 to 0.55: 4 of 4
ed5f: seeds with no pair >= 0.3: 0, with one pair: 4, with two or more: 0 of 4

Mean R2 by area and variable, range over the 4 seeds:
                   L            Ldot           theta               I             tau
S1TJ  0.568 to 0.568  0.361 to 0.362  0.566 to 0.569  0.270 to 0.270  0.197 to 0.198
M1TJ  0.422 to 0.423  0.325 to 0.328  0.533 to 0.534  0.412 to 0.412  0.438 to 0.440
ALM   0.117 to 0.118  0.080 to 0.083  0.520 to 0.521  0.693 to 0.693  0.637 to 0.637
S1BF  0.032 to 0.032  0.013 to 0.019  0.107 to 0.111  0.113 to 0.114  0.086 to 0.087

Units per session do not depend on the seed: SD over the 35 sessions 14.4 with seed 11, as with seed 3.

cc = pd.DataFrame([{"area": o["area"], "g1": o["cca"][1], "g2": o["cca"][2], "g3": o["cca"][3]} for o in f3])
print("Mean canonical correlation by group (mean over sessions):")
print(cc.groupby("area")[["g1", "g2", "g3"]].mean().round(3).loc[FIG3].to_string())
tests = [("ALM", "g2", "g1"), ("ALM", "g2", "g3"), ("S1TJ", "g1", "g2"), ("S1TJ", "g3", "g2"), ("M1TJ", "g3", "g2")]
held = 0
print("\nComparisons the figure marks (paired two-tailed sign-flip permutation, exact):")
for a, hi, lo in tests:
    d = cc[cc.area == a]
    p = perm_paired(d[hi], d[lo])
    ok = d[hi].mean() > d[lo].mean() and p < 0.05
    held += ok
    print(f"  {a}: {hi} {d[hi].mean():.3f} vs {lo} {d[lo].mean():.3f}, difference {d[hi].mean() - d[lo].mean():+.3f}, "
          f"p = {p:.4f}, n = {len(d)} -> {'holds' if ok else 'does not hold'}")
print(f"{held} of 5 comparisons hold")

Mean canonical correlation by group (mean over sessions):
         g1     g2     g3
area                     
S1TJ  0.746  0.589  0.735
M1TJ  0.712  0.662  0.809
ALM   0.607  0.900  0.743
S1BF  0.661  0.605  0.690

Comparisons the figure marks (paired two-tailed sign-flip permutation, exact):
  ALM: g2 0.900 vs g1 0.607, difference +0.293, p = 0.0002, n = 13 -> holds
  ALM: g2 0.900 vs g3 0.743, difference +0.157, p = 0.0002, n = 13 -> holds
  S1TJ: g1 0.746 vs g2 0.589, difference +0.156, p = 0.0234, n = 8 -> holds
  S1TJ: g3 0.735 vs g2 0.589, difference +0.146, p = 0.0078, n = 8 -> holds
  M1TJ: g3 0.809 vs g2 0.662, difference +0.147, p = 0.0195, n = 9 -> holds
5 of 5 comparisons hold

print("Sessions with backtracking in both directions, by area:", Counter(o["area"] for o in bt))
print("Trials per class per session (min, max): standard",
      (min(o["n"]["standard"] for o in bt), max(o["n"]["standard"] for o in bt)), "backtracking",
      (min(o["n"]["backtrack"] for o in bt), max(o["n"]["backtrack"] for o in bt)))
fig, axes = plt.subplots(1, 4, figsize=(12, 2.5), sharey=True)
rows = []
for ax, a in zip(axes, FIG3):
    os_ = [o for o in bt if o["area"] == a]
    acc = np.mean([o["acc"] for o in os_], 0)
    real = np.mean([o["acc"][o["null_bins"]].mean() for o in os_])
    null = np.mean([o["null"].mean(1) for o in os_], 0)
    rows.append({"area": a, "sessions": len(os_), "mice": len({o["subject"] for o in os_}),
                 "acc_0.2_0.4": real, "null_mean": null.mean(), "null_95": np.percentile(null, 95),
                 "decodes": real > np.percentile(null, 95),
                 "acc_before_0.1": np.mean([o["acc"][o["t"] < 0.1].mean() for o in os_])})
    ax.plot(os_[0]["t"], acc)
    ax.axhline(0.5, color="k", lw=0.5)
    ax.set_title(a)
    ax.set_xlabel("time from backtracking trigger (s)")
axes[0].set_ylabel("accuracy")
plt.show()
res = pd.DataFrame(rows)
print(res.round(3).to_string(index=False))
print(f"Plotted mean accuracies range from {min(np.mean([o['acc'] for o in bt if o['area'] == a], 0).min() for a in FIG3):.3f} "
      f"to {max(np.mean([o['acc'] for o in bt if o['area'] == a], 0).max() for a in FIG3):.3f}; none outside the axes.")
print("\nPer session, mean accuracy in 0.2 to 0.4 s and its own null 95th percentile:")
for o in bt:
    print(f"  {o['area']:5s} {o['path'].split('/')[1][:32]}  {o['acc'][o['null_bins']].mean():.3f}  "
          f"{np.percentile(o['null'].mean(1), 95):.3f}")

Sessions with backtracking in both directions, by area: Counter({'M1TJ': 9, 'S1TJ': 8, 'ALM': 6, 'S1BF': 5})
Trials per class per session (min, max): standard (32, 84) backtracking (20, 38)

area  sessions  mice  acc_0.2_0.4  null_mean  null_95  decodes  acc_before_0.1
S1TJ         8     4        0.595      0.500    0.516     True           0.487
M1TJ         9     4        0.734      0.504    0.523     True           0.513
 ALM         6     4        0.740      0.500    0.526     True           0.487
S1BF         5     3        0.556      0.499    0.515     True           0.519
Plotted mean accuracies range from 0.439 to 0.830; none outside the axes.

Per session, mean accuracy in 0.2 to 0.4 s and its own null 95th percentile:
  S1BF  sub-MX180803_ses-20181128T130139  0.564  0.527
  ALM   sub-MX180804_ses-20181125T151623  0.798  0.566
  S1BF  sub-MX180804_ses-20181128T151110  0.603  0.542
  S1TJ  sub-MX180804_ses-20181202T124714  0.653  0.547
  S1TJ  sub-MX180804_ses-20181201T130007  0.588  0.555
  S1BF  sub-MX180804_ses-20181129T150336  0.520  0.541
  ALM   sub-MX181002_ses-20181224T154648  0.680  0.554
  S1TJ  sub-MX181002_ses-20181229T151032  0.495  0.544
  S1TJ  sub-MX181002_ses-20181228T150236  0.541  0.526
  S1TJ  sub-MX181002_ses-20181230T141519  0.605  0.561
  S1BF  sub-MX181002_ses-20181219T132039  0.585  0.551
  S1BF  sub-MX181002_ses-20181218T140101  0.507  0.544
  ALM   sub-MX181002_ses-20181223T145116  0.680  0.550
  S1TJ  sub-MX181101_ses-20190329T105021  0.723  0.553
  M1TJ  sub-MX181101_ses-20190401T114039  0.945  0.592
  S1TJ  sub-MX181101_ses-20190328T114227  0.634  0.553
  ALM   sub-MX181101_ses-20190331T103702  0.705  0.550
  ALM   sub-MX181101_ses-20190330T112923  0.892  0.550
  M1TJ  sub-MX181101_ses-20190402T123529  0.839  0.574
  M1TJ  sub-MX181301_ses-20190522T105041  0.673  0.575
  M1TJ  sub-MX181301_ses-20190524T135431  0.676  0.533
  M1TJ  sub-MX181301_ses-20190521T122250  0.606  0.546
  M1TJ  sub-MX181302_ses-20190210T163258  0.776  0.552
  M1TJ  sub-MX181302_ses-20190212T154624  0.777  0.553
  M1TJ  sub-MX190103_ses-20190513T110307  0.669  0.585
  M1TJ  sub-MX190103_ses-20190512T114257  0.645  0.540
  ALM   sub-MX190103_ses-20190515T110643  0.684  0.560
  S1TJ  sub-MX190103_ses-20190514T112816  0.522  0.543

print(f"Zigzag sessions {len(zz)}, mice {len({o['subject'] for o in zz})}")
print("Perfectly executed trials per session (A starts at position 0, B at position 4):",
      [(o["subject"], o["n"]["A"], o["n"]["B"]) for o in zz])
zrows = []
fig, axes = plt.subplots(1, 4, figsize=(12, 2.5), sharey=True)
for k, ax in enumerate(axes):
    accs = [o["comps"][k]["acc"] for o in zz]
    real = np.mean([a.mean() for a in accs])
    null = np.mean([o["comps"][k]["null"].mean(1) for o in zz], 0)
    zrows.append({"comparison": zz[0]["comps"][k]["label"], "mean_acc": real, "null_mean": null.mean(),
                  "null_95": np.percentile(null, 95), "decodes": real > np.percentile(null, 95),
                  "lowest_session": min(a.mean() for a in accs)})
    ax.plot(zz[0]["comps"][k]["t"], np.mean(accs, 0))
    ax.axhline(0.5, color="k", lw=0.5)
    ax.set_title(zz[0]["comps"][k]["label"], fontsize=9)
    ax.set_xlabel("time from aligned touch (s)")
axes[0].set_ylabel("accuracy")
plt.show()
print(pd.DataFrame(zrows).round(3).to_string(index=False))
print("By mouse (mean over its sessions and the four comparisons):")
print(pd.DataFrame([{"subject": o["subject"], "acc": np.mean([c["acc"].mean() for c in o["comps"]])} for o in zz])
      .groupby("subject").acc.mean().round(3).to_string())

Zigzag sessions 6, mice 3
Perfectly executed trials per session (A starts at position 0, B at position 4): [('MX200101', 36, 44), ('MX200101', 18, 54), ('MX200105', 51, 62), ('MX200105', 42, 29), ('MX210301', 63, 72), ('MX210301', 79, 78)]

       comparison  mean_acc  null_mean  null_95  decodes  lowest_session
 fix B6, shift A2     0.819      0.504    0.518     True           0.739
fix B6, shift A10     0.809      0.500    0.511     True           0.654
 fix A6, shift B2     0.858      0.501    0.510     True           0.730
fix A6, shift B10     0.785      0.501    0.513     True           0.651
By mouse (mean over its sessions and the four comparisons):
subject
MX200101    0.807
MX200105    0.730
MX210301    0.916

it = pd.DataFrame(iti)
print(it.groupby("area")[["TP", "TP_trialfolds", "n_trials"]].agg(["mean", "std"]).round(3).loc[FIG3].to_string())
alm = it[it.area == "ALM"].TP
for a in ["M1TJ", "S1TJ", "S1BF"]:
    o = it[it.area == a].TP
    p = stats.mannwhitneyu(alm, o, alternative="greater").pvalue
    pt = stats.mannwhitneyu(it[it.area == "ALM"].TP_trialfolds, it[it.area == a].TP_trialfolds, alternative="greater").pvalue
    print(f"ALM > {a}: difference {alm.mean() - o.mean():+.3f}, Mann-Whitney p = {p:.4f} (trial-fold variant p = {pt:.4f})")
mouse = it.groupby(["area", "subject"]).TP.mean().reset_index()
print("\nVariant, one value per mouse (mean over its sessions):")
for a in ["M1TJ", "S1TJ", "S1BF"]:
    x, y = mouse[mouse.area == "ALM"].TP, mouse[mouse.area == a].TP
    print(f"  ALM ({len(x)} mice) > {a} ({len(y)} mice): difference {x.mean() - y.mean():+.3f}, "
          f"Mann-Whitney p = {stats.mannwhitneyu(x, y, alternative='greater').pvalue:.4f}")
print("\nBy subject:")
print(it.groupby(["area", "subject"]).TP.mean().round(3).to_string())

         TP        TP_trialfolds        n_trials        
       mean    std          mean    std     mean     std
area                                                    
S1TJ  0.075  0.036         0.034  0.030  209.625  35.079
M1TJ  0.097  0.066         0.028  0.048  162.333  31.468
ALM   0.271  0.161         0.159  0.172  220.692  49.663
S1BF  0.051  0.041         0.022  0.035  215.000  25.971
ALM > M1TJ: difference +0.173, Mann-Whitney p = 0.0025 (trial-fold variant p = 0.0226)
ALM > S1TJ: difference +0.196, Mann-Whitney p = 0.0007 (trial-fold variant p = 0.0319)
ALM > S1BF: difference +0.220, Mann-Whitney p = 0.0014 (trial-fold variant p = 0.0473)

Variant, one value per mouse (mean over its sessions):
  ALM (8 mice) > M1TJ (4 mice): difference +0.195, Mann-Whitney p = 0.0242
  ALM (8 mice) > S1TJ (4 mice): difference +0.210, Mann-Whitney p = 0.0141
  ALM (8 mice) > S1BF (3 mice): difference +0.248, Mann-Whitney p = 0.0121

By subject:
area  subject 
ALM   MX170903    0.642
      MX180202    0.191
      MX180203    0.336
      MX180701    0.274
      MX180804    0.400
      MX181002    0.055
      MX181101    0.263
      MX190103    0.180
M1TJ  MX181101    0.090
      MX181301    0.099
      MX181302    0.062
      MX190103    0.138
S1BF  MX180803    0.015
      MX180804    0.092
      MX181002    0.027
S1TJ  MX180804    0.105
      MX181002    0.056
      MX181101    0.052
      MX190103    0.117

print(it.groupby("area")[["time", "time_trialfolds"]].agg(["mean", "max"]).round(3).loc[FIG3].to_string())
print("Regions with mean R2 of the time to cue >= 0.1:", list(it.groupby("area").time.mean().loc[lambda x: x >= 0.1].index))

       time        time_trialfolds       
       mean    max            mean    max
area                                     
S1TJ  0.024  0.043           0.008  0.023
M1TJ  0.033  0.063           0.013  0.028
ALM   0.034  0.091           0.018  0.066
S1BF  0.008  0.019           0.004  0.010
Regions with mean R2 of the time to cue >= 0.1: []

units = []
for o in pe:
    q = o["quality"]
    for k in range(q.shape[0]):
        p6 = {key: o["peth"][key][1][k] for key in o["peth"]}
        mx = max(v.max() for v in p6.values())
        peak = max(p6, key=lambda key: p6[key].max())
        r_avg = np.nanmean([sl.rhythmicity(o["peth"][("mid", d)][0], p6[("mid", d)]) for d in ("LR", "RL")])
        units.append({"area": o["area"], "subject": o["subject"], "path": o["path"], "unit": k,
                      "single": (q[k, 0] <= sl.RPV_MAX) and (q[k, 1] <= sl.CONTAM_MAX),
                      "max_rate": mx, "peak_period": peak[0], "r_avg": r_avg})
units = pd.DataFrame(units)
print("Selected trials per PETH (min over sessions):",
      {f"{p}-{d}": min(o["n_trials"][(p, d)] for o in pe) for p in sl.PETH_WINDOWS for d in ("LR", "RL")})
paper = {"ALM": 324, "M1TJ": 233, "S1TJ": 119}
for thr in (5, 10, 15):
    act = units[units.max_rate >= thr]
    a = act[act.single].groupby("area").size()
    b = act.groupby("area").size()
    print(f"Max PETH rate >= {thr} spikes/s: (a) single units {a.reindex(FIG3).to_dict()}; (b) all units {b.reindex(FIG3).to_dict()}")
print("Paper:", paper)
act = units[(units.max_rate >= 10) & units.single]
print("\nReading (a) at 10 spikes/s by subject:")
print(act.groupby(["area", "subject"]).size().to_string())

Selected trials per PETH (min over sessions): {'init-LR': 16, 'init-RL': 13, 'mid-LR': 16, 'mid-RL': 13, 'term-LR': 16, 'term-RL': 13}
Max PETH rate >= 5 spikes/s: (a) single units {'S1TJ': 164, 'M1TJ': 259, 'ALM': 387, 'S1BF': 67}; (b) all units {'S1TJ': 264, 'M1TJ': 303, 'ALM': 426, 'S1BF': 74}
Max PETH rate >= 10 spikes/s: (a) single units {'S1TJ': 155, 'M1TJ': 223, 'ALM': 351, 'S1BF': 55}; (b) all units {'S1TJ': 251, 'M1TJ': 261, 'ALM': 387, 'S1BF': 59}
Max PETH rate >= 15 spikes/s: (a) single units {'S1TJ': 140, 'M1TJ': 200, 'ALM': 321, 'S1BF': 44}; (b) all units {'S1TJ': 226, 'M1TJ': 234, 'ALM': 353, 'S1BF': 48}
Paper: {'ALM': 324, 'M1TJ': 233, 'S1TJ': 119}

Reading (a) at 10 spikes/s by subject:
area  subject 
ALM   MX170903    39
      MX180202    46
      MX180203    38
      MX180701    47
      MX180804    27
      MX181002    52
      MX181101    68
      MX190103    34
M1TJ  MX181101    56
      MX181301    75
      MX181302    36
      MX190103    56
S1BF  MX180803    15
      MX180804    25
      MX181002    15
S1TJ  MX180804    54
      MX181002    46
      MX181101    38
      MX190103    17

tab = pd.crosstab(act.area, act.peak_period).reindex(FIG3)
tab["frac_init"] = tab["init"] / tab[["init", "mid", "term"]].sum(1)
print(tab.round(3).to_string())
for a in ["M1TJ", "S1TJ", "S1BF"]:
    t = [[tab.loc["ALM", "init"], tab.loc["ALM", ["mid", "term"]].sum()], [tab.loc[a, "init"], tab.loc[a, ["mid", "term"]].sum()]]
    print(f"ALM against {a}: fraction {tab.loc['ALM', 'frac_init']:.3f} vs {tab.loc[a, 'frac_init']:.3f}, "
          f"Fisher one-sided p = {stats.fisher_exact(t, alternative='greater')[1]:.4f}")
print("S1BF is outside the criterion; it is shown because the claim says 'the other regions'.")
print("Window lengths (s):", {p: w[1] - w[0] for p, w in sl.PETH_WINDOWS.items()},
      "; fraction peaking in each window per second of window, by area:")
print((tab[["init", "mid", "term"]].div(tab[["init", "mid", "term"]].sum(1), axis=0)
       / pd.Series({p: w[1] - w[0] for p, w in sl.PETH_WINDOWS.items()})).round(3).to_string())
allu = units[units.max_rate >= 10]
print("\nAll active units (reading b):")
print(pd.crosstab(allu.area, allu.peak_period).reindex(FIG3).to_string())

peak_period  init  mid  term  frac_init
area                                   
S1TJ           20   94    41      0.129
M1TJ           48  115    60      0.215
ALM            91  194    66      0.259
S1BF           17   14    24      0.309
ALM against M1TJ: fraction 0.259 vs 0.215, Fisher one-sided p = 0.1355
ALM against S1TJ: fraction 0.259 vs 0.129, Fisher one-sided p = 0.0006
ALM against S1BF: fraction 0.259 vs 0.309, Fisher one-sided p = 0.8275
S1BF is outside the criterion; it is shown because the claim says 'the other regions'.
Window lengths (s): {'init': 0.5, 'mid': 2.0, 'term': 1.2} ; fraction peaking in each window per second of window, by area:
peak_period   init    mid   term
area                            
S1TJ         0.258  0.303  0.220
M1TJ         0.430  0.258  0.224
ALM          0.519  0.276  0.157
S1BF         0.618  0.127  0.364

All active units (reading b):
peak_period  init  mid  term
area                        
S1TJ           32  137    82
M1TJ           53  136    72
ALM            98  216    73
S1BF           19   15    25

/tmp/ipykernel_2421/2350594777.py:2: Pandas4Warning: Starting with pandas version 4.0 all arguments of sum will be keyword-only.
  tab["frac_init"] = tab["init"] / tab[["init", "mid", "term"]].sum(1)
/tmp/ipykernel_2421/2350594777.py:11: Pandas4Warning: Starting with pandas version 4.0 all arguments of sum will be keyword-only.
  print((tab[["init", "mid", "term"]].div(tab[["init", "mid", "term"]].sum(1), axis=0)

print(act.groupby("area").r_avg.describe().round(3).loc[FIG3].to_string())
for a in ["S1TJ", "M1TJ"]:
    p = stats.mannwhitneyu(act[act.area == a].r_avg.dropna(), act[act.area == "ALM"].r_avg.dropna(), alternative="greater").pvalue
    print(f"{a} > ALM: medians {act[act.area == a].r_avg.median():.3f} vs {act[act.area == 'ALM'].r_avg.median():.3f}, "
          f"one-sided Mann-Whitney p = {p:.2g}")
print("Units without r_avg (flat PETH):", int(act.r_avg.isna().sum()))
fig, ax = plt.subplots(figsize=(4, 2.5))
for a in ["S1TJ", "M1TJ", "ALM"]:
    x = np.sort(act[act.area == a].r_avg.dropna())
    ax.plot(x, np.linspace(0, 1, len(x)), label=a)
ax.set_xlabel("r_avg")
ax.set_ylabel("cumulative fraction")
ax.legend()
plt.show()
print(f"r_avg ranges from {act.r_avg.min():.3f} to {act.r_avg.max():.3f}; all within the axes.")

      count   mean    std    min    25%    50%    75%    max
area                                                        
S1TJ  155.0  0.223  0.113  0.028  0.142  0.201  0.285  0.561
M1TJ  223.0  0.166  0.080  0.028  0.106  0.156  0.208  0.487
ALM   351.0  0.102  0.052  0.014  0.063  0.093  0.133  0.300
S1BF   55.0  0.149  0.072  0.031  0.094  0.136  0.184  0.403
S1TJ > ALM: medians 0.201 vs 0.093, one-sided Mann-Whitney p = 2.5e-36
M1TJ > ALM: medians 0.156 vs 0.093, one-sided Mann-Whitney p = 4.1e-24
Units without r_avg (flat PETH): 0

r_avg ranges from 0.014 to 0.561; all within the axes.

rw = []
for o in pe:
    for k in range(o["reward_pre"].shape[1]):
        a, b = o["reward_pre"][:, k], o["reward_post"][:, k]
        p = stats.wilcoxon(a, b).pvalue if np.any(a != b) else 1.0
        rw.append({"area": o["area"], "path": o["path"], "unit": k, "pre": a.mean(), "post": b.mean(), "p": p,
                   "trials": len(a)})
rw = pd.DataFrame(rw).merge(units[["path", "unit", "single", "max_rate"]], on=["path", "unit"])
rw = rw[(rw.max_rate >= 10) & rw.single]
for a in FIG3:
    d = rw[rw.area == a]
    dec = int(((d.p < 0.01) & (d.post < d.pre)).sum())
    inc = int(((d.p < 0.01) & (d.post > d.pre)).sum())
    pb = stats.binomtest(dec, dec + inc, 0.5, alternative="greater").pvalue if dec + inc else np.nan
    print(f"{a}: active single units {len(d)}, decrease {dec}, increase {inc}, binomial p = {pb:.3g}; "
          f"median rate before {d.pre.median():.1f}, after {d.post.median():.1f} spikes/s")
for thr in (0.001, 0.05):
    d = rw[rw.area == "ALM"]
    print(f"ALM at p < {thr}: decrease {int(((d.p < thr) & (d.post < d.pre)).sum())}, increase {int(((d.p < thr) & (d.post > d.pre)).sum())}")

S1TJ: active single units 155, decrease 50, increase 51, binomial p = 0.579; median rate before 7.2, after 8.6 spikes/s
M1TJ: active single units 223, decrease 96, increase 74, binomial p = 0.0535; median rate before 9.2, after 6.9 spikes/s
ALM: active single units 351, decrease 188, increase 92, binomial p = 4.95e-09; median rate before 8.5, after 7.0 spikes/s
S1BF: active single units 55, decrease 16, increase 15, binomial p = 0.5; median rate before 5.7, after 5.9 spikes/s
ALM at p < 0.001: decrease 175, increase 79
ALM at p < 0.05: decrease 198, increase 102

def last_position_touch(ev):
    """The first touch at the last position: the first touch after the last port-moving touch,
    leaving 60 ms (the lick detector merges contacts closer than this)."""
    t = ev["touch"][ev["touch"] > ev["pos"][-1] + 0.06] if len(ev["pos"]) else []
    return t[0] if len(t) else np.nan


dur = []
for t in tables:
    if not ((t["kind"] == "ephys" and t["area"] != "ZZ") or
            (t["kind"] == "opto VGAT-CRE Ai32" and t["area"] in ("2s 5V", "2s 2.5V"))):
        continue
    for i in range(t["trials"]["n"]):
        ev = sl.trial_events(t, i)
        if sl.trial_type(ev["kind"])[0] != "standard" or not np.isnan(t["trials"]["opto_times"][i]):
            continue
        end = last_position_touch(ev)
        dur.append({"subject": t["subject"], "kind": t["kind"], "first_to_last": end - ev["pos"][0],
                    "cue_to_last": end - ev["cue"]})
dur = pd.DataFrame(dur).dropna()
per_mouse = dur.groupby("subject")[["first_to_last", "cue_to_last"]].median()
print(f"Standard trials {len(dur)} from {dur.subject.nunique()} mice (paper, Extended Data Fig. 1h: 17 mice)")
print(f"(a) first port-moving touch to the first touch at the last position: median over mice "
      f"{per_mouse.first_to_last.median():.3f} s (mice from {per_mouse.first_to_last.min():.3f} to {per_mouse.first_to_last.max():.3f})")
print(f"(b) cue to that touch: median over mice {per_mouse.cue_to_last.median():.3f} s "
      f"(from {per_mouse.cue_to_last.min():.3f} to {per_mouse.cue_to_last.max():.3f})")
print("\nPer mouse:")
print(per_mouse.round(3).to_string())
print("\nRecording sessions only (no opto sessions):",
      dur[dur.kind == "ephys"].groupby("subject").first_to_last.median().median().round(3), "s for reading (a)")

Standard trials 25438 from 17 mice (paper, Extended Data Fig. 1h: 17 mice)
(a) first port-moving touch to the first touch at the last position: median over mice 0.992 s (mice from 0.890 to 1.138)
(b) cue to that touch: median over mice 1.324 s (from 1.150 to 1.673)

Per mouse:
          first_to_last  cue_to_last
subject                             
MX170903          1.082        1.318
MX180201          1.047        1.457
MX180202          1.091        1.673
MX180203          0.916        1.324
MX180601          0.954        1.246
MX180602          1.040        1.382
MX180701          0.985        1.439
MX180803          0.890        1.150
MX180804          1.015        1.289
MX181002          0.952        1.233
MX181003          1.039        1.323
MX181101          0.992        1.327
MX181301          1.138        1.607
MX181302          0.942        1.284
MX190103          0.992        1.428
MX200101          0.898        1.205
MX200102          1.013        1.369

Recording sessions only (no opto sessions): 1.004 s for reading (a)

learn = sorted([t for t in tables if t["kind"] == "learning" and t["path"] not in dup_paths],
               key=lambda t: (t["subject"], t["session_start"]))
print(f"Learning files {len(learn)} (left out as deposited twice: {sorted(dup_paths)})")
lt = []
for t in learn:
    for i in range(t["trials"]["n"]):
        ev = sl.trial_events(t, i)
        ty = sl.trial_type(ev["kind"])[0]
        row = {"subject": t["subject"], "session": t["session_start"], "type": ty}
        if ty == "standard":
            end = last_position_touch(ev)
            touches = np.concatenate([ev["pos"], [end]])
            row["init"] = sl.first_touch(ev) - ev["cue"]
            row["speed"] = np.mean(1 / np.diff(touches)) if np.all(np.diff(touches) > 0) else np.nan
        lt.append(row)
lt = pd.DataFrame(lt)
lt["trial_number"] = lt.groupby("subject").cumcount() + 1
print(f"Mice with learning files: {lt.subject.nunique()}; trials {len(lt)}, of which standard {int((lt.type == 'standard').sum())}")
rows = []
for m, d in lt[lt.type == "standard"].groupby("subject"):
    a, b = d.iloc[:100], d.iloc[-100:]
    rows.append({"subject": m, "sessions": d.session.nunique(), "standard_trials": len(d),
                 "init_first": a.init.median(), "init_last": b.init.median(),
                 "speed_first": a.speed.median(), "speed_last": b.speed.median()})
lc = pd.DataFrame(rows).set_index("subject")
lc["init_down"] = lc.init_last < lc.init_first
lc["speed_up"] = lc.speed_last > lc.speed_first
print(lc.round(3).to_string())
for col, (x, y) in {"initiation time": ("init_first", "init_last"), "sequence speed": ("speed_first", "speed_last")}.items():
    print(f"{col}: first 100 median over mice {lc[x].median():.3f}, last 100 {lc[y].median():.3f}; "
          f"Wilcoxon p = {stats.wilcoxon(lc[x], lc[y]).pvalue:.2g}")
print(f"Mice with lower initiation time: {int(lc.init_down.sum())}; higher speed: {int(lc.speed_up.sum())}; "
      f"both: {int((lc.init_down & lc.speed_up).sum())} of {len(lc)}")
l14 = lc.drop(index="MX181003")
print(f"Variant without MX181003 (its other file is deposited twice): {int((l14.init_down & l14.speed_up).sum())} of {len(l14)} improve "
      f"in both; Wilcoxon p = {stats.wilcoxon(l14.init_first, l14.init_last).pvalue:.2g} and "
      f"{stats.wilcoxon(l14.speed_first, l14.speed_last).pvalue:.2g}")

Learning files 572 (left out as deposited twice: ['sub-MX181002/sub-MX181002_ses-20181002T153057.nwb', 'sub-MX181003/sub-MX181003_ses-20181002T153057.nwb'])

Mice with learning files: 15; trials 140474, of which standard 114357
          sessions  standard_trials  init_first  init_last  speed_first  speed_last  init_down  speed_up
subject                                                                                                 
MX170903        14             3665       0.468      0.230        5.189       6.005       True      True
MX180201        18             4416       0.367      0.252        6.033       6.584       True      True
MX180601        32             5929       4.335      0.427        4.380       6.733       True      True
MX180602        23             4526       2.863      0.285        4.532       5.884       True      True
MX180701        29             5057       2.573      0.787        5.264       5.749       True      True
MX180803        49            14444       1.988      0.227        5.138       6.798       True      True
MX180804        59            12986       6.207      0.320        4.600       5.780       True      True
MX181002        50            12847       4.639      0.261        5.022       6.429       True      True
MX181003        34             7716       0.323      1.594        5.288       3.337      False     False
MX181101        43             8915       2.855      0.270        4.537       6.182       True      True
MX181301        40             6892       3.313      0.492        3.503       5.188       True      True
MX181302        23             4302       4.796      0.193        4.190       6.620       True      True
MX190103        32             5363       2.567      0.386        3.659       6.038       True      True
MX200101        25             6384       1.789      0.291        4.902       6.847       True      True
MX200102        44            10915       1.323      0.523        5.346       5.550       True      True
initiation time: first 100 median over mice 2.573, last 100 0.291; Wilcoxon p = 0.00043
sequence speed: first 100 median over mice 4.902, last 100 6.038; Wilcoxon p = 0.0043
Mice with lower initiation time: 14; higher speed: 14; both: 14 of 15
Variant without MX181003 (its other file is deposited twice): 14 of 14 improve in both; Wilcoxon p = 0.00012 and 0.00012

prof = []
for m, d in lt.groupby("subject"):
    n = len(d)
    final = d.iloc[-500:]
    fin_speed = final[final.type == "standard"].speed.median()
    fin_init = final[final.type == "standard"].init.median()
    a = b = np.nan
    for start in range(0, n - 99, 50):
        w = d.iloc[start:start + 100]
        w = w[w.type == "standard"]
        if len(w) < 10:
            continue
        if np.isnan(a) and w.speed.median() >= 0.9 * fin_speed:
            a = start + 100
        if np.isnan(b) and w.init.median() <= 1.1 * fin_init:
            b = start + 100
    prof.append({"subject": m, "trials_in_files": n, "a_speed": a, "b_initiation": b})
prof = pd.DataFrame(prof).set_index("subject")
print(prof.to_string())
print(f"\nMedian over mice: (a) speed reaches 90% of final at trial {prof.a_speed.median():.0f}; "
      f"(b) initiation within 110% of final at trial {prof.b_initiation.median():.0f}  (paper: about 1500)")
for f in (0.8, 0.95):
    v = []
    for m, d in lt.groupby("subject"):
        fin = d.iloc[-500:]
        fs = fin[fin.type == "standard"].speed.median()
        hit = [s + 100 for s in range(0, len(d) - 99, 50)
               if d.iloc[s:s + 100][lambda x: x.type == "standard"].speed.median() >= f * fs]
        v.append(hit[0] if hit else np.nan)
    print(f"Variant, speed reaches {int(f * 100)}% of final: median over mice {np.nanmedian(v):.0f}")
first_dates = lt.groupby("subject").session.min()
print("\nFirst learning file of each mouse:", first_dates.str[:10].to_dict())

          trials_in_files  a_speed  b_initiation
subject                                         
MX170903             3669    150.0        1200.0
MX180201             4416    100.0         150.0
MX180601             6318   1300.0        1350.0
MX180602             4527    200.0        1150.0
MX180701             5765    150.0        1100.0
MX180803            16133    300.0        1100.0
MX180804            14931    200.0        1050.0
MX181002            14187    500.0         300.0
MX181003             7716    100.0         100.0
MX181101            11112    200.0        3900.0
MX181301             8309    200.0         950.0
MX181302             4915    850.0        3250.0
MX190103             6992    650.0         950.0
MX200101            16927      NaN           NaN
MX200102            14557    100.0        1100.0

Median over mice: (a) speed reaches 90% of final at trial 200; (b) initiation within 110% of final at trial 1100  (paper: about 1500)

Variant, speed reaches 80% of final: median over mice 150

Variant, speed reaches 95% of final: median over mice 250

First learning file of each mouse: {'MX170903': '2018-02-16', 'MX180201': '2018-03-23', 'MX180601': '2018-06-17', 'MX180602': '2018-06-17', 'MX180701': '2018-08-24', 'MX180803': '2018-07-16', 'MX180804': '2018-07-15', 'MX181002': '2018-09-15', 'MX181003': '2018-09-16', 'MX181101': '2019-01-06', 'MX181301': '2019-01-08', 'MX181302': '2019-01-08', 'MX190103': '2019-03-30', 'MX200101': '2020-09-11', 'MX200102': '2020-09-18'}

# Do the learning files start at the start of training? Per mouse: the dates of the first and last
# learning file and of the first recording or opto session, the trials and trial types of the
# first file, and the initiation time and speed in its first 100 standard trials.
other_first = (pd.DataFrame([{"subject": t["subject"], "date": t["session_start"]} for t in tables if t["kind"] != "learning"])
               .groupby("subject").date.min())
ff = []
for m, d in lt.groupby("subject"):
    first = d[d.session == d.session.min()]
    std = d[d.type == "standard"]
    ff.append({"subject": m, "learning_files": d.session.nunique(), "first_file": d.session.min()[:10],
               "last_file": d.session.max()[:10], "first_other_session": str(other_first.get(m, ""))[:10],
               "first_file_trials": len(first), "first_file_types": dict(Counter(first.type)),
               "init_first100": std.init.iloc[:100].median(), "speed_first100": std.speed.iloc[:100].median(),
               "standard_in_last500": int((d.iloc[-500:].type == "standard").sum())})
ff = pd.DataFrame(ff).set_index("subject")
print(ff.round(3).to_string())
print(f"\nMice whose first file has no standard trials: {int(ff.first_file_types.map(lambda x: 'standard' not in x).sum())}; "
      f"mice whose first 100 standard trials already have a median initiation time under 1 s: {int((ff.init_first100 < 1).sum())} of {len(ff)}")
print("Mice with fewer than 10 standard trials among their last 500 trials (no final value):",
      list(ff.index[ff.standard_in_last500 < 10]))
print("\nposIndex patterns of the first learning file of each mouse (top 3):")
for m, d in lt.groupby("subject"):
    t = [t for t in learn if t["subject"] == m and t["session_start"] == d.session.min()][0]
    print(f"  {m}: {Counter(t['trials']['posIndex']).most_common(3)}")

          learning_files  first_file   last_file first_other_session  first_file_trials               first_file_types  init_first100  speed_first100  standard_in_last500
subject                                                                                                                                                                   
MX170903              14  2018-02-16  2018-03-05          2018-03-04                132  {'standard': 130, 'other': 2}          0.468           5.189                  500
MX180201              18  2018-03-23  2018-04-09          2018-04-03                 95               {'standard': 95}          0.367           6.033                  500
MX180601              32  2018-06-17  2018-07-29          2018-07-20                100              {'standard': 100}          4.335           4.380                  340
MX180602              23  2018-06-17  2018-07-14          2018-07-07                104              {'standard': 104}          2.863           4.532                  500
MX180701              29  2018-08-24  2018-09-24          2018-09-16                141              {'standard': 141}          2.573           5.264                  340
MX180803              49  2018-07-16  2018-11-28          2018-08-06                128  {'other': 1, 'standard': 127}          1.988           5.138                  341
MX180804              59  2018-07-15  2018-12-06          2018-08-06                 82               {'standard': 82}          6.207           4.600                  334
MX181002              50  2018-09-15  2019-01-01          2018-10-10                 87   {'other': 1, 'standard': 86}          4.639           5.022                  388
MX181003              34  2018-09-16  2019-04-03          2018-10-10                129              {'standard': 129}          0.323           5.288                  500
MX181101              43  2019-01-06  2019-04-02          2019-03-28                135              {'standard': 135}          2.855           4.537                  358
MX181301              40  2019-01-08  2019-05-31          2019-05-21                142              {'standard': 142}          3.313           3.503                  347
MX181302              23  2019-01-08  2019-02-12          2019-02-01                130              {'standard': 130}          4.796           4.190                  360
MX190103              32  2019-03-30  2019-05-15          2019-04-23                140              {'standard': 140}          2.567           3.659                  349
MX200101              64  2020-09-11  2020-12-03          2020-09-27                108              {'standard': 108}          1.789           4.902                    0
MX200102              62  2020-09-18  2021-01-20          2020-10-05                158              {'standard': 158}          1.323           5.346                  500

Mice whose first file has no standard trials: 0; mice whose first 100 standard trials already have a median initiation time under 1 s: 3 of 15
Mice with fewer than 10 standard trials among their last 500 trials (no final value): ['MX200101']

posIndex patterns of the first learning file of each mouse (top 3):
  MX170903: [(np.str_('[1 2 3 4 5 6]'), 66), (np.str_('[5 4 3 2 1 0]'), 64), (np.str_('[1 2 3 4 5 6 0 1 2 3 4 5 6 0]'), 1)]
  MX180201: [(np.str_('[1 2 3 4 5 6]'), 48), (np.str_('[5 4 3 2 1 0]'), 47)]
  MX180601: [(np.str_('[1 2 3 4 5 6]'), 50), (np.str_('[5 4 3 2 1 0]'), 50)]
  MX180602: [(np.str_('[1 2 3 4 5 6]'), 52), (np.str_('[5 4 3 2 1 0]'), 52)]
  MX180701: [(np.str_('[1 2 3 4 5 6]'), 71), (np.str_('[5 4 3 2 1 0]'), 70)]
  MX180803: [(np.str_('[1 2 3 4 5 6]'), 64), (np.str_('[5 4 3 2 1 0]'), 63), (np.str_('nan'), 1)]
  MX180804: [(np.str_('[1 2 3 4 5 6]'), 41), (np.str_('[5 4 3 2 1 0]'), 41)]

  MX181002: [(np.str_('[1 2 3 4 5 6]'), 43), (np.str_('[5 4 3 2 1 0]'), 43), (np.str_('nan'), 1)]
  MX181003: [(np.str_('[1 2 3 4 5 6]'), 65), (np.str_('[5 4 3 2 1 0]'), 64)]

  MX181101: [(np.str_('[1 2 3 4 5 6]'), 68), (np.str_('[5 4 3 2 1 0]'), 67)]
  MX181301: [(np.str_('[1 2 3 4 5 6]'), 71), (np.str_('[5 4 3 2 1 0]'), 71)]
  MX181302: [(np.str_('[1 2 3 4 5 6]'), 65), (np.str_('[5 4 3 2 1 0]'), 65)]
  MX190103: [(np.str_('[1 2 3 4 5 6]'), 70), (np.str_('[5 4 3 2 1 0]'), 70)]
  MX200101: [(np.str_('[1 2 3 4 5 6]'), 54), (np.str_('[5 4 3 2 1 0]'), 54)]
  MX200102: [(np.str_('[1 2 3 4 5 6]'), 79), (np.str_('[5 4 3 2 1 0]'), 79)]

bk = []
for t in tables:
    if t["kind"] == "learning":
        continue
    types = Counter(sl.trial_type(k)[0] for k in t["trials"]["posIndex"])
    if types["backtrack"]:
        bk.append({"path": t["path"], "subject": t["subject"], "kind": t["kind"], "area": t["area"], "standard": types["standard"],
                   "backtrack": types["backtrack"]})
bk = pd.DataFrame(bk)
readings = {"(a) all non-learning files with backtracking": bk,
            "(b) without air and earplug sessions": bk[~bk.area.isin(["air", "earplug"])],
            "(c) as (b), without opto efficiency sessions": bk[~bk.area.isin(["air", "earplug", "efficiency"])]}
for name, d in readings.items():
    print(f"{name}: sessions {len(d)}, mice {d.subject.nunique()}, standard {d.standard.sum()}, backtracking {d.backtrack.sum()}")
print("Paper: 47 sessions, 9 mice, 7458 standard, 2695 backtracking")
print("\nReading (c) by kind of file:")
print(readings["(c) as (b), without opto efficiency sessions"].groupby(["kind", "area"])[["standard", "backtrack"]]
      .agg(["size", "sum"]).to_string())
print("Mice:", sorted(readings["(c) as (b), without opto efficiency sessions"].subject.unique()))
# Variant, not part of the criterion: the Methods leave out 8 sessions "where only touch onsets were
# correctly registered". Sessions whose touches are mostly very short (5th percentile of the touch
# duration under 10 ms) are candidates for these.
short = []
for t in tables:
    if t["kind"] == "learning" or t["area"] in ("air", "earplug", "efficiency"):
        continue
    if not any(sl.trial_type(k)[0] == "backtrack" for k in t["trials"]["posIndex"]):
        continue
    d = np.concatenate(t["trials"]["lickOff_times"]) - np.concatenate(t["trials"]["lickOn_times"])
    short.append({"path": t["path"], "subject": t["subject"], "p5_ms": 1000 * np.percentile(d, 5),
                  "median_ms": 1000 * np.median(d)})
short = pd.DataFrame(short)
print("\nTouch duration (ms) per session of reading (c), 5th percentile, lowest 10:")
print(short.sort_values("p5_ms").head(10).round(1).to_string(index=False))
keep = set(short[short.p5_ms >= 10].path)
d = bk[~bk.area.isin(["air", "earplug", "efficiency"]) & bk.path.isin(keep)]
print(f"Variant, (c) without sessions whose 5th percentile touch is under 10 ms: sessions {len(d)}, mice {d.subject.nunique()}, "
      f"standard {d.standard.sum()}, backtracking {d.backtrack.sum()}")

(a) all non-learning files with backtracking: sessions 74, mice 15, standard 13248, backtracking 4849
(b) without air and earplug sessions: sessions 59, mice 10, standard 10890, backtracking 4030
(c) as (b), without opto efficiency sessions: sessions 55, mice 9, standard 10365, backtracking 3839
Paper: 47 sessions, 9 mice, 7458 standard, 2695 backtracking

Reading (c) by kind of file:
                           standard       backtrack      
                               size   sum      size   sum
kind               area                                  
ephys              ALM            7  1437         7   513
                   M1TJ           9  1470         9   467
                   S1BF           5  1080         5   512
                   S1L            6   969         6   331
                   S1TJ           8  1685         8   627
                   others        14  2639        14  1034
opto VGAT-CRE Ai32 2s 2.5V        1   292         1    95
                   2s 5V          5   793         5   260
Mice: ['MX180601', 'MX180701', 'MX180803', 'MX180804', 'MX181002', 'MX181101', 'MX181301', 'MX181302', 'MX190103']

Touch duration (ms) per session of reading (c), 5th percentile, lowest 10:
                                                      path  subject  p5_ms  median_ms
sub-MX181302/sub-MX181302_ses-20190210T163258_behavior.nwb MX181302    3.0        5.0
sub-MX180803/sub-MX180803_ses-20181128T130139_behavior.nwb MX180803    3.0        4.0
sub-MX180804/sub-MX180804_ses-20181206T135700_behavior.nwb MX180804    4.0       40.0
sub-MX180804/sub-MX180804_ses-20181205T150915_behavior.nwb MX180804    4.0       39.0
sub-MX180803/sub-MX180803_ses-20181122T135902_behavior.nwb MX180803    4.0       10.0
sub-MX180803/sub-MX180803_ses-20181123T122105_behavior.nwb MX180803    5.0       36.0
sub-MX181302/sub-MX181302_ses-20190201T153422_behavior.nwb MX181302    7.0       46.0
sub-MX181302/sub-MX181302_ses-20190212T154624_behavior.nwb MX181302    7.0       41.0
sub-MX181002/sub-MX181002_ses-20181224T154648_behavior.nwb MX181002   14.0       34.0
sub-MX181002/sub-MX181002_ses-20181219T132039_behavior.nwb MX181002   14.0       33.0
Variant, (c) without sessions whose 5th percentile touch is under 10 ms: sessions 47, mice 8, standard 8891, backtracking 3207

# Variant, not part of the criterion: the paper's trial exclusions on the 47 sessions above. The
# Methods remove the first and last trial of every session and, from analyses of tongue kinematics
# (Extended Data Fig. 2c), trials not finished before the video stopped. Initiation under 1 s is
# the Methods' trial selection for Figure 3 and is shown as a further step.
v47 = [a for a in assets if a["path"] in set(d.path)]
spans = {r["path"]: r["spans"] for r in sl.run_parallel(sl.video_spans, v47, workers=sl.cpu_workers())}
steps = Counter()
for t in tables:
    if t["path"] not in spans:
        continue
    n = t["trials"]["n"]
    for i in range(n):
        ev = sl.trial_events(t, i)
        ty = sl.trial_type(ev["kind"])[0]
        if ty not in ("standard", "backtrack"):
            continue
        steps[("all", ty)] += 1
        if i in (0, n - 1):
            continue
        steps[("without first and last trial", ty)] += 1
        sp = spans[t["path"]].get(i)
        if sp is None or not len(ev["pos"]) or ev["pos"][-1] > sp[1]:
            continue
        steps[("and sequence finished within the video", ty)] += 1
        if not sl.first_touch(ev) - ev["cue"] < 1.0:
            continue
        steps[("and initiation under 1 s", ty)] += 1
st = pd.Series(steps).unstack()
print(f"47-session variant ({len(spans)} files, {sum(len(v) > 0 for v in spans.values())} with video):")
print(st.loc[["all", "without first and last trial", "and sequence finished within the video", "and initiation under 1 s"]].to_string())
print("Paper: 7458 standard, 2695 backtracking (tolerance 6712 to 8204 and 2426 to 2965)")

47-session variant (47 files, 47 with video):
                                        backtrack  standard
all                                          3207      8891
without first and last trial                 3196      8809
and sequence finished within the video       3104      8627
and initiation under 1 s                     2894      8058
Paper: 7458 standard, 2695 backtracking (tolerance 6712 to 8204 and 2426 to 2965)

ml = []
for s in fig3_sessions:
    for i in range(min(s["trials"]["n"], len(s["video"]))):
        ev = sl.trial_events(s, i)
        ty = sl.trial_type(ev["kind"])[0]
        if ty not in ("standard", "backtrack") or len(ev["pos"]) < 5:
            continue
        v = s["video"][i]
        if ev["pos"][4] > v["start"] + len(v["out"]) / v["rate"]:
            continue
        lk = s["licks"][i]
        inside = lk[(lk[:, 0] > ev["pos"][3]) & (lk[:, 0] < ev["pos"][4])]
        miss = {tol: sum(1 for on, off, *_ in inside if not np.any((ev["touch"] >= on - tol) & (ev["touch"] <= off + tol)))
                for tol in (0.0, 0.0125, 0.025)}
        ml.append({"subject": s["subject"], "path": s["path"], "type": ty, "missed": miss[0.0125],
                   "missed_tol0": miss[0.0], "missed_tol25": miss[0.025], "time": ev["pos"][4] - ev["pos"][3]})
ml = pd.DataFrame(ml)
ml = ml[ml.path.isin(ml[ml.type == "backtrack"].path.unique())]
print(f"Sessions with backtracking trials: {ml.path.nunique()}, mice {ml.subject.nunique()}")
print(ml.groupby("type").agg(trials=("missed", "size"), mean_missed=("missed", "mean"), median_missed=("missed", "median"),
                             mean_time_s=("time", "mean")).round(3).to_string())
p = stats.mannwhitneyu(ml[ml.type == "backtrack"].missed, ml[ml.type == "standard"].missed, alternative="greater").pvalue
print(f"Backtracking > standard: Mann-Whitney p = {p:.2g}")
print("\nBy mouse (mean missed licks):")
print(ml.groupby(["subject", "type"]).missed.mean().unstack().round(2).to_string())
mb = ml[ml.type == "backtrack"]
print("\nVariant, touch tolerance, same trials: mean missed licks in backtracking trials "
      f"{mb.missed_tol0.mean():.3f} (0 ms), {mb.missed.mean():.3f} (12.5 ms), {mb.missed_tol25.mean():.3f} (25 ms)")

Sessions with backtracking trials: 29, mice 8
           trials  mean_missed  median_missed  mean_time_s
type                                                      
backtrack    2090        1.854            1.0        0.525
standard     5590        0.113            0.0        0.189
Backtracking > standard: Mann-Whitney p = 0

By mouse (mean missed licks):
type      backtrack  standard
subject                      
MX180701       2.00      0.03
MX180803       2.11      0.12
MX180804       0.97      0.03
MX181002       1.63      0.19
MX181101       2.08      0.08
MX181301       2.26      0.19
MX181302       5.39      0.05
MX190103       2.14      0.11

Variant, touch tolerance, same trials: mean missed licks in backtracking trials 1.854 (0 ms), 1.854 (12.5 ms), 1.854 (25 ms)

cl = []
for s in fig3_sessions:
    vid_ili = np.concatenate([np.diff(l[:, 2]) for l in s["licks"] if len(l) > 1])
    vmax = 1.5 * np.median(vid_ili[vid_ili < 1])
    tmax = sl.session_ili_max(s)
    for i in range(min(s["trials"]["n"], len(s["video"]))):
        ev = sl.trial_events(s, i)
        if sl.trial_type(ev["kind"])[0] != "standard" or np.isnan(ev["water"]):
            continue
        tw = ev["touch"][ev["touch"] >= ev["water"]]
        if not len(tw):
            continue
        w = tw[0]
        k = 0                                          # (b) touches after the first water touch
        while k + 1 < len(tw) and tw[k + 1] - tw[k] < tmax:
            k += 1
        k2 = 0                                         # (b) variant: video lick threshold for touches
        while k2 + 1 < len(tw) and tw[k2 + 1] - tw[k2] < vmax:
            k2 += 1
        row = {"subject": s["subject"], "touches_after": k, "touches_after_video_thr": k2,
               "thr_touch_ms": 1000 * tmax, "thr_video_ms": 1000 * vmax}
        lk = s["licks"][i]
        v = s["video"][i]
        end = v["start"] + len(v["out"]) / v["rate"]
        j0 = np.where((lk[:, 0] - 0.0125 <= w) & (lk[:, 1] + 0.0125 >= w))[0]
        if len(j0):
            j = j0[0]
            while j + 1 < len(lk) and lk[j + 1, 2] - lk[j, 2] < vmax:
                j += 1
            censored = end - lk[j, 2] < vmax
            row.update(video_after=np.nan if censored else j - j0[0], censored=censored)
        cl.append(row)
cl = pd.DataFrame(cl)
print(f"Standard trials with a water touch: {len(cl)}; with the water lick in the video: {int(cl.censored.notna().sum())}; "
      f"video ended during the bout: {int(cl.censored.fillna(False).astype(bool).sum())}")
va = cl.video_after.dropna()
print(f"(a) video licks after the first water lick: median {va.median():.1f}, mean {va.mean():.2f} (trials {len(va)})")
print(f"(b) touches after the first water touch: median {cl.touches_after.median():.1f}, mean {cl.touches_after.mean():.2f}")
print(f"(c) video licks counting the first water lick: median {va.median() + 1:.1f}")
print("\nBy mouse, median of (a) and (b):")
print(cl.groupby("subject")[["video_after", "touches_after"]].median().to_string())
print("\nDistribution over trials (fraction of trials at each count):")
dist = pd.DataFrame({"(a) video licks": va.value_counts(normalize=True), "(b) touches": cl.touches_after.value_counts(normalize=True),
                     "(b) with the video lick threshold": cl.touches_after_video_thr.value_counts(normalize=True)}).fillna(0).sort_index()
print(dist.loc[0:12].round(3).to_string())
print(f"Bout thresholds over sessions (ms): from touch intervals median {cl.thr_touch_ms.median():.1f}, "
      f"from video lick intervals median {cl.thr_video_ms.median():.1f}")
print(f"(b) with the video lick threshold: median {cl.touches_after_video_thr.median():.1f}, mean {cl.touches_after_video_thr.mean():.2f}")
print(f"Mean of (a) {va.mean():.2f}; trials with (a) of 4, 5 or 6: {va.isin([4, 5, 6]).mean():.3f}")

Standard trials with a water touch: 7106; with the water lick in the video: 6952; video ended during the bout: 88
(a) video licks after the first water lick: median 6.0, mean 6.38 (trials 6864)
(b) touches after the first water touch: median 4.0, mean 3.87
(c) video licks counting the first water lick: median 7.0

By mouse, median of (a) and (b):
          video_after  touches_after
subject                             
MX170903          6.0            3.0
MX180202          6.0            2.0
MX180203          7.0            5.0
MX180701          5.0            2.0
MX180803          7.0            5.0
MX180804          5.0            3.0
MX181002          8.0            6.0
MX181101          5.0            2.0
MX181301          4.0            3.0
MX181302         10.0            7.0
MX190103          6.0            4.0

Distribution over trials (fraction of trials at each count):
      (a) video licks  (b) touches  (b) with the video lick threshold
0.0             0.000        0.148                              0.174
1.0             0.000        0.100                              0.101
2.0             0.005        0.066                              0.062
3.0             0.040        0.132                              0.126
4.0             0.128        0.137                              0.131
5.0             0.202        0.144                              0.138
6.0             0.196        0.105                              0.103
7.0             0.167        0.078                              0.076
8.0             0.113        0.045                              0.045
9.0             0.068        0.025                              0.024
10.0            0.041        0.013                              0.013
11.0            0.021        0.004                              0.004
12.0            0.013        0.001                              0.001
Bout thresholds over sessions (ms): from touch intervals median 243.0, from video lick intervals median 238.0
(b) with the video lick threshold: median 4.0, mean 3.76
Mean of (a) 6.38; trials with (a) of 4, 5 or 6: 0.526

cs_rows = []
for s in fig3_sessions:
    for i in range(min(s["trials"]["n"], len(s["video"]))):
        ev = sl.trial_events(s, i)
        if sl.trial_type(ev["kind"])[0] != "standard":
            continue
        ft = sl.first_touch(ev)
        lk = s["licks"][i]
        j = np.where((lk[:, 0] - 0.0125 <= ft) & (lk[:, 1] + 0.0125 >= ft))[0]
        if not len(j) or np.isnan(lk[j[0], 4]):
            continue
        side = np.sign(lk[j[0], 4])
        before = lk[(lk[:, 0] >= ev["cue"]) & (lk[:, 0] < lk[j[0], 0]), 4]
        cs_rows.append({"subject": s["subject"], "licks_before": len(before),
                        "explored": bool(np.any(np.sign(before) == -side)),
                        "explored_5deg": bool(np.any(before * side < -5))})
csd = pd.DataFrame(cs_rows)
fr = csd.groupby("subject").agg(trials=("explored", "size"), without_exploring=("explored", lambda x: 1 - x.mean()),
                                without_5deg=("explored_5deg", lambda x: 1 - x.mean()),
                                licks_before_mean=("licks_before", "mean"))
print(fr.round(3).to_string())
print(f"\nMice with fraction above 0.5: {int((fr.without_exploring > 0.5).sum())} of {len(fr)}; "
      f"with the 5 degree margin: {int((fr.without_5deg > 0.5).sum())} of {len(fr)}; "
      f"pooled fraction {1 - csd.explored.mean():.3f}")

          trials  without_exploring  without_5deg  licks_before_mean
subject                                                             
MX170903     261              0.923         0.973              0.966
MX180202     382              0.626         0.788              1.717
MX180203     750              0.900         0.980              1.665
MX180701     166              0.741         0.849              0.723
MX180803     245              0.763         0.845              1.016
MX180804     949              0.920         0.952              1.065
MX181002    1668              0.562         0.733              1.610
MX181101    1163              0.735         0.835              1.727
MX181301     394              0.731         0.802              1.165
MX181302     320              0.569         0.647              1.262
MX190103     709              0.568         0.705              1.444

Mice with fraction above 0.5: 11 of 11; with the 5 degree margin: 11 of 11; pooled fraction 0.714

ctrl_assets = [a for a in assets if sl.group_label(a["description"]) in
               (("behav", "earplug"), ("behav", "numbing"), ("behav", "air"))]
ctrl = sl.load_many(ctrl_assets, workers=sl.cpu_workers(), spikes=False)
for s in ctrl:
    s["kind"], s["area"] = sl.group_label(s["description"])
cm = sl.run_parallel(sl.sequence_missed_licks, ctrl, workers=sl.cpu_workers())
rows = []
for o, s in zip(cm, ctrl):
    r = pd.DataFrame(o["rows"])
    rows.append({"condition": o["area"], "subject": o["subject"], "session": o["session_start"][:16],
                 "standard_trials": sum(sl.trial_type(k)[0] == "standard" for k in s["trials"]["posIndex"]),
                 "with_video": len(r), "missed_mean": r.missed.mean() if len(r) else np.nan,
                 "init_median": r.init.median() if len(r) else np.nan})
cs_ = pd.DataFrame(rows).sort_values(["condition", "subject", "session"])
print(cs_.round(3).to_string(index=False))


def boot_p(diff, rng, n=10000):
    """Paired one-tailed bootstrap: fraction of resampled mean differences at or below 0."""
    diff = np.asarray(diff)
    b = rng.choice(diff, (n, len(diff))).mean(1)
    return float(np.mean(b <= 0))

condition  subject          session  standard_trials  with_video  missed_mean  init_median
      air MX190201 2019-10-08T11:28              249         246        0.980        0.263
      air MX190602 2019-12-02T13:24              197         195        0.554        0.356
      air VC010102 2019-11-16T11:24              158         149        0.584        0.366
      air VC010103 2019-11-14T14:25              194         192        0.266        0.239
      air WO010401 2019-11-28T13:26              184         182        0.934        0.539
      air WO010402 2019-11-28T14:13              194         187        0.209        0.349
  earplug MX190602 2019-12-03T14:29              147         132        0.227        1.691
  earplug MX190602 2019-12-04T12:36              169         163        0.147        0.253
  earplug VC010102 2019-11-17T12:19              113         106        0.123        1.875
  earplug VC010102 2019-11-18T13:34              129         128        0.109        0.387
  earplug VC010103 2019-11-15T12:27              149         134        0.142        0.644
  earplug VC010103 2019-11-16T14:08              176         175        0.143        0.249
  earplug WO010401 2019-11-29T12:41              150         146        0.486        1.405
  earplug WO010401 2019-12-01T14:00              160         158        0.114        0.359
  earplug WO010402 2019-11-29T13:48              100          76        0.276        1.568
  earplug WO010402 2019-11-30T17:10              138         112        0.580        0.378
  numbing MX200102 2021-01-17T16:18              143         131        2.687        0.572
  numbing MX200102 2021-01-20T17:00              202         183        0.956        0.607
  numbing MX200103 2021-01-17T17:29              109          96        2.667        0.484
  numbing MX200103 2021-01-20T18:02              186         181        0.232        0.350
  numbing MX210203 2021-02-26T17:35              190         186        2.586        0.314
  numbing MX210203 2021-02-28T17:16              209         204        1.049        0.332
  numbing MX210502 2021-05-03T14:52              205         196        0.224        0.355
  numbing MX210601 2021-05-30T13:40               54          42        1.857        0.331
  numbing MX210601 2021-05-31T14:10               84          58        0.207        0.249
  numbing MX210602 2021-05-29T12:19              120         109        0.991        0.614
  numbing MX210602 2021-05-30T12:30              131         128        0.664        0.341

# ed1i: under each of the 2^5 assignments of 'earplug' to one of each mouse's two sessions, the
# paired one-tailed bootstrap of earplug > sham in missed licks. The paper's top panel (time to
# first touch, p < 0.001) suggests one assignment; it is shown too.
rng = np.random.default_rng(31)
ep = cs_[cs_.condition == "earplug"]
pairs = [g.sort_values("session") for _, g in ep.groupby("subject")]
print(f"Mice with two earplug sessions: {sum(len(g) == 2 for g in pairs)} of {len(pairs)}")
d0 = np.array([g.missed_mean.iloc[1] - g.missed_mean.iloc[0] for g in pairs if len(g) == 2])
ps = [boot_p(d0 * np.array(sg), rng) for sg in itertools.product([1, -1], repeat=len(d0))]
print(f"Missed licks, later minus earlier session per mouse: {np.round(d0, 3).tolist()}")
print(f"Over all {len(ps)} assignments: smallest p = {min(ps):.4f}; assignments with p < 0.05: {sum(p < 0.05 for p in ps)}")
di = np.array([g.init_median.iloc[1] - g.init_median.iloc[0] for g in pairs if len(g) == 2])
sign = np.sign(di)
print(f"Median initiation time, later minus earlier session: {np.round(di, 3).tolist()}")
print(f"Assignment by the slower initiation (earplug = the session with the longer time to first touch): "
      f"missed licks earplug minus sham {np.round(d0 * sign, 3).tolist()}, mean {np.mean(d0 * sign):+.3f}, "
      f"p = {boot_p(d0 * sign, rng):.4f}")

Mice with two earplug sessions: 5 of 5
Missed licks, later minus earlier session per mouse: [-0.08, -0.013, 0.001, -0.372, 0.304]
Over all 32 assignments: smallest p = 0.0000; assignments with p < 0.05: 4
Median initiation time, later minus earlier session: [-1.438, -1.488, -0.395, -1.045, -1.19]
Assignment by the slower initiation (earplug = the session with the longer time to first touch): missed licks earplug minus sham [0.08, 0.013, -0.001, 0.372, -0.304], mean +0.032, p = 0.3586

# ed1j: the air flow changes after the 100th trial, on first in two mice and off first in four.
# Missed licks in the standard trials before and after the 100th trial of each session; under each
# of the 15 ways to choose the two on-first mice, the paired one-tailed bootstrap of air on > off.
rng = np.random.default_rng(32)
air = [(o, s) for o, s in zip(cm, ctrl) if o["area"] == "air"]
halves = []
for o, s in air:
    r = pd.DataFrame(o["rows"])
    halves.append({"subject": o["subject"], "first_100": r[r.trial < 100].missed.mean(), "after_100": r[r.trial >= 100].missed.mean(),
                   "trials_first": int((r.trial < 100).sum()), "trials_after": int((r.trial >= 100).sum())})
halves = pd.DataFrame(halves).set_index("subject")
print(halves.round(3).to_string())
dd = (halves.after_100 - halves.first_100).to_numpy()
ps = []
for on_first in itertools.combinations(range(len(dd)), 2):
    sgn = np.array([-1 if k in on_first else 1 for k in range(len(dd))])
    ps.append(boot_p(dd * sgn, rng))
print(f"Over all {len(ps)} assignments: smallest p = {min(ps):.4f}; assignments with p < 0.05: {sum(p < 0.05 for p in ps)}")

          first_100  after_100  trials_first  trials_after
subject                                                   
MX190201      1.293      0.769            99           147
MX190602      0.488      0.600            80           115
VC010102      1.051      0.070            78            71
VC010103      0.267      0.265            75           117
WO010401      0.800      1.039            80           102
WO010402      0.295      0.147            78           109
Over all 15 assignments: smallest p = 0.0179; assignments with p < 0.05: 1

# ed1k: lidocaine should increase missed licks. Without labels, only whether each mouse's two
# sessions differ can be shown, over the trials of the first 30 min of each session (the Methods
# limit the analysis to about 30 min).
rows = []
for (m, g) in cs_[cs_.condition == "numbing"].groupby("subject"):
    per_s = []
    for _, r in g.sort_values("session").iterrows():
        o = [x for x in cm if x["subject"] == m and x["session_start"][:16] == r.session][0]
        t = pd.DataFrame(o["rows"])
        t = t[t.cue <= t.cue.min() + 1800] if len(t) else t
        per_s.append(t.missed.to_numpy() if len(t) else np.array([]))
    row = {"subject": m, "sessions": len(per_s), "missed_session1": per_s[0].mean() if len(per_s[0]) else np.nan,
           "trials1": len(per_s[0])}
    if len(per_s) == 2:
        row.update(missed_session2=per_s[1].mean() if len(per_s[1]) else np.nan, trials2=len(per_s[1]),
                   p_two_sided=stats.mannwhitneyu(per_s[0], per_s[1]).pvalue if len(per_s[0]) and len(per_s[1]) else np.nan)
    rows.append(row)
nk = pd.DataFrame(rows).set_index("subject")
print(nk.round(4).to_string())
print(f"Mice with two sessions: {int((nk.sessions == 2).sum())}; whose sessions differ at p < 0.05: {int((nk.p_two_sided < 0.05).sum())}; "
      f"ratio of the larger to the smaller mean, median {np.nanmedian(np.maximum(nk.missed_session1, nk.missed_session2) / np.minimum(nk.missed_session1, nk.missed_session2)):.2f}")
print(f"\nWhole notebook up to here: {(time.time() - started) / 60:.1f} min")

          sessions  missed_session1  trials1  missed_session2  trials2  p_two_sided
subject                                                                            
MX200102         2           2.7381      126           0.9931    145.0       0.0018
MX200103         2           2.6667       96           0.2320    181.0       0.0000
MX210203         2           2.5860      186           1.0663    196.0       0.0000
MX210502         1           0.2222      117              NaN      NaN          NaN
MX210601         2           1.8571       42           0.2703     37.0       0.0002
MX210602         2           0.9908      109           0.6641    128.0       0.2259
Mice with two sessions: 5; whose sessions differ at p < 0.05: 4; ratio of the larger to the smaller mean, median 2.76

Whole notebook up to here: 21.0 min

Claim	Tested or blocked
ephys-unit-counts	tested
fig3-session-counts	tested
fig3-units-per-session	tested
fig3a-trial-count	tested
fig3c-length	tested
fig3c-velocity	tested
fig3c-angle	tested
fig3c-identity	tested
fig3c-time	tested
ed5b-variables-uncorrelated	tested
ed5e-variance-explained	tested
ed5f-axes-orthogonal	tested
fig3e-cca	tested
fig3f-backtracking-svm	tested
fig4g-zigzag-context	tested
ed7f-target-position-alm	tested
ed7h-no-time-coding	tested
ed4m-active-neuron-counts	tested
ed4l-alm-initiation	tested
ed4o-rhythmicity	tested
ed6-alm-reward-drop	tested
ed1h-sequence-duration	tested
ed1l-learning-curves	tested
ed1l-proficiency-1500	tested
ed2-backtracking-counts	tested
ed2e-missed-licks	tested
ed6b-consummatory-licks	tested
ed7b-correct-side	tested
fig2b-s1tj-angle-variability, fig2b-almm1tj-angle-length, fig2c-almm1tj-initiation, fig2c-termination, ed3lm-half-power	blocked: metadata insufficient (inhibited region not stated)
ed3f-inhibition-efficiency	blocked: data missing (no spike waveforms, no fibre distance)
ed1i-earplug, ed1j-odor-masking, ed1k-tongue-numbing	blocked: metadata insufficient (condition not stated)
ed1f-lick-force	blocked: metadata insufficient (no force calibration)
ed4j-nnmf-clusters	blocked: methods underspecified

Claim	Outcome
ephys-unit-counts	partially reproduced
fig3-session-counts	reproduced
fig3-units-per-session	reproduced
fig3a-trial-count	reproduced
fig3c-length, fig3c-velocity, fig3c-angle, fig3c-identity, fig3c-time	reproduced
ed5b-variables-uncorrelated	reproduced
ed5e-variance-explained	reproduced
ed5f-axes-orthogonal	partially reproduced
fig3e-cca	reproduced
fig3f-backtracking-svm	partially reproduced
fig4g-zigzag-context	reproduced
ed7f-target-position-alm	reproduced
ed7h-no-time-coding	reproduced
ed4m-active-neuron-counts	partially reproduced
ed4l-alm-initiation	partially reproduced
ed4o-rhythmicity	reproduced
ed6-alm-reward-drop	reproduced
ed1h-sequence-duration	reproduced
ed1l-learning-curves	reproduced
ed1l-proficiency-1500	not reproduced
ed2-backtracking-counts	not reproduced
ed2e-missed-licks	reproduced
ed6b-consummatory-licks	not reproduced
ed7b-correct-side	reproduced
Opto claims (Figure 2b,c; Extended Data Fig. 3l,m), five claims	blocked: metadata insufficient
ed3f-inhibition-efficiency	blocked: data missing
ed1i-earplug, ed1j-odor-masking, ed1k-tongue-numbing	blocked: metadata insufficient
ed1f-lick-force	blocked: metadata insufficient
ed4j-nnmf-clusters	blocked: methods underspecified

Reproducing Xu et al. (2022) from dandiset 000239¶

Data check¶

ephys-unit-counts¶

fig3-session-counts¶

fig3-units-per-session¶

fig3a-trial-count¶

Figure 3c: decoding of the five variables¶

fig3c-length¶

fig3c-velocity¶

fig3c-angle¶

fig3c-identity¶

fig3c-time¶

ed5b-variables-uncorrelated¶

ed5e-variance-explained¶

ed5f-axes-orthogonal¶

Seed variants of the Figure 3 fits¶

fig3e-cca¶

fig3f-backtracking-svm¶

fig4g-zigzag-context¶

ed7f-target-position-alm¶

ed7h-no-time-coding¶

ed4m-active-neuron-counts¶

ed4l-alm-initiation¶

ed4o-rhythmicity¶

ed6-alm-reward-drop¶

ed1h-sequence-duration¶

ed1l-learning-curves¶

ed1l-proficiency-1500¶

ed2-backtracking-counts¶

ed2e-missed-licks¶

ed6b-consummatory-licks¶

ed7b-correct-side¶

Blocked control claims: what can be checked without the labels¶

Summary¶