import itertools
import time
import warnings

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

import prepspace as ps

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

started = time.time()
data = ps.load_train()
trials, units, spikes, hand = data["trials"], data["units"], data["spikes"], data["hand"]
print(f"loaded dandiset {ps.DANDISET_ID} version {ps.DANDISET_VERSION}, {ps.TRAIN_PATH}, "
      f"in {time.time() - started:.0f} s")

loaded dandiset 000138 version 0.220113.0407, sub-Jenkins/sub-Jenkins_ses-large_desc-train_behavior+ecephys.nwb, in 4 s

api = "https://api.dandiarchive.org/api/dandisets/000138"
versions = requests.get(f"{api}/versions/").json()["results"]
print("Versions:", [(v["version"], v["asset_count"], v["size"]) for v in versions])
for v in ["0.220113.0407", "draft"]:
    assets = requests.get(f"{api}/versions/{v}/assets/").json()["results"]
    print(v, [(a["path"], a["size"], a["blob"]) for a in assets])
meta = requests.get(f"{api}/versions/{ps.DANDISET_VERSION}/").json()
print("\nRelated resources:", [(r["identifier"], r["relation"]) for r in meta["relatedResource"]])
print("Description:", meta["description"])
print("\nFile metadata:", data["meta"])

Versions: [('draft', 2, 149392888), ('0.220113.0407', 2, 149392888)]

0.220113.0407 [('sub-Jenkins/sub-Jenkins_ses-large_desc-test_ecephys.nwb', 802352, 'd78a051b-669f-4b7c-bcb6-769f106afb1c'), ('sub-Jenkins/sub-Jenkins_ses-large_desc-train_behavior+ecephys.nwb', 148590536, 'c7549202-993a-42a0-b8b3-27635e053d60')]

draft [('sub-Jenkins/sub-Jenkins_ses-large_desc-test_ecephys.nwb', 802352, 'd78a051b-669f-4b7c-bcb6-769f106afb1c'), ('sub-Jenkins/sub-Jenkins_ses-large_desc-train_behavior+ecephys.nwb', 148590536, 'c7549202-993a-42a0-b8b3-27635e053d60')]

Related resources: [('10.1016/j.neuron.2010.09.015', 'dcite:IsDescribedBy'), ('10.1038/nature11129', 'dcite:IsDescribedBy')]
Description: This dataset contains sorted unit spiking times and behavioral data from a macaque performing a delayed reaching task. The experimental task was a center-out reaching task with obstructing barriers forming a maze, resulting in a variety of straight and curved reaches. Neural activity was recorded from electrode arrays implanted in the motor cortex (M1) and dorsal premotor cortex (PMd). Cursor position, hand position, and eye position were also recorded during the experiment, and hand velocity was calculated offline from hand position. The provided data has been limited to 500 train trials and 100 test trials. Provided as part of the Neural Latents Benchmark: https://neurallatents.github.io.

File metadata: {'session_start': '2009-10-06T00:00:00-07:00', 'subject': 'Jenkins', 'species': 'Macaca mulatta', 'session_description': 'Data from monkey Jenkins performing center-out delayed reaching task. This file contains trials from the full session on 2009-10-06 that can be used for training models for the Neural Latents Benchmark.'}

print("Trials:", len(trials), "; success:", trials.success.value_counts().to_dict(),
      "; split:", trials.split.value_counts().to_dict())
print("Conditions (trial_type x trial_version):", trials.condition.nunique(),
      "; maze_id values:", trials.maze_id.nunique())
print(trials.groupby(["trial_type", "trial_version"]).size().unstack().to_string())
check = trials.groupby("condition").agg(maze_ids=("maze_id", "nunique"),
                                        barriers=("num_barriers", lambda x: sorted(set(x))),
                                        targets=("num_targets", lambda x: sorted(set(x))),
                                        active=("active_target", lambda x: sorted(set(x))))
print("\nPer condition, distinct maze_id, number of barriers, targets, active target index:")
print(check.to_string())
print("\ndelay field minus (go cue - target on) in ms, max abs:",
      float(np.abs(trials.delay - (trials.go_cue_time - trials.target_on_time) * 1000).max()))
print("rt field minus (move onset - go cue) in ms, max abs:",
      float(np.abs(trials.rt - (trials.move_onset_time - trials.go_cue_time) * 1000).max()))
print("Delay (ms) range:", trials.delay.min(), "to", trials.delay.max(),
      "; trials with delay > 400 ms:", int((trials.delay > 400).sum()))
print("Trials with delay of 0 ms or less:", int((trials.delay <= 0).sum()))
print("RT (ms) range:", trials.rt.min(), "to", trials.rt.max())
print("Target on after trial start (s):", round((trials.target_on_time - trials.start_time).min(), 3),
      "to", round((trials.target_on_time - trials.start_time).max(), 3))
print("Trial stop after movement onset (s):", round((trials.stop_time - trials.move_onset_time).min(), 3),
      "to", round((trials.stop_time - trials.move_onset_time).max(), 3))

Trials: 500 ; success: {True: 500} ; split: {'train': 375, 'val': 125}
Conditions (trial_type x trial_version): 27 ; maze_id values: 9
trial_version   0   1   2
trial_type               
2              18  19  19
3              18  18  19
4              18  18  18
6              18  18  19
7              19  18  19
8              18  18  19
10             18  19  19
11             19  19  19
12             19  18  19

Per condition, distinct maze_id, number of barriers, targets, active target index:
           maze_ids barriers targets active
condition                                  
(2, 0)            1      [0]     [1]    [0]
(2, 1)            1      [6]     [1]    [0]
(2, 2)            1      [6]     [3]    [1]
(3, 0)            1      [0]     [1]    [0]
(3, 1)            1      [6]     [1]    [0]
(3, 2)            1      [6]     [3]    [0]
(4, 0)            1      [0]     [1]    [0]
(4, 1)            1      [6]     [1]    [0]
(4, 2)            1      [6]     [3]    [2]
(6, 0)            1      [0]     [1]    [0]
(6, 1)            1      [9]     [1]    [0]
(6, 2)            1      [9]     [3]    [0]
(7, 0)            1      [0]     [1]    [0]
(7, 1)            1      [9]     [1]    [0]
(7, 2)            1      [9]     [3]    [1]
(8, 0)            1      [0]     [1]    [0]
(8, 1)            1      [9]     [1]    [0]
(8, 2)            1      [9]     [3]    [2]
(10, 0)           1      [0]     [1]    [0]
(10, 1)           1      [8]     [1]    [0]
(10, 2)           1      [8]     [3]    [0]
(11, 0)           1      [0]     [1]    [0]
(11, 1)           1      [8]     [1]    [0]
(11, 2)           1      [8]     [3]    [2]
(12, 0)           1      [0]     [1]    [0]
(12, 1)           1      [8]     [1]    [0]
(12, 2)           1      [8]     [3]    [1]

delay field minus (go cue - target on) in ms, max abs: 2.1100277081131935e-10
rt field minus (move onset - go cue) in ms, max abs: 2.219167072325945e-10
Delay (ms) range: -2 to 999 ; trials with delay > 400 ms: 401
Trials with delay of 0 ms or less: 8
RT (ms) range: 218 to 1200
Target on after trial start (s): 0.65 to 0.961
Trial stop after movement onset (s): 0.963 to 1.54

el = data["electrodes"]
print("Electrodes:", len(el), "; by group:", el.group_name.value_counts().to_dict(),
      "; x, y, z all NaN:", bool(el[["x", "y", "z"]].isna().all().all()))
print("Units:", len(units), "; heldout:", units.heldout.value_counts().to_dict())
print("Location of the electrode each unit references:", units.electrode_location.value_counts().to_dict())
print("Unit id = array * 1000 + channel * 10 + sort: id prefix counts",
      units.id_array.value_counts().sort_index().to_dict())
print("Channel in the id equals the referenced electrode index + 1 for every unit:",
      bool((units.id_channel == units.electrode + 1).all()))
print("Heldout units by id prefix:", units.groupby("id_array").heldout.sum().to_dict())
n_spk = np.array([len(s) for s in spikes])
print("Spikes per unit:", n_spk.min(), "to", n_spk.max(), "; total", n_spk.sum())
print("Spike times sorted in every unit:", all(np.all(np.diff(s) >= 0) for s in spikes))
starts, stops = trials.start_time.values, trials.stop_time.values
allsp = np.concatenate(spikes)
k = np.searchsorted(starts, allsp, side="right") - 1
print("Spikes outside the trial intervals:", int(np.sum((k < 0) | (allsp > stops[np.clip(k, 0, None)]))))
dur = (stops - starts).sum()
rates = n_spk / dur
print(f"Mean rate per unit over the trial intervals (spikes/s): {rates.min():.2f} to {rates.max():.2f}, "
      f"median {np.median(rates):.2f}")

Electrodes: 192 ; by group: {'electrode_group_PMd': 96, 'electrode_group_M1': 96} ; x, y, z all NaN: True
Units: 162 ; heldout: {False: 122, True: 40}
Location of the electrode each unit references: {'PMd': 162}
Unit id = array * 1000 + channel * 10 + sort: id prefix counts {1: 83, 2: 79}
Channel in the id equals the referenced electrode index + 1 for every unit: True
Heldout units by id prefix: {1: 20, 2: 20}
Spikes per unit: 197 to 26207 ; total 916892
Spike times sorted in every unit: True

Spikes outside the trial intervals: 0
Mean rate per unit over the trial intervals (spikes/s): 0.14 to 18.03, median 2.97

hp, hv, ht = hand["hand_pos"], hand["hand_vel"], hand["t"]
for name in ["hand_pos", "hand_vel", "cursor_pos", "eye_pos"]:
    x = hand[name]
    print(f"{name}: stored unit {hand[name + '_unit']!r}, conversion {hand[name + '_conversion']}, "
          f"description {hand[name + '_description']!r}; x {x[:, 0].min():.1f} to {x[:, 0].max():.1f}, "
          f"y {x[:, 1].min():.1f} to {x[:, 1].max():.1f}")
dt = np.diff(ht)
print("Hand samples:", len(ht), "; interval within trials", np.median(dt) * 1000, "ms; samples outside trials:",
      int(np.sum(ht > stops[np.searchsorted(starts, ht, side='right') - 1] + 1e-9)))
# velocity against the derivative of position, within trials
within = np.r_[False, np.isclose(dt, 0.001)]
dpos = np.r_[np.zeros((1, 2)), np.diff(hp, axis=0)] * 1000
fast = within & (np.hypot(hv[:, 0], hv[:, 1]) > 100)
for d, ax in enumerate("xy"):
    r = np.corrcoef(dpos[fast, d], hv[fast, d])[0, 1]
    slope = np.polyfit(dpos[fast, d], hv[fast, d], 1)[0]
    print(f"hand_vel {ax} against the derivative of hand_pos {ax} (samples above 100 mm/s): r = {r:.3f}, slope {slope:.3f}")
# reach endpoint against the active target
end_err = []
for _, tr in trials.iterrows():
    i = np.searchsorted(ht, tr.move_onset_time + 0.75)
    target = np.asarray(tr.target_pos)[int(tr.active_target)]
    end_err.append(np.hypot(*(hp[i] - target)))
end_err = np.array(end_err)
print(f"Hand 750 ms after movement onset to the active target (mm): median {np.median(end_err):.1f}, "
      f"95th percentile {np.percentile(end_err, 95):.1f}, max {end_err.max():.1f}")

