import time
import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats
from statsmodels.formula.api import ols

import stretch

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

started = time.time()
S = stretch.load_all(stretch.asset_urls())
print(f"loaded {len(S)} files of dandiset {stretch.DANDISET_ID} version {stretch.DANDISET_VERSION} "
      f"in {time.time() - started:.0f} s")

loaded 447 files of dandiset 001636 version 0.260518.1819 in 170 s

sessions = pd.DataFrame([{
    "si": si, "session": s["session_id"], "code": s["session_id"].split("++")[1],
    "subject": s["subject"], "post": s["post"], "date": s["date"], "n_trials": len(s["trials"]),
    "n_flex": int((s["pert_direction"] == "flexion").sum()),
    "n_ext": int((s["pert_direction"] == "extension").sum()), "n_units": len(s["units"]),
    "analog_dt_ms": s.get("analog_dt_ms", np.nan), "has_analog": "windows" in s or "analog_dt_ms" in s,
    "angle_all_nan": s.get("angle_all_nan", True), "torque_all_nan": s.get("torque_all_nan", True),
    "emg_channels": len(s.get("emg_channels", [])), "emg_same_clock": s.get("emg_same_clock", np.nan),
} for si, s in enumerate(S)])
units = pd.DataFrame([{
    "si": si, "ui": ui, "session": s["session_id"], "subject": s["subject"], "post": s["post"],
    "is_post_mptp": u["is_post_mptp"], "projection": u["type"], "ptype": stretch.short_type(u["type"]),
    "unit_name": u["unit_name"], "also_in": u["also_in"], "date": s["date"],
    "n_flex": int((s["pert_direction"] == "flexion").sum()),
    "n_ext": int((s["pert_direction"] == "extension").sum()),
    "antidromic_latency_ms": u["antidromic_latency_ms"], "n_spikes": len(u["spikes"]),
    "isi_below_2ms": u["isi_below_2ms"], "outside_trials": u["outside_trials"], "sorted": u["sorted"],
} for si, s in enumerate(S) for ui, u in enumerate(s["units"])])
units["state"] = np.where(units.post, "post", "pre")
units["studied"] = (units.n_flex >= stretch.MIN_REPS) & (units.n_ext >= stretch.MIN_REPS)

print("Files and units by subject and MPTP state (from the session id):")
print(pd.concat([sessions.groupby(["subject", "post"]).size().rename("files"),
                 units.groupby(["subject", "post"]).size().rename("units")], axis=1).to_string())
print("\nUnits whose is_post_mptp differs from the state in the session id:",
      int((units.post != units.is_post_mptp).sum()))
print("\nRecording dates by subject:")
print(sessions.groupby("subject").date.agg(["min", "max"]).to_string())
print("\nProjection types as stored, by subject and state:")
print(pd.crosstab(units.projection, [units.subject, units.state]).to_string())
print("\nTorque perturbations per file (flexion / extension), by subject and state:")
print(sessions.groupby(["subject", "post"])[["n_flex", "n_ext"]].describe()
      .loc[:, (slice(None), ["min", "50%", "max"])].to_string())
print("\nUnits studied (at least 5 perturbations of each direction), by subject, state and type:")
print(pd.crosstab([units.subject, units.state], [units.studied, units.ptype]).to_string())

Files and units by subject and MPTP state (from the session id):
               files  units
subject post               
Leu     False     98    138
        True      51     72
Venus   False    157    194
        True     141    164

Units whose is_post_mptp differs from the state in the session id: 0

Recording dates by subject:
               min        max
subject                      
Leu     1992-08-06 1993-12-01
Venus   1999-04-06 2000-04-03

Projection types as stored, by subject and state:
subject                                     Leu     Venus    
state                                      post pre  post pre
projection                                                   
corticostriatal_neuron                       10  25    49  51
no_response                                  46  84    52  63
not_tested                                    0   0     4   3
pyramidal_tract_and_corticostriatal_neuron    0   0     1   2
pyramidal_tract_and_subthalamic_neuron        0   1     0   0
pyramidal_tract_and_thalamic_neuron           0   0     1   0
pyramidal_tract_neuron                       10  19    55  75
subthalamic_and_corticostriatal_neuron        0   1     0   0
subthalamic_and_thalamic_neuron               0   1     0   0
subthalamic_projection_neuron                 6   7     0   0
thalamic_projection_neuron                    0   0     2   0

Torque perturbations per file (flexion / extension), by subject and state:
              n_flex             n_ext            
                 min   50%   max   min   50%   max
subject post                                      
Leu     False    0.0   0.0   0.0   0.0   0.0   0.0
        True     0.0   0.0   0.0   0.0   0.0   0.0
Venus   False    0.0   6.0  12.0   1.0   6.0  12.0
        True     1.0  13.0  28.0   3.0  13.0  28.0

Units studied (at least 5 perturbations of each direction), by subject, state and type:
studied       False         True              
ptype           CSN  NA PTN   CSN  NA PTN dual
subject state                                 
Leu     post     10  52  10     0   0   0    0
        pre      25  94  19     0   0   0    0
Venus   post      1   3   1    48  56  54    1
        pre      20  36  26    31  30  49    2

hold = []
for s in S:
    if "hold_excursion_deg" not in s["trials"]:
        continue
    tr = s["trials"]
    hold.append(pd.DataFrame({"subject": s["subject"], "post": s["post"], "perturbed": tr.torque != "none",
                              "excursion": tr.hold_excursion_deg, "movement": tr.movement,
                              "active_change": tr.active_angle_change, "end_position": tr.end_position}))
hold = pd.concat(hold, ignore_index=True)
print("Largest angle excursion from its median, 1.0 to 2.5 s after the center target appears (deg):")
print(hold.groupby(["subject", "post", "perturbed"]).excursion.describe(percentiles=[.1, .5, .9]).round(2).to_string())
print("\nFraction of trials with an excursion above 2 deg:")
print(hold.assign(big=hold.excursion > 2).groupby(["subject", "post", "perturbed"]).big.mean().round(3).to_string())
leu = sessions[sessions.subject == "Leu"]
print(f"\nLeu files with no analog data: {int((~leu.has_analog).sum())} of {len(leu)} "
      f"({int((~leu[leu.post].has_analog).sum())} of {int(leu.post.sum())} post-MPTP)")
print("Leu files with analog data whose torque command is all NaN, by state:")
print(leu[leu.has_analog].groupby("post").torque_all_nan.agg(["sum", "count"]).to_string())
print("Analog sample interval (ms) by subject:")
print(sessions.groupby("subject").analog_dt_ms.describe()[["count", "min", "50%", "max"]].round(4).to_string())

Largest angle excursion from its median, 1.0 to 2.5 s after the center target appears (deg):
                          count  mean   std   min   10%   50%   90%    max
subject post  perturbed                                                   
Leu     False False      2997.0  1.06  3.35  0.00  0.02  0.08  2.66  45.24
        True  False      2091.0  0.17  1.37  0.00  0.02  0.05  0.16  39.62
Venus   False False       958.0  0.20  0.64  0.01  0.01  0.03  0.48  17.27
              True       1862.0  2.79  1.28  0.35  1.27  2.72  4.40  19.56
        True  False      1829.0  0.89  0.99  0.01  0.07  0.57  2.18  12.34
              True       3649.0  3.31  1.44  0.11  1.70  3.19  4.77  31.75

Fraction of trials with an excursion above 2 deg:
subject  post   perturbed
Leu      False  False        0.148
         True   False        0.011
Venus    False  False        0.007
                True         0.657
         True   False        0.108
                True         0.863

Leu files with no analog data: 15 of 149 (15 of 51 post-MPTP)
Leu files with analog data whose torque command is all NaN, by state:
       sum  count
post             
False   98     98
True     1     36
Analog sample interval (ms) by subject:
         count      min      50%      max
subject                                  
Leu      134.0  10.0156  10.0176  10.0192
Venus    298.0   1.0001   1.0001   1.0001

# Every excursion of the torque command (more than 20 units from rest, lasting at least 10 ms),
# in every file with a torque command; a hold-period impulse starts in the hold window.
exc_rows = []
for s in S:
    if "command_excursions" not in s:
        continue
    e = s["command_excursions"]
    e = e[e.duration_ms >= 10]
    h = e[e.in_hold]
    exc_rows.append({"subject": s["subject"], "state": "post" if s["post"] else "pre", "code": s["session_id"].split("++")[1],
                     "trials": len(s["trials"]), "excursions": len(e),
                     "labelled_torque_trials": int((s["trials"].torque != "none").sum()),
                     "trials_hold_excursion": h.trial.nunique(),
                     "trials_hold_40_70ms": h[(h.duration_ms >= 40) & (h.duration_ms <= 70)].trial.nunique(),
                     "hold_median_ms": h.duration_ms.median(), "hold_median_abs_peak": h.peak.abs().median(),
                     "all_median_from_center_s": e.from_center_s.median(),
                     "all_median_from_lateral_s": e.from_lateral_s.median()})
exc = pd.DataFrame(exc_rows)
print("Files with a torque command, by subject and state (trials with an excursion starting in the hold window,"
      " 1.0 to 2.5 s after the center target and before the lateral target):")
print(exc.groupby(["subject", "state"]).agg(files=("code", "size"), trials=("trials", "sum"),
      labelled_torque_trials=("labelled_torque_trials", "sum"), trials_hold_excursion=("trials_hold_excursion", "sum"),
      trials_hold_40_70ms=("trials_hold_40_70ms", "sum")).to_string())
print("\nEvery monkey L file with a torque command:")
print(exc[exc.subject == "Leu"].drop(columns="subject").round(2).to_string(index=False))

Files with a torque command, by subject and state (trials with an excursion starting in the hold window, 1.0 to 2.5 s after the center target and before the lateral target):
               files  trials  labelled_torque_trials  trials_hold_excursion  trials_hold_40_70ms
subject state                                                                                   
Leu     post      35    2048                       0                      1                    0
Venus   post     141    5573                    3649                   3563                 2596
        pre      157    2820                    1862                   1811                 1811

Every monkey L file with a torque command:
state    code  trials  excursions  labelled_torque_trials  trials_hold_excursion  trials_hold_40_70ms  hold_median_ms  hold_median_abs_peak  all_median_from_center_s  all_median_from_lateral_s
 post  l14001      61         349                       0                      0                    0             NaN                   NaN                      3.74                       0.54
 post  l14002      60         355                       0                      0                    0             NaN                   NaN                      3.69                       0.54
 post  l14101      60         353                       0                      0                    0             NaN                   NaN                      3.52                       0.42
 post  l14102      60         352                       0                      0                    0             NaN                   NaN                      3.67                       0.49
 post  l14103      60         359                       0                      0                    0             NaN                   NaN                      3.70                       0.54
 post  l14201      60         359                       0                      0                    0             NaN                   NaN                      3.68                       0.49
 post  l14501      85         381                       0                      0                    0             NaN                   NaN                      3.73                       0.50
 post  l14601      60         352                       0                      0                    0             NaN                   NaN                      3.66                       0.41
 post  l14602      61         354                       0                      0                    0             NaN                   NaN                      3.54                       0.43
 post  l14701      62         304                       0                      0                    0             NaN                   NaN                      3.93                       0.84
 post  l14801      88         352                       0                      0                    0             NaN                   NaN                      3.33                       0.35
 post  l14802      18          91                       0                      0                    0             NaN                   NaN                      3.85                       0.67
 post  l14803      60         215                       0                      0                    0             NaN                   NaN                      3.02                       0.31
 post  l14804      60         239                       0                      0                    0             NaN                   NaN                      3.31                       0.34
 post  l15001      62         229                       0                      0                    0             NaN                   NaN                      3.86                       0.73
 post  l15601      60         410                       0                      0                    0             NaN                   NaN                      3.45                       0.34
 post  l15602      37         260                       0                      0                    0             NaN                   NaN                      3.52                       0.43
 post  l15701      60         413                       0                      0                    0             NaN                   NaN                      3.56                       0.39
 post  l15702      62         412                       0                      0                    0             NaN                   NaN                      3.41                       0.39
 post  l15703      60         400                       0                      0                    0             NaN                   NaN                      3.48                       0.33
 post  l15801      60         395                       0                      0                    0             NaN                   NaN                      3.41                       0.39
 post  l15901      61         399                       0                      0                    0             NaN                   NaN                      3.57                       0.38
 post  l16001      60         385                       0                      0                    0             NaN                   NaN                      3.48                       0.39
 post  l16002      60         393                       0                      0                    0             NaN                   NaN                      3.43                       0.39
 post  l16101      60         385                       0                      0                    0             NaN                   NaN                      3.27                       0.32
 post  l16102      43         284                       0                      0                    0             NaN                   NaN                      3.38                       0.36
 post  l16103      47         297                       0                      1                    0           70.12                 41.12                      3.20                       0.34
 post  l16201      60         388                       0                      0                    0             NaN                   NaN                      3.41                       0.32
 post  l16301      60         381                       0                      0                    0             NaN                   NaN                      3.27                       0.37
 post l16301b      60         391                       0                      0                    0             NaN                   NaN                      3.33                       0.36
 post  l16302      60         382                       0                      0                    0             NaN                   NaN                      3.12                       0.38
 post  l16401      60         402                       0                      0                    0             NaN                   NaN                      3.54                       0.36
 post  l16402      60         380                       0                      0                    0             NaN                   NaN                      3.26                       0.35
 post  l16403      61         381                       0                      0                    0             NaN                   NaN                      3.30                       0.34
 post  l16501      40         266                       0                      0                    0             NaN                   NaN                      3.38                       0.34

venus = [s for s in S if s["subject"] == "Venus" and "windows" in s]
AX = np.arange(stretch.ANALOG_WINDOW_MS[0], stretch.ANALOG_WINDOW_MS[1] + 1)


def command_onset(w):
    # first ms where the torque command leaves its pre-onset level by more than 6 SD (and 10 units)
    base = np.nanmedian(w[:, :400], axis=1, keepdims=True)
    sd = np.nanstd(w[:, :400], axis=1, keepdims=True)
    dev = np.abs(w - base) > np.maximum(10, 6 * sd)
    return np.where(dev.any(1), AX[np.argmax(dev, axis=1)], np.nan)


