import time
import io
import warnings
import json
import urllib.request

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import io as sio
from scipy import signal, stats

import chamber

warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
started = time.time()
f = chamber.open_file()
events = chamber.load_events(f)
trials = chamber.trial_table(f, events)
fscv = chamber.load_fscv_table(f)
lfp_t, lfp_v = chamber.load_lfp(f)
thr = chamber.spike_units(f, "processing/ecephys/thresholded_units")
plx = chamber.spike_units(f, "units")
print(f"dandiset {chamber.DANDISET_ID} version {chamber.DANDISET_VERSION}, {chamber.ASSET_PATH}; "
      f"loaded in {time.time() - started:.0f} s")

dandiset 002014 version 0.261012.0921, sub-Monkey-T/sub-Monkey-T_ses-09272024_ecephys.nwb; loaded in 10 s

def text(path):
    x = f[path][()]
    return x.decode() if isinstance(x, bytes) else x

print("Subject:", {k: text(f"general/subject/{k}") for k in ("subject_id", "species", "sex", "age")})
print("Subject description:", text("general/subject/description"))
print("Session start:", text("session_start_time"), "; session id:", text("general/session_id"))
print("Experiment description:", text("general/experiment_description")[:160], "...")
print("Keywords:", [k.decode() for k in f["general/keywords"][()]])
print("Related publication:", [k.decode() for k in f["general/related_publications"][()]])

el = f["general/extracellular_ephys/electrodes"]
cols = ["location", "channel_name", "group_name", "recording_opened", "recording_closed", "DspLowCutFrequency", "DspHighCutFrequency", "InputRange"]
print("\nElectrodes table:")
print(pd.DataFrame({c: [x.decode() if isinstance(x, bytes) else x for x in el[c][()]] for c in cols}).to_string())

series = {
    "raw ephys": "acquisition/electrical_series",
    "bipolar LFP": "processing/ecephys/FilteredEphys/bipolar_lfp_series",
    "raw FSCV current": "acquisition/fscv_response_series",
    "FSCV applied voltage": "stimulus/presentation/fscv_excitation_series",
}
rows = []
for name, path in series.items():
    g = f[path]
    if "timestamps" in g:
        ts = g["timestamps"]
        span = (ts[0], ts[-1])
        step = np.median(np.diff(ts[:5000]))
    else:
        rate = g["starting_time"].attrs["rate"]
        span = (g["starting_time"][()], g["starting_time"][()] + g["data"].shape[0] / rate)
        step = 1 / rate
    rows.append({"series": name, "shape": g["data"].shape, "dtype": str(g["data"].dtype),
                 "unit": g["data"].attrs["unit"], "conversion": g["data"].attrs["conversion"],
                 "electrodes": f[path + "/electrodes"][()].tolist() if "electrodes" in f[path] else None, "first_s": round(span[0], 3),
                 "last_s": round(span[1], 3), "sample_step_s": step})
print("\nTime series (times in s from session start):")
print(pd.DataFrame(rows).to_string(index=False))

Subject: {'subject_id': 'Monkey T', 'species': 'Macaca mulatta', 'sex': 'F', 'age': 'P11Y'}
Subject description: Two adult female rhesus monkeys were used for in vivo validation of the s-µIP device, with chronic implants targeting the striatum for simultaneous FSCV and EPhys recordings during behavioral tasks.
Session start: 2024-09-27T08:49:52-04:00 ; session id: 09272024
Experiment description: This experiment contains in vivo recordings of reward- and movement-related dopamine, spikes, and local field potentials obtained using single- and dual-channel ...
Keywords: ['fast scan cyclic voltammetry', 'micro-invasive neural implant', 'dopamine', 'striatum']
Related publication: ['https://doi.org/10.1016/j.jneumeth.2025.110467']

Electrodes table:
  location channel_name          group_name     recording_opened     recording_closed  DspLowCutFrequency DspHighCutFrequency InputRange
0      c7b        CSC47      ElectrodeGroup  2024-09-27 08:49:52  2024-09-27 15:12:24                 0.1                8000       1000
1      c5c         CSC7      ElectrodeGroup  2024-09-27 08:49:52  2024-09-27 15:12:22                 0.1                8000       1000
2      c5d        CSC76      ElectrodeGroup  2024-09-27 08:49:53  2024-09-27 15:12:27                 0.1                8000       1000
3     c8ds            6  FSCVElectrodeGroup                                                           NaN

Time series (times in s from session start):
              series          shape   dtype    unit   conversion electrodes   first_s    last_s  sample_step_s
           raw ephys (586751679, 3)   int16   volts 3.051758e-08  [0, 1, 2]     0.000 18335.990       0.000031
         bipolar LFP  (18144393, 1) float64   volts 3.051758e-08        [0] 11824.993 29968.682       0.001000
    raw FSCV current    (42115200,) float64 amperes 1.000000e+00        [3] 11854.993 31534.902       0.000040
FSCV applied voltage    (42115200,) float64   volts 1.000000e+00       None 11854.993 31534.902       0.000040

counts = chamber.raw_event_counts(f)
print("Events as stored and after removing repeats:")
print(counts.to_string(index=False))
print("\nTrials:", len(trials), "; trial windows", np.unique(np.round(trials.stop_time - trials.start_time, 6)), "s;",
      "trial-start events matching the trials' midpoint_time:",
      int(np.sum([np.min(np.abs(events["intended trial start"] - m)) < 1e-6 for m in trials.trial_start])))
tags = f["intervals/trials/tags"][()]
print("Trial tags:", pd.Series([t.decode() for t in tags]).value_counts().to_dict())
print("Trials with more than one V:", int((trials.n_V > 1).sum()))
print("\nTrial outcomes (first error event, or the reward size):")
print(pd.crosstab(trials.error.fillna("none"), trials.reward.fillna("no reward")).to_string())
succ = trials[trials.success]
print("\nSuccessful trials: large", int((succ.reward == "large").sum()), ", small", int((succ.reward == "small").sum()))
print("Successful trials by side and reward:")
print(pd.crosstab(succ.side, succ.reward).to_string())
lat = pd.DataFrame({"V - trial start": succ.V - succ.trial_start, "RW - V": succ.RW - succ.V})
print("\nLatencies in successful trials (s):")
print(lat.describe(percentiles=[0.05, 0.5, 0.95]).round(3).to_string())

Events as stored and after removing repeats:
                                   label                     description  stored  unique  sorted
                    intended trial start The event times for code '128'.    7221    1412   False
                          left condition The event times for code '117'.    3763     733   False
            central cue fixation started  The event times for code '50'.    6478    1285   False
                      value object start  The event times for code '23'.    6076    1249   False
              value cue fixation started  The event times for code '51'.    5785    1204   False
                                feedback  The event times for code '21'.    4209     952   False
                              big reward The event times for code '101'.    1951     453   False
                     reward delivery end The event times for code '105'.    4206     952   False
                               end trial  The event times for code '18'.    6960    1412   False
                         right condition The event times for code '118'.    3458     679   False
                            small reward The event times for code '100'.    2249     499   False
  error value cue initial fixation break  The event times for code '65'.     237      45   False
          error value cue fixation break  The event times for code '63'.    1300     252   False
        error central cue fixation break  The event times for code '60'.     197      36   False
            error fixation never started  The event times for code '62'.     692     122   False
error central cue initial fixation break  The event times for code '61'.      30       5    True

Trials: 1411 ; trial windows [60.] s; trial-start events matching the trials' midpoint_time: 1411
Trial tags: {'initial cue display': 1411}
Trials with more than one V: 0

Trial outcomes (first error event, or the reward size):
reward                                    large  no reward  small
error                                                            
error central cue fixation break              0         36      0
error central cue initial fixation break      0          5      0
error fixation never started                  0        122      0
error value cue fixation break                0        252      0
error value cue initial fixation break        0         45      0
none                                        452          0    499

Successful trials: large 452 , small 499
Successful trials by side and reward:
reward  large  small
side                
left      224    258
right     228    241