hand_pos: stored unit 'meters', conversion 0.001, description 'Monkey hand position in mm. First column x, second column y'; x -144.7 to 165.3, y -149.8 to 78.1
hand_vel: stored unit 'm/s', conversion 0.001, description 'Monkey hand velocity in mm/s, calculated from monkey hand position. First column x, second column y'; x -960.2 to 1080.2, y -984.6 to 1091.2
cursor_pos: stored unit 'meters', conversion 0.001, description 'Cursor position in mm. First column x, second column y'; x -146.4 to 167.5, y -116.9 to 115.0
eye_pos: stored unit 'meters', conversion 0.001, description 'Monkey eye position in mm. First column x, second column y'; x -290.1 to 220.3, y -430.6 to 350.7

Hand samples: 1453340 ; interval within trials 0.9999999999763531 ms; samples outside trials: 0

hand_vel x against the derivative of hand_pos x (samples above 100 mm/s): r = 1.000, slope 1.000
hand_vel y against the derivative of hand_pos y (samples above 100 mm/s): r = 1.000, slope 1.000
Hand 750 ms after movement onset to the active target (mm): median 36.3, 95th percentile 45.8, max 50.4

def decode(t0_s, window):
    i0 = np.round((t0_s + window[0] / 1000) * 1000).astype(int)
    counts = np.zeros((len(trials), len(spikes)))
    for u, s in enumerate(spikes):
        ms = np.floor(s * 1000).astype(int)
        cs = np.r_[0, np.cumsum(np.bincount(ms, minlength=int(stops.max() * 1000) + 3000))]
        counts[:, u] = cs[i0 + window[1] - window[0]] - cs[i0]
    X = np.sqrt(counts)
    y = trials.trial_type.values
    labels = np.unique(y)
    sums = np.stack([X[y == l].sum(0) for l in labels])
    ns = np.array([(y == l).sum() for l in labels])
    correct = 0
    for i in range(len(y)):
        li = np.searchsorted(labels, y[i])
        cent = sums.copy()
        cent[li] -= X[i]
        nn = ns.copy()
        nn[li] -= 1
        cent = cent / nn[:, None]
        correct += labels[np.argmin(((cent - X[i]) ** 2).sum(1))] == y[i]
    return correct / len(y)

clock = []
for name, ev, win in [("target onset", trials.target_on_time.values, (-250, -50)),
                      ("target onset", trials.target_on_time.values, (100, 300)),
                      ("go cue", trials.go_cue_time.values, (-200, 0)),
                      ("movement onset", trials.move_onset_time.values, (0, 200))]:
    clock.append({"event": name, "window_ms": win, "accuracy": decode(ev, win),
                  "accuracy_window_1s_earlier": decode(ev - 1.0, win) if name == "target onset" else np.nan})
print("Chance for 9 mazes about", round(1 / 9, 3))
print(pd.DataFrame(clock).round(3).to_string(index=False))
speed_on = []
for tm in trials.move_onset_time.values:
    i = np.searchsorted(ht, tm - 0.0005)
    speed_on.append(np.hypot(*hv[i - 100:i + 101].T))
speed_on = np.array(speed_on)
print("Median hand speed (mm/s) at -100, -50, 0, +50, +100 ms from the table's movement onset:",
      np.round(np.median(speed_on, 0)[[0, 50, 100, 150, 200]], 1))

Chance for 9 mazes about 0.111
         event   window_ms  accuracy  accuracy_window_1s_earlier
  target onset (-250, -50)     0.112                       0.088
  target onset  (100, 300)     0.548                       0.106
        go cue   (-200, 0)     0.562                         NaN
movement onset    (0, 200)     0.774                         NaN
Median hand speed (mm/s) at -100, -50, 0, +50, +100 ms from the table's movement onset: [  2.3   5.1  44.8 397.7 723.9]

test = ps.load_test_summary()
print("Test file trial columns:", test["trial_columns"], "; trials:", test["n_trials"],
      "; units:", len(test["unit_ids"]),
      "; unit ids all among the training file's held-in units:",
      set(test["unit_ids"]) == set(units.id[~units.heldout]))
print("Segment durations (s):", np.unique(np.round(test["stops"] - test["starts"], 3)))
hits, tried = 0, 0
for uid, st in zip(test["unit_ids"], test["spikes"]):
    train_isi = np.round(np.diff(spikes[int(np.nonzero(units.id.values == uid)[0][0])]) * 1e4).astype(int)
    train_keys = {tuple(train_isi[i:i + 5]) for i in range(len(train_isi) - 4)}
    isi = np.round(np.diff(st) * 1e4).astype(int)
    for i in range(0, len(isi) - 4, 5):
        tried += 1
        hits += tuple(isi[i:i + 5]) in train_keys
print(f"Units checked: {len(test['unit_ids'])}; runs of 5 test intervals found in the training intervals "
      f"of the same unit: {hits} of {tried}")

Test file trial columns: ['id', 'move_onset_time', 'split', 'start_time', 'stop_time'] ; trials: 100 ; units: 122 ; unit ids all among the training file's held-in units: True
Segment durations (s): [0.7]

Units checked: 122; runs of 5 test intervals found in the training intervals of the same unit: 0 of 8734

onsets = {"a": ps.onset_times(data, "a"), "b": ps.onset_times(data, "b")}
d_on = (onsets["b"] - onsets["a"]) * 1000
print(f"Onset (b) minus onset (a), ms: median {np.median(d_on):.0f}, 5th to 95th percentile "
      f"{np.percentile(d_on, 5):.0f} to {np.percentile(d_on, 95):.0f}, range {d_on.min():.0f} to {d_on.max():.0f}")
t0 = time.time()
R = ps.trial_rates(data, onsets)
print(f"Smoothed rates computed in {time.time() - t0:.0f} s")
ctxs = {}
rows = []
for proto, combine, onset in itertools.product(ps.READINGS["proto"], ps.READINGS["combine"], ps.READINGS["onset"]):
    ctx = ps.prepare(data, R, onsets, proto, combine, onset)
    ctxs[(proto, combine, onset)] = ctx
    per = np.bincount(trials.cond_idx[ctx["mask"]], minlength=27)
    row = {"proto": proto, "combine": combine, "onset": onset, "kept": int(ctx["mask"].sum()),
           "of_delay_over_400": int((trials.delay > 400).sum()),
           "fraction": ctx["mask"].sum() / (trials.delay > 400).sum(),
           "min_per_condition": per.min(), "max_per_condition": per.max()}
    for tun in ps.READINGS["tuning"]:
        tp = ps.tuning(ctx["prep"], tun)
        row[f"included_tuning_{tun}"] = int((tp > ps.ARRAY_MIN_TUNING).sum())
    rows.append(row)
selection = pd.DataFrame(rows)
print(selection.round(3).to_string(index=False))
per_cond = pd.DataFrame({k: np.bincount(trials.cond_idx[c["mask"]], minlength=27) for k, c in ctxs.items()}).T
per_cond.index.names = ["proto", "combine", "onset"]
print("\nTrials kept per condition (columns: condition index):")
print(per_cond.to_string())

Onset (b) minus onset (a), ms: median 14, 5th to 95th percentile 7 to 24, range -74 to 39

Smoothed rates computed in 17 s

proto combine onset  kept  of_delay_over_400  fraction  min_per_condition  max_per_condition  included_tuning_a  included_tuning_b
    a  concat     a   386                401     0.963                 11                 18                 92                  6
    a  concat     b   386                401     0.963                 11                 18                 92                  7
    a    each     a   325                401     0.810                  3                 17                101                  7
    a    each     b   332                401     0.828                  3                 18                 95                  7
    b  concat     a   380                401     0.948                 10                 18                 91                  6
    b  concat     b   382                401     0.953                 11                 18                 92                  6
    b    each     a   317                401     0.791                  2                 16                103                  7
    b    each     b   326                401     0.813                  3                 18                 97                  7

Trials kept per condition (columns: condition index):
                     0   1   2   3   4   5   6   7   8   9   10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26
proto combine onset                                                                                                            
a     concat  a      15  14  17  13  13  14  14  15  14  15  17  12  15  13  13  14  11  16  13  14  18  14  11  16  12  16  17
              b      15  13  17  13  13  13  14  16  14  15  17  12  15  14  13  14  12  16  13  15  18  14  11  15  12  16  16
      each    a       9  13  14  12  13  10   6  14  14  13  14  11  14  13   8  14  10  14  12  12  17  11  10  15   3  15  14
              b       8  11  15  12  13   9   7  15  14  14  15  11  14  13  10  14  10  15  13  11  18  12   9  14   3  16  16
b     concat  a      15  14  17  13  13  14  14  15  14  15  17  12  14  14  13  14  10  16  13  14  18  14  10  16  11  15  15
              b      15  12  17  13  13  13  14  16  14  15  17  12  15  14  13  14  11  16  13  14  18  14  11  15  11  16  16
      each    a       7  13  14  12  13  10   5  14  14  13  14  11  14  12   8  14  10  14  13  11  16  11  10  15   2  13  14
              b       7  11  14  12  13   9   6  15  14  14  15  11  14  13   9  14  10  13  13  10  18  13   9  15   3  16  15

ctx = ctxs[("a", "concat", "a")]
tp_a, tp_b = ps.tuning(ctx["prep"], "a"), ps.tuning(ctx["prep"], "b")
print(f"Preparatory condition means (spikes/s): {ctx['prep'].min():.2f} to {ctx['prep'].max():.2f}")
print(f"Preparatory tuning, range (spikes/s): {tp_a.min():.2f} to {tp_a.max():.2f}, median {np.median(tp_a):.2f}")
print(f"Preparatory tuning, SD (spikes/s): {tp_b.min():.2f} to {tp_b.max():.2f}, median {np.median(tp_b):.2f}")
peri = ps.epoch_mean(ctx["move"], ps.PERI_WINDOW)
print(f"Peri-movement condition means (spikes/s): {peri.min():.2f} to {peri.max():.2f}")
inc_a = tp_a > ps.ARRAY_MIN_TUNING
print("Included units (range reading) by id prefix:", pd.Series(units.id_array[inc_a]).value_counts().sort_index().to_dict())
print("Median preparatory tuning (range) by id prefix:",
      pd.Series(tp_a).groupby(units.id_array.values).median().round(2).to_dict())
print("Median peri-movement tuning (range) by id prefix:",
      pd.Series(peri.max(1) - peri.min(1)).groupby(units.id_array.values).median().round(2).to_dict())
print(f"Hand position, condition means (mm): x {ctx['pos'][..., 0].min():.0f} to {ctx['pos'][..., 0].max():.0f}, "
      f"y {ctx['pos'][..., 1].min():.0f} to {ctx['pos'][..., 1].max():.0f}; "
      f"peak speed {np.hypot(ctx['vel'][..., 0], ctx['vel'][..., 1]).max():.0f} mm/s")

Preparatory condition means (spikes/s): 0.00 to 35.11
Preparatory tuning, range (spikes/s): 0.24 to 34.89, median 5.64
Preparatory tuning, SD (spikes/s): 0.07 to 10.30, median 1.45
Peri-movement condition means (spikes/s): 0.00 to 36.44
Included units (range reading) by id prefix: {1: 55, 2: 37}
Median preparatory tuning (range) by id prefix: {1: 6.71, 2: 4.79}
Median peri-movement tuning (range) by id prefix: {1: 10.079999923706055, 2: 6.0}
Hand position, condition means (mm): x -129 to 145, y -132 to 68; peak speed 960 mm/s

print(ps.MAIN)