realigned = np.concatenate([command_onset(s["windows"]["torque_cmd"]) for s in venus])
stored = np.concatenate([command_onset(s["windows"]["torque_cmd_stored"]) for s in venus])
lag = np.concatenate([s["analog_lag_ms"][s["pert_trial"]] for s in venus])
ok = np.isfinite(stored) & np.isfinite(realigned)
print(f"Monkey V perturbations: {len(realigned)}; command step found in {int(np.isfinite(realigned).sum())}")
print("Torque command step relative to the labelled onset, realigned clock (ms), percentiles 1/50/99:",
      np.nanpercentile(realigned, [1, 50, 99]))
print("Same on the stored analog clock (ms), percentiles 1/50/99:", np.nanpercentile(stored, [1, 50, 99]))
print(f"Stored analog clock behind the trials table at each perturbation: {lag.min():.0f} to {lag.max():.0f} ms; "
      f"correlation of the stored-clock step time with this lag {np.corrcoef(stored[ok], lag[ok])[0, 1]:.3f}")
per_trial = np.concatenate([np.diff(s["analog_lag_ms"]) for s in S if "analog_lag_ms" in s])
print(f"Lag added per trial: median {np.median(per_trial):.1f} ms (range {per_trial.min():.1f} to {per_trial.max():.1f})")
print("EMG series share the angle's timestamps in every monkey V file:",
      bool(sessions[sessions.subject == "Venus"].emg_same_clock.astype(bool).all()))
print("\nStored angle change over active movements (lateral target to reward), median deg:")
print(hold.groupby(["subject", "movement"]).active_change.median().round(1).to_string())
print("derived_end_position, median deg:")
print(hold.groupby(["subject", "movement"]).end_position.median().round(1).to_string())
chg = {d: np.concatenate([s["windows"]["angle"][s["pert_direction"] == d][:, AX == 100][:, 0]
                          - s["windows"]["angle"][s["pert_direction"] == d][:, (AX >= -50) & (AX < 0)].mean(1)
                          for s in venus]) for d in ["flexion", "extension"]}
for d, v in chg.items():
    print(f"Angle change 100 ms after a {d} torque: median {np.nanmedian(v):.2f} deg, positive in "
          f"{np.nanmean(v > 0):.1%} of {np.isfinite(v).sum()} perturbations")

Monkey V perturbations: 5511; command step found in 5474
Torque command step relative to the labelled onset, realigned clock (ms), percentiles 1/50/99: [4. 5. 7.]
Same on the stored analog clock (ms), percentiles 1/50/99: [-108.  -42.    6.]
Stored analog clock behind the trials table at each perturbation: 0 to 365 ms; correlation of the stored-clock step time with this lag -0.998
Lag added per trial: median 5.0 ms (range 4.0 to 14.0)
EMG series share the angle's timestamps in every monkey V file: True

Stored angle change over active movements (lateral target to reward), median deg:
subject  movement 
Leu      extension   -16.7
         flexion      17.7
Venus    extension   -20.9
         flexion      18.3
derived_end_position, median deg:
subject  movement 
Leu      extension   -18.6
         flexion      20.1
Venus    extension   -20.2
         flexion      19.1
Angle change 100 ms after a flexion torque: median 3.40 deg, positive in 99.8% of 2745 perturbations
Angle change 100 ms after a extension torque: median -2.48 deg, positive in 1.1% of 2766 perturbations

# Angle scale: active movements measured here (same filter as the torque-evoked movements, around
# derived_movement_onset_time) against the files' derived columns and target amplitudes.
act = pd.concat([s["trials"].assign(subject=s["subject"]) for s in S if "active_peak_velocity" in s["trials"]])
act = act[act.active_peak_velocity.notna() & act.derived_peak_velocity.notna()]
act["ratio_velocity"] = act.active_peak_velocity / act.derived_peak_velocity
act["ratio_amplitude"] = act.active_amplitude / act.derived_amplitude.abs()
print(act.groupby("subject").agg(movements=("ratio_velocity", "size"),
      target_amplitude_median=("target_amplitude", "median"),
      derived_amplitude_median=("derived_amplitude", lambda x: x.abs().median()),
      amplitude_here_median=("active_amplitude", "median"),
      derived_peak_velocity_median=("derived_peak_velocity", "median"),
      peak_velocity_here_median=("active_peak_velocity", "median"),
      ratio_velocity_median=("ratio_velocity", "median"),
      ratio_amplitude_median=("ratio_amplitude", "median")).round(2).to_string())
for subj, g in act.groupby("subject"):
    print(f"{subj}: Spearman rho of peak velocity here against derived_peak_velocity "
          f"{stats.spearmanr(g.active_peak_velocity, g.derived_peak_velocity)[0]:.3f}; "
          f"velocity ratio within 10% of 1 in {((g.ratio_velocity - 1).abs() <= 0.1).mean():.1%}")

         movements  target_amplitude_median  derived_amplitude_median  amplitude_here_median  derived_peak_velocity_median  peak_velocity_here_median  ratio_velocity_median  ratio_amplitude_median
subject                                                                                                                                                                                             
Leu           4793                     20.0                     19.18                  19.86                        100.52                     101.43                   1.00                    1.01
Venus         7491                     20.0                     19.46                  20.40                         99.86                      99.33                   0.99                    1.02
Leu: Spearman rho of peak velocity here against derived_peak_velocity 0.989; velocity ratio within 10% of 1 in 93.5%
Venus: Spearman rho of peak velocity here against derived_peak_velocity 0.993; velocity ratio within 10% of 1 in 99.5%

# Spikes against the trials table: does the torque direction become decodable from single-cell
# spike counts only after the labelled onset? 20 ms windows, Mann-Whitney flexion vs extension.
edges = np.arange(-200, 201, 20)
studied_v = units[units.studied & (units.subject == "Venus")]
sig = np.zeros(len(edges) - 1)
for r in studied_v.itertuples():
    s, u = S[r.si], S[r.si]["units"][r.ui]
    counts = np.array([np.histogram(x, edges)[0] for x in u["pert"]])
    f, e = counts[s["pert_direction"] == "flexion"], counts[s["pert_direction"] == "extension"]
    for j in range(len(edges) - 1):
        if np.ptp(np.r_[f[:, j], e[:, j]]) > 0:
            sig[j] += stats.mannwhitneyu(f[:, j], e[:, j]).pvalue < 0.01
frac = sig / len(studied_v)
print(f"Fraction of the {len(studied_v)} monkey V cells studied whose spike counts differ between "
      "flexion and extension trials (Mann-Whitney p < 0.01), by window (ms from labelled onset):")
print(pd.Series(frac.round(3), index=[f"{a} to {b}" for a, b in zip(edges[:-1], edges[1:])]).to_string())

Fraction of the 271 monkey V cells studied whose spike counts differ between flexion and extension trials (Mann-Whitney p < 0.01), by window (ms from labelled onset):
-200 to -180    0.004
-180 to -160    0.000
-160 to -140    0.000
-140 to -120    0.000
-120 to -100    0.011
-100 to -80     0.000
-80 to -60      0.004
-60 to -40      0.000
-40 to -20      0.000
-20 to 0        0.007
0 to 20         0.011
20 to 40        0.070
40 to 60        0.089
60 to 80        0.048
80 to 100       0.063
100 to 120      0.048
120 to 140      0.030
140 to 160      0.033
160 to 180      0.052
180 to 200      0.018

# Files that repeat another file's recording: same subject and date, and the trial durations of one
# are a contiguous block of the other's.
dur = {si: np.round((s["trials"].stop - s["trials"].start).values, 3) for si, s in enumerate(S)}
repeats = []
for (subj, date), g in sessions.groupby(["subject", "date"]):
    for a in g.si:
        for b in g.si:
            if a != b and len(dur[b]) <= len(dur[a]) and (len(dur[b]) < len(dur[a]) or a < b):
                n = len(dur[b])
                if any(np.array_equal(dur[a][i:i + n], dur[b]) for i in range(len(dur[a]) - n + 1)):
                    repeats.append({"file": sessions.code[a], "repeated_in": sessions.code[b],
                                    "trials": len(dur[a]), "repeated_trials": n,
                                    "units_file": [u["unit_name"] for u in S[a]["units"]],
                                    "units_repeat": [u["unit_name"] for u in S[b]["units"]],
                                    "types_file": [stretch.short_type(u["type"]) for u in S[a]["units"]],
                                    "types_repeat": [stretch.short_type(u["type"]) for u in S[b]["units"]]})
repeats = pd.DataFrame(repeats)
print("Files whose trials repeat a block of another file's trials:")
print(repeats.to_string(index=False))


def trial_spikes(s, u):
    # spike times relative to the start of their trial, one array per trial
    tr = s["trials"]
    k = np.searchsorted(tr.start.values, u["spikes"], side="right") - 1
    return [np.round(u["spikes"][k == i] - tr.start.values[i], 4) for i in range(len(tr))]


same_unit = []
for r in repeats.itertuples():
    a = S[int(sessions[sessions.code == r.file].si.iloc[0])]
    b = S[int(sessions[sessions.code == r.repeated_in].si.iloc[0])]
    da, db = dur[int(sessions[sessions.code == r.file].si.iloc[0])], dur[int(sessions[sessions.code == r.repeated_in].si.iloc[0])]
    off = [i for i in range(len(da) - len(db) + 1) if np.array_equal(da[i:i + len(db)], db)][0]
    for ub in b["units"]:
        tb = trial_spikes(b, ub)
        for ua in a["units"]:
            ta = trial_spikes(a, ua)[off:off + len(db)]
            n_same = sum(np.array_equal(x, y) for x, y in zip(ta, tb))
            same_unit.append({"file": r.file, "unit": ua["unit_name"], "type": ua["type"],
                              "repeat": r.repeated_in, "repeat_unit": ub["unit_name"], "repeat_type": ub["type"],
                              "identical_trials": f"{n_same}/{len(db)}",
                              "spikes_in_block": (sum(map(len, ta)), sum(map(len, tb)))})
same_unit = pd.DataFrame(same_unit)
print("\nUnits of each repeated recording compared trial by trial over the repeated block:")
print(same_unit.to_string(index=False))
dup_units = same_unit[same_unit.identical_trials.str.split("/", regex=False).str[0].astype(int) > 0]
# the copy in the repeated file is the one dropped in the variant that counts each unit once
dup_keys = {(r.repeat, r.repeat_unit) for r in dup_units.itertuples()}
units["duplicate_copy"] = [(sess.split("++")[1], un) in dup_keys for sess, un in zip(units.session, units.unit_name)]
print("Unit copies in repeated files (dropped in the variant that counts each unit once):", sorted(dup_keys))

codes = set(sessions.code)
xref = units[units.also_in != ""].copy()
xref["ref_code"] = xref.also_in.str.split("++", regex=False).str[1]
xref["ref_prefix"] = xref.also_in.str.split("++", regex=False).str[0]
xref["ref_in_dandiset"] = xref.ref_code.isin(codes)
print(f"\nUnits with unit_also_in_session_id: {len(xref)}; prefixes used {xref.ref_prefix.value_counts().to_dict()} "
      f"(session ids use {sorted(sessions.subject.unique())}); referenced file in the dandiset: "
      f"{int(xref.ref_in_dandiset.sum())}")
rows = []
for r in xref[xref.ref_in_dandiset].itertuples():
    other = S[int(sessions[sessions.code == r.ref_code].si.iloc[0])]
    rows.append({"unit": f"{r.session.split('++')[1]}/{r.unit_name}", "type": r.ptype, "spikes": r.n_spikes,
                 "referenced": r.ref_code,
                 "its units": [(u["unit_name"], stretch.short_type(u["type"]), len(u["spikes"])) for u in other["units"]]})
print(pd.DataFrame(rows).to_string(index=False))
print("\nSpike times sorted in every unit:", bool(units["sorted"].all()),
      "| spikes outside the trials:", int(units.outside_trials.sum()))
print("Within-trial ISIs below 2 ms, fraction per unit: median "
      f"{units.isi_below_2ms.median():.4f}, units above 1%: {int((units.isi_below_2ms > 0.01).sum())} of {len(units)} "
      f"({int((units[units.studied & (units.subject == 'Venus')].isi_below_2ms > 0.01).sum())} of the cells studied)")

Files whose trials repeat a block of another file's trials:
 file repeated_in  trials  repeated_trials units_file units_repeat types_file types_repeat
v4527      v4527b      47               27     [1, 2]          [2]  [PTN, NA]        [PTN]
v4810      v4810b      80               40        [1]          [2]      [PTN]         [NA]
v5604      v5604b      80               40        [1]          [2]      [PTN]        [CSN]
v5811       v5812      80               40     [1, 2]          [2] [PTN, PTN]         [NA]

Units of each repeated recording compared trial by trial over the repeated block:
 file  unit                   type repeat  repeat_unit            repeat_type identical_trials spikes_in_block
v4527     1 pyramidal_tract_neuron v4527b            2 pyramidal_tract_neuron             0/27    (3389, 4403)
v4527     2             not_tested v4527b            2 pyramidal_tract_neuron            27/27    (4403, 4403)
v4810     1 pyramidal_tract_neuron v4810b            2            no_response             0/40     (2155, 595)
v5604     1 pyramidal_tract_neuron v5604b            2 corticostriatal_neuron             0/40      (3327, 29)
v5811     1 pyramidal_tract_neuron  v5812            2            no_response             0/40    (3693, 2265)
v5811     2 pyramidal_tract_neuron  v5812            2            no_response            40/40    (2265, 2265)
Unit copies in repeated files (dropped in the variant that counts each unit once): [('v4527b', 2), ('v5812', 2)]

Units with unit_also_in_session_id: 43; prefixes used {'V': 42, 'L': 1} (session ids use ['Leu', 'Venus']); referenced file in the dandiset: 8
    unit type  spikes referenced        its units
l16301/1  CSN    3444    l16301b   [(2, NA, 861)]
 v2701/1  CSN     536      v2702   [(1, NA, 181)]
 v3204/1  PTN    1940      v3205  [(2, CSN, 661)]
 v3607/1  PTN    3375     v3607b  [(2, CSN, 141)]
 v4527/1  PTN    5578     v4527b [(2, PTN, 4403)]
 v4810/1  PTN    4229     v4810b   [(2, NA, 595)]
 v5604/1  PTN    4008     v5604b   [(2, CSN, 29)]
 v5811/1  PTN    6933      v5812  [(2, NA, 2265)]

Spike times sorted in every unit: True | spikes outside the trials: 0
Within-trial ISIs below 2 ms, fraction per unit: median 0.0006, units above 1%: 69 of 568 (37 of the cells studied)