Latencies in successful trials (s):
       V - trial start   RW - V
count          951.000  951.000
mean             1.916    4.780
std              0.299    0.066
min              1.549    4.716
5%               1.707    4.733
50%              1.816    4.766
95%              2.583    4.849
max              3.549    5.499

# The paper's trial lists on Zenodo (doi:10.5281/zenodo.14825992), 1-based trial ids.
ZENODO = "https://zenodo.org/api/records/14825992/files/"
lists = {}
for size, name in (("large", "big"), ("small", "small")):
    url = ZENODO + urllib.request.quote(f"09272024_nlx_c8c_value object start_{name} rew and success_tr1-1411.mat") + "/content"
    with urllib.request.urlopen(url) as r:
        lists[size] = set(sio.loadmat(io.BytesIO(r.read()))["trids_out"].ravel().tolist())
for size in ("large", "small"):
    mine = set(succ.trial[succ.reward == size])
    print(f"{size}: Zenodo list {len(lists[size])} trials, dandiset events {len(mine)}, both {len(mine & lists[size])}, "
          f"only in the dandiset {sorted(mine - lists[size])}, only in the list {sorted(lists[size] - mine)}")
print("Last trial:", trials.iloc[-1][["trial", "reward", "error", "trial_start", "stop_time"]].to_dict(),
      "; last LFP sample", round(lfp_t[-1], 3))

large: Zenodo list 451 trials, dandiset events 452, both 451, only in the dandiset [1411], only in the list []
small: Zenodo list 499 trials, dandiset events 499, both 499, only in the dandiset [], only in the list []
Last trial: {'trial': 1411, 'reward': 'large', 'error': nan, 'trial_start': 29938.681817, 'stop_time': 29968.681817} ; last LFP sample 29968.682

# Units and ranges of every variable used.
print("Trial-aligned FSCV table: rows", fscv["da"].shape[0], ", samples per row", fscv["da"].shape[1],
      ", rate", np.unique(fscv["rate"]), "Hz, good", int(fscv["good"].sum()), "of", len(fscv["good"]))
print("Rows start at the trials' start_time:", bool(np.allclose(fscv["start_time"], trials.start_time)),
      "; row time axis", fscv["t_rel_start"][[0, -1]], "s from the trial start")
for k in ("da", "ph", "m", "iox"):
    x = fscv[k]
    print(f"  {k:3s}: unit not stated; range {np.nanmin(x):.3f} to {np.nanmax(x):.3f}; NaN {np.isnan(x).mean():.3%} of samples, "
          f"in {np.isnan(x).any(1).sum()} rows")
print("Columns holding NaN in more than 5% of rows (s from trial start):",
      np.round(fscv["t_rel_start"][np.isnan(fscv["da"]).mean(0) > 0.05][[0, -1]], 1))
print("Dopamine NaN in successful trials within V-1 s to V+6 s:",
      f"{np.isnan(chamber.dopamine_aligned(fscv, succ, np.arange(-1, 6.01, 0.1))[1]).mean():.2%} of samples")

ok = ~np.isnan(lfp_v)
print(f"\nBipolar LFP: values {np.nanmin(lfp_v):.2e} to {np.nanmax(lfp_v):.2e} as stored, SD {np.nanstd(lfp_v):.2e};"
      f" NaN {1 - ok.mean():.2%}; sample step {np.median(np.diff(lfp_t[:100000])):.6f} s; gaps over 1.5 ms: "
      f"{int(np.sum(np.diff(lfp_t) > 0.0015))}")
print("The stored conversion factor would scale them to", f"{np.nanstd(lfp_v) * f['processing/ecephys/FilteredEphys/bipolar_lfp_series/data'].attrs['conversion']:.1e}", "V SD")
seg_raw = chamber.raw_segment(f, 3000, 300, slice(0, 3))
raw_bip = signal.decimate(signal.decimate(seg_raw[:, 0] - seg_raw[:, 1], 8), 4)
k0 = np.searchsorted(lfp_t, 20000)
fr, p_lfp = signal.welch(np.nan_to_num(lfp_v[k0:k0 + 300000]), fs=1000, nperseg=4096)
_, p_raw = signal.welch(raw_bip, fs=1000, nperseg=4096)
bands = [(1, 4), (4, 8), (8, 13), (13, 30), (30, 58), (58, 62), (62, 100), (100, 300)]
print("\nPower spectral density (V^2/Hz): stored LFP as stored, and raw c7b minus c5c in volts (300 s each)")
print(pd.DataFrame([{"band_Hz": f"{a}-{b}", "stored_LFP": p_lfp[(fr >= a) & (fr < b)].mean(),
                     "raw_bipolar": p_raw[(fr >= a) & (fr < b)].mean()} for a, b in bands]).to_string(index=False))
print("\nRaw ephys in volts (300 s from 3000 s): range per channel", np.round(seg_raw.min(0) * 1e6), "to", np.round(seg_raw.max(0) * 1e6), "uV")
full_scale = f["acquisition/electrical_series/data"].attrs["conversion"] * 32767
clip = []
for start in np.linspace(100, 18000, 6):
    x = chamber.raw_segment(f, start, 30, slice(0, 3))
    clip.append(np.mean(np.abs(x) >= 0.999 * full_scale, axis=0))
print("Fraction of samples at the input range (+-1 mV), 30 s at six points, per channel (c7b, c5c, c5d):", np.round(np.array(clip), 4).tolist())
peak10 = p_raw[(fr >= 9.5) & (fr <= 10.5)].max() / np.median(p_raw[(fr >= 4) & (fr <= 8)])
print(f"Raw bipolar: peak density at 9.5-10.5 Hz is {peak10:.0f} times the median at 4-8 Hz; stored LFP: "
      f"{p_lfp[(fr >= 9.5) & (fr <= 10.5)].max() / np.median(p_lfp[(fr >= 4) & (fr <= 8)]):.1f} times")

Trial-aligned FSCV table: rows 1411 , samples per row 601 , rate [10.] Hz, good 1411 of 1411
Rows start at the trials' start_time: True ; row time axis [-30.  30.] s from the trial start
  da : unit not stated; range -53.560 to 42.437; NaN 5.019% of samples, in 1053 rows
  ph : unit not stated; range -2.372 to 2.372; NaN 0.000% of samples, in 0 rows
  m  : unit not stated; range -0.474 to 0.484; NaN 0.000% of samples, in 0 rows
  iox: unit not stated; range -1.040 to 0.805; NaN 0.000% of samples, in 0 rows
Columns holding NaN in more than 5% of rows (s from trial start): [-8.9 30. ]
Dopamine NaN in successful trials within V-1 s to V+6 s: 9.71% of samples

Bipolar LFP: values -6.38e-04 to 9.69e-04 as stored, SD 4.25e-05; NaN 1.60%; sample step 0.001000 s; gaps over 1.5 ms: 0
The stored conversion factor would scale them to 1.3e-12 V SD

Power spectral density (V^2/Hz): stored LFP as stored, and raw c7b minus c5c in volts (300 s each)
band_Hz   stored_LFP  raw_bipolar
    1-4 4.208983e-11 3.271436e-11
    4-8 5.819154e-12 5.004912e-12
   8-13 2.106545e-12 3.111888e-10
  13-30 2.996336e-12 1.233671e-10
  30-58 8.276712e-13 1.536662e-10
  58-62 4.163779e-13 3.405444e-10
 62-100 4.138574e-13 8.882068e-11
100-300 1.110503e-13 1.023431e-11

Raw ephys in volts (300 s from 3000 s): range per channel [-954. -624. -253.] to [1000. 1000.  141.] uV

Fraction of samples at the input range (+-1 mV), 30 s at six points, per channel (c7b, c5c, c5d): [[0.0002, 0.0161, 0.0], [0.0, 0.0167, 0.0], [0.0, 0.0168, 0.0], [0.0, 0.0173, 0.0], [0.0, 0.0176, 0.0], [0.0, 0.0175, 0.0]]
Raw bipolar: peak density at 9.5-10.5 Hz is 911 times the median at 4-8 Hz; stored LFP: 0.1 times