{'proto': 'a', 'combine': 'concat', 'onset': 'a', 'tuning': 'a', 't_units': 'all', 'fold': 'once', 'norm': 'neuron', 'kin_end': 'fixed', 'task_end': 'speed', 'split': 'length', 'kmax': 14}

def fig4_selection(ctx, tun, onset):
    prep = ctx["prep"]
    peri = ps.epoch_mean(ctx["move"], ps.PERI_WINDOW)
    tp, tq = ps.tuning(prep, tun), ps.tuning(peri, tun)
    inc = (tp > ps.ARRAY_MIN_TUNING) & (tq >= ps.PERI_RATIO * tp)
    return inc, prep, peri

rows = []
fig4 = {}
for (proto, combine, onset), ctx in ctxs.items():
    for tun in ps.READINGS["tuning"]:
        inc, prep, peri = fig4_selection(ctx, tun, onset)
        row = {"proto": proto, "combine": combine, "onset": onset, "tuning": tun,
               "n_units": int(inc.sum())}
        if inc.sum() >= 3:
            r, rs = ps.corr_pearson_spearman(prep[inc], peri[inc])
            fig4[(proto, combine, onset, tun)] = (inc, r, rs)
            row.update({"mean_r": r.mean(), "mean_spearman": rs.mean(),
                        "p_t": stats.ttest_1samp(r, 0).pvalue,
                        "p_t_fisher": stats.ttest_1samp(np.arctanh(np.clip(r, -0.999999, 0.999999)), 0).pvalue,
                        "p_wilcoxon": stats.wilcoxon(r).pvalue})
            row["r_within_tol"] = abs(row["mean_r"] - 0.13) <= 0.013
            row["spearman_within_tol"] = abs(row["mean_spearman"] - 0.11) <= 0.011
            row["positive_p05"] = row["mean_r"] > 0 and min(row["p_t"], row["p_t_fisher"], row["p_wilcoxon"]) < 0.05
        rows.append(row)
fig4d = pd.DataFrame(rows)
print("Expected: mean r 0.13 +/- 0.013, mean Spearman 0.11 +/- 0.011, mean r > 0 at p < 0.05")
print(fig4d.round(4).to_string(index=False))
main4 = fig4d[(fig4d.proto == "a") & (fig4d["combine"] == "concat") & (fig4d.onset == "a") & (fig4d.tuning == "a")].iloc[0]
print("\nMain reading:", {k: (round(v, 4) if isinstance(v, float) else v) for k, v in main4.to_dict().items()})
ok = fig4d.dropna(subset=["mean_r"]).astype({"r_within_tol": bool, "spearman_within_tol": bool, "positive_p05": bool})
print("Readings with a result:", len(ok), "; both means within tolerance:",
      int((ok.r_within_tol & ok.spearman_within_tol).sum()), "; mean r positive at p < 0.05:", int(ok.positive_p05.sum()))
print("Range of mean r over readings with a result:", round(ok.mean_r.min(), 3), "to", round(ok.mean_r.max(), 3),
      "; mean Spearman:", round(ok.mean_spearman.min(), 3), "to", round(ok.mean_spearman.max(), 3))

Expected: mean r 0.13 +/- 0.013, mean Spearman 0.11 +/- 0.011, mean r > 0 at p < 0.05
proto combine onset tuning  n_units  mean_r  mean_spearman    p_t  p_t_fisher  p_wilcoxon r_within_tol spearman_within_tol positive_p05
    a  concat     a      a       20  0.1954         0.2119 0.0233      0.0225      0.0240        False               False         True
    a  concat     a      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    a  concat     b      a       19  0.1115         0.1528 0.1802      0.1659      0.1688        False               False        False
    a  concat     b      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    a    each     a      a       22  0.1339         0.1651 0.0845      0.0795      0.0738         True               False        False
    a    each     a      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    a    each     b      a       20  0.1240         0.1609 0.1403      0.1254      0.1429         True               False        False
    a    each     b      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    b  concat     a      a       22  0.1353         0.1601 0.1072      0.0961      0.0854         True               False        False
    b  concat     a      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    b  concat     b      a       19  0.1045         0.1465 0.2046      0.1909      0.1688        False               False        False
    b  concat     b      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    b    each     a      a       24  0.1392         0.1774 0.0772      0.0692      0.0491         True               False         True
    b    each     a      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN
    b    each     b      a       24  0.1224         0.1508 0.1086      0.0958      0.1140         True               False        False
    b    each     b      b        0     NaN            NaN    NaN         NaN         NaN          NaN                 NaN          NaN

Main reading: {'proto': 'a', 'combine': 'concat', 'onset': 'a', 'tuning': 'a', 'n_units': 20, 'mean_r': 0.1954, 'mean_spearman': 0.2119, 'p_t': 0.0233, 'p_t_fisher': 0.0225, 'p_wilcoxon': 0.024, 'r_within_tol': False, 'spearman_within_tol': False, 'positive_p05': True}
Readings with a result: 8 ; both means within tolerance: 0 ; mean r positive at p < 0.05: 2
Range of mean r over readings with a result: 0.105 to 0.195 ; mean Spearman: 0.146 to 0.212

inc, r, rs = fig4[("a", "concat", "a", "a")]
print("Main reading, per-unit Pearson r: range", round(r.min(), 3), "to", round(r.max(), 3),
      "; by id prefix:", pd.Series(r).groupby(units.id_array.values[inc]).agg(["count", "mean"]).round(3).to_dict())
# how the result moves with the peri-movement criterion
for ratio in [1.0, 1.25, 1.5, 2.0]:
    ctx = ctxs[("a", "concat", "a")]
    prep = ctx["prep"]; peri = ps.epoch_mean(ctx["move"], ps.PERI_WINDOW)
    tp, tq = ps.tuning(prep, "a"), ps.tuning(peri, "a")
    sel = (tp > 5) & (tq >= ratio * tp)
    rr, rrs = ps.corr_pearson_spearman(prep[sel], peri[sel])
    print(f"peri/prep tuning ratio >= {ratio}: {sel.sum()} units, mean r {rr.mean():.3f}, "
          f"mean Spearman {rrs.mean():.3f}, t-test p {stats.ttest_1samp(rr, 0).pvalue:.4f}")
fig, ax = plt.subplots(figsize=(4, 3))
ax.hist(r, bins=np.linspace(-1, 1, 21), color="0.5")
ax.axvline(r.mean(), color="r")
ax.set_xlabel("Pearson r, preparatory vs peri-movement")
ax.set_ylabel("units")
ax.set_title("Main reading (all values inside the axes)")
plt.show()

Main reading, per-unit Pearson r: range -0.459 to 0.794 ; by id prefix: {'count': {1: 12, 2: 8}, 'mean': {1: 0.209, 2: 0.176}}
peri/prep tuning ratio >= 1.0: 52 units, mean r 0.137, mean Spearman 0.132, t-test p 0.0041
peri/prep tuning ratio >= 1.25: 35 units, mean r 0.106, mean Spearman 0.114, t-test p 0.0687
peri/prep tuning ratio >= 1.5: 20 units, mean r 0.195, mean Spearman 0.212, t-test p 0.0233
peri/prep tuning ratio >= 2.0: 12 units, mean r 0.202, mean Spearman 0.231, t-test p 0.0450

tt = np.arange(*ps.MOVE_TRACE)
rows = []
curves = {}
for key, (inc, r, rs) in fig4.items():
    ctx = ctxs[key[:3]]
    prep = ctx["prep"][inc]
    trace = ctx["move"][inc]
    pc = prep - prep.mean(1, keepdims=True)
    tc = trace - trace.mean(1, keepdims=True)
    curve = np.nanmean((pc[:, :, None] * tc).sum(1) / np.sqrt((pc ** 2).sum(1)[:, None] * (tc ** 2).sum(1)), 0)
    curves[key] = curve
    at = {ms: curve[ms - ps.MOVE_TRACE[0]] for ms in [-400, -300, -250, -200, -100, 0]}
    rows.append({"proto": key[0], "combine": key[1], "onset": key[2], "tuning": key[3], "n_units": int(inc.sum()),
                 **{f"r_at_{ms}": v for ms, v in at.items()},
                 "near_zero_at_0": abs(at[0]) < 0.1, "below_half_of_-250": at[0] < 0.5 * at[-250],
                 "no_early_decline": at[-300] >= 0.75 * at[-400]})
fig4e = pd.DataFrame(rows)
fig4e["criterion_met"] = fig4e["near_zero_at_0"] & fig4e["below_half_of_-250"] & fig4e["no_early_decline"]
print(fig4e.round(3).to_string(index=False))
print("Readings meeting the criterion:", int(fig4e.criterion_met.sum()), "of", len(fig4e))
main_curve = curves[("a", "concat", "a", "a")]
print("Main reading, mean correlation every 50 ms from -600 to +400 ms:")
print(pd.Series(main_curve[::50], index=tt[::50]).round(3).to_string())

proto combine onset tuning  n_units  r_at_-400  r_at_-300  r_at_-250  r_at_-200  r_at_-100  r_at_0  near_zero_at_0  below_half_of_-250  no_early_decline  criterion_met
    a  concat     a      a       20      0.773      0.764      0.713      0.572      0.374   0.218           False                True              True          False
    a  concat     b      a       19      0.762      0.734      0.653      0.533      0.304   0.133           False                True              True          False
    a    each     a      a       22      0.735      0.718      0.634      0.521      0.309   0.163           False                True              True          False
    a    each     b      a       20      0.733      0.706      0.603      0.510      0.324   0.146           False                True              True          False
    b  concat     a      a       22      0.760      0.743      0.687      0.549      0.326   0.170           False                True              True          False
    b  concat     b      a       19      0.761      0.741      0.649      0.533      0.299   0.133           False                True              True          False
    b    each     a      a       24      0.716      0.700      0.639      0.492      0.283   0.151           False                True              True          False
    b    each     b      a       24      0.722      0.683      0.559      0.483      0.286   0.129           False                True              True          False
Readings meeting the criterion: 0 of 8
Main reading, mean correlation every 50 ms from -600 to +400 ms:
-600    0.609
-550    0.572
-500    0.679
-450    0.748
-400    0.773
-350    0.762
-300    0.764
-250    0.713
-200    0.572
-150    0.523
-100    0.374
-50     0.278
 0      0.218
 50     0.161
 100    0.149
 150    0.110
 200    0.064
 250    0.010
 300    0.039
 350    0.033
 400   -0.001
 450   -0.022
 500   -0.049
 550   -0.032
 600   -0.050
 650   -0.005
 700   -0.000
 750   -0.018
 800   -0.035

inc = fig4[("a", "concat", "a", "a")][0]
ctx = ctxs[("a", "concat", "a")]
gt = np.arange(-200, 500)
gtrace = ps.condition_means(R["go_trace"], trials, ctx["mask"])[inc]
prep = ctx["prep"][inc]
pc = prep - prep.mean(1, keepdims=True)
gc = gtrace - gtrace.mean(1, keepdims=True)
go_curve = np.nanmean((pc[:, :, None] * gc).sum(1) / np.sqrt((pc ** 2).sum(1)[:, None] * (gc ** 2).sum(1)), 0)
print("Go-locked mean correlation every 50 ms from -200 to +450 ms after the go cue:")
print(pd.Series(go_curve[::50], index=gt[::50]).round(3).to_string())
fig, axs = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
axs[0].plot(gt, go_curve, "k"); axs[0].set_xlabel("ms from go cue")
axs[1].plot(tt, main_curve, "k"); axs[1].set_xlabel("ms from movement onset")
for a in axs:
    a.axhline(0, color="0.7"); a.set_ylim(-0.4, 1.05)
