OSCR

Shallow recurrent decoders for neural and behavioural dynamics.

Code ↔ Paper

5 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 5 matches
  1. [1] § Model systems › Mouse › SHRED reconstruction of neural response to visual stimuli ↔ mice/shred_mice_expB_i_flash.ipynb, lines 560–629 · score 0.59 · SCig, MB, POL, VISp, APN, POST
  2. [2] § Model systems › Mouse › SHRED reconstruction of neural response to visual stimuli ↔ mice/shred_mice_expB_i_dg.ipynb, lines 589–658 · score 0.59 · SCig, MB, POL, VISp, APN, POST
  3. [3] § Model systems › Mouse › SHRED reconstruction of local field potential from pupil area response to a visual stimulus ↔ mice/shred_mice_expB_iv.ipynb, lines 125–150 · score 0.56 · cubic spline, pupil area, stimuli
  4. [4] § Model systems › Control comparison ↔ mice/shred_mice_expB_i_dg.ipynb, lines 756–826 · score 0.55 · power spectral density, linear regression, PSD, neurons, SHRED
  5. [5] § Model systems › Control comparison ↔ mice/shred_mice_expB_i_dg.ipynb, lines 756–826 · score 0.52 · Power spectral densities, linear regression, ground truth, PSDs, neurons, SHRED

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