# Dopamine is referenced at V. Where the sample at V is NaN, `reference_at` interpolates between the
# nearest valid samples on each side. Check on a linear ramp with a NaN gap at every position of V
# relative to the 10 Hz grid: the referenced trace must be 0 at V.
worst = 0.0
for phase in np.linspace(0, 0.1, 101):
    t_syn = np.arange(-30, 30.01, 0.1) + phase
    for t0 in (phase, phase + 0.03, phase + 0.05, 0.0):
        ramp = 2 * t_syn
        ramp[np.abs(t_syn - t0) < 0.06] = np.nan
        ref = chamber.reference_at(ramp, t_syn, t0)
        worst = max(worst, np.nanmax(np.abs(ref - 2 * (t_syn - t0))))
print(f"Synthetic ramp: largest error of the reference {worst:.2e}")
assert worst < 1e-9
on_nan = sum(np.isnan(np.interp(row.V, row.trial_start + fscv["t_rel_start"], fscv["da"][row.trial - 1])) for row in succ.itertuples())
print("Successful trials whose dopamine is NaN at V (reference interpolated across the gap, or trial dropped):", int(on_nan), "of", len(succ))

Synthetic ramp: largest error of the reference 3.55e-15
Successful trials whose dopamine is NaN at V (reference interpolated across the gap, or trial dropped): 62 of 951

# Spikes: the thresholded table (the paper's Zenodo file) and the Plexon unit.
g = f["processing/ecephys/thresholded_units"]
print("Thresholded table description:", g.attrs["description"])
print("Electrodes of its units:", g["electrodes"][()].tolist())
for name, st in thr.items():
    print(f"  thresholded unit {name}: {len(st)} spikes, {st.min():.3f} to {st.max():.3f} s, sorted {bool(np.all(np.diff(st) >= 0))}")
for name, st in plx.items():
    print(f"  Plexon unit {name}: {len(st)} spikes, {st.min():.3f} to {st.max():.3f} s, "
          f"sorted {bool(np.all(np.diff(st) >= 0))}, decreasing steps {int(np.sum(np.diff(st) < 0))}")
main_unit = max(thr, key=lambda k: len(thr[k]))
spk = thr[main_unit]
plx_spk = np.sort(next(iter(plx.values())))
nearest = np.minimum(np.abs(spk[np.clip(np.searchsorted(spk, plx_spk), 0, len(spk) - 1)] - plx_spk),
                     np.abs(spk[np.clip(np.searchsorted(spk, plx_spk) - 1, 0, len(spk) - 1)] - plx_spk))
print(f"Plexon spikes within 0.2 ms of a spike of thresholded unit {main_unit}: {np.sum(nearest < 2e-4)} of {len(plx_spk)}; "
      f"within 1 us: {np.sum(nearest < 1e-6)}")
wf = chamber.thresholded_waveforms(f, list(thr).index(main_unit), max_spikes=20000)
print(f"Waveforms: {wf.shape[1]} samples; unit attribute '{g['waveforms'].attrs['unit']}'; mean waveform {wf.mean(0).min():.1f} to "
      f"{wf.mean(0).max():.1f} as stored; fraction of spikes whose largest |value| exceeds 40: {np.mean(np.abs(wf).max(1) > 40):.3f}")
print("Datasets: waveforms", g["waveforms"].shape, ", waveforms_index", g["waveforms_index"].shape, ", spike_times",
      g["spike_times"].shape, ", spike_times_index", g["spike_times_index"].shape,
      "; waveform samples per spike (index steps)", np.unique(np.diff(np.r_[0, g["waveforms_index"][()]])).tolist())
es_dur = f["acquisition/electrical_series/data"].shape[0] / chamber.RAW_RATE
print(f"Raw ephys duration {es_dur:.2f} s; spikes {spk.min():.3f} to {spk.max():.3f} s")

Thresholded table description: EPhys captured and thresholded spike data from site c5d. Each row contains the initial time stamp for the waveform, the unit ID (-1 indicates a poor-quality spike), followed by the detected waveform samples (48) points long sampled at 32000.0 Hz).
Electrodes of its units: [2, 2]
  thresholded unit 0: 3 spikes, 9589.894 to 13183.096 s, sorted True
  thresholded unit 1: 101733 spikes, 0.011 to 18335.951 s, sorted True
  Plexon unit 0: 90912 spikes, 0.011 to 18335.951 s, sorted False, decreasing steps 3
Plexon spikes within 0.2 ms of a spike of thresholded unit 1: 54411 of 90912; within 1 us: 52625

Waveforms: 48 samples; unit attribute 'volts'; mean waveform -52.9 to 37.4 as stored; fraction of spikes whose largest |value| exceeds 40: 1.000
Datasets: waveforms (4883328, 1) , waveforms_index (101736,) , spike_times (101736,) , spike_times_index (2,) ; waveform samples per spike (index steps) [48]
Raw ephys duration 18335.99 s; spikes 0.011 to 18335.951 s

# Clock 1: where the raw ephys falls on the trial clock. Segments of raw c7b minus c5c are located in the
# stored LFP, which is placed on the trial clock by the trials' start times.
candidate, clock = chamber.locate_raw_in_lfp(f, lfp_t, lfp_v)
print(f"Raw c7b minus c5c at {candidate['raw_time']:.0f}-{candidate['raw_time'] + chamber.CLOCK_SEGMENT_S:.0f} s, searched in the whole stored LFP: "
      f"best at offset {candidate['offset']:.4f} s, r = {candidate['r']:.4f}")
print("Stored LFP segments searched in the raw signal within 60 s of that offset:")
print(clock.round(4).to_string(index=False))
OFFSET = float(np.median(clock.offset))
print(f"\nOffset of the trial clock from the raw ephys clock: median {OFFSET:.4f} s; range {clock.offset.min():.4f} to "
      f"{clock.offset.max():.4f} s; correlation {clock.r.min():.3f} to {clock.r.max():.3f}, best 1 s or more away at most {clock.best_r_1s_away.max():.3f}")
drift = np.polyfit(clock.raw_time, clock.offset, 1)[0]
print(f"Drift {drift * 3600 * 1000:.2f} ms per hour")

# Clock 2: the spike times against raw c5d (300-6000 Hz), spike-triggered, over 300 s.
lag_ms, snip = chamber.spike_triggered_raw(f, spk, 6000.0, 300.0, channel=2)
sta_c5d = snip.mean(0)
k_min = int(np.argmin(sta_c5d))
per_spike = lag_ms[np.argmin(snip, 1)]
print(f"\nSpike-triggered average of raw c5d over {len(snip)} spikes: trough {sta_c5d[k_min] * 1e6:.1f} uV at {lag_ms[k_min]:.2f} ms "
      f"from the stored spike time; peak-to-peak {np.ptp(sta_c5d) * 1e6:.1f} uV; mean |signal| far from the spikes "
      f"{np.median(np.abs(snip)) * 1e6:.1f} uV")
print("Per-spike trough lag (ms), 10th/50th/90th percentile:", np.round(np.percentile(per_spike, [10, 50, 90]), 2))
print(f"Per-spike trough lag within 0.5 ms of {lag_ms[k_min]:.2f} ms: {np.mean(np.abs(per_spike - lag_ms[k_min]) <= 0.5):.1%} of {len(per_spike)} spikes")
rng_clock = np.random.default_rng(0)
null_lag, null_snip = chamber.spike_triggered_raw(f, spk + rng_clock.uniform(0.05, 0.5, len(spk)), 6000.0, 300.0, channel=2)
print(f"Spike times moved by a random 50-500 ms: trough {null_snip.mean(0).min() * 1e6:.1f} uV")
for ch, site in ((0, "c7b"), (1, "c5c")):
    _, sn = chamber.spike_triggered_raw(f, spk, 6000.0, 300.0, channel=ch)
    print(f"Same at raw {site}: trough {sn.mean(0).min() * 1e6:.1f} uV, peak-to-peak {np.ptp(sn.mean(0)) * 1e6:.1f} uV")
SPIKE_LAG = lag_ms[k_min] / 1000
print("Mean stored waveform trough at sample", int(np.argmin(wf.mean(0))), f"({np.argmin(wf.mean(0)) / 32:.2f} ms after the waveform start)")

# Absolute times.
opened = pd.Timestamp(text("session_start_time"))
closed_hdr = el["recording_closed"][()][0].decode()
raw_end_trial_clock = OFFSET + es_dur
print(f"\nSession start {opened}; Neuralynx header: recording closed {closed_hdr}")
print(f"The raw ephys on the trial clock runs from {opened + pd.Timedelta(seconds=OFFSET)} to {opened + pd.Timedelta(seconds=raw_end_trial_clock)}")
print(f"Header open-to-close span {(pd.Timestamp(closed_hdr) - pd.Timestamp(opened.tz_localize(None))).total_seconds():.0f} s; raw samples {es_dur:.0f} s")

Raw c7b minus c5c at 60-90 s, searched in the whole stored LFP: best at offset 11812.1030 s, r = 0.9973
Stored LFP segments searched in the raw signal within 60 s of that offset:
 trial_clock_time   raw_time    offset      r  best_r_1s_away
       12184.9931   372.8901 11812.103 0.9887          0.3275
       14631.2349  2819.1289 11812.106 0.9871          0.3509
       17137.4762  5325.3662 11812.110 0.9767          0.2829
       19643.7176  7831.6056 11812.112 0.9978          0.2475
       22239.9585 10427.8444 11812.114 0.9910          0.2571
       24656.1997 12844.0816 11812.118 0.9954          0.2756
       27162.4413 15350.3203 11812.121 0.9976          0.2408
       29668.6828 17856.5588 11812.124 0.9956          0.2524

Offset of the trial clock from the raw ephys clock: median 11812.1130 s; range 11812.1030 to 11812.1240 s; correlation 0.977 to 0.998, best 1 s or more away at most 0.351
Drift 4.25 ms per hour

Spike-triggered average of raw c5d over 1633 spikes: trough -26.1 uV at -7.12 ms from the stored spike time; peak-to-peak 40.0 uV; mean |signal| far from the spikes 2.6 uV
Per-spike trough lag (ms), 10th/50th/90th percentile: [-7.41 -7.19 -6.88]
Per-spike trough lag within 0.5 ms of -7.12 ms: 98.8% of 1633 spikes

Spike times moved by a random 50-500 ms: trough -0.3 uV

Same at raw c7b: trough -2.4 uV, peak-to-peak 4.8 uV

Same at raw c5c: trough -1.4 uV, peak-to-peak 2.6 uV
Mean stored waveform trough at sample 10 (0.31 ms after the waveform start)

Session start 2024-09-27 08:49:52-04:00; Neuralynx header: recording closed 2024-09-27 15:12:24
The raw ephys on the trial clock runs from 2024-09-27 12:06:44.113006-04:00 to 2024-09-27 17:12:20.102974750-04:00
Header open-to-close span 22952 s; raw samples 18336 s

# Contrasts: does each stream respond to the task at V on the trial clock?
rng_check = np.random.default_rng(1)
V = succ.V.to_numpy()
spk_trial = spk + OFFSET
def null_alignments(n, rng):
    """V replaced by a uniform time within each successful trial's span (its start to the next trial's start)."""
    return [rng.uniform(succ.trial_start.to_numpy(), succ.trial_end.to_numpy()) for _ in range(n)]
nulls = null_alignments(100, rng_check)
null_mod = np.array([chamber.modulation(chamber.psth(spk_trial, a)[1]) for a in nulls])
rows = []
for label, extra in [("measured offset", 0.0), ("offset - 5 s", -5.0), ("offset + 5 s", 5.0)]:
    _, c = chamber.psth(spk_trial + extra, V)
    rows.append({"spike clock": label, "modulation": chamber.modulation(c)})
_, c0 = chamber.psth(spk, V)
rows.append({"spike clock": "as stored (no offset)", "modulation": chamber.modulation(c0),
             })
tab = pd.DataFrame(rows)
tab["null 95th pct"] = np.percentile(null_mod, 95)
print("Spikes: variance across 100 ms bins of the trial-averaged rate from V-1 s to V+3 s (spikes/s squared)")
print(tab.round(2).to_string(index=False))
print("Successful trials with V inside the stored spike times' range:", int(np.sum((V > spk.min()) & (V < spk.max()))), "of", len(V))

# LFP: evoked potential at V against random times in the same trials.
def evoked_raw(times, pad=1.0):
    out = []
    for t0 in times:
        k = np.searchsorted(lfp_t, t0 - pad)
        seg = lfp_v[k:k + 2001]
        if len(seg) == 2001 and not np.isnan(seg).any():
            out.append(seg - seg[:1000].mean())
    return np.array(out)
ev_true = evoked_raw(V)
ev_null = [evoked_raw(a) for a in nulls[:30]]
amp = lambda e: np.ptp(e.mean(0)[1000:1500])
print(f"\nLFP: peak-to-peak of the mean LFP 0-0.5 s after V, {len(ev_true)} trials: {amp(ev_true) * 1e6:.1f} uV; "
      f"random times in the same trials, 95th percentile {np.percentile([amp(e) for e in ev_null], 95) * 1e6:.1f} uV")

# Dopamine: when do large and small trials diverge, relative to V?
t_da, da_al = chamber.dopamine_aligned(fscv, succ)
big = (succ.reward == "large").to_numpy()
tt = stats.ttest_ind(da_al[big], da_al[~big], nan_policy="omit").statistic
tt = np.asarray(tt)
run = np.convolve(np.abs(tt) > 3, np.ones(5, int), "valid") == 5
first = t_da[np.argmax(run)] if run.any() else None
print(f"\nDopamine: first time from V at which |t| (large against small) exceeds 3 for 5 samples running: {first:.1f} s;"
      f" largest |t| before V: {np.nanmax(np.abs(tt[t_da < 0])):.2f}")

Spikes: variance across 100 ms bins of the trial-averaged rate from V-1 s to V+3 s (spikes/s squared)
          spike clock  modulation  null 95th pct
      measured offset        0.49           0.06
         offset - 5 s        0.23           0.06
         offset + 5 s        0.04           0.06
as stored (no offset)        0.01           0.06
Successful trials with V inside the stored spike times' range: 399 of 951

LFP: peak-to-peak of the mean LFP 0-0.5 s after V, 951 trials: 18.6 uV; random times in the same trials, 95th percentile 3.7 uV

Dopamine: first time from V at which |t| (large against small) exceeds 3 for 5 samples running: 2.5 s; largest |t| before V: 2.06

# Raw FSCV: can the trial-aligned oxidation current be found in it?
rs = f["acquisition/fscv_response_series"]
ex = f["stimulus/presentation/fscv_excitation_series"]
n_scan = rs["data"].shape[0] // chamber.FSCV_SAMPLES_PER_SCAN
scans = rs["data"][()].reshape(n_scan, chamber.FSCV_SAMPLES_PER_SCAN)
scan_t = rs["timestamps"][::chamber.FSCV_SAMPLES_PER_SCAN]
e0 = ex["data"][:chamber.FSCV_SAMPLES_PER_SCAN]
print(f"Raw FSCV: {n_scan} scans of {chamber.FSCV_SAMPLES_PER_SCAN} samples; scan interval {np.unique(np.round(np.diff(scan_t), 4))} s; "
      f"span {scan_t[0]:.2f} to {scan_t[-1]:.2f} s ({(scan_t[-1] - scan_t[0]) / 60:.1f} min, the trials span "
      f"{(trials.stop_time.iloc[-1] - trials.start_time.iloc[0]) / 60:.1f} min)")
print(f"Applied voltage per scan {e0.min():.3f} to {e0.max():.3f} V (paper: -0.4 to 1.3 V); "
      f"current stated in {rs['data'].attrs['unit']}, values {np.percentile(scans, 1):.0f} to {np.percentile(scans, 99):.0f} (1st to 99th percentile)")
j = int(np.argmin(np.abs(e0[:len(e0) // 2] - 0.6)))
sos = signal.butter(2, [0.02, 0.5], "bandpass", fs=10, output="sos")
cur = signal.sosfiltfilt(sos, scans[:, j - 1:j + 2].mean(1))
energy = np.cumsum(np.r_[0, cur ** 2])
rows = []
for k in np.linspace(0, len(trials) - 1, 7).astype(int):
    y = signal.sosfiltfilt(sos, fscv["iox"][k])
    r = signal.correlate(cur, y, "valid") / np.sqrt((energy[len(y):] - energy[:-len(y)]) * np.sum(y ** 2))
    i = int(np.argmax(np.abs(r)))
    rows.append({"trial": k + 1, "row start s": round(fscv["start_time"][k], 2), "best lag s": round(scan_t[i] - fscv["start_time"][k], 1),
                 "best r": round(r[i], 3), "30th best |r|": round(np.sort(np.abs(r))[-30], 3)})
print(f"\nThe 'iox' row of a trial searched in the raw current at {e0[j]:.2f} V (both 0.02-0.5 Hz):")
print(pd.DataFrame(rows).to_string(index=False))
del scans

Raw FSCV: 196800 scans of 214 samples; scan interval [0.1] s; span 11854.99 to 31534.89 s (328.0 min, the trials span 302.4 min)

Applied voltage per scan -0.347 to 1.340 V (paper: -0.4 to 1.3 V); current stated in amperes, values -1267 to 1197 (1st to 99th percentile)

The 'iox' row of a trial searched in the raw current at 0.59 V (both 0.02-0.5 Hz):
 trial  row start s  best lag s  best r  30th best |r|
     1     11824.99     11787.2  -0.744          0.729
   236     15205.38      6868.9   0.623          0.580
   471     18298.77     -4811.0   0.741          0.708
   706     21183.17      9476.7   0.726          0.671
   941     24234.81      3658.1   0.722          0.650
  1176     27088.46    -10050.4   0.623          0.580
  1411     29908.68     -7502.9   0.706          0.680

t_full, da_full = chamber.dopamine_aligned(fscv, trials)
def sd_ratio(row):
    after = da_full[row, (t_full >= 0) & (t_full <= 5)]
    before = da_full[row, (t_full >= -5) & (t_full < 0)]
    return np.nanstd(after) / np.nanstd(before), np.isnan(after).mean(), np.isnan(before).mean()
for label, trial_id in (("trial id 133 (1-based)", 133), ("index 133 (trial id 134)", 134)):
    row = trial_id - 1
    r, na, nb = sd_ratio(row)
    tr = trials.iloc[row]
    print(f"{label}: reward {tr.reward}, error {tr.error}, V at {tr.V - tr.trial_start:.2f} s, RW at {tr.RW - tr.V:.2f} s from V; "
          f"SD ratio after/before {r:.2f} (NaN after {na:.0%}, before {nb:.0%})")
ratios = np.array([sd_ratio(i - 1)[0] for i in succ.trial])
valid = ratios[~np.isnan(ratios)]
print(f"Reference, {len(valid)} of {len(ratios)} successful trials with a ratio (the others have no reference at V): median ratio "
      f"{np.median(valid):.2f}, 5th-95th percentile {np.percentile(valid, 5):.2f} to {np.percentile(valid, 95):.2f}, above 1 in {np.mean(valid > 1):.0%}")
r133 = sd_ratio(132)[0]
print(f"Trial 133's ratio {r133:.2f} is above {np.mean(valid < r133):.1%} of these {len(valid)} ratios")
fig, ax = plt.subplots(figsize=(7, 3))
for trial_id in (133, 134):
    ax.plot(t_full, da_full[trial_id - 1], label=f"trial {trial_id}")
ax.axvline(0, color="k", lw=0.5)
ax.set(xlabel="time from V (s)", ylabel="dopamine referenced at V (stated unit: none)", xlim=(-10, 10))
ax.legend()
print("Plotted range:", np.round(np.nanmin(da_full[[132, 133]]), 2), "to", np.round(np.nanmax(da_full[[132, 133]]), 2))

trial id 133 (1-based): reward large, error nan, V at 1.90 s, RW at 4.78 s from V; SD ratio after/before 2.02 (NaN after 0%, before 0%)
index 133 (trial id 134): reward small, error nan, V at 1.72 s, RW at 4.73 s from V; SD ratio after/before 2.04 (NaN after 36%, before 46%)

Reference, 889 of 951 successful trials with a ratio (the others have no reference at V): median ratio 1.27, 5th-95th percentile 0.43 to 3.30, above 1 in 64%
Trial 133's ratio 2.02 is above 76.9% of these 889 ratios
Plotted range: -18.32 to 13.59

def window_by_trial(da, t_axis, lo, hi):
    """Mean of each row between per-trial bounds lo[i] and hi[i] (s from V)."""
    out = np.full(len(da), np.nan)
    for i in range(len(da)):
        m = (t_axis >= lo[i]) & (t_axis <= hi[i])
        if np.any(~np.isnan(da[i, m])):
            out[i] = np.nanmean(da[i, m])
    return out
rw = (succ.RW - succ.V).to_numpy()
main = window_by_trial(da_al, t_da, np.zeros(len(succ)), rw)
def compare(values, mask_big, label):
    a, b = values[mask_big], values[~mask_big]
    a, b = a[~np.isnan(a)], b[~np.isnan(b)]
    p = stats.mannwhitneyu(a, b, alternative="two-sided").pvalue
    d = (a.mean() - b.mean()) / np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
    return {"window": label, "n_large": len(a), "n_small": len(b), "mean_large": a.mean(), "mean_small": b.mean(),
            "difference": a.mean() - b.mean(), "cohen_d": d, "mann_whitney_p": p}
rows = [compare(main, big, "V to RW (main)")]
print("Dropped (no dopamine sample in V to RW): large", int(np.isnan(main[big]).sum()), ", small", int(np.isnan(main[~big]).sum()))
rows.append(compare(window_by_trial(da_al, t_da, np.zeros(len(succ)), np.full(len(succ), 2.0)), big, "0 to 2 s"))
rows.append(compare(window_by_trial(da_al, t_da, rw, rw + 3), big, "RW to RW + 3 s"))
rows.append(compare(window_by_trial(da_al, t_da, np.full(len(succ), -3.0), np.zeros(len(succ))), big, "-3 to 0 s (before V)"))
half = (succ.trial <= len(trials) / 2).to_numpy()
for name, m in (("first half of session", half), ("second half of session", ~half)):
    r = compare(main[m], big[m], "V to RW")
    r["window"] = f"V to RW, {name}"
    rows.append(r)
keep = (succ.trial != 1411).to_numpy()
r = compare(main[keep], big[keep], "V to RW")
r["window"] = "V to RW, without trial 1411"
rows.append(r)
print(pd.DataFrame(rows).round(4).to_string(index=False))
# Trials are not independent: dopamine drifts across the session. Labels shuffled within blocks of 50
# consecutive successful trials keep the drift in the null.
rng_8b = np.random.default_rng(8)
blocks = np.arange(len(succ)) // 50
diff_means = lambda a, b: a.mean() - b.mean()
for label, vals in (("V to RW (main)", main), ("RW to RW + 3 s", window_by_trial(da_al, t_da, rw, rw + 3))):
    obs, p_blk = chamber.block_permutation_p(vals, big, blocks, diff_means, rng=rng_8b)
    print(f"{label}: difference of means {obs:.3f}, permutation p within blocks of 50 trials = {p_blk:.4f} (2000 permutations, {blocks.max() + 1} blocks)")
fig, ax = plt.subplots(figsize=(7, 3))
for m, lab, col in ((big, "large", "tab:red"), (~big, "small", "tab:blue")):
    mu = np.nanmean(da_al[m], 0)
    se = np.nanstd(da_al[m], 0) / np.sqrt(np.sum(~np.isnan(da_al[m]), 0))
    ax.plot(t_da, mu, color=col, label=lab)
    ax.fill_between(t_da, mu - se, mu + se, color=col, alpha=0.3)
ax.axvline(0, color="k", lw=0.5)
ax.axvline(np.median(rw), color="k", ls="--", lw=0.5)
ax.set(xlabel="time from V (s); dashed: median RW", ylabel="dopamine referenced at V", xlim=(-10, 10))
ax.legend()
print("Plotted means range:", np.round(np.nanmin([np.nanmean(da_al[big], 0), np.nanmean(da_al[~big], 0)]), 2), "to",
      np.round(np.nanmax([np.nanmean(da_al[big], 0), np.nanmean(da_al[~big], 0)]), 2), "; median RW", round(np.median(rw), 2), "s")

Dropped (no dopamine sample in V to RW): large 23 , small 39
                         window  n_large  n_small  mean_large  mean_small  difference  cohen_d  mann_whitney_p
                 V to RW (main)      429      460      1.8795      1.0955      0.7841   0.1882          0.0035
                       0 to 2 s      429      460      1.3220      1.1808      0.1412   0.0463          0.4395
                 RW to RW + 3 s      426      454      1.0290     -0.5575      1.5864   0.2657          0.0001
           -3 to 0 s (before V)      429      460     -0.1558     -0.3057      0.1499   0.0474          0.4952
 V to RW, first half of session      243      262      2.3919      1.6363      0.7556   0.1760          0.0747
V to RW, second half of session      186      198      1.2102      0.3798      0.8304   0.2131          0.0125
    V to RW, without trial 1411      428      460      1.8871      1.0955      0.7917   0.1900          0.0031

V to RW (main): difference of means 0.784, permutation p within blocks of 50 trials = 0.0045 (2000 permutations, 20 blocks)

RW to RW + 3 s: difference of means 1.586, permutation p within blocks of 50 trials = 0.0010 (2000 permutations, 20 blocks)
Plotted means range: -2.51 to 3.82 ; median RW 4.77 s

rng_9a = np.random.default_rng(9)
def fig9a_test(spike_times, label):
    centers, counts = chamber.psth(spike_times, V)
    has_spike = (counts.sum(1) > 0).mean()
    obs = chamber.modulation(counts)
    null = np.array([chamber.modulation(chamber.psth(spike_times, a)[1]) for a in null_alignments(200, rng_9a)])
    return {"spikes": label, "trials with a spike": round(has_spike, 4), "modulation": round(obs, 2),
            "null 95th pct": round(np.percentile(null, 95), 2), "null max": round(null.max(), 2),
            "p (random alignments)": (np.sum(null >= obs) + 1) / (len(null) + 1)}, centers, counts
res_main, centers, counts = fig9a_test(spk + OFFSET, f"thresholded unit {main_unit}, measured offset (main)")
res_plx, _, counts_plx = fig9a_test(plx_spk + OFFSET, "Plexon unit, measured offset")
res_raw, _, _ = fig9a_test(spk, f"thresholded unit {main_unit}, as stored")
res_lag, _, _ = fig9a_test(spk + OFFSET + SPIKE_LAG, f"thresholded unit {main_unit}, offset and 7 ms lag")
print(pd.DataFrame([res_main, res_plx, res_raw, res_lag]).to_string(index=False))
# Variants: other events and a longer window, 250 ms bins.
rows = []
fixation = events["value cue fixation started"]
v_fix = fixation[np.searchsorted(fixation, V)]
for label, align, window in (("V, -1 to 6 s", V, (-1.0, 6.0)), ("RW, -2 to 3 s", succ.RW.to_numpy(), (-2.0, 3.0)),
                             ("value cue fixation, -2 to 3 s", v_fix, (-2.0, 3.0))):
    _, c = chamber.psth(spk + OFFSET, align, window=window, bin_s=0.25)
    # null: the same random times as for V, moved by each trial's event-to-V interval
    nm = [chamber.modulation(chamber.psth(spk + OFFSET, a + (align - V), window=window, bin_s=0.25)[1], 0.25)
          for a in null_alignments(100, rng_9a)]
    rows.append({"aligned to": label, "modulation": chamber.modulation(c, 0.25), "null 95th pct": np.percentile(nm, 95)})
print("\nVariants, 250 ms bins:")
print(pd.DataFrame(rows).round(3).to_string(index=False))
rate = counts / 0.1
for m, lab in ((big, "large"), (~big, "small")):
    r = rate[m].mean(0)
    print(f"{lab}: {m.sum()} trials; mean rate before V (-1 to 0 s) {rate[m][:, centers < 0].mean():.2f} spikes/s, "
          f"0 to 0.5 s {rate[m][:, (centers > 0) & (centers < 0.5)].mean():.2f}, peak {r.max():.2f} at {centers[r.argmax()]:.2f} s")
dip = rate[:, (centers > 0) & (centers < 0.2)].mean()
print(f"All successful trials: mean rate -1 to 0 s {rate[:, centers < 0].mean():.2f}, 0 to 0.2 s {dip:.2f}, 0.3 to 0.5 s "
      f"{rate[:, (centers > 0.3) & (centers < 0.5)].mean():.2f} spikes/s")
post = rate[:, (centers > 0) & (centers < 1)].mean(1)
print(f"Large against small, mean rate 0 to 1 s after V: {post[big].mean():.2f} against {post[~big].mean():.2f} spikes/s, "
      f"Mann-Whitney p = {stats.mannwhitneyu(post[big], post[~big]).pvalue:.3g} (descriptive)")
psth_tab = pd.DataFrame({"bin centre (s)": centers, "large (spikes/s)": rate[big].mean(0), "small (spikes/s)": rate[~big].mean(0)})
print("\nPSTH per condition, 100 ms bins:")
print(psth_tab.round(2).T.to_string(header=False))
fig, axs = plt.subplots(2, 1, figsize=(7, 6), sharex=True)
order = np.r_[np.where(big)[0], np.where(~big)[0]]
for row, i in enumerate(order):
    lo, hi = np.searchsorted(spk + OFFSET, [V[i] - 1, V[i] + 3])
    axs[0].plot((spk + OFFSET)[lo:hi] - V[i], np.full(hi - lo, row), "|", ms=1, color="tab:red" if big[i] else "tab:blue")
axs[0].set(ylabel="trial (large, then small)")
for m, lab, col in ((big, "large", "tab:red"), (~big, "small", "tab:blue")):
    axs[1].plot(centers, rate[m].mean(0), color=col, label=lab)
axs[1].axvline(0, color="k", lw=0.5)
axs[1].set(xlabel="time from V (s)", ylabel="spikes/s")
axs[1].legend()
print("PSTH range:", round(min(rate[big].mean(0).min(), rate[~big].mean(0).min()), 2), "to", round(max(rate[big].mean(0).max(), rate[~big].mean(0).max()), 2), "spikes/s")

                                    spikes  trials with a spike  modulation  null 95th pct  null max  p (random alignments)
thresholded unit 1, measured offset (main)               1.0000        0.49           0.06      0.07               0.004975
              Plexon unit, measured offset               0.9085        0.17           0.06      0.07               0.004975
             thresholded unit 1, as stored               0.4196        0.01           0.03      0.04               0.905473
   thresholded unit 1, offset and 7 ms lag               1.0000        0.48           0.06      0.08               0.004975

Variants, 250 ms bins:
                   aligned to  modulation  null 95th pct
                 V, -1 to 6 s       0.258          0.020
                RW, -2 to 3 s       0.123          0.020
value cue fixation, -2 to 3 s       0.299          0.023
large: 452 trials; mean rate before V (-1 to 0 s) 5.66 spikes/s, 0 to 0.5 s 5.68, peak 8.50 at 0.45 s
small: 499 trials; mean rate before V (-1 to 0 s) 5.59 spikes/s, 0 to 0.5 s 5.90, peak 8.04 at 0.35 s
All successful trials: mean rate -1 to 0 s 5.63, 0 to 0.2 s 4.58, 0.3 to 0.5 s 7.70 spikes/s
Large against small, mean rate 0 to 1 s after V: 5.69 against 5.66 spikes/s, Mann-Whitney p = 0.896 (descriptive)

PSTH per condition, 100 ms bins:
bin centre (s)   -0.95 -0.85 -0.75 -0.65 -0.55 -0.45 -0.35 -0.25 -0.15 -0.05  0.05  0.15  0.25  0.35  0.45  0.55  0.65  0.75  0.85  0.95  1.05  1.15  1.25  1.35  1.45  1.55  1.65  1.75  1.85  1.95  2.05  2.15  2.25  2.35  2.45  2.55  2.65  2.75  2.85  2.95
large (spikes/s)  5.84  5.91  5.80  5.02  5.55  5.64  6.00  5.73  5.62  5.53  4.96  4.00  3.54  7.41  8.50  8.27  5.77  4.93  4.29  5.27  5.46  4.67  5.18  4.96  4.96  4.93  5.40  4.89  4.69  4.82  6.13  5.18  5.64  5.60  5.18  5.42  5.44  5.71  5.53  5.40
small (spikes/s)  5.07  6.27  5.55  5.71  5.47  5.71  5.49  5.41  5.51  5.73  5.25  4.09  5.21  8.04  6.89  6.23  5.21  5.41  4.79  5.45  4.99  5.43  5.47  5.01  5.25  5.59  5.79  5.23  5.13  5.37  5.75  5.55  4.85  5.99  5.53  5.19  5.61  5.57  5.13  5.61

PSTH range: 3.54 to 8.5 spikes/s

wf_all = chamber.thresholded_waveforms(f, list(thr).index(main_unit))
isi = np.diff(spk)
mean_wf, sd_wf = wf_all.mean(0), wf_all.std(0)
k = int(np.argmax(np.abs(mean_wf)))
rows = []
for label, times, w in ((f"thresholded unit {main_unit} (main)", spk, wf_all), ("Plexon unit", plx_spk, None)):
    i = np.diff(times)
    row = {"spikes": label, "n": len(times), "n ISI < 1 ms": int(np.sum(i < 0.001)), "frac ISI < 1 ms": np.mean(i < 0.001),
           "n ISI < 2 ms": int(np.sum(i < 0.002)), "frac ISI < 2 ms": np.mean(i < 0.002), "median ISI ms": np.median(i) * 1000}
    if w is not None:
        row.update({"mean wf extreme (uV)": mean_wf[k], "SD there (uV)": sd_wf[k], "ratio": abs(mean_wf[k]) / sd_wf[k]})
    row["shortest ISI ms"] = i.min() * 1000
    rows.append(row)
print(pd.DataFrame(rows).round(6).to_string(index=False))
dead = isi.min()
amp_all = np.abs(wf_all).max(1)
print(f"All {len(wf_all)} spikes: fraction whose largest |value| exceeds 40: {np.mean(amp_all > 40):.6f}; smallest largest |value| {amp_all.min():.1f}")
print(f"\nDetector dead time: shortest ISI of thresholded unit {main_unit} {dead * 1000:.3f} ms; the stored waveform spans "
      f"{wf_all.shape[1]} samples ({wf_all.shape[1] / 32:.2f} ms)")
plx_isi = np.diff(plx_spk)
print(f"Plexon unit: ISIs under {dead * 1000:.2f} ms {np.sum(plx_isi < dead)} ({np.mean(plx_isi < dead):.4%})")
# Refractoriness beyond the dead time: ISI count per ms just above it against a stretch where it is flat.
for label, i in ((f"thresholded unit {main_unit}", isi), ("Plexon unit", plx_isi)):
    short = np.sum((i >= dead) & (i < 0.003)) / ((0.003 - dead) * 1000)
    flat = np.sum((i >= 0.01) & (i < 0.03)) / 20
    print(f"{label}: ISIs per ms in {dead * 1000:.2f}-3 ms {short:.1f}, in 10-30 ms {flat:.1f}, ratio {short / flat:.3f}")
in_trials = np.zeros(len(spk), bool)
for a in V:
    lo, hi = np.searchsorted(spk + OFFSET, [a - 1, a + 3])
    in_trials[lo:hi] = True
print(f"ISI < 1 ms among spikes within V-1 to V+3 s of successful trials: {np.mean(np.diff(spk[in_trials]) < 0.001):.4f} "
      "(pairs across trial windows included)")
fig, axs = plt.subplots(1, 2, figsize=(9, 3))
t_wf = np.arange(wf_all.shape[1]) / 32.0
axs[0].plot(t_wf, mean_wf, color="k")
axs[0].fill_between(t_wf, mean_wf - sd_wf, mean_wf + sd_wf, color="0.7")
axs[0].set(xlabel="ms", ylabel="uV")
axs[1].hist(isi[isi < 0.1] * 1000, bins=np.arange(0, 101, 1), color="0.4")
axs[1].set(xlabel="ISI (ms), 1 ms bins", ylabel="count")
print(f"Waveform plot range {np.min(mean_wf - sd_wf):.1f} to {np.max(mean_wf + sd_wf):.1f} uV; "
      f"ISI histogram shows {np.mean(isi < 0.1):.1%} of ISIs (those under 100 ms)")

                   spikes      n  n ISI < 1 ms  frac ISI < 1 ms  n ISI < 2 ms  frac ISI < 2 ms  median ISI ms  mean wf extreme (uV)  SD there (uV)    ratio  shortest ISI ms
thresholded unit 1 (main) 101733             0         0.000000             9         0.000088        141.938             -56.35675      11.763006 4.791016            1.562
              Plexon unit  90912             1         0.000011            39         0.000429        139.988                   NaN            NaN      NaN            0.250
All 101733 spikes: fraction whose largest |value| exceeds 40: 1.000000; smallest largest |value| 40.3

Detector dead time: shortest ISI of thresholded unit 1 1.562 ms; the stored waveform spans 48 samples (1.50 ms)
Plexon unit: ISIs under 1.56 ms 14 (0.0154%)
thresholded unit 1: ISIs per ms in 1.56-3 ms 15.3, in 10-30 ms 144.7, ratio 0.106
Plexon unit: ISIs per ms in 1.56-3 ms 27.1, in 10-30 ms 132.9, ratio 0.204
ISI < 1 ms among spikes within V-1 to V+3 s of successful trials: 0.0000 (pairs across trial windows included)
Waveform plot range -68.1 to 44.9 uV; ISI histogram shows 36.4% of ISIs (those under 100 ms)

def beta_test(band, post, base, label):
    t_ax, pw = chamber.band_power_around(lfp_t, lfp_v, V, band=band)
    ratio = np.log10(chamber.window_mean(t_ax, pw, *post) / chamber.window_mean(t_ax, pw, *base))
    a, b = ratio[big], ratio[~big]
    dropped = (int(np.isnan(a).sum()), int(np.isnan(b).sum()))
    a, b = a[~np.isnan(a)], b[~np.isnan(b)]
    p = stats.mannwhitneyu(a, b, alternative="two-sided").pvalue
    w_a, w_b = stats.wilcoxon(a).pvalue, stats.wilcoxon(b).pvalue
    return {"variant": label, "n_large": len(a), "n_small": len(b), "dropped (l, s)": dropped,
            "median_large": np.median(a), "median_small": np.median(b), "diff_medians": np.median(a) - np.median(b),
            "mann_whitney_p": p, "p large != 0": w_a, "p small != 0": w_b}, t_ax, pw
rows = []
r, t_beta, pw_beta = beta_test((13, 30), (0, 0.5), (-0.5, 0), "13-30 Hz, 0-0.5 s vs 0.5 s before (main)")
rows.append(r)
for band, post, base, label in [((13, 30), (0, 0.3), (-0.5, 0), "post 0-0.3 s"), ((13, 30), (0, 1.0), (-0.5, 0), "post 0-1 s"),
                                ((13, 30), (0, 0.5), (-1.0, 0), "1 s baseline"), ((15, 30), (0, 0.5), (-0.5, 0), "15-30 Hz"),
                                ((12, 25), (0, 0.5), (-0.5, 0), "12-25 Hz")]:
    rows.append(beta_test(band, post, base, label)[0])
print(pd.DataFrame(rows).round(4).to_string(index=False))
raw_post = chamber.window_mean(t_beta, pw_beta, 0, 0.5)
ok = ~np.isnan(raw_post)
print(f"Variant without baseline: mean beta power 0-0.5 s, large {np.median(raw_post[big & ok]):.3e} against small {np.median(raw_post[~big & ok]):.3e} V^2 (medians), "
      f"Mann-Whitney p = {stats.mannwhitneyu(raw_post[big & ok], raw_post[~big & ok]).pvalue:.3g}")
for name, m in (("first half of session", half), ("second half of session", ~half)):
    t_ax, pw = t_beta, pw_beta[m]
    ratio = np.log10(chamber.window_mean(t_ax, pw, 0, 0.5) / chamber.window_mean(t_ax, pw, -0.5, 0))
    a, b = ratio[big[m]], ratio[~big[m]]
    a, b = a[~np.isnan(a)], b[~np.isnan(b)]
    print(f"{name}: medians large {np.median(a):.4f}, small {np.median(b):.4f}, p = {stats.mannwhitneyu(a, b).pvalue:.3g}")
# Trial-averaged normalized power over 0 to 0.5 s (the paper's figure shows trial averages), with
# a bootstrap 95% CI over trials.
rng_9d = np.random.default_rng(4)
norm_post = chamber.window_mean(t_beta, pw_beta, 0, 0.5) / chamber.window_mean(t_beta, pw_beta, -0.5, 0)
okn = ~np.isnan(norm_post)
na, nb_ = norm_post[big & okn], norm_post[~big & okn]
boot = np.array([[rng_9d.choice(na, len(na)).mean(), rng_9d.choice(nb_, len(nb_)).mean()] for _ in range(2000)])
for j, (lab, x) in enumerate((("large", na), ("small", nb_))):
    print(f"Trial-averaged beta power 0-0.5 s / 0.5 s before, {lab}: {x.mean():.3f} (95% CI {np.percentile(boot[:, j], 2.5):.3f} to "
          f"{np.percentile(boot[:, j], 97.5):.3f})")
dd = boot[:, 0] - boot[:, 1]
print(f"Difference large - small: {na.mean() - nb_.mean():.3f} (95% CI {np.percentile(dd, 2.5):.3f} to {np.percentile(dd, 97.5):.3f})")
lr = np.log10(norm_post)
obs, p_blk = chamber.block_permutation_p(lr, big, np.arange(len(succ)) // 50, lambda a, b: np.median(a) - np.median(b), rng=rng_9d)
print(f"Main measure, difference of medians {obs:.4f}, permutation p within blocks of 50 trials = {p_blk:.4f}")
fig, ax = plt.subplots(figsize=(7, 3))
for m, lab, col in ((big, "large", "tab:red"), (~big, "small", "tab:blue")):
    base = chamber.window_mean(t_beta, pw_beta[m], -0.5, 0)
    norm = pw_beta[m] / base[:, None]
    ax.plot(t_beta, np.nanmean(norm, 0), color=col, label=lab)
ax.axvline(0, color="k", lw=0.5)
ax.set(xlabel="time from V (s)", ylabel="beta power / 0.5 s before V", xlim=(-1.5, 1.5))
ax.legend()
lines = [np.nanmean(pw_beta[m] / chamber.window_mean(t_beta, pw_beta[m], -0.5, 0)[:, None], 0) for m in (big, ~big)]
sel = (t_beta >= -1.5) & (t_beta <= 1.5)
print("Plotted range (-1.5 to 1.5 s):", round(min(l[sel].min() for l in lines), 3), "to", round(max(l[sel].max() for l in lines), 3))
print("At the segment edges, -3 s and +3 s (filter edge effects, outside the plot):", [np.round(l[[0, -1]], 3).tolist() for l in lines])

                                 variant  n_large  n_small dropped (l, s)  median_large  median_small  diff_medians  mann_whitney_p  p large != 0  p small != 0
13-30 Hz, 0-0.5 s vs 0.5 s before (main)      452      499         (0, 0)       -0.1577       -0.1247       -0.0330          0.1080           0.0           0.0
                            post 0-0.3 s      452      499         (0, 0)       -0.0923       -0.0934        0.0011          0.5872           0.0           0.0
                              post 0-1 s      452      499         (0, 0)       -0.1112       -0.1031       -0.0082          0.8371           0.0           0.0
                            1 s baseline      452      499         (0, 0)       -0.1526       -0.1274       -0.0251          0.1081           0.0           0.0
                                15-30 Hz      452      499         (0, 0)       -0.1281       -0.1012       -0.0269          0.0627           0.0           0.0
                                12-25 Hz      452      499         (0, 0)       -0.1949       -0.1763       -0.0186          0.1936           0.0           0.0
Variant without baseline: mean beta power 0-0.5 s, large 7.303e-11 against small 7.884e-11 V^2 (medians), Mann-Whitney p = 0.169
first half of session: medians large -0.1832, small -0.1331, p = 0.117
second half of session: medians large -0.1266, small -0.1149, p = 0.52
Trial-averaged beta power 0-0.5 s / 0.5 s before, large: 0.839 (95% CI 0.782 to 0.903)
Trial-averaged beta power 0-0.5 s / 0.5 s before, small: 0.885 (95% CI 0.831 to 0.940)
Difference large - small: -0.046 (95% CI -0.125 to 0.036)

Main measure, difference of medians -0.0330, permutation p within blocks of 50 trials = 0.1864
Plotted range (-1.5 to 1.5 s): 0.531 to 1.295
At the segment edges, -3 s and +3 s (filter edge effects, outside the plot): [[0.262, 0.159], [0.317, 0.194]]

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

Notebook ran in 1.5 min

Claim	Figure	Tested
fig8a-single-trial-dopamine	8A	yes
fig8b-dopamine-large-vs-small	8B	yes
fig9a-spikes-task-aligned	9A	yes
fig9bc-unit-waveform-isi	9B, 9C	yes
fig9d-beta-suppression-large	9D	yes
recordings-two-monkeys	Abstract, 3.5	blocked: data missing (monkey J)
scfe-functional-yield	3.5	blocked: data missing
scfe-background-yield	3.5	blocked: data missing
mcfe-functional-yield	3.5	blocked: data missing
dopamine-signal-yield	3.5	blocked: data missing
asepsis-cultures	3.2	blocked: data missing
fig6-drainage-volume	6	blocked: data missing
fig7-bone-growth-margins	7	blocked: data missing

Claim	Outcome
fig8a-single-trial-dopamine	reproduced (weak test)
fig8b-dopamine-large-vs-small	reproduced
fig9a-spikes-task-aligned	reproduced (after correcting the spike clock)
fig9bc-unit-waveform-isi	reproduced by the criterion (weak evidence for a single unit)
fig9d-beta-suppression-large	partially reproduced
recordings-two-monkeys	blocked: data missing
scfe-functional-yield	blocked: data missing
scfe-background-yield	blocked: data missing
mcfe-functional-yield	blocked: data missing
dopamine-signal-yield	blocked: data missing
asepsis-cultures	blocked: data missing
fig6-drainage-volume	blocked: data missing
fig7-bone-growth-margins	blocked: data missing

Reproducing Choi et al. (2025) from dandiset 002014¶

Data check¶

fig8a-single-trial-dopamine¶

fig8b-dopamine-large-vs-small¶

fig9a-spikes-task-aligned¶

fig9bc-unit-waveform-isi¶

fig9d-beta-suppression-large¶

Summary¶