emg_rows = []
for s in venus:
    for ch in stretch.EMG_CHANNELS[:3]:
        if ch not in s["windows"]:
            state = "absent"
        else:
            w = s["windows"][ch]
            if np.isnan(w).any():
                state = "NaN"
            elif np.nanstd(w[:, AX < 0]) < 1:
                state = "flat"
            else:
                state = "used"
        emg_rows.append({"code": s["session_id"].split("++")[1], "post": s["post"], "channel": ch, "state": state})
emg_avail = pd.DataFrame(emg_rows)
emg_avail["repeat"] = emg_avail.code.isin(repeats.repeated_in)
print("Monkey V triceps channels per file (repeated files excluded), by state:")
print(pd.crosstab([emg_avail[~emg_avail.repeat].channel], [emg_avail[~emg_avail.repeat].post, emg_avail[~emg_avail.repeat].state]).to_string())
vals = np.concatenate([s["windows"][ch].ravel() for s in venus for ch in stretch.EMG_CHANNELS[:3] if ch in s["windows"]])
print(f"\nTriceps EMG values (A/D units, as stored in a.u.): {np.nanmin(vals):.0f} to {np.nanmax(vals):.0f}, "
      f"median {np.nanmedian(vals):.0f}")
ang = np.concatenate([s["windows"]["angle"].ravel() for s in venus])
print(f"Monkey V angle around perturbations: {np.nanmin(ang):.1f} to {np.nanmax(ang):.1f} deg")
tq = np.concatenate([s["windows"]["torque_cmd"].ravel() for s in venus])
print(f"Monkey V torque command (arbitrary units): {np.nanmin(tq):.0f} to {np.nanmax(tq):.0f}")
iv = np.concatenate([(s["trials"].onset - s["trials"].center_on).dropna().values for s in venus])
# size of the torque command step (sign-corrected, 10 to 40 ms after onset, less the 100 ms before)
command = []
for s in venus:
    w = s["windows"]["torque_cmd"]
    sign = np.where(s["pert_direction"] == "flexion", 1.0, -1.0)[:, None]
    rel = (w - np.nanmedian(w[:, (AX >= -100) & (AX < 0)], axis=1, keepdims=True)) * sign
    m = np.nanmean(rel, axis=0)
    command.append({"session": s["session_id"], "code": s["session_id"].split("++")[1], "post": s["post"],
                    "date": s["date"], "step": float(np.nanmedian(rel[:, (AX >= 10) & (AX < 40)])),
                    "width_ms": int((m > m.max() / 2).sum()), "after_pulse": float(m[(AX >= 70) & (AX < 150)].mean())})
command = pd.DataFrame(command)
command["level"] = np.round(command.step / command.step[command.step > 10].min()).astype(int)
print("\nTorque command step (units above the pre-onset level) by state; level = step in multiples of the smallest:")
print(pd.crosstab(command.level, command.post).to_string())
print(command.groupby(["post", "level"])[["step", "width_ms", "after_pulse"]].median().round(1).to_string())
print("Dates of the level-2 sessions by state:", command[command.level == 2].groupby("post").date.agg(["min", "max", "count"]).to_string())
print("Torque onset after the center target appears (s), percentiles 0/5/50/95/100:",
      np.percentile(iv, [0, 5, 50, 95, 100]).round(3))

Monkey V triceps channels per file (repeated files excluded), by state:
post                False                  True                  
state                 NaN absent flat used   NaN absent flat used
channel                                                          
EMGTricepsLateralis     6    115    0   36     0     97    4   36
EMGTricepsLongus1       2      0    9  146     0      4    3  130
EMGTricepsLongus2       0      0   27  130    13      4    4  116

Triceps EMG values (A/D units, as stored in a.u.): 1705 to 4294, median 2072
Monkey V angle around perturbations: -28.9 to 29.8 deg
Monkey V torque command (arbitrary units): 1846 to 2233

Torque command step (units above the pre-onset level) by state; level = step in multiples of the smallest:
post   False  True 
level              
0          3      0
1        108      0
2         43     12
4          3    129
              step  width_ms  after_pulse
post  level                              
False 0       -0.1      15.0         -0.4
      1       45.7      50.0         -0.3
      2       91.3      50.0         -0.0
      4      182.7      50.0         -0.2
True  2       91.2      50.0        -30.3
      4      182.5      50.0        -37.1
Dates of the level-2 sessions by state:              min        max  count
post                              
False 1999-04-06 1999-06-23     43
True  2000-01-13 2000-02-16     12
Torque onset after the center target appears (s), percentiles 0/5/50/95/100: [ 1.257  1.696  1.696  1.732 11.841]

cells = units[units.studied & (units.subject == "Venus")].copy().reset_index(drop=True)
resp_rows, base, mean_sdf = [], [], {}
for c in cells.itertuples():
    s, u = S[c.si], S[c.si]["units"][c.ui]
    X = stretch.sdfs(u["pert"])
    base.append(X[:, stretch._mask(*stretch.BASELINE_MS)].mean())
    for d in ("flexion", "extension"):
        Xd = X[s["pert_direction"] == d]
        mean_sdf[(c.Index, d)] = Xd.mean(0)
        for method, shift in (("mean", 0), ("trials", 0), ("mean", stretch.PSEUDO_ONSET_MS),
                              ("trials", stretch.PSEUDO_ONSET_MS)):
            r = stretch.detect(Xd, shift_ms=shift, method=method)
            r.update(cell=c.Index, direction=d, method=method, pseudo=shift != 0)
            resp_rows.append(r)
cells["baseline_rate"] = base
R = pd.DataFrame(resp_rows).merge(cells[["state", "ptype", "si", "date", "duplicate_copy"]],
                                  left_on="cell", right_index=True)
main = R[(R.method == "mean") & ~R.pseudo]
by_cell = main.pivot(index="cell", columns="direction")
cells["torque_related"] = by_cell["responds"].any(axis=1).values
cells["n_directions"] = by_cell["responds"].sum(axis=1).values
cells["latency"] = by_cell["latency"].min(axis=1).values
dsis = [stretch.dsi(main[(main.cell == i) & (main.direction == "flexion")].iloc[0].to_dict(),
                    main[(main.cell == i) & (main.direction == "extension")].iloc[0].to_dict())
        for i in cells.index]
cells["dsi"] = [d for d, _ in dsis]
cells["category"] = [c for _, c in dsis]
resp = main[main.responds].copy()
print(f"{len(cells)} cells studied, {len(main)} cell-directions; responses found: {len(resp)} "
      f"in {int(cells.torque_related.sum())} cells")
flat = main[main.baseline_sd == 0]
print("Cell-directions without a usable baseline (SD of the baseline SDF 0; they can never respond):",
      len(flat), "| by type and state:", flat.groupby(["ptype", "state"]).size().to_dict())
print("Trials per cell and direction: ", main.n_trials.min(), "to", main.n_trials.max(),
      "| baseline rate range", round(cells.baseline_rate.min(), 2), "to", round(cells.baseline_rate.max(), 2), "spikes/s")
print("Response latencies (ms):", resp.latency.min(), "to", resp.latency.max(),
      "| offsets:", resp.offset.min(), "to", resp.offset.max())

271 cells studied, 542 cell-directions; responses found: 202 in 136 cells
Cell-directions without a usable baseline (SD of the baseline SDF 0; they can never respond): 14 | by type and state: {('CSN', 'post'): 3, ('CSN', 'pre'): 7, ('NA', 'post'): 3, ('NA', 'pre'): 1}
Trials per cell and direction:  5 to 28 | baseline rate range 0.0 to 48.96 spikes/s
Response latencies (ms): 20.0 to 60.0 | offsets: 24.0 to 600.0

counts = pd.crosstab(cells.ptype, cells.state).reindex(columns=["pre", "post"])
counts.loc["all"] = counts.sum()
print("Monkey V cells studied by type and state:")
print(counts.to_string())
dual_all = units[units.ptype == "dual"]
print(f"\nUnits activated from both the peduncle and the putamen, all files: {len(dual_all)} "
      f"({dual_all.subject.value_counts().to_dict()}); among the cells studied: {int((cells.ptype == 'dual').sum())}")
print("Other dual projections (with thalamus or STN):", int(units.projection.str.contains('_and_').sum() - len(dual_all)))
print("\nVariant, the two unit copies in repeated files dropped:")
print(pd.crosstab(cells[~cells.duplicate_copy].ptype, cells[~cells.duplicate_copy].state)
      .reindex(columns=["pre", "post"]).to_string())

Monkey V cells studied by type and state:
state  pre  post
ptype           
CSN     31    48
NA      30    56
PTN     49    54
dual     2     1
all    112   159

Units activated from both the peduncle and the putamen, all files: 3 ({'Venus': 3}); among the cells studied: 3
Other dual projections (with thalamus or STN): 4

Variant, the two unit copies in repeated files dropped:
state  pre  post
ptype           
CSN     31    48
NA      30    55
PTN     49    53
dual     2     1

def mwu(a, b):
    return stats.mannwhitneyu(a, b, alternative="two-sided").pvalue


expected_rate = {"M1": (12.6, 8.9, 1.98, 1.33), "PTN": (20.3, 14.8, 2.46, 2.72), "CSN": (3.5, 3.4, 1.47, 1.01)}
rows = []
for g, (e_pre, e_post, t_pre, t_post) in expected_rate.items():
    sub = cells if g == "M1" else cells[cells.ptype == g]
    pre, post = sub[sub.state == "pre"].baseline_rate, sub[sub.state == "post"].baseline_rate
    rows.append({"group": g, "n_pre": len(pre), "pre": pre.mean(), "sd_pre": pre.std(), "expected_pre": e_pre,
                 "in_tol_pre": abs(pre.mean() - e_pre) <= t_pre, "n_post": len(post), "post": post.mean(),
                 "sd_post": post.std(), "expected_post": e_post, "in_tol_post": abs(post.mean() - e_post) <= t_post,
                 "change_%": 100 * (post.mean() / pre.mean() - 1), "mwu_p": mwu(pre, post)})
rate_table = pd.DataFrame(rows)
print(rate_table.round(4).to_string(index=False))
nodup = cells[~cells.duplicate_copy]
p = nodup[nodup.ptype == "PTN"]
print(f"\nVariant without the unit copies, PTN: {p[p.state == 'pre'].baseline_rate.mean():.2f} and "
      f"{p[p.state == 'post'].baseline_rate.mean():.2f} spikes/s, p = {mwu(p[p.state == 'pre'].baseline_rate, p[p.state == 'post'].baseline_rate):.4f}")

group  n_pre     pre  sd_pre  expected_pre  in_tol_pre  n_post    post  sd_post  expected_post  in_tol_post  change_%  mwu_p
   M1    112 12.4222 10.5828          12.6        True     159  8.8714   9.4321            8.9         True  -28.5840 0.0033
  PTN     49 18.0317  9.2008          20.3        True      54 14.4832  10.1348           14.8         True  -19.6793 0.0343
  CSN     31  3.4201  4.0269           3.5        True      48  3.3647   4.1454            3.4         True   -1.6173 0.7253

Variant without the unit copies, PTN: 18.03 and 14.46 spikes/s, p = 0.0343

pre = cells[cells.state == "pre"]
ptn_pre, csn_pre = pre[pre.ptype == "PTN"].baseline_rate, pre[pre.ptype == "CSN"].baseline_rate
print(f"PTN {ptn_pre.mean():.2f} spikes/s (range {ptn_pre.min():.2f} to {ptn_pre.max():.2f}, n = {len(ptn_pre)}); "
      f"CSN {csn_pre.mean():.2f} (range {csn_pre.min():.2f} to {csn_pre.max():.2f}, n = {len(csn_pre)}); "
      f"Mann-Whitney p = {mwu(ptn_pre, csn_pre):.2e}")
print("Paper (both monkeys): text 16.3 and 3.6, Table 1 18.3 and 3.2 spikes/s")

PTN 18.03 spikes/s (range 2.56 to 48.96, n = 49); CSN 3.42 (range 0.00 to 20.73, n = 31); Mann-Whitney p = 1.25e-11
Paper (both monkeys): text 16.3 and 3.6, Table 1 18.3 and 3.2 spikes/s

def antidromic(df, label):
    rows = []
    for g, e in (("PTN", (1.9, 1.9, 0.22, 0.2)), ("CSN", (5.1, 6.1, 0.61, 0.89))):
        sub = df[df.ptype == g]
        a, b = sub[sub.state == "pre"].antidromic_latency_ms.dropna(), sub[sub.state == "post"].antidromic_latency_ms.dropna()
        rows.append({"cells": label, "group": g, "n_pre": len(a), "pre": a.mean(), "sd_pre": a.std(),
                     "in_tol_pre": abs(a.mean() - e[0]) <= e[2], "n_post": len(b), "post": b.mean(),
                     "sd_post": b.std(), "in_tol_post": abs(b.mean() - e[1]) <= e[3], "mwu_p": mwu(a, b)})
    return rows


leu_units = units[(units.subject == "Leu") & units.ptype.isin(["PTN", "CSN"])]
ad = pd.DataFrame(antidromic(cells, "monkey V studied")
                  + antidromic(pd.concat([cells, leu_units]), "variant: plus all monkey L units"))
print(ad.round(4).to_string(index=False))
print("Antidromic latency range (ms), cells studied:", cells.antidromic_latency_ms.min(), "to",
      cells.antidromic_latency_ms.max(), "| NaN for", int(cells[cells.ptype.isin(['PTN', 'CSN'])].antidromic_latency_ms.isna().sum()), "PTNs and CSNs")

                           cells group  n_pre    pre  sd_pre  in_tol_pre  n_post   post  sd_post  in_tol_post  mwu_p
                monkey V studied   PTN     47 1.9553  0.8285        True      48 1.8594   0.8072         True 0.5366
                monkey V studied   CSN     31 4.8919  1.7768        True      38 6.8474   3.7243         True 0.0180
variant: plus all monkey L units   PTN     63 1.8254  0.8021        True      55 1.8318   0.7730         True 0.9140
variant: plus all monkey L units   CSN     55 5.3827  2.3482        True      47 7.0149   3.4964        False 0.0071
Antidromic latency range (ms), cells studied: 0.65 to 20.5 | NaN for 18 PTNs and CSNs

fp = R[R.pseudo].groupby("method").agg(cell_directions=("responds", "size"),
                                       any_crossing_0_200=("any_crossing_0_200", "mean"),
                                       response_20_60=("responds", "mean"))