axs[0].set_ylabel("mean r with preparatory rates")
plt.show()
print("Curve values outside the axes:", int(np.sum((main_curve < -0.4) | (main_curve > 1.05)) + np.sum((go_curve < -0.4) | (go_curve > 1.05))))

Go-locked mean correlation every 50 ms from -200 to +450 ms after the go cue:
-200    0.634
-150    0.787
-100    0.784
-50     0.780
 0      0.799
 50     0.778
 100    0.676
 150    0.553
 200    0.447
 250    0.355
 300    0.289
 350    0.214
 400    0.160
 450    0.159

Curve values outside the axes: 0

rows = []
for key, (inc, r, rs) in fig4.items():
    ap = (units.id_array.values[inc] == 1).astype(float)
    if len(set(ap)) < 2:
        continue
    lr = stats.linregress(ap, r)
    rows.append({"proto": key[0], "combine": key[1], "onset": key[2], "tuning": key[3],
                 "n_anterior": int(ap.sum()), "n_posterior": int((1 - ap).sum()),
                 "mean_r_anterior": r[ap == 1].mean(), "mean_r_posterior": r[ap == 0].mean(),
                 "regression_r": lr.rvalue, "p": lr.pvalue, "no_effect_p_ge_05": lr.pvalue >= 0.05})
ap_tab = pd.DataFrame(rows)
print(ap_tab.round(3).to_string(index=False))
print("Readings with no effect (p >= 0.05):", int(ap_tab.no_effect_p_ge_05.sum()), "of", len(ap_tab))

proto combine onset tuning  n_anterior  n_posterior  mean_r_anterior  mean_r_posterior  regression_r     p  no_effect_p_ge_05
    a  concat     a      a          12            8            0.209             0.176         0.047 0.844               True
    a  concat     b      a          11            8            0.140             0.073         0.098 0.690               True
    a    each     a      a          13            9            0.173             0.078         0.138 0.541               True
    a    each     b      a          12            8            0.167             0.060         0.148 0.532               True
    b  concat     a      a          12           10            0.195             0.064         0.177 0.432               True
    b  concat     b      a          11            8            0.128             0.072         0.083 0.736               True
    b    each     a      a          13           11            0.178             0.093         0.117 0.585               True
    b    each     b      a          13           11            0.168             0.069         0.140 0.513               True
Readings with no effect (p >= 0.05): 8 of 8

def n_for_90(M):
    Mc = M - M.mean(0)
    s = np.linalg.svd(Mc, compute_uv=False) ** 2
    return int(np.searchsorted(np.cumsum(s) / s.sum(), 0.9) + 1), np.cumsum(s) / s.sum()

rows = []
for (proto, combine, onset), ctx in ctxs.items():
    M = ctx["prep"].T
    n_all, cum_all = n_for_90(M)
    row = {"proto": proto, "combine": combine, "onset": onset, "units_all": M.shape[1], "k90_all": n_all}
    for tun in ps.READINGS["tuning"]:
        inc = ps.tuning(ctx["prep"], tun) > ps.ARRAY_MIN_TUNING
        row[f"units_included_{tun}"] = int(inc.sum())
        row[f"k90_included_{tun}"] = n_for_90(M[:, inc])[0]
    rows.append(row)
dim = pd.DataFrame(rows)
print("Expected: 7 to 10, widened to 6.3 to 11")
print(dim.to_string(index=False))
_, cum = n_for_90(ctxs[("a", "concat", "a")]["prep"].T)
print("Main reading, all units, cumulative variance of the first 15 PCs:", np.round(cum[:15], 3))
# the same on a split of trials, to see how much of the 90% is trial noise
rng = np.random.default_rng(4)
ctx = ctxs[("a", "concat", "a")]
half = ctx["mask"] & (rng.random(len(trials)) < 0.5)
other = ctx["mask"] & ~half
A = ps.condition_means(R["prep"], trials, half).T
B = ps.condition_means(R["prep"], trials, other).T
print("Split halves of trials (seed 4): k90 =", n_for_90(A)[0], "and", n_for_90(B)[0])
print("Trial counts per condition in each half:", np.bincount(trials.cond_idx[half], minlength=27).min(), "to",
      np.bincount(trials.cond_idx[half], minlength=27).max())

Expected: 7 to 10, widened to 6.3 to 11
proto combine onset  units_all  k90_all  units_included_a  k90_included_a  units_included_b  k90_included_b
    a  concat     a        162       11                92              10                 6               4
    a  concat     b        162       11                92              10                 7               4
    a    each     a        162       12               101              11                 7               4
    a    each     b        162       12                95              11                 7               4
    b  concat     a        162       11                91              10                 6               4
    b  concat     b        162       11                92              10                 6               4
    b    each     a        162       12               103              11                 7               4
    b    each     b        162       12                97              11                 7               4
Main reading, all units, cumulative variance of the first 15 PCs: [0.297 0.49  0.601 0.682 0.738 0.79  0.829 0.856 0.877 0.896 0.91  0.924
 0.934 0.944 0.951]
Split halves of trials (seed 4): k90 = 14 and 13
Trial counts per condition in each half: 2 to 11

t0 = time.time()
peri_rows, kin_rows, task_rows = [], [], []
for (proto, combine, onset), ctx in ctxs.items():
    for tun in ps.READINGS["tuning"]:
        inc = ps.tuning(ctx["prep"], tun) > ps.ARRAY_MIN_TUNING
        X = ctx["prep"][inc]
        base = {"proto": proto, "combine": combine, "onset": onset, "tuning": tun, "n_units": int(inc.sum())}
        for t_units in ps.READINGS["t_units"]:
            for fold in ps.READINGS["fold"]:
                for start in [-150, 0]:
                    err = ps.peri_errors(ctx, inc, t_units, fold, 14, start)
                    for k in range(1, 15):
                        pf = ps.perf_both(err[:, k - 1], X)
                        for norm in ps.READINGS["norm"]:
                            peri_rows.append({**base, "t_units": t_units, "fold": fold, "start": start,
                                              "k": k, "norm": norm, "perf": pf[norm]})
        for fold in ps.READINGS["fold"]:
            for kin_end in ps.READINGS["kin_end"]:
                err = ps.kinematic_errors(ctx, inc, fold, kin_end)
                for k in range(1, 15):
                    pf = ps.perf_both(err[:, k - 1], X)
                    for norm in ps.READINGS["norm"]:
                        kin_rows.append({**base, "fold": fold, "kin_end": kin_end, "k": k, "norm": norm, "perf": pf[norm]})
        for task_end in ps.READINGS["task_end"]:
            for split in ps.READINGS["split"]:
                errs, dims = ps.task_errors(ctx, inc, task_end, split)
                for name, err in errs.items():
                    pf = ps.perf_both(err, X)
                    for norm in ps.READINGS["norm"]:
                        task_rows.append({**base, "task_end": task_end, "split": split, "space": name,
                                          "k": dims[name], "norm": norm, "perf": pf[norm]})
peri_df, kin_df, task_df = pd.DataFrame(peri_rows), pd.DataFrame(kin_rows), pd.DataFrame(task_rows)
print(f"Grid computed in {time.time() - t0:.0f} s: {len(peri_df)} peri-movement, {len(kin_df)} kinematic, "
      f"{len(task_df)} task-space values")

Grid computed in 120 s: 3584 peri-movement, 1792 kinematic, 1024 task-space values

PREP_KEYS = ["proto", "combine", "onset", "tuning"]

def combos():
    names = ["proto", "combine", "onset", "tuning", "t_units", "fold", "norm", "kin_end", "task_end", "split", "kmax"]
    for vals in itertools.product(*[ps.READINGS[n] for n in names]):
        yield dict(zip(names, vals))

pi = peri_df[peri_df.start == -150].set_index(PREP_KEYS + ["t_units", "fold", "norm", "k"]).perf.sort_index()
pi0 = peri_df[peri_df.start == 0].set_index(PREP_KEYS + ["t_units", "fold", "norm", "k"]).perf.sort_index()
ki = kin_df.set_index(PREP_KEYS + ["fold", "kin_end", "norm", "k"]).perf.sort_index()
ti = task_df.set_index(PREP_KEYS + ["task_end", "split", "norm", "space"]).sort_index()

def curves_for(c):
    pk = tuple(c[k] for k in PREP_KEYS)
    ks = np.arange(1, c["kmax"] + 1)
    peri = pi.loc[pk + (c["t_units"], c["fold"], c["norm"])].loc[ks]
    peri0 = pi0.loc[pk + (c["t_units"], c["fold"], c["norm"])].loc[ks]
    kin = ki.loc[pk + (c["fold"], c["kin_end"], c["norm"])].loc[ks]
    task = ti.loc[pk + (c["task_end"], c["split"], c["norm"])]
    return peri, peri0, kin, task

rows = []
for c in combos():
    peri, peri0, kin, task = curves_for(c)
    best_other = max(kin.max(), task.perf.max())
    rows.append({**c, "peri_best": peri.max(), "peri_best_k": int(peri.idxmax()), "peri_k1": peri.loc[1],
                 "peri_kmax": peri.loc[c["kmax"]],
                 "kin_best": kin.max(), "kin_best_k": int(kin.idxmax()), "kin_k1": kin.loc[1], "kin_kmax": kin.loc[c["kmax"]],
                 "task_best": task.perf.max(), "task_best_space": task.perf.idxmax(),
                 "n_task_above_0": int((task.perf > 0).sum()),
                 "peri0_best": peri0.max(), "peri0_best_k": int(peri0.idxmax()),
                 "peri_beats_all": peri.max() > best_other, "margin": peri.max() - best_other})
grid = pd.DataFrame(rows)
MAIN = ps.MAIN
is_main = (grid[list(MAIN)] == pd.Series(MAIN)).all(1)
print("Combinations of readings:", len(grid), "; peri-movement space best in", int(grid.peri_beats_all.sum()))
print("\nMain reading:")
print(grid[is_main].T.to_string())
print("\nShare of combinations where the peri-movement space is best, by each reading:")
for n in ["proto", "combine", "onset", "tuning", "t_units", "fold", "norm", "kin_end", "task_end", "split", "kmax"]:
    print(f"  {n}: " + ", ".join(f"{v} {grid[grid[n] == v].peri_beats_all.mean():.2f}" for v in ps.READINGS[n]))
print("\nMargin (peri best minus best other) over combinations: min", round(grid.margin.min(), 3),
      "median", round(grid.margin.median(), 3), "max", round(grid.margin.max(), 3))
print("Peri-movement best k over combinations:", grid.peri_best_k.value_counts().sort_index().to_dict())

Combinations of readings: 2048 ; peri-movement space best in 1468

Main reading:
                            0
proto                       a
combine                concat
onset                       a
tuning                      a
t_units                   all
fold                     once
norm                   neuron
kin_end                 fixed
task_end                speed
split                  length
kmax                       14
peri_best            0.350618
peri_best_k                10
peri_k1              0.056274
peri_kmax           -0.455551
kin_best             0.103079
kin_best_k                  4
kin_k1               0.061186
kin_kmax            -0.451454
task_best            0.250235
task_best_space  kitchen sink
n_task_above_0              7
peri0_best           0.318372
peri0_best_k                9
peri_beats_all           True
margin               0.100384