Jupyter notebook · 1,010 lines · 35 KB · no license · 3 matches

  1. # %% [markdown]
  2. # # Applying SHRED to the Mice Dataset
  4. # %% [markdown]
  5. # Dataset obtained from Allen Brain Atlas https://allensdk.readthedocs.io/en/latest/visual_coding_neuropixels.html
  7. # %% [markdown]
  8. # SHRED (SHallow REcurrent Decoder) models are a network architecture that merges a recurrent layer (LSTM) with a shallow decoder network (SDN) to reconstruct high-dimensional spatio-temporal fields from a trajectory of sensor measurements of the field. More formally, the SHRED architecture can be written as
  9. # $$ \mathcal {H} \left( \{ y_i \} _{i=t-k}^t \right) = \mathcal {F} \left( \mathcal {G} \left( \{ y_i \} _{i=t-k}^t \right) ; W_{RN}) ; W_{SD} \right)$$
  10. # where $\mathcal F$ is a feed forward network parameterized by weights $W_{SD}$, $\mathcal G$ is a LSTM network parameterized by weights $W_{RN}$, and $\{ y_i \} _{i=t-k}^t$ is a trajectory of sensor measurements of a high-dimensional spatio-temporal field $\{ x_i \} _{i=t-k}^t$.
  12. # %% [markdown]
  13. # Goal: Reconstruct LFP potential of different brain regions from only one session 756029989 and 180 degree orientation drifting gratings
  14. #
  15. # LFP band:
  16. # - ~2.5 kHz original sample rate
  17. # - 1000 Hz analog lo-pass
  18. # - 625 Hz digital lo-pass
  19. # - NWB includes every 2nd sample and every 4th channel
  21. # %%
  22. # Importing all dependencies
  23. import numpy as np
  24. import torch
  25. import subprocess
  26. import os
  27. import torch
  28. import matplotlib.pyplot as plt
  29. from sklearn.preprocessing import MinMaxScaler
  30. import random
  31. import mne
  34. # %%
  35. import pandas as pd
  36. from scipy.ndimage import gaussian_filter
  37. from pathlib import Path
  38. import json
  39. from IPython.display import display
  40. from PIL import Image
  44. from allensdk.brain_observatory.ecephys.visualization import plot_mean_waveforms, plot_spike_counts, raster_plot
  45. from allensdk.brain_observatory.visualization import plot_running_speed
  47. # tell pandas to show all columns when we display a DataFrame
  48. pd.set_option("display.max_columns", None)
  50. # %%
  51. #Set up data cache
  52. from allensdk.brain_observatory.ecephys.ecephys_project_cache import EcephysProjectCache
  53. HOME_DIR = "/Users/amyrude/Downloads/Kutz_Research/SHRED_neuro/neuralSHRED/"
  54. data_directory = os.path.join(HOME_DIR, "mice/data") #where data will be stored
  55. manifest_path = os.path.join(data_directory, 'manifest.json')
  56. cache = EcephysProjectCache.from_warehouse(manifest = manifest_path)
  58. # %%
  59. sessions = cache.get_session_table()
  61. print('Total number of sessions: ' + str(len(sessions)))
  63. sessions.head()
  66. # Can also filter which mice to select
  67. filtered_sessions = sessions[(sessions.index == 756029989)]
  69. filtered_sessions.head()
  71. # %%
  72. #Loading dataset for specific session
  73. session_id = 756029989 # for example
  74. session = cache.get_session_data(session_id)
  76. print([attr_or_method for attr_or_method in dir(session) if attr_or_method[0] != '_'])
  79. # %% [markdown]
  80. # **Importing local field potential (LFP) data**
  82. # %%
  83. # list the probes recorded from in this session
  84. session.probes.head()
  86. # %%
  87. {session.probes.loc[probe_id].description :
  88. list(session.channels[session.channels.probe_id == probe_id].ecephys_structure_acronym.unique())
  89. for probe_id in session.probes.index.values}
  91. # %%
  92. probe_id = session.probes.index.values[2] ### Select 0 to probe A, 2 for probe C
  93. lfp = session.get_lfp(probe_id)
  95. # %%
  96. print(lfp)
  98. # %%
  99. # now use a utility to associate intervals of /rows with structures
  101. ## For probe C, mismatch between channels in lfp["channel"] and session.channels.index
  102. lfp_channels = lfp["channel"].values.tolist()
  103. valid_channels = [ch for ch in lfp_channels if ch in session.channels.index]
  104. structure_acronyms, intervals = session.channel_structure_intervals(valid_channels)
  105. interval_midpoints = [aa + (bb - aa) / 2 for aa, bb in zip(intervals[:-1], intervals[1:])]
  107. print(structure_acronyms)
  108. print(intervals)
  109. print(interval_midpoints)
  111. # %% [markdown]
  112. # Aligning the LFP data to a particular stimulus
  114. # %%
  115. stim_table = session.get_stimulus_table('drifting_gratings')
  116. stim_table[stim_table['stimulus_condition_id'] == 246].head(10)
  118. # %%
  119. # presentation_table = session.stimulus_presentations[session.stimulus_presentations.stimulus_name == 'drifting_gratings']
  120. presentation_table = session.stimulus_presentations[session.stimulus_presentations.stimulus_condition_id == 246]
  122. presentation_times = presentation_table.start_time.values
  123. presentation_ids = presentation_table.index.values
  125. # %%
  126. print(presentation_table)
  128. # %%
  131. # %%
  132. ### Probe A and C
  133. dt = 1/500
  134. sr = 1/dt
  135. trial_window = np.arange(0, 2, dt)
  139. time_selection = np.concatenate([trial_window + t for t in presentation_times])
  141. inds = pd.MultiIndex.from_product((presentation_ids, trial_window),
  142. names=('presentation_id', 'time_from_presentation_onset'))
  144. ds = lfp.sel(time = time_selection, method='nearest').to_dataset(name = 'aligned_lfp')
  145. ds = ds.assign(time=inds).unstack('time')
  147. aligned_lfp = ds['aligned_lfp']
  149. # %%
  150. print(aligned_lfp.shape)
  152. # %%
  153. fig, ax = plt.subplots()
  154. im = ax.imshow(aligned_lfp.mean(dim='presentation_id'), aspect='auto', origin='lower', vmin=-1e-4, vmax=1e-4)
  156. ax.set_yticks(intervals)
  157. ax.set_yticks(interval_midpoints, minor=True)
  158. ax.set_yticklabels(structure_acronyms, minor=True)
  159. plt.tick_params("y", which="major", labelleft=False, length=40)
  161. num_time_labels = 8
  162. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  163. # time_labels = [ f"{val:1.3}" for val in lfp["time"].values[window][time_label_indices]]
  164. # ax.set_xticks(time_label_indices + 0.5)
  165. # ax.set_xticklabels(time_labels)
  166. ax.set_xlabel("Sample Number", fontsize=20)
  167. plt.colorbar(im, fraction=0.036, pad=0.04)
  169. plt.show()
  172. ## just one
  173. fig, ax = plt.subplots()
  174. im = ax.imshow(aligned_lfp[:,0,:], aspect='auto', origin='lower', vmin=-0.5e-3, vmax=0.5e-3)
  176. ax.set_yticks(intervals)
  177. ax.set_yticks(interval_midpoints, minor=True)
  178. ax.set_yticklabels(structure_acronyms, minor=True)
  179. plt.tick_params("y", which="major", labelleft=False, length=40)
  181. num_time_labels = 8
  182. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  183. # time_labels = [ f"{val:1.3}" for val in lfp["time"].values[window][time_label_indices]]
  184. # ax.set_xticks(time_label_indices + 0.5)
  185. # ax.set_xticklabels(time_labels)
  186. ax.set_xlabel("Sample Number", fontsize=20)
  187. plt.colorbar(im, fraction=0.036, pad=0.04)
  189. plt.show()
  193. # %%
  194. ### save data as np arrays
  195. sens_max = 76 #74
  196. probe = 'c'
  198. lfp_np = aligned_lfp.data
  199. data = np.hstack((lfp_np[:sens_max,0,:], lfp_np[:sens_max,1,:], lfp_np[:sens_max,2,:]))
  200. print(data.shape)
  202. n_t = data.shape[1]
  203. n_s = data.shape[0]
  206. ## Plot the numpy array
  207. fig, ax = plt.subplots(figsize = (5,6))
  208. im = ax.imshow(data, aspect='auto', origin='lower', vmin=-0.4e-3, vmax=0.4e-3, cmap = 'Reds')
  210. ax.set_yticks(intervals)
  211. ax.set_yticks(interval_midpoints, minor=True)
  212. ax.set_yticklabels(structure_acronyms, minor=True)
  213. plt.tick_params("y", which="major", labelleft=False, length=40)
  214. for y in [0,1000, 2000]:
  215. plt.axvline(y, color='white', linestyle='--', linewidth = '1.25')
  216. ax.set_xlim((-15,data.shape[1]))
  217. ax.set_ylim((0,sens_max))
  219. num_time_labels = 8
  220. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  222. ax.set_xlabel("Sample Number", fontsize=20)
  223. plt.colorbar(im, fraction=0.036, pad=0.04)
  224. plt.show()
  228. # %% [markdown]
  229. # **Applying SHRED on Dataset**
  231. # %%
  232. #Importing packages
  234. os.environ["CUDA_VISIBLE_DEVICES"]="0"
  235. SHRED_DIR = os.path.join(HOME_DIR, "neuralSHRED/sindy-shred-main")
  236. import sys
  237. sys.path.append(SHRED_DIR)
  239. from processdata import load_data
  240. from processdata import TimeSeriesDataset
  241. import sindy
  242. import sindy_shred
  245. os.environ["CUDA_VISIBLE_DEVICES"]="0"
  247. # %%
  248. load_X = data.T
  250. latent_dim = 32
  251. poly_order = 1
  252. num_neurons = 3
  253. lags = 100
  254. test_val_size = 1000 - lags
  257. ## Ignore below
  258. include_sine = False
  259. library_dim = sindy.library_size(latent_dim, poly_order, include_sine, True)
  263. train_indices = np.arange(0, n_t - lags - test_val_size)
  264. mask = np.ones(n_t - lags)
  265. mask[train_indices] = 0
  266. valid_test_indices = np.arange(0, n_t - lags)[np.where(mask!=0)[0]]
  267. valid_indices = valid_test_indices[int(test_val_size/2):test_val_size]
  268. test_indices = valid_test_indices[:int(test_val_size/2)]
  271. # %%
  272. ## Plot the lag
  273. fig, ax = plt.subplots(figsize = (7,3))
  274. im = ax.imshow(data[:, :lags], aspect='auto', origin='lower', vmin=-0.4e-3, vmax=0.4e-3, cmap = 'Reds')
  276. ax.set_yticks(intervals)
  277. ax.set_yticks(interval_midpoints, minor=True)
  278. ax.set_yticklabels(structure_acronyms, minor=True)
  279. plt.tick_params("y", which="major", labelleft=False, length=40)
  280. # for y in [0,1000, 2000]:
  281. # plt.axvline(y, color='white', linestyle='--', linewidth = '1.25')
  282. # ax.set_xlim((-15,data.shape[1]))
  283. ax.set_ylim((0,sens_max))
  285. num_time_labels = 8
  286. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  288. ax.set_xlabel("Sample Number", fontsize=20)
  289. plt.colorbar(im, fraction=0.036, pad=0.04)
  290. plt.show()
  292. # %%
  293. # Randomly select the sensors from the CA1 region for probe A (35 - 51) and SUB for probe C (31 - 54)
  294. print(intervals)
  295. print(structure_acronyms)
  297. for k in range(4,5):
  298. neuron_locations = np.random.choice(np.arange(31, 54), size=num_neurons, replace=False)
  299. file_path = os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/neuron_loc_{probe}_{k}.npy")
  300. np.save(file_path, neuron_locations)
  301. sc = MinMaxScaler()
  302. sc = sc.fit(load_X[train_indices])
  303. transformed_X = sc.transform(load_X)
  305. ### Generate input sequences to a SHRED model
  306. all_data_in = np.zeros((n_t - lags, lags, num_neurons))
  307. for i in range(len(all_data_in)):
  308. all_data_in[i] = transformed_X[i:i+lags, neuron_locations]
  310. ### Generate training validation and test datasets both for reconstruction of states and forecasting sensors
  311. device = 'cuda' if torch.cuda.is_available() else 'cpu'
  313. train_data_in = torch.tensor(all_data_in[train_indices], dtype=torch.float32).to(device)
  314. valid_data_in = torch.tensor(all_data_in[valid_indices], dtype=torch.float32).to(device)
  315. test_data_in = torch.tensor(all_data_in[test_indices], dtype=torch.float32).to(device)
  317. ### -1 to have output be at the same time as final sensor measurements
  318. train_data_out = torch.tensor(transformed_X[train_indices + lags - 1], dtype=torch.float32).to(device)
  319. valid_data_out = torch.tensor(transformed_X[valid_indices + lags - 1], dtype=torch.float32).to(device)
  320. test_data_out = torch.tensor(transformed_X[test_indices + lags - 1], dtype=torch.float32).to(device)
  322. train_dataset = TimeSeriesDataset(train_data_in, train_data_out)
  323. valid_dataset = TimeSeriesDataset(valid_data_in, valid_data_out)
  324. test_dataset = TimeSeriesDataset(test_data_in, test_data_out)
  326. shred = sindy_shred.SINDy_SHRED(num_neurons, n_s, hidden_size=latent_dim, hidden_layers=2, l1=350, l2=400, dropout=0.1,
  327. library_dim=library_dim, poly_order=poly_order,
  328. include_sine=include_sine, dt=dt, layer_norm=False, sindy = False).to(device)
  330. validation_errors = sindy_shred.fit(shred, train_dataset, valid_dataset, batch_size=128, num_epochs=200,
  331. lr=1e-3, verbose=True, threshold=0.25, patience=5, sindy_regularization= 0.0,
  332. optimizer="AdamW", thres_epoch=100)
  334. test_recons = sc.inverse_transform(shred(test_dataset.X).detach().cpu().numpy())
  335. test_ground_truth = sc.inverse_transform(test_dataset.Y.detach().cpu().numpy())
  336. train_recons = sc.inverse_transform(shred(train_dataset.X).detach().cpu().numpy())
  337. train_ground_truth = sc.inverse_transform(train_dataset.Y.detach().cpu().numpy())
  339. mask = np.ones(train_ground_truth.shape[1], dtype=bool)
  340. mask[neuron_locations] = False
  341. train_ex = train_recons.T[mask]
  342. train_gt_ex = train_ground_truth.T[mask]
  343. mask = np.ones(test_ground_truth.shape[1], dtype=bool)
  344. mask[neuron_locations] = False
  345. test_ex = test_recons.T[mask]
  346. test_gt_ex = test_ground_truth.T[mask]
  347. mse_train = np.linalg.norm(train_ex - train_gt_ex) / np.linalg.norm(train_gt_ex)
  348. mse_test = np.linalg.norm(test_ex - test_gt_ex) / np.linalg.norm(test_gt_ex)
  350. print('mse test', mse_test)
  351. print('mse train', mse_train)
  354. file_path = os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_recons_{probe}_{k}.npy")
  355. np.save(file_path, train_recons)
  356. file_path = os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_recons_{probe}_{k}.npy")
  357. np.save(file_path, test_recons)
  358. file_path = os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_gt_{probe}_{k}.npy")
  359. np.save(file_path, train_ground_truth)
  360. file_path = os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_gt_{probe}_{k}.npy")
  361. np.save(file_path, test_ground_truth)
  362. file_path = os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/mse_{probe}_{k}.npy")
  363. np.save(file_path, [mse_train, mse_test])
  365. # %% [markdown]
  366. # **Plotting the Data**
  368. # %%
  369. ## Probe A
  370. structure_acronyms = ['APN', 'DG', 'CA1', 'VISam', 'nan']
  371. intervals = [0, 27, 35, 51, 74 ,87]
  372. interval_midpoints = [13.5, 31.0, 43.0, 62.5, 80.5]
  374. ### Load the data
  375. sens_max = 74
  376. probe = 'a'
  377. trial = 2
  378. train_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_recons_a_{trial}.npy"))
  379. test_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_recons_a_{trial}.npy"))
  380. data_recon = np.vstack((train_recon, test_recon))
  382. train_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_gt_a_{trial}.npy"))
  383. test_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_gt_a_{trial}.npy"))
  384. data_gt= np.vstack((train_gt, test_gt))
  386. neuron_loc = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/neuron_loc_a_{trial}.npy"))
  387. mse = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/mse_a_{trial}.npy"))
  389. ## Plot the Reconstruction
  390. fig, ax = plt.subplots(figsize = (6.5,3))
  391. im = ax.imshow(data_recon.T, aspect='auto', origin='lower', vmin=-0.4e-3, vmax=0.4e-3, cmap = 'YlOrRd')
  393. ax.set_yticks(intervals)
  394. ax.set_yticks(interval_midpoints, minor=True)
  395. ax.set_yticklabels(structure_acronyms, minor=True)
  396. for tick in ax.yaxis.get_minor_ticks():
  397. tick.label1.set_fontsize(12)
  399. plt.tick_params("y", which="major", labelleft=False, length=40)
  400. for y in [0,1000, 2000]:
  401. plt.axvline(y, color='white', linestyle='--', linewidth = '1.25')
  402. for k in range(len(neuron_loc)):
  403. ax.plot(-5, neuron_loc[k], 'o', color='red', zorder = 10, clip_on=False, markersize = 5)
  404. ax.set_clip_on(False)
  405. ax.set_xlim((-15,data_recon.shape[0]))
  406. ax.set_ylim((0,sens_max))
  407. plt.tick_params(axis='both', which='major', labelsize=13)
  409. num_time_labels = 8
  410. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  412. ax.set_xlabel("Sample Number", fontsize=15)
  413. cbar = plt.colorbar(im, fraction=0.036, pad=0.04)
  414. cbar.set_label("V", fontsize = 12)
  415. cbar.ax.tick_params(labelsize=11)
  417. plt.tight_layout()
  418. plt.savefig(os.path.join(HOME_DIR, "mice/data_output/figs/lfp_dg_a_recon.png"), transparent = True, dpi = 400)
  419. plt.show()
  421. ## Plot the Ground Truth
  422. fig, ax = plt.subplots(figsize = (6.5,3))
  423. im = ax.imshow(data_gt.T, aspect='auto', origin='lower', vmin=-0.4e-3, vmax=0.4e-3, cmap = 'YlOrRd')
  425. ax.set_yticks(intervals)
  426. ax.set_yticks(interval_midpoints, minor=True)
  427. ax.set_yticklabels(structure_acronyms, minor=True)
  428. for tick in ax.yaxis.get_minor_ticks():
  429. tick.label1.set_fontsize(12)
  431. plt.tick_params("y", which="major", labelleft=False, length=40)
  432. for y in [0,1000, 2000]:
  433. plt.axvline(y, color='white', linestyle='--', linewidth = '1.25')
  434. for k in range(len(neuron_loc)):
  435. ax.plot(-5, neuron_loc[k], 'o', color='red', zorder = 10, clip_on=False, markersize = 5)
  436. ax.set_clip_on(False)
  437. ax.set_xlim((-15,data_gt.shape[0]))
  438. ax.set_ylim((0,sens_max))
  439. plt.tick_params(axis='both', which='major', labelsize=13)
  441. num_time_labels = 8
  442. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  444. ax.set_xlabel("Sample Number", fontsize=15)
  445. cbar = plt.colorbar(im, fraction=0.036, pad=0.04)
  446. cbar.set_label("V", fontsize = 12)
  447. cbar.ax.tick_params(labelsize=11)
  448. plt.tight_layout()
  449. plt.savefig(os.path.join(HOME_DIR, "mice/data_output/figs/lfp_dg_a_gt.png"), transparent = True, dpi = 400)
  451. plt.show()
  452. print('mse', mse)
  456. # %%
  457. structure_acronyms = ['APN', 'DG', 'CA1', 'VISam', 'nan']
  458. intervals = [0, 27, 35, 51, 74 ,87]
  460. APN_train_gt = train_gt[0:8,:]
  461. DG_train_gt = train_gt[27:35,:]
  462. CA1_train_gt = train_gt[35:43,:]
  463. VISam_train_gt = train_gt[51:59,:]
  465. APN_test_gt = test_gt[0:8,:]
  466. DG_test_gt = test_gt[27:35,:]
  467. CA1_test_gt = test_gt[35:43,:]
  468. VISam_test_gt = test_gt[51:59,:]
  470. APN_train_recon = train_recon[0:8,:]
  471. DG_train_recon = train_recon[27:35,:]
  472. CA1_train_recon = train_recon[35:43,:]
  473. VISam_train_recon = train_recon[51:59,:]
  475. APN_test_recon = test_recon[0:8,:]
  476. DG_test_recon = test_recon[27:35,:]
  477. CA1_test_recon = test_recon[35:43,:]
  478. VISam_test_recon = test_recon[51:59,:]
  480. ### CA1 to CA1
  481. mse_ca1_train = np.linalg.norm(CA1_train_recon - CA1_train_gt) / np.linalg.norm(CA1_train_gt)
  482. mse_ca1_test = np.linalg.norm(CA1_test_recon - CA1_test_gt) / np.linalg.norm(CA1_test_gt)
  484. ### CA1 to VISam
  485. mse_visam_train = np.linalg.norm(VISam_train_recon - VISam_train_gt) / np.linalg.norm(VISam_train_gt)
  486. mse_visam_test = np.linalg.norm(VISam_test_recon - VISam_test_gt) / np.linalg.norm(VISam_test_gt)
  488. ### CA1 to DG
  489. mse_dg_train = np.linalg.norm(DG_train_recon - DG_train_gt) / np.linalg.norm(DG_train_gt)
  490. mse_dg_test = np.linalg.norm(DG_test_recon - DG_test_gt) / np.linalg.norm(DG_test_gt)
  492. ### CA1 to APN
  493. mse_apn_train = np.linalg.norm(APN_train_recon - APN_train_gt) / np.linalg.norm(APN_train_gt)
  494. mse_apn_test = np.linalg.norm(APN_test_recon - APN_test_gt) / np.linalg.norm(APN_test_gt)
  496. print('mse visam train', mse_visam_train, mse_visam_test)
  497. print('mse dg train', mse_dg_train, mse_dg_test)
  498. print('mse apn train', mse_apn_train, mse_apn_test)
  499. print('mse ca1 train', mse_ca1_train, mse_ca1_test)
  504. # %%
  505. ## Probe C
  506. structure_acronyms = ['POL', '', '', 'SCig', '', 'SUB', 'VISp', 'nan']
  507. intervals = [ 0, 5 , 9 ,12 ,29, 31 ,54, 76 ,81]
  508. interval_midpoints = [2.5, 7.0, 10.5, 20.5, 30.0, 42.5, 65.0, 78.5]
  510. ### Load the data
  511. sens_max = 76
  512. probe = 'c'
  513. trial = 3 ## trial 4- 64 latent dimension
  514. train_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_recons_c_{trial}.npy"))
  515. test_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_recons_c_{trial}.npy"))
  516. data_recon = np.vstack((train_recon, test_recon))
  518. train_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_gt_c_{trial}.npy"))
  519. test_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_gt_c_{trial}.npy"))
  520. data_gt= np.vstack((train_gt, test_gt))
  522. neuron_loc = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/neuron_loc_c_{trial}.npy"))
  523. mse = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/mse_c_{trial}.npy"))
  525. ## Plot the Reconstruction
  526. fig, ax = plt.subplots(figsize = (6.5,3))
  527. im = ax.imshow(data_recon.T, aspect='auto', origin='lower', vmin=-0.4e-3, vmax=0.4e-3, cmap = 'GnBu')
  529. ax.set_yticks(intervals)
  530. ax.set_yticks(interval_midpoints, minor=True)
  531. ax.set_yticklabels(structure_acronyms, minor=True)
  532. for tick in ax.yaxis.get_minor_ticks():
  533. tick.label1.set_fontsize(12)
  535. plt.tick_params("y", which="major", labelleft=False, length=40)
  536. for y in [0,1000, 2000]:
  537. plt.axvline(y, color='white', linestyle='--', linewidth = '1.25')
  538. for k in range(len(neuron_loc)):
  539. ax.plot(-5, neuron_loc[k], 'o', color='red', zorder = 10, clip_on=False, markersize = 5)
  540. ax.set_clip_on(False)
  541. ax.set_xlim((-15,data_recon.shape[0]))
  542. ax.set_ylim((0,sens_max))
  543. plt.tick_params(axis='both', which='major', labelsize=13)
  545. num_time_labels = 8
  546. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  548. ax.set_xlabel("Sample Number", fontsize=15)
  549. cbar = plt.colorbar(im, fraction=0.036, pad=0.04)
  550. cbar.ax.tick_params(labelsize=11)
  551. cbar.set_label("V", fontsize = 12)
  552. plt.tight_layout()
  553. plt.savefig(os.path.join(HOME_DIR, "mice/data_output/figs/lfp_dg_c_recon.png"), transparent = True, dpi = 400)
  554. plt.show()
  556. ## Plot the Ground Truth
  557. fig, ax = plt.subplots(figsize = (6.5,3))
  558. im = ax.imshow(data_gt.T, aspect='auto', origin='lower', vmin=-0.4e-3, vmax=0.4e-3, cmap = 'GnBu')
  560. ax.set_yticks(intervals)
  561. ax.set_yticks(interval_midpoints, minor=True)
  562. ax.set_yticklabels(structure_acronyms, minor=True)
  563. for tick in ax.yaxis.get_minor_ticks():
  564. tick.label1.set_fontsize(12)
  566. plt.tick_params("y", which="major", labelleft=False, length=40)
  567. for y in [0,1000, 2000]:
  568. plt.axvline(y, color='white', linestyle='--', linewidth = '1.25')
  569. for k in range(len(neuron_loc)):
  570. ax.plot(-5, neuron_loc[k], 'o', color='red', zorder = 10, clip_on=False, markersize = 5)
  571. ax.set_clip_on(False)
  572. ax.set_xlim((-15,data_gt.shape[0]))
  573. ax.set_ylim((0,sens_max))
  574. plt.tick_params(axis='both', which='major', labelsize=13)
  576. num_time_labels = 8
  577. time_label_indices = np.around(np.linspace(1, len(trial_window), num_time_labels)).astype(int) - 1
  579. ax.set_xlabel("Sample Number", fontsize=15)
  580. cbar = plt.colorbar(im, fraction=0.036, pad=0.04)
  581. cbar.ax.tick_params(labelsize=11)
  582. cbar.set_label("V", fontsize = 12)
  583. plt.tight_layout()
  584. plt.savefig(os.path.join(HOME_DIR, "mice/data_output/figs/lfp_dg_c_gt.png"), transparent = True, dpi = 400)
  586. plt.show()
  587. print('mse', mse)
  589. # %%
  590. structure_acronyms = ['POL', 'APN', 'MB', 'SCig', 'POST', 'SUB', 'VISp', 'nan']
  591. intervals = [ 0, 5 , 9 ,12 ,29, 31 ,54, 76 ,81]
  593. POL_train_gt = train_gt[0:2,:]
  594. APN_train_gt = train_gt[5:7,:]
  595. MB_train_gt = train_gt[9:11,:]
  596. SCig_train_gt = train_gt[12:14,:]
  597. POST_train_gt = train_gt[29:31,:]
  598. SUB_train_gt = train_gt[31:33,:]
  599. VISp_train_gt = train_gt[54:57,:]
  601. POL_train_recon = train_recon[0:2,:]
  602. APN_train_recon = train_recon[5:7,:]
  603. MB_train_recon = train_recon[9:11,:]
  604. SCig_train_recon = train_recon[12:14,:]
  605. POST_train_recon = train_recon[29:31,:]
  606. SUB_train_recon = train_recon[31:33,:]
  607. VISp_train_recon = train_recon[54:57,:]
  609. POL_test_gt = test_gt[0:2,:]
  610. APN_test_gt = test_gt[5:7,:]
  611. MB_test_gt = test_gt[9:11,:]
  612. SCig_test_gt = test_gt[12:14,:]
  613. POST_test_gt = test_gt[29:31,:]
  614. SUB_test_gt = test_gt[31:33,:]
  615. VISp_test_gt = test_gt[54:57,:]
  617. POL_test_recon = test_recon[0:2,:]
  618. APN_test_recon = test_recon[5:7,:]
  619. MB_test_recon = test_recon[9:11,:]
  620. SCig_test_recon = test_recon[12:14,:]
  621. POST_test_recon = test_recon[29:31,:]
  622. SUB_test_recon = test_recon[31:33,:]
  623. VISp_test_recon = test_recon[54:57,:]
  625. ### APN
  626. mse_apn_train = np.linalg.norm(APN_train_recon - APN_train_gt) / np.linalg.norm(APN_train_gt)
  627. mse_apn_test = np.linalg.norm(APN_test_recon - APN_test_gt) / np.linalg.norm(APN_test_gt)
  628. print('apn', mse_apn_train, mse_apn_test)
  630. ### VISp
  631. mse_visp_train = np.linalg.norm(VISp_train_recon - VISp_train_gt) / np.linalg.norm(VISp_train_gt)
  632. mse_visp_test = np.linalg.norm(VISp_test_recon - VISp_test_gt) / np.linalg.norm(VISp_test_gt)
  633. print('visp', mse_visp_train, mse_visp_test)
  635. ### MB
  636. mse_mb_train = np.linalg.norm(MB_train_recon - MB_train_gt) / np.linalg.norm(MB_train_gt)
  637. mse_mb_test = np.linalg.norm(MB_test_recon - MB_test_gt) / np.linalg.norm(MB_test_gt)
  638. print('mb', mse_mb_train, mse_mb_test)
  640. ### SCig
  641. mse_scig_train = np.linalg.norm(SCig_train_recon - SCig_train_gt) / np.linalg.norm(SCig_train_gt)
  642. mse_scig_test = np.linalg.norm(SCig_test_recon - SCig_test_gt) / np.linalg.norm(SCig_test_gt)
  643. print('scig', mse_scig_train, mse_scig_test)
  645. ### POL
  646. mse_pol_train = np.linalg.norm(POL_train_recon - POL_train_gt) / np.linalg.norm(POL_train_gt)
  647. mse_pol_test = np.linalg.norm(POL_test_recon - POL_test_gt) / np.linalg.norm(POL_test_gt)
  648. print('pol', mse_pol_train, mse_pol_test)
  650. ##SUB
  651. mse_sub_train = np.linalg.norm(SUB_train_recon - SUB_train_gt) / np.linalg.norm(SUB_train_gt)
  652. mse_sub_test = np.linalg.norm(SUB_test_recon - SUB_test_gt) / np.linalg.norm(SUB_test_gt)
  653. print('sub', mse_sub_train, mse_sub_test)
  655. ## POST
  656. mse_post_train = np.linalg.norm(POST_train_recon - POST_train_gt) / np.linalg.norm(POST_train_gt)
  657. mse_post_test = np.linalg.norm(POST_test_recon - POST_test_gt) / np.linalg.norm(POST_test_gt)
  658. print('post', mse_post_train, mse_post_test)
  660. # %%
  661. #### Linear Regression
  662. ## Probe A
  664. trial = 2
  665. train_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_recons_a_{trial}.npy"))
  666. test_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_recons_a_{trial}.npy"))
  667. data_recon = np.vstack((train_recon, test_recon))
  669. train_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_gt_a_{trial}.npy"))
  670. test_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_gt_a_{trial}.npy"))
  671. data_gt= np.vstack((train_gt, test_gt))
  673. neuron_loc = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/neuron_loc_a_{trial}.npy"))
  674. mse = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/mse_a_{trial}.npy"))
  677. ##### Linear Regression
  679. from sklearn.linear_model import LinearRegression
  683. train_in = train_gt[:, neuron_loc]
  684. test_in = test_gt[:, neuron_loc]
  685. data_in = data_gt[:, neuron_loc]
  686. train_out = train_gt
  690. model = LinearRegression()
  691. model.fit(train_in, train_out)
  693. test_out = model.predict(test_in)
  694. print('test out shape',test_out.shape)
  695. data_out = model.predict(data_in)
  696. print('data out shape',data_out.shape)
  697. plt.figure(figsize = (6,2))
  698. plt.plot(test_out)
  699. plt.title('Linear Regression')
  700. plt.show()
  702. plt.figure(figsize = (6,2))
  703. plt.plot(test_gt)
  704. plt.title('Ground Truth')
  705. plt.show()
  707. ######### Plotting the reconstruction
  708. fig, ax = plt.subplots(figsize = (6.5,2))
  709. plt.imshow(data_gt.T, aspect = "auto")
  710. plt.xlabel('Time (s)', fontsize = 12)
  711. plt.ylabel('Neuron', fontsize =12)
  712. # Add red dots for neuron locations on y-axis
  713. for loc in neuron_loc:
  714. ax.plot(0.1, loc, 'ro', markersize=5, clip_on=False, zorder=10)
  716. plt.colorbar()
  717. plt.tight_layout()
  718. # plt.savefig(os.path.join(HOME_DIR, f'worms/data_output/figs/{name}_recon.png'), transparent = True, dpi = 400)
  720. plt.show()
  722. print('out shape',data_out.shape)
  723. fig, ax = plt.subplots(figsize = (6.5,2))
  724. plt.imshow(data_out.T, aspect = "auto")
  725. plt.xlabel('Time (s)', fontsize = 12)
  726. plt.ylabel('Neuron', fontsize = 12)
  727. plt.colorbar()
  728. # Add red dots for neuron locations on y-axis
  729. for loc in neuron_loc:
  730. ax.plot(0.1, loc, 'ro', markersize=5, clip_on=False, zorder=10)
  731. plt.tight_layout()
  732. # plt.savefig(os.path.join(HOME_DIR, f'worms/data_output/figs/{name}_linear.png'), transparent = True, dpi = 400)
  733. plt.show()
  735. n_train = train_out.shape[0]
  736. mask = np.ones(train_out.shape[1], dtype=bool)
  737. mask[neuron_loc] = False
  738. train_ex = data_out[:n_train,:].T[mask]
  739. train_gt_ex = train_gt.T[mask]
  740. mask = np.ones(test_out.shape[1], dtype=bool)
  741. mask[neuron_loc] = False
  742. test_ex = test_out.T[mask]
  743. test_gt_ex = test_gt.T[mask]
  744. mse_train = np.linalg.norm(train_ex - train_gt_ex) / np.linalg.norm(train_gt_ex)
  745. mse_test = np.linalg.norm(test_ex - test_gt_ex) / np.linalg.norm(test_gt_ex)
  746. print('mse train, test', [mse_train, mse_test])
  749. mask = np.ones(train_out.shape[1], dtype=bool)
  750. mask[neuron_loc] = False
  751. data_ex = data_out.T[mask]
  752. data_gt_ex = data_gt.T[mask]
  753. mse_whole = np.linalg.norm(data_ex - data_gt_ex) / np.linalg.norm(data_gt_ex)
  754. print('mse whole', mse_whole)
  756. # %%
  757. # Compute power spectral density for each dataset
  758. from scipy import signal
  759. import numpy as np
  760. name = 'expB_i_dg'
  761. # Parameters for PSD computation
  762. nperseg = min(256, data_out.shape[0] // 4) # Window length for Welch's method
  763. noverlap = nperseg // 2
  764. time = np.linspace(0, 2*data_out.shape[0]/3000, data_out.shape[0])
  765. # Initialize lists to store PSDs
  766. psds_linear = []
  767. psds_gt = []
  768. psds_shred = []
  769. freqs = None
  770. # Compute PSD for each neuron/channel
  771. for i in range(data_out.shape[1]):
  772. # Linear regression PSD
  773. f, psd_linear = signal.welch(data_out[:, i], fs=1/(time[1]-time[0]),
  774. nperseg=nperseg, noverlap=noverlap)
  775. psds_linear.append(psd_linear)
  777. # Ground truth PSD
  778. f, psd_gt = signal.welch(data_gt[:, i], fs=1/(time[1]-time[0]),
  779. nperseg=nperseg, noverlap=noverlap)
  780. psds_gt.append(psd_gt)
  782. # SHRED PSD
  783. f, psd_shred = signal.welch(data_recon[:, i], fs=1/(time[1]-time[0]),
  784. nperseg=nperseg, noverlap=noverlap)
  785. psds_shred.append(psd_shred)
  787. if freqs is None:
  788. freqs = f
  790. # Convert to arrays
  791. psds_linear = np.array(psds_linear)
  792. psds_gt = np.array(psds_gt)
  793. psds_shred = np.array(psds_shred)
  794. # Compute mean and standard deviation across neurons
  795. mean_psd_linear = np.mean(psds_linear, axis=0)
  796. std_psd_linear = np.std(psds_linear, axis=0)
  798. mean_psd_gt = np.mean(psds_gt, axis=0)
  799. std_psd_gt = np.std(psds_gt, axis=0)
  801. mean_psd_shred = np.mean(psds_shred, axis=0)
  802. std_psd_shred = np.std(psds_shred, axis=0)
  804. # Plot the power spectral densities with smooth distributions
  805. fig, ax = plt.subplots(figsize=(3.5, 2))
  806. neuron = 10
  807. # Plot mean with shaded regions for standard deviation
  808. ax.loglog(freqs, mean_psd_linear, 'b-', linewidth=2, label='Linear Regression', alpha=0.8)
  809. ax.fill_between(freqs, mean_psd_linear - std_psd_linear, mean_psd_linear + std_psd_linear,
  810. alpha=0.3, color='blue')
  812. ax.loglog(freqs, mean_psd_gt, 'k-', linewidth=2, label='Ground Truth', alpha=0.8)
  813. ax.fill_between(freqs, mean_psd_gt - std_psd_gt, mean_psd_gt + std_psd_gt,
  814. alpha=0.3, color='black')
  816. ax.loglog(freqs, mean_psd_shred, 'r-', linewidth=2, label='SHRED', alpha=0.8)
  817. ax.fill_between(freqs, mean_psd_shred - std_psd_shred, mean_psd_shred + std_psd_shred,
  818. alpha=0.3, color='red')
  820. ax.set_xlabel('Frequency (Hz)', fontsize=12)
  821. ax.set_ylabel('PSD', fontsize=12)
  822. # ax.legend(fontsize=12)
  824. plt.tight_layout()
  825. plt.savefig(os.path.join(HOME_DIR, f'mice/data_output/figs/{name}_psd.png'), transparent = True, dpi = 400)
  826. plt.show()
  829. # %%
  830. mse_linear = np.linalg.norm(mean_psd_linear - mean_psd_gt) / np.linalg.norm(mean_psd_gt)
  831. mse_shred = np.linalg.norm(mean_psd_shred - mean_psd_gt) / np.linalg.norm(mean_psd_gt)
  832. print(mse_shred)
  833. print(mse_linear)
  835. # %%
  836. #### Linear Regression
  837. ## Probe C
  839. trial = 3
  840. train_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_recons_c_{trial}.npy"))
  841. test_recon = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_recons_c_{trial}.npy"))
  842. data_recon = np.vstack((train_recon, test_recon))
  844. train_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/train_gt_c_{trial}.npy"))
  845. test_gt = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/test_gt_c_{trial}.npy"))
  846. data_gt= np.vstack((train_gt, test_gt))
  848. neuron_loc = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/neuron_loc_c_{trial}.npy"))
  849. mse = np.load(os.path.join(HOME_DIR, f"mice/data_output/exp_Bi/dg/mse_c_{trial}.npy"))
  852. ##### Linear Regression
  854. from sklearn.linear_model import LinearRegression
  858. train_in = train_gt[:, neuron_loc]
  859. test_in = test_gt[:, neuron_loc]
  860. data_in = data_gt[:, neuron_loc]
  861. train_out = train_gt
  865. model = LinearRegression()
  866. model.fit(train_in, train_out)
  868. test_out = model.predict(test_in)
  869. print('test out shape',test_out.shape)
  870. data_out = model.predict(data_in)
  871. print('data out shape',data_out.shape)
  872. plt.figure(figsize = (6,2))
  873. plt.plot(test_out)
  874. plt.title('Linear Regression')
  875. plt.show()
  877. plt.figure(figsize = (6,2))
  878. plt.plot(test_gt)
  879. plt.title('Ground Truth')
  880. plt.show()
  882. ######### Plotting the reconstruction
  883. fig, ax = plt.subplots(figsize = (6.5,2))
  884. plt.imshow(data_gt.T, aspect = "auto")
  885. plt.xlabel('Time (s)', fontsize = 12)
  886. plt.ylabel('Neuron', fontsize =12)
  887. # Add red dots for neuron locations on y-axis
  888. for loc in neuron_loc:
  889. ax.plot(0.1, loc, 'ro', markersize=5, clip_on=False, zorder=10)
  891. plt.colorbar()
  892. plt.tight_layout()
  893. # plt.savefig(os.path.join(HOME_DIR, f'worms/data_output/figs/{name}_recon.png'), transparent = True, dpi = 400)
  895. plt.show()
  897. print('out shape',data_out.shape)
  898. fig, ax = plt.subplots(figsize = (6.5,2))
  899. plt.imshow(data_out.T, aspect = "auto")
  900. plt.xlabel('Time (s)', fontsize = 12)
  901. plt.ylabel('Neuron', fontsize = 12)
  902. plt.colorbar()
  903. # Add red dots for neuron locations on y-axis
  904. for loc in neuron_loc:
  905. ax.plot(0.1, loc, 'ro', markersize=5, clip_on=False, zorder=10)
  906. plt.tight_layout()
  907. # plt.savefig(os.path.join(HOME_DIR, f'worms/data_output/figs/{name}_linear.png'), transparent = True, dpi = 400)
  908. plt.show()
  910. n_train = train_out.shape[0]
  911. mask = np.ones(train_out.shape[1], dtype=bool)
  912. mask[neuron_loc] = False
  913. train_ex = data_out[:n_train,:].T[mask]
  914. train_gt_ex = train_gt.T[mask]
  915. mask = np.ones(test_out.shape[1], dtype=bool)
  916. mask[neuron_loc] = False
  917. test_ex = test_out.T[mask]
  918. test_gt_ex = test_gt.T[mask]
  919. mse_train = np.linalg.norm(train_ex - train_gt_ex) / np.linalg.norm(train_gt_ex)
  920. mse_test = np.linalg.norm(test_ex - test_gt_ex) / np.linalg.norm(test_gt_ex)
  921. print('mse train, test', [mse_train, mse_test])
  924. mask = np.ones(train_out.shape[1], dtype=bool)
  925. mask[neuron_loc] = False
  926. data_ex = data_out.T[mask]
  927. data_gt_ex = data_gt.T[mask]
  928. mse_whole = np.linalg.norm(data_ex - data_gt_ex) / np.linalg.norm(data_gt_ex)
  929. print('mse whole', mse_whole)
  931. # %%
  932. #Compute power spectral density for each dataset
  933. from scipy import signal
  934. import numpy as np
  935. name = 'expB_i_dg_c'
  936. # Parameters for PSD computation
  937. nperseg = min(256, data_out.shape[0] // 4) # Window length for Welch's method
  938. noverlap = nperseg // 2
  939. time = np.linspace(0, 2*data_out.shape[0]/3000, data_out.shape[0])
  940. # Initialize lists to store PSDs
  941. psds_linear = []
  942. psds_gt = []
  943. psds_shred = []
  944. freqs = None
  945. # Compute PSD for each neuron/channel
  946. for i in range(data_out.shape[1]):
  947. # Linear regression PSD
  948. f, psd_linear = signal.welch(data_out[:, i], fs=1/(time[1]-time[0]),
  949. nperseg=nperseg, noverlap=noverlap)
  950. psds_linear.append(psd_linear)
  952. # Ground truth PSD
  953. f, psd_gt = signal.welch(data_gt[:, i], fs=1/(time[1]-time[0]),
  954. nperseg=nperseg, noverlap=noverlap)
  955. psds_gt.append(psd_gt)
  957. # SHRED PSD
  958. f, psd_shred = signal.welch(data_recon[:, i], fs=1/(time[1]-time[0]),
  959. nperseg=nperseg, noverlap=noverlap)
  960. psds_shred.append(psd_shred)
  962. if freqs is None:
  963. freqs = f
  965. # Convert to arrays
  966. psds_linear = np.array(psds_linear)
  967. psds_gt = np.array(psds_gt)
  968. psds_shred = np.array(psds_shred)
  969. # Compute mean and standard deviation across neurons
  970. mean_psd_linear = np.mean(psds_linear, axis=0)
  971. std_psd_linear = np.std(psds_linear, axis=0)
  973. mean_psd_gt = np.mean(psds_gt, axis=0)
  974. std_psd_gt = np.std(psds_gt, axis=0)
  976. mean_psd_shred = np.mean(psds_shred, axis=0)
  977. std_psd_shred = np.std(psds_shred, axis=0)
  979. # Plot the power spectral densities with smooth distributions
  980. fig, ax = plt.subplots(figsize=(3.5, 2))
  981. neuron = 10
  982. # Plot mean with shaded regions for standard deviation
  983. ax.loglog(freqs, mean_psd_linear, 'b-', linewidth=2, label='Linear Regression', alpha=0.8)
  984. ax.fill_between(freqs, mean_psd_linear - std_psd_linear, mean_psd_linear + std_psd_linear,
  985. alpha=0.3, color='blue')
  987. ax.loglog(freqs, mean_psd_gt, 'k-', linewidth=2, label='Ground Truth', alpha=0.8)
  988. ax.fill_between(freqs, mean_psd_gt - std_psd_gt, mean_psd_gt + std_psd_gt,
  989. alpha=0.3, color='black')
  991. ax.loglog(freqs, mean_psd_shred, 'r-', linewidth=2, label='SHRED', alpha=0.8)
  992. ax.fill_between(freqs, mean_psd_shred - std_psd_shred, mean_psd_shred + std_psd_shred,
  993. alpha=0.3, color='red')
  995. ax.set_xlabel('Frequency (Hz)', fontsize=12)
  996. ax.set_ylabel('PSD', fontsize=12)
  997. # ax.legend(fontsize=12)
  999. plt.tight_layout()
  1000. plt.savefig(os.path.join(HOME_DIR, f'mice/data_output/figs/{name}_psd.png'), transparent = True, dpi = 400)
  1001. plt.show()
  1003. # %%
  1004. mse_linear = np.linalg.norm(mean_psd_linear - mean_psd_gt) / np.linalg.norm(mean_psd_gt)
  1005. mse_shred = np.linalg.norm(mean_psd_shred - mean_psd_gt) / np.linalg.norm(mean_psd_gt)
  1006. print(mse_shred)
  1007. print(mse_linear)
  1009. # %%

shred_mice_expB_i_dg.ipynb at commit 6b86875, no license · at the source

Overview

Authors: Amy Rude1, J Nathan Kutz1,2
ORCID iDs: J Nathan Kutz
  1. Department of Applied Mathematics, University of Washington, Seattle, WA, USA
  2. Autodesk Research, London, UK
Institutions: University of Washington (United States); Autodesk (United Kingdom) (United Kingdom)
Journal: Philosophical transactions of the Royal Society of London. Series B, Biological sciences, volume 381, issue 1958, article 20240461
Dates: received 12 June 2025; accepted 17 October 2025; published online 17 September 2026; in print September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1098/rstb.2024.0461 · PMID 42750452 · PMCID PMC13583487 · OpenAlex W7213464516
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), mouse (organism), C. elegans (organism), computational (subfield)
Methods: Spectral & time-frequency, Connectivity, Statistics, Graphs, Machine learning, Single-unit activity, calcium imaging, Physiology & signal measures
Keywords: computational neuroscience, machine learning, population codes, shallow recurrent decoding, sensing
MeSH: Brain*, Machine Learning*, Models, Neurological*, Neurons*, Animals, Caenorhabditis elegans, Humans, Mice, Recurrent Neural Networks (* major topic)
Topic: Genetics, Aging, and Longevity in Model Organisms (Aging, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Air Force Office of Scientific Research (FA9550-24-1-0141); National Science Foundation AI Institute in Dynamic Systems (2112085)
Citations: not cited yet (Europe PMC); 68 references in the paper

Abstract

Machine learning algorithms are affording new opportunities for building bio-inspired and data-driven models characterizing neural activity. Critical to understanding decision-making and behaviour is quantifying the relationship between the activity of neuronal population codes and individual neurons. We leverage a SHallow REcurrent Decoder (SHRED) architecture for mapping the dynamics of population codes to individual neurons and other proxy measures of neural activity and behaviour. SHRED is constructed from a temporal sequence model, which encodes the temporal dynamics of limited sensor data in multiple scenarios, and a shallow decoder, which reconstructs the corresponding high-dimensional neuronal and/or behavioural states. It is a robust and flexible sensing strategy which allows for decoding the diversity of neural measurements with only a few sensor measurements. Thus, estimates of whole-brain activity, behaviour and individual neurons can be constructed with only a few neural time-series recordings. Several examples in this article further highlight the potential of leveraging non-invasive or minimally invasive measurements to estimate large-scale brain dynamics. We empirically demonstrate the capabilities of the method on a number of model organisms including Caenorhabditis elegans, mouse, zebrafish and human biolocomotion.

This article is part of the discussion meeting issue ‘Digital healthcare for the management of functional neurological disorders’.

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 5 matches between paragraphs and lines of code.

amysrude/neuralSHRED

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 6b8687501a637ec64ba2ffe65313c6acced4c83a, 15 September 2025
Languages: Jupyter (16), Python (8)
Size: 1,819 files, 24 scripts
Software Heritage: not archived
Found in: “Data accessibility”
Holds: README, environment (sindy-shred-main/requirements.txt), 16 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (19 files), PyTorch (19 files), SciPy (16 files), scikit-learn (14 files), Matplotlib (12 files), MNE-Python (11 files), AllenSDK (7 files), pandas (7 files), Pillow (7 files), xarray (6 files)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
21 files
At the source: github.com/amysrude/neuralSHRED

meganebers/mobileSHRED

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 66e85893bea3a463e18aaaf8b087387c64561018, 7 March 2024
Languages: Python (5)
Size: 12 files, 5 scripts
Software Heritage: not archived
Found in: “Data accessibility”
Holds: README, environment (SHRED_base/requirements.txt)
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (5 files), PyTorch (5 files), scikit-learn (3 files), h5py (1 file), imageio (1 file), Matplotlib (1 file), SciPy (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
6 files
At the source: github.com/meganebers/mobileSHRED

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;
  • 25 scripts, each with its path and the digest of its content;
  • 5 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 accessibility

All code and data are publicly available. Figures and experiments can be reproduced from the following GitHub: http://github.com/amysrude/neuralSHRED. The datasets used in the analyses are from the following sites. C. elegans: http://osf.io/2395t/. Mouse: download data: http://allensdk.readthedocs.io/en/latest/_static/examples/nb/ecephys_data_access.html. Documentation: http://allensdk.readthedocs.io/en/latest/visual_coding_neuropixels.html. Zebrafish: http://www.nature.com/articles/nmeth.2434. Biolocomotion: http://github.com/meganebers/mobileSHRED, http://simtk.org/projects/ankleexopred and [65,66,67].

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, 5 keywords, 9 MeSH terms, 2 funders, 58 references.

Cite

This paper

Rude, A., & Kutz, J. N. (2026). Shallow recurrent decoders for neural and behavioural dynamics. Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 381(1958), 20240461. https://doi.org/10.1098/rstb.2024.0461

BibTeX

@article{rude2026shallow,
author = {Rude, Amy and Kutz, J Nathan},
title = {{Shallow recurrent decoders for neural and behavioural dynamics}},
journal = {Philosophical transactions of the Royal Society of London. Series B, Biological sciences},
year = {2026},
month = sep,
volume = {381},
number = {1958},
pages = {20240461},
publisher = {Royal Society},
issn = {0962-8436},
doi = {10.1098/rstb.2024.0461},
url = {https://doi.org/10.1098/rstb.2024.0461},
pmid = {42750452},
pmcid = {PMC13583487}
}

RIS

TY - JOUR
AU - Rude, Amy
AU - Kutz, J Nathan
TI - Shallow recurrent decoders for neural and behavioural dynamics
T2 - Philosophical transactions of the Royal Society of London. Series B, Biological sciences
J2 - Philos Trans R Soc Lond B Biol Sci
PY - 2026
DA - 2026/09/01
VL - 381
IS - 1958
SP - 20240461
SN - 0962-8436
PB - Royal Society
DO - 10.1098/rstb.2024.0461
UR - https://doi.org/10.1098/rstb.2024.0461
LA - en
ER -

CSL-JSON

{
"id": "10.1098/rstb.2024.0461",
"type": "article-journal",
"title": "Shallow recurrent decoders for neural and behavioural dynamics",
"container-title": "Philosophical transactions of the Royal Society of London. Series B, Biological sciences",
"author": [
{
"family": "Rude",
"given": "Amy"
},
{
"family": "Kutz",
"given": "J Nathan"
}
],
"container-title-short": "Philos Trans R Soc Lond B Biol Sci",
"volume": "381",
"issue": "1958",
"page": "20240461",
"DOI": "10.1098/rstb.2024.0461",
"PMID": "42750452",
"PMCID": "PMC13583487",
"ISSN": "0962-8436",
"publisher": "Royal Society",
"URL": "https://doi.org/10.1098/rstb.2024.0461",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
1
]
]
}
}

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.1038/s41593-026-02257-5 [code]
Neural sequences underlying directed turning in Caenorhabditis elegans.
Journal: Nature neuroscience
In common: h5py, Pillow, PyTorch, 4 other tools, C. elegans, 3 references
[2] doi:10.1038/s41467-026-72709-w [code]
An epifluorescence microscope design for naturalistic behavior and cellular activity in freely moving Caenorhabditis elegans.
Journal: Nature communications
In common: imageio, h5py, Pillow, 5 other tools, C. elegans, 1 reference
[3] doi:10.3389/frai.2026.1771088 [code]
Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: imageio, h5py, Pillow, 6 other tools, 1 reference
[4] doi:10.1038/s41467-026-76939-w [code]
HIPPIE: a generative model for electrophysiological analysis across species, technologies, and modalities.
Journal: Nature communications
In common: AllenSDK, h5py, Pillow, 6 other tools, mouse
[5] doi:10.3389/fncom.2026.1786996 [code]
Schumann-anchored golden ratio organization of human neural oscillations.
Journal: Frontiers in computational neuroscience
In common: imageio, MNE-Python, h5py, 6 other tools
[6] 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, pandas, 3 other tools, computational, 3 references
[7] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: xarray, h5py, Pillow, 6 other tools, mouse
[8] doi:10.1093/bioinformatics/btag540 [code]
Deciphering spatial heterogeneity by multimodal spatial transcriptomics modelling with SpatialModal.
Journal: Bioinformatics (Oxford, England)
In common: imageio, h5py, Pillow, 6 other tools, mouse
[9] doi:10.1038/s41467-026-76098-y [code]
A single computational objective can produce specialization of streams in visual cortex.
Journal: Nature communications
In common: xarray, h5py, Pillow, 6 other tools
[10] doi:10.1371/journal.pcbi.1014152 [code]
Synchronization properties in C. elegans: Relating behavioral circuits to structural and functional neuronal connectivity.
Journal: PLoS computational biology
In common: pandas, SciPy, Matplotlib, 1 other tool, C. elegans, 4 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.