OSCR

Protocol for analyzing slow cortical dynamics in mouse neuronal recordings.

Code ↔ Paper

30 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 30 matches
  1. [1] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.ipynb, lines 36–62 · score 0.78 · 0–1.2 s, Responsive cells, stimulus onset, stimulus triggered, threshold, post
  2. [2] § Step-by-step method details › Trial-to-trial variability analysis ↔ functions/sd_utils.py, lines 380–427 · score 0.71 · stimulus response window, peak exceeds, Responsive cells, threshold, 2 s
  3. [3] § Step-by-step method details › Intrinsic timescale analysis ↔ functions/sd_utils.py, lines 782–809 · score 0.65 · tau_cell, tau_net, network intrinsic timescales, sd
  4. [4] § Step-by-step method details › Intrinsic timescale analysis ↔ functions/sd_utils.py, lines 782–809 · score 0.65 · tau_cell, tau_net, network intrinsic timescales, sd
  5. [5] § Step-by-step method details › Intrinsic timescale analysis ↔ slow_dynamics_net_int.ipynb, lines 55–66 · score 0.64 · sd.get_network_intrinsic_timescales, tau_cell, tau_net
  6. [6] § Step-by-step method details › Intrinsic timescale analysis ↔ slow_dynamics_net_int.ipynb, lines 55–66 · score 0.64 · sd.get_network_intrinsic_timescales, tau_cell, tau_net
  7. [7] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.ipynb, lines 36–62 · score 0.64 · plot_t_tuning, trial_frames_tuning, sd.get_frames
  8. [8] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.py, lines 41–65 · score 0.64 · plot_t_tuning, trial_frames_tuning, sd.get_frames
  9. [9] § Step-by-step method details › Intrinsic timescale analysis ↔ functions/sd_utils.py, lines 505–511 · score 0.64 · Euclidean distance, population vector, baseline, Firing rates, network, frame
  10. [10] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation_comb.ipynb, lines 167–190 · score 0.63 · mouse_idx2, trial_idx, n_freq, echo
  11. [11] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.ipynb, lines 114–132 · score 0.63 · mouse_idx2, trial_idx, n_freq, echo
  12. [12] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.ipynb, lines 114–132 · score 0.62 · mouse_list, mouse_idx2, mouse_tag
  13. [13] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.py, lines 125–139 · score 0.62 · mouse_list, mouse_idx2, mouse_tag
  14. [14] § Step-by-step method details › Intrinsic timescale analysis ↔ slow_dynamics_net_int.ipynb, lines 34–53 · score 0.61 · num_initial_trials_skip, flatten_runs, untrained, firing rate, oddball, network
  15. [15] § Step-by-step method details › Intrinsic timescale analysis ↔ slow_dynamics_net_int.py, lines 39–58 · score 0.61 · num_initial_trials_skip, flatten_runs, untrained, firing rate, oddball, network
  16. [16] § Step-by-step method details › Data loading and preprocessing ↔ slow_dynamics_decoder_comb_data_rnn.ipynb, lines 86–101 · score 0.60 · trial_frames_rnn, plot_t_rnn, sd.get_frames
  17. [17] § Step-by-step method details › Data loading and preprocessing ↔ slow_dynamics_decoder_comb.py, lines 91–96 · score 0.60 · trial_frames_rnn, plot_t_rnn, sd.get_frames
  18. [18] § Step-by-step method details › Data loading and preprocessing ↔ slow_dynamics_decoder.ipynb, lines 1–14 · score 0.59 · import functions, pipeline_dir, slow dynamics
  19. [19] § Step-by-step method details › Data loading and preprocessing ↔ slow_dynamics_decoder_comb_data_rnn.ipynb, lines 1–12 · score 0.59 · import functions, pipeline_dir, slow dynamics
  20. [20] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.ipynb, lines 134–143 · score 0.58 · n_freq, mouse_tag, title_tag
  21. [21] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.py, lines 143–143 · score 0.58 · n_freq, mouse_tag, title_tag
  22. [22] § Step-by-step method details › Trial-to-trial variability analysis ↔ slow_dynamics_isi_correlation.py, lines 125–139 · score 0.58 · mouse_idx2, trial_idx, stim_trig_resp
  23. [23] § Step-by-step method details › Intrinsic timescale analysis ↔ slow_dynamics_net_int.ipynb, lines 82–110 · score 0.57 · Freq neuron, sd.stat_compare, CaIm, Intrinsic, timescale
  24. [24] § Step-by-step method details › Intrinsic timescale analysis ↔ slow_dynamics_net_int.ipynb, lines 82–110 · score 0.57 · Ob neuron, sd.stat_compare, CaIm, Intrinsic, timescale
  25. [25] § Step-by-step method details › Data loading and preprocessing ↔ functions/sd_utils.py, lines 57–183 · score 0.56 · CaImAn, derivative, rectified, oasis, firing rate, smoothdfdt
  26. [26] § Step-by-step method details › Intrinsic timescale analysis ↔ functions/sd_utils.py, lines 513–563 · score 0.55 · zero lag, fit, autocorrelation, trace, Intrinsic, timescale
  27. [27] § Step-by-step method details › Binwise decoder ↔ functions/sd_decoder.py, lines 90–141 · score 0.55 · temporal generalization, decoder trained, cross, Binwise
  28. [28] § Step-by-step method details › Binwise decoder ↔ slow_dynamics_decoder.ipynb, lines 71–92 · score 0.51 · run_binwise_dec, cross validation, diagonal, decoder
  29. [29] § Step-by-step method details › Binwise decoder ↔ slow_dynamics_decoder.py, lines 67–88 · score 0.51 · run_binwise_dec, cross validation, diagonal, decoder
  30. [30] § Before you begin › Data preparation ↔ slow_dynamics_decoder_comb_data_rnn.ipynb, lines 36–55 · score 0.50 · neuronal activity, RNNs trained, untrained, oddball, network, stimulus

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,678 lines · 69 KB · no license · 6 matches

  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Wed Nov 26 13:32:23 2025
  5. @author: ys2605
  6. """
  8. import os
  9. import h5py
  10. import numpy as np
  11. from scipy.stats import norm, wilcoxon, mannwhitneyu, ttest_ind, ttest_rel, f_oneway
  12. from scipy.spatial.distance import pdist, squareform, cdist
  13. from scipy.signal import correlate #, correlation_lags
  15. import matplotlib.pyplot as plt
  16. from matplotlib.ticker import FuncFormatter, NullFormatter
  17. from datetime import datetime
  18. from sd_decoder import plot_diag_binwise_dec, plot_full_binwise_dec
  20. #%%
  21. def get_fnames_from_dir(data_dir, ext_list = [], tags = None, f_list=None):
  22. # list files in data_dir, keeping those with an ext_list extension AND containing every tag; f_list overrides the listing
  24. if f_list is None:
  25. if not os.path.isdir(data_dir):
  26. raise FileNotFoundError('Data directory not found: %s' % data_dir)
  27. f_list = os.listdir(data_dir)
  28. f_list2 = []
  29. for fil1 in f_list:
  30. if len(ext_list):
  31. for ext1 in ext_list:
  32. if fil1.endswith(ext1):
  33. f_list2.append(fil1)
  34. else:
  35. f_list2.append(fil1)
  37. if tags is not None:
  38. if type(tags) is str:
  39. tags = [tags]
  41. f_list_out = []
  42. for fil1 in f_list2:
  43. has_tag = True
  44. for tag in tags:
  45. if tag not in fil1:
  46. has_tag = False
  47. if has_tag:
  48. f_list_out.append(fil1)
  50. else:
  51. f_list_out = f_list2
  53. return f_list_out
  55. #%%
  57. def load_caim_data_mat(data_dir, ext_list = [], tags = None, num_files=None, data_tag = 'results_cnmf_sort.mat', proc_tag = 'processed_data.mat', deconvolution='oasis', smooth_std_duration=0.1, norm_first=False):
  59. # load matched caiman .mat sessions from data_dir into a list of dataset dicts (one per session);
  60. # main fields: firing_rates (neurons x time), trial_types, stim_times, volume_period, isi.
  61. # deconvolution methods are either oasis (caiman default) or smoothdftd - smoothed, rectified first derivative
  62. # norm_first: False (default) = smooth then peak-normalize (original); True = peak-normalize raw S then smooth (MATLAB order)
  63. # smooth_std_duration in sec
  66. dir_files = get_fnames_from_dir(data_dir)
  67. flist = get_fnames_from_dir(data_dir, ext_list=ext_list, tags=tags, f_list=dir_files)
  69. if len(flist) == 0:
  70. print('Warning: no files found in %s matching ext %s and tags %s' % (data_dir, ext_list, tags))
  71. return []
  73. if num_files is None:
  74. num_files = len(flist)
  75. else:
  76. num_files = np.min([len(flist), num_files])
  78. data_out = []
  80. for n_fl in range(num_files):
  82. fname_core = flist[n_fl]
  83. if data_tag in fname_core:
  84. fname_core = fname_core.removesuffix(data_tag)
  85. if proc_tag in fname_core:
  86. fname_core = fname_core.removesuffix(proc_tag)
  88. flist_data = get_fnames_from_dir(data_dir, ext_list = ['.mat'], tags = [fname_core, data_tag], f_list=dir_files)
  89. flist_proc = get_fnames_from_dir(data_dir, ext_list = ['.mat'], tags = [fname_core, proc_tag], f_list=dir_files)
  91. do_load = False
  92. if len(flist_proc):
  93. if len(flist_data):
  94. do_load = True
  95. else:
  96. print(fname_core + " data file with " + data_tag + " tag not found, skipping")
  97. else:
  98. print(fname_core + " proc file with " + proc_tag + " tag not found, skipping")
  100. if do_load:
  102. data_slice = {'flist_data': flist_data,
  103. 'flist_proc': flist_proc}
  105. f_proc = h5py.File(data_dir + '/' + flist_proc[0], 'r')
  106. vid_cuts_trace = f_proc[f_proc['data']['file_cuts_params'][0][0]]['vid_cuts_trace'][()].flatten().astype(bool)
  107. trial_types = f_proc['data']['trial_types'][()].flatten().astype(int)
  108. stim_times = f_proc[f_proc['data']['stim_times_frame'][0][0]][()].flatten().astype(int)
  110. if 'volume_period' in f_proc['data']['frame_data'].keys():
  111. data_slice['volume_period'] = f_proc['data']['frame_data']['volume_period'][()].flatten()[0]
  112. if data_slice.get('volume_period', 0) <= 0:
  113. print('Warning: %s missing or invalid volume_period in _processed_data; using default 33.3 ms (~30 Hz).' % fname_core)
  114. data_slice['volume_period'] = 33.3 # ms, default ~30 Hz
  115. if 'isi' in f_proc['data']['stim_params'].keys():
  116. data_slice['isi'] = f_proc['data']['stim_params']['isi'][()].flatten()[0]
  117. if 'MMN_orientations' in f_proc['data'].keys():
  118. data_slice['MMN_ori'] = f_proc['data']['MMN_orientations'][()].flatten().astype(int)
  119. if 'MMN_freq' in f_proc['data']['stim_params'].keys():
  120. data_slice['MMN_ori'] = f_proc['data']['stim_params']['MMN_freq'][()].flatten().astype(int)
  122. f_proc.close()
  124. firing_rates_all = []
  125. dset_idx_all = []
  127. for n_fl in range(len(flist_data)):
  128. fname_data = flist_data[n_fl]
  130. f = h5py.File(os.path.join(data_dir, fname_data), 'r')
  132. d_est = f['est']
  133. d_proc = f['proc']
  135. comp_acc = d_proc['comp_accepted'][()].flatten().astype(bool)
  137. if deconvolution == 'oasis':
  138. firing_rates_cut = d_est['S'][()][:,comp_acc].T
  140. elif deconvolution == 'smoothdfdt':
  141. C = d_est['C'][()]
  142. YrA = d_est['YrA'][()]
  143. ca_traces_cut = (C + YrA)[:,comp_acc].T
  144. firing_rates_cut = smooth_dfdt(ca_traces_cut, sigma_frames=1000/data_slice['volume_period']*0.1, do_smooth=True)
  146. sigma_fr = 1000/data_slice['volume_period']*smooth_std_duration
  147. if norm_first:
  148. # MATLAB order: peak-normalize the raw deconvolved S first, then smooth
  149. peak_rate = np.max(firing_rates_cut, axis=1)[:,None]
  150. peak_rate[peak_rate == 0] = 1
  151. firing_rates_cutn = gauss_smooth(firing_rates_cut/peak_rate, sigma_frames=sigma_fr)
  152. else:
  153. # original order: smooth, then peak-normalize the smoothed trace
  154. firing_rates_cut = gauss_smooth(firing_rates_cut, sigma_frames=sigma_fr)
  155. peak_rate = np.max(firing_rates_cut, axis=1)[:,None]
  156. peak_rate[peak_rate == 0] = 1
  157. firing_rates_cutn = firing_rates_cut/peak_rate
  159. firing_rates = np.zeros((firing_rates_cutn.shape[0], vid_cuts_trace.shape[0]))
  160. firing_rates[:, vid_cuts_trace] = firing_rates_cutn
  162. firing_rates_all.append(firing_rates)
  163. dset_idx_all.append(np.ones(firing_rates.shape[0], dtype=int)*(n_fl))
  164. f.close()
  166. firing_rates = np.vstack(firing_rates_all)
  167. dset_idx = np.hstack(dset_idx_all)
  169. data_slice['fname_core'] = fname_core
  170. data_slice['firing_rates'] = firing_rates
  171. data_slice['trial_types'] = trial_types
  172. data_slice['stim_times'] = stim_times
  173. data_slice['vid_cuts_trace'] = vid_cuts_trace
  174. data_slice['files_loaded'] = do_load
  175. data_slice['dset_idx'] = dset_idx
  177. data_out.append(data_slice)
  179. if len(data_out):
  180. mouse_ids = list(np.unique(get_mouse_id(data_out)))
  181. print('Loaded %d datasets from %d mice: %s' % (len(data_out), len(mouse_ids), mouse_ids))
  183. return data_out
  185. def load_caim_data_mat2(*args, **kwargs):
  186. # same as load_caim_data_mat but with the MATLAB preprocessing order:
  187. # peak-normalize the raw deconvolved S first, then smooth (norm_first=True)
  188. kwargs['norm_first'] = True
  189. return load_caim_data_mat(*args, **kwargs)
  191. def h5_load_group(group, keys=None):
  192. # load the given keys (default: all) of an open h5py group into a plain dict
  193. if keys is None:
  194. keys = group.keys()
  196. data = {}
  197. for key1 in keys:
  198. data[key1] = group[key1][()]
  200. return data
  202. def get_values(data_out, key):
  203. # collect data_slice[key] across the datasets in the list (skips slices missing key)
  205. values = []
  207. for data_slice in data_out:
  208. if key in data_slice.keys():
  209. values.append(data_slice[key])
  211. return values
  213. def get_mouse_id(data_echo):
  214. # mouse id per dataset = text before the first underscore of fname_core (e.g. 'M4372')
  216. fnames = get_values(data_echo, 'fname_core')
  218. mouse_id = []
  219. for fname in fnames:
  220. mouse_id.append(fname.split('_')[0])
  222. return mouse_id
  224. def get_example_mouse(data_echo, example_mouse='M4372'):
  225. # return (per-dataset mouse-id array, example_mouse); error if example_mouse is not among the loaded datasets
  227. mouse_list = np.array(get_mouse_id(data_echo))
  228. mouse_uq = np.unique(mouse_list)
  230. if example_mouse not in mouse_uq:
  231. raise ValueError('Example mouse %s not found among loaded datasets. The trial-to-trial '
  232. 'similarity demo requires all six variable-ISI (echo) datasets: '
  233. 'M226, M4264, M4265, M4266, M4371, M4372. Loaded mice: %s'
  234. % (example_mouse, list(mouse_uq)))
  236. return mouse_list, example_mouse
  238. def smooth_dfdt(data, do_smooth=True, sigma_frames=1, rectify=True, normalize=True):
  239. # smoothdfdt firing-rate estimate: per-neuron smoothed first derivative, optionally rectified and peak-normalized
  241. num_cells, num_frames = data.shape
  243. firing_rates = np.zeros((num_cells, num_frames));
  245. if sigma_frames == 0:
  246. do_smooth=False
  248. if do_smooth:
  249. s_fr = np.ceil(sigma_frames).astype(int)
  250. x = np.linspace(-3*s_fr, 3*s_fr, s_fr*6+1)
  251. gauss_kernel = np.exp(-x**2 / (2 * sigma_frames**2))
  253. for n_cell in range(num_cells):
  254. temp_data = np.diff(data[n_cell,:], prepend=0)
  256. if do_smooth:
  257. temp_data = np.convolve(temp_data, gauss_kernel, mode='same')
  259. if rectify:
  260. temp_data = np.maximum(temp_data, 0)
  262. if normalize:
  263. temp_data = temp_data - np.mean(temp_data)
  264. temp_data = temp_data/np.max(temp_data)
  266. firing_rates[n_cell,:] = temp_data;
  268. return firing_rates
  270. def gauss_smooth(firing_rates, sigma_frames=1):
  271. # gaussian smoothing of each neuron's trace (sigma in frames; 0 = no smoothing)
  273. if sigma_frames:
  274. num_cells, num_frames = firing_rates.shape
  276. s_fr = np.ceil(sigma_frames).astype(int)
  277. x = np.linspace(-3*s_fr, 3*s_fr, s_fr*6+1)
  278. gauss_kernel = np.exp(-x**2 / (2 * sigma_frames**2))
  280. firing_rates_sm = np.zeros((num_cells, num_frames));
  281. for n_cell in range(num_cells):
  282. firing_rates_sm[n_cell,:] = np.convolve(firing_rates[n_cell,:], gauss_kernel, mode='same')
  283. else:
  284. firing_rates_sm = firing_rates
  286. return firing_rates_sm
  288. def normalize(rates):
  289. # assumes neurons x time
  291. rates_n = rates - np.min(rates, axis=1)[:, None]
  293. # ignoring cells that are always zero
  294. max_rates = np.max(rates_n, axis=1)
  295. has_max = max_rates > 0
  297. rates_n = rates_n[has_max,:] / max_rates[has_max][:,None]
  299. return rates_n
  302. def get_frames(trial_win = [-0.05, .95], frame_rate = 30):
  303. # anchor at 0
  304. frame_start = np.ceil(trial_win[0] * frame_rate)
  305. frame_end = np.ceil(trial_win[1] * frame_rate)
  306. trial_frames = [int(frame_start), int(frame_end)]
  307. plot_t = np.round(np.arange(frame_start/frame_rate, frame_end/frame_rate, 1/frame_rate), decimals=4)
  308. return trial_frames, plot_t
  310. def get_stim_trig_resp(firing_rates, stim_times, trial_frames = [-29, 85]):
  311. # input: cells x time
  313. num_cells, T = firing_rates.shape
  314. num_trials = len(stim_times)
  316. win_size = trial_frames[1] - trial_frames[0]
  318. stim_trig_resp = np.zeros((num_cells, win_size, num_trials))
  320. for n_tr in range(num_trials):
  321. cur_frame = round(stim_times[n_tr]-1) # correct for matlab to python
  323. raw_start = cur_frame + trial_frames[0]
  324. raw_end = cur_frame + trial_frames[1]
  326. # clamp the window to the recording bounds; trials near the edges stay zero-padded
  327. src_start = np.max([raw_start, 0])
  328. src_end = np.min([raw_end, T])
  330. if src_end > src_start:
  331. dst_start = src_start - raw_start
  332. stim_trig_resp[:, dst_start:dst_start + (src_end - src_start), n_tr] = firing_rates[:, src_start:src_end]
  334. return stim_trig_resp
  336. def save_fig(fig, path='/', name_tag=''):
  337. # save the figure as svg + png (1200 dpi), named from the first-axis title + date + name_tag
  339. plt.rcParams['svg.fonttype'] = 'none'
  340. name1 = fig.axes[0].title.get_text()
  341. now1 = datetime.now()
  343. date_tag = '%d_%d_%d_%dh_%dm' % (now1.year, now1.month, now1.day, now1.hour, now1.minute)
  345. fig.savefig('%s/%s_%s%s.svg' % (path, name1, date_tag, name_tag))
  346. fig.savefig('%s/%s_%s%s.png' % (path, name1, date_tag, name_tag), dpi=1200)
  349. def get_trial_peak(trial_ave, peak_size=3):
  350. # per-neuron peak response: mean over a peak_size-frame window centred on each neuron's argmax; returns (peak_vals, peak_locs)
  351. num_cells, num_bins = trial_ave.shape
  353. pad_left = np.floor((peak_size-1)/2).astype(int)
  354. pad_right = np.ceil((peak_size-1)/2).astype(int)
  356. peak_locs = np.argmax(trial_ave,axis=1)
  358. peak_start = peak_locs - pad_left
  359. peak_end = peak_locs + pad_right + 1
  361. idx_fix_start = peak_start < 0
  362. if sum(idx_fix_start):
  363. peak_start_to_fix = peak_start[idx_fix_start]
  364. peak_start[idx_fix_start] = peak_start[idx_fix_start] - peak_start_to_fix
  365. peak_end[idx_fix_start] = peak_end[idx_fix_start] - peak_start_to_fix
  367. idx_fix_end = peak_end > num_bins
  368. if sum(idx_fix_end):
  369. peak_end_to_fix = peak_end[idx_fix_end]
  370. peak_end[idx_fix_end] = peak_end[idx_fix_end] - peak_end_to_fix + num_bins
  371. peak_start[idx_fix_end] = peak_start[idx_fix_end] - peak_end_to_fix + num_bins
  373. peak_vals = np.zeros(num_cells)
  374. for n_cell in range(num_cells):
  375. peak_vals[n_cell] = np.mean(trial_ave[n_cell,peak_start[n_cell]:peak_end[n_cell]])
  377. return peak_vals, peak_locs
  380. def compute_tuning(stim_trig_resp, trial_types, trials_analyze, plot_t, num_samp=2000, z_thresh = 3, sig_resp_win = [0, 1.5], seed=None):
  381. # find stimulus-responsive cells per trial type by comparing each cell's peak response to a shuffled null
  382. # (z_thresh over num_samp shuffles, within sig_resp_win); returns a (cells x trial types) responsive-cell mask
  384. tt_use_idx = np.sum(trial_types == trials_analyze[:,None],axis=0).astype(bool)
  385. trial_types_use = trial_types[tt_use_idx]
  386. stim_trig_resp_use = stim_trig_resp[:,:,tt_use_idx]
  388. num_cells, _, num_trials = stim_trig_resp_use.shape
  389. num_tt = len(trials_analyze)
  391. # get data
  392. peak_vals = np.full((num_cells, num_tt), np.nan)
  393. peak_locs = np.full((num_cells, num_tt), np.nan)
  394. for n_tt in range(num_tt):
  395. tt1_idx = trial_types_use == trials_analyze[n_tt]
  396. if np.any(tt1_idx):
  397. trial_ave1 = np.mean(stim_trig_resp_use[:,:,tt1_idx], axis=2)
  398. peak_vals[:,n_tt], peak_locs[:,n_tt] = get_trial_peak(trial_ave1, peak_size=3)
  400. # null distribution of peaks from random trial resamples (with replacement)
  401. trials_per_stim = np.zeros(num_tt, dtype=int)
  402. for n_tt in range(num_tt):
  403. trials_per_stim[n_tt] = np.sum(trial_types_use == trials_analyze[n_tt])
  404. trials_per_stim_ave = np.round(np.mean(trials_per_stim)).astype(int)
  406. samp_peak_vals = np.full((num_cells, num_samp), np.nan)
  407. samp_peak_locs = np.full((num_cells, num_samp), np.nan)
  408. rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed)
  409. for n_cell in range(num_cells):
  410. random_integers = rng.integers(low=0, high=num_trials, size=(trials_per_stim_ave, num_samp))
  411. samp_trial_ave = np.mean(stim_trig_resp_use[n_cell,:,random_integers], axis=0)
  412. samp_peak_vals[n_cell,:], samp_peak_locs[n_cell,:] = get_trial_peak(samp_trial_ave, peak_size=3)
  414. idx1 = ~np.isnan(peak_locs[0,:])
  415. peak_locs_t = np.full((num_cells, num_tt), np.nan)
  416. peak_locs_t[:,idx1] = plot_t[peak_locs[:,idx1].astype(int)]
  418. peak_in_resp_win = np.logical_and(peak_locs_t >= sig_resp_win[0], peak_locs_t <= sig_resp_win[1])
  420. # responsive if the observed peak exceeds the z-threshold percentile of the null and peaks within the response window
  421. peak_prcntle = norm.cdf(z_thresh)*100
  422. prc_thresh = np.percentile(samp_peak_vals, peak_prcntle, axis=1)
  423. resp_cells_peak = np.zeros((num_cells, num_tt), dtype=bool)
  425. resp_cells_peak[:,idx1] = np.logical_and(peak_vals[:,idx1] > prc_thresh[:,None], peak_in_resp_win[:,idx1])
  427. return resp_cells_peak
  430. def compute_correlation(stim_trig_resp, trial_types, trials_analyze, resp_cells=None, min_resp_cells=5, subtract_mean=False, add_noise_sigma=1e-5, metric='correlation', cell_select='resp_marg', drop_zero_cells=True, seed=None):
  431. # cell_select: 'resp_marg' (union of responsive cells; = MATLAB 'Resp marg'),
  432. # 'resp_split' (per-frequency responsive cells; = MATLAB 'Resp split'), or 'all'
  433. # drop_zero_cells: True drops cells with non-positive mean across trials (Python-only; MATLAB keeps all)
  434. rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed)
  435. corr_vals = np.full((len(trials_analyze)), np.nan)
  437. if subtract_mean:
  438. stim_trig_resp = stim_trig_resp - np.mean(stim_trig_resp)
  440. if add_noise_sigma: # add some uncorrelated noise to add stability
  441. stim_trig_resp = stim_trig_resp + rng.normal(0, add_noise_sigma, size=stim_trig_resp.shape)
  443. if resp_cells is not None:
  444. resp_marg = np.sum(resp_cells, axis=1).astype(bool)
  445. for n_tn in range(len(trials_analyze)):
  447. if resp_cells is not None:
  448. trig = sum(resp_cells[:,n_tn]) > min_resp_cells
  449. else:
  450. trig = 1
  452. if trig:
  454. tn1 = trials_analyze[n_tn]
  455. tr_idx = trial_types == tn1
  457. stim_trig_resp2 = stim_trig_resp[:,:,tr_idx]
  458. # cell selection mode
  459. if resp_cells is None or cell_select == 'all':
  460. cell_mask = np.ones(stim_trig_resp2.shape[0], dtype=bool)
  461. elif cell_select == 'resp_split':
  462. cell_mask = resp_cells[:,n_tn].astype(bool) # per-frequency responsive cells
  463. else: # 'resp_marg'
  464. cell_mask = resp_marg # union of responsive cells
  465. stim_trig_resp3 = stim_trig_resp2[cell_mask,:,:]
  466. stim_trig_resp4 = np.mean(stim_trig_resp3, axis=1)
  468. if drop_zero_cells:
  469. act_idx = np.mean(stim_trig_resp4, axis=1) > 0
  470. stim_trig_resp5 = stim_trig_resp4[act_idx,:]
  471. else:
  472. stim_trig_resp5 = stim_trig_resp4
  474. distances = squareform(pdist(stim_trig_resp5.T, metric=metric)) # cosine, correlation
  475. # trial-to-trial similarity = 1 - distance; average over the unique trial pairs (lower triangle)
  476. SI = 1 - distances
  477. SI2 = np.tril(SI, k=-1)
  478. SI2_vals = SI2[SI2.astype(bool)]
  479. corr_vals[n_tn] = np.mean(SI2_vals) if len(SI2_vals) else np.nan
  481. # if 0:
  482. # if tn1==4:
  483. # plt.figure()
  484. # plt.imshow(SI)
  485. # plt.title('isi = ' + str(data_echo[n_fl]['isi']))
  487. return corr_vals
  489. def compute_correlation_mat(stim_trig_resp, subtract_mean=False, add_noise_sigma=1e-5, metric='correlation', seed=None):
  490. # trial-by-trial similarity matrix (1 - pairwise distance) of each trial's mean stimulus response
  491. rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed)
  492. if subtract_mean:
  493. stim_trig_resp = stim_trig_resp - np.mean(stim_trig_resp)
  495. if add_noise_sigma: # add some uncorrelated noise to add stability
  496. stim_trig_resp = stim_trig_resp + rng.normal(0, add_noise_sigma, size=stim_trig_resp.shape)
  498. stim_trig_resp2 = np.mean(stim_trig_resp, axis=1)
  500. distances = squareform(pdist(stim_trig_resp2.T, metric=metric)) # cosine, correlation
  501. SI = 1 - distances
  503. return SI
  505. #%%
  506. def get_network_distance(firing_rates):
  507. # euclidean distance of the population vector from its time-averaged baseline at each frame: d(t) = ||r(t) - mean_t r||
  508. base_pop_vec = np.mean(firing_rates, axis=1)
  509. tr_dist = cdist(np.reshape(base_pop_vec, (1,len(base_pop_vec))), firing_rates.T, 'euclidean')[0]
  511. return tr_dist
  513. def get_trace_tau(trace, sm_bin = 0):
  514. # intrinsic timescale of a trace = first autocorrelation lag (in frames) where it drops below 0.5; returns (tau, autocorr)
  516. #sm_bin = 10#round(1/params['dt'])*50;
  517. #trial_len = out_temp_all.shape[1]
  520. # z-score the trace so the zero-lag autocorrelation equals 1 (half-max crossing = 0.5)
  521. tracen = trace - np.mean(trace)
  522. trace_std = np.std(tracen)
  523. if trace_std == 0:
  524. return np.nan, np.full(len(trace), np.nan)
  525. tracen = tracen/trace_std
  527. corr1 = correlate(tracen, tracen)/len(tracen)
  529. #lags = correlation_lags(len(tracen), len(tracen))
  531. if sm_bin:
  532. kernel = np.ones(sm_bin)/sm_bin
  533. corr1_sm = np.convolve(corr1, kernel, mode='same')
  535. corr1_smn = corr1_sm - np.mean(corr1_sm)
  536. corr1_smn = corr1_smn/np.max(corr1_smn)
  537. else:
  538. corr1_smn = corr1
  540. corr1_smn2 = corr1_smn[len(trace)-1:] # positive lags only (zero lag is at index len-1)
  542. # plt.figure(); plt.plot(corr1)
  544. below = np.where(corr1_smn2 < 0.5)[0]
  545. tau_corr = below[0] if len(below) else np.nan
  547. # x = np.arange(corr_len)+1
  548. # y = corr1[num_trials2*num_run:num_trials2*num_run+corr_len]
  550. # yn = y - np.min(y)+0.01
  551. # yn = yn/np.max(yn)
  553. # fit = np.polyfit(x, np.log(yn), 1)
  555. # y_fit = np.exp(x*fit[0]+fit[1])
  557. # tau_corr = np.log(1/2)/fit[0]*params['dt']
  559. # x = np.random.rand(1000)
  560. # corrx = correlate(x, x)
  561. # plt.figure(); plt.plot(corrx)
  563. return tau_corr, corr1_smn2
  566. #%%
  568. def load_rnn_test(data_dir, fname_data, fname_params, max_net_load = 999, limit_network_types=[], flatten_runs = False, max_trial_types = 10, max_trials = 500, cut_zero_trials = False, num_initial_trials_skip = 0, seed=None):
  569. rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed)
  570. # input is the spectrogram inputs
  571. # target is the index when oddball trial happens
  572. # loaded rates shape is (time, run, neurons) - inside converted to
  573. # output rates are converted to (runs, neurons, time)
  574. # flatten runs makes output rates 2D (neurons, time)
  575. # max trial types and max trials only apply during flattened runs
  576. # returns a list of dicts (one per network), matching the calcium-data format;
  577. # main field 'firing_rates' is (neurons, time) if flatten_runs else (runs, neurons, time)
  580. data_path = os.path.join(data_dir, fname_data)
  581. params_path = os.path.join(data_dir, fname_params)
  583. for p in (data_path, params_path):
  584. if not os.path.isfile(p):
  585. raise FileNotFoundError('RNN test file not found: %s' % p)
  587. if '_ob_data' in fname_data:
  588. print('Note: the oddball (_ob_data) RNN file is large (tens of GB) and is loaded fully into memory '
  589. 'before subsetting, so max_net_load does not reduce the peak memory. Ensure sufficient RAM, '
  590. 'or use the smaller control file (_cont_data).')
  592. test_data_load = np.load(data_path, allow_pickle=True).item()
  593. data_key = list(test_data_load.keys())[0]
  594. test_data_load_all = test_data_load[data_key]
  596. # for n1 in range(len(test_data_load_all)):
  597. # for n2 in range(len(test_data_load_all[n1])):
  598. # del test_data_load_all[n1][n2]['input']
  599. # del test_data_load_all[n1][n2]['target']
  600. # del test_data_load_all[n1][n2]['output']
  601. # test_data_load_all[0][0].keys()
  602. # test_data_load_all[0][0]['rates'].shape
  603. # np.save(data_dir + 'RNN_test_data_2024_5_24_9h_42m_ob_data2.npy', test_data_load)
  605. deets_load = np.load(params_path, allow_pickle=True).item()
  606. params_all = deets_load['params_all']
  607. params_test = deets_load['params_test']
  608. # deets_load['ob_data'].keys()
  609. # del deets_load['ob_data']['input_oddball']
  610. # del deets_load['ob_data']['target_oddball_freq']
  611. # del deets_load['cont_data']['input_control']
  612. # del deets_load['cont_data']['target_control']
  613. # np.save(data_dir + 'RNN_test_data_2024_5_24_9h_42m2_params.npy', deets_load)
  615. data_all = []
  617. if len(limit_network_types) == 0:
  618. limit_network_types = list(np.unique(deets_load['rnn_leg']))
  619. use_net_type_lim = True
  622. if 'ob' in data_key:
  623. is_ob = True
  624. else:
  625. is_ob = False
  627. # oddball tained, cont trained, and untrained
  628. for n_net in range(len(test_data_load_all)):
  630. num_rnn = np.min([len(test_data_load_all[n_net]), max_net_load])
  632. if use_net_type_lim:
  633. if deets_load['rnn_leg'][n_net] in limit_network_types:
  635. for n_rnn in range(num_rnn):
  637. net1 = test_data_load_all[n_net][n_rnn]
  639. T, num_runs, num_neurons = net1['rates'].shape
  641. # get stim times
  642. stim_times_runs = []
  643. trial_types_runs = []
  644. if is_ob:
  645. trial_types_red_dd_runs = []
  646. for n_run in range(num_runs):
  648. # deets_load['ob_data'].keys()
  649. # deets_load['ob_data']['trials_oddball_freq']
  650. if 1:
  651. stim_on_trace = 1 - deets_load['ob_data']['target_oddball_ctx3'][:,n_run,0]
  652. else:
  653. stim_trace = np.max(net1['input'][:,n_run], axis=1)
  654. stim_trace_n = stim_trace - np.percentile(stim_trace, 20)
  655. stim_trace_n = stim_trace_n / np.max(stim_trace_n)
  656. stim_on_trace = (stim_trace_n > 0.5).astype(int)
  658. stim_onset_trace = np.diff(stim_on_trace, prepend=[0]) > 0
  659. stim_times = np.where(stim_onset_trace)[0]
  660. stim_times_runs.append(stim_times)
  662. if is_ob:
  663. ob_deets = deets_load['ob_data']
  664. trial_idx = ob_deets['trials_oddball_ctx3'][:,n_run] > 0
  665. trials_types_red_dd = ob_deets['trials_oddball_ctx3'][:,n_run][trial_idx] - 1
  666. trials_types = ob_deets['trials_oddball_freq'][:, n_run][trial_idx]
  667. red_dd_seq = ob_deets['red_dd_seq']
  669. trial_types_red_dd_runs.append(trials_types_red_dd)
  671. else:
  672. cont_deets = deets_load['cont_data']
  673. trial_idx = cont_deets['trials_control_freq'][:,n_run] > 0
  674. trials_types = cont_deets['trials_control_freq'][:, n_run][trial_idx]
  676. trial_types_runs.append(trials_types)
  678. stim_times_runs2 = np.vstack(stim_times_runs)
  679. trial_types_runs2 = np.vstack(trial_types_runs)
  680. if is_ob:
  681. trial_types_red_dd_runs2 = np.vstack(trial_types_red_dd_runs)
  683. if cut_zero_trials or num_initial_trials_skip > 0:
  685. trial_len = round((params_test['stim_duration'] + params_test['isi_duration']) / params_test['dt'])
  686. if T % trial_len != 0:
  687. raise ValueError('RNN trace length T=%d is not divisible by trial_len=%d; cannot reshape into whole trials (check stim_duration/isi_duration/dt).' % (T, trial_len))
  688. rates3d = net1['rates']
  689. rates4d = np.reshape(rates3d, (round(T/trial_len), trial_len, num_runs, num_neurons), order='C')
  691. if cut_zero_trials:
  692. num_skip = params_test['num_prepend_zeros'] + num_initial_trials_skip
  693. else:
  694. num_skip = num_initial_trials_skip
  696. if num_skip >= rates4d.shape[0]:
  697. raise ValueError('num_initial_trials_skip (+ prepended zeros) = %d exceeds the %d available trials; reduce num_initial_trials_skip.' % (num_skip, rates4d.shape[0]))
  699. rates4d_cut = rates4d[num_skip:,:,:]
  700. rates2 = np.reshape(rates4d_cut, ((round(T/trial_len) - num_skip) * trial_len, num_runs, num_neurons), order='C')
  702. stim_times_runs2 = stim_times_runs2[:,num_skip:]
  703. trial_types_runs2 = trial_types_runs2[:,num_skip:]
  705. else:
  706. rates2 = net1['rates']
  708. T, num_runs, num_neurons = rates2.shape
  710. if flatten_runs:
  711. stim_times_runs_flat = (stim_times_runs2 + np.arange(num_runs)[:,None] * T).flatten()
  712. trial_types_runs_flat = trial_types_runs2.flatten()
  713. if is_ob:
  714. trial_types_red_dd_runs2 = trial_types_red_dd_runs2.flatten()
  715. rates = np.reshape(rates2, (T * num_runs, num_neurons), order='F').T
  717. trials_uq = np.unique(trial_types_runs_flat)
  718. num_tt = len(trials_uq)
  720. num_stim_all = len(stim_times_runs_flat)
  722. if max_trial_types < num_tt:
  723. step = num_tt/10
  724. trial_sel = np.arange(step/2, num_tt, step=step, dtype=int)
  725. sel_trial_idx = np.sum(trial_types_runs_flat[:,None] == trials_uq[trial_sel][None,:], axis=1).astype(bool)
  726. sel_trials = np.where(sel_trial_idx)[0]
  727. else:
  728. sel_trials = np.arange(num_stim_all)
  730. num_stim_all2 = len(sel_trials)
  732. # limit trials to max number
  733. if max_trials < num_stim_all:
  734. trial_idx = rng.choice(sel_trials, size=np.min([max_trials, num_stim_all2]), replace=False)
  735. trial_idx.sort()
  737. else:
  738. trial_idx = sel_trials
  740. stim_times_runs3 = stim_times_runs_flat[trial_idx]
  741. trial_types_runs3 = trial_types_runs_flat[trial_idx]
  742. if is_ob:
  743. trial_types_red_dd_runs3 = trial_types_red_dd_runs2[trial_idx]
  745. else:
  746. stim_times_runs3 = stim_times_runs2
  747. trial_types_runs3 = trial_types_runs2
  748. if is_ob:
  749. trial_types_red_dd_runs3 = trial_types_red_dd_runs2
  750. rates = np.transpose(rates2, (1, 2, 0))
  753. data_slice = {}
  754. data_slice['training'] = deets_load['rnn_leg'][n_net]
  755. data_slice['test_data_key'] = data_key
  756. data_slice['firing_rates'] = rates
  757. data_slice['stim_times'] = stim_times_runs3
  758. data_slice['trial_types'] = trial_types_runs3
  759. if is_ob:
  760. data_slice['trial_types_red_dd'] = trial_types_red_dd_runs3
  761. data_slice['red_dd_seq'] = red_dd_seq
  762. data_slice['num_runs'] = num_runs
  763. data_slice['volume_period'] = params_test['dt']*1000 # for frame rate equivalent
  764. data_slice['params_train'] = params_all[n_net][n_rnn]
  765. data_slice['params_test'] = params_test
  767. data_all.append(data_slice)
  769. del net1
  771. if len(data_all):
  772. trainings = [d['training'] for d in data_all]
  773. uq, cnt = np.unique(trainings, return_counts=True)
  774. print('Loaded %d RNN networks: %s' % (len(data_all), dict(zip(list(uq), [int(c) for c in cnt]))))
  775. else:
  776. print('Warning: no RNN networks loaded (file %s, limit_network_types=%s)' % (fname_data, limit_network_types))
  778. return data_all
  780. #%%
  782. def get_network_intrinsic_timescales(firing_rates, frame_rate):
  783. # per recording (or RNN run): network tau from the baseline-distance autocorrelation and per-neuron tau; both in seconds
  785. if len(firing_rates.shape) == 2:
  786. rates3d = firing_rates[None,:,:]
  787. else:
  788. rates3d = firing_rates
  790. num_runs, num_neurons, _ = rates3d.shape
  792. tau_net = np.zeros(num_runs)
  793. tau_cell = np.full(shape=(num_runs, num_neurons), fill_value=np.nan)
  795. for n_run in range(num_runs):
  797. tr_dist = get_network_distance(rates3d[n_run,:,:])
  799. tau_net1, _ = get_trace_tau(tr_dist, sm_bin = 0)
  800. tau_net[n_run] = tau_net1/frame_rate
  802. for n_nr in range(num_neurons):
  804. neur = rates3d[n_run,n_nr,:]
  805. if np.sum(neur) > 0.1:
  806. tau_neur1, _ = get_trace_tau(neur, sm_bin = 0)
  807. tau_cell[n_run, n_nr] = tau_neur1/frame_rate
  809. return tau_net, tau_cell
  811. #%%
  814. def plot_fig_raster(firing_rates_ob, firing_rates_rnn, frame_rate=None, frame_rate_rnn=None):
  815. # ---- plot example raster ----
  816. fig, ax = plt.subplots(1,2, figsize=(12,5), layout='constrained')
  817. fig.set_constrained_layout_pads(w_pad=0.2, h_pad=0.3, hspace=0, wspace=0)
  818. fig.text(0.015, .89, 'A', fontsize=18)
  819. fig.text(0.515, .89, 'B', fontsize=18)
  821. num_cells, num_t = firing_rates_ob.shape
  823. if frame_rate is not None:
  824. x_end = num_t/frame_rate
  825. x_lab = 'Time (sec)'
  826. else:
  827. x_end = num_t
  828. x_lab = 'Frames'
  830. ax[0].imshow(firing_rates_ob,
  831. aspect='auto',
  832. cmap='gist_yarg',
  833. vmin=0,
  834. vmax=.5,
  835. extent=[0, x_end, 1, num_cells],
  836. interpolation='none')
  837. ax[0].set_title('CaIm control data')
  838. ax[0].set_ylabel('Neurons')
  839. ax[0].set_xlabel(x_lab)
  841. num_cells_rnn, num_t_rnn = firing_rates_rnn.shape
  843. if frame_rate_rnn is not None:
  844. x_end = num_t_rnn/frame_rate_rnn
  845. x_lab = 'Time (sec)'
  846. else:
  847. x_end = num_t_rnn
  848. x_lab = 'Frames'
  850. ax[1].imshow(normalize(firing_rates_rnn),
  851. aspect='auto',
  852. cmap='gist_yarg',
  853. vmin=0,
  854. vmax=1,
  855. extent=[0, x_end, 1, num_cells_rnn],
  856. interpolation='none')
  857. ax[1].set_title('RNN control test data')
  858. ax[1].set_ylabel('Neurons')
  859. ax[1].set_xlabel(x_lab)
  861. return fig
  863. def plot_fig_diag_decoder(dec_data_all, dec_data_all_rnn, training_type, plot_t=None, plot_t_rnn=None, add_sig=True):
  864. # ---- Plotting diagonal decoder results ----
  865. # add_sig=True overlays binwise data-vs-shuffle significance (see plot_diag_binwise_dec)
  866. fig_diag, ax = plt.subplots(1,2, figsize=(12,5), layout='constrained')
  867. fig_diag.set_constrained_layout_pads(w_pad=0.2, h_pad=0.1, hspace=0, wspace=0)
  868. fig_diag.text(0.015, .87, 'A', fontsize=18)
  869. fig_diag.text(0.515, .87, 'B', fontsize=18)
  870. fig_diag.suptitle('Binwise decoder')
  872. plot_diag_binwise_dec(
  873. dec_data_all,
  874. plot_t=plot_t,
  875. plot_legend=('Data', 'Shuff'),
  876. plot_start=-1, # plot window start
  877. plot_end=3, # plot window end
  878. axis = ax[0],
  879. title_tag='CaIm data',
  880. add_sig=add_sig)
  882. plot_diag_binwise_dec(
  883. np.array(dec_data_all_rnn)[training_type == 'ob trained'],
  884. plot_t=plot_t_rnn,
  885. plot_start=-1,
  886. plot_end=8,
  887. axis=ax[1],
  888. title_tag='RNN test data',
  889. colors = ['limegreen', 'darkgreen'],
  890. add_sig=add_sig)
  892. plot_diag_binwise_dec(
  893. np.array(dec_data_all_rnn)[training_type == 'freq trained'],
  894. plot_t=plot_t_rnn,
  895. plot_start=-1,
  896. plot_end=8,
  897. axis=ax[1],
  898. colors = ['orange', 'saddlebrown'],
  899. add_sig=add_sig)
  901. # legend from the mean-trace lines only (exclude the '.-' binwise-significance lines)
  902. main_lines = [l for l in ax[1].get_lines() if l.get_marker() in ('', 'None', None)]
  903. ax[1].legend(main_lines, ['Oddball trained', 'Oddball trained shuff', 'Freq trained', 'Freq trained shuff'])
  905. return fig_diag
  907. def plot_fig_full_decoder_caim(dec_data_full, plot_t, clim=[0, 0.5]):
  909. # ---- Plotting full decoder results CaIm data ----
  910. fig_full, ax = plt.subplots(1, 3, gridspec_kw={'width_ratios': [20, 20, 1]}, figsize=(12, 5.7), layout='constrained')
  911. fig_full.set_constrained_layout_pads(w_pad=0.01, h_pad=0.1, hspace=0, wspace=0)
  912. fig_full.text(0.01, .88, 'A', fontsize=16)
  913. fig_full.text(0.47, .88, 'B', fontsize=16)
  914. fig_full.suptitle('CaIm data full space decoder', fontsize=14)
  916. plot_full_binwise_dec(
  917. dec_data_full,
  918. plot_t=plot_t,
  919. plot_legend=('Data', 'Shuff'),
  920. plot_start=-1,
  921. plot_end=2,
  922. clim=clim,
  923. axis=ax,
  924. title_tag='')
  926. return fig_full
  928. def plot_fig_full_decoder_RNN(dec_data_full_ob, dec_data_full_freq, plot_t):
  930. # ---- Plotting full decoder results RNN data ----
  931. fig_full, ax = plt.subplots(2, 3, gridspec_kw={'width_ratios': [20, 20, 1]}, figsize=(12,11), layout='constrained')
  932. fig_full.set_constrained_layout_pads(w_pad=0.05, h_pad=0.1, hspace=0, wspace=0)
  933. fig_full.text(0.005, .94, 'A', fontsize=16)
  934. fig_full.text(0.465, .94, 'B', fontsize=16)
  935. fig_full.text(0.005, .46, 'C', fontsize=16)
  936. fig_full.text(0.465, .46, 'D', fontsize=16)
  937. fig_full.suptitle('RNN neurons full space decoder', fontsize=14)
  939. plot_full_binwise_dec(
  940. dec_data_full_ob,
  941. plot_t=plot_t,
  942. plot_legend=('data', 'shuff'),
  943. plot_start=-1,
  944. plot_end=10,
  945. clim=[0, 1],
  946. axis=ax[0,:],
  947. title_tag='RNN Oddball trained')
  949. plot_full_binwise_dec(
  950. dec_data_full_freq,
  951. plot_t=plot_t,
  952. plot_legend=('data', 'shuff'),
  953. plot_start=-1,
  954. plot_end=10,
  955. clim=[0, 1],
  956. axis=ax[1,:],
  957. title_tag='RNN Control freq trained')
  959. return fig_full
  961. def plot_fig_full_decoder_comb(dec_data_full, dec_data_full_ob, dec_data_full_freq, plot_t, plot_t_rnn):
  962. # ---- plotting full decoder results ----
  963. fig_full, ax = plt.subplots(3, 3, gridspec_kw={'width_ratios': [10, 10, 1]}, figsize=(15, 12))
  964. fig_full.text(0.09, .88, 'A', fontsize=16)
  965. fig_full.text(0.51, .88, 'B', fontsize=16)
  967. plot_full_binwise_dec(
  968. dec_data_full,
  969. plot_t=plot_t,
  970. plot_legend=('Data', 'Shuff'),
  971. plot_start=-1,
  972. plot_end=2,
  973. clim=[0, 0.5],
  974. axis=ax[0,:],
  975. title_tag='CaIm data')
  978. plot_full_binwise_dec(
  979. dec_data_full_ob,
  980. plot_t=plot_t_rnn,
  981. plot_legend=('Oddball trained', 'Shuff'),
  982. plot_start=-1,
  983. plot_end=10,
  984. clim=[0, 1],
  985. axis=ax[1,:],
  986. title_tag='RNN Ob trained')
  989. plot_full_binwise_dec(
  990. dec_data_full_freq,
  991. plot_t=plot_t_rnn,
  992. plot_legend=('Freq trained', 'Shuff'),
  993. plot_start=-1,
  994. plot_end=10,
  995. clim=[0, 1],
  996. axis=ax[2,:],
  997. title_tag='RNN Freq trained')
  999. return fig_full
  1001. #%%
  1003. def plot_fig_isi_corr_trials(corr_vals, isi_list, colormap='jet', metric_tag = None):
  1004. # trial-to-trial correlation vs ISI: panel A = mean over frequencies, panel B = per-frequency lines; returns (fig, ax, groups)
  1006. fig, ax = plt.subplots(1,2, figsize=(12,5))
  1007. fig.text(0.07, .88, 'A', fontsize=16)
  1008. fig.text(0.51, .88, 'B', fontsize=16)
  1010. idx_uq = np.unique(isi_list)
  1011. num_trials = corr_vals.shape[1]
  1013. col1 = plt.colormaps[colormap](np.linspace(0, 1, num_trials))
  1015. corr_tn_all = np.zeros((num_trials, len(idx_uq)))
  1016. for n_tn in range(num_trials):
  1017. corr_tn = np.full(len(idx_uq), np.nan)
  1018. for n_isi in range(len(idx_uq)):
  1019. idx1 = (idx_uq[n_isi] == np.array(isi_list)).flatten()
  1020. if np.sum(~np.isnan(corr_vals[idx1,n_tn])):
  1021. corr_tn[n_isi] = np.nanmean(corr_vals[idx1,n_tn])
  1022. corr_tn_all[n_tn, n_isi] = np.nanmean(corr_vals[idx1,n_tn])
  1023. if np.sum(~np.isnan(corr_tn)):
  1024. ax[1].plot(idx_uq, corr_tn, '-o', color=col1[n_tn])
  1026. if metric_tag is not None:
  1027. ax[1].set_ylabel(metric_tag)
  1028. else:
  1029. ax[1].set_ylabel('Correlation')
  1030. ax[1].set_xlabel('ISI duration (sec)')
  1031. freqs = np.logspace(np.log10(2), np.log10(76.9), num_trials) # kHz, log-spaced 2 -> 76.9 (x1.5/step)
  1032. ax[1].legend(['%g' % round(f, 1) for f in freqs], loc='upper right', title='Freq (kHz)')
  1033. ax[1].set_title('Individual freqs.')
  1035. mean1 = np.nanmean(corr_tn_all, axis=0)
  1036. sem1 = np.nanstd(corr_tn_all, axis=0)/np.sqrt(np.sum(~np.isnan(corr_tn_all), axis=0) - 1)
  1038. ax[0].errorbar(idx_uq, mean1, sem1, fmt='-o', color='k', capsize=4)
  1040. ax[0].set_ylabel('Correlation')
  1041. ax[0].set_xlabel('ISI duration (sec)')
  1042. ax[0].set_title('Average over freqs.')
  1044. # per-ISI groups for stat_compare; each entry is (data_array, x-position = ISI value).
  1045. # here we include every individual (dataset x frequency) correlation value at each ISI
  1046. # (not the per-dataset mean), so all the individual points enter the test. Pass ax[0], e.g.
  1047. # sd.stat_compare(ax[0], groups, 'ISI 0.5', 'ISI 1', test='mannwhitney', alternative='two-sided')
  1048. groups = {}
  1049. for isi in idx_uq:
  1050. sel = (np.array(isi_list) == isi).flatten()
  1051. d = corr_vals[sel, :].flatten() # all datasets x frequencies at this ISI
  1052. groups['ISI %g' % isi] = (d[np.isfinite(d)], float(isi))
  1054. return fig, ax, groups
  1056. def plot_fig_isi_corr_trials2(corr_vals, isi_list, colormap='jet', metric_tag=None):
  1057. # nicer version: panel A = individual (dataset x freq) points + mean line + shaded SEM;
  1058. # panel B = per-frequency lines. Returns (fig, ax, groups) like plot_fig_isi_corr_trials.
  1059. fig, ax = plt.subplots(1, 2, figsize=(12, 5))
  1060. fig.text(0.07, .88, 'A', fontsize=16)
  1061. fig.text(0.51, .88, 'B', fontsize=16)
  1063. idx_uq = np.unique(isi_list)
  1064. num_trials = corr_vals.shape[1]
  1065. isi_arr = np.array(isi_list).flatten()
  1066. rng = np.random.default_rng(0) # reproducible jitter
  1068. # ---- panel B: per-frequency correlation vs ISI ----
  1069. col1 = plt.colormaps[colormap](np.linspace(0, 1, num_trials))
  1070. corr_tn_all = np.full((num_trials, len(idx_uq)), np.nan)
  1071. for n_tn in range(num_trials):
  1072. for n_isi in range(len(idx_uq)):
  1073. vals = corr_vals[idx_uq[n_isi] == isi_arr, n_tn]
  1074. if np.sum(~np.isnan(vals)):
  1075. corr_tn_all[n_tn, n_isi] = np.nanmean(vals)
  1076. ax[1].plot(idx_uq, corr_tn_all[n_tn, :], '-o', color=col1[n_tn], markersize=4, linewidth=1.2)
  1077. ax[1].set_xlabel('ISI duration (sec)')
  1078. ax[1].set_ylabel(metric_tag if metric_tag is not None else 'Correlation')
  1079. ax[1].set_title('Individual freqs.')
  1080. ax[1].legend([str(n + 1) for n in range(num_trials)], loc='upper right', fontsize=8, ncol=2, frameon=False)
  1082. # ---- panel A: individual points + mean line + shaded SEM ----
  1083. means = np.full(len(idx_uq), np.nan)
  1084. sems = np.full(len(idx_uq), np.nan)
  1085. for n_isi in range(len(idx_uq)):
  1086. pts = corr_vals[idx_uq[n_isi] == isi_arr, :].flatten()
  1087. pts = pts[np.isfinite(pts)]
  1088. if len(pts) == 0:
  1089. continue
  1090. jit = (rng.random(len(pts)) - 0.5) * 0.08
  1091. ax[0].plot(np.full(len(pts), idx_uq[n_isi]) + jit, pts, '.', color='0.7',
  1092. markersize=4, alpha=0.5, zorder=1)
  1093. means[n_isi] = np.mean(pts)
  1094. sems[n_isi] = np.std(pts) / np.sqrt(len(pts) - 1) if len(pts) > 1 else np.nan
  1096. ax[0].fill_between(idx_uq, means - sems, means + sems, color='steelblue', alpha=0.3, zorder=2)
  1097. ax[0].plot(idx_uq, means, '-o', color='steelblue', markersize=6, linewidth=2, zorder=3)
  1098. ax[0].set_xlabel('ISI duration (sec)')
  1099. ax[0].set_ylabel(metric_tag if metric_tag is not None else 'Correlation')
  1100. ax[0].set_title('Average over freqs.')
  1102. for a in ax:
  1103. a.spines['top'].set_visible(False)
  1104. a.spines['right'].set_visible(False)
  1106. # ---- groups (all individual dataset x freq points per ISI) for stat_compare ----
  1107. groups = {}
  1108. for isi in idx_uq:
  1109. d = corr_vals[isi_arr == isi, :].flatten()
  1110. groups['ISI %g' % isi] = (d[np.isfinite(d)], float(isi))
  1112. return fig, ax, groups
  1114. def plot_fig_SI_mat(SI_list, isi_uq, title_tag = ''):
  1115. # plot the per-ISI trial-by-trial similarity matrices side by side with a shared colorbar
  1116. fig, ax = plt.subplots(1,len(isi_uq)+1, gridspec_kw={'width_ratios': list(np.ones(len(isi_uq))*10) + [1]}, figsize=(12, 2.8))
  1117. for n_isi in range(len(SI_list)):
  1118. ax1 = ax.flatten()[n_isi]
  1119. im = ax1.imshow(SI_list[n_isi], vmin=0, vmax=.7)
  1121. ax1.set_title('isi %s sec' % (isi_uq[n_isi]))
  1122. if title_tag is not None:
  1123. fig.suptitle(title_tag)
  1124. ax.flatten()[0].set_ylabel('trials')
  1125. ax.flatten()[0].set_xlabel('trials')
  1126. fig.colorbar(im, cax=ax.flatten()[-1])
  1127. ax.flatten()[-1].set_ylabel('cosine similarity')
  1129. return fig
  1132. def plot_fig_tau_networks_comb3(tau_ob_net_all, tau_rnn_net_all, training_type, tau_ob_cell_all, tau_rnn_cell_all, do_log=True):
  1133. # legacy two-panel network/neuron tau figure (superseded by plot_fig_tau_networks_comb)
  1135. fig, ax, = plt.subplots(1, 2, sharey=True, figsize=(12,5))
  1137. fig.text(0.075, .88, 'A', fontsize=16)
  1138. fig.text(0.495, .88, 'B', fontsize=16)
  1140. data_all = [np.array(tau_ob_net_all).flatten()] + tau_rnn_net_all
  1141. labels_all = np.array(['Caim data'] + list(training_type))
  1143. plot_int_violin2(data_all,
  1144. net_labels = labels_all,
  1145. axis=ax[0],
  1146. points=1000,
  1147. mean_std=True,
  1148. showmeans=True,
  1149. showmedians=False,
  1150. quantile = [0.05, 0.95],
  1151. colors=['blue', 'green', 'orange', 'gray'],
  1152. do_log=True)
  1154. tau_rnn_cell_all2 = []
  1155. for n_net in range(len(tau_rnn_cell_all)):
  1156. tau_rnn_cell_all2.append(np.nanmean(tau_rnn_cell_all[n_net], axis=0))
  1158. data_cell_all = [np.hstack(tau_ob_cell_all).flatten()] + tau_rnn_cell_all2
  1160. plot_int_violin2(data_cell_all,
  1161. net_labels = labels_all,
  1162. axis=ax[1],
  1163. points=1000,
  1164. mean_std=True,
  1165. showmeans=True,
  1166. showmedians=False,
  1167. quantile = [0.05, 0.95],
  1168. colors=['blue', 'green', 'orange', 'gray'],
  1169. do_log=True)
  1171. ax[1].yaxis.set_tick_params(labelleft=True)
  1172. ax[0].set_title('Network tau')
  1173. ax[1].set_title('Neuron tau')
  1175. return fig
  1177. def get_tau_groups(tau_ob_net_all, tau_rnn_net_all, training_type, tau_ob_cell_all, tau_rnn_cell_all, gap=1.0):
  1178. # build the 8 tau groups (network block + neuron block) shared by the plot and the stats.
  1179. # returns:
  1180. # groups: dict {label: (data_array, x_position)} e.g. 'CaIm net', 'CaIm neuron', ...
  1181. # meta: dict with 'labels', 'x_labels', 'positions', 'colors', 'n_net', 'data'
  1182. tt = np.array(training_type)
  1183. uq_types = list(dict.fromkeys(tt.tolist())) # RNN types in order of appearance
  1184. short = {'ob trained': 'Ob', 'freq trained': 'Freq', 'untrained': 'Untr'}
  1185. sh = lambda l: short.get(l, l)
  1187. # network tau, one group per dataset type
  1188. net_groups = [np.array(tau_ob_net_all).flatten()]
  1189. net_names = ['CaIm']
  1190. for u in uq_types:
  1191. net_groups.append(np.hstack([np.asarray(tau_rnn_net_all[i]).flatten()
  1192. for i in range(len(tau_rnn_net_all)) if tt[i] == u]))
  1193. net_names.append(sh(u))
  1195. # neuron tau (mean over runs per network, then pooled per type)
  1196. tau_rnn_cell_mean = [np.nanmean(tau_rnn_cell_all[n], axis=0) for n in range(len(tau_rnn_cell_all))]
  1197. neuron_groups = [np.hstack(tau_ob_cell_all).flatten()]
  1198. for u in uq_types:
  1199. neuron_groups.append(np.hstack([np.asarray(tau_rnn_cell_mean[i]).flatten()
  1200. for i in range(len(tau_rnn_cell_mean)) if tt[i] == u]))
  1202. data_all = [g[np.isfinite(g)] for g in (net_groups + neuron_groups)]
  1203. labels = [n + ' net' for n in net_names] + [n + ' neuron' for n in net_names]
  1204. x_labels = net_names + net_names
  1205. base_colors = ['blue', 'green', 'orange', 'gray']
  1206. colors = base_colors[:len(net_names)] + base_colors[:len(net_names)]
  1208. n_net = len(net_names)
  1209. positions = list(range(n_net)) + [p + n_net + gap for p in range(n_net)] # e.g. [0,1,2,3, 5,6,7,8]
  1211. groups = {labels[i]: (data_all[i], positions[i]) for i in range(len(labels))}
  1212. meta = {'labels': labels, 'x_labels': x_labels, 'positions': positions,
  1213. 'colors': colors, 'n_net': n_net, 'data': data_all}
  1214. return groups, meta
  1217. def plot_fig_tau_networks_comb2(tau_ob_net_all, tau_rnn_net_all, training_type, tau_ob_cell_all, tau_rnn_cell_all, do_log=True):
  1218. # all 8 groups (network block + neuron block) on one shared-y axis.
  1219. # returns (fig, ax, groups); pass ax + groups to stat_compare() to add significance brackets.
  1220. groups, meta = get_tau_groups(tau_ob_net_all, tau_rnn_net_all, training_type, tau_ob_cell_all, tau_rnn_cell_all)
  1221. data_all = meta['data']
  1222. positions = meta['positions']
  1223. colors_all = meta['colors']
  1224. x_labels = meta['x_labels']
  1225. n_net = meta['n_net']
  1226. num = len(data_all)
  1228. fig, ax1 = plt.subplots(1, 1, figsize=(12, 5))
  1230. parts = ax1.violinplot(data_all, positions=positions, showmeans=False, showextrema=False,
  1231. quantiles=[[0.05, 0.95]] * num, points=1000)
  1232. if 'cquantiles' in parts:
  1233. parts['cquantiles'].set_color('k')
  1234. for i in range(num):
  1235. parts['bodies'][i].set_facecolor(colors_all[i])
  1236. parts['bodies'][i].set_edgecolor(colors_all[i])
  1237. for i in range(num):
  1238. y = data_all[i]
  1239. ax1.plot(positions[i], np.mean(y), '_', color='black', mew=2, markersize=30)
  1240. ax1.errorbar(positions[i], np.mean(y), np.std(y), fmt='o', color='black', markersize=4, linewidth=2, capsize=8)
  1242. ax1.set_xticks(positions)
  1243. ax1.set_xticklabels(x_labels, rotation=0)
  1244. if do_log:
  1245. ax1.set_yscale('log')
  1246. ax1.set_ylabel('Tau (sec)')
  1247. ax1.set_title('Network and neuron intrinsic timescales')
  1249. # group underlines + labels beneath (x in data coords, y in axes fraction)
  1250. trans = ax1.get_xaxis_transform()
  1251. net_pos = positions[:n_net]
  1252. neuron_pos = positions[n_net:]
  1253. for xs, name, col in [(net_pos, 'Networks', 'black'), (neuron_pos, 'Neurons', 'red')]:
  1254. ax1.plot([xs[0] - 0.45, xs[-1] + 0.45], [-0.13, -0.13], transform=trans,
  1255. color=col, lw=4, clip_on=False)
  1256. ax1.text(np.mean(xs), -0.19, name, transform=trans, ha='center', va='top',
  1257. fontsize=13, fontweight='bold', color=col)
  1259. fig.subplots_adjust(bottom=0.26)
  1261. return fig, ax1, groups
  1264. def plot_fig_tau_networks_comb(tau_ob_net_all, tau_rnn_net_all, training_type, tau_ob_cell_all, tau_rnn_cell_all, do_log=True):
  1265. # cleaner version of comb2: violin body + a slim inner boxplot (median + IQR + whiskers),
  1266. # instead of the mean/SD dash + 5/95 quantile lines. Returns (fig, ax, groups) for stat_compare().
  1267. groups, meta = get_tau_groups(tau_ob_net_all, tau_rnn_net_all, training_type, tau_ob_cell_all, tau_rnn_cell_all)
  1268. data_all = meta['data']
  1269. positions = meta['positions']
  1270. colors_all = meta['colors']
  1271. x_labels = meta['x_labels']
  1272. n_net = meta['n_net']
  1273. num = len(data_all)
  1275. fig, ax1 = plt.subplots(1, 1, figsize=(12, 5))
  1277. # violin bodies only (no extrema / quantile lines)
  1278. parts = ax1.violinplot(data_all, positions=positions, showmeans=False, showextrema=False,
  1279. showmedians=False, points=1000)
  1280. for i in range(num):
  1281. parts['bodies'][i].set_facecolor(colors_all[i])
  1282. parts['bodies'][i].set_edgecolor(colors_all[i])
  1283. parts['bodies'][i].set_alpha(0.55)
  1285. # slim inner boxplot: median + IQR + whiskers, no outliers
  1286. ax1.boxplot(data_all, positions=positions, widths=0.12, showfliers=False, patch_artist=True,
  1287. medianprops=dict(color='black', linewidth=1.5),
  1288. boxprops=dict(facecolor='white', edgecolor='black', linewidth=1.0),
  1289. whiskerprops=dict(color='black', linewidth=1.0),
  1290. capprops=dict(color='black', linewidth=1.0))
  1292. ax1.set_xticks(positions)
  1293. ax1.set_xticklabels(x_labels, rotation=0)
  1294. if do_log:
  1295. ax1.set_yscale('log')
  1296. # plain-number tick labels (0.1, 1, 10, ...) instead of 10^0 exponential form
  1297. ax1.yaxis.set_major_formatter(FuncFormatter(lambda v, _: ('%g' % v)))
  1298. ax1.yaxis.set_minor_formatter(NullFormatter())
  1299. ax1.set_ylabel('Tau (sec)')
  1300. ax1.set_title('Network and neuron intrinsic timescales')
  1302. # group underlines + labels beneath (x in data coords, y in axes fraction)
  1303. trans = ax1.get_xaxis_transform()
  1304. net_pos = positions[:n_net]
  1305. neuron_pos = positions[n_net:]
  1306. for xs, name, col in [(net_pos, 'Networks', 'black'), (neuron_pos, 'Neurons', 'red')]:
  1307. ax1.plot([xs[0] - 0.45, xs[-1] + 0.45], [-0.13, -0.13], transform=trans,
  1308. color=col, lw=4, clip_on=False)
  1309. ax1.text(np.mean(xs), -0.19, name, transform=trans, ha='center', va='top',
  1310. fontsize=13, fontweight='bold', color=col)
  1312. fig.subplots_adjust(bottom=0.26)
  1314. return fig, ax1, groups
  1317. def p_to_star(p):
  1318. # p-value to significance stars: *** <0.001, ** <0.01, * <0.05, else n.s.
  1319. if p < 0.001: return '***'
  1320. if p < 0.01: return '**'
  1321. if p < 0.05: return '*'
  1322. return 'n.s.'
  1325. def draw_sig_bracket(ax, x1, x2, p, star=True, color='black', y=None):
  1326. # draw a significance bracket + label between x1 and x2 in the space above the data,
  1327. # extending the y-axis to fit. Stacked calls space evenly (pixel gaps) on log or linear.
  1328. sig = p_to_star(p)
  1329. def _yoff(y_data, pts):
  1330. x0 = ax.get_xlim()[0]
  1331. yd = ax.transData.transform((x0, y_data))[1] + pts * ax.figure.dpi / 72.0
  1332. return ax.transData.inverted().transform((x0, yd))[1]
  1334. step_pts = 15 # spacing between stacked brackets
  1335. margin_pts = 24 # space above the topmost bracket/star and the plot edge
  1336. last = getattr(ax, '_sig_last_y', None)
  1337. if y is not None:
  1338. y_line = y
  1339. elif last is None:
  1340. y_line = _yoff(ax.get_ylim()[1], step_pts)
  1341. else:
  1342. y_line = _yoff(last, step_pts)
  1343. ax._sig_last_y = y_line
  1344. y_tick = _yoff(y_line, -5)
  1345. ax.set_ylim(top=_yoff(y_line, margin_pts))
  1346. ax.plot([x1, x1, x2, x2], [y_tick, y_line, y_line, y_tick], color=color, lw=1.5, clip_on=False)
  1347. # '*' renders high in its box (va='top'); text like 'n.s.'/'p=..' sits low (va='bottom')
  1348. label = sig if star else ('p=%.3g' % p)
  1349. if set(label) <= set('*'):
  1350. va, dy = 'top', 7
  1351. else:
  1352. va, dy = 'bottom', 3
  1353. ax.annotate(label, xy=((x1 + x2) / 2.0, y_line), xytext=(0, dy), textcoords='offset points',
  1354. ha='center', va=va, color=color, fontsize=12, annotation_clip=False)
  1357. def stat_compare(ax, groups, name1, name2, test='mannwhitney', alternative='two-sided',
  1358. stat_file=None, y=None, star=True, color='black'):
  1359. # compare two groups from get_tau_groups(); print result, append to a stats file,
  1360. # and draw a significance bracket on the same panel (ax) between the two groups.
  1361. # test: 'mannwhitney' (unpaired) or 'wilcoxon' (paired; needs equal-length data)
  1362. # alternative: 'two-sided', 'greater', or 'less' (name1 vs name2)
  1363. d1, x1 = groups[name1]
  1364. d2, x2 = groups[name2]
  1365. d1 = np.asarray(d1, dtype=float); d1 = d1[np.isfinite(d1)]
  1366. d2 = np.asarray(d2, dtype=float); d2 = d2[np.isfinite(d2)]
  1368. if test in ('wilcoxon', 'paired'):
  1369. if len(d1) != len(d2):
  1370. raise ValueError('wilcoxon needs paired equal-length data (got %d and %d); '
  1371. 'use test="mannwhitney" or pass paired per-dataset values' % (len(d1), len(d2)))
  1372. stat, p = wilcoxon(d1, d2, alternative=alternative)
  1373. n_str = 'n=%d pairs' % len(d1)
  1374. elif test in ('mannwhitney', 'mannwhitneyu', 'mwu'):
  1375. stat, p = mannwhitneyu(d1, d2, alternative=alternative)
  1376. n_str = 'n1=%d; n2=%d' % (len(d1), len(d2))
  1377. elif test in ('ttest', 'ttest_ind', 't'):
  1378. stat, p = ttest_ind(d1, d2, equal_var=False, alternative=alternative) # Welch's t-test
  1379. n_str = 'n1=%d; n2=%d' % (len(d1), len(d2))
  1380. elif test in ('ttest_rel', 'ttest_paired', 'paired_t'):
  1381. if len(d1) != len(d2):
  1382. raise ValueError('paired t-test needs equal-length data (got %d and %d)' % (len(d1), len(d2)))
  1383. stat, p = ttest_rel(d1, d2, alternative=alternative)
  1384. n_str = 'n=%d pairs' % len(d1)
  1385. else:
  1386. raise ValueError('unknown test: %s (use "mannwhitney", "wilcoxon", "ttest", or "ttest_rel")' % test)
  1388. sig = p_to_star(p)
  1389. print('%s vs %s | %s (%s) | %s | stat=%.4g | p=%.4g | %s'
  1390. % (name1, name2, test, alternative, n_str, stat, p, sig))
  1392. if stat_file is not None:
  1393. new = not os.path.isfile(stat_file)
  1394. with open(stat_file, 'a') as f:
  1395. if new:
  1396. f.write('group1,group2,test,alternative,n,statistic,p_value,significance\n')
  1397. f.write('%s,%s,%s,%s,%s,%.6g,%.6g,%s\n'
  1398. % (name1, name2, test, alternative, n_str, stat, p, sig))
  1400. # draw the significance bracket on the panel
  1401. draw_sig_bracket(ax, x1, x2, p, star=star, color=color, y=y)
  1403. return {'name1': name1, 'name2': name2, 'test': test, 'alternative': alternative,
  1404. 'statistic': stat, 'p': p, 'significance': sig}
  1407. def _bh_fdr(p):
  1408. # Benjamini-Hochberg FDR-adjusted p-values (matches MATLAB f_FDR_correction.m)
  1409. p = np.asarray(p, dtype=float)
  1410. n = p.size
  1411. order = np.argsort(p) # ascending
  1412. ranked = p[order] * n / (np.arange(n) + 1) # p * n / rank
  1413. ranked = np.minimum.accumulate(ranked[::-1])[::-1] # monotone from the top
  1414. ranked = np.minimum(ranked, 1.0)
  1415. out = np.empty(n)
  1416. out[order] = ranked
  1417. return out
  1420. def isi_stats(groups, ax=None, draw='consecutive', color='black', posthoc='tukey'):
  1421. # Stats for the ISI correlation figure: one-way ANOVA across ALL ISI groups + post-hoc.
  1422. # groups: dict {label: (data, x_position)} from plot_fig_isi_corr_trials.
  1423. # posthoc: 'tukey' (Tukey HSD) OR 'fisher_fdr' = MATLAB f_dv_plot_anova1: Fisher LSD using the
  1424. # pooled ANOVA error MSE, t with df=n1+n2-1, two-tailed, then Benjamini-Hochberg FDR.
  1425. # draw: None, 'consecutive' (0.5-1,1-2,2-4), 'first' (vs first group), or 'all'.
  1426. labels = list(groups)
  1427. data = [np.asarray(groups[k][0], dtype=float) for k in labels]
  1428. data = [d[np.isfinite(d)] for d in data]
  1429. ncat = len(data)
  1431. F, p_anova = f_oneway(*data)
  1432. print('One-way ANOVA across %d groups: F=%.4g, p=%.4g' % (ncat, F, p_anova))
  1434. if posthoc in ('fisher_fdr', 'fisher', 'matlab'):
  1435. from scipy.stats import t as _tdist
  1436. ns = np.array([len(d) for d in data])
  1437. means = np.array([d.mean() for d in data])
  1438. N = int(ns.sum())
  1439. MSE = float(sum(((d - d.mean()) ** 2).sum() for d in data)) / (N - ncat) # pooled ANOVA error
  1440. pairs = [(i, j) for i in range(ncat) for j in range(i + 1, ncat)]
  1441. praw = np.array([2 * _tdist.sf(abs((means[i] - means[j]) / np.sqrt(MSE / ns[i] + MSE / ns[j])),
  1442. ns[i] + ns[j] - 1) for (i, j) in pairs])
  1443. padj = _bh_fdr(praw)
  1444. pmat = np.ones((ncat, ncat))
  1445. for k, (i, j) in enumerate(pairs):
  1446. pmat[i, j] = pmat[j, i] = padj[k]
  1447. method = 'Fisher LSD (pooled MSE) + Benjamini-Hochberg FDR'
  1448. else:
  1449. try:
  1450. from scipy.stats import tukey_hsd
  1451. pmat = np.asarray(tukey_hsd(*data).pvalue)
  1452. method = 'Tukey HSD'
  1453. except Exception:
  1454. m = ncat * (ncat - 1) // 2
  1455. pmat = np.ones((ncat, ncat))
  1456. for i in range(ncat):
  1457. for j in range(i + 1, ncat):
  1458. _, pp = ttest_ind(data[i], data[j], equal_var=False)
  1459. pmat[i, j] = pmat[j, i] = min(pp * m, 1.0)
  1460. method = 'pairwise Welch t-test + Bonferroni'
  1462. print('Post-hoc (%s):' % method)
  1463. for i in range(len(labels)):
  1464. for j in range(i + 1, len(labels)):
  1465. print(' %s vs %s: p=%.4g %s' % (labels[i], labels[j], pmat[i, j], p_to_star(pmat[i, j])))
  1467. if ax is not None and draw:
  1468. if draw == 'consecutive':
  1469. pairs = [(i, i + 1) for i in range(len(labels) - 1)]
  1470. elif draw == 'first':
  1471. pairs = [(0, j) for j in range(1, len(labels))]
  1472. elif draw == 'all':
  1473. pairs = [(i, j) for i in range(len(labels)) for j in range(i + 1, len(labels))]
  1474. else:
  1475. pairs = []
  1476. for (i, j) in pairs:
  1477. draw_sig_bracket(ax, groups[labels[i]][1], groups[labels[j]][1], pmat[i, j], color=color)
  1479. return {'labels': labels, 'F': F, 'p_anova': p_anova, 'pvalue': pmat, 'method': method}
  1481. def plot_fig_tau_networks(tau_ob_net_all, tau_rnn_net_all, training_type):
  1482. # two-panel network-tau figure (CaIm | RNN by training type) drawn as violins
  1484. fig, ax, = plt.subplots(1, 2, sharey=True, gridspec_kw={'width_ratios': [1, 3]}, figsize=(6,5))
  1486. fig.text(0.02, .89, 'A', fontsize=16)
  1487. fig.text(0.32, .89, 'B', fontsize=16)
  1488. fig.suptitle('Network Tau')
  1490. plot_int_violin2(tau_ob_net_all,
  1491. net_labels = None,
  1492. title_tag = 'CaIm',
  1493. axis=ax[0],
  1494. points=1000,
  1495. mean_std=True,
  1496. showmeans=True,
  1497. showmedians=False,
  1498. quantile = [0.05, 0.95],
  1499. colors=['blue', 'magenta', 'green'],
  1500. do_log=True)
  1503. plot_int_violin2(tau_rnn_net_all,
  1504. net_labels = training_type,
  1505. title_tag = 'RNN',
  1506. axis=ax[1],
  1507. points=1000,
  1508. mean_std=True,
  1509. showmeans=True,
  1510. showmedians=False,
  1511. quantile = [0.05, 0.95],
  1512. colors=['blue', 'magenta', 'green'],
  1513. do_log=True)
  1515. ax[1].set_ylabel(None)
  1517. return fig
  1519. def plot_fig_tau_networks2(tau_ob_net_all, tau_rnn_net_all, training_type):
  1520. # single-panel network-tau violins (CaIm data + each RNN training type on one axis)
  1522. fig, ax, = plt.subplots(1, 1, figsize=(6,5))
  1524. fig.text(0.02, .89, 'A', fontsize=16)
  1525. fig.text(0.32, .89, 'B', fontsize=16)
  1526. fig.suptitle('Network Tau')
  1528. data_all = [np.array(tau_ob_net_all).flatten()] + tau_rnn_net_all
  1529. labels_all = np.array(['Caim data'] + list(training_type))
  1531. plot_int_violin2(data_all,
  1532. net_labels = labels_all,
  1533. axis=ax,
  1534. points=1000,
  1535. mean_std=True,
  1536. showmeans=True,
  1537. showmedians=False,
  1538. quantile = [0.05, 0.95],
  1539. colors=['blue', 'green', 'orange', 'gray'],
  1540. do_log=True)
  1542. return fig
  1544. #%%
  1546. def plot_cat_data2(y_data_in, rnn_leg, title_tag = '', do_log=False):
  1547. # scatter of each category's values with a mean +/- std marker
  1549. num_cat = len(y_data_in)
  1551. plt.figure()
  1552. ax1 = plt.subplot(111)
  1553. ax1.bar(rnn_leg, np.zeros(num_cat))
  1554. for n_net in range(num_cat):
  1556. y_data = np.concatenate(y_data_in[n_net])
  1558. x_data = ((np.random.rand(len(y_data)))-0.5)/5+n_net
  1560. ax1.plot(x_data, y_data, '.', color='gray')
  1561. ax1.plot(n_net, np.mean(y_data), '_', color='black', mew=2, markersize=40)
  1562. ax1.errorbar(n_net, np.mean(y_data), np.std(y_data), fmt='o', color='black', mew=2, markersize=5, linewidth=2, capsize=10)
  1563. ax1.set_title(title_tag)
  1564. if do_log:
  1565. ax1.set_yscale('log')
  1567. def plot_cat_data_violin(y_data_in, rnn_leg, title_tag = '', points=100, mean_std=True, showmeans=False, showmedians=False, quantile = [], colors=[], do_log=False):
  1568. # violin plot of categorical groups, optional mean +/- std overlay and log y-axis
  1570. num_cat = len(y_data_in)
  1572. plt.figure()
  1573. ax1 = plt.subplot(111)
  1574. ax1.bar(rnn_leg, np.zeros(num_cat))
  1575. parts = ax1.violinplot(y_data_in, positions=range(num_cat), showmeans=showmeans, showextrema=False, showmedians=showmedians, quantiles=[quantile, quantile, quantile], points=points)
  1576. for key in ['cmeans', 'cmedians', 'cquantiles']:
  1577. if key in parts:
  1578. parts[key].set_color('k')
  1580. if len(colors):
  1581. for n_net in range(num_cat):
  1582. pc = parts['bodies'][n_net]
  1583. pc.set_facecolor(colors[n_net])
  1584. pc.set_edgecolor(colors[n_net])
  1585. if mean_std:
  1586. for n_net in range(num_cat):
  1587. y_data = y_data_in[n_net]
  1588. ax1.plot(n_net, np.mean(y_data), '_', color='black', mew=2, markersize=40)
  1589. ax1.errorbar(n_net, np.mean(y_data), np.std(y_data), fmt='o', color='black', mew=2, markersize=5, linewidth=2, capsize=10)
  1590. ax1.set_title(title_tag)
  1591. if do_log:
  1592. ax1.set_yscale('log')
  1594. def plot_cat_data_bar(y_data_in, rnn_leg, title_tag = '', do_sem=True, colors=[]):
  1595. # bar plot (mean +/- sem or std) of categorical groups
  1597. num_cat = len(y_data_in)
  1599. plt.figure()
  1600. plt.bar(rnn_leg, np.zeros(num_cat))
  1601. for n_net in range(num_cat):
  1602. y_data = y_data_in[n_net]
  1603. if do_sem:
  1604. stds = np.std(y_data)/np.sqrt(len(y_data)-1)
  1605. else:
  1606. stds = np.std(y_data)
  1607. if len(colors):
  1608. plt.bar(n_net, np.mean(y_data), color=colors[n_net], alpha=0.5, edgecolor=colors[n_net])
  1609. else:
  1610. plt.bar(n_net, np.mean(y_data))
  1611. plt.plot(n_net, np.mean(y_data), '_', color='black', mew=2, markersize=40)
  1612. plt.errorbar(n_net, np.mean(y_data), stds, fmt='o', color='black', mew=2, markersize=5, linewidth=2, capsize=10)
  1613. plt.title(title_tag)
  1615. def plot_int_violin2(tau_net_list, net_labels = None, data_lab = ['CaIm data'], title_tag = '', axis=None, points=100, mean_std=True, showmeans=False, showmedians=False, quantile = [0.05, 0.95], colors=['blue', 'green', 'orange', 'gray'], do_log=False):
  1616. # violin plot of tau distributions grouped by net_labels, with optional mean +/- std overlay and log y-axis
  1618. num_net = len(tau_net_list)
  1620. net_type_list = []
  1621. net_idx = np.zeros(num_net, dtype=int)
  1622. rnn_leg = []
  1623. quantiles_all = []
  1624. if net_labels is not None:
  1625. _, idx1 = np.unique(net_labels, return_index=True)
  1626. idx1.sort()
  1627. net_uq = net_labels[idx1]
  1628. for n_net in range(len(net_uq)):
  1629. net_idx[net_labels == net_uq[n_net]] = n_net
  1630. temp_bin = []
  1631. for n_net2 in range(num_net):
  1632. if net_uq[n_net] == net_labels[n_net2]:
  1633. temp_bin.append(tau_net_list[n_net2].flatten())
  1635. net_type_list.append(np.hstack(temp_bin).flatten())
  1636. rnn_leg.append(net_uq[n_net].capitalize())
  1637. quantiles_all.append(quantile)
  1639. else:
  1640. net_uq = data_lab
  1641. net_type_list.append(np.array(tau_net_list).flatten())
  1642. quantiles_all.append(quantile)
  1643. rnn_leg = data_lab
  1646. num_net_types = len(net_uq)
  1648. if axis is None:
  1649. plt.figure()
  1650. ax1 = plt.subplot(111)
  1651. else:
  1652. ax1 = axis
  1654. ax1.bar(rnn_leg, np.zeros(num_net_types))
  1655. parts = ax1.violinplot(net_type_list, positions=range(num_net_types), showmeans=showmeans, showextrema=False, showmedians=showmedians, quantiles=quantiles_all, points=points)
  1656. for key in ['cmeans', 'cmedians', 'cquantiles']:
  1657. if key in parts:
  1658. parts[key].set_color('k')
  1660. if len(colors):
  1661. for n_net in range(num_net_types):
  1662. pc = parts['bodies'][n_net]
  1663. pc.set_facecolor(colors[n_net])
  1664. pc.set_edgecolor(colors[n_net])
  1665. if mean_std:
  1666. for n_net in range(num_net_types):
  1667. y_data = net_type_list[n_net]
  1668. ax1.plot(n_net, np.mean(y_data), '_', color='black', mew=2, markersize=40)
  1669. ax1.errorbar(n_net, np.mean(y_data), np.std(y_data), fmt='o', color='black', mew=2, markersize=5, linewidth=2, capsize=10)
  1670. ax1.set_title(title_tag)
  1672. if do_log:
  1673. ax1.set_yscale('log')
  1674. ax1.set_ylabel('Tau (sec)')
  1675. else:
  1676. ax1.set_ylabel('Tau (sec)')

sd_utils.py at commit 53e1dfd, no license · at the source

Overview

Authors: Yuriy Shymkiv1, Rafael Yuste1
ORCID iDs: Yuriy Shymkiv
  1. Neurotechnology Center, Department of Biological Sciences, Columbia University, New York, NY, USA
Institutions: Columbia University (United States)
Journal: STAR protocols, volume 7, issue 3, article 104804
Dates: published online 27 August 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.xpro.2026.104804 · PMID 42658683 · PMCID PMC13544393 · OpenAlex W7204494996
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: mouse (organism)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, Evoked potentials, Connectivity, Single-unit activity, calcium imaging
Keywords: Neuroscience, Cognitive Neuroscience, Computer sciences
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: HHS | NIH | National Institute of Neurological Disorders and Stroke (NINDS) (RM1NS132981); HHS | NIH | National Eye Institute (NEI) (R01EY035248); HHS | NIH | National Institute of Mental Health (NIMH) (R01MH115900, F31MH122137, 5T32MH018870)
Citations: not cited yet (Europe PMC); 31 references in the paper
Research resources: Scikit-learn (v1.7.1) RRID:SCR_002577, SciPy (v1.16.3) RRID:SCR_008058, Python 3 (v3.12.12) RRID:SCR_008394, Matplotlib (v3.10.7) RRID:SCR_008624, NumPy (v2.3.5) RRID:SCR_008633, Jupyter Notebook (v7.5.0) RRID:SCR_018315, Miniconda RRID:SCR_018317, h5py (v3.15.1) RRID:SCR_024812

Abstract

Temporal information processing is critical for brain function, supporting neural computations such as novelty detection, adaptation, and temporal normalization. Its disruption is implicated in schizophrenia. We present a protocol for analyzing ongoing neuronal network activity using binwise decoding, trial-to-trial variability analysis, and estimation of network-intrinsic timescales (INTs). We apply these techniques to identify slow dynamics that encode the memory of recent stimuli in neuronal populations in the mouse auditory cortex and in artificial neural networks trained on a novelty-detection task.

For complete details on the use and execution of this protocol, please refer to Shymkiv et al.1

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 30 matches between paragraphs and lines of code.

shymkivy/slow_dynamics_protocol

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 53e1dfd5b9a9b053e0e8ad5c13f815a5861b789b, 10 July 2026
Languages: Jupyter (11), Python (7)
Size: 22 files, 18 scripts
Software Heritage: not archived
Found in: “Data and code availability”
Holds: README, 6 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (18 files), Matplotlib (12 files), SciPy (3 files), h5py (1 file), scikit-learn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
19 files
At the source: github.com/shymkivy/slow_dynamics_protocol

doi:10.5061/dryad.xsj3tx9q6

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: the link answers
Software Heritage: not checked
Found in: “Data and code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
At the source: doi.org/10.5061/dryad.xsj3tx9q6

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 18 scripts, each with its path and the digest of its content;
  • 30 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data and code availability

• The calcium imaging and RNN datasets analyzed in this protocol are publicly available at Dryad: https://doi.org/10.5061/dryad.xsj3tx9q6. • All original analysis code is publicly available on GitHub: https://github.com/shymkivy/slow_dynamics_protocol. • Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Reproduced under the paper's license (CC BY), from the paper cited above.

Materials availability

This study did not generate new unique reagents or materials. All software and datasets used are listed in the key resources table.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 3 keywords, 3 funders, 28 references, 8 RRIDs.

Cite

This paper

Shymkiv, Y., & Yuste, R. (2026). Protocol for analyzing slow cortical dynamics in mouse neuronal recordings. STAR protocols, 7(3), 104804. https://doi.org/10.1016/j.xpro.2026.104804

BibTeX

@article{shymkiv2026protocol,
author = {Shymkiv, Yuriy and Yuste, Rafael},
title = {{Protocol for analyzing slow cortical dynamics in mouse neuronal recordings}},
journal = {STAR protocols},
year = {2026},
month = aug,
volume = {7},
number = {3},
pages = {104804},
publisher = {Elsevier},
issn = {2666-1667},
doi = {10.1016/j.xpro.2026.104804},
url = {https://doi.org/10.1016/j.xpro.2026.104804},
pmid = {42658683},
pmcid = {PMC13544393}
}

RIS

TY - JOUR
AU - Shymkiv, Yuriy
AU - Yuste, Rafael
TI - Protocol for analyzing slow cortical dynamics in mouse neuronal recordings
T2 - STAR protocols
J2 - STAR Protoc
PY - 2026
DA - 2026/08/27
VL - 7
IS - 3
SP - 104804
SN - 2666-1667
PB - Elsevier
DO - 10.1016/j.xpro.2026.104804
UR - https://doi.org/10.1016/j.xpro.2026.104804
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.xpro.2026.104804",
"type": "article-journal",
"title": "Protocol for analyzing slow cortical dynamics in mouse neuronal recordings",
"container-title": "STAR protocols",
"author": [
{
"family": "Shymkiv",
"given": "Yuriy"
},
{
"family": "Yuste",
"given": "Rafael"
}
],
"container-title-short": "STAR Protoc",
"volume": "7",
"issue": "3",
"page": "104804",
"DOI": "10.1016/j.xpro.2026.104804",
"PMID": "42658683",
"PMCID": "PMC13544393",
"ISSN": "2666-1667",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.xpro.2026.104804",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
27
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1126/sciadv.aed6417 [code]
Intrinsic timing, not temporal prediction, underlies ramping dynamics in visual and parietal cortex during passive behavior.
Journal: Science advances
In common: h5py, scikit-learn, SciPy, 2 other tools, mouse, 3 references
[2] doi:10.1371/journal.pcbi.1014162 [code]
Exploring neural manifolds across a wide range of intrinsic dimensions.
Journal: PLoS computational biology
In common: h5py, scikit-learn, SciPy, 2 other tools, 3 references
[3] doi:10.1016/j.crmeth.2026.101421 [code]
EthoPy provides an accessible platform for reproducible behavioral neuroscience.
Journal: Cell reports methods
In common: h5py, scikit-learn, SciPy, 2 other tools, mouse, 3 references
[4] doi:10.1038/s41467-026-74816-0 [code]
Hierarchical optimization predicts plasticity in the macaque inferior temporal cortex following object training.
Journal: Nature communications
In common: h5py, scikit-learn, SciPy, 2 other tools, 3 references
[5] doi:10.3390/biomimetics11080569 [code]
Pretraining of Embodied Recurrent Networks Bridges the Gap Between Artificial and Cortical Neural Activities.
Journal: Biomimetics (Basel, Switzerland)
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 4 references
[6] doi:10.1038/s41593-026-02333-w [code]
Learning shapes neural geometry in the primate prefrontal cortex.
Journal: Nature neuroscience
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 3 references
[7] doi:10.1038/s41467-026-76106-1 [code]
Facial expression discrimination emerges from partially overlapping neural subspaces of detection and identity.
Journal: Nature communications
In common: h5py, scikit-learn, SciPy, 2 other tools, 2 references
[8] doi:10.1038/s41467-026-75924-7 [code]
Data-driven reduced modeling of neural dynamics.
Journal: Nature communications
In common: h5py, scikit-learn, SciPy, 2 other tools, 2 references
[9] doi:10.1016/j.isci.2026.117492 [code]
Neural subspace reorganization reflects value-based decision-making.
Journal: iScience
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 3 references
[10] doi:10.1038/s41467-026-74347-8 [code]
Compositionality of social gaze in the prefrontal-amygdala circuits.
Journal: Nature communications
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.