Share of combinations where the peri-movement space is best, by each reading:
  proto: a 0.76, b 0.67
  combine: concat 0.76, each 0.68
  onset: a 0.67, b 0.76
  tuning: a 0.95, b 0.49
  t_units: all 0.94, included 0.50
  fold: once 0.66, refit 0.77
  norm: neuron 0.70, pooled 0.73
  kin_end: fixed 0.72, speed 0.72
  task_end: speed 0.70, fixed 0.73
  split: length 0.70, time 0.73
  kmax: 14 0.72, 12 0.71

Margin (peri best minus best other) over combinations: min -0.543 median 0.129 max 0.338
Peri-movement best k over combinations: {2: 16, 3: 32, 4: 32, 5: 80, 6: 112, 7: 304, 8: 112, 9: 144, 10: 160, 11: 128, 12: 472, 13: 72, 14: 384}

/tmp/ipykernel_2482/1839405699.py:35: Pandas4Warning: Starting with pandas version 4.0 all arguments of all will be keyword-only.
  is_main = (grid[list(MAIN)] == pd.Series(MAIN)).all(1)

c = dict(MAIN)
peri, peri0, kin, task = curves_for(c)
pk = tuple(c[k] for k in PREP_KEYS)
n_inc = int(peri_df.set_index(PREP_KEYS).loc[pk].n_units.iloc[0])
print(f"Main reading, {n_inc} included units. Mean performance by k:")
print(pd.DataFrame({"peri-movement": peri, "kinematic": kin}).round(3).T.to_string())
print("\nTask spaces (dimensionality, mean performance):")
print(task[["k", "perf"]].round(4).to_string())
alt = dict(MAIN, fold="refit")
p2, _, k2, t2 = curves_for(alt)
print("\nVariant fold = refit:")
print(pd.DataFrame({"peri-movement": p2, "kinematic": k2}).round(3).T.to_string())
fig, axs = plt.subplots(1, 2, figsize=(9, 3.2), sharey=True)
for ax, (pp, kk, tk, title) in zip(axs, [(peri, kin, task, "space computed once"), (p2, k2, t2, "space refit per fold")]):
    ax.plot(pp.index, pp.values, "k-o", ms=3, label="peri-movement")
    ax.plot(kk.index, kk.values, "b-o", ms=3, label="kinematic")
    ax.plot(tk.k, tk.perf, "gs", label="task spaces")
    ax.axhline(0, color="0.7"); ax.set_ylim(-0.5, 0.5); ax.set_xlabel("dimensionality"); ax.set_title(title)
axs[0].set_ylabel("mean performance"); axs[0].legend(fontsize=7)
plt.show()
allv = np.r_[peri.values, kin.values, task.perf.values, p2.values, k2.values, t2.perf.values]
print("Values below -0.5 (outside the axes):", int((allv < -0.5).sum()), "; min", round(allv.min(), 3))

Main reading, 92 included units. Mean performance by k:
k                 1      2      3      4      5      6      7      8      9      10     11     12     13     14
peri-movement  0.056  0.124  0.126  0.163  0.238  0.277  0.336  0.350  0.348  0.351  0.322  0.251 -0.047 -0.456
kinematic      0.061  0.061  0.097  0.103  0.079  0.046  0.073  0.095  0.051  0.043 -0.128 -0.235 -0.312 -0.451

Task spaces (dimensionality, mean performance):
                                   k    perf
space                                       
endpoints + halfpoint              4  0.1021
endpoints + max speed              3  0.0003
endpoints + max velocity           4 -0.1246
endpoints + orthogonal excursion   3  0.0419
endpoints + viapoint               4  0.1267
initial velocity                   2  0.0927
kitchen sink                      11  0.2502
segment directions                 6  0.1164

Variant fold = refit:
k                 1      2      3      4      5      6      7      8      9      10     11     12     13     14
peri-movement  0.054  0.128  0.132  0.140  0.153  0.205  0.355  0.382  0.403  0.419  0.447  0.454  0.458  0.461
kinematic      0.062  0.061  0.098  0.103  0.102  0.100  0.103  0.129  0.131  0.160  0.128  0.092  0.108  0.085

Values below -0.5 (outside the axes): 0 ; min -0.456

rows = []
for c in combos():
    peri, peri0, kin, task = curves_for(c)
    wins = [peri.loc[k] > kin.loc[k] for k in kin.index]
    wins += [peri.loc[int(r.k)] > r.perf for _, r in task.iterrows() if r.k <= c["kmax"]]
    e3 = task.loc["endpoints + max speed"].perf
    rows.append({**c, "comparisons": len(wins), "won": int(sum(wins)), "share": np.mean(wins),
                 "peri_k3": peri.loc[3], "endpoints_max_speed": e3, "beats_e3": peri.loc[3] > e3})
each = pd.DataFrame(rows)
each["criterion_met"] = (each.share >= 0.8) & each.beats_e3
print("Main reading:", each[is_main.values][["comparisons", "won", "share", "peri_k3", "endpoints_max_speed", "criterion_met"]].round(3).to_dict("records"))
print("Combinations meeting the criterion:", int(each.criterion_met.sum()), "of", len(each))
print("Share won over combinations: min", round(each.share.min(), 3), "median", round(each.share.median(), 3), "max", round(each.share.max(), 3))
print("By reading, fraction of combinations meeting the criterion:")
for n in ["proto", "combine", "onset", "tuning", "t_units", "fold", "norm", "kin_end", "task_end", "split", "kmax"]:
    print(f"  {n}: " + ", ".join(f"{v} {each[each[n] == v].criterion_met.mean():.2f}" for v in ps.READINGS[n]))
peri, _, kin, task = curves_for(dict(MAIN))
print("\nMain reading, kinematic minus peri-movement by k:", (kin - peri).round(3).to_dict())
# the main reading with the space refit without the left-out condition
cr = dict(MAIN, fold="refit")
peri_r, _, kin_r, task_r = curves_for(cr)
print("\nMain reading with fold=refit, mean performance by k:")
print(pd.DataFrame({"peri-movement": peri_r, "kinematic": kin_r}).round(3).T.to_string())
print("Task spaces (refit does not apply):", task_r.perf.round(3).to_dict())
print("Main reading with fold=refit:", each[(each[list(cr)] == pd.Series(cr)).all(axis=1)][
    ["comparisons", "won", "share", "peri_k3", "endpoints_max_speed", "criterion_met"]].round(3).to_dict("records"))
print("Comparisons lost under fold=refit, main reading otherwise:",
      [f"kinematic k={k}" for k in kin_r.index if not peri_r.loc[k] > kin_r.loc[k]]
      + [f"{s} k={int(r.k)}" for s, r in task_r.iterrows() if r.k <= cr["kmax"] and not peri_r.loc[int(r.k)] > r.perf])

Main reading: [{'comparisons': 22, 'won': 20, 'share': 0.909, 'peri_k3': 0.126, 'endpoints_max_speed': 0.0, 'criterion_met': True}]
Combinations meeting the criterion: 1074 of 2048
Share won over combinations: min 0.0 median 0.85 max 1.0
By reading, fraction of combinations meeting the criterion:
  proto: a 0.53, b 0.52
  combine: concat 0.69, each 0.36
  onset: a 0.52, b 0.53
  tuning: a 0.75, b 0.30
  t_units: all 0.67, included 0.38
  fold: once 0.37, refit 0.68
  norm: neuron 0.49, pooled 0.56
  kin_end: fixed 0.54, speed 0.51
  task_end: speed 0.52, fixed 0.53
  split: length 0.52, time 0.53
  kmax: 14 0.52, 12 0.53

Main reading, kinematic minus peri-movement by k: {1: 0.005, 2: -0.063, 3: -0.029, 4: -0.06, 5: -0.16, 6: -0.231, 7: -0.264, 8: -0.255, 9: -0.296, 10: -0.308, 11: -0.45, 12: -0.486, 13: -0.265, 14: 0.004}

Main reading with fold=refit, mean performance by k:
k                 1      2      3      4      5      6      7      8      9      10     11     12     13     14
peri-movement  0.054  0.128  0.132  0.140  0.153  0.205  0.355  0.382  0.403  0.419  0.447  0.454  0.458  0.461
kinematic      0.062  0.061  0.098  0.103  0.102  0.100  0.103  0.129  0.131  0.160  0.128  0.092  0.108  0.085
Task spaces (refit does not apply): {'endpoints + halfpoint': 0.102, 'endpoints + max speed': 0.0, 'endpoints + max velocity': -0.125, 'endpoints + orthogonal excursion': 0.042, 'endpoints + viapoint': 0.127, 'initial velocity': 0.093, 'kitchen sink': 0.25, 'segment directions': 0.116}
Main reading with fold=refit: [{'comparisons': 22, 'won': 21, 'share': 0.955, 'peri_k3': 0.132, 'endpoints_max_speed': 0.0, 'criterion_met': True}]
Comparisons lost under fold=refit, main reading otherwise: ['kinematic k=1']

trend = grid.copy()
trend["peri_rises"] = trend.peri_best > trend.peri_k1
trend["peri_best_k_le_10"] = trend.peri_best_k <= 10
trend["peri_falls"] = trend.peri_kmax < trend.peri_best
trend["criterion_met"] = trend.peri_rises & trend.peri_best_k_le_10 & trend.peri_falls
trend["kin_met"] = (trend.kin_best > trend.kin_k1) & (trend.kin_best_k <= 10) & (trend.kin_kmax < trend.kin_best)
print("Main reading:", trend[is_main][["peri_k1", "peri_best", "peri_best_k", "peri_kmax", "criterion_met",
                                        "kin_k1", "kin_best", "kin_best_k", "kin_kmax", "kin_met"]].round(3).to_dict("records"))
print("Combinations meeting the criterion (peri-movement):", int(trend.criterion_met.sum()), "of", len(trend),
      "; kinematic:", int(trend.kin_met.sum()))
print("By reading, fraction meeting the criterion:")
for n in ["proto", "combine", "onset", "tuning", "t_units", "fold", "norm", "kmax"]:
    print(f"  {n}: " + ", ".join(f"{v} {trend[trend[n] == v].criterion_met.mean():.2f}" for v in ps.READINGS[n]))
print("\nBest k of the peri-movement space, by fold reading:")
print(trend.groupby("fold").peri_best_k.describe().to_string())

Main reading: [{'peri_k1': 0.056, 'peri_best': 0.351, 'peri_best_k': 10, 'peri_kmax': -0.456, 'criterion_met': True, 'kin_k1': 0.061, 'kin_best': 0.103, 'kin_best_k': 4, 'kin_kmax': -0.451, 'kin_met': True}]
Combinations meeting the criterion (peri-movement): 992 of 2048 ; kinematic: 1776
By reading, fraction meeting the criterion:
  proto: a 0.47, b 0.50
  combine: concat 0.46, each 0.51
  onset: a 0.48, b 0.49
  tuning: a 0.50, b 0.47
  t_units: all 0.42, included 0.55
  fold: once 0.91, refit 0.06
  norm: neuron 0.52, pooled 0.45
  kmax: 14 0.47, 12 0.50

Best k of the peri-movement space, by fold reading:
        count       mean       std   min   25%   50%   75%   max
fold                                                            
once   1024.0   7.484375  2.092614   2.0   6.0   7.0   9.0  11.0
refit  1024.0  12.664062  1.188260  10.0  12.0  12.0  14.0  14.0