real = R[~R.pseudo].groupby("method").agg(any_crossing_0_200=("any_crossing_0_200", "mean"),
                                          response_20_60=("responds", "mean"))
print("At the pseudo onset (600 ms before the torque), fraction of cell-directions:")
print(fp.round(4).to_string())
print("\nAt the true onset, for reference:")
print(real.round(4).to_string())
print("\nPseudo onset, main method, by state:")
print(R[R.pseudo & (R.method == "mean")].groupby("state")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())
print(f"\nPseudo-onset over true-onset rate of responses (main): {fp.loc['mean', 'response_20_60'] / real.loc['mean', 'response_20_60']:.3f}")

At the pseudo onset (600 ms before the torque), fraction of cell-directions:
        cell_directions  any_crossing_0_200  response_20_60
method                                                     
mean                542              0.2196          0.0443
trials              542              0.4631          0.0923

At the true onset, for reference:
        any_crossing_0_200  response_20_60
method                                    
mean                0.6642          0.3727
trials              0.8653          0.3506

Pseudo onset, main method, by state:
       any_crossing_0_200  responds
state                              
post               0.2075    0.0472
pre                0.2366    0.0402

Pseudo-onset over true-onset rate of responses (main): 0.119

# Variant null: unperturbed trials, at the time a torque usually comes (the session's median delay
# from center target to torque onset), so the epoch matches the true one but has no torque.
null_rows = []
for c in cells.itertuples():
    s, u = S[c.si], S[c.si]["units"][c.ui]
    tr = s["trials"]
    delay = np.nanmedian(tr.onset - tr.center_on)
    t0 = (tr.center_on + delay)[tr.torque == "none"].dropna().values
    if len(t0) < stretch.MIN_REPS:
        continue
    rel = [np.round((u["spikes"][(u["spikes"] >= t - 1.8) & (u["spikes"] < t + 0.7)] - t) * 1000) for t in t0]
    for method in ("mean", "trials"):
        r = stretch.detect(stretch.sdfs(rel), method=method)
        null_rows.append({"cell": c.Index, "method": method, "n_trials": len(t0), "state": c.state,
                          "any_crossing_0_200": r["any_crossing_0_200"], "responds": r["responds"]})
null = pd.DataFrame(null_rows)
print(f"Unperturbed trials at the usual torque time: {null.cell.nunique()} cells with at least 5 such trials "
      f"({null.n_trials.min()} to {null.n_trials.max()} trials)")
print(null.groupby("method")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())

Unperturbed trials at the usual torque time: 253 cells with at least 5 such trials (5 to 27 trials)
        any_crossing_0_200  responds
method                              
mean                0.1937    0.0474
trials              0.4506    0.0751

# Where the false crossings come from. (1) The pseudo-onset rate by baseline rate and by number of
# trials. (2) Stationary Poisson spike trains at each cell-direction's baseline rate, with its number
# of trials, through the same detection: no response exists, so every crossing is false.
pm = R[R.pseudo & (R.method == "mean")].copy()
pm["rate_group"] = pd.qcut(pm.cell.map(cells.baseline_rate), 4)
pm["trials_group"] = pd.qcut(pm.n_trials, 3, duplicates="drop")
print("Pseudo onset, main detection, by baseline rate quartile (spikes/s):")
print(pm.groupby("rate_group", observed=True)[["any_crossing_0_200", "responds"]].agg(["mean", "size"]).round(3).to_string())
print("\nBy number of trials per direction:")
print(pm.groupby("trials_group", observed=True)[["any_crossing_0_200", "responds"]].agg(["mean", "size"]).round(3).to_string())

rng_fp = np.random.default_rng(1636)
span = int((stretch.SPIKE_WINDOW_S[1] - stretch.SPIKE_WINDOW_S[0]) * 1000)
sim = []
for r in main.itertuples():
    rate = cells.baseline_rate[r.cell]
    trains = [np.flatnonzero(rng_fp.random(span) < rate / 1000) + int(stretch.SPIKE_WINDOW_S[0] * 1000)
              for _ in range(r.n_trials)]
    X = stretch.sdfs(trains)
    for method in ("mean", "trials"):
        q = stretch.detect(X, method=method)
        sim.append({"method": method, "any_crossing_0_200": q["any_crossing_0_200"], "responds": q["responds"]})
sim = pd.DataFrame(sim)
print("\nStationary Poisson trains at each cell-direction's baseline rate (seeded), same detection:")
print(sim.groupby("method")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())

Pseudo onset, main detection, by baseline rate quartile (spikes/s):
                 any_crossing_0_200      responds     
                               mean size     mean size
rate_group                                            
(-0.001, 2.419]               0.169  136    0.044  136
(2.419, 6.968]                0.272  136    0.051  136
(6.968, 15.744]               0.194  134    0.060  134
(15.744, 48.965]              0.243  136    0.022  136

By number of trials per direction:
             any_crossing_0_200      responds     
                           mean size     mean size
trials_group                                      
(4.999, 7.0]              0.211  190    0.032  190
(7.0, 13.0]               0.220  264    0.045  264
(13.0, 28.0]              0.239   88    0.068   88

Stationary Poisson trains at each cell-direction's baseline rate (seeded), same detection:
        any_crossing_0_200  responds
method                              
mean                0.1531    0.0387
trials              0.4908    0.0978

# How the result changes with the threshold. The main detection crosses when the largest change over
# 0 to 200 ms exceeds t times the baseline SD of the mean SDF; the t that gives a pseudo-onset rate
# of 0.01 is read off the pseudo-onset distribution of that ratio and applied at the true onset.
z_pseudo = pm.max_abs_z.fillna(-np.inf).values
t_main = stats.t.ppf(1 - stretch.P_THRESHOLD / 2, int(np.diff(stretch.BASELINE_MS)[0]) - 1)
t_cal = float(np.quantile(z_pseudo, 0.99))
print(f"t of the main detection {t_main:.2f}; t giving a pseudo-onset rate of 0.01: {t_cal:.2f}")
cal_rows = []
for (cell, d), m in mean_sdf.items():
    for shift in (0, stretch.PSEUDO_ONSET_MS):
        q = stretch.detect(m[None, :], shift_ms=shift, tcrit=t_cal)
        cal_rows.append({"cell": cell, "direction": d, "pseudo": shift != 0, "responds": q["responds"],
                         "any_crossing_0_200": q["any_crossing_0_200"]})
cal = pd.DataFrame(cal_rows)
print(cal.groupby("pseudo")[["any_crossing_0_200", "responds"]].mean().round(4).to_string())
cal_cell = cal[~cal.pseudo].groupby("cell").responds.any()
first_overall = main.assign(early=main.first_crossing < stretch.ONSET_RANGE_MS[0])
print(f"\nDetected responses (main) with an earlier crossing at 0 to {stretch.ONSET_RANGE_MS[0] - 1} ms: "
      f"{int((first_overall.responds & first_overall.early).sum())} of {int(main.responds.sum())}")
fo_cell = first_overall.assign(r=first_overall.responds & ~first_overall.early).groupby("cell").r.any()
rows = []
for label, rel in (("main", cells.torque_related), ("calibrated threshold", cal_cell.reindex(cells.index)),
                   ("first crossing overall in 20 to 60 ms", fo_cell.reindex(cells.index))):
    for g in ("M1", "PTN", "CSN"):
        sub = cells.index if g == "M1" else cells.index[cells.ptype == g]
        for st in ("pre", "post"):
            idx = [k for k in sub if cells.state[k] == st]
            rows.append({"detection": label, "group": g, "state": st, "percent": 100 * rel[idx].mean(),
                         "expected": {('M1', 'pre'): 75, ('M1', 'post'): 68, ('PTN', 'pre'): 55, ('PTN', 'post'): 68, ('CSN', 'pre'): 45, ('CSN', 'post'): 21}[(g, st)]})
print("\nTorque-related percent by detection (paper's monkey V values for reference):")
print(pd.DataFrame(rows).pivot(index=["group", "state"], columns="detection", values="percent").round(1)
      .join(pd.DataFrame(rows).groupby(["group", "state"]).expected.first()).to_string())

# Example mean SDFs with their thresholds, as Figure 3 draws them: the cell-directions at the
# 10th, 50th and 90th percentile of baseline rate among detected responses and among non-responses.
pick = []
for flag in (True, False):
    g = main[main.responds == flag].sort_values("baseline")
    pick += [g.iloc[int(q * (len(g) - 1))] for q in (0.1, 0.5, 0.9)]
fig, axes = plt.subplots(2, 3, figsize=(10, 5), sharex=True)
for ax, r in zip(axes.ravel(), pick):
    m = mean_sdf[(r.cell, r.direction)]
    thr = t_main * r.baseline_sd
    ax.plot(stretch.GRID, m, color="k", lw=0.8)
    for y in (r.baseline - thr, r.baseline + thr):
        ax.axhline(y, color="0.5", lw=0.8, ls="--")
    ax.axvspan(stretch.PSEUDO_ONSET_MS, stretch.PSEUDO_ONSET_MS + 200, color="0.92")
    ax.axvspan(0, 200, color="0.85")
    ax.set_xlim(stretch.GRID[0], stretch.GRID[-1])
    ax.set_title(f"{cells.ptype[r.cell]} {r.state} {r.direction}, {r.baseline:.1f} sp/s, "
                 f"{'response' if r.responds else 'none'}", fontsize=8)
axes[1, 0].set_xlabel("ms from torque onset")
axes[0, 0].set_ylabel("mean SDF (spikes/s)")
plt.tight_layout()
plt.show()
for r in pick:
    m = mean_sdf[(r.cell, r.direction)]
    print(f"cell {r.cell} {r.direction}: SDF {m.min():.1f} to {m.max():.1f} spikes/s, threshold "
          f"{r.baseline - t_main * r.baseline_sd:.1f} to {r.baseline + t_main * r.baseline_sd:.1f}; "
          f"pseudo-onset crossing {bool(pm[(pm.cell == r.cell) & (pm.direction == r.direction)].any_crossing_0_200.iloc[0])}")
print("Shaded: 200 ms after the pseudo onset (light) and the torque (dark). Every SDF is within its axes.")

t of the main detection 3.11; t giving a pseudo-onset rate of 0.01: 6.67
        any_crossing_0_200  responds
pseudo                              
False               0.3376    0.2140
True                0.0111    0.0018

Detected responses (main) with an earlier crossing at 0 to 19 ms: 6 of 202

Torque-related percent by detection (paper's monkey V values for reference):
             calibrated threshold  first crossing overall in 20 to 60 ms  main  expected
group state                                                                             
CSN   post                    0.0                                   16.7  16.7        21
      pre                     3.2                                   16.1  16.1        45
M1    post                   32.7                                   46.5  47.2        68
      pre                    30.4                                   52.7  54.5        75
PTN   post                   42.6                                   64.8  66.7        68
      pre                    32.7                                   71.4  71.4        55

cell 189 flexion: SDF 0.0 to 35.1 spikes/s, threshold -2.2 to 4.5; pseudo-onset crossing False
cell 265 flexion: SDF 2.7 to 19.0 spikes/s, threshold 7.6 to 18.3; pseudo-onset crossing True
cell 77 extension: SDF 6.0 to 138.6 spikes/s, threshold 6.1 to 50.4; pseudo-onset crossing True
cell 232 extension: SDF 0.0 to 3.4 spikes/s, threshold -1.3 to 2.0; pseudo-onset crossing False
cell 4 flexion: SDF 0.0 to 14.8 spikes/s, threshold -5.6 to 15.5; pseudo-onset crossing False
cell 49 flexion: SDF 0.1 to 31.2 spikes/s, threshold 7.7 to 30.5; pseudo-onset crossing True
Shaded: 200 ms after the pseudo onset (light) and the torque (dark). Every SDF is within its axes.

expected_frac = {("M1", "pre"): 75, ("M1", "post"): 68, ("PTN", "pre"): 55, ("PTN", "post"): 68,
                 ("CSN", "pre"): 45, ("CSN", "post"): 21}


def related_table(df, label):
    rows = []
    for (g, st), e in expected_frac.items():
        sub = df if g == "M1" else df[df.ptype == g]
        sub = sub[sub.state == st]
        k, n = int(sub.torque_related.sum()), len(sub)
        rows.append({"detection": label, "group": g, "state": st, "related": k, "n": n,
                     "percent": 100 * k / n, "expected": e, "in_tol": abs(100 * k / n - e) <= 0.1 * e})
    return rows


# variant detection: the t-test across trials
var = R[(R.method == "trials") & ~R.pseudo].pivot(index="cell", columns="direction")
cells["torque_related_trials"] = var["responds"].any(axis=1).values
frac_table = pd.DataFrame(related_table(cells, "main")
                          + related_table(cells.assign(torque_related=cells.torque_related_trials), "t-test variant")
                          + related_table(cells[~cells.duplicate_copy], "main, unit copies dropped"))
print(frac_table.round(1).to_string(index=False))

                detection group state  related   n  percent  expected  in_tol
                     main    M1   pre       61 112     54.5        75   False
                     main    M1  post       75 159     47.2        68   False
                     main   PTN   pre       35  49     71.4        55   False
                     main   PTN  post       36  54     66.7        68    True
                     main   CSN   pre        5  31     16.1        45   False
                     main   CSN  post        8  48     16.7        21   False
           t-test variant    M1   pre       44 112     39.3        75   False
           t-test variant    M1  post       95 159     59.7        68   False
           t-test variant   PTN   pre       23  49     46.9        55   False
           t-test variant   PTN  post       39  54     72.2        68    True
           t-test variant   CSN   pre        8  31     25.8        45   False
           t-test variant   CSN  post       18  48     37.5        21   False
main, unit copies dropped    M1   pre       61 112     54.5        75   False
main, unit copies dropped    M1  post       74 157     47.1        68   False
main, unit copies dropped   PTN   pre       35  49     71.4        55   False
main, unit copies dropped   PTN  post       35  53     66.0        68    True
main, unit copies dropped   CSN   pre        5  31     16.1        45   False
main, unit copies dropped   CSN  post        8  48     16.7        21   False

def chi2_yates(a_k, a_n, b_k, b_n):
    table = np.array([[a_k, a_n - a_k], [b_k, b_n - b_k]])
    chi2, p, _, _ = stats.chi2_contingency(table, correction=True)
    return chi2, p


rows = []
for label, col in (("main", "torque_related"), ("t-test variant", "torque_related_trials")):
    for g in ("M1", "PTN", "CSN"):
        sub = cells if g == "M1" else cells[cells.ptype == g]
        a, b = sub[sub.state == "pre"][col], sub[sub.state == "post"][col]
        chi2, p = chi2_yates(a.sum(), len(a), b.sum(), len(b))
        rows.append({"detection": label, "group": g, "pre_%": 100 * a.mean(), "post_%": 100 * b.mean(),
                     "change_points": 100 * (b.mean() - a.mean()), "chi2": chi2, "p": p, "p_below_0.05/3": p < 0.05 / 3})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("\nThe paper's monkey V counts, through the same test:",
      {g: round(chi2_yates(*k)[1], 4) for g, k in {"M1": (84, 112, 108, 159), "PTN": (27, 49, 37, 54), "CSN": (14, 31, 10, 48)}.items()})

     detection group   pre_%  post_%  change_points    chi2      p  p_below_0.05/3
          main    M1 54.4643 47.1698        -7.2945  1.1221 0.2895           False
          main   PTN 71.4286 66.6667        -4.7619  0.0951 0.7578           False
          main   CSN 16.1290 16.6667         0.5376  0.0000 1.0000           False
t-test variant    M1 39.2857 59.7484        20.4627 10.2096 0.0014            True
t-test variant   PTN 46.9388 72.2222        25.2834  5.8391 0.0157            True
t-test variant   CSN 25.8065 37.5000        11.6935  0.6970 0.4038           False

The paper's monkey V counts, through the same test: {'M1': np.float64(0.26), 'PTN': np.float64(0.2307), 'CSN': np.float64(0.0408)}

for label, col in (("main", "torque_related"), ("t-test variant", "torque_related_trials")):
    p_, c_ = cells[cells.ptype == "PTN"][col], cells[cells.ptype == "CSN"][col]
    chi2, p = chi2_yates(p_.sum(), len(p_), c_.sum(), len(c_))
    print(f"{label}: PTN {int(p_.sum())}/{len(p_)} ({100 * p_.mean():.1f}%), CSN {int(c_.sum())}/{len(c_)} "
          f"({100 * c_.mean():.1f}%), chi-square {chi2:.1f}, p = {p:.2e}; expected 62.1% and 30.4%")

main: PTN 71/103 (68.9%), CSN 13/79 (16.5%), chi-square 47.5, p = 5.64e-12; expected 62.1% and 30.4%
t-test variant: PTN 62/103 (60.2%), CSN 26/79 (32.9%), chi-square 12.3, p = 4.64e-04; expected 62.1% and 30.4%

tr_cells = cells[cells.torque_related]
one = (tr_cells.n_directions == 1)
print(f"Torque-related cells responding in one direction only: {int(one.sum())}/{len(tr_cells)} "
      f"({100 * one.mean():.1f}%); by state: {tr_cells.groupby('state').apply(lambda d: round(100 * (d.n_directions == 1).mean(), 1)).to_dict()}")
vt = cells[cells.torque_related_trials]
n_dir_t = var["responds"].sum(axis=1).values[cells.torque_related_trials.values]
print(f"t-test variant: {int((n_dir_t == 1).sum())}/{len(vt)} ({100 * (n_dir_t == 1).mean():.1f}%)")

Torque-related cells responding in one direction only: 70/136 (51.5%); by state: {'post': 45.3, 'pre': 59.0}
t-test variant: 88/139 (63.3%)

rows = []
for g in ("M1", "PTN", "CSN"):
    sub = resp if g == "M1" else resp[resp.ptype == g]
    for st in ("pre", "post"):
        x = sub[sub.state == st]
        rows.append({"group": g, "state": st, "increases": int((x.sign > 0).sum()), "responses": len(x),
                     "percent": 100 * (x.sign > 0).mean(),
                     "flexion": f"{int(((x.sign > 0) & (x.direction == 'flexion')).sum())}/{int((x.direction == 'flexion').sum())}",
                     "extension": f"{int(((x.sign > 0) & (x.direction == 'extension')).sum())}/{int((x.direction == 'extension').sum())}"})
inc = pd.DataFrame(rows)
print(inc.round(1).to_string(index=False))
for g in ("PTN", "CSN"):
    a, b = inc[(inc.group == g) & (inc.state == "pre")].iloc[0], inc[(inc.group == g) & (inc.state == "post")].iloc[0]
    chi2, p = chi2_yates(a.increases, a.responses, b.increases, b.responses)
    print(f"{g} increases pre against post: chi-square {chi2:.2f}, p = {p:.3f}")
print("All responses, both states:", f"{100 * (resp.sign > 0).mean():.1f}% increases of {len(resp)}")
chi2, p = chi2_yates(23, 24, 19, 22)
print(f"Reference, the paper's own CSN counts in Table 1 (23/24 against 19/22): chi-square {chi2:.2f}, p = {p:.2f} "
      "(Table 1 marks this p < 0.01)")

group state  increases  responses  percent flexion extension
   M1   pre         71         86     82.6   39/47     32/39
   M1  post         89        116     76.7   48/64     41/52
  PTN   pre         32         47     68.1   17/25     15/22
  PTN  post         35         59     59.3   19/33     16/26
  CSN   pre          6          6    100.0     4/4       2/2
  CSN  post          8         10     80.0     5/6       3/4
PTN increases pre against post: chi-square 0.53, p = 0.467
CSN increases pre against post: chi-square 0.15, p = 0.696
All responses, both states: 79.2% increases of 202
Reference, the paper's own CSN counts in Table 1 (23/24 against 19/22): chi-square 0.38, p = 0.54 (Table 1 marks this p < 0.01)

measures = ["mean_change", "max_change", "auc"]
rows = []
for g in ("M1", "PTN"):
    sub = resp if g == "M1" else resp[resp.ptype == g]
    for d in ("both", "flexion", "extension"):
        x = sub if d == "both" else sub[sub.direction == d]
        a, b = x[x.state == "pre"], x[x.state == "post"]
        for m in measures:
            rows.append({"group": g, "direction": d, "measure": m, "n_pre": len(a), "n_post": len(b),
                         "mean_pre": a[m].mean(), "mean_post": b[m].mean(),
                         "ks_p": stats.ks_2samp(a[m], b[m]).pvalue, "mwu_p": mwu(a[m], b[m])})
mag = pd.DataFrame(rows)
print(mag.round(4).to_string(index=False))
main_tests = mag[mag.direction == "both"]
print(f"\nMain (directions pooled): smallest p {main_tests[['ks_p', 'mwu_p']].min().min():.4f}; "
      f"tests with p < 0.05: {int((main_tests[['ks_p', 'mwu_p']] < 0.05).sum().sum())} of 12")
for m in measures:
    print(f"{m} range: {resp[m].min():.2f} to {resp[m].max():.2f}")

group direction     measure  n_pre  n_post  mean_pre  mean_post   ks_p  mwu_p
   M1      both mean_change     86     116   27.3703    20.6158 0.0004 0.0005
   M1      both  max_change     86     116   43.6552    32.3955 0.0017 0.0024
   M1      both         auc     86     116    4.3185     1.9977 0.0254 0.0189
   M1   flexion mean_change     47      64   26.3613    20.5170 0.0060 0.0090
   M1   flexion  max_change     47      64   40.5192    31.4474 0.0105 0.0243
   M1   flexion         auc     47      64    4.7054     1.8180 0.1129 0.0605
   M1 extension mean_change     39      52   28.5863    20.7373 0.0287 0.0178
   M1 extension  max_change     39      52   47.4345    33.5623 0.0786 0.0313
   M1 extension         auc     39      52    3.8523     2.2189 0.2415 0.1368
  PTN      both mean_change     47      59   23.4471    15.3504 0.0007 0.0001
  PTN      both  max_change     47      59   33.3657    21.9125 0.0047 0.0013
  PTN      both         auc     47      59    3.3663     1.1540 0.0286 0.0932
  PTN   flexion mean_change     25      33   22.2380    14.2063 0.0102 0.0012
  PTN   flexion  max_change     25      33   31.6402    19.8632 0.0162 0.0041
  PTN   flexion         auc     25      33    4.5609     0.9730 0.0276 0.0662
  PTN extension mean_change     22      26   24.8211    16.8025 0.0394 0.0199
  PTN extension  max_change     22      26   35.3265    24.5135 0.1716 0.1000
  PTN extension         auc     22      26    2.0088     1.3837 0.7502 0.6415

Main (directions pooled): smallest p 0.0001; tests with p < 0.05: 11 of 12
mean_change range: 1.33 to 130.45
max_change range: 1.48 to 211.86
auc range: 0.01 to 25.46

expected_mag = {("M1", "flexion"): (42.4, 37.9), ("M1", "extension"): (50.9, 42.6),
                ("PTN", "flexion"): (33.4, 26.1), ("PTN", "extension"): (38.8, 29.3)}
rows = []
for reading, col in (("a max change", "max_change"), ("b mean rate", "mean_rate"), ("c mean change", "mean_change")):
    within = 0
    for (g, d), (e_pre, e_post) in expected_mag.items():
        sub = resp if g == "M1" else resp[resp.ptype == g]
        x = sub[sub.direction == d]
        for st, e in (("pre", e_pre), ("post", e_post)):
            v = x[x.state == st][col].mean()
            within += abs(v - e) <= 0.1 * e
            rows.append({"reading": reading, "group": g, "direction": d, "state": st, "value": v, "expected": e,
                         "in_tol": abs(v - e) <= 0.1 * e, "n": int((x.state == st).sum())})
mags = pd.DataFrame(rows)
print(mags.round(2).to_string(index=False))
print("\nWithin tolerance per reading (of 8):", mags.groupby("reading").in_tol.sum().to_dict())

      reading group direction state  value  expected  in_tol  n
 a max change    M1   flexion   pre  40.52      42.4    True 47
 a max change    M1   flexion  post  31.45      37.9   False 64
 a max change    M1 extension   pre  47.43      50.9    True 39
 a max change    M1 extension  post  33.56      42.6   False 52
 a max change   PTN   flexion   pre  31.64      33.4    True 25
 a max change   PTN   flexion  post  19.86      26.1   False 33
 a max change   PTN extension   pre  35.33      38.8    True 22
 a max change   PTN extension  post  24.51      29.3   False 26
  b mean rate    M1   flexion   pre  36.48      42.4   False 47
  b mean rate    M1   flexion  post  28.41      37.9   False 64
  b mean rate    M1 extension   pre  37.25      50.9   False 39
  b mean rate    M1 extension  post  28.80      42.6   False 52
  b mean rate   PTN   flexion   pre  32.30      33.4    True 25
  b mean rate   PTN   flexion  post  21.12      26.1   False 33
  b mean rate   PTN extension   pre  34.05      38.8   False 22
  b mean rate   PTN extension  post  26.05      29.3   False 26
c mean change    M1   flexion   pre  26.36      42.4   False 47
c mean change    M1   flexion  post  20.52      37.9   False 64
c mean change    M1 extension   pre  28.59      50.9   False 39
c mean change    M1 extension  post  20.74      42.6   False 52
c mean change   PTN   flexion   pre  22.24      33.4   False 25
c mean change   PTN   flexion  post  14.21      26.1   False 33
c mean change   PTN extension   pre  24.82      38.8   False 22
c mean change   PTN extension  post  16.80      29.3   False 26

Within tolerance per reading (of 8): {'a max change': 4, 'b mean rate': 1, 'c mean change': 0}

expected_lat = {("M1", "flexion"): (39.2, 34.56), ("M1", "extension"): (40.4, 33.4),
                ("PTN", "flexion"): (43.9, 35.5), ("PTN", "extension"): (44.4, 36.7)}


def latency_table(rsp):
    rows = []
    for (g, d), (e_pre, e_post) in expected_lat.items():
        sub = rsp if g == "M1" else rsp[rsp.ptype == g]
        x = sub[sub.direction == d]
        a, b = x[x.state == "pre"].latency, x[x.state == "post"].latency
        rows.append({"group": g, "direction": d, "n_pre": len(a), "pre": a.mean(), "expected_pre": e_pre,
                     "in_tol_pre": abs(a.mean() - e_pre) <= 0.1 * e_pre, "n_post": len(b), "post": b.mean(),
                     "expected_post": e_post, "in_tol_post": abs(b.mean() - e_post) <= 0.1 * e_post,
                     "ks_p": stats.ks_2samp(a, b).pvalue, "mwu_p": mwu(a, b)})
    return pd.DataFrame(rows)


lat = latency_table(resp)
print(lat.round(4).to_string(index=False))
print("\nPooled means (paper, both monkeys: M1 41 and 35 ms, PTN 46 and 38 ms):")
for g in ("M1", "PTN"):
    sub = resp if g == "M1" else resp[resp.ptype == g]
    cs = cells[cells.torque_related] if g == "M1" else cells[cells.torque_related & (cells.ptype == g)]
    print(f"  {g}: responses {sub[sub.state == 'pre'].latency.mean():.1f} and {sub[sub.state == 'post'].latency.mean():.1f} ms; "
          f"earliest per cell {cs[cs.state == 'pre'].latency.mean():.1f} and {cs[cs.state == 'post'].latency.mean():.1f} ms "
          f"(KS p = {stats.ks_2samp(cs[cs.state == 'pre'].latency, cs[cs.state == 'post'].latency).pvalue:.4f})")
resp_t = R[(R.method == "trials") & ~R.pseudo & R.responds]
print("\nt-test variant:")
print(latency_table(resp_t)[["group", "direction", "n_pre", "pre", "n_post", "post", "ks_p"]].round(4).to_string(index=False))

fig, axes = plt.subplots(1, 2, figsize=(9, 3.2), sharey=True)
for ax, g in zip(axes, ("M1", "PTN")):
    sub = resp if g == "M1" else resp[resp.ptype == g]
    for st, col in (("pre", "k"), ("post", "0.6")):
        for d, ls in (("flexion", "-"), ("extension", "--")):
            v = np.sort(sub[(sub.state == st) & (sub.direction == d)].latency)
            ax.step(v, np.arange(1, len(v) + 1) / len(v), color=col, ls=ls, where="post", label=f"{st} {d}")
    ax.set_title(g)
    ax.set_xlabel("response onset (ms)")
axes[0].set_ylabel("cumulative fraction")
axes[1].legend(fontsize=7)
plt.tight_layout()
plt.show()
print("Plotted latencies span", resp.latency.min(), "to", resp.latency.max(), "ms; none outside the axes.")

group direction  n_pre     pre  expected_pre  in_tol_pre  n_post    post  expected_post  in_tol_post   ks_p  mwu_p
   M1   flexion     47 40.7660          39.2        True      64 35.9375          34.56         True 0.0086 0.0144
   M1 extension     39 42.9487          40.4        True      52 35.8462          33.40         True 0.0048 0.0031
  PTN   flexion     25 43.5600          43.9        True      33 37.8485          35.50         True 0.0217 0.0366
  PTN extension     22 46.6364          44.4        True      26 37.7692          36.70         True 0.0023 0.0019

Pooled means (paper, both monkeys: M1 41 and 35 ms, PTN 46 and 38 ms):
  M1: responses 41.8 and 35.9 ms; earliest per cell 40.9 and 34.7 ms (KS p = 0.0032)
  PTN: responses 45.0 and 37.8 ms; earliest per cell 44.2 and 35.4 ms (KS p = 0.0026)

t-test variant:
group direction  n_pre     pre  n_post    post   ks_p
   M1   flexion     30 43.3667      67 40.3881 0.2666
   M1 extension     28 42.8929      65 41.3538 0.8709
  PTN   flexion     18 46.5000      27 42.0741 0.2647
  PTN extension     12 48.3333      25 44.0000 0.2789

Plotted latencies span 20.0 to 60.0 ms; none outside the axes.

tr_cells = cells[cells.torque_related].copy()
first_day = sessions[sessions.subject == "Venus"].date.min()
tr_cells["day"] = (tr_cells.date - first_day).dt.days
rho, p = stats.spearmanr(tr_cells.day, tr_cells.latency)
print(f"All torque-related cells (n = {len(tr_cells)}): Spearman rho = {rho:.3f}, p = {p:.4f}; days {tr_cells.day.min()} to {tr_cells.day.max()}")
for st in ("pre", "post"):
    x = tr_cells[tr_cells.state == st]
    rho, p = stats.spearmanr(x.day, x.latency)
    print(f"  {st}: n = {len(x)}, rho = {rho:.3f}, p = {p:.4f}")

All torque-related cells (n = 136): Spearman rho = -0.261, p = 0.0022; days 0 to 360
  pre: n = 61, rho = 0.105, p = 0.4195
  post: n = 75, rho = -0.145, p = 0.2147

expected_fwhm = {"flexion": (97, 52), "extension": (57, 41)}
ptn = resp[resp.ptype == "PTN"]
rows = []
for subset, x in (("increases", ptn[ptn.sign > 0]), ("all PTN responses", ptn)):
    for d, (e_pre, e_post) in expected_fwhm.items():
        y = x[x.direction == d]
        a, b = y[y.state == "pre"].fwhm, y[y.state == "post"].fwhm
        rows.append({"responses": subset, "direction": d, "n_pre": len(a), "pre": a.mean(),
                     "in_tol_pre": abs(a.mean() - e_pre) <= 0.1 * e_pre, "n_post": len(b), "post": b.mean(),
                     "in_tol_post": abs(b.mean() - e_post) <= 0.1 * e_post, "mwu_p": mwu(a, b)})
print(pd.DataFrame(rows).round(4).to_string(index=False))
inc = ptn[ptn.sign > 0]
print(f"\nPTN increases, directions pooled: {inc[inc.state == 'pre'].fwhm.mean():.1f} and {inc[inc.state == 'post'].fwhm.mean():.1f} ms "
      f"(paper 73 and 50), p = {mwu(inc[inc.state == 'pre'].fwhm, inc[inc.state == 'post'].fwhm):.4f}; FWHM range {ptn.fwhm.min():.0f} to {ptn.fwhm.max():.0f} ms")

        responses direction  n_pre     pre  in_tol_pre  n_post    post  in_tol_post  mwu_p
        increases   flexion     17 93.8235        True      19 40.7895        False 0.0394
        increases extension     15 56.9333        True      16 41.5625         True 0.3948
all PTN responses   flexion     25 97.4000        True      33 52.3333         True 0.0173
all PTN responses extension     22 56.8636        True      26 60.7308        False 0.8118

PTN increases, directions pooled: 76.5 and 41.1 ms (paper 73 and 50), p = 0.0486; FWHM range 17 to 382 ms

cats = ["reciprocal", "unidirectional", "bidirectional", "non-directional"]
tr_cells = cells[cells.torque_related]
tab = pd.crosstab(tr_cells.category, tr_cells.state).reindex(index=cats, columns=["pre", "post"]).fillna(0)
print(tab.to_string())
chi2, p, dof, _ = stats.chi2_contingency(tab.values)
high = tr_cells.category.isin(["reciprocal", "unidirectional"])
fr = high.groupby(tr_cells.state).mean()
print(f"\nReciprocal or unidirectional: pre {100 * fr['pre']:.1f}%, post {100 * fr['post']:.1f}% "
      f"(change {100 * (fr['post'] - fr['pre']):.1f} points); chi-square (2 x 4) {chi2:.2f}, dof {dof}, p = {p:.4f}")
for g in ("PTN", "CSN"):
    x = tr_cells[tr_cells.ptype == g]
    t = pd.crosstab(x.category, x.state).reindex(index=cats, columns=["pre", "post"]).fillna(0)
    t = t[t.sum(axis=1) > 0]
    c2, pp, _, _ = stats.chi2_contingency(t.values)
    h = x.category.isin(["reciprocal", "unidirectional"]).groupby(x.state).mean()
    print(f"  {g}: reciprocal or unidirectional pre {100 * h.get('pre', np.nan):.1f}%, post {100 * h.get('post', np.nan):.1f}%; chi-square p = {pp:.3f}")

# Variant: categories from the DSI value as the paper words them: reciprocal DSI > 1, unidirectional
# a response in one direction only (DSI <= 1), bidirectional 0.5 < DSI <= 1, non-directional DSI <= 0.5.
def by_value(r):
    if r.dsi > 1:
        return "reciprocal"
    if r.n_directions == 1:
        return "unidirectional"
    return "bidirectional" if r.dsi > 0.5 else "non-directional"


tr_v = tr_cells.assign(category_value=[by_value(r) for r in tr_cells.itertuples()])
tab_v = pd.crosstab(tr_v.category_value, tr_v.state).reindex(index=cats, columns=["pre", "post"]).fillna(0)
chi2_v, p_v, _, _ = stats.chi2_contingency(tab_v.values)
fr_v = tr_v.category_value.isin(["reciprocal", "unidirectional"]).groupby(tr_v.state).mean()
print("\nVariant, categories by DSI value:")
print(tab_v.to_string())
print(f"Reciprocal or unidirectional: pre {100 * fr_v['pre']:.1f}%, post {100 * fr_v['post']:.1f}%; chi-square {chi2_v:.2f}, p = {p_v:.4f}")
print("Smallest expected count in the main 2 x 4 table:", round(stats.chi2_contingency(tab.values)[3].min(), 1))

state            pre  post
category                  
reciprocal         0     5
unidirectional    36    34
bidirectional      7     7
non-directional   18    29

Reciprocal or unidirectional: pre 59.0%, post 52.0% (change -7.0 points); chi-square (2 x 4) 6.26, dof 3, p = 0.0998
  PTN: reciprocal or unidirectional pre 65.7%, post 50.0%; chi-square p = 0.015
  CSN: reciprocal or unidirectional pre 80.0%, post 75.0%; chi-square p = 1.000

Variant, categories by DSI value:
state            pre  post
category_value            
reciprocal        14    10
unidirectional    24    28
bidirectional      5     8
non-directional   18    29
Reciprocal or unidirectional: pre 62.3%, post 50.7%; chi-square 2.83, p = 0.4186
Smallest expected count in the main 2 x 4 table: 2.2

expected_dsi = {("M1", "pre"): (0.72, 0.087), ("M1", "post"): (0.66, 0.077), ("PTN", "pre"): (0.73, 0.115),
                ("PTN", "post"): (0.75, 0.132), ("CSN", "pre"): (0.81, 0.214), ("CSN", "post"): (0.83, 0.253)}
rows = []
for (g, st), (e, tol) in expected_dsi.items():
    x = tr_cells if g == "M1" else tr_cells[tr_cells.ptype == g]
    v = x[x.state == st].dsi
    rows.append({"group": g, "state": st, "n": len(v), "mean_dsi": v.mean(), "sd": v.std(), "expected": e,
                 "in_tol": abs(v.mean() - e) <= tol})
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("DSI range:", round(tr_cells.dsi.min(), 3), "to", round(tr_cells.dsi.max(), 3))

group state  n  mean_dsi    sd  expected  in_tol
   M1   pre 61     0.628 0.537      0.72   False
   M1  post 75     0.475 0.447      0.66   False
  PTN   pre 35     0.780 0.580      0.73    True
  PTN  post 36     0.496 0.487      0.75   False
  CSN   pre  5     0.662 0.516      0.81    True
  CSN  post  8     0.695 0.415      0.83    True
DSI range: 0.0 to 1.905

groups = [tr_cells[tr_cells.category == c].baseline_rate for c in cats]
f4, p4 = stats.f_oneway(*groups)
merged = tr_cells.category.replace({"reciprocal": "highly directional", "unidirectional": "highly directional"})
f3, p3 = stats.f_oneway(*[tr_cells[merged == c].baseline_rate for c in merged.unique()])
print("Baseline rate by category (mean, n):", {c: (round(g.mean(), 2), len(g)) for c, g in zip(cats, groups)})
print(f"One-way ANOVA, four categories: F = {f4:.2f}, p = {p4:.4f}; three groups: F = {f3:.2f}, p = {p3:.4f}")

Baseline rate by category (mean, n): {'reciprocal': (np.float64(15.93), 5), 'unidirectional': (np.float64(12.84), 70), 'bidirectional': (np.float64(12.71), 14), 'non-directional': (np.float64(15.35), 47)}
One-way ANOVA, four categories: F = 0.58, p = 0.6295; three groups: F = 0.69, p = 0.5016

kin_rows = []
for s in venus:
    if s["session_id"].split("++")[1] in set(repeats.repeated_in):
        continue
    for k, d in enumerate(s["pert_direction"]):
        r = stretch.kinematics(s["windows"]["angle"][k], d, s["analog_step_ms"])
        if r["valid"]:
            kin_rows.append({"session": s["session_id"], "post": s["post"], "direction": d, **r})
kin = pd.DataFrame(kin_rows)
kin["state"] = np.where(kin.post, "post", "pre")
print(f"{len(kin)} perturbations in {kin.session.nunique()} sessions; moving before onset: {int(kin.moving_before.sum())}; "
      f"no velocity crossing within 150 ms: {int(kin.latency_ms.isna().sum())}")
print("Ranges: latency", kin.latency_ms.min(), "to", kin.latency_ms.max(), "ms; peak velocity",
      round(kin.peak_velocity.min(), 1), "to", round(kin.peak_velocity.max(), 1), "deg/s; amplitude",
      round(kin.amplitude.min(), 2), "to", round(kin.amplitude.max(), 2), "deg")
sess_kin = kin.groupby(["session", "state", "direction"])[["latency_ms", "peak_velocity", "amplitude"]].mean().reset_index()
sess_kin["level"] = sess_kin.session.map(command.set_index("session").level)
print("\nSession means, mean over sessions:")
print(sess_kin.groupby(["direction", "state"])[["latency_ms", "peak_velocity", "amplitude"]].agg(["mean", "count"]).round(2).to_string())

5413 perturbations in 294 sessions; moving before onset: 490; no velocity crossing within 150 ms: 4
Ranges: latency 0.0 to 41.0 ms; peak velocity -5.0 to 572.4 deg/s; amplitude -1.66 to 13.88 deg

Session means, mean over sessions:
                latency_ms       peak_velocity       amplitude      
                      mean count          mean count      mean count
direction state                                                     
extension post        6.01   137         58.04   137      2.86   137
          pre        11.03   157         45.82   157      2.39   157
flexion   post        3.20   137         75.30   137      3.78   137
          pre         9.79   156         57.78   156      3.38   156

def anova(df, measure):
    model = ols(f"{measure} ~ C(state) * C(direction)", data=df).fit()
    return sm.stats.anova_lm(model, typ=2)


for m in ("peak_velocity", "amplitude"):
    print(f"{m}:")
    print(anova(sess_kin, m).round(4).to_string())
    print()
rows = []
for reading, x in (("a both directions", sess_kin), ("b flexion", sess_kin[sess_kin.direction == "flexion"]),
                   ("extension, for reference", sess_kin[sess_kin.direction == "extension"])):
    g = x.groupby("state")[["peak_velocity", "amplitude"]].mean()
    rows.append({"reading": reading, "velocity_change_%": 100 * (g.loc["post", "peak_velocity"] / g.loc["pre", "peak_velocity"] - 1),
                 "amplitude_change_%": 100 * (g.loc["post", "amplitude"] / g.loc["pre", "amplitude"] - 1)})
kin_change = pd.DataFrame(rows)
kin_change["velocity_in_tol"] = (kin_change["velocity_change_%"] + 16.2).abs() <= 1.62
kin_change["amplitude_in_tol"] = (kin_change["amplitude_change_%"] + 6.6).abs() <= 0.66
print(kin_change.round(2).to_string(index=False))
lvl2 = sess_kin[sess_kin.level == 2]
g = lvl2.groupby("state")[["peak_velocity", "amplitude"]].mean()
print(f"\nVariant, sessions with the middle command size only ({lvl2[lvl2.state == 'pre'].session.nunique()} pre, "
      f"{lvl2[lvl2.state == 'post'].session.nunique()} post): peak velocity {g.loc['pre', 'peak_velocity']:.1f} and "
      f"{g.loc['post', 'peak_velocity']:.1f} deg/s ({100 * (g.loc['post', 'peak_velocity'] / g.loc['pre', 'peak_velocity'] - 1):.1f}%), "
      f"amplitude {g.loc['pre', 'amplitude']:.2f} and {g.loc['post', 'amplitude']:.2f} deg "
      f"({100 * (g.loc['post', 'amplitude'] / g.loc['pre', 'amplitude'] - 1):.1f}%)")
for m in ("peak_velocity", "amplitude"):
    print(f"  {m}: MPTP F = {anova(lvl2, m).loc['C(state)', 'F']:.2f}, p = {anova(lvl2, m).loc['C(state)', 'PR(>F)']:.4f}")
print("\nBy command size (session means, both directions):")
print(sess_kin.groupby(["level", "state"])[["peak_velocity", "amplitude", "latency_ms"]].agg(["mean", "count"]).round(1).to_string())
fig, axes = plt.subplots(1, 3, figsize=(10, 3))
for ax, m in zip(axes, ("latency_ms", "peak_velocity", "amplitude")):
    for i, (d, st) in enumerate([(d, st) for d in ("flexion", "extension") for st in ("pre", "post")]):
        v = sess_kin[(sess_kin.direction == d) & (sess_kin.state == st)][m]
        ax.scatter(np.full(len(v), i) + np.random.default_rng(0).uniform(-0.15, 0.15, len(v)), v, s=6,
                   color="k" if st == "pre" else "0.6")
    ax.set_xticks(range(4), ["flex pre", "flex post", "ext pre", "ext post"], fontsize=7)
    ax.set_title(m)
plt.tight_layout()
plt.show()
print("All session means are plotted (ranges printed above).")
print(f"\nAbsolute scale: torque-induced amplitude, mean of session means {sess_kin.amplitude.mean():.2f} deg "
      f"(paper: about 10 deg); peak velocity {sess_kin.peak_velocity.mean():.1f} deg/s, highest session mean "
      f"{sess_kin.peak_velocity.max():.1f} (paper's matched sessions 130 to 200 deg/s). Active movements in the "
      f"same files, same filter: median peak velocity {act[act.subject == 'Venus'].active_peak_velocity.median():.1f} "
      f"deg/s and amplitude {act[act.subject == 'Venus'].active_amplitude.median():.1f} deg.")

peak_velocity:
                            sum_sq     df         F  PR(>F)
C(state)                32306.1073    1.0  103.3281  0.0000
C(direction)            30568.6125    1.0   97.7709  0.0000
C(state):C(direction)    1029.8107    1.0    3.2938  0.0701
Residual               182278.2433  583.0       NaN     NaN

amplitude:
                         sum_sq     df         F  PR(>F)
C(state)                27.6754    1.0   26.0982  0.0000
C(direction)           133.1410    1.0  125.5534  0.0000
C(state):C(direction)    0.1437    1.0    0.1355  0.7129
Residual               618.2328  583.0       NaN     NaN

                 reading  velocity_change_%  amplitude_change_%  velocity_in_tol  amplitude_in_tol
       a both directions              28.76               15.15            False             False
               b flexion              30.34               11.96            False             False
extension, for reference              26.67               19.49            False             False

Variant, sessions with the middle command size only (43 pre, 12 post): peak velocity 41.6 and 83.9 deg/s (101.9%), amplitude 2.14 and 4.58 deg (113.4%)
  peak_velocity: MPTP F = 42.80, p = 0.0000
  amplitude: MPTP F = 46.11, p = 0.0000

By command size (session means, both directions):
            peak_velocity       amplitude       latency_ms      
                     mean count      mean count       mean count
level state                                                     
0     pre            96.2     6       5.2     6        7.9     6
1     pre            54.2   216       3.1   216       11.7   216
2     post           83.9    24       4.6    24        7.7    24
      pre            41.6    85       2.1    85        7.8    85
4     post           65.0   250       3.2   250        4.3   250
      pre            66.4     6       3.8     6        3.8     6

All session means are plotted (ranges printed above).

Absolute scale: torque-induced amplitude, mean of session means 3.09 deg (paper: about 10 deg); peak velocity 58.7 deg/s, highest session mean 145.3 (paper's matched sessions 130 to 200 deg/s). Active movements in the same files, same filter: median peak velocity 99.3 deg/s and amplitude 20.4 deg.

print(anova(sess_kin, "latency_ms").round(4).to_string())
print(sess_kin.groupby(["direction", "state"]).latency_ms.mean().round(2).to_string())
a = anova(lvl2, "latency_ms")
print(f"Variant, middle command size only: MPTP F = {a.loc['C(state)', 'F']:.2f}, p = {a.loc['C(state)', 'PR(>F)']:.4f}; "
      f"means {lvl2.groupby('state').latency_ms.mean().round(2).to_dict()}")
print(f"\nPerturbations with a latency of 0 ms: {int((kin.latency_ms == 0).sum())} of {len(kin)}")
still = kin[~kin.moving_before]
for label, df, col in (("not moving before onset", still, "latency_ms"),
                       ("causal 10 ms average of the unfiltered velocity", kin, "latency_causal_ms"),
                       ("causal, not moving before onset", still, "latency_causal_ms")):
    sk = df.groupby(["session", "state", "direction"])[col].mean().reset_index()
    sk["level"] = sk.session.map(command.set_index("session").level)
    a = anova(sk.rename(columns={col: "lat"}), "lat")
    a2 = anova(sk[sk.level == 2].rename(columns={col: "lat"}), "lat")
    print(f"Variant, {label}: means {sk.groupby(['direction', 'state'])[col].mean().round(2).to_dict()}; "
          f"MPTP F = {a.loc['C(state)', 'F']:.1f}, p = {a.loc['C(state)', 'PR(>F)']:.2g}; "
          f"middle command size only F = {a2.loc['C(state)', 'F']:.2f}, p = {a2.loc['C(state)', 'PR(>F)']:.2g}")

                          sum_sq     df          F  PR(>F)
C(state)               4919.5365    1.0  1032.2752     0.0
C(direction)            571.0072    1.0   119.8155     0.0
C(state):C(direction)    91.0699    1.0    19.1094     0.0
Residual               2778.4159  583.0        NaN     NaN
direction  state
extension  post      6.01
           pre      11.03
flexion    post      3.20
           pre       9.79
Variant, middle command size only: MPTP F = 0.07, p = 0.7855; means {'post': 7.73, 'pre': 7.81}

Perturbations with a latency of 0 ms: 219 of 5413
Variant, not moving before onset: means {('extension', 'post'): 5.77, ('extension', 'pre'): 11.03, ('flexion', 'post'): 3.3, ('flexion', 'pre'): 9.79}; MPTP F = 1036.7, p = 1.8e-131; middle command size only F = 0.00, p = 0.97
Variant, causal 10 ms average of the unfiltered velocity: means {('extension', 'post'): 15.5, ('extension', 'pre'): 18.63, ('flexion', 'post'): 12.62, ('flexion', 'pre'): 17.61}; MPTP F = 781.3, p = 1e-109; middle command size only F = 1.33, p = 0.25
Variant, causal, not moving before onset: means {('extension', 'post'): 15.53, ('extension', 'pre'): 18.64, ('flexion', 'post'): 13.36, ('flexion', 'pre'): 17.6}; MPTP F = 742.8, p = 4.3e-106; middle command size only F = 5.12, p = 0.026

emg_rows, traces = [], {}
for s in venus:
    code_ = s["session_id"].split("++")[1]
    if code_ in set(repeats.repeated_in):
        continue
    for ch in stretch.EMG_CHANNELS[:3]:
        if ch not in s["windows"]:
            continue
        w = s["windows"][ch]
        if np.isnan(w).any() or np.nanstd(w[:, AX < 0]) < 1:
            continue
        for d in ("flexion", "extension"):
            r = stretch.emg_response(w[s["pert_direction"] == d], s["analog_step_ms"])
            emg_rows.append({"session": code_, "post": s["post"], "channel": ch, "direction": d,
                             **{k: r[k] for k in ("m1", "m2", "onset_ms", "duration_ms")}})
            traces[(code_, ch, d)] = r["trace"]
emg = pd.DataFrame(emg_rows)
emg["state"] = np.where(emg.post, "post", "pre")
emg["level"] = emg.session.map(command.set_index("code").level)
print("Channel-direction-sessions without perturbations of that direction (dropped):", int(emg.m2.isna().sum()))
emg = emg.dropna(subset=["m2"])
print("Channel-sessions used:")
print(pd.crosstab([emg.channel, emg.direction], emg.state).to_string())
print("\nMedian m2 (A/D units) by channel, direction and state:")
print(emg.groupby(["channel", "direction", "state"]).m2.median().unstack().round(2).to_string())
print("\nRanges: m1", round(emg.m1.min(), 1), "to", round(emg.m1.max(), 1), "; m2", round(emg.m2.min(), 1), "to",
      round(emg.m2.max(), 1), "; responses found in", int(emg.onset_ms.notna().sum()), "of", len(emg))
fig, axes = plt.subplots(1, 3, figsize=(11, 3), sharex=True)
for ax, ch in zip(axes, stretch.EMG_CHANNELS[:3]):
    for st, col in (("pre", "k"), ("post", "0.6")):
        keys = [(r.session, ch, "flexion") for r in emg[(emg.channel == ch) & (emg.state == st) & (emg.direction == "flexion")].itertuples()]
        if keys:
            m = np.mean([traces[k] for k in keys], axis=0)
            ax.plot(AX, m, color=col, label=f"{st} (n = {len(keys)})")
    ax.set_xlim(-50, 200)
    ax.axvspan(41, 80, color="0.9")
    ax.set_title(f"{ch}, flexion torque")
    ax.set_xlabel("ms from torque onset")
    ax.legend(fontsize=7)
axes[0].set_ylabel("EMG less baseline (A/D units)")
plt.tight_layout()
plt.show()
print("Grand averages are plotted from -50 to 200 ms (shaded: m2 window); y range is set by the data.")

Channel-direction-sessions without perturbations of that direction (dropped): 2
Channel-sessions used:
state                          post  pre
channel             direction           
EMGTricepsLateralis extension    36   36
                    flexion      36   36
EMGTricepsLongus1   extension   130  146
                    flexion     130  145
EMGTricepsLongus2   extension   116  130
                    flexion     116  129

Median m2 (A/D units) by channel, direction and state:
state                            post   pre
channel             direction              
EMGTricepsLateralis extension    5.73  1.47
                    flexion     76.18  5.20
EMGTricepsLongus1   extension   -2.55 -0.51
                    flexion    164.25  3.01
EMGTricepsLongus2   extension   -3.06 -0.19
                    flexion    297.52  3.14

Ranges: m1 -76.7 to 265.1 ; m2 -96.9 to 508.7 ; responses found in 1049 of 1186

Grand averages are plotted from -50 to 200 ms (shaded: m2 window); y range is set by the data.

rows = []
flex = emg[emg.direction == "flexion"]
for ch in stretch.EMG_CHANNELS[:3]:
    x = flex[flex.channel == ch]
    a, b = x[x.state == "pre"].m2, x[x.state == "post"].m2
    rows.append({"channel": ch, "n_pre": len(a), "median_pre": a.median(), "mean_pre": a.mean(), "n_post": len(b),
                 "median_post": b.median(), "mean_post": b.mean(), "ks_p": stats.ks_2samp(a, b).pvalue,
                 "larger_and_p<0.05": (b.median() > a.median()) and stats.ks_2samp(a, b).pvalue < 0.05})
m2_table = pd.DataFrame(rows)
print(m2_table.round({c: 2 for c in ("median_pre", "mean_pre", "median_post", "mean_post")})
      .to_string(index=False, formatters={"ks_p": "{:.2e}".format}))
print("\nVariant, middle command size only:")
for ch in stretch.EMG_CHANNELS[:3]:
    x = flex[(flex.channel == ch) & (flex.level == 2)]
    a, b = x[x.state == "pre"].m2, x[x.state == "post"].m2
    if len(a) and len(b):
        print(f"  {ch}: n {len(a)} and {len(b)}, median m2 {a.median():.2f} and {b.median():.2f}, KS p = {stats.ks_2samp(a, b).pvalue:.4f}")
print("\nMedian m2 by command size and state (flexion, all triceps channels):")
print(flex.groupby(["level", "state"]).m2.agg(["median", "count"]).round(2).to_string())

            channel  n_pre  median_pre  mean_pre  n_post  median_post  mean_post     ks_p  larger_and_p<0.05
  EMGTricepsLongus1    145        3.01     14.41     130       164.25     157.36 1.14e-58               True
  EMGTricepsLongus2    129        3.14     28.27     116       297.52     285.61 1.86e-50               True
EMGTricepsLateralis     36        5.20      9.54      36        76.18      76.35 4.65e-15               True

Variant, middle command size only:
  EMGTricepsLongus1: n 40 and 6, median m2 46.71 and 1.33, KS p = 0.0132
  EMGTricepsLongus2: n 36 and 4, median m2 49.52 and 124.79, KS p = 0.3200
  EMGTricepsLateralis: n 16 and 1, median m2 19.66 and 0.52, KS p = 0.1176

Median m2 by command size and state (flexion, all triceps channels):
             median  count
level state               
0     pre      3.68      9
1     pre      0.61    204
2     post     1.81     11
      pre     40.81     92
4     post   194.62    271
      pre     17.09      5

fr = flex.dropna(subset=["duration_ms"])
a, b = fr[fr.state == "pre"], fr[fr.state == "post"]
print(f"Responses found: pre {len(a)} of {int((flex.state == 'pre').sum())}, post {len(b)} of {int((flex.state == 'post').sum())} channel-sessions")
print(f"Duration: pre {a.duration_ms.mean():.1f} ms (SEM {a.duration_ms.sem():.1f}, median {a.duration_ms.median():.0f}), "
      f"post {b.duration_ms.mean():.1f} ms (SEM {b.duration_ms.sem():.1f}, median {b.duration_ms.median():.0f}); "
      f"KS p = {stats.ks_2samp(a.duration_ms, b.duration_ms).pvalue:.4g}; range {fr.duration_ms.min():.0f} to {fr.duration_ms.max():.0f} ms")
print(f"Onset: pre {a.onset_ms.mean():.1f} ms, post {b.onset_ms.mean():.1f} ms; KS p = {stats.ks_2samp(a.onset_ms, b.onset_ms).pvalue:.4g}; "
      f"range {fr.onset_ms.min():.0f} to {fr.onset_ms.max():.0f} ms")
l2 = fr[fr.level == 2]
print(f"Variant, middle command size only: duration {l2[l2.state == 'pre'].duration_ms.mean():.1f} and "
      f"{l2[l2.state == 'post'].duration_ms.mean():.1f} ms (n {int((l2.state == 'pre').sum())} and {int((l2.state == 'post').sum())}); "
      f"onset {l2[l2.state == 'pre'].onset_ms.mean():.1f} and {l2[l2.state == 'post'].onset_ms.mean():.1f} ms; "
      f"KS p duration {stats.ks_2samp(l2[l2.state == 'pre'].duration_ms, l2[l2.state == 'post'].duration_ms).pvalue:.4f}, "
      f"onset {stats.ks_2samp(l2[l2.state == 'pre'].onset_ms, l2[l2.state == 'post'].onset_ms).pvalue:.4f}")
print(f"Durations of 10 ms or less: pre {(a.duration_ms <= 10).mean():.1%}, post {(b.duration_ms <= 10).mean():.1%}; "
      f"median duration pre {a.duration_ms.median():.0f} ms, post {b.duration_ms.median():.0f} ms")
print("\nBy channel:")
print(fr.groupby(["channel", "state"])[["duration_ms", "onset_ms"]].agg(["mean", "count"]).round(1).to_string())
print("\nThreshold variants (pooled channels):")
for sd in (1.5, 3.0):
    rr = []
    post_of = dict(zip(sessions.code, sessions.post))
    for k, tr in traces.items():
        if k[2] != "flexion":
            continue
        st = "post" if post_of[k[0]] else "pre"
        b0 = (AX >= -200) & (AX < 0)
        thr = sd * tr[b0].std(ddof=1)
        above = (tr > thr) & (AX >= 10)
        cand = np.flatnonzero(above & (AX <= 150))
        if len(cand):
            rest = np.flatnonzero(~above[cand[0]:])
            rr.append({"state": st, "duration": int(rest[0]) if len(rest) else np.nan, "onset": AX[cand[0]]})
    rr = pd.DataFrame(rr).dropna()
    pa, pb = rr[rr.state == "pre"], rr[rr.state == "post"]
    print(f"  {sd} SD: duration {pa.duration.mean():.1f} and {pb.duration.mean():.1f} ms (KS p = {stats.ks_2samp(pa.duration, pb.duration).pvalue:.4g}); "
          f"onset {pa.onset.mean():.1f} and {pb.onset.mean():.1f} ms (KS p = {stats.ks_2samp(pa.onset, pb.onset).pvalue:.4g})")

Responses found: pre 283 of 310, post 281 of 282 channel-sessions
Duration: pre 24.5 ms (SEM 3.8, median 9), post 50.4 ms (SEM 2.4, median 59); KS p = 1.757e-39; range 1 to 491 ms
Onset: pre 40.8 ms, post 24.0 ms; KS p = 5.527e-21; range 10 to 147 ms
Variant, middle command size only: duration 28.4 and 27.9 ms (n 92 and 10); onset 28.8 and 45.6 ms; KS p duration 0.9631, onset 0.1828
Durations of 10 ms or less: pre 63.6%, post 25.6%; median duration pre 9 ms, post 59 ms

By channel:
                          duration_ms       onset_ms      
                                 mean count     mean count
channel             state                                 
EMGTricepsLateralis post         72.9    36     24.8    36
                    pre          16.3    33     44.5    33
EMGTricepsLongus1   post         47.3   130     23.3   130
                    pre          27.8   128     40.8   128
EMGTricepsLongus2   post         46.7   115     24.5   115
                    pre          23.2   122     39.7   122

Threshold variants (pooled channels):
  1.5 SD: duration 16.4 and 48.5 ms (KS p = 2.279e-36); onset 33.6 and 22.3 ms (KS p = 5.383e-14)
  3.0 SD: duration 17.3 and 49.5 ms (KS p = 1.988e-43); onset 54.4 and 25.9 ms (KS p = 7.121e-38)

for st in ("pre", "post"):
    x = flex[flex.state == st]
    ratio = x.m1 / x.m2
    print(f"{st}: m2 > m1 in {100 * (x.m2 > x.m1).mean():.1f}% of {len(x)} channel-sessions; median m1 {x.m1.median():.2f}, "
          f"median m2 {x.m2.median():.2f}, median m1/m2 {ratio.median():.2f} (over channel-sessions with m2 > 0: "
          f"{(x.m1 / x.m2)[x.m2 > 0].median():.2f})")

pre: m2 > m1 in 73.5% of 310 channel-sessions; median m1 0.38, median m2 3.18, median m1/m2 0.10 (over channel-sessions with m2 > 0: 0.09)
post: m2 > m1 in 96.5% of 282 channel-sessions; median m1 40.65, median m2 191.67, median m1/m2 0.22 (over channel-sessions with m2 > 0: 0.23)

sv = sess_kin.set_index(["session", "direction"]).peak_velocity
resp_v = resp.copy()
resp_v["session"] = [S[si]["session_id"] for si in resp_v.si]
resp_v["session_velocity"] = [sv.get((s, d), np.nan) for s, d in zip(resp_v.session, resp_v.direction)]
print("Responses without a session velocity (repeated files):", int(resp_v.session_velocity.isna().sum()))
matched = resp_v[(resp_v.session_velocity >= 130) & (resp_v.session_velocity <= 200)]
rows = []
for d in ("flexion", "extension"):
    x = matched[matched.direction == d]
    a, b = x[x.state == "pre"], x[x.state == "post"]
    row = {"direction": d, "n_pre": len(a), "n_post": len(b),
           "velocity_pre": a.session_velocity.mean(), "velocity_post": b.session_velocity.mean(),
           "velocity_p": mwu(a.session_velocity, b.session_velocity)}
    for m in ("max_change", "mean_change", "auc"):
        row[f"{m}_pre"], row[f"{m}_post"], row[f"{m}_p"] = a[m].mean(), b[m].mean(), mwu(a[m], b[m])
    rows.append(row)
print(pd.DataFrame(rows).round(4).T.to_string())
print("Session mean peak velocity range:", round(sess_kin.peak_velocity.min(), 1), "to", round(sess_kin.peak_velocity.max(), 1), "deg/s")

Responses without a session velocity (repeated files): 2
                         0          1
direction          flexion  extension
n_pre                    0          0
n_post                   1          0
velocity_pre           NaN        NaN
velocity_post     133.1665        NaN
velocity_p             NaN        NaN
max_change_pre         NaN        NaN
max_change_post    93.5684        NaN
max_change_p           NaN        NaN
mean_change_pre        NaN        NaN
mean_change_post   56.2528        NaN
mean_change_p          NaN        NaN
auc_pre                NaN        NaN
auc_post           10.8005        NaN
auc_p                  NaN        NaN
Session mean peak velocity range: 8.5 to 145.3 deg/s

ptn_cells = cells[cells.torque_related & (cells.ptype == "PTN")].copy()
ptn_cells["code"] = [S[si]["session_id"].split("++")[1] for si in ptn_cells.si]
ptn_cells["session"] = [S[si]["session_id"] for si in ptn_cells.si]
per_session = ptn_cells.groupby(["session", "code"]).latency.mean().reset_index()
vel = sess_kin.groupby("session").peak_velocity.mean()
m2 = flex[flex.channel == "EMGTricepsLongus1"].set_index("session").m2
per_session["velocity"] = per_session.session.map(vel)
per_session["m2"] = per_session.code.map(m2)
print(f"Sessions with a torque-related PTN: {len(per_session)}; with velocity {int(per_session.velocity.notna().sum())}; "
      f"with triceps longus 1 m2 {int(per_session.m2.notna().sum())}")
for col in ("velocity", "m2"):
    x = per_session.dropna(subset=[col])
    rho, p = stats.spearmanr(x.latency, x[col])
    print(f"PTN latency against {col}: Spearman rho = {rho:.3f}, p = {p:.4f} (n = {len(x)})")
    for st, sub in x.groupby(x.session.str.contains("PostMPTP")):
        r2, p2 = stats.spearmanr(sub.latency, sub[col])
        print(f"   {'post' if st else 'pre'} only: rho = {r2:.3f}, p = {p2:.4f} (n = {len(sub)})")

Sessions with a torque-related PTN: 68; with velocity 67; with triceps longus 1 m2 64
PTN latency against velocity: Spearman rho = -0.216, p = 0.0795 (n = 67)
   pre only: rho = 0.021, p = 0.9080 (n = 33)
   post only: rho = 0.001, p = 0.9973 (n = 34)
PTN latency against m2: Spearman rho = -0.316, p = 0.0109 (n = 64)
   pre only: rho = 0.108, p = 0.5549 (n = 32)
   post only: rho = 0.166, p = 0.3627 (n = 32)

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

Notebook run time: 3.1 min

Claim	Tested
`cell-counts`: 112 and 159 monkey V cells, 49 and 54 PTNs, 31 and 48 CSNs, 3 dual	yes
`monkey-l-perturbations`: monkey L torque responses	blocked (data missing)
`baseline-rate`: lower resting rate after MPTP in PTNs, not CSNs	yes
`ptn-faster-than-csn`: PTNs fire faster than CSNs before MPTP	yes
`antidromic-latency`: CSN antidromic latency longer after MPTP, PTN unchanged	yes
`false-positive-rate`: the detection threshold's false-positive rate near 0.01	yes
`torque-related-fractions`: monkey V fractions of torque-related cells	yes
`prevalence-unchanged`: MPTP does not change the fraction of torque-related cells	yes
`ptn-more-responsive`: PTNs respond more often than CSNs	yes
`one-direction-half`: 44% of torque-related cells respond to one direction only	yes
`excitatory-majority`: most responses are increases; PTN balance unchanged	yes
`magnitude-unchanged`: response magnitudes unchanged by MPTP	yes
`venus-response-magnitudes`: monkey V response magnitudes	yes
`latency-shorter`: response latencies shorter after MPTP	yes
`latency-not-training`: latencies do not follow days of training	yes
`fwhm-shorter-ptn`: PTN increases briefer after MPTP	yes
`dsi-reduced`: fewer directional responses after MPTP	yes
`venus-dsi`: monkey V mean DSI	yes
`dsi-not-baseline`: resting rate does not differ between directional categories	yes
`kinematics-slower`: torque-induced movements slower and smaller after MPTP	yes
`kinematic-latency-unchanged`: movement latency unchanged	yes
`emg-m2-larger`: larger long-latency (m2) triceps response after MPTP	yes
`emg-duration-longer`: longer EMG responses after MPTP	yes
`emg-onset-unchanged`: EMG onset unchanged	yes
`emg-m1-weak`: m1 much smaller than m2	yes
`velocity-matched-magnitude`: magnitudes unchanged in velocity-matched sessions	yes
`session-correlations`: no session-level correlation of PTN latency with velocity or m2	yes
`trial-correlations`: no trial-level correlation	blocked (methods underspecified)

Claim	Outcome
`cell-counts`	reproduced
`monkey-l-perturbations`	blocked: data missing
`baseline-rate`	reproduced
`ptn-faster-than-csn`	reproduced
`antidromic-latency`	partially reproduced
`false-positive-rate`	blocked: methods underspecified
`torque-related-fractions`	blocked: methods underspecified
`prevalence-unchanged`	blocked: methods underspecified
`ptn-more-responsive`	blocked: methods underspecified
`one-direction-half`	blocked: methods underspecified
`excitatory-majority`	blocked: methods underspecified
`magnitude-unchanged`	blocked: methods underspecified
`venus-response-magnitudes`	blocked: methods underspecified
`latency-shorter`	blocked: methods underspecified
`latency-not-training`	blocked: methods underspecified
`fwhm-shorter-ptn`	blocked: methods underspecified
`dsi-reduced`	blocked: methods underspecified
`venus-dsi`	blocked: methods underspecified
`dsi-not-baseline`	blocked: methods underspecified
`kinematics-slower`	not reproduced
`kinematic-latency-unchanged`	not reproduced
`emg-m2-larger`	reproduced
`emg-duration-longer`	partially reproduced
`emg-onset-unchanged`	not reproduced
`emg-m1-weak`	partially reproduced
`velocity-matched-magnitude`	blocked: methods underspecified
`session-correlations`	blocked: methods underspecified
`trial-correlations`	blocked: methods underspecified

Reproducing Pasquereau & Turner (2013) from dandiset 001636¶

Data check¶

Responses of the cells studied¶

`cell-counts`¶

`monkey-l-perturbations`¶

`baseline-rate`¶

`ptn-faster-than-csn`¶

`antidromic-latency`¶

`false-positive-rate`¶

`prevalence-unchanged`¶

`ptn-more-responsive`¶

`one-direction-half`¶

`excitatory-majority`¶

`magnitude-unchanged` and `venus-response-magnitudes`¶

`latency-shorter`¶

`latency-not-training`¶

`fwhm-shorter-ptn`¶

`dsi-reduced`, `venus-dsi` and `dsi-not-baseline`¶

Kinematics and EMG¶

`kinematics-slower`¶

`kinematic-latency-unchanged`¶

EMG¶

`emg-m2-larger`¶

`emg-duration-longer` and `emg-onset-unchanged`¶

`emg-m1-weak`¶

`velocity-matched-magnitude`¶

`session-correlations`¶

`trial-correlations`¶

Summary¶

Reproducing Pasquereau & Turner (2013) from dandiset 001636¶

Data check¶

Responses of the cells studied¶

cell-counts¶

monkey-l-perturbations¶

baseline-rate¶

ptn-faster-than-csn¶

antidromic-latency¶

false-positive-rate¶

torque-related-fractions¶

prevalence-unchanged¶

ptn-more-responsive¶

one-direction-half¶

excitatory-majority¶

magnitude-unchanged and venus-response-magnitudes¶

latency-shorter¶

latency-not-training¶

fwhm-shorter-ptn¶

dsi-reduced, venus-dsi and dsi-not-baseline¶

Kinematics and EMG¶

kinematics-slower¶

kinematic-latency-unchanged¶

EMG¶

emg-m2-larger¶

emg-duration-longer and emg-onset-unchanged¶

emg-m1-weak¶

velocity-matched-magnitude¶

session-correlations¶

trial-correlations¶

Summary¶

`cell-counts`¶

`monkey-l-perturbations`¶

`baseline-rate`¶

`ptn-faster-than-csn`¶

`antidromic-latency`¶

`false-positive-rate`¶

`torque-related-fractions`¶

`prevalence-unchanged`¶

`ptn-more-responsive`¶

`one-direction-half`¶

`excitatory-majority`¶

`magnitude-unchanged` and `venus-response-magnitudes`¶

`latency-shorter`¶

`latency-not-training`¶

`fwhm-shorter-ptn`¶

`dsi-reduced`, `venus-dsi` and `dsi-not-baseline`¶

`kinematics-slower`¶

`kinematic-latency-unchanged`¶

`emg-m2-larger`¶

`emg-duration-longer` and `emg-onset-unchanged`¶

`emg-m1-weak`¶

`velocity-matched-magnitude`¶

`session-correlations`¶

`trial-correlations`¶