import collections
import re
import time

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

import rodgers as R

t_start = time.time()
meta = R.version_metadata()
assets = R.list_assets()
nwb_assets = sorted([a for a in assets if a["path"].endswith(".nwb")], key=lambda a: a["path"])
mkv_assets = sorted([a for a in assets if a["path"].endswith(".mkv")], key=lambda a: a["path"])
S = R.summarize_all(nwb_assets)
V = pd.DataFrame(R.video_headers(mkv_assets))
V["duration_tag"] = V.stream_tags.map(R.tag_duration)
V["frames_stream"] = V.frames   # frame count of the stream as the container reports it
V["frames"] = np.round(V.duration_tag * 30).astype(int)   # nominal 30 frames/s in the file
V["nwb_name"] = V.path.str.rsplit("/", n=1).str[0]
print(f"version {meta['version']} (doi {meta['metadata'].get('doi')}), {meta['asset_count']} assets, {meta['size'] / 1e12:.2f} TB, status {meta['status']}")
print(f"{len(nwb_assets)} NWB files and {len(mkv_assets)} videos read in {time.time() - t_start:.0f} s")

version 0.220904.1554 (doi 10.48324/dandi.000231/0.220904.1554), 4228 assets, 2.00 TB, status Valid
115 NWB files and 115 videos read in 352 s

# per-session table and pooled trials
def task_by_rule(s):
    """Task-rules criterion: detection when any trial has stimulus "nothing" or the metadata
    names detection; otherwise discrimination."""
    if (s["trials"].stimulus == "nothing").any() or R.task_from_metadata(s) == "detection":
        return "detection"
    return "discrimination"


def task_from_description(s):
    d = s["experiment_description"].lower()
    return "detection" if "shape detection" in d else "discrimination" if "shape discrimination" in d else None


sess = pd.DataFrame({
    "path": [s["path"] for s in S],
    "folder": [s["folder"] for s in S],
    "subject": [s["subject"] for s in S],
    "session": [s["session"] for s in S],
    "types": [s["types"] for s in S],
    "ecephys": [s["ecephys"] is not None for s in S],
    "task_description": [task_from_description(s) for s in S],
    "task_keywords": [";".join(k for k in s["keywords"] if "shape" in k) for s in S],
    "task": [task_by_rule(s) for s in S],
    "n_trials": [len(s["trials"]) for s in S],
    "start": pd.to_datetime([s["session_start_time"] for s in S], utc=True),
})
T = pd.concat([s["trials"].assign(path=s["path"], subject=s["subject"], task=t)
               for s, t in zip(S, sess.task)], ignore_index=True)
by_path = {s["path"]: s for s in S}

print("Counts")
print(f"  subjects (folders): {sess.folder.nunique()}, subject_id in metadata: "
      f"{len({s['subject_meta']['subject_id'] for s in S})}, folder equals subject_id in "
      f"{sum(s['folder'] == s['subject_meta']['subject_id'] for s in S)} of {len(S)} files")
print(f"  sessions (NWB files): {len(sess)}, with ecephys: {sess.ecephys.sum()}")
print(sess.groupby(["subject", "task"]).agg(sessions=("path", "size"), ecephys=("ecephys", "sum"),
                                            trials=("n_trials", "sum")).to_string())
print("\nTrials by outcome, stimulus, servo_position, choice, rewarded_side")
for c in ["outcome", "stimulus", "servo_position", "choice", "rewarded_side"]:
    print(f"  {c}: {T[c].value_counts().to_dict()}")
print(f"  flags: direct_delivery {T.direct_delivery.sum()}, stim_is_random False {(~T.stim_is_random).sum()}, "
      f"optogenetic {T.optogenetic.sum()}, ignore_trial {T.ignore_trial.sum()} of {len(T)}")
print("\nRanges (s, common video timebase)")
for c in ["start_time", "stop_time", "response_window_open_time", "choice_time"]:
    print(f"  {c}: {T[c].min():.3f} to {T[c].max():.3f}")
w_rows = [dict(path=s["path"], whisker=w, **{k: d.get(k) for k in
          ["n_processed", "t_first", "t_last", "n_pose", "n_contacts", "n_interpolated"]})
          for s in S for w, d in s["whiskers"].items()]
W = pd.DataFrame(w_rows)
print(f"  processed whisker series: {W.n_processed.min()} to {W.n_processed.max()} frames, "
      f"spanning {W.t_first.min():.3f} to {W.t_last.max():.3f} s")
units_proc = collections.Counter(u for s in S for d in s["whiskers"].values() for u in d.get("units", {}).items())
print(f"  units of processed series: {dict(units_proc)}")
print(f"  pose units {collections.Counter(u for s in S for d in s['whiskers'].values() for u in d.get('pose_unit', []))}, "
      f"conversions {sorted({c for s in S for d in s['whiskers'].values() for c in d.get('pose_conversion', [])})}")
E = [s for s in S if s["ecephys"]]
print(f"  ecephys: channels {min(e['ecephys']['n_channels'] for e in E)} to {max(e['ecephys']['n_channels'] for e in E)}, "
      f"rate {min(e['ecephys']['rate'] for e in E):.4f} to {max(e['ecephys']['rate'] for e in E):.4f} Hz, "
      f"dtype {collections.Counter(e['ecephys']['dtype'] for e in E)}, conversion {collections.Counter(e['ecephys']['conversion'] for e in E)} V")
print(f"  units: {sum(len(e['ecephys']['unit_id']) for e in E)} in {len(E)} sessions, "
      f"{sum(e['ecephys']['unit_inhibitory'].sum() for e in E)} inhibitory; layers "
      f"{dict(collections.Counter(l for e in E for l in e['ecephys']['unit_layer']))}")
print(f"  electrode depth y: {min(e['ecephys']['el_y'].min() for e in E)} to {max(e['ecephys']['el_y'].max() for e in E)} um")
print(f"  videos: {V.format.unique()}, codec {V.codec.unique()}, nominal rate {V.average_rate.unique()}, "
      f"sizes {V.groupby(['width', 'height']).size().to_dict()}")

print("\nChecks of the metadata against the data")
print("  task named by experiment_description:", sess.task_description.value_counts().to_dict())
print("  task keyword in general/keywords:", sess.task_keywords.value_counts().to_dict())
print("  sessions whose keyword disagrees with the description:",
      int((sess.task_description != sess.task_keywords.str.replace("shape ", "")).sum()))
print("  stimulus values outside Table 6 {convex, concave, nothing}:",
      T.loc[~T.stimulus.isin(["convex", "concave", "nothing"]), "stimulus"].value_counts().to_dict())
same_r = sum(np.array_equal(s["events"]["rewards"]["right"], s["events"]["licks"]["right"]) for s in S)
same_l = sum(np.array_equal(s["events"]["rewards"]["left"], s["events"]["licks"]["left"]) for s in S)
print(f"  rewards/right identical to licks/right in {same_r} of {len(S)} sessions; "
      f"rewards/left identical to licks/left in {same_l}")
n_rr = sum(len(s["events"]["rewards"]["right"]) for s in S)
n_rl = sum(len(s["events"]["rewards"]["left"]) for s in S)
print(f"  reward times: right {n_rr}, left {n_rl}; correct trials rewarded right {((T.outcome == 'correct') & (T.rewarded_side == 'right')).sum()}, left {((T.outcome == 'correct') & (T.rewarded_side == 'left')).sum()}")
print(f"  neural recording stored with timestamps array: {sum(e['ecephys']['has_timestamps'] for e in E)} of {len(E)} "
      "(Table 2 describes a timestamps dataset)")

Counts
  subjects (folders): 15, subject_id in metadata: 15, folder equals subject_id in 115 of 115 files
  sessions (NWB files): 115, with ecephys: 44
                        sessions  ecephys  trials
subject task                                     
219CR   discrimination         9        0    2024
221CR   discrimination        13        2    4137
228CR   detection              7        5    1651
229CR   discrimination         8        6    2012
231CR   discrimination        12        4    3054
234CR   detection              9        0    2690
245CR   detection              6        5    1456
255CR   detection              1        1     248
267CR   detection              7        6    1932
KF119   discrimination         6        1    1710
KF132   discrimination         7        2    2275
KF134   discrimination         9        7    3168
KM101   discrimination         7        0    1688
KM131   discrimination         9        5    2845
KM91    discrimination         5        0    1460

Trials by outcome, stimulus, servo_position, choice, rewarded_side
  outcome: {'correct': 24962, 'error': 7123, 'spoil': 265}
  stimulus: {'concave': 13979, 'convex': 13089, 'nothing': 3749, 'flatter convex': 1058, 'flatter concave': 475}
  servo_position: {'close': 10854, 'far': 10811, 'medium': 10685}
  choice: {'left': 16218, 'right': 15867, 'nogo': 265}
  rewarded_side: {'right': 16278, 'left': 16072}
  flags: direct_delivery 46, stim_is_random False 7574, optogenetic 2482, ignore_trial 9445 of 32350

Ranges (s, common video timebase)
  start_time: 0.035 to 4945.925
  stop_time: 6.710 to 4953.650
  response_window_open_time: 3.057 to 4949.360
  choice_time: 3.321 to 4949.361
  processed whisker series: 326070 to 979096 frames, spanning 0.000 to 4958.995 s
  units of processed series: {('angle', 'degrees'): 460, ('base_x', 'mm'): 460, ('base_y', 'mm'): 460, ('bend', '1 / mm'): 460, ('interpolated', 'bool'): 460, ('tip_x', 'mm'): 460, ('tip_y', 'mm'): 460}
  pose units Counter({'mm': 460}), conversions [0.04026712407335282, 0.04046107493905031, 0.04079220087387068, 0.04106968604119641, 0.04134806269863899, 0.04153846153846154, 0.04161239537603088, 0.04163712785689264, 0.04169169491508716, 0.04199823714803065, 0.04216691570992364, 0.04224943421489739, 0.04238978953521744, 0.042452620285415264, 0.04247888234810981, 0.04249992712638771, 0.04251573125229859, 0.04263748482520342, 0.04280863447390447, 0.04281401616158352, 0.042900400140783466, 0.04299275848264754, 0.0430966917612516, 0.04314067921006985, 0.04329570268419228, 0.04365641250653994, 0.043707863806076885, 0.04373078951644225, 0.04374800749316417, 0.0441261304060914, 0.04416743872845641, 0.04419109563359587, 0.0443517802364878, 0.04438772657244647, 0.04445988223086794, 0.04452633497595642, 0.04459916981467541, 0.04467236324982891, 0.04469684125892593, 0.04484455912101959, 0.044875518433731026, 0.044975020814061954, 0.04506890795433316, 0.04514444248993608, 0.04516970519056461, 0.04529028487803037, 0.045392578141888344, 0.045405413594751094, 0.04557976482527636, 0.04560576583413677, 0.045861629866987386, 0.0459013770619681, 0.04600119894947903, 0.04607481697700146, 0.04613531351670509, 0.04630461798847739, 0.046407099747155, 0.0466002335770648, 0.04671168422013739, 0.04688743803618518, 0.04696540072601875, 0.04705089606509522, 0.04711532486128489, 0.04713686001508462, 0.047361133042579585, 0.0476924625369598, 0.047699904600286204, 0.04796259192376326, 0.048175895481497026, 0.04823736395368353, 0.048682245103286774, 0.0487695248891345, 0.04905848444464303, 0.049254002261427635, 0.04992301766027062, 0.050137741307804004, 0.05031152949374527, 0.05051363757817357, 0.050718201081660766, 0.051162430111725186, 0.05203867575318974, 0.05244943656729227]
  ecephys: channels 56 to 64, rate 29999.9354 to 30000.0199 Hz, dtype Counter({'int16': 44}), conversion Counter({1.9499999999999999e-07: 44}) V
  units: 999 in 44 sessions, 172 inhibitory; layers {'5b': 249, '5a': 186, '4': 269, '6': 158, '2/3': 114, '1': 12, '6b': 11}
  electrode depth y: -53.200000000000045 to 1437.5000000000057 um
  videos: <StringArray>
['matroska,webm']
Length: 1, dtype: str, codec <StringArray>
['h264']
Length: 1, dtype: str, nominal rate [30.], sizes {(640, 550): 111, (640, 600): 4}

Checks of the metadata against the data
  task named by experiment_description: {'discrimination': 85, 'detection': 30}
  task keyword in general/keywords: {'shape discrimination': 115}
  sessions whose keyword disagrees with the description: 30
  stimulus values outside Table 6 {convex, concave, nothing}: {'flatter convex': 1058, 'flatter concave': 475}
  rewards/right identical to licks/right in 115 of 115 sessions; rewards/left identical to licks/left in 0
  reward times: right 315782, left 13036; correct trials rewarded right 12447, left 12515
  neural recording stored with timestamps array: 0 of 44 (Table 2 describes a timestamps dataset)

keys = collections.defaultdict(list)
for s in S:
    keys[("trials", tuple(np.round(s["trials"].start_time.to_numpy()[:20], 3)))].append(s["path"])
    keys[("whisker frames", tuple(d.get("n_processed") for d in s["whiskers"].values()))].append(s["path"])
for v in V.itertuples():
    keys[("video frames", v.frames)].append(v.path)
dups = {k: v for k, v in keys.items() if len(v) > 1}
print(f"groups sharing a key: {len(dups)}")
for k, v in dups.items():
    print(" ", k[0], k[1] if k[0] != "trials" else "", v)
    if k[0] == "video frames":
        sz = V.set_index("path").loc[v, "size"].tolist()
        tr_n = [len(by_path[p.rsplit("/", 1)[0] + ".nwb"]["trials"]) for p in v]
        print(f"    file sizes {sz}, trials {tr_n}")
print(f"videos whose duration is a whole number of seconds: {(V.duration_tag % 1 == 0).sum()} of {len(V)}")
print("pairs matching on trial times or whisker frame counts:", sum(k[0] != "video frames" for k in dups))
print(f"sessions per subject and date: max {sess.groupby(['subject', sess.start.dt.date]).size().max()}")

groups sharing a key: 7
  video frames 419400 ['sub-219CR/sub-219CR_ses-20190403T123013_behavior+image/a969df11-1d5e-4088-822e-e1a6730d74a2_external_file_0.mkv', 'sub-231CR/sub-231CR_ses-20190921T144923_behavior+ecephys+image/f80f5af1-6140-4d88-9734-52fcfe0dca47_external_file_0.mkv']
    file sizes [3658207689, 2862700447], trials [251, 246]
  video frames 625200 ['sub-229CR/sub-229CR_ses-20191012T141045_behavior+ecephys+image/3bb4fe4e-123a-4aaf-832e-dc5282a5ede8_external_file_0.mkv', 'sub-234CR/sub-234CR_ses-20190827T122436_behavior+image/257e02e1-7588-4bc4-9cd4-060b5e94b361_external_file_0.mkv']
    file sizes [4355957504, 4843524424], trials [280, 326]
  video frames 516000 ['sub-229CR/sub-229CR_ses-20191013T152728_behavior+ecephys+image/3e9479d0-ac03-4831-b7ed-7ed2e97aeae1_external_file_0.mkv', 'sub-267CR/sub-267CR_ses-20200326T151949_behavior+ecephys+image/dae16b81-9404-4436-bb17-887b8d422d5d_external_file_0.mkv']
    file sizes [3689304865, 4013962062], trials [257, 261]
  video frames 525600 ['sub-229CR/sub-229CR_ses-20191014T171024_behavior+ecephys+image/93ba6e76-2f50-4cde-b5cd-5c821a8506e5_external_file_0.mkv', 'sub-267CR/sub-267CR_ses-20200330T154638_behavior+ecephys+image/32d77d21-7071-4fd0-86f4-38ad13b09134_external_file_0.mkv']
    file sizes [3728118514, 4222049627], trials [269, 331]
  video frames 507600 ['sub-231CR/sub-231CR_ses-20190927T145645_behavior+ecephys+image/ee6d5168-ffea-4956-9b24-5c5405e554c4_external_file_0.mkv', 'sub-KM101/sub-KM101_ses-20170712T104720_behavior+image/2d79fd63-697b-455a-8b7e-3d6d565af385_external_file_0.mkv']
    file sizes [3405991144, 6828821267], trials [311, 231]
  video frames 558000 ['sub-234CR/sub-234CR_ses-20190905T123046_behavior+image/6078df1e-52cd-4398-bd53-6d4368543a81_external_file_0.mkv', 'sub-267CR/sub-267CR_ses-20200325T173201_behavior+image/c23b425e-4854-496a-a17b-6f49a5591849_external_file_0.mkv']
    file sizes [4261203332, 4137342755], trials [289, 308]
  video frames 394800 ['sub-KM101/sub-KM101_ses-20170717T115845_behavior+image/74fabd19-3e85-4666-962b-2fbbe875d40d_external_file_0.mkv', 'sub-KM131/sub-KM131_ses-20180116T184757_behavior+ecephys+image/672f9fce-72c4-4e68-9f2a-3017fc1cbe99_external_file_0.mkv']
    file sizes [4963420399, 5103315829], trials [204, 161]
videos whose duration is a whole number of seconds: 115 of 115
pairs matching on trial times or whisker frame counts: 0
sessions per subject and date: max 2

def all_licks(s):
    ev = s["events"]["licks"]
    t = np.concatenate([ev[k] for k in ev])
    side = np.concatenate([[k] * len(ev[k]) for k in ev])
    o = np.argsort(t, kind="stable")
    return t[o], side[o]


rows = []
for s in S:
    tr = s["trials"]
    lt, _ = all_licks(s)
    ns = tr[tr.outcome != "spoil"]
    def frac_lick_at(times):
        i = np.searchsorted(lt, times - 0.005)
        ok = (i < len(lt)) & (np.abs(lt[np.clip(i, 0, len(lt) - 1)] - times) <= 0.005)
        return ok.mean()
    c2 = s["whiskers"]["C2"]
    gb, gm = c2["gap_begin"], c2["gap_missing"]
    def frac_gap(starts):
        return np.mean([((np.abs(gb - x) <= 0.2) & (gm >= 0.1)).any() for x in starts])
    rows.append(dict(path=s["path"], lick_at_choice=frac_lick_at(ns.choice_time.to_numpy()),
                     lick_at_choice_shift=frac_lick_at(ns.choice_time.to_numpy() + 1.0),
                     gap_at_start=frac_gap(tr.start_time.to_numpy()),
                     gap_at_start_shift=frac_gap(tr.start_time.to_numpy() + 2.0)))
clock = pd.DataFrame(rows)
print("Fraction of non-spoiled trials with a lick within 5 ms of choice_time (true, times + 1 s):")
print(f"  median over sessions {clock.lick_at_choice.median():.3f}, min {clock.lick_at_choice.min():.3f}; "
      f"shifted median {clock.lick_at_choice_shift.median():.3f}, max {clock.lick_at_choice_shift.max():.3f}")
print("Fraction of trials with a C2 gap of 100 ms or more beginning within 0.2 s of start_time (true, starts + 2 s):")
print(f"  median {clock.gap_at_start.median():.3f}, min {clock.gap_at_start.min():.3f}; "
      f"shifted median {clock.gap_at_start_shift.median():.3f}, max {clock.gap_at_start_shift.max():.3f}")

# contacts against the stimulus, detection sessions
rows = []
for s, task in zip(S, sess.task):
    if task != "detection":
        continue
    tr = s["trials"].reset_index(drop=True)
    st = np.concatenate([d["contact_start"] for d in s["whiskers"].values() if "contact_start" in d])
    lo = tr.response_window_open_time.to_numpy() - 2.0
    hi = tr.response_window_open_time.to_numpy()
    st.sort()
    n = np.searchsorted(st, hi) - np.searchsorted(st, lo)
    shape = (tr.stimulus != "nothing").to_numpy()
    shifted = np.roll(shape, -1)[:-1]
    rows.append(dict(path=s["path"], shape=n[shape].mean(), nothing=n[~shape].mean(),
                     shape_shift=n[:-1][shifted].mean(), nothing_shift=n[:-1][~shifted].mean()))
cc = pd.DataFrame(rows)
print(f"Detection sessions ({len(cc)}): contacts per trial in the 2 s before the response window,")
print(f"  shape {cc['shape'].median():.2f} against nothing {cc.nothing.median():.2f} (medians over sessions); "
      f"shape > nothing in {(cc['shape'] > cc.nothing).sum()} of {len(cc)}")
print(f"  labels shifted by one trial: shape {cc.shape_shift.median():.2f} against nothing {cc.nothing_shift.median():.2f}; "
      f"shape > nothing in {(cc.shape_shift > cc.nothing_shift).sum()} of {len(cc)}")

# spikes against contact onsets, ecephys sessions
a0, a1 = R.CONTACT_WINDOW
rows = []
for s in E:
    on = np.sort(np.concatenate([d["contact_start"] for d in s["whiskers"].values() if "contact_start" in d]))
    sp = np.sort(np.concatenate([u["spikes"] for u in s["ecephys"]["units"]]))
    def ratio(x):
        after = (np.searchsorted(sp, x + a1) - np.searchsorted(sp, x + a0)).sum()
        before = (np.searchsorted(sp, x - a0) - np.searchsorted(sp, x - a1)).sum()
        return after / max(before, 1)
    rows.append(dict(path=s["path"], ratio=ratio(on), ratio_shift=ratio(on + 1.0)))
nc = pd.DataFrame(rows)
print(f"Ecephys sessions ({len(nc)}): spikes 5 to 25 ms after contact onsets over 25 to 5 ms before,")
print(f"  median ratio {nc.ratio.median():.2f} (range {nc.ratio.min():.2f} to {nc.ratio.max():.2f}), above 1.1 in "
      f"{(nc.ratio > 1.1).sum()} of {len(nc)}; onsets + 1 s: median {nc.ratio_shift.median():.2f} "
      f"(range {nc.ratio_shift.min():.2f} to {nc.ratio_shift.max():.2f})")
print("Sessions where the neural clock check does not show a contact response (ratio <= 1.1):",
      nc.loc[nc.ratio <= 1.1, "path"].tolist())

Fraction of non-spoiled trials with a lick within 5 ms of choice_time (true, times + 1 s):
  median over sessions 0.802, min 0.525; shifted median 0.026, max 0.080
Fraction of trials with a C2 gap of 100 ms or more beginning within 0.2 s of start_time (true, starts + 2 s):
  median 1.000, min 1.000; shifted median 0.000, max 0.004
Detection sessions (30): contacts per trial in the 2 s before the response window,
  shape 4.01 against nothing 0.00 (medians over sessions); shape > nothing in 30 of 30
  labels shifted by one trial: shape 1.83 against nothing 2.19; shape > nothing in 8 of 30

Ecephys sessions (44): spikes 5 to 25 ms after contact onsets over 25 to 5 ms before,
  median ratio 1.44 (range 1.01 to 2.86), above 1.1 in 43 of 44; onsets + 1 s: median 0.99 (range 0.66 to 1.18)
Sessions where the neural clock check does not show a contact response (ratio <= 1.1): ['sub-KF134/sub-KF134_ses-20180209T181601_behavior+ecephys+image.nwb']

mice, nsess = sess.folder.nunique(), len(sess)
trials = {"a, all rows": len(T), "b, ignore_trial False": int((~T.ignore_trial).sum())}
fr_b = sum(max(d["n_pose"] for d in s["whiskers"].values()) for s in S)
fr_span = sum(int(round(max(d["t_last"] for d in s["whiskers"].values()) / R.FRAME)) + 1 for s in S)
frames = {"a, video files": int(V.frames.sum()), "b, raw pose series": fr_b}
print(f"mice {mice} (expected 15), sessions {nsess} (expected 115)")
for k, v in trials.items():
    print(f"trials {k}: {v}, {100 * (v - 18514) / 18514:+.1f}% from 18514")
for k, v in frames.items():
    print(f"frames {k}: {v}, {100 * (v - 63979800) / 63979800:+.2f}% from 63979800, hours {v / 200 / 3600:.2f}")
print(f"variant, span of the processed series (first to last frame): {fr_span} frames, {fr_span / 720000:.3f} h")
print(f"videos whose duration tag is not a whole number of 1/30 s frames: "
      f"{(np.abs(V.duration_tag * 30 - V.frames) > 1e-3).sum()}; container duration equals the tag in "
      f"{(np.abs(V.container_duration - V.duration_tag) < 0.01).sum()} of {len(V)}")
print(f"frame count reported by the video streams: {(V.frames_stream > 0).sum()} of {len(V)} report one; "
      f"equal to duration x 30 in {(V.frames_stream == V.frames).sum()}; sum of those reported {int(V.frames_stream.sum())}")
# other trial counts, for information
for name, m in {"stim_is_random True": T.stim_is_random,
                "random and not spoiled": T.stim_is_random & (T.outcome != "spoil"),
                "ignore_trial False, stimulus in Table 6": ~T.ignore_trial & T.stimulus.isin(["convex", "concave", "nothing"]),
                "ignore_trial False in sessions without optogenetics": ~T.ignore_trial & ~T.groupby("path").optogenetic.transform("any")}.items():
    print(f"  other count, {name}: {int(m.sum())}")
exact_tr = any(v == 18514 for v in trials.values())
exact_fr = any(v == 63979800 for v in frames.values())
near_tr = any(abs(v - 18514) <= 0.1 * 18514 for v in trials.values())
near_fr = any(abs(v - 63979800) <= 0.1 * 63979800 for v in frames.values())
ok_ms = mice == 15 and nsess == 115
verdict = ("reproduced" if ok_ms and exact_tr and exact_fr else
           "partially_reproduced" if ok_ms and near_tr and near_fr else "not_reproduced")
print("criterion:", verdict)

mice 15 (expected 15), sessions 115 (expected 115)
trials a, all rows: 32350, +74.7% from 18514
trials b, ignore_trial False: 22905, +23.7% from 18514
frames a, video files: 63979800, +0.00% from 63979800, hours 88.86
frames b, raw pose series: 63082162, -1.40% from 63979800, hours 87.61
variant, span of the processed series (first to last frame): 63979788 frames, 88.861 h
videos whose duration tag is not a whole number of 1/30 s frames: 0; container duration equals the tag in 115 of 115
frame count reported by the video streams: 0 of 115 report one; equal to duration x 30 in 0; sum of those reported 0
  other count, stim_is_random True: 24776
  other count, random and not spoiled: 24590
  other count, ignore_trial False, stimulus in Table 6: 21752
  other count, ignore_trial False in sessions without optogenetics: 17287
criterion: not_reproduced

n_name = int(sess.types.str.contains("ecephys").sum())
n_data = sum(("extracellular array recording" in s["acquisition"]) and s["ecephys"] is not None
             and len(s["ecephys"]["unit_id"]) > 0 for s in S)
print(f"files named with ecephys: {n_name}; files with the ElectricalSeries and a non-empty units table: {n_data}")
print(f"name and content agree in every file: {all(('ecephys' in s['types']) == (s['ecephys'] is not None) for s in S)}")
print("criterion:", "reproduced" if n_name == 44 and n_data == 44 else "not_reproduced")

files named with ecephys: 44; files with the ElectricalSeries and a non-empty units table: 44
name and content agree in every file: True
criterion: reproduced

span_a = np.array([e["ecephys"]["ts_last"] - e["ecephys"]["ts_first"] for e in E])
src = [a for a in assets if a["path"].startswith("sourcedata/") and a["path"].endswith(".continuous")]
src_df = pd.DataFrame({"dir": [a["path"].rsplit("/", 1)[0] for a in src],
                       "file": [a["path"].rsplit("/", 1)[1] for a in src],
                       "samples": [R.openephys_samples(a["size"]) for a in src]})
per_dir = src_df.groupby("dir").agg(n_files=("file", "size"), smin=("samples", "min"), smax=("samples", "max"))
OE = pd.DataFrame(R.openephys_records_all(sorted([a for a in src if a["path"].endswith("/100_CH1.continuous")], key=lambda a: a["path"])))
OE["dir"] = OE.path.str.rsplit("/", n=1).str[0]
OE["n_segments"] = OE.segments.map(len)
OE["main_records"] = OE.segments.map(lambda x: max(n for _, n in x))
OE["other_records"] = OE.n_records - OE.main_records
print(f"(a) NWB recordings: {span_a.sum() / 3600:.3f} h over {len(span_a)} sessions "
      f"({span_a.min() / 60:.1f} to {span_a.max() / 60:.1f} min)")
print(f"Open Ephys files: {len(per_dir)} sessions; channel files per session {per_dir.n_files.min()} to {per_dir.n_files.max()}; "
      f"all channels of a session equal in length: {(per_dir.smin == per_dir.smax).all()}")
print(f"record headers of CH1: recording numbers seen {OE.recordings.map(tuple).value_counts().to_dict()}; "
      f"contiguous segments per file {OE.n_segments.value_counts().to_dict()}")
for r in OE[OE.n_segments > 1].itertuples():
    print(f"  {r.dir.split('/')[-1]}: recordings {r.recordings}, segments (first record, records) {r.segments}, "
          f"left out {r.other_records * 1024 / 30000:.2f} s")
h_b = OE.main_records.sum() * 1024 / 30000 / 3600
h_whole = per_dir.smax.sum() / 30000 / 3600
print(f"(b) main recordings: {h_b:.3f} h over {len(OE)} sessions; the other segments hold "
      f"{OE.other_records.sum() * 1024 / 30000:.1f} s in {(OE.n_segments > 1).sum()} sessions")
print(f"variant, whole Open Ephys files: {h_whole:.3f} h")
for name, h in [("a", span_a.sum() / 3600), ("b", h_b), ("whole files", h_whole)]:
    print(f"  {name}: {h:.3f} h, {h - 33.8:+.3f} h from 33.8 ({100 * (h - 33.8) / 33.8:+.2f}%)")
nwb_ch = [e["ecephys"]["n_channels"] for e in E]
print(f"NWB channels per session {min(nwb_ch)} to {max(nwb_ch)} (broken channels excluded)")
hs = [span_a.sum() / 3600, h_b]
print("criterion:", "reproduced" if any(abs(h - 33.8) <= 0.05 for h in hs) else
      "partially_reproduced" if any(abs(h - 33.8) <= 3.38 for h in hs) else "not_reproduced")

(a) NWB recordings: 33.656 h over 44 sessions (32.6 to 77.0 min)
Open Ephys files: 44 sessions; channel files per session 64 to 64; all channels of a session equal in length: True
record headers of CH1: recording numbers seen {(0,): 43, (0, 1): 1}; contiguous segments per file {1: 37, 2: 7}
  sub-KF132_ses-20180223T152924_behavior+ecephys+image: recordings [0], segments (first record, records) [(0, 15), (15, 119980)], left out 0.51 s
  sub-KF134_ses-20180207T200156_behavior+ecephys+image: recordings [0], segments (first record, records) [(0, 15), (15, 112248)], left out 0.51 s
  sub-KF134_ses-20180208T171905_behavior+ecephys+image: recordings [0], segments (first record, records) [(0, 14), (14, 135418)], left out 0.48 s
  sub-KF134_ses-20180210T191835_behavior+ecephys+image: recordings [0, 1], segments (first record, records) [(0, 1164), (1164, 96069)], left out 39.73 s
  sub-KF134_ses-20180212T183227_behavior+ecephys+image: recordings [0], segments (first record, records) [(0, 15), (15, 104022)], left out 0.51 s
  sub-KM131_ses-20180116T184757_behavior+ecephys+image: recordings [0], segments (first record, records) [(0, 15), (15, 57264)], left out 0.51 s
  sub-KM131_ses-20180117T194023_behavior+ecephys+image: recordings [0], segments (first record, records) [(0, 15), (15, 73345)], left out 0.51 s
(b) main recordings: 33.770 h over 44 sessions; the other segments hold 42.8 s in 7 sessions
variant, whole Open Ephys files: 33.782 h
  a: 33.656 h, -0.144 h from 33.8 (-0.42%)
  b: 33.770 h, -0.030 h from 33.8 (-0.09%)
  whole files: 33.782 h, -0.018 h from 33.8 (-0.05%)
NWB channels per session 56 to 64 (broken channels excluded)
criterion: reproduced

print("task by rule against task in the experiment description:")
print(pd.crosstab(sess.task, sess.task_description).to_string())
sess["task_nothing"] = np.where([(s["trials"].stimulus == "nothing").any() for s in S], "detection", "discrimination")
print("task by \"nothing\" alone against task in the experiment description:")
print(pd.crosstab(sess.task_nothing, sess.task_description).to_string())


def expected_side(task, stim):
    base = stim.replace("flatter ", "")
    if task == "discrimination":
        return {"convex": "right", "concave": "left"}.get(base)
    return {"convex": "right", "concave": "right", "nothing": "left"}.get(base)


T["expected_side"] = [expected_side(t, s) for t, s in zip(T.task, T.stimulus)]
T["table6_stim"] = T.stimulus.isin(["convex", "concave", "nothing"])
T["follows"] = T.rewarded_side == T.expected_side
for name, m in [("trials to analyse, Table 6 stimuli", ~T.ignore_trial & T.table6_stim),
                ("trials to analyse, flatter shapes by their shape", ~T.ignore_trial & ~T.table6_stim),
                ("all trials, Table 6 stimuli", T.table6_stim),
                ("all trials, every stimulus", pd.Series(True, index=T.index))]:
    sub = T[m]
    per = sub.groupby("path").follows.mean()
    print(f"{name}: {sub.follows.sum()} of {len(sub)} follow the rule ({sub.follows.mean():.4f}); "
          f"sessions where every trial follows {int((per == 1).sum())} of {len(per)}; lowest session {per.min():.4f}")
print(pd.crosstab([T.task, T.stimulus], T.rewarded_side).to_string())
main = T[~T.ignore_trial & T.table6_stim].groupby("path").follows.mean()
print("criterion:", "reproduced" if (main == 1).all() else "partially_reproduced" if (main >= 0.95).all() else "not_reproduced")

task by rule against task in the experiment description:
task_description  detection  discrimination
task                                       
detection                30               0
discrimination            0              85
task by "nothing" alone against task in the experiment description:
task_description  detection  discrimination
task_nothing                               
detection                30               0
discrimination            0              85
trials to analyse, Table 6 stimuli: 21752 of 21752 follow the rule (1.0000); sessions where every trial follows 115 of 115; lowest session 1.0000
trials to analyse, flatter shapes by their shape: 1153 of 1153 follow the rule (1.0000); sessions where every trial follows 6 of 6; lowest session 1.0000
all trials, Table 6 stimuli: 30817 of 30817 follow the rule (1.0000); sessions where every trial follows 115 of 115; lowest session 1.0000
all trials, every stimulus: 32350 of 32350 follow the rule (1.0000); sessions where every trial follows 115 of 115; lowest session 1.0000
rewarded_side                    left  right
task           stimulus                     
detection      concave              0   2131
               convex               0   2097
               nothing           3749      0
discrimination concave          11848      0
               convex               0  10992
               flatter concave    475      0
               flatter convex       0   1058
criterion: reproduced

n_nothing = T.groupby("path").apply(lambda d: int((d.stimulus == "nothing").sum()), include_groups=False)
tab = sess.set_index("path").assign(n_nothing=n_nothing)
print(tab.groupby("task_description").n_nothing.agg(["size", lambda x: int((x > 0).sum()), "sum"])
      .rename(columns={"size": "sessions", "<lambda_0>": "sessions with nothing", "sum": "nothing trials"}).to_string())
bad = tab[(tab.task_description == "discrimination") & (tab.n_nothing > 0)]
print("discrimination sessions with a nothing trial:", len(bad))
print("keywords naming a task, by the task of the experiment description:")
print(pd.crosstab(sess.task_description, sess.task_keywords).to_string())
print("criterion:", "reproduced" if len(bad) == 0 else "not_reproduced")

                  sessions  sessions with nothing  nothing trials
task_description                                                 
detection               30                     30            3749
discrimination          85                      0               0
discrimination sessions with a nothing trial: 0
keywords naming a task, by the task of the experiment description:
task_keywords     shape discrimination
task_description                      
detection                           30
discrimination                      85
criterion: reproduced

per_mouse = sess.groupby("subject").agg(tasks=("task", lambda x: sorted(set(x))),
                                        description=("task_description", lambda x: sorted(set(x))),
                                        sessions=("path", "size"))
print(per_mouse.to_string())
print("mice per task:", sess.groupby("task").subject.nunique().to_dict())
print("criterion:", "reproduced" if (per_mouse.tasks.map(len) == 1).all() else "not_reproduced")

                    tasks       description  sessions
subject                                              
219CR    [discrimination]  [discrimination]         9
221CR    [discrimination]  [discrimination]        13
228CR         [detection]       [detection]         7
229CR    [discrimination]  [discrimination]         8
231CR    [discrimination]  [discrimination]        12
234CR         [detection]       [detection]         9
245CR         [detection]       [detection]         6
255CR         [detection]       [detection]         1
267CR         [detection]       [detection]         7
KF119    [discrimination]  [discrimination]         6
KF132    [discrimination]  [discrimination]         7
KF134    [discrimination]  [discrimination]         9
KM101    [discrimination]  [discrimination]         7
KM131    [discrimination]  [discrimination]         9
KM91     [discrimination]  [discrimination]         5
mice per task: {'detection': 5, 'discrimination': 10}
criterion: reproduced

shape_tr = T[T.stimulus != "nothing"]
vals = set(shape_tr.servo_position)
print("servo_position over trials with a shape:", shape_tr.servo_position.value_counts().to_dict())
print("servo_position on nothing trials:", T.loc[T.stimulus == "nothing", "servo_position"].value_counts().to_dict())
pos_per = shape_tr.groupby("path").servo_position.apply(lambda x: ",".join(sorted(set(x))))
print("positions present per session:", pos_per.value_counts().to_dict())
print("criterion:", "reproduced" if vals == {"close", "medium", "far"} else "not_reproduced")

servo_position over trials with a shape:

 {'close': 9564, 'far': 9543, 'medium': 9494}
servo_position on nothing trials: {'close': 1290, 'far': 1268, 'medium': 1191}

positions present per session: {'close,far,medium': 115}
criterion: reproduced

d = T.response_window_open_time - T.start_time
print(f"response window open minus start_time: median {d.median():.3f} s, 5th to 95th percentile "
      f"{d.quantile(0.05):.3f} to {d.quantile(0.95):.3f} s, range {d.min():.3f} to {d.max():.3f} s")
print(d.groupby(T.servo_position).describe()[["count", "50%", "min", "max"]].round(3).to_string())
motor_cols = [c for c in T.columns if "motor" in c or "servo" in c]
print("trials columns naming the motor or servo:", motor_cols)

response window open minus start_time: median 3.023 s, 5th to 95th percentile 3.020 to 3.443 s, range 3.009 to 4.258 s
                  count    50%    min    max
servo_position                              
close           10854.0  3.023  3.009  4.257
far             10811.0  3.023  3.011  4.257
medium          10685.0  3.023  3.010  4.258
trials columns naming the motor or servo: ['servo_position']

rows = []
for s in S:
    lt, ls = all_licks(s)
    tr = s["trials"]
    for r in tr.itertuples():
        if r.outcome == "spoil":
            rows.append(dict(path=s["path"], spoil=True, ok_time=abs(r.choice_time - r.response_window_open_time - 45) <= 0.005,
                             ok_side=r.choice == "nogo", ok_tol=None, d=r.choice_time - r.response_window_open_time))
            continue
        i = np.searchsorted(lt, r.response_window_open_time)
        j = np.searchsorted(lt, r.response_window_open_time - 0.005)
        def check(k):
            if k >= len(lt):
                return False, False, np.nan
            return abs(lt[k] - r.choice_time) <= 0.005, ls[k] == r.choice or ls[k] == "both", lt[k] - r.choice_time
        ok_t, ok_s, dd = check(i)
        ok_t2, _, _ = check(j)
        rows.append(dict(path=s["path"], spoil=False, ok_time=ok_t, ok_side=ok_s, ok_tol=ok_t2, d=dd))
CT = pd.DataFrame(rows)
ns, sp = CT[~CT.spoil], CT[CT.spoil]
print(f"non-spoiled trials: {len(ns)}; first lick at or after the window within 5 ms of choice_time: "
      f"{ns.ok_time.mean():.4f}; side matches choice: {ns.ok_side.mean():.4f}; both: {(ns.ok_time & ns.ok_side).mean():.4f}")
print(f"  variant, first lick from 5 ms before the window: {ns.ok_tol.astype(bool).mean():.4f}")
print(f"  first lick minus choice_time: median {ns.d.median() * 1000:.2f} ms, 1st to 99th percentile "
      f"{ns.d.quantile(0.01) * 1000:.1f} to {ns.d.quantile(0.99) * 1000:.1f} ms; no lick after the window on {ns.d.isna().sum()} trials")
print(f"spoiled trials: {len(sp)}; choice_time 45 s after the window (within 5 ms): {sp.ok_time.mean():.4f}; "
      f"choice nogo: {sp.ok_side.mean():.4f}; choice_time NaN on {sp.d.isna().sum()} of them")
# where the first lick in the window is not at choice_time: where is choice_time, and was the mouse already licking?
mis = []
for s in S:
    lt, _ = all_licks(s)
    for r in s["trials"][s["trials"].outcome != "spoil"].itertuples():
        i = np.searchsorted(lt, r.response_window_open_time)
        if i < len(lt) and abs(lt[i] - r.choice_time) > 0.005:
            mis.append(dict(c_minus_w=r.choice_time - r.response_window_open_time,
                            lick_before=bool(((lt > r.response_window_open_time - 0.1) & (lt < r.response_window_open_time)).any()),
                            any_lick_at_choice=bool((np.abs(lt - r.choice_time) <= 0.005).any())))
mis = pd.DataFrame(mis)
print(f"mismatched trials {len(mis)}: choice_time within 5 ms after the window opening on {(mis.c_minus_w.between(0, 0.005)).mean():.3f}; "
      f"a lick in the 100 ms before the opening on {mis.lick_before.mean():.3f}; any lick within 5 ms of choice_time on {mis.any_lick_at_choice.mean():.3f}")
print(f"  choice_time minus window opening on mismatched trials: median {mis.c_minus_w.median() * 1000:.1f} ms, "
      f"95th percentile {mis.c_minus_w.quantile(0.95) * 1000:.1f} ms")
nsw = np.concatenate([(d.choice_time - d.response_window_open_time).to_numpy()
                      for d in (s["trials"][s["trials"].outcome != "spoil"] for s in S)])
print(f"all non-spoiled trials {len(nsw)}: choice_time 0 to 5 ms after the window opening on {int(((nsw >= 0) & (nsw <= 0.005)).sum())} "
      f"({((nsw >= 0) & (nsw <= 0.005)).mean():.3f}), of them {int(mis.c_minus_w.between(0, 0.005).sum())} mismatched")
per = ns.groupby("path").apply(lambda d: (d.ok_time & d.ok_side).mean(), include_groups=False)
print(f"per session, non-spoiled: lowest {per.min():.4f}, sessions below 0.99 {int((per < 0.99).sum())}")
f_ns, f_sp = (ns.ok_time & ns.ok_side).mean(), (sp.ok_time & sp.ok_side).mean() if len(sp) else 1.0
print("criterion:", "reproduced" if min(f_ns, f_sp) >= 0.99 else "partially_reproduced" if min(f_ns, f_sp) >= 0.95 else "not_reproduced")

non-spoiled trials: 32085; first lick at or after the window within 5 ms of choice_time: 0.7761; side matches choice: 0.9906; both: 0.7761
  variant, first lick from 5 ms before the window: 0.7800
  first lick minus choice_time: median -1.00 ms, 1st to 99th percentile -3.0 to 1767.4 ms; no lick after the window on 0 trials
spoiled trials: 265; choice_time 45 s after the window (within 5 ms): 0.0000; choice nogo: 1.0000; choice_time NaN on 265 of them

mismatched trials 7185: choice_time within 5 ms after the window opening on 1.000; a lick in the 100 ms before the opening on 0.848; any lick within 5 ms of choice_time on 0.057
  choice_time minus window opening on mismatched trials: median 1.0 ms, 95th percentile 2.0 ms
all non-spoiled trials 32085: choice_time 0 to 5 ms after the window opening on 7848 (0.245), of them 7183 mismatched
per session, non-spoiled: lowest 0.5245, sessions below 0.99 115
criterion: not_reproduced

exp_out = np.where(T.choice == "nogo", "spoil", np.where(T.choice == T.rewarded_side, "correct", "error"))
T["outcome_ok"] = T.outcome.to_numpy() == exp_out
for name, m in [("all trials", np.ones(len(T), bool)), ("direct delivery", T.direct_delivery),
                ("not direct delivery", ~T.direct_delivery)]:
    sub = T[m]
    print(f"{name}: {sub.outcome_ok.sum()} of {len(sub)} consistent ({sub.outcome_ok.mean():.4f})")
print(pd.crosstab([T.outcome], [T.choice == T.rewarded_side, T.choice == "nogo"]).to_string())
f = T.outcome_ok.mean()
print("criterion:", "reproduced" if f == 1 else "partially_reproduced" if f >= 0.99 else "not_reproduced")

all trials: 32350 of 32350 consistent (1.0000)
direct delivery: 46 of 46 consistent (1.0000)
not direct delivery: 32304 of 32304 consistent (1.0000)
col_0   False        True 
choice  False True   False
outcome                   
correct     0     0  24962
error    7123     0      0
spoil       0   265      0
criterion: reproduced

rows = []
for s in S:
    rw = s["events"]["rewards"]
    for r in s["trials"][~s["trials"].ignore_trial].itertuples():
        t = rw[r.rewarded_side]
        k = ((t >= r.choice_time - 0.005) & (t < r.stop_time)).sum()
        rows.append(dict(path=s["path"], outcome=r.outcome, side=r.rewarded_side, rewarded=k > 0,
                         after=r.stop_time - r.choice_time))
RW = pd.DataFrame(rows)
for side in ["left", "right", "both sides"]:
    sub = RW if side == "both sides" else RW[RW.side == side]
    c, e = sub[sub.outcome == "correct"], sub[sub.outcome == "error"]
    print(f"rewarded side {side}: reward on {c.rewarded.mean():.4f} of {len(c)} correct trials, "
          f"on {e.rewarded.mean():.4f} of {len(e)} error trials")
c, e = RW[RW.outcome == "correct"], RW[RW.outcome == "error"]
ra = e.after.median()
rb = ra - c.after.median()
print(f"timeout (a) median stop_time - choice_time on error trials: {ra:.3f} s (expected 9 +- 0.9)")
print(f"timeout (b) error median minus correct median ({c.after.median():.3f} s): {rb:.3f} s")
print(f"  error trials, 25th to 75th percentile {e.after.quantile(0.25):.3f} to {e.after.quantile(0.75):.3f} s")
per = RW.groupby(["path", "outcome"]).after.median().unstack()
print(f"  per session, error minus correct median: median {(per.error - per.correct).median():.3f} s, "
      f"range {(per.error - per.correct).min():.3f} to {(per.error - per.correct).max():.3f} s")
rew_ok = c.rewarded.mean() == 1 and e.rewarded.mean() == 0
rew_part = c.rewarded.mean() >= 0.95 and e.rewarded.mean() <= 0.05
ta, tb = abs(ra - 9) <= 0.9, abs(rb - 9) <= 0.9
print(f"reward part holds: {rew_ok} (partially: {rew_part}); timeout (a) holds: {ta}; (b) holds: {tb}")
print("criterion:", "reproduced" if rew_ok and ta and tb else
      "not_reproduced" if not ta and not tb else "not_reproduced (reward part)" if ta and tb else
      "blocked (timeout readings disagree)" + ("" if rew_ok else
      "; the reward part also fails as computed, on the right side only, where the reward series is the right lick series"))

rewarded side left: reward on 1.0000 of 8917 correct trials, on 0.0000 of 2497 error trials
rewarded side right: reward on 0.9990 of 8698 correct trials, on 0.1664 of 2633 error trials
rewarded side both sides: reward on 0.9995 of 17615 correct trials, on 0.0854 of 5130 error trials
timeout (a) median stop_time - choice_time on error trials: 12.212 s (expected 9 +- 0.9)
timeout (b) error median minus correct median (3.876 s): 8.336 s
  error trials, 25th to 75th percentile 12.072 to 12.387 s
  per session, error minus correct median: median 8.377 s, range 7.044 to 8.535 s
reward part holds: False (partially: False); timeout (a) holds: False; (b) holds: True
criterion: blocked (timeout readings disagree); the reward part also fails as computed, on the right side only, where the reward series is the right lick series

rows = []
for s in S:
    tr = s["trials"]
    early = tr[tr.trial <= 44]
    first_random = tr.loc[tr.stim_is_random, "trial"].min() if tr.stim_is_random.any() else np.nan
    rows.append(dict(path=s["path"], first_trial=tr.trial.min(), n_early=len(early),
                     n_early_random=int(early.stim_is_random.sum()), first_random=first_random))
FA = pd.DataFrame(rows)
has_early = FA[FA.n_early > 0]
ok = has_early.n_early_random == 0
print(f"sessions with trials numbered 0 to 44 present: {len(has_early)} of {len(FA)}; "
      f"first trial number in the table: {FA.first_trial.min()} to {FA.first_trial.max()}")
print(f"sessions where every such trial is non-random: {int(ok.sum())} of {len(has_early)}")
print(f"random trials among trials 0 to 44: {has_early.n_early_random.sum()} of {has_early.n_early.sum()}")
print(f"first random trial number per session: {FA.first_random.value_counts().sort_index().to_dict()}")
for r in FA[FA.n_early == 0].itertuples():
    print(f"  no trials 0 to 44: {r.path.split('/')[-1]}: first trial in the table {r.first_trial}, first random trial {r.first_random}")
print("criterion:", "reproduced" if ok.all() else "partially_reproduced" if ok.mean() > 0.5 else "not_reproduced")

sessions with trials numbered 0 to 44 present: 110 of 115; first trial number in the table: 0 to 124
sessions where every such trial is non-random: 110 of 110
random trials among trials 0 to 44: 0 of 4229
first random trial number per session: {45: 110, 61: 1, 76: 1, 88: 1, 102: 1, 124: 1}
  no trials 0 to 44: sub-219CR_ses-20190409T144717_behavior+image.nwb: first trial in the table 61, first random trial 61
  no trials 0 to 44: sub-KF119_ses-20180110T161350_behavior+image.nwb: first trial in the table 124, first random trial 124
  no trials 0 to 44: sub-KF132_ses-20180220T135608_behavior+image.nwb: first trial in the table 76, first random trial 76
  no trials 0 to 44: sub-KF134_ses-20180211T154918_behavior+ecephys+image.nwb: first trial in the table 88, first random trial 88
  no trials 0 to 44: sub-KM101_ses-20170714T132037_behavior+image.nwb: first trial in the table 102, first random trial 102
criterion: reproduced

flags = {"direct_delivery": T.direct_delivery, "stim_is_random False": ~T.stim_is_random,
         "outcome spoil": T.outcome == "spoil", "optogenetic": T.optogenetic}
any_flag = np.zeros(len(T), bool)
for k, m in flags.items():
    m = m.to_numpy()
    any_flag |= m
    print(f"{k}: {m.sum()} trials, ignore_trial True on {T.ignore_trial[m].sum()}, False on {(~T.ignore_trial[m]).sum()}")
viol = any_flag & ~T.ignore_trial.to_numpy()
print(f"flagged trials with ignore_trial False: {viol.sum()} of {any_flag.sum()}, in {T.path[viol].nunique()} sessions")
print(f"ignore_trial True with none of the four flags (munged): {(T.ignore_trial.to_numpy() & ~any_flag).sum()}")
print("ignore_trial equals (stim_is_random False or optogenetic) on every trial:",
      bool((T.ignore_trial == (~T.stim_is_random | T.optogenetic)).all()))
print("flag combinations among the violating trials:")
print(T[viol].groupby(["direct_delivery", "stim_is_random", "outcome", "optogenetic"]).size().to_string())
print("criterion:", "reproduced" if viol.sum() == 0 else "not_reproduced")

direct_delivery: 46 trials, ignore_trial True on 46, False on 0
stim_is_random False: 7574 trials, ignore_trial True on 7574, False on 0
outcome spoil: 265 trials, ignore_trial True on 105, False on 160
optogenetic: 2482 trials, ignore_trial True on 2482, False on 0
flagged trials with ignore_trial False: 160 of 9605, in 50 sessions
ignore_trial True with none of the four flags (munged): 0
ignore_trial equals (stim_is_random False or optogenetic) on every trial: True
flag combinations among the violating trials:
direct_delivery  stim_is_random  outcome  optogenetic
False            True            spoil    False          160
criterion: not_reproduced

n_pairs = n_bad = n_id = n_bad_gap = n_skip = 0
worst = 0.0
for s in S:
    tr = s["trials"]
    d = tr.start_time.to_numpy()[1:] - tr.stop_time.to_numpy()[:-1]
    gap = np.diff(tr.trial.to_numpy()) != 1
    n_pairs += len(d)
    bad = np.abs(d) > 0.001
    n_skip += int(gap.sum())
    n_bad += int(bad.sum())
    n_bad_gap += int((bad & gap).sum())
    worst = max(worst, float(np.abs(d).max()))
    n_id += int((tr.id != tr.trial).sum())
    if bad.any():
        print(f"  {s['path']}: {bad.sum()} pairs, next start minus stop {np.round(d[bad], 3).tolist()} s, at a missing trial number: {(bad & gap).sum()}")
print(f"consecutive pairs {n_pairs}; stop_time differs from the next start_time by more than 1 ms: {n_bad} "
      f"({n_bad_gap} of them where trial numbers skip); largest {worst * 1000:.3f} ms")
print(f"rows whose id differs from trial: {n_id} of {len(T)}")
print(f"trial numbers consecutive within each session: {sum(bool((np.diff(s['trials'].trial) == 1).all()) for s in S)} of {len(S)}")
print(f"variant, rows whose trial numbers are consecutive: {n_pairs - n_skip} pairs, {n_bad - n_bad_gap} differ by more than 1 ms")
print("criterion:", "reproduced" if n_bad == 0 and n_id == 0 else "not_reproduced")

  sub-219CR/sub-219CR_ses-20190404T152625_behavior+image.nwb: 1 pairs, next start minus stop [7.02] s, at a missing trial number: 1
  sub-KF119/sub-KF119_ses-20180112T160910_behavior+image.nwb: 1 pairs, next start minus stop [7.36] s, at a missing trial number: 1
  sub-KF132/sub-KF132_ses-20180220T135608_behavior+image.nwb: 2 pairs, next start minus stop [15.845, 15.205] s, at a missing trial number: 2
  sub-KF132/sub-KF132_ses-20180221T155455_behavior+image.nwb: 1 pairs, next start minus stop [68.045] s, at a missing trial number: 1
consecutive pairs 32235; stop_time differs from the next start_time by more than 1 ms: 5 (5 of them where trial numbers skip); largest 68045.000 ms
rows whose id differs from trial: 0 of 32350
trial numbers consecutive within each session: 111 of 115
variant, rows whose trial numbers are consecutive: 32230 pairs, 0 differ by more than 1 ms
criterion: not_reproduced

WW = pd.DataFrame([dict(path=s["path"], **{k: v for k, v in d.items() if np.isscalar(v)})
                   for s in S for w, d in s["whiskers"].items()])
print(f"processed: {len(WW)} whisker series; largest distance from the grid {WW.grid_offset.max():.2e} s; "
      f"steps shorter than a frame {WW.n_steps_short.sum()}; steps of one frame "
      f"{WW.n_steps_frame.sum()} of {WW.n_steps.sum()}; gaps {WW.n_gaps.sum()}")
print(f"  the seven series of a whisker share length and ends: {bool(WW.processed_lengths_equal.all() and WW.processed_ends_equal.all())}")
print(f"pose: largest distance from the grid {WW.pose_grid_offset.max():.2e} s; steps shorter than a frame "
      f"{WW.pose_n_steps_short.sum()}; steps other than one frame {WW.pose_n_steps_not_frame.sum()} of {WW.pose_n_steps.sum()}")
print(f"contacts: largest distance of start or stop from the grid {WW.contact_grid_offset.max():.2e} s; "
      f"frame_start x 0.005 against start_time, largest difference {WW.contact_frame_vs_time.max():.2e} s")
ok = (WW.grid_offset.max() <= 1e-6 and WW.pose_grid_offset.max() <= 1e-6 and WW.n_steps_short.sum() == 0
      and WW.pose_n_steps_short.sum() == 0)
print(f"first processed time per session: {WW.groupby('path').t_first.min().value_counts().to_dict()}")
print("criterion:", "reproduced" if ok else "not_reproduced")

processed: 460 whisker series; largest distance from the grid 9.09e-13 s; steps shorter than a frame 0; steps of one frame 252019206 of 252170504; gaps 151298
  the seven series of a whisker share length and ends: True
pose: largest distance from the grid 9.09e-13 s; steps shorter than a frame 0; steps other than one frame 397036 of 251801577
contacts: largest distance of start or stop from the grid 9.09e-13 s; frame_start x 0.005 against start_time, largest difference 9.09e-13 s

first processed time per session: {0.0: 114, 0.03: 1}
criterion: reproduced

print(f"(1) longest run of interpolated frames: {WW.max_interp_run.max()}; series with a run over 4: {(WW.max_interp_run > 4).sum()}")
rows = []
for s in S:
    for w, d in s["whiskers"].items():
        g = d["trial_gap"]
        rows.append(dict(path=s["path"], whisker=w, n=len(g), with_gap=np.isfinite(g).mean(),
                         within=np.mean(np.abs(g - 0.133) <= 0.0133), median=np.nanmedian(g) if np.isfinite(g).any() else np.nan))
G = pd.DataFrame(rows)
print("(2) per whisker, over sessions: fraction of trials whose blackout gap is within 10% of 133 ms")
print(G.groupby("whisker").agg(sessions=("path", "size"), lowest=("within", "min"), median=("within", "median"),
                               sessions_over_half=("within", lambda x: int((x > 0.5).sum())),
                               median_gap_s=("median", "median")).round(3).to_string())
allg = np.concatenate([d["trial_gap"] for s in S for w, d in s["whiskers"].items() if w != "C0"])
allg = allg[np.isfinite(allg)]
print(f"  C1 to C3 pooled: gap per trial median {np.median(allg):.4f} s, 5th to 95th percentile "
      f"{np.quantile(allg, 0.05):.4f} to {np.quantile(allg, 0.95):.4f} s; within 10%: {np.mean(np.abs(allg - 0.133) <= 0.0133):.3f}")
vals, cnt = np.unique(np.round(allg, 3), return_counts=True)
print("  most common gap values (s):", dict(zip(vals[np.argsort(-cnt)][:6], np.sort(cnt)[::-1][:6])))
p1 = WW.max_interp_run.max() <= 4
p2 = bool((G[G.whisker != "C0"].within > 0.5).all())
print(f"part (1) holds: {p1}; part (2) holds in every session for C1 to C3: {p2}")
print("criterion:", "reproduced" if p1 and p2 else "partially_reproduced" if p1 or p2 else "not_reproduced")
fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(allg * 1000, bins=np.arange(0, 400, 5))
ax.axvspan(119.7, 146.3, color="orange", alpha=0.3, label="133 ms +- 10%")
ax.set_xlabel("gap at trial start (ms, missing time)"); ax.set_ylabel("trials, C1 to C3"); ax.legend()
print(f"plotted range 0 to 400 ms; {np.sum(allg * 1000 >= 400)} values above it")
plt.show()

(1) longest run of interpolated frames: 4; series with a run over 4: 0
(2) per whisker, over sessions: fraction of trials whose blackout gap is within 10% of 133 ms
         sessions  lowest  median  sessions_over_half  median_gap_s
whisker                                                            
C0            115   0.903     1.0                 115         0.135
C1            115   0.984     1.0                 115         0.135
C2            115   0.994     1.0                 115         0.135
C3            115   0.987     1.0                 115         0.135
  C1 to C3 pooled: gap per trial median 0.1350 s, 5th to 95th percentile 0.1350 to 0.1400 s; within 10%: 0.999
  most common gap values (s): {np.float64(0.135): np.int64(89298), np.float64(0.14): np.int64(6959), np.float64(0.13): np.int64(715), np.float64(0.145): np.int64(22), np.float64(0.15): np.int64(10), np.float64(0.155): np.int64(6)}
part (1) holds: True; part (2) holds in every session for C1 to C3: True
criterion: reproduced
plotted range 0 to 400 ms; 3 values above it

print(f"interpolated frames: {WW.n_interpolated.sum()}; with a bend value that is not NaN: {WW.bend_notnan_on_interp.sum()}")
print(f"non-interpolated frames with bend NaN: {WW.bend_nan_on_noninterp.sum()} in "
      f"{(WW.bend_nan_on_noninterp > 0).sum()} whisker series")
print("criterion:", "reproduced" if WW.bend_notnan_on_interp.sum() == 0 else "not_reproduced")

interpolated frames: 368927; with a bend value that is not NaN: 0
non-interpolated frames with bend NaN: 0 in 0 whisker series
criterion: reproduced

NT = pd.DataFrame([dict(path=e["path"], rate=e["ecephys"]["rate"], has_ts=e["ecephys"]["has_timestamps"],
                        rate_span=(e["ecephys"]["n_samples"] - 1) / (e["ecephys"]["ts_last"] - e["ecephys"]["ts_first"]))
                   for e in E])
NT["dev_pct"] = 100 * (NT.rate - 30000) / 30000
print(f"sessions {len(NT)}; with a timestamps array {NT.has_ts.sum()}")
print(f"rate {NT.rate.min():.4f} to {NT.rate.max():.4f} Hz; deviation from 30 kHz {NT.dev_pct.min():.5f}% to {NT.dev_pct.max():.5f}%")
print(f"sessions at exactly 30000 Hz: {(NT.rate == 30000).sum()}; rate from span equals stored rate (1e-6 Hz): {bool((abs(NT.rate_span - NT.rate) < 1e-6).all())}")
ok = (NT.dev_pct.abs() < 1).all() and (NT.rate != 30000).any()
print("criterion:", "reproduced" if ok else "not_reproduced")

sessions 44; with a timestamps array 0
rate 29999.9354 to 30000.0199 Hz; deviation from 30 kHz -0.00022% to 0.00007%
sessions at exactly 30000 Hz: 0; rate from span equals stored rate (1e-6 Hz): True
criterion: reproduced

idx = [a for a in assets if a["path"].endswith("Continuous_Data.openephys")]
listed = collections.Counter()
for a in idx:
    txt = requests.get(R.s3_url(a["asset_id"]), timeout=60).text
    for name in re.findall(r'filename="100_([A-Z]+)\d+\.continuous"', txt):
        listed[name] += 1
present = collections.Counter(re.sub(r"\d+$", "", a["path"].rsplit("_", 1)[1].split(".")[0])
                              for a in assets if a["path"].endswith(".continuous"))
print(f"channel files listed in the {len(idx)} .openephys indexes, by kind: {dict(listed)}")
print(f"channel files deposited, by kind: {dict(present)}")
ev_sizes = collections.Counter(a["size"] for a in assets if a["path"].endswith("all_channels.events"))
print(f"all_channels.events sizes (bytes; the header alone is 1024): {dict(ev_sizes)}")

channel files listed in the 44 .openephys indexes, by kind: {'CH': 2880, 'AUX': 270, 'ADC': 360}
channel files deposited, by kind: {'CH': 2816}
all_channels.events sizes (bytes; the header alone is 1024): {1024: 44}

kinds = collections.Counter(a["path"].rsplit(".", 1)[-1] for a in assets)
print("file kinds in the dandiset:", dict(kinds))
print("trials columns:", list(T.columns[:15]))
logs = [a["path"] for a in assets if a["path"].endswith(".log")]
print(f".log files: {len(logs)}, in {len({p.rsplit('/', 1)[0] for p in logs})} folders, all under sourcedata/: "
      f"{all(p.startswith('sourcedata/') for p in logs)}; e.g. {logs[:2]}")
print("  names:", dict(collections.Counter(p.rsplit("/", 1)[1] for p in logs)))

file kinds in the dandiset: {'nwb': 115, 'mkv': 115, 'continuous': 2816, 'openephys': 44, 'events': 88, 'xml': 2, 'npy': 660, 'mat': 88, 'tsv': 80, 'm': 88, 'py': 44, 'log': 44, 'dat': 44}
trials columns: ['id', 'choice', 'choice_time', 'direct_delivery', 'ignore_trial', 'optogenetic', 'outcome', 'response_window_open_time', 'rewarded_side', 'servo_position', 'start_time', 'stim_is_random', 'stimulus', 'stop_time', 'trial']
.log files: 44, in 44 folders, all under sourcedata/: True; e.g. ['sourcedata/sub-KM131/sub-KM131_ses-20180116T184757_behavior+ecephys+image/sort20180116202010/phy.log', 'sourcedata/sub-KM131/sub-KM131_ses-20180117T194023_behavior+ecephys+image/sort20180509034311/phy.log']
  names: {'phy.log': 44}

rows = []
for e in E:
    tr = e["trials"]
    for k, u in enumerate(e["ecephys"]["units"]):
        a, b = u["intervals"][0]
        rows.append(dict(path=e["path"], subject=e["subject"], unit=int(e["ecephys"]["unit_id"][k]), n_spikes=len(u["spikes"]),
                         isi=R.isi_fraction(u["spikes"]), low60=R.low_rate_fraction_bins(u["spikes"], a, b, 60.0),
                         lowtrial=R.low_rate_fraction_trials(u["spikes"], a, b, tr)))
U = pd.DataFrame(rows)
print(f"units {len(U)} in {U.path.nunique()} sessions; spikes per unit {U.n_spikes.min()} to {U.n_spikes.max()}")
print(f"ISI below 3 ms: median {U.isi.median():.4f}, max {U.isi.max():.4f}; units above 5%: {(U.isi > 0.05).sum()}")
for c, name in [("low60", "(a) 60 s bins"), ("lowtrial", "(b) trial bins")]:
    print(f"{name}: fraction of bins below 5% of mean rate, median {U[c].median():.4f}, max {U[c].max():.4f}; "
          f"units above 20%: {(U[c] > 0.2).sum()}; units with no bin: {U[c].isna().sum()}")
print("units over a threshold, by subject:")
print(U.assign(isi_bad=U.isi > 0.05, a_bad=U.low60 > 0.2, b_bad=U.lowtrial > 0.2)
      .groupby("subject")[["isi_bad", "a_bad", "b_bad"]].sum().join(U.groupby("subject").size().rename("units")).to_string())
for thr in [0.05, 0.1, 0.2, 0.3]:
    print(f"  firing-rate threshold {thr:.2f}: units above, (a) {(U.low60 > thr).sum()}, (b) {(U.lowtrial > thr).sum()}")
isi_ok = (U.isi <= 0.05).all()
fa, fb = (U.low60.fillna(0) <= 0.2).all(), (U.lowtrial.fillna(0) <= 0.2).all()
print(f"ISI part holds: {isi_ok}; firing-rate part (a): {fa}, (b): {fb}")
kf = U.path.str.contains("sub-KF134_ses-20180209")
print(f"variant without KF134 2018-02-09 ({kf.sum()} units): units above 5% ISI {(U[~kf].isi > 0.05).sum()}, "
      f"above 20% (a) {(U[~kf].low60 > 0.2).sum()}, (b) {(U[~kf].lowtrial > 0.2).sum()}")
print("criterion:", "reproduced" if isi_ok and fa and fb else "not_reproduced" if (not isi_ok) or (not fa and not fb)
      else "blocked (readings disagree)")

# amplitude drift from the Kilosort outputs (outside the outcome)
sort_files = collections.defaultdict(dict)
for a in assets:
    parts = a["path"].split("/")
    if a["path"].startswith("sourcedata/") and len(parts) == 5 and parts[3].startswith("sort"):
        sort_files["/".join(parts[1:3]) + ".nwb", parts[3]][parts[4]] = a
print(f"sort folders: {len(sort_files)} for {len({k[0] for k in sort_files})} ecephys sessions")
jobs = [(sort_files[k], e["ecephys"]["unit_id"], [u["spikes"] for u in e["ecephys"]["units"]], e["ecephys"]["rate"])
        for e in E for k in sort_files if k[0] == e["path"]]
AD = pd.DataFrame([r for res in R.amplitude_drifts(jobs) for r in res])
print(f"units matched by cluster id with equal spike counts: {(AD.n_ks == AD.n_nwb).sum()} of {len(AD)}; "
      f"largest spread of NWB time minus Kilosort sample / rate within a unit: {AD.offset_spread.max():.2e} s")
print(f"units not matched: {(AD.n_ks != AD.n_nwb).sum()}, of which no Kilosort spike with the id {(AD.n_ks == 0).sum()}, "
      f"more Kilosort spikes {(AD.n_ks > AD.n_nwb).sum()}, fewer {((AD.n_ks > 0) & (AD.n_ks < AD.n_nwb)).sum()}")
ad = AD.drift_pct_per_min.abs()
print(f"amplitude change per minute (|slope| / mean, %): median {ad.median():.3f}, 95th percentile {ad.quantile(0.95):.3f}, "
      f"max {ad.max():.3f}; units above 1.5%: {(ad > 1.5).sum()} of {ad.notna().sum()}")

units 999 in 44 sessions; spikes per unit 745 to 107009
ISI below 3 ms: median 0.0066, max 0.0493; units above 5%: 0
(a) 60 s bins: fraction of bins below 5% of mean rate, median 0.0000, max 0.2326; units above 20%: 6; units with no bin: 0
(b) trial bins: fraction of bins below 5% of mean rate, median 0.0444, max 0.4015; units above 20%: 78; units with no bin: 0
units over a threshold, by subject:
         isi_bad  a_bad  b_bad  units
subject                              
221CR          0      0      6     66
228CR          0      1      7     81
229CR          0      1     18    168
231CR          0      1      6     64
245CR          0      0      5     55
255CR          0      0      3     19
267CR          0      1      7    148
KF119          0      0      2     20
KF132          0      0      7     92
KF134          0      2     13    144
KM131          0      0      4    142
  firing-rate threshold 0.05: units above, (a) 306, (b) 477
  firing-rate threshold 0.10: units above, (a) 180, (b) 319
  firing-rate threshold 0.20: units above, (a) 6, (b) 78
  firing-rate threshold 0.30: units above, (a) 0, (b) 6
ISI part holds: True; firing-rate part (a): False, (b): False
variant without KF134 2018-02-09 (16 units): units above 5% ISI 0, above 20% (a) 6, (b) 78
criterion: not_reproduced
sort folders: 44 for 44 ecephys sessions

units matched by cluster id with equal spike counts: 893 of 999; largest spread of NWB time minus Kilosort sample / rate within a unit: 2.27e-12 s
units not matched: 106, of which no Kilosort spike with the id 0, more Kilosort spikes 106, fewer 0
amplitude change per minute (|slope| / mean, %): median 0.328, 95th percentile 1.011, max 1.471; units above 1.5%: 0 of 893

rows = []
for e in E:
    x = e["ecephys"]
    for k in range(len(x["unit_id"])):
        ch = int(x["unit_channel"][k])
        rows.append(dict(path=e["path"], depth=x["unit_depth"][k], layer=x["unit_layer"][k], ch=ch,
                         el_y=x["el_y"][ch] if 0 <= ch < len(x["el_y"]) else np.nan,
                         el_layer=re.search(r"layer L?(\S+),", x["el_location"][ch]).group(1) if 0 <= ch < len(x["el_y"]) else None))
UD = pd.DataFrame(rows)
UD["depth_ok"] = np.isclose(UD.depth, UD.el_y)
print(f"units {len(UD)}; depth equals the y of the electrodes row neuron_channel: {UD.depth_ok.sum()}")
print(f"electrode ids equal their row index in every session: {all((e['ecephys']['el_id'] == np.arange(len(e['ecephys']['el_id']))).all() for e in E)}")
rng = UD.groupby("layer").depth.agg(["size", "min", "max"]).sort_values("min")
print(rng.to_string())
# the files name layers "1", "2/3", "4", "5a", "5b", "6", "6b": L1 to L6 with sublayers of L5 and L6
order = ["1", "2/3", "4", "5a", "5b", "6", "6b"]
present = [l for l in order if l in rng.index]
mono = all(rng.loc[a, "max"] <= rng.loc[b, "min"] for a, b in zip(present, present[1:]))
print(f"layers present {list(rng.index)}; in order {present}; ranges do not overlap and run superficial to deep: {mono}")
print(f"unit layer equals the layer in its electrode's location string: {(UD.layer == UD.el_layer).sum()} of {len(UD)} as spelled, "
      f"{(UD.layer.str.lower() == UD.el_layer.str.lower()).sum()} ignoring case (the strings write 5A, 5B, 6B where the units write 5a, 5b, 6b)")
print("unit layer (rows) against the layer in its electrode's location string (columns):")
print(pd.crosstab(UD.layer, UD.el_layer).to_string())
print("depth ranges of the location-string layers over all electrodes of all sessions (um):")
EL = pd.DataFrame([dict(y=y, el_layer=re.search(r"layer L?(\S+),", loc).group(1)) for e in E
                   for y, loc in zip(e["ecephys"]["el_y"], e["ecephys"]["el_location"])])
print(EL.groupby("el_layer").y.agg(["size", "min", "max"]).sort_values("min").to_string())
for lab in sorted(set(zip(UD.layer, UD.el_layer)) - {(l, l) for l in UD.layer}, key=str)[:3]:
    r = UD[(UD.layer == lab[0]) & (UD.el_layer == lab[1])].iloc[0]
    print(f"  example, unit layer {r.layer}, depth {r.depth:.0f} um: {by_path[r.path]['ecephys']['el_location'][r.ch][:60]}...")
print("criterion:", "reproduced" if UD.depth_ok.all() and mono and set(rng.index) <= set(order) else "not_reproduced")

units 999; depth equals the y of the electrodes row neuron_channel: 999
electrode ids equal their row index in every session: True
       size      min     max
layer                       
1        12    10.00   126.9
2/3     114   130.20   416.6
4       269   421.50   625.2
5a      186   626.80   732.5
5b      249   733.50  1005.0
6       158  1008.05  1364.3
6b       11  1371.60  1424.3
layers present ['1', '2/3', '4', '5a', '5b', '6', '6b']; in order ['1', '2/3', '4', '5a', '5b', '6', '6b']; ranges do not overlap and run superficial to deep: True
unit layer equals the layer in its electrode's location string: 553 of 999 as spelled, 999 ignoring case (the strings write 5A, 5B, 6B where the units write 5a, 5b, 6b)
unit layer (rows) against the layer in its electrode's location string (columns):
el_layer   1  2/3    4   5A   5B    6  6B
layer                                    
1         12    0    0    0    0    0   0
2/3        0  114    0    0    0    0   0
4          0    0  269    0    0    0   0
5a         0    0    0  186    0    0   0
5b         0    0    0    0  249    0   0
6          0    0    0    0    0  158   0
6b         0    0    0    0    0    0  11
depth ranges of the location-string layers over all electrodes of all sessions (um):
          size      min     max
el_layer                       
1          220   -53.20   127.3
2/3        627   128.00   418.8
4          448   419.05   625.5
5A         240   626.80   732.5
5B         584   733.00  1005.5
6          615  1006.80  1364.3
6B          13  1371.60  1437.5
  example, unit layer 5a, depth 629 um: barrel cortex, (approx) layer L5A, Intan channel 46, OpenEph...
  example, unit layer 5b, depth 789 um: barrel cortex, (approx) layer L5B, Intan channel 14, OpenEph...
  example, unit layer 6b, depth 1373 um: barrel cortex, (approx) layer L6B, Intan channel 30, OpenEph...
criterion: reproduced

rows = []
pat = re.compile(r"layer (\S+), Intan channel (\d+), OpenEphys channel (\d+).*channel (\d+) starting with 0 = most superficial and including.*channel (\d+) starting with 0 = most superficial and excluding")
for e in E:
    x = e["ecephys"]
    y = x["el_y"]
    dy = np.diff(y)
    parsed = [pat.search(l) for l in x["el_location"]]
    incl = np.array([int(p.group(4)) for p in parsed if p])
    excl = np.array([int(p.group(5)) for p in parsed if p])
    rows.append(dict(path=e["path"], n=len(y), n_data=x["n_channels"], increasing=bool((dy > 0).all()),
                     multiple20=bool(np.allclose(dy / 20, np.round(dy / 20))), steps_over_20=int((dy > 20.001).sum()),
                     region_ok=bool((x["electrodes_region"] == np.arange(len(y))).all()),
                     parsed=sum(p is not None for p in parsed), excl_ok=bool((excl == np.arange(len(y))).all()),
                     y_vs_incl=float(np.ptp(y - 20 * incl)) if len(incl) == len(y) else np.nan,
                     oe_intan=all(int(p.group(3)) == int(p.group(2)) + 1 for p in parsed if p),
                     span_full=(63 * 20)))
EO = pd.DataFrame(rows)
print(f"sessions {len(EO)}; rows {EO.n.min()} to {EO.n.max()}; rows equal data columns: {(EO.n == EO.n_data).sum()}")
print(f"y strictly increasing: {EO.increasing.sum()}; steps multiples of 20 um: {EO.multiple20.sum()}; "
      f"sessions with a larger step (a broken channel): {(EO.steps_over_20 > 0).sum()}")
print(f"ElectricalSeries electrodes region is 0..n-1: {EO.region_ok.sum()}")
print(f"location strings with all five items: {EO.parsed.sum()} of {EO.n.sum()}; 'excluding' index equals the row: {EO.excl_ok.sum()}; "
      f"OpenEphys = Intan + 1 in every string: {EO.oe_intan.all()}")
print(f"y minus 20 x ('including' index) constant within a session (range 0): {(EO.y_vs_incl == 0).sum()} of {len(EO)}; "
      f"so a full probe spans 63 x 20 = {63 * 20} um between end sites")
ok = (EO.increasing & EO.multiple20 & (EO.n <= 64) & (EO.n == EO.n_data) & (EO.parsed == EO.n)).all()
print("criterion:", "reproduced" if ok else "not_reproduced")

sessions 44; rows 56 to 64; rows equal data columns: 44
y strictly increasing: 44; steps multiples of 20 um: 44; sessions with a larger step (a broken channel): 26
ElectricalSeries electrodes region is 0..n-1: 44
location strings with all five items: 2747 of 2747; 'excluding' index equals the row: 44; OpenEphys = Intan + 1 in every string: True
y minus 20 x ('including' index) constant within a session (range 0): 44 of 44; so a full probe spans 63 x 20 = 1260 um between end sites
criterion: reproduced

locs = pd.Series({e["path"]: list(e["ecephys"]["groups"].values()) for e in E})
print(f"electrode groups per session: {locs.map(len).value_counts().to_dict()}")
col = locs.map(lambda l: re.match(r"^\s*([A-E]\d|alpha|beta|gamma|delta|[A-E]\d?/[A-E]?\d)", l[0]) if l else None)
print("group locations:", locs.map(lambda l: l[0]).value_counts().to_dict())
names = col.map(lambda m: m.group(1) if m else None)
print("column per session:", names.value_counts(dropna=False).to_dict())
print("criterion:", "reproduced" if names.notna().all() else "not_reproduced")

electrode groups per session: {1: 44}
group locations: {'C3 column of the S1 barrel field': 13, 'C2 column of the S1 barrel field': 10, 'C1 column of the S1 barrel field': 10, 'D2 column of the S1 barrel field': 3, 'C4 column of the S1 barrel field': 2, 'D1 column of the S1 barrel field': 2, 'B2 column of the S1 barrel field': 1, 'B3 column of the S1 barrel field': 1, 'E1 column of the S1 barrel field': 1, 'E2 column of the S1 barrel field': 1}
column per session: {'C3': 13, 'C2': 10, 'C1': 10, 'D2': 3, 'C4': 2, 'D1': 2, 'B2': 1, 'B3': 1, 'E1': 1, 'E2': 1}
criterion: reproduced

rows = []
for e in E:
    x = e["ecephys"]
    one = 1 / x["rate"]
    for k, u in enumerate(x["units"]):
        iv = u["intervals"]
        a, b = iv[0]
        sp = u["spikes"]
        rows.append(dict(path=e["path"], n_int=len(iv), inside=bool(((sp >= a) & (sp <= b)).all()),
                         d_start=a - x["ts_first"], d_stop=b - x["ts_last"], one=one,
                         within=(a >= x["ts_first"] - one * (1 + 1e-6)) and (b <= x["ts_last"] + one * (1 + 1e-6))))
OI = pd.DataFrame(rows)
OI["ends_match"] = (OI.d_start.abs() <= OI.one * (1 + 1e-6)) & (OI.d_stop.abs() <= OI.one * (1 + 1e-6))
print(f"units {len(OI)}; with exactly one interval {(OI.n_int == 1).sum()}; all spikes inside {OI.inside.sum()}")
print(f"interval start minus first sample time: {(OI.d_start / OI.one).min():.4f} to {(OI.d_start / OI.one).max():.4f} samples; "
      f"stop minus last sample time: {(OI.d_stop / OI.one).min():.4f} to {(OI.d_stop / OI.one).max():.4f} samples")
print(f"ends match within one sample: {OI.ends_match.sum()}; interval within the recording: {OI.within.sum()}")
print("criterion:", "reproduced" if (OI.n_int == 1).all() and OI.inside.all() and OI.ends_match.all() else
      "partially_reproduced" if (OI.n_int == 1).all() and OI.inside.all() and OI.within.all() else "not_reproduced")

units 999; with exactly one interval 999; all spikes inside 999
interval start minus first sample time: 0.0000 to 0.0000 samples; stop minus last sample time: 1.0000 to 1.0000 samples
ends match within one sample: 999; interval within the recording: 999
criterion: reproduced

rows = []
for s in S:
    ws = [d for d in s["whiskers"].values() if "pose_markers" in d]
    conv = sorted({c for d in ws for c in d["pose_conversion"]})
    rows.append(dict(path=s["path"], n_pose=len(ws), eight=all(len(d["pose_markers"]) == 8 for d in ws),
                     shape_ok=all(len({sh for sh in d["pose_shapes"]}) == 1 and d["pose_shapes"][0][1] == 2 for d in ws),
                     n_conv=len(conv), conv=conv[0] if conv else np.nan))
PS = pd.DataFrame(rows)
print(f"sessions with four PoseEstimation objects: {(PS.n_pose == 4).sum()} of {len(PS)}; eight series each: {PS.eight.sum()}; "
      f"series of one shape (N, 2) within each whisker: {PS.shape_ok.sum()}; one conversion value per session: {(PS.n_conv == 1).sum()}")
print(f"conversion (mm per pixel) across sessions: {PS.conv.round(6).value_counts().to_dict()}")
print("criterion:", "reproduced" if ((PS.n_pose == 4) & PS.eight & PS.shape_ok & (PS.n_conv == 1)).all() else "not_reproduced")

sessions with four PoseEstimation objects: 115 of 115; eight series each: 115; series of one shape (N, 2) within each whisker: 115; one conversion value per session: 115
conversion (mm per pixel) across sessions: {0.0466: 4, 0.04529: 4, 0.042167: 4, 0.041348: 3, 0.044697: 3, 0.046305: 3, 0.048176: 3, 0.0477: 2, 0.0425: 2, 0.043748: 2, 0.043731: 2, 0.045862: 2, 0.045901: 2, 0.046001: 2, 0.044876: 2, 0.04558: 2, 0.04517: 2, 0.045606: 2, 0.048237: 2, 0.049254: 2, 0.044126: 2, 0.042637: 2, 0.041538: 2, 0.048682: 1, 0.050718: 1, 0.050312: 1, 0.051162: 1, 0.049058: 1, 0.052039: 1, 0.052449: 1, 0.044975: 1, 0.043708: 1, 0.04107: 1, 0.040792: 1, 0.042516: 1, 0.044388: 1, 0.041612: 1, 0.045393: 1, 0.043097: 1, 0.042479: 1, 0.042809: 1, 0.043141: 1, 0.04446: 1, 0.044167: 1, 0.045405: 1, 0.046135: 1, 0.044191: 1, 0.044599: 1, 0.047051: 1, 0.045144: 1, 0.049923: 1, 0.047692: 1, 0.04877: 1, 0.050514: 1, 0.047963: 1, 0.044526: 1, 0.042453: 1, 0.041637: 1, 0.04239: 1, 0.042993: 1, 0.043296: 1, 0.050138: 1, 0.046965: 1, 0.044845: 1, 0.041692: 1, 0.042814: 1, 0.041998: 1, 0.044672: 1, 0.043656: 1, 0.045069: 1, 0.044352: 1, 0.047361: 1, 0.046887: 1, 0.046075: 1, 0.046712: 1, 0.046407: 1, 0.047115: 1, 0.040267: 1, 0.040461: 1, 0.0429: 1, 0.042249: 1, 0.047137: 1}
criterion: reproduced

rows = []
for s in S:
    tr = s["trials"].reset_index(drop=True)
    keep = ~tr.ignore_trial.to_numpy()
    for w, d in s["whiskers"].items():
        st = d.get("contact_start", np.array([]))
        i = R.assign_trials(st, tr)
        rows.append(dict(subject=s["subject"], whisker=w, n=int(((i >= 0) & keep[np.clip(i, 0, None)]).sum())))
CM = pd.DataFrame(rows).pivot_table(index="subject", columns="whisker", values="n", aggfunc="sum")
CM["C0_fraction"] = CM.C0 / CM[["C0", "C1", "C2", "C3"]].sum(axis=1)
CM["C0_fewest"] = (CM.C0 < CM.C1) & (CM.C0 < CM.C2) & (CM.C0 < CM.C3)
print(CM.round(3).to_string())
print(f"mice where C0 makes fewer contacts than each of C1 to C3: {CM.C0_fewest.sum()} of {len(CM)}")
print("criterion:", "reproduced" if CM.C0_fewest.mean() > 0.5 else "not_reproduced")

whisker    C0     C1     C2     C3  C0_fraction  C0_fewest
subject                                                   
219CR    4099   4834   8248   2952        0.204      False
221CR       3  14194  21929  11689        0.000       True
228CR       7   1025   1264    469        0.003       True
229CR      26   6068   6457   2625        0.002       True
231CR     502   2674  13171   7617        0.021       True
234CR    3142   2334    996      2        0.485      False
245CR    1144    746    539    227        0.431      False
255CR     221    266    149     71        0.313      False
267CR       1   2963   2896    659        0.000       True
KF119    6572   4415   2793    803        0.451      False
KF132     406   1446   2665   2400        0.059       True
KF134      67  10140  11694   6113        0.002       True
KM101      37   2254   3924   1503        0.005       True
KM131    1156   5073   4833   4826        0.073       True
KM91       47   6690   4294    229        0.004       True
mice where C0 makes fewer contacts than each of C1 to C3: 10 of 15
criterion: reproduced

cnt = T.groupby("path").agg(a=("ignore_trial", lambda x: int((~x).sum())))
cnt["b"] = T.assign(m=T.stim_is_random & (T.outcome != "spoil")).groupby("path").m.sum()
print(f"(a) ignore_trial False: smallest {cnt.a.min()}, sessions below 90: {(cnt.a < 90).sum()}")
print(f"(b) random and not spoiled: smallest {cnt.b.min()}, sessions below 90: {(cnt.b < 90).sum()}")
print("sessions below 90 under (a):", cnt[cnt.a < 90].to_dict("index"))
pa, pb = (cnt.a >= 90).all(), (cnt.b >= 90).all()
print("criterion:", "reproduced" if pa and pb else "not_reproduced" if (cnt[["a", "b"]].min(axis=1) < 90).any() and not pa and not pb
      else "blocked (readings disagree)")

(a) ignore_trial False: smallest 84, sessions below 90: 1
(b) random and not spoiled: smallest 90, sessions below 90: 0
sessions below 90 under (a): {'sub-228CR/sub-228CR_ses-20190715T182135_behavior+image.nwb': {'a': 84, 'b': 101}}
criterion: blocked (readings disagree)

nvid = V.set_index("nwb_name").frames
rows = []
for s in S:
    nf = nvid[s["path"][:-4]]
    for w, d in s["whiskers"].items():
        cm = d.get("contact_missing", np.array([], bool))
        rows.append(dict(path=s["path"], whisker=w, fa=1 - d["n_interpolated"] / d["n_processed"],
                         fb=d["n_pose"] / nf, n_contacts=len(cm), miss=cm.mean() if len(cm) else 0.0))
LQ = pd.DataFrame(rows)
res = {}
for wset, ws in [("i", ["C0", "C1", "C2", "C3"]), ("ii", ["C1", "C2", "C3"])]:
    sub = LQ[LQ.whisker.isin(ws)]
    for r in ["a", "b"]:
        fmin = sub["f" + r].min()
        cmax = sub.miss.max()
        res[(wset, r)] = (fmin >= 0.95) and (cmax <= 0.02)
        print(f"whiskers ({wset}), frames ({r}): smallest labeled fraction {fmin:.4f} "
              f"(sessions below 0.95: {sub[sub['f' + r] < 0.95].path.nunique()}); largest fraction of contacts with a missing frame "
              f"{cmax:.4f} (sessions above 0.02: {sub[sub.miss > 0.02].path.nunique()})")
print(LQ.groupby("whisker")[["fa", "fb", "miss"]].describe().loc[:, (slice(None), ["min", "50%", "max"])].round(4).to_string())
CT = pd.DataFrame([dict(path=s["path"], whisker=w, n=d["n_contacts"], incl=d["contact_duration_inclusive"],
                        excl=d["contact_duration_exclusive"], off=d["contact_stop_vs_frame_stop"],
                        miss=int(d["contact_missing"].sum()), miss_excl=int(d["contact_missing_excl"].sum()))
                   for s in S for w, d in s["whiskers"].items() if d.get("n_contacts")])
print(f"contacts {CT.n.sum()}: duration_frames = frame_stop - frame_start + 1 on {CT.incl.sum()}, "
      f"= frame_stop - frame_start on {CT.excl.sum()}; largest |stop_time - frame_stop x 0.005| {CT.off.max():.1e} s")
CT["frac_excl"] = CT.miss_excl / CT.n
print(f"variant, stop exclusive: largest fraction with a missing frame {CT.frac_excl.max():.4f} "
      f"(C1 to C3 {CT[CT.whisker != 'C0'].frac_excl.max():.4f}); sessions above 0.02: "
      f"{CT[CT.frac_excl > 0.02].path.nunique()} (C1 to C3 {CT[(CT.whisker != 'C0') & (CT.frac_excl > 0.02)].path.nunique()})")
for wset, ws in [("i", ["C0", "C1", "C2", "C3"]), ("ii", ["C1", "C2", "C3"])]:
    pool = CT[CT.whisker.isin(ws)].groupby("path")[["n", "miss", "miss_excl"]].sum()
    for col in ["miss", "miss_excl"]:
        fr = pool[col] / pool.n
        print(f"variant, contacts pooled over whiskers ({wset}), stop {'inclusive' if col == 'miss' else 'exclusive'}: "
              f"largest fraction {fr.max():.5f}, sessions above 0.02: {(fr > 0.02).sum()} of {len(fr)}")
        for path in fr[fr > 0.02].index:
            print(f"    {path}: {pool.loc[path, col]} of {pool.loc[path, 'n']} contacts")
print("readings holding:", {f"{k[0]}{k[1]}": v for k, v in res.items()})
print("criterion:", "reproduced" if all(res.values()) else "not_reproduced" if not any(res.values()) else "blocked (readings disagree)")

whiskers (i), frames (a): smallest labeled fraction 0.9309 (sessions below 0.95: 3); largest fraction of contacts with a missing frame 0.0769 (sessions above 0.02: 27)
whiskers (i), frames (b): smallest labeled fraction 0.8970 (sessions below 0.95: 7); largest fraction of contacts with a missing frame 0.0769 (sessions above 0.02: 27)
whiskers (ii), frames (a): smallest labeled fraction 0.9935 (sessions below 0.95: 0); largest fraction of contacts with a missing frame 0.0769 (sessions above 0.02: 16)
whiskers (ii), frames (b): smallest labeled fraction 0.9760 (sessions below 0.95: 0); largest fraction of contacts with a missing frame 0.0769 (sessions above 0.02: 16)
             fa                   fb                 miss                
            min     50%  max     min     50%     max  min     50%     max
whisker                                                                  
C0       0.9309  0.9996  1.0  0.8970  0.9849  0.9893  0.0  0.0000  0.0769
C1       0.9935  0.9994  1.0  0.9760  0.9848  0.9891  0.0  0.0012  0.0370
C2       0.9945  0.9999  1.0  0.9794  0.9854  0.9890  0.0  0.0000  0.0399
C3       0.9936  0.9998  1.0  0.9779  0.9852  0.9891  0.0  0.0000  0.0769

contacts 290345: duration_frames = frame_stop - frame_start + 1 on 290345, = frame_stop - frame_start on 0; largest |stop_time - frame_stop x 0.005| 9.1e-13 s
variant, stop exclusive: largest fraction with a missing frame 0.0769 (C1 to C3 0.0375); sessions above 0.02: 24 (C1 to C3 14)

variant, contacts pooled over whiskers (i), stop inclusive: largest fraction 0.02007, sessions above 0.02: 1 of 115
    sub-KM91/sub-KM91_ses-20161216T143049_behavior+image.nwb: 54 of 2691 contacts
variant, contacts pooled over whiskers (i), stop exclusive: largest fraction 0.01942, sessions above 0.02: 0 of 115
variant, contacts pooled over whiskers (ii), stop inclusive: largest fraction 0.02008, sessions above 0.02: 2 of 115
    sub-KM91/sub-KM91_ses-20161216T143049_behavior+image.nwb: 54 of 2689 contacts
    sub-KM91/sub-KM91_ses-20170118T170550_behavior+image.nwb: 70 of 3497 contacts
variant, contacts pooled over whiskers (ii), stop exclusive: largest fraction 0.01945, sessions above 0.02: 0 of 115
readings holding: {'ia': np.False_, 'ib': np.False_, 'iia': np.False_, 'iib': np.False_}
criterion: not_reproduced

rows = []
for s in S:
    for w, d in s["whiskers"].items():
        cols = d.get("contact_columns")
        rows.append(dict(path=s["path"], whisker=w, has=cols is not None, n=d.get("n_contacts", 0),
                         bad=d.get("contact_start_after_stop", 0),
                         has_loc=bool(cols) and {"tip_x", "tip_y"} <= set(cols), angle=bool(cols) and "angle" in cols,
                         cols=",".join(cols) if cols else ""))
CS = pd.DataFrame(rows)
print(f"contacts tables {CS.has.sum()} of {len(CS)} whisker-sessions; missing by whisker: {CS[~CS.has].whisker.value_counts().to_dict()}")
print(f"columns: {CS[CS.has].cols.value_counts().to_dict()}")
print(f"tables with no rows: {(CS.has & (CS.n == 0)).sum()}; rows with start after stop: {CS.bad.sum()} of {CS.n.sum()}")
h = CS[CS.has]
print("criterion:", "reproduced" if (h.has_loc.all() and h.angle.all() and h.bad.sum() == 0) else "not_reproduced")

contacts tables 431 of 460 whisker-sessions; missing by whisker: {'C0': 21, 'C3': 8}
columns: {'start_time,stop_time,base_x,base_y,tip_x,tip_y,angle,frame_start,frame_stop,duration_frames': 431}
tables with no rows: 0; rows with start after stop: 0 of 290345
criterion: reproduced

proc_names = collections.Counter(k for s in S for k in s["processing"])
print("processing modules:", dict(proc_names))
need_series = {"angle", "base_x", "base_y", "tip_x", "tip_y", "bend", "interpolated"}
rows = []
for s in S:
    b = set(s["processing"].get("behavior", []))
    ok = dict(
        licks=bool({"left", "right"} <= set(s["events"]["licks"])),
        rewards=bool({"left", "right"} <= set(s["events"]["rewards"])),
        pose=all(f"raw_position_whisker_{w}" in b for w in R.WHISKERS),
        processed=all(f"processed_position_whisker_{w}" in b and set(s["whiskers"][w]["processed_series"]) == need_series for w in R.WHISKERS),
        contacts_any=any(f"contacts_by_whisker_{w}" in b for w in R.WHISKERS),
        contacts_all4=all(f"contacts_by_whisker_{w}" in b for w in R.WHISKERS),
        trials=set(R.TRIAL_COLUMNS) <= set(s["trials"].columns),
        video=s.get("video_external_file") is not None,
        acquisition_ok=("behavioral video" in s["acquisition"]) and (("extracellular array recording" in s["acquisition"]) == ("ecephys" in s["types"])),
    )
    if s["ecephys"]:
        ok["units_cols"] = {"depth", "layer", "neuron_channel", "inhibitory", "spike_times", "obs_intervals"} <= set(s["ecephys"]["unit_columns"])
        x = s["ecephys"]
        ok["electrodes"] = len(x["el_y"]) == x["n_channels"]
        # Table 4: y, location and the electrode group (one per session); x and z left null
        ok["electrodes_table4"] = ({"y", "location", "group", "group_name"} <= set(x["el_columns"])
                                   and len(set(x["el_group_name"])) == 1 and len(x["groups"]) == 1
                                   and bool(x["el_x_all_nan"]) and bool(x["el_z_all_nan"]))
    ok["path"] = s["path"]
    rows.append(ok)
NS = pd.DataFrame(rows).set_index("path")
print("files holding each section (behavior naming):")
print(NS.sum().to_string())
print("extra trials columns:", sorted(set(T.columns) - set(R.TRIAL_COLUMNS) - {"id", "path", "subject", "task", "expected_side", "table6_stim", "follows", "outcome_ok"}))
# the external video file exists at the relative path
vid_paths = {a["path"] for a in mkv_assets}
rel_ok = sum((s["path"].rsplit("/", 1)[0] + "/" + s["video_external_file"][0]) in vid_paths for s in S)
print(f"behavioral video external_file resolves to a deposited video: {rel_ok} of {len(S)}; starting_time {set(s['video_starting_time'] for s in S)}, rate {set(s['video_rate'] for s in S)}")
core = ["licks", "rewards", "pose", "processed", "contacts_any", "trials", "video", "acquisition_ok"]
eph = ["units_cols", "electrodes", "electrodes_table4"]
ok_all = NS[core].all().all() and NS.loc[NS.units_cols.notna(), eph].all().all() and rel_ok == len(S)
print("electrodes table columns:", collections.Counter(tuple(e["ecephys"]["el_columns"]) for e in E))
print(f"files with all four contacts tables: {NS.contacts_all4.sum()} of {len(NS)}")
print("criterion:", "reproduced" if ok_all else "not_reproduced")

processing modules: {'behavior': 115}
files holding each section (behavior naming):
licks                115
rewards              115
pose                 115
processed            115
contacts_any         115
contacts_all4         86
trials               115
video                115
acquisition_ok       115
units_cols            44
electrodes            44
electrodes_table4     44
extra trials columns: []
behavioral video external_file resolves to a deposited video: 115 of 115; starting_time {0.0}, rate {200.0}
electrodes table columns: Counter({('filtering', 'group', 'group_name', 'id', 'imp', 'location', 'reference', 'x', 'y', 'z'): 44})
files with all four contacts tables: 86 of 115
criterion: reproduced

name_ok = sess.types.isin(["behavior+image", "behavior+ecephys+image"]) & (sess.folder == sess.subject)
print(f"NWB files in their subject's folder and named as described: {name_ok.sum()} of {len(sess)}; types {sess.types.value_counts().to_dict()}")
vid_by = V.groupby("nwb_name").size()
has_vid = sess.path.str[:-4].map(vid_by).fillna(0)
print(f"sessions with exactly one video in a folder named as the NWB file: {(has_vid == 1).sum()} of {len(sess)}")
src_dirs = {a["path"].split("/")[1] + "/" + a["path"].split("/")[2] for a in assets if a["path"].startswith("sourcedata/")}
src_files = collections.defaultdict(set)
for a in assets:
    if a["path"].startswith("sourcedata/"):
        p = a["path"].split("/")
        src_files[p[1] + "/" + p[2]].add(p[-1])
eph_names = set(sess.loc[sess.ecephys, "path"].str[:-4])
print(f"sourcedata session folders: {len(src_dirs)}; equal to the ecephys NWB names: {src_dirs == eph_names}")
oe = sum(any(f.endswith(".continuous") for f in v) and "Continuous_Data.openephys" in v for v in src_files.values())
ks = sum({"spike_times.npy", "spike_clusters.npy", "cluster_group.tsv"} <= v for v in src_files.values())
print(f"folders with Open Ephys files: {oe}; with Kilosort/Phy outputs (spike_times, spike_clusters, cluster_group): {ks}")
other_top = {a["path"].split("/")[0] for a in assets} - set("sub-" + sess.folder) - {"sourcedata"}
print("other top-level entries:", other_top)
ok = name_ok.all() and (has_vid == 1).all() and src_dirs == eph_names and oe == ks == len(eph_names)
print("criterion:", "reproduced" if ok else "partially_reproduced")

NWB files in their subject's folder and named as described: 115 of 115; types {'behavior+image': 71, 'behavior+ecephys+image': 44}
sessions with exactly one video in a folder named as the NWB file: 115 of 115
sourcedata session folders: 44; equal to the ecephys NWB names: True
folders with Open Ephys files: 44; with Kilosort/Phy outputs (spike_times, spike_clusters, cluster_group): 44
other top-level entries: set()
criterion: reproduced

print(f"containers {V.format.value_counts().to_dict()}; codecs {V.codec.value_counts().to_dict()}; "
      f"nominal rate {V.average_rate.value_counts().to_dict()}; extension mkv {V.path.str.endswith('.mkv').sum()} of {len(V)}")
print(f"frame sizes (width x height): {V.groupby(['width', 'height']).size().to_dict()}")
span = pd.Series({s["path"][:-4]: max(d["t_last"] for d in s["whiskers"].values()) for s in S})
V["vid_s"] = V.frames / 200
V["span_s"] = V.nwb_name.map(span)
dd = V.vid_s - V.span_s
print(f"video frames / 200 minus last whisker time: median {dd.median():.3f} s, range {dd.min():.3f} to {dd.max():.3f} s "
      f"(one frame = 0.005 s)")
n550 = ((V.width == 640) & (V.height == 550) | (V.width == 550) & (V.height == 640)).sum()
ok = (V.format.str.contains("matroska")).all() and (V.codec == "h264").all() and (V.average_rate == 30).all() and n550 > len(V) / 2
print("criterion:", "reproduced" if ok else "not_reproduced")

containers {'matroska,webm': 115}; codecs {'h264': 115}; nominal rate {30.0: 115}; extension mkv 115 of 115
frame sizes (width x height): {(640, 550): 111, (640, 600): 4}
video frames / 200 minus last whisker time: median 0.005 s, range 0.005 to 0.065 s (one frame = 0.005 s)
criterion: reproduced

subj = pd.DataFrame([dict(subject=s["subject"], **s["subject_meta"], start=pd.Timestamp(s["session_start_time"])) for s in S])
per = subj.groupby("subject").agg(sex=("sex", "first"), strain=("strain", "first"), dob=("date_of_birth", "first"),
                                  first=("start", "min"), n_sex=("sex", "nunique"))
per["age_days"] = [(f.tz_convert(None) - pd.Timestamp(d).tz_convert(None)).days for f, d in zip(per["first"], per.dob)]
print(per[["sex", "strain", "dob", "age_days"]].to_string())
norm = per.strain.str.upper().str.replace(r"[^A-Z0-9]", "", regex=True)
strain_ok = (norm == "C57BL6J").all()
age_ok = (per.age_days >= 90).all()
counts = per.sex.value_counts().to_dict()
k = min(counts.values()) if len(counts) == 2 else 0
p = stats.binomtest(k, len(per), 0.5).pvalue
ra, rb = p >= 0.05, k / len(per) >= 1 / 3
print(f"strain all C57BL/6J: {strain_ok}; sex {counts}; binomial p = {p:.3f} (reading a holds: {ra}); minority fraction "
      f"{k / len(per):.3f} (reading b holds: {rb}); first session age {per.age_days.min()} to {per.age_days.max()} days (all >= 90: {age_ok})")
print("criterion:", "reproduced" if strain_ok and age_ok and ra and rb else
      "blocked (sex readings disagree)" if strain_ok and age_ok and ra != rb else "not_reproduced")

        sex    strain                        dob  age_days
subject                                                   
219CR     F  C57BL6/J  2018-04-30T00:00:00-04:00       338
221CR     F  C57BL6/J  2018-05-20T00:00:00-04:00       341
228CR     F  C57BL6/J  2018-07-07T00:00:00-04:00       356
229CR     F  C57BL6/J  2018-09-14T00:00:00-04:00       384
231CR     M  C57BL6/J  2018-09-14T00:00:00-04:00       350
234CR     M  C57BL6/J  2018-10-01T00:00:00-04:00       330
245CR     F  C57BL6/J  2019-01-15T00:00:00-05:00       436
255CR     F  C57BL6/J  2019-03-11T00:00:00-04:00       377
267CR     F  C57BL6/J  2019-09-10T00:00:00-04:00       197
KF119     F  C57BL6/J  2017-01-01T00:00:00-05:00       373
KF132     F  C57BL6/J  2017-02-24T00:00:00-05:00       360
KF134     F  C57BL6/J  2017-02-25T00:00:00-05:00       345
KM101     M  C57BL6/J  2016-07-07T00:00:00-04:00       370
KM131     M  C57BL6/J  2017-02-06T00:00:00-05:00       336
KM91      M  C57BL6/J  2016-05-22T00:00:00-04:00       207
strain all C57BL/6J: True; sex {'F': 10, 'M': 5}; binomial p = 0.302 (reading a holds: True); minority fraction 0.333 (reading b holds: True); first session age 197 to 436 days (all >= 90: True)
criterion: reproduced

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

run time 6.8 min

Claim	Tested or blocked
dataset-size, neural-sessions, neural-hours	tested
task-rules, nothing-detection-only, one-task-per-mouse, servo-positions	tested
response-window-timing	blocked: no motor times in the dataset (the distribution of the window opening after the trial start is shown)
choice-time-first-lick, outcome-consistency, reward-and-timeout, forced-alternation-start, ignore-trial-definition, trials-table-structure	tested
video-timebase, sync-gap-and-interpolation, bend-null-interpolated, neural-timestamps	tested
neural-video-alignment	blocked: the analog input that carries the pulse is not deposited
behavior-clock-drift	blocked: the behavioral (Arduino) times are not in the dataset
unit-quality-isi, unit-depth-layer, electrodes-order, recording-columns, units-obs-intervals	tested
pose-structure, c0-rare-contacts, session-min-trials, session-labeling-quality, contacts-structure, nwb-structure, file-layout, video-format, subjects-sex-age-strain	tested
tracking-mean-error, tracking-error-dataset1, tracking-error-challenge	blocked: the curated training sets and predictions are not in the dataset
prior-task-strategies, prior-orthogonal-subspaces	blocked: findings of earlier papers, no method given here

Claim	Outcome
dataset-size	not reproduced: mice, sessions and frames match exactly; trials 32,350 or 22,905 against 18,514
neural-sessions	reproduced
neural-hours	reproduced (33.77 h in the main Open Ephys recordings)
task-rules	reproduced
nothing-detection-only	reproduced
one-task-per-mouse	reproduced
servo-positions	reproduced
response-window-timing	blocked: data missing (no motor times)
choice-time-first-lick	not reproduced: 77.6% of choice times are the first lick in the window; spoiled trials have NaN
outcome-consistency	reproduced
reward-and-timeout	blocked: methods underspecified (timeout readings 12.21 and 8.34 s disagree; right rewards untestable)
forced-alternation-start	reproduced
ignore-trial-definition	not reproduced: 160 spoiled trials not ignored
trials-table-structure	not reproduced: 5 of 32,235 pairs, all at missing trials
video-timebase	reproduced
sync-gap-and-interpolation	reproduced
bend-null-interpolated	reproduced
neural-timestamps	reproduced
neural-video-alignment	blocked: data missing (analog input not deposited)
behavior-clock-drift	blocked: data missing (no Arduino times)
unit-quality-isi	not reproduced: ISI part holds; firing-rate part fails under both bin readings
unit-depth-layer	reproduced
electrodes-order	reproduced
recording-columns	reproduced
units-obs-intervals	reproduced
pose-structure	reproduced
c0-rare-contacts	reproduced (10 of 15 mice)
session-min-trials	blocked: methods underspecified (one session has 84 or 101 trials by reading)
session-labeling-quality	not reproduced: up to 7.7% of a whisker's contacts hold a missing frame (2.007% pooled over whiskers)
contacts-structure	reproduced
nwb-structure	reproduced, under the naming of Table 1
file-layout	reproduced
video-format	reproduced
subjects-sex-age-strain	reproduced
tracking-mean-error, tracking-error-dataset1, tracking-error-challenge	blocked: data missing
prior-task-strategies, prior-orthogonal-subspaces	blocked: methods underspecified

Reproducing Rodgers (2022) from dandiset 000231¶

Recordings deposited more than once¶

dataset-size¶

neural-sessions¶

neural-hours¶

task-rules¶

nothing-detection-only¶

one-task-per-mouse¶

servo-positions¶

response-window-timing (blocked)¶

choice-time-first-lick¶

outcome-consistency¶

reward-and-timeout¶

forced-alternation-start¶

ignore-trial-definition¶

trials-table-structure¶

video-timebase¶

sync-gap-and-interpolation¶

bend-null-interpolated¶

neural-timestamps¶

neural-video-alignment (blocked)¶

behavior-clock-drift (blocked)¶

unit-quality-isi¶

unit-depth-layer¶

electrodes-order¶

recording-columns¶

units-obs-intervals¶

pose-structure¶

c0-rare-contacts¶

session-min-trials¶

session-labeling-quality¶

contacts-structure¶

nwb-structure¶

file-layout¶

video-format¶

subjects-sex-age-strain¶

Blocked: tracking-mean-error, tracking-error-dataset1, tracking-error-challenge¶

Blocked: prior-task-strategies, prior-orthogonal-subspaces¶

Summary¶

Reproducing Rodgers (2022) from dandiset 000231¶

Recordings deposited more than once¶

Do the streams share the clock of the trials table?¶

dataset-size¶

neural-sessions¶

neural-hours¶

task-rules¶

nothing-detection-only¶

one-task-per-mouse¶

servo-positions¶

response-window-timing (blocked)¶

choice-time-first-lick¶

outcome-consistency¶

reward-and-timeout¶

forced-alternation-start¶

ignore-trial-definition¶

trials-table-structure¶

video-timebase¶

sync-gap-and-interpolation¶

bend-null-interpolated¶

neural-timestamps¶

neural-video-alignment (blocked)¶

behavior-clock-drift (blocked)¶

unit-quality-isi¶

unit-depth-layer¶

electrodes-order¶

recording-columns¶

units-obs-intervals¶

pose-structure¶

c0-rare-contacts¶

session-min-trials¶

session-labeling-quality¶

contacts-structure¶

nwb-structure¶

file-layout¶

video-format¶

subjects-sex-age-strain¶

Blocked: tracking-mean-error, tracking-error-dataset1, tracking-error-challenge¶

Blocked: prior-task-strategies, prior-orthogonal-subspaces¶

Summary¶