ad = grid.copy()
ad["criterion_met"] = (ad.peri_best > 0) & (ad.kin_best > 0) & (ad.n_task_above_0 > 4)
print("Main reading:", ad[is_main][["peri_best", "kin_best", "n_task_above_0", "criterion_met"]].round(3).to_dict("records"))
print("Combinations meeting the criterion:", int(ad.criterion_met.sum()), "of", len(ad))
print("Task spaces above zero over combinations:", ad.n_task_above_0.value_counts().sort_index().to_dict(), "(of 8)")
print("By reading, fraction meeting the criterion:")
for n in ["proto", "combine", "onset", "tuning", "fold", "norm", "kin_end", "task_end", "split"]:
    print(f"  {n}: " + ", ".join(f"{v} {ad[ad[n] == v].criterion_met.mean():.2f}" for v in ps.READINGS[n]))
print("\nPer task space, share of combinations with mean performance above zero:")
print(task_df.groupby("space").perf.apply(lambda x: (x > 0).mean()).round(2).to_string())
print("\nPer task space, main reading: performance per unit above zero (count of units):")
ctx = ctxs[("a", "concat", "a")]
inc = ps.tuning(ctx["prep"], "a") > 5
errs, dims = ps.task_errors(ctx, inc, "speed", "length")
print({k: f"{int((ps.performance(e, ctx['prep'][inc]) > 0).sum())} of {int(inc.sum())}" for k, e in errs.items()})

Main reading: [{'peri_best': 0.351, 'kin_best': 0.103, 'n_task_above_0': 7, 'criterion_met': True}]
Combinations meeting the criterion: 1936 of 2048
Task spaces above zero over combinations: {3: 80, 4: 32, 5: 112, 6: 192, 7: 256, 8: 1376} (of 8)
By reading, fraction meeting the criterion:
  proto: a 0.95, b 0.94
  combine: concat 0.95, each 0.94
  onset: a 0.95, b 0.94
  tuning: a 0.89, b 1.00
  fold: once 0.95, refit 0.95
  norm: neuron 0.89, pooled 1.00
  kin_end: fixed 0.95, speed 0.95
  task_end: speed 1.00, fixed 0.89
  split: length 1.00, time 0.89

Per task space, share of combinations with mean performance above zero:
space
endpoints + halfpoint               0.94
endpoints + max speed               0.81
endpoints + max velocity            0.67
endpoints + orthogonal excursion    1.00
endpoints + viapoint                1.00
initial velocity                    1.00
kitchen sink                        0.91
segment directions                  0.94

Per task space, main reading: performance per unit above zero (count of units):
{'initial velocity': '56 of 92', 'endpoints + max speed': '42 of 92', 'endpoints + orthogonal excursion': '43 of 92', 'endpoints + viapoint': '60 of 92', 'endpoints + max velocity': '30 of 92', 'endpoints + halfpoint': '58 of 92', 'segment directions': '62 of 92', 'kitchen sink': '69 of 92'}

def shuffled_move(move, rng):
    i0 = ps.T_WINDOW[0] - ps.MOVE_TRACE[0]
    flip = rng.random(move.shape[:2]) < 0.5
    out = move.copy()
    ref = move[:, :, i0:i0 + 1]
    out[:, :, i0:] = np.where(flip[:, :, None], 2 * ref - move[:, :, i0:], move[:, :, i0:])
    return out

variants = [dict(MAIN)] + [dict(MAIN, **{n: ps.READINGS[n][1]}) for n in ["combine", "proto", "onset", "t_units", "fold", "norm", "kin_end"]]
rng7a = np.random.default_rng(7)
rows = []
t0 = time.time()
for v in variants:
    ctx = ctxs[(v["proto"], v["combine"], v["onset"])]
    inc = ps.tuning(ctx["prep"], v["tuning"]) > ps.ARRAY_MIN_TUNING
    X = ctx["prep"][inc]
    n_draws = 20
    draws = np.zeros((n_draws, 14))
    for d in range(n_draws):
        err = ps.peri_errors(ctx, inc, v["t_units"], v["fold"], 14, move=shuffled_move(ctx["move"], rng7a))
        draws[d] = [ps.perf_both(err[:, k], X)[v["norm"]] for k in range(14)]
    shuf = draws.mean(0)
    peri, _, kin, task = curves_for(v)
    rows.append({"variant": ", ".join(f"{n}={v[n]}" for n in v if v[n] != MAIN[n]) or "main", "draws": n_draws,
                 "shuffled_best": shuf.max(), "shuffled_best_k": int(shuf.argmax() + 1),
                 "shuffled_best_se": draws[:, shuf.argmax()].std(ddof=1) / np.sqrt(n_draws),
                 "peri_best": peri.max(), "kin_best": kin.max(),
                 "below_peri_every_k": bool(np.all(shuf < peri.values[:14])),
                 "below_kin_best": shuf.max() < kin.max(),
                 "shuffled_by_k": np.round(shuf, 3).tolist(), "peri_by_k": np.round(peri.values[:14], 3).tolist()})
fig7a = pd.DataFrame(rows)
fig7a["criterion_met"] = fig7a.below_peri_every_k & fig7a.below_kin_best
print(f"computed in {time.time() - t0:.0f} s")
print(fig7a.drop(columns=["shuffled_by_k", "peri_by_k"]).round(3).to_string(index=False))
print("\nShuffled-space performance by k:")
for _, r in fig7a.iterrows():
    print(f"  {r.variant}: {r.shuffled_by_k}")
print("\nPeri-movement performance by k, for the variants where the shuffled space is not below it at every k:")
for _, r in fig7a[~fig7a.below_peri_every_k].iterrows():
    print(f"  {r.variant}: {r.peri_by_k}")

computed in 65 s
         variant  draws  shuffled_best  shuffled_best_k  shuffled_best_se  peri_best  kin_best  below_peri_every_k  below_kin_best  criterion_met
            main     20          0.014                3             0.008      0.351     0.103                True            True           True
    combine=each     20         -0.014                2             0.012      0.150     0.068               False            True          False
         proto=b     20          0.023                3             0.009      0.338     0.098                True            True           True
         onset=b     20          0.029                4             0.011      0.351     0.099                True            True           True
t_units=included     20          0.017                3             0.010      0.353     0.103                True            True           True
      fold=refit     20          0.159               14             0.009      0.461     0.160                True            True           True
     norm=pooled     20          0.100                3             0.016      0.562     0.289                True            True           True
   kin_end=speed     20          0.027                3             0.008      0.351     0.109                True            True           True

Shuffled-space performance by k:
  main: [-0.043, 0.009, 0.014, -0.0, -0.036, -0.085, -0.144, -0.179, -0.283, -0.378, -0.48, -0.605, -0.829, -1.061]
  combine=each: [-0.055, -0.014, -0.037, -0.07, -0.112, -0.164, -0.238, -0.3, -0.391, -0.482, -0.61, -0.758, -0.952, -1.242]
  proto=b: [-0.037, 0.009, 0.023, -0.0, -0.026, -0.057, -0.084, -0.156, -0.233, -0.333, -0.44, -0.582, -0.753, -0.941]
  onset=b: [-0.04, 0.008, 0.022, 0.029, -0.003, -0.028, -0.081, -0.143, -0.197, -0.286, -0.396, -0.502, -0.634, -0.827]
  t_units=included: [-0.04, 0.013, 0.017, -0.013, -0.037, -0.072, -0.123, -0.18, -0.254, -0.342, -0.438, -0.55, -0.698, -0.85]
  fold=refit: [-0.043, 0.025, 0.064, 0.083, 0.098, 0.117, 0.129, 0.135, 0.142, 0.149, 0.15, 0.155, 0.156, 0.159]
  norm=pooled: [-0.018, 0.071, 0.1, 0.089, 0.065, 0.057, 0.045, 0.019, -0.028, -0.092, -0.164, -0.265, -0.375, -0.53]
  kin_end=speed: [-0.042, 0.006, 0.027, 0.017, -0.001, -0.036, -0.086, -0.137, -0.225, -0.33, -0.422, -0.535, -0.681, -0.861]

Peri-movement performance by k, for the variants where the shuffled space is not below it at every k:
  combine=each: [0.026, 0.071, 0.08, 0.05, 0.094, 0.15, 0.132, -0.183, -0.425, -1.666, -1.756, -1.71, -1.838, -3.102]

fb = grid.copy()
fb["ratio"] = fb.peri0_best / fb.peri_best
fb["within_tol"] = (fb.ratio - 0.97).abs() <= 0.097
print("Expected 0.97 +/- 0.097 (0.873 to 1.067)")
print("Main reading:", fb[is_main][["peri_best", "peri_best_k", "peri0_best", "peri0_best_k", "ratio", "within_tol"]].round(3).to_dict("records"))
print("Combinations within tolerance:", int(fb.within_tol.sum()), "of", len(fb))
print("Ratio over combinations: min", round(fb.ratio.min(), 3), "median", round(fb.ratio.median(), 3), "max", round(fb.ratio.max(), 3))
print("By reading, median ratio and fraction within tolerance:")
for n in ["proto", "combine", "onset", "tuning", "t_units", "fold", "norm", "kmax"]:
    print(f"  {n}: " + ", ".join(f"{v} {fb[fb[n] == v].ratio.median():.3f}/{fb[fb[n] == v].within_tol.mean():.2f}" for v in ps.READINGS[n]))
# finer sweep of the start under the main reading
ctx = ctxs[("a", "concat", "a")]
inc = ps.tuning(ctx["prep"], "a") > 5
X = ctx["prep"][inc]
sweep = {}
for start in [-150, -100, -50, 0, 50, 100]:
    err = ps.peri_errors(ctx, inc, "all", "once", 14, start)
    pf = [ps.perf_both(err[:, k], X)["neuron"] for k in range(14)]
    sweep[start] = (round(max(pf), 3), int(np.argmax(pf) + 1))
print("\nMain reading, best performance (and k) by window start (ms from onset):", sweep)

Expected 0.97 +/- 0.097 (0.873 to 1.067)
Main reading: [{'peri_best': 0.351, 'peri_best_k': 10, 'peri0_best': 0.318, 'peri0_best_k': 9, 'ratio': 0.908, 'within_tol': True}]
Combinations within tolerance: 1336 of 2048
Ratio over combinations: min 0.367 median 0.917 max 2.004
By reading, median ratio and fraction within tolerance:
  proto: a 0.915/0.69, b 0.919/0.62
  combine: concat 0.930/0.82, each 0.886/0.48
  onset: a 0.914/0.62, b 0.920/0.68
  tuning: a 0.912/0.70, b 0.932/0.60
  t_units: all 0.934/0.75, included 0.899/0.55
  fold: once 0.872/0.47, refit 0.930/0.84
  norm: neuron 0.905/0.63, pooled 0.929/0.67
  kmax: 14 0.920/0.66, 12 0.912/0.65

Main reading, best performance (and k) by window start (ms from onset): {-150: (0.351, 10), -100: (0.34, 8), -50: (0.33, 9), 0: (0.318, 9), 50: (0.303, 9), 100: (0.279, 9)}

def subset_scores(ctx, members, inc, fold, t_units, self_in_upper):
    """Mean performance at k = 6 in the peri-movement space of `members`, and of the upper-limit
    space (k = 1..14), over included units of the subset."""
    fit = members & inc
    in_t = members if t_units == "all" else fit
    X = ctx["prep"][fit]
    Gp = ps.neuron_grams(ps.t_matrix(ctx["move"]))
    G = Gp[in_t].sum(0)
    err = ps.loo_pca_space(G, Gp[fit] * in_t[fit][:, None, None], X, 6, fold)
    peri6 = np.mean(ps.performance(err[:, 5], X))
    P = ctx["prep"][:, :, None]
    Pp = ps.neuron_grams(P)
    Gu = Pp[members].sum(0)
    parts = None if self_in_upper else Pp[fit]
    if self_in_upper:
        erru = ps.loo_pca_space(Gu, None, X, 14, fold)
    else:
        erru = ps.loo_pca_space(Gu, parts, X, 14, fold)
    upper = np.array([np.mean(ps.performance(erru[:, k], X)) for k in range(14)])
    return peri6, upper, int(fit.sum())

rng7c = np.random.default_rng(17)
rows = []
for v in [dict(MAIN), dict(MAIN, fold="refit"), dict(MAIN, combine="each"), dict(MAIN, t_units="included")]:
    ctx = ctxs[(v["proto"], v["combine"], v["onset"])]
    inc = ps.tuning(ctx["prep"], v["tuning"]) > ps.ARRAY_MIN_TUNING
    for self_in in [False, True]:
        res = {}
        for name, members in [("anterior", units.id_array.values == 1), ("posterior", units.id_array.values == 2),
                              ("all", np.ones(len(units), bool))]:
            res[name] = subset_scores(ctx, members, inc, v["fold"], v["t_units"], self_in)
        rand = []
        for d in range(20):
            perm = rng7c.permutation(len(units))
            for half in (perm[: len(units) // 2], perm[len(units) // 2:]):
                m = np.zeros(len(units), bool); m[half] = True
                rand.append(subset_scores(ctx, m, inc, v["fold"], v["t_units"], self_in))
        for up in ["k6", "best"]:
            def norm(p, u):
                return p / (u[5] if up == "k6" else u.max())
            r_all = np.array([norm(p, u) for p, u, _ in rand])
            r_norm = r_all.mean()
            a_norm, p_norm = norm(*res["anterior"][:2]), norm(*res["posterior"][:2])
            rows.append({"variant": ", ".join(f"{n}={v[n]}" for n in v if v[n] != MAIN[n]) or "main",
                         "self_in_upper": self_in, "upper": up,
                         "n_anterior": res["anterior"][2], "n_posterior": res["posterior"][2],
                         "all_norm": norm(*res["all"][:2]), "anterior_norm": a_norm, "posterior_norm": p_norm,
                         "random_half_norm": r_norm, "random_half_sd": r_all.std(ddof=1),
                         "random_half_above_1": np.mean(r_all > 1),
                         "random_half_p5": np.percentile(r_all, 5), "random_half_p95": np.percentile(r_all, 95),
                         "anterior_raw": res["anterior"][0], "posterior_raw": res["posterior"][0],
                         "random_raw": np.mean([p for p, _, _ in rand]),
                         "ant_vs_random": a_norm / r_norm, "post_vs_random": p_norm / r_norm})
fig7c = pd.DataFrame(rows)
fig7c["criterion_met"] = ((fig7c.anterior_norm > 0) & (fig7c.posterior_norm > 0)
                          & ((fig7c.ant_vs_random - 1).abs() <= 0.1) & ((fig7c.post_vs_random - 1).abs() <= 0.1))
fig7c["anterior_better"] = fig7c.anterior_norm > fig7c.posterior_norm
print(fig7c.round(3).to_string(index=False))
print("Rows meeting the criterion:", int(fig7c.criterion_met.sum()), "of", len(fig7c),
      "; with the unit being fit in its upper-limit space (main):",
      int(fig7c[fig7c.self_in_upper].criterion_met.sum()), "of", int(fig7c.self_in_upper.sum()),
      "; with it left out:", int(fig7c[~fig7c.self_in_upper].criterion_met.sum()), "of", int((~fig7c.self_in_upper).sum()))
print("Main reading:", fig7c[(fig7c.variant == "main") & fig7c.self_in_upper & (fig7c.upper == "k6")]
      .round(3).to_dict("records"))
print("Anterior and posterior normalized performance against the 5th to 95th percentile of the 40 random halves:")
print(fig7c[["variant", "self_in_upper", "upper", "anterior_norm", "posterior_norm", "random_half_p5", "random_half_p95"]]
      .assign(anterior_inside=lambda d: d.anterior_norm.between(d.random_half_p5, d.random_half_p95),
              posterior_inside=lambda d: d.posterior_norm.between(d.random_half_p5, d.random_half_p95))
      .round(3).to_string(index=False))

         variant  self_in_upper upper  n_anterior  n_posterior  all_norm  anterior_norm  posterior_norm  random_half_norm  random_half_sd  random_half_above_1  random_half_p5  random_half_p95  anterior_raw  posterior_raw  random_raw  ant_vs_random  post_vs_random  criterion_met  anterior_better
            main          False    k6          55           37     0.878          0.785           0.814             0.961           0.190                0.275           0.746            1.252         0.291          0.143       0.262          0.817           0.847          False            False
            main          False  best          55           37     0.755          0.576           0.814             0.855           0.123                0.125           0.684            1.091         0.291          0.143       0.262          0.674           0.952          False            False
            main           True    k6          55           37     0.722          0.597           0.382             0.620           0.079                0.000           0.512            0.751         0.291          0.143       0.258          0.963           0.615          False             True
            main           True  best          55           37     0.617          0.479           0.322             0.495           0.067                0.000           0.390            0.591         0.291          0.143       0.258          0.967           0.650          False             True
      fold=refit          False    k6          55           37     0.687          0.708           0.994             0.781           0.210                0.150           0.545            1.172         0.227          0.196       0.211          0.907           1.273          False            False
      fold=refit          False  best          55           37     0.463          0.418           0.772             0.547           0.105                0.000           0.375            0.686         0.227          0.196       0.211          0.764           1.410          False            False
      fold=refit           True    k6          55           37     0.567          0.527           0.539             0.543           0.125                0.000           0.363            0.740         0.227          0.196       0.211          0.970           0.993           True            False
      fold=refit           True  best          55           37     0.378          0.331           0.383             0.357           0.072                0.000           0.252            0.475         0.227          0.196       0.211          0.928           1.075           True            False
    combine=each          False    k6          56           45     0.708          0.665           0.763            -0.826          10.811                0.275           0.212            1.543         0.216         -0.014       0.119         -0.805          -0.924          False            False
    combine=each          False  best          56           45     0.691          0.508          -0.526             0.720           0.259                0.075           0.315            1.156         0.216         -0.014       0.119          0.706          -0.731          False             True
    combine=each           True    k6          56           45     0.534          0.485          -0.072             0.374           0.135                0.000           0.166            0.541         0.216         -0.014       0.114          1.295          -0.192          False             True
    combine=each           True  best          56           45     0.496          0.395          -0.067             0.334           0.126                0.000           0.137            0.511         0.216         -0.014       0.114          1.183          -0.200          False             True
t_units=included          False    k6          55           37     0.856          0.808           0.802             0.955           0.239                0.350           0.678            1.506         0.300          0.141       0.246          0.847           0.840          False             True
t_units=included          False  best          55           37     0.736          0.593           0.802             0.800           0.126                0.100           0.642            1.041         0.300          0.141       0.246          0.742           1.002          False            False
t_units=included           True    k6          55           37     0.704          0.615           0.376             0.599           0.087                0.000           0.450            0.702         0.300          0.141       0.251          1.027           0.628          False             True
t_units=included           True  best          55           37     0.601          0.493           0.317             0.483           0.059                0.000           0.404            0.549         0.300          0.141       0.251          1.020           0.656          False             True
Rows meeting the criterion: 2 of 16 ; with the unit being fit in its upper-limit space (main): 2 of 8 ; with it left out: 0 of 8
Main reading: [{'variant': 'main', 'self_in_upper': True, 'upper': 'k6', 'n_anterior': 55, 'n_posterior': 37, 'all_norm': 0.722, 'anterior_norm': 0.597, 'posterior_norm': 0.382, 'random_half_norm': 0.62, 'random_half_sd': 0.079, 'random_half_above_1': 0.0, 'random_half_p5': 0.512, 'random_half_p95': 0.751, 'anterior_raw': 0.291, 'posterior_raw': 0.143, 'random_raw': 0.258, 'ant_vs_random': 0.963, 'post_vs_random': 0.615, 'criterion_met': False, 'anterior_better': True}]
Anterior and posterior normalized performance against the 5th to 95th percentile of the 40 random halves:
         variant  self_in_upper upper  anterior_norm  posterior_norm  random_half_p5  random_half_p95  anterior_inside  posterior_inside
            main          False    k6          0.785           0.814           0.746            1.252             True              True
            main          False  best          0.576           0.814           0.684            1.091            False              True
            main           True    k6          0.597           0.382           0.512            0.751             True             False
            main           True  best          0.479           0.322           0.390            0.591             True             False
      fold=refit          False    k6          0.708           0.994           0.545            1.172             True              True
      fold=refit          False  best          0.418           0.772           0.375            0.686             True             False
      fold=refit           True    k6          0.527           0.539           0.363            0.740             True              True
      fold=refit           True  best          0.331           0.383           0.252            0.475             True              True
    combine=each          False    k6          0.665           0.763           0.212            1.543             True              True
    combine=each          False  best          0.508          -0.526           0.315            1.156             True             False
    combine=each           True    k6          0.485          -0.072           0.166            0.541             True             False
    combine=each           True  best          0.395          -0.067           0.137            0.511             True             False
t_units=included          False    k6          0.808           0.802           0.678            1.506             True              True
t_units=included          False  best          0.593           0.802           0.642            1.041            False              True
t_units=included           True    k6          0.615           0.376           0.450            0.702             True             False
t_units=included           True  best          0.493           0.317           0.404            0.549             True             False

def consistency(S, X, rng, n_draws=100):
    c = S.shape[0]
    D = np.column_stack([np.ones(c), S])
    vals = np.zeros(X.shape[0])
    for _ in range(n_draws):
        i1, i2 = rng.integers(0, c, 25), rng.integers(0, c, 25)
        b1 = np.linalg.lstsq(D[i1], X[:, i1].T, rcond=None)[0][1:]
        b2 = np.linalg.lstsq(D[i2], X[:, i2].T, rcond=None)[0][1:]
        b1 /= np.linalg.norm(b1, axis=0, keepdims=True)
        b2 /= np.linalg.norm(b2, axis=0, keepdims=True)
        vals += (b1 * b2).sum(0) / n_draws
    return vals

rng7d = np.random.default_rng(23)
rows = []
for v in [dict(MAIN), dict(MAIN, combine="each"), dict(MAIN, t_units="included"), dict(MAIN, task_end="fixed"),
          dict(MAIN, split="time"), dict(MAIN, kin_end="speed"), dict(MAIN, onset="b"), dict(MAIN, proto="b")]:
    ctx = ctxs[(v["proto"], v["combine"], v["onset"])]
    inc = ps.tuning(ctx["prep"], v["tuning"]) > ps.ARRAY_MIN_TUNING
    X = ctx["prep"][inc]
    spaces, _ = ps.task_spaces(ctx["pos"], ctx["vel"], ctx["times"], trials, v["task_end"], v["split"])
    _, _, kin, _ = curves_for(v)
    kbest = int(kin.idxmax())
    KT = ps.kinematic_matrix(ctx["pos"], ctx["vel"], ctx["times"], ps.kinematic_end_ms(ctx, v["kin_end"]))
    Skin = ps._scores_from_gram(ps.centre_gram(KT @ KT.T), 14)[0][:, :kbest]
    Gp = ps.neuron_grams(ps.t_matrix(ctx["move"]))
    in_t = np.ones(len(inc), bool) if v["t_units"] == "all" else inc
    G = Gp[in_t].sum(0)
    for scale in ["raw", "unit SD"]:
        def sc(S):
            return S / S.std(0) if scale == "unit SD" else S
        comps = [("endpoints + max speed", sc(spaces["endpoints + max speed"]), 3),
                 ("segment directions", sc(spaces["segment directions"]), 6),
                 (f"kinematic k={kbest}", sc(Skin), kbest)]
        for name, S, k in comps:
            other = consistency(S, X, rng7d)
            peri_c = np.zeros(X.shape[0])
            for j, u in enumerate(np.nonzero(inc)[0]):
                Gj = G - Gp[u] * in_t[u]
                Sp = sc(ps._scores_from_gram(ps.centre_gram(Gj), k)[0])
                peri_c[j] = consistency(Sp, X[j:j + 1], rng7d)[0]
            rows.append({"variant": ", ".join(f"{n}={v[n]}" for n in v if v[n] != MAIN[n]) or "main",
                         "scale": scale, "space": name, "k": k, "other_consistency": other.mean(),
                         "peri_consistency": peri_c.mean(), "peri_higher": peri_c.mean() > other.mean(),
                         "units_peri_higher": int((peri_c > other).sum()), "n_units": len(other)})
fig7d = pd.DataFrame(rows)
# how much two draws of 25 of the 27 conditions with replacement overlap
ov = np.random.default_rng(24)
d1 = [set(ov.integers(0, 27, 25)) for _ in range(10000)]
d2 = [set(ov.integers(0, 27, 25)) for _ in range(10000)]
print(f"Distinct conditions per draw: mean {np.mean([len(a) for a in d1]):.1f} of 27; "
      f"shared by two draws: mean {np.mean([len(a & b) for a, b in zip(d1, d2)]):.1f}")
print(fig7d.round(3).to_string(index=False))
by = fig7d.groupby(["variant", "scale"]).peri_higher.all()
print("\nPeri-movement higher in every comparison, by variant and scaling:")
print(by.to_string())

Distinct conditions per draw: mean 16.5 of 27; shared by two draws: mean 10.1
         variant   scale                 space  k  other_consistency  peri_consistency  peri_higher  units_peri_higher  n_units
            main     raw endpoints + max speed  3              0.610             0.756         True                 74       92
            main     raw    segment directions  6              0.570             0.684         True                 64       92
            main     raw         kinematic k=4  4              0.666             0.748         True                 69       92
            main unit SD endpoints + max speed  3              0.641             0.788         True                 81       92
            main unit SD    segment directions  6              0.562             0.722         True                 71       92
            main unit SD         kinematic k=4  4              0.755             0.781         True                 56       92
    combine=each     raw endpoints + max speed  3              0.619             0.732         True                 78      101
    combine=each     raw    segment directions  6              0.554             0.668         True                 68      101
    combine=each     raw         kinematic k=4  4              0.642             0.681         True                 65      101
    combine=each unit SD endpoints + max speed  3              0.623             0.750         True                 83      101
    combine=each unit SD    segment directions  6              0.555             0.700         True                 67      101
    combine=each unit SD         kinematic k=4  4              0.713             0.726         True                 58      101
t_units=included     raw endpoints + max speed  3              0.607             0.750         True                 77       92
t_units=included     raw    segment directions  6              0.593             0.684         True                 59       92
t_units=included     raw         kinematic k=4  4              0.621             0.746         True                 71       92
t_units=included unit SD endpoints + max speed  3              0.622             0.786         True                 80       92
t_units=included unit SD    segment directions  6              0.575             0.728         True                 66       92
t_units=included unit SD         kinematic k=4  4              0.742             0.783         True                 61       92
  task_end=fixed     raw endpoints + max speed  3              0.629             0.752         True                 68       92
  task_end=fixed     raw    segment directions  6              0.588             0.681         True                 59       92
  task_end=fixed     raw         kinematic k=4  4              0.665             0.744         True                 66       92
  task_end=fixed unit SD endpoints + max speed  3              0.650             0.789         True                 82       92
  task_end=fixed unit SD    segment directions  6              0.557             0.724         True                 69       92
  task_end=fixed unit SD         kinematic k=4  4              0.758             0.778         True                 52       92
      split=time     raw endpoints + max speed  3              0.609             0.752         True                 78       92
      split=time     raw    segment directions  6              0.610             0.682         True                 58       92
      split=time     raw         kinematic k=4  4              0.655             0.741         True                 69       92
      split=time unit SD endpoints + max speed  3              0.622             0.786         True                 81       92
      split=time unit SD    segment directions  6              0.619             0.726         True                 65       92
      split=time unit SD         kinematic k=4  4              0.746             0.784         True                 56       92
   kin_end=speed     raw endpoints + max speed  3              0.614             0.757         True                 80       92
   kin_end=speed     raw    segment directions  6              0.567             0.683         True                 60       92
   kin_end=speed     raw         kinematic k=4  4              0.675             0.743         True                 69       92
   kin_end=speed unit SD endpoints + max speed  3              0.624             0.786         True                 82       92
   kin_end=speed unit SD    segment directions  6              0.606             0.723         True                 66       92
   kin_end=speed unit SD         kinematic k=4  4              0.744             0.779         True                 55       92
         onset=b     raw endpoints + max speed  3              0.622             0.758         True                 73       92
         onset=b     raw    segment directions  6              0.567             0.684         True                 64       92
         onset=b     raw         kinematic k=4  4              0.646             0.744         True                 69       92
         onset=b unit SD endpoints + max speed  3              0.615             0.789         True                 83       92
         onset=b unit SD    segment directions  6              0.579             0.723         True                 65       92
         onset=b unit SD         kinematic k=4  4              0.741             0.777         True                 57       92
         proto=b     raw endpoints + max speed  3              0.594             0.755         True                 79       91
         proto=b     raw    segment directions  6              0.559             0.676         True                 63       91
         proto=b     raw         kinematic k=4  4              0.645             0.739         True                 65       91
         proto=b unit SD endpoints + max speed  3              0.630             0.786         True                 78       91
         proto=b unit SD    segment directions  6              0.569             0.721         True                 67       91
         proto=b unit SD         kinematic k=4  4              0.740             0.781         True                 59       91

Peri-movement higher in every comparison, by variant and scaling:
variant           scale  
combine=each      raw        True
                  unit SD    True
kin_end=speed     raw        True
                  unit SD    True
main              raw        True
                  unit SD    True
onset=b           raw        True
                  unit SD    True
proto=b           raw        True
                  unit SD    True
split=time        raw        True
                  unit SD    True
t_units=included  raw        True
                  unit SD    True
task_end=fixed    raw        True
                  unit SD    True

rng8 = np.random.default_rng(8)
rows = []
for init in ["uniform", "gaussian"]:
    for release in [360, 180]:
        rs, perfs = [], []
        for d in range(1000):
            act, x0 = ps.oscillator(rng8, 25, 50, release, init)
            rel = act[:, 50:, :]
            rs.append(np.corrcoef(x0[:, 0], rel[:, :, 0].mean(1))[0, 1])
            T = rel.transpose(0, 2, 1).reshape(25, -1)
            S = ps._scores_from_gram(ps.centre_gram(T @ T.T), 2)[0]
            err = ps.loo_errors_fixed(S, x0[None, :, 0])
            perfs.append(ps.performance(err, x0[None, :, 0])[0])
        rs, perfs = np.array(rs), np.array(perfs)
        lo, hi = np.percentile(rs, [2.5, 97.5])
        rows.append({"init": init, "interval_steps": release, "median_r": np.median(rs), "r_2.5%": lo, "r_97.5%": hi,
                     "0.08_inside": lo <= 0.08 <= hi, "min_pd_perf": perfs.min(), "median_pd_perf": np.median(perfs),
                     "pd_ge_0.99": perfs.min() >= 0.99})
fig8 = pd.DataFrame(rows)
print(fig8.round(4).to_string(index=False))

    init  interval_steps  median_r  r_2.5%  r_97.5%  0.08_inside  min_pd_perf  median_pd_perf  pd_ge_0.99
 uniform             360    0.0119 -0.3968   0.3870         True          1.0             1.0        True
 uniform             180   -0.0193 -0.4114   0.3799         True          1.0             1.0        True
gaussian             360    0.0004 -0.4000   0.4037         True          1.0             1.0        True
gaussian             180   -0.0065 -0.4087   0.3835         True          1.0             1.0        True

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

Notebook run time: 7.2 min

Claim	Tested	Outcome
fig4abc-tuning-correlation-single-electrode	no	blocked, data missing
fig4d-tuning-correlation-j-array	yes	partially reproduced (assumption)
fig4e-correlation-time-course	yes	not reproduced
results-ap-location-no-effect	J-array only	reproduced (assumption)
supp-preparatory-dimensionality	J-array only	reproduced (assumption)
fig6-peri-movement-space-best	J-array only	reproduced
fig6-each-dimensionality	J-array only	reproduced (assumption)
fig6-dimensionality-trend	J-array only	reproduced (assumption)
fig6-pd-adequate-all-spaces	J-array only	reproduced (assumption)
fig6e-summary-rank	no	blocked, data missing
fig7a-shuffled-space-poor	J-array only	reproduced (assumption)
fig7b-window-start-onset	J-array only	reproduced (assumption)
fig7c-anterior-posterior	J-array only	not reproduced (assumption)
fig7d-pd-consistency	J-array only	reproduced
fig8-oscillator-model	simulation	reproduced

Claim	Outcome	Main numbers
fig4abc-tuning-correlation-single-electrode	blocked (data missing)	single-electrode data not in the dandiset
fig4d-tuning-correlation-j-array	partially reproduced (assumption)	mean r 0.195, Spearman 0.212, p = 0.023
fig4e-correlation-time-course	not reproduced	r 0.218 at movement onset
results-ap-location-no-effect	reproduced (assumption)	r = 0.047, p = 0.844 (J-array)
supp-preparatory-dimensionality	reproduced (assumption)	11 components
fig6-peri-movement-space-best	reproduced	0.351 against 0.250 and 0.103
fig6-each-dimensionality	reproduced (assumption)	20 of 22 comparisons
fig6-dimensionality-trend	reproduced (assumption)	best k 10, falls to -0.456 at 14
fig6-pd-adequate-all-spaces	reproduced (assumption)	7 of 8 task spaces above 0
fig6e-summary-rank	blocked (data missing)	monkeys A, B, J and EMG not in the dandiset
fig7a-shuffled-space-poor	reproduced (assumption)	shuffled best 0.014
fig7b-window-start-onset	reproduced (assumption)	ratio 0.908
fig7c-anterior-posterior	not reproduced (assumption)	anterior 0.963, posterior 0.615 of random halves
fig7d-pd-consistency	reproduced	higher in 3 of 3 comparisons
fig8-oscillator-model	reproduced	0.08 inside the central 95%

Reproducing Churchland et al. (2010) from dandiset 000138¶

Data check¶

Readings¶

fig4abc-tuning-correlation-single-electrode¶

fig4d-tuning-correlation-j-array¶

fig4e-correlation-time-course¶

results-ap-location-no-effect¶

supp-preparatory-dimensionality¶

PD performance in each space (Figure 6, J-array)¶

fig6-peri-movement-space-best¶

fig6-each-dimensionality¶

fig6-dimensionality-trend¶

fig6-pd-adequate-all-spaces¶

fig6e-summary-rank¶

fig7a-shuffled-space-poor¶

fig7b-window-start-onset¶

fig7c-anterior-posterior¶

fig7d-pd-consistency¶

fig8-oscillator-model¶

Summary¶