OSCR

Locomotion optimizes sensory representations through a computational principle shared by rodents and primates.

Code ↔ Paper

15 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 15 matches
  1. [1] § METHODS › Nonlinearity optimization ↔ fig2/optimize_all_gabor.ipynb, lines 655–702 · score 0.86 · mutual information, response entropy, Gaussian distribution, stimulus bin, response bin, optimal parameters
  2. [2] § METHODS › Optimization of inhibitory connections ↔ fig4/network_interaction_analysis.ipynb, lines 341–483 · score 0.83 · enforce positivity, Pearson correlation, loss function, network interaction, gradient, subtracted
  3. [3] § METHODS › Nonlinearity optimization ↔ fig2/compute_ori_tuning.py, lines 24–135 · score 0.74 · logistic nonlinearity, Gaussian distribution, response bin, optimal parameters, entropy, utility
  4. [4] § METHODS › Temporal dynamics ↔ fig3/temporal_filtering_analysis.ipynb, lines 1170–1208 · score 0.73 · inverse Fourier transform, white noise, autocorrelation function, temporal filters, spectra, power
  5. [5] § METHODS › Temporal filtering ↔ fig3/temporal_filtering_analysis.ipynb, lines 1170–1208 · score 0.69 · filtered spectra, white noise, temporal filtering, autocorrelation function, spectral, power
  6. [6] § METHODS › Minimal model of sensing during locomotion ↔ minimal_model/agent_sim.py, lines 74–129 · score 0.69 · starting position, minimal model, agent, circle, uniform, velocity
  7. [7] § METHODS › Orientation tuning curves ↔ fig2/optimize_all_gabor.ipynb, lines 1327–1389 · score 0.64 · stationary tuning curve, moving tuning curve, regressing, additive, Gabor filters, fit
  8. [8] § METHODS › Minimal model of sensing during locomotion ↔ minimal_model/agent_sim_analysis.ipynb, lines 78–177 · score 0.63 · starting position, minimal model, trajectories, velocity, radius, agent
  9. [9] § METHODS › Population coding fidelity ↔ fig2/decoding_error_decreases_gabor.ipynb, lines 415–443 · score 0.62 · decoding error, Gaussian filter, decoder, firing rate, smoothed, concatenated
  10. [10] § METHODS › Orientation tuning curves ↔ fig2/compute_ori_tuning.py, lines 177–226 · score 0.61 · stationary tuning curve, moving tuning curve, regressing, additive, Gabor filters, fit
  11. [11] § RESULTS › Modulation of temporal filtering ensures efficiency of sensory coding during locomotion ↔ fig3/temporal_filtering_analysis.ipynb, lines 228–257 · score 0.59 · frequency domain, temporal frequencies, temporal filter, linear, Figure 3
  12. [12] § METHODS › Analysis of data from freely moving mice ↔ supp2/visualize_data_for_anti_mod.ipynb, lines 48–142 · score 0.58 · visual angle, Gabor filter bank, cpd, resolution
  13. [13] § METHODS › Analysis of data from freely moving mice ↔ supp2/visualize_data_for_anti_mod_no_eye_movement_correction.ipynb, lines 35–128 · score 0.58 · visual angle, Gabor filter bank, cpd, resolution
  14. [14] § METHODS › Surround suppression ↔ fig4/network_interaction_analysis.ipynb, lines 1678–1720 · score 0.53 · circular mask, radii, gratings, radius, Gabor filters, fitting
  15. [15] § RESULTS › Stimulus statistics explain differential modulation of sensory coding by locomotion in rodents and primates ↔ fig5/optimize_all_gabor_species_comp.ipynb, lines 1229–1363 · score 0.51 · confidence interval, firing rate, foveal, species, optimized, modulated

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 · 2,590 lines · 123 KB · no license · 3 matches

  1. # %% [markdown]
  2. # # Load Data
  4. # %%
  5. '''
  6. This script optimizes logistic nonlinearities for a variety of Gaussian distribtuions and plots the optimal parameters.
  7. Author: Jonathan Gant
  8. Date: 29.08.2024
  9. '''
  11. import numpy as np
  12. import matplotlib.pyplot as plt
  13. from utilities import logistic_func, calc_MI, calc_entropy
  14. from tqdm import tqdm
  15. import h5py
  16. import os
  17. import bottleneck as bn
  19. # Set random seed for reproducibility
  20. np.random.seed(0)
  22. # load in the data
  23. # all_gabor_responses = h5py.File('../results/new_nat_videos_gabor_responses_full_res_z_score.h5', 'r')
  24. all_gabor_responses = h5py.File('../results/new_nat_videos_gabor_responses_full_res_z_score_more_freq.h5', 'r')
  25. # all_gabor_responses_low_freq = h5py.File('../results/new_nat_videos_gabor_responses_full_res_more_low_freq_z_score.h5', 'r')
  26. # all_gabor_responses_eye_movements = h5py.File('../results/new_nat_videos_gabor_responses_full_res_more_low_freq_z_score_eye_movements.h5', 'r')
  27. # all_gabor_responses_eye_movements_long = h5py.File('../results/new_nat_videos_gabor_responses_full_res_more_low_freq_z_score_eye_movements_2s_interval_stat_only.h5', 'r')
  29. # all_gabor_responses = h5py.File('../results/new_nat_vids_gabor_responses_full_res.h5', 'r')
  31. # video size
  32. resolution_height = 1080
  33. resolution_width = 1920
  35. # fov
  36. horizontal_fov = 92
  37. vertical_fov = 61
  39. # conversion factor of pixels to degrees
  40. horizontal_pixels_per_degree = resolution_width / horizontal_fov
  41. vertical_pixels_per_degree = resolution_height / vertical_fov
  43. # average of the conversion factors to the nearest integer
  44. pixels_per_degree = np.ceil((horizontal_pixels_per_degree + vertical_pixels_per_degree) / 2)
  45. print(pixels_per_degree)
  47. # data hyperparameters
  48. orientation_arr = all_gabor_responses['orientation_arr'][()]
  49. phase_arr = all_gabor_responses['phase_arr'][()]
  50. position_arr = all_gabor_responses['position_arr'][()]
  51. wavelength_arr = all_gabor_responses['wavelength_arr'][()]
  52. # wavelength_arr_low_freq = all_gabor_responses_low_freq['wavelength_arr'][()]
  53. freq_arr = pixels_per_degree / wavelength_arr
  54. # freq_arr_low_freq = pixels_per_degree / wavelength_arr_low_freq
  55. filter_size = (resolution_height, resolution_width)
  56. print(freq_arr)
  57. low_spatial_freq_idx = np.arange(0, 31)
  58. high_spatial_freq_idx = np.arange(35, 70)
  60. # %%
  61. plt.plot(all_gabor_responses['field']['stationary_1'][0, 0, 0, 4, :], color='tab:red')
  62. plt.plot(all_gabor_responses['field']['moving_1'][0, 0, 0, 4, :], color='tab:blue')
  64. # %% [markdown]
  65. # # Temporal filtering
  67. # %% [markdown]
  68. # ## Load and average PSD
  70. # %%
  71. # parse the data and compute fourier transform for each filter and video
  72. environments = ['field', 'forest', 'orchard', 'tall_grass', 'pond']
  74. num_videos = 10
  75. vid_length = 50*30
  76. low_spatial_freq_idx = np.arange(0, 31)
  77. # low_spatial_freq_idx = np.arange(0, 19)
  78. # high_spatial_freq_idx = np.arange(35, 70)
  79. # low_spatial_freq_idx = np.arange(35, 70)
  81. stationary_stim = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), vid_length))
  82. moving_stim = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), vid_length))
  84. # get the responses for each environment
  85. for i, env_key in enumerate(environments):
  86. stationary_count = 0
  87. moving_count = 0
  88. all_gabor_responses_env = all_gabor_responses[env_key]
  89. print(all_gabor_responses_env.keys())
  90. for vid_key in all_gabor_responses_env.keys():
  91. if 'stationary' in vid_key:
  92. stationary_stim[i, stationary_count, :, :, :, :, :] = all_gabor_responses_env[vid_key][()][:, :, low_spatial_freq_idx, :, :vid_length]
  93. stationary_count += 1
  94. if 'moving' in vid_key and 'free_moving' not in vid_key:
  95. moving_stim[i, moving_count, :, :, :, :, :] = all_gabor_responses_env[vid_key][()][:, :, low_spatial_freq_idx, :, :vid_length]
  96. moving_count += 1
  98. # %%
  99. def compute_psd_rfft(x, fs):
  100. """
  101. Compute one-sided PSD using rfft for a real-valued signal.
  103. Parameters:
  104. - x: 1D real-valued input signal
  105. - fs: sampling frequency (Hz)
  107. Returns:
  108. - freqs: frequency bins (Hz)
  109. - psd: power spectral density (power/Hz)
  110. """
  111. N = len(x)
  112. X = np.fft.rfft(x, axis=-1)
  113. X_mag2 = np.abs(X) ** 2
  115. psd = X_mag2 / (N * fs)
  117. # Multiply by 2 to account for negative freqs (except DC and Nyquist if N even)
  118. if N % 2 == 0:
  119. psd[1:-1] *= 2
  120. else:
  121. psd[1:] *= 2
  123. return psd
  125. # %%
  126. # compute the fourier transform for each filter and video using scipy
  127. from scipy.fft import rfftn, rfftfreq
  129. sampling_rate = 30
  131. stationary_stim_psd = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), vid_length//2+1))
  132. moving_stim_psd = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), vid_length//2+1))
  133. for i, env_key in enumerate(environments):
  134. for j in range(num_videos):
  135. stationary_stim_psd[i, j] = compute_psd_rfft(stationary_stim[i, j], sampling_rate)
  136. moving_stim_psd[i, j] = compute_psd_rfft(moving_stim[i, j], sampling_rate)
  138. # %%
  139. # compute the temporal frequencies
  140. sampling_rate = 30
  141. num_samples = vid_length
  142. # frequencies = fftfreq(num_samples, 1/sampling_rate)[1:vid_length//2]
  143. frequencies = rfftfreq(num_samples, 1/sampling_rate)
  145. # pick a cutoff frequency
  146. rf_window = vid_length/sampling_rate # .5 # seconds
  147. cutoff_freq = 1/rf_window
  149. # find the first index wheter the frequency is greater than the cutoff frequency
  150. cutoff_idx = np.where(frequencies >= cutoff_freq)[0][0]
  152. # %% [markdown]
  153. # ## Welch method
  155. # %%
  156. from scipy.signal import welch
  158. def compute_psd_welch(x, fs, window_size, overlap):
  159. """
  160. Compute one-sided PSD using Welch's method for a real-valued signal with variable window size and full overlap.
  162. Parameters:
  163. - x: 1D real-valued input signal
  164. - fs: sampling frequency (Hz)
  165. - window_size: length of each segment (samples)
  167. Returns:
  168. - freqs: frequency bins (Hz)
  169. - psd: power spectral density (power/Hz)
  170. """
  171. freqs, psd = welch(
  172. x,
  173. fs=fs,
  174. window='boxcar',
  175. nperseg=window_size,
  176. noverlap=overlap, # full overlap
  177. return_onesided=True,
  178. scaling='density',
  179. detrend='constant'
  180. )
  181. return freqs, psd
  183. # %%
  184. # compute the fourier transform for each filter and video using scipy
  185. from scipy.fft import rfftn, rfftfreq
  187. sampling_rate = 30
  188. window_size = int(5 * sampling_rate) # 5 seconds window size
  190. stationary_stim_psd = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), window_size//2+1))
  191. moving_stim_psd = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), window_size//2+1))
  192. for i, env_key in enumerate(environments):
  193. for j in range(num_videos):
  194. _, stationary_stim_psd[i, j] = compute_psd_welch(stationary_stim[i, j], sampling_rate, window_size, window_size//2)
  195. frequencies, moving_stim_psd[i, j] = compute_psd_welch(moving_stim[i, j], sampling_rate, window_size, window_size//2)
  197. # %%
  198. # compute the PSD as the square of the absolute value of the fourier transform, also truncate the DC component and the negative frequencies
  199. # set the DC to zero
  200. # stationary_stim_freq_no_dc = np.zeros_like(stationary_stim_psd)
  201. # stationary_stim_freq_no_dc[:, :, :, :, :, :, cutoff_idx:] = stationary_stim_psd[:, :, :, :, :, :, cutoff_idx:]
  202. # moving_stim_freq_no_dc = np.zeros_like(moving_stim_psd)
  203. # moving_stim_freq_no_dc[:, :, :, :, :, :, cutoff_idx:] = moving_stim_psd[:, :, :, :, :, :, cutoff_idx:]
  205. # stationary_stim_psd_avg = np.mean(stationary_stim_freq_no_dc, axis=(0, 1))
  206. # moving_stim_psd_avg = np.mean(moving_stim_freq_no_dc, axis=(0, 1))
  209. stationary_stim_psd_avg = np.mean(stationary_stim_psd, axis=(0, 1))
  210. moving_stim_psd_avg = np.mean(moving_stim_psd, axis=(0, 1))
  212. # %%
  213. ori_idx = 1
  214. phase_idx = 3
  215. freq_idx = 10
  216. pos_idx = 4
  218. plt.plot(frequencies[1:-1], moving_stim_psd_avg[ori_idx, phase_idx, freq_idx, pos_idx, 1:-1], label='Moving', color='tab:orange')
  219. plt.plot(frequencies[1:-1], stationary_stim_psd_avg[ori_idx, phase_idx, freq_idx, pos_idx, 1:-1], label='Stationary', color='tab:gray')
  220. # plt.plot(frequencies[cutoff_idx:], moving_stim_psd_avg[ori_idx, phase_idx, freq_idx, pos_idx, cutoff_idx:], label='Moving', color='tab:orange')
  221. # plt.plot(frequencies[cutoff_idx:], stationary_stim_psd_avg[ori_idx, phase_idx, freq_idx, pos_idx, cutoff_idx:], label='Stationary', color='tab:gray')
  222. plt.xlabel('Frequency (Hz)')
  223. plt.ylabel('Power')
  224. plt.yscale('log')
  226. plt.savefig(f'../manuscript_figures/fig3_psd_example_ori_{ori_idx}_phase_{phase_idx}_freq_{freq_idx}_pos_{pos_idx}.pdf', format='pdf', bbox_inches='tight')
  228. # %% [markdown]
  229. # ## Apply filter
  231. # %%
  232. # definte the temporal filter in the frequency domain with a linear dropoff at a specific frequency determined by a free parameter theta. Define the filter as a function of a set of temporal frequencies
  233. def temporal_filter(temporal_freqs, theta, low_freq_cutoff=2):
  234. # create the filter
  235. filter = np.zeros(temporal_freqs.shape)
  236. # set the filter to 1 for frequencies below theta
  237. filter[temporal_freqs <= theta] = 1
  238. # filter[temporal_freqs < low_freq_cutoff] = 1 - (low_freq_cutoff - temporal_freqs[temporal_freqs < low_freq_cutoff])
  239. # create a linear dropoff for frequencies above theta with a fixed slope
  240. filter[temporal_freqs > theta] = 1 - (temporal_freqs[temporal_freqs > theta] - theta)
  241. # ensure that the filter is non-negative
  242. filter[filter < 0] = 0
  243. return np.square(filter)
  245. def temporal_filter_gauss(temporal_freqs, theta_u, theta_s):
  246. # create the filter as a Gaussian in the frequency domain with mean theta_u and standard deviation theta_s
  247. filter = np.exp(-0.5 * ((temporal_freqs - theta_u) / theta_s)**2)
  248. # ensure that the filter is non-negative
  249. filter[filter < 0] = 0
  250. return np.square(filter)
  252. def temporal_filter_exp(temporal_freq, theta):
  253. # create the filter as an exponential in the frequency domain with mean theta_u and standard deviation theta_s
  254. filter = np.exp(-temporal_freq/theta)
  255. # ensure that the filter is non-negative
  256. filter[filter < 0] = 0
  257. return np.square(filter)
  259. # %%
  260. # take the inverse fourier transform of the sqrt of the temporal filter
  261. from scipy.fft import irfftn, ifftshift
  262. # temporal_filter_freq = temporal_filter_gauss(frequencies, 0, 5)
  263. temporal_filter_freq = temporal_filter(frequencies, 6)
  264. temporal_filter_freq[0] = 0
  265. # temporal_filter_freq = temporal_filter(freq, 3)
  266. # temporal_filter_freq = temporal_filter_exp(frequencies, 5)
  268. inverse_temporal_filter = irfftn(np.sqrt(temporal_filter_freq), axes=-1)[:len(frequencies)]
  271. # %%
  272. # plot the filter in the frequency domain
  273. plt.figure(figsize=(8, 6))
  274. plt.plot(frequencies[1:], np.sqrt(temporal_filter_freq[1:]), label='Temporal filter', color='tab:blue')
  275. plt.xlabel('Frequency (Hz)')
  276. plt.ylabel('Filter amplitude')
  277. plt.savefig('../manuscript_figures/fig3_temporal_filter_freq.pdf', format='pdf', bbox_inches='tight')
  279. # %%
  280. plt.plot(np.linspace(0, 5000, len(frequencies)), inverse_temporal_filter)
  281. plt.axhline(0, color='k', linestyle='--')
  282. plt.xlim(0, 2000)
  284. # %%
  285. # simulate white noise
  286. noise_level = 100
  287. # noise = np.ones(frequencies.shape) * noise_level
  288. noise = np.ones(freq.shape) * noise_level
  290. # %%
  291. # repeat the process again for a range of theta values
  292. num_theta_values = 100
  293. array_shape = moving_stim_psd_avg.shape
  294. theta_values = np.linspace(2, 15, num_theta_values)
  295. moving_stim_information = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  296. stationary_stim_information = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  297. moving_stim_energy = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  298. stationary_stim_energy = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  299. moving_stim_energy_linear = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  300. stationary_stim_energy_linear = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  301. moving_stim_energy_quadratic = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  302. stationary_stim_energy_quadratic = np.zeros((num_theta_values, array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  303. for i in tqdm(range(num_theta_values)):
  304. # compute the filter
  305. filter = temporal_filter(frequencies, theta_values[i])
  306. # compute the filtered spectra
  307. filtered_stationary_stim_psd_avg = np.multiply(stationary_stim_psd_avg, filter)
  308. filtered_moving_stim_psd_avg = np.multiply(moving_stim_psd_avg, filter)
  309. # compute the information
  310. moving_stim_information[i] = np.sum(np.log2(1 + filtered_moving_stim_psd_avg / noise), axis=-1)
  311. stationary_stim_information[i] = np.sum(np.log2(1 + filtered_stationary_stim_psd_avg / noise), axis=-1)
  312. # compute the energy
  313. moving_stim_energy[i] = np.sum(filtered_moving_stim_psd_avg, axis=-1)
  314. stationary_stim_energy[i] = np.sum(filtered_stationary_stim_psd_avg, axis=-1)
  316. # compute the energy weighted by the square of the frequencies
  317. moving_stim_energy_quadratic[i] = np.sum(filtered_moving_stim_psd_avg*frequencies**2, axis=-1)
  318. stationary_stim_energy_quadratic[i] = np.sum(filtered_stationary_stim_psd_avg*frequencies**2, axis=-1)
  320. # compute the energy linearly weighted by frequency
  321. moving_stim_energy_linear[i] = np.sum(filtered_moving_stim_psd_avg*frequencies, axis=-1)
  322. stationary_stim_energy_linear[i] = np.sum(filtered_stationary_stim_psd_avg*frequencies, axis=-1)
  325. # %%
  326. # plot an example of the psd before and after filtering
  327. fig, ax = plt.subplots(1, 2, figsize=(16, 6), sharey=True)
  329. filter = temporal_filter(frequencies, 6)
  331. ori_idx = 0
  332. phase_idx = 0
  333. wavelength_idx = 8
  334. position_idx = 4
  336. ax[0].plot(frequencies, moving_stim_psd_avg[ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:blue')
  337. ax[0].plot(frequencies, stationary_stim_psd_avg[ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Stationary', color='tab:red')
  338. ax[1].plot(frequencies, np.multiply(moving_stim_psd_avg[ori_idx, phase_idx, wavelength_idx, position_idx, :], filter), label='Moving', color='tab:blue')
  339. ax[1].plot(frequencies, np.multiply(stationary_stim_psd_avg[ori_idx, phase_idx, wavelength_idx, position_idx, :], filter), label='Stationary', color='tab:red')
  340. ax[1].set_yscale('log')
  342. # %%
  343. moving_test_informations = moving_stim_information[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  344. stationary_test_informations = stationary_stim_information[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  345. moving_test_energy = moving_stim_energy[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  346. stationary_test_energy = stationary_stim_energy[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  348. # create a subplot with 1 row and 3 columns. First column plot the information vs. theta, second column plot the energy vs. theta, third column plot the information vs. energy
  349. fig, axs = plt.subplots(1, 3, figsize=(15, 5))
  350. axs[0].plot(theta_values, moving_test_informations, color='tab:blue', label='Moving test')
  351. axs[0].plot(theta_values, stationary_test_informations, color='tab:red', label='Stationary test')
  352. axs[0].set_xlabel('Theta')
  353. axs[0].set_ylabel('Information')
  354. axs[0].set_title('Information vs. theta')
  355. axs[0].legend()
  356. axs[1].plot(theta_values, moving_test_energy, color='tab:blue', label='Moving test')
  357. axs[1].plot(theta_values, stationary_test_energy, color='tab:red', label='Stationary test')
  358. axs[1].set_xlabel('Theta')
  359. axs[1].set_ylabel('Energy')
  360. axs[1].set_title('Energy vs. theta')
  361. axs[1].legend()
  362. axs[2].plot(moving_test_energy, moving_test_informations, color='tab:blue', label='Moving test')
  363. axs[2].plot(stationary_test_energy, stationary_test_informations, color='tab:red', label='Stationary test')
  364. axs[2].set_xlabel('Energy')
  365. axs[2].set_ylabel('Information')
  366. axs[2].set_title('Information vs. energy')
  367. axs[2].legend()
  368. plt.tight_layout()
  369. plt.show()
  371. # %%
  372. lmbd_arr = np.logspace(-5, -3, 16)
  374. fig, ax = plt.subplots(4, 4, figsize=(16, 16))
  376. for i, lmbd in enumerate(lmbd_arr):
  377. moving_utility = moving_test_informations - lmbd*moving_test_energy
  378. stationary_utility = stationary_test_informations - lmbd*stationary_test_energy
  379. ax[i//4, i%4].plot(theta_values, moving_utility, color='tab:blue', label='Moving test')
  380. ax[i//4, i%4].plot(theta_values, stationary_utility, color='tab:red', label='Stationary test')
  381. ax[i//4, i%4].set_title('Lambda = ' + str(np.round(lmbd, 7)))
  383. # %%
  384. # for each filter, for a given energy find the closest value of theta which has that energy
  385. energy_arr = np.logspace(6, 11, 1000)
  387. optimal_theta_val_mov = np.zeros((len(energy_arr), len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  388. optimal_theta_val_stat = np.zeros((len(energy_arr), len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  389. optimal_theta_val_mov_linear = np.zeros((len(energy_arr), len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  390. optimal_theta_val_stat_linear = np.zeros((len(energy_arr), len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  391. optimal_theta_val_mov_quadratic = np.zeros((len(energy_arr), len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  392. optimal_theta_val_stat_quadratic = np.zeros((len(energy_arr), len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  394. # loop over orientation, phase, wavelength, position
  395. for ori_idx in range(len(orientation_arr)):
  396. for phase_idx in range(len(phase_arr)):
  397. for wavelength_idx in range(len(wavelength_arr[low_spatial_freq_idx])):
  398. for position_idx in range(len(position_arr)):
  399. test_moving_energy = moving_stim_energy[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  400. test_stationary_energy = stationary_stim_energy[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  401. # interpolate the energy values to get a smooth curve
  402. test_moving_theta_interp = np.interp(energy_arr, test_moving_energy, theta_values)
  403. test_stationary_theta_interp = np.interp(energy_arr, test_stationary_energy, theta_values)
  404. optimal_theta_val_mov[:, ori_idx, phase_idx, wavelength_idx, position_idx] = test_moving_theta_interp
  405. optimal_theta_val_stat[:, ori_idx, phase_idx, wavelength_idx, position_idx] = test_stationary_theta_interp
  407. # repeat for the linear energy
  408. test_moving_energy_linear = moving_stim_energy_linear[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  409. test_stationary_energy_linear = stationary_stim_energy_linear[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  410. # interpolate the energy values to get a smooth curve
  411. test_moving_theta_interp_linear = np.interp(energy_arr, test_moving_energy_linear, theta_values)
  412. test_stationary_theta_interp_linear = np.interp(energy_arr, test_stationary_energy_linear, theta_values)
  413. optimal_theta_val_mov_linear[:, ori_idx, phase_idx, wavelength_idx, position_idx] = test_moving_theta_interp_linear
  414. optimal_theta_val_stat_linear[:, ori_idx, phase_idx, wavelength_idx, position_idx] = test_stationary_theta_interp_linear
  416. # repeat for the quadratic energy
  417. test_moving_energy_quadratic = moving_stim_energy_quadratic[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  418. test_stationary_energy_quadratic = stationary_stim_energy_quadratic[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  419. # interpolate the energy values to get a smooth curve
  420. test_moving_theta_interp_quadratic = np.interp(energy_arr, test_moving_energy_quadratic, theta_values)
  421. test_stationary_theta_interp_quadratic = np.interp(energy_arr, test_stationary_energy_quadratic, theta_values)
  422. optimal_theta_val_mov_quadratic[:, ori_idx, phase_idx, wavelength_idx, position_idx] = test_moving_theta_interp_quadratic
  423. optimal_theta_val_stat_quadratic[:, ori_idx, phase_idx, wavelength_idx, position_idx] = test_stationary_theta_interp_quadratic
  425. # %%
  426. ori_idx = 0
  427. phase_idx = 0
  428. wavelength_idx = -1
  429. position_idx = 4
  431. # plot the difference in theta values vs. energy
  433. plt.plot(energy_arr, optimal_theta_val_mov[:, ori_idx, phase_idx, wavelength_idx, position_idx]-optimal_theta_val_stat[:, ori_idx, phase_idx, wavelength_idx, position_idx], color='tab:purple')
  435. # %%
  436. optimal_theta_diff = optimal_theta_val_mov - optimal_theta_val_stat
  437. optimal_theta_diff_avg = np.mean(optimal_theta_diff, axis=(1, 2, 4))
  438. optimal_theta_diff_linear = optimal_theta_val_mov_linear - optimal_theta_val_stat_linear
  439. optimal_theta_diff_linear_avg = np.mean(optimal_theta_diff_linear, axis=(1, 2, 4))
  440. optimal_theta_diff_quadratic = optimal_theta_val_mov_quadratic - optimal_theta_val_stat_quadratic
  441. optimal_theta_diff_quadratic_avg = np.mean(optimal_theta_diff_quadratic, axis=(1, 2, 4))
  443. # %%
  444. # plot the average difference in optimal theta, color them with viridis
  446. from matplotlib import cm
  447. from matplotlib.colors import Normalize
  448. from mpl_toolkits.axes_grid1.inset_locator import inset_axes
  450. colors_map = cm.get_cmap('viridis', len(wavelength_arr[low_spatial_freq_idx]))
  452. fig, ax = plt.subplots(1, 3, figsize=(4*3*1.5, 3*1.5), sharey=True)
  453. for i in range(wavelength_arr[low_spatial_freq_idx].shape[0]):
  454. ax[0].plot(energy_arr, optimal_theta_diff_avg[:, i], color=colors_map(i))
  455. ax[1].plot(energy_arr, optimal_theta_diff_linear_avg[:, i], color=colors_map(i))
  456. ax[2].plot(energy_arr, optimal_theta_diff_quadratic_avg[:, i], color=colors_map(i))
  457. ax[0].set_xscale('log')
  458. ax[1].set_xscale('log')
  459. ax[2].set_xscale('log')
  460. ax[0].set_title('Unweighted energy')
  461. ax[1].set_title('Linearly weighted energy')
  462. ax[2].set_title('Quadratically weighted energy')
  463. ax[0].set_xlabel('Energy constraint')
  464. ax[1].set_xlabel('Energy constraint')
  465. ax[2].set_xlabel('Energy constraint')
  466. ax[0].set_ylabel('$\\Delta\\theta$ [moving - stationary]')
  467. # add something to indicate that the darker the color, the lower the spatial frequency
  468. # add a colorbar
  469. norm = Normalize(vmin=min(freq_arr[low_spatial_freq_idx]), vmax=max(freq_arr[low_spatial_freq_idx]))
  470. sm = plt.cm.ScalarMappable(cmap=colors_map, norm=norm)
  471. sm.set_array([])
  474. axins1 = inset_axes(
  475. ax[0],
  476. width="25%", # width: 50% of parent_bbox width
  477. height="5%", # height: 5%
  478. loc="upper right",
  479. )
  480. axins1.xaxis.set_ticks_position("bottom")
  481. cbar = fig.colorbar(sm, cax=axins1, orientation="horizontal")
  482. cbar.set_label('spatial freq.\n[cycles/degree]')
  484. axins2 = inset_axes(
  485. ax[1],
  486. width="25%", # width: 50% of parent_bbox width
  487. height="5%", # height: 5%
  488. loc="upper right",
  489. )
  490. axins2.xaxis.set_ticks_position("bottom")
  491. cbar = fig.colorbar(sm, cax=axins2, orientation="horizontal")
  492. cbar.set_label('spatial freq.\n[cycles/degree]')
  493. axins2 = inset_axes(
  494. ax[2],
  495. width="25%", # width: 50% of parent_bbox width
  496. height="5%", # height: 5%
  497. loc="upper right",
  498. )
  499. axins2.xaxis.set_ticks_position("bottom")
  500. cbar = fig.colorbar(sm, cax=axins2, orientation="horizontal")
  501. cbar.set_label('spatial freq.\n[cycles/degree]')
  503. # %%
  504. ori_idx = 0
  505. phase_idx = 0
  506. position_idx = 4
  508. wavelength_idices = [0, 15, -1]
  510. moving_test_energy = moving_stim_energy[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  511. stationary_test_energy = stationary_stim_energy[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  513. # create a subplot with 1 row and 3 columns. First column plot the information vs. theta, second column plot the energy vs. theta, third column plot the information vs. energy
  514. fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharex=True)
  515. axs[0].plot(theta_values, moving_stim_energy[:, ori_idx, phase_idx, wavelength_idices[0], position_idx], color='tab:blue', label='Moving test')
  516. axs[0].plot(theta_values, stationary_stim_energy[:, ori_idx, phase_idx, wavelength_idices[0], position_idx], color='tab:red', label='Stationary test')
  517. axs[0].set_xlabel('Theta')
  518. axs[0].set_ylabel('Energy')
  519. axs[0].set_title('Low frequency')
  520. axs[0].legend()
  521. axs[1].plot(theta_values, moving_stim_energy[:, ori_idx, phase_idx, wavelength_idices[1], position_idx], color='tab:blue', label='Moving test')
  522. axs[1].plot(theta_values, stationary_stim_energy[:, ori_idx, phase_idx, wavelength_idices[1], position_idx], color='tab:red', label='Stationary test')
  523. axs[1].set_xlabel('Theta')
  524. axs[1].set_ylabel('Energy')
  525. axs[1].set_title('Mid frequency')
  526. axs[1].legend()
  527. axs[2].plot(theta_values, moving_stim_energy[:, ori_idx, phase_idx, wavelength_idices[2], position_idx], color='tab:blue', label='Moving test')
  528. axs[2].plot(theta_values, stationary_stim_energy[:, ori_idx, phase_idx, wavelength_idices[2], position_idx], color='tab:red', label='Stationary test')
  529. axs[2].set_xlabel('Theta')
  530. axs[2].set_ylabel('Energy')
  531. axs[2].set_title('High frequency')
  532. axs[2].legend()
  533. plt.tight_layout()
  534. plt.show()
  536. # %%
  537. # pick a contraint and then see how a stimulus filtered by the optimal filter looks
  538. constraint_energy = 1e9
  539. constraint_idx = np.argmin(np.abs(energy_arr-constraint_energy))
  541. optimal_theta_mov = optimal_theta_val_mov[constraint_idx, :, :, :, :]
  542. optimal_theta_stat = optimal_theta_val_stat[constraint_idx, :, :, :, :]
  544. ori_idx = 0
  545. phase_idx = 0
  546. wavelength_idx = 3
  547. position_idx = 4
  549. theta_mov = optimal_theta_mov[ori_idx, phase_idx, wavelength_idx, position_idx]
  550. theta_stat = optimal_theta_stat[ori_idx, phase_idx, wavelength_idx, position_idx]
  552. # compute the filter
  553. filter_mov = temporal_filter(frequencies, theta_mov)
  554. filter_stat = temporal_filter(frequencies, theta_stat)
  555. # compute the filtered spectra
  556. filtered_moving_stim_psd = np.multiply(noise, filter_mov)
  557. filtered_stationary_stim_psd = np.multiply(noise, filter_stat)
  558. # compute the energy
  559. test_moving_stim_energy = np.sum(filtered_moving_stim_psd, axis=-1)
  560. test_stationary_stim_energy = np.sum(filtered_stationary_stim_psd, axis=-1)
  562. # compute the energy weighted by the square of the frequencies
  563. test_moving_stim_energy_quadratic = np.sum(filtered_moving_stim_psd*frequencies**2, axis=-1)
  564. test_stationary_stim_energy_quadratic = np.sum(filtered_stationary_stim_psd*frequencies**2, axis=-1)
  566. # compute the energy linearly weighted by frequency
  567. test_moving_stim_energy_linear = np.sum(filtered_moving_stim_psd*frequencies, axis=-1)
  568. test_stationary_stim_energy_linear = np.sum(filtered_stationary_stim_psd*frequencies, axis=-1)
  571. # %%
  572. # create a figure which shows the spectra before filtering, the filter, and the spectra after filtering
  573. fig, ax = plt.subplots(1, 3, figsize=(15, 5))
  574. ax[0].plot(frequencies, noise, label='Noise', color='tab:gray')
  575. ax[1].plot(frequencies, filter_mov, label='Moving', color='tab:blue')
  576. ax[1].plot(frequencies, filter_stat, label='Stationary', color='tab:red')
  577. ax[2].plot(frequencies, filtered_moving_stim_psd, label='Moving', color='tab:blue')
  578. ax[2].plot(frequencies, filtered_stationary_stim_psd, label='Stationary', color='tab:red')
  580. # %%
  581. # compute the autocorrelation function by taking the inverse fourier transform of the PSD
  582. from scipy.fft import irfft
  583. # compute the autocorrelation function
  584. moving_autocorr = irfft(filtered_moving_stim_psd)[:len(filtered_moving_stim_psd)//2]
  585. stationary_autocorr = irfft(filtered_stationary_stim_psd)[:len(filtered_stationary_stim_psd)//2]
  587. # normalize the autocorrelation function
  588. moving_autocorr /= np.max(moving_autocorr)
  589. stationary_autocorr /= np.max(stationary_autocorr)
  591. time = np.arange(0, len(moving_autocorr)) / 30
  593. # plot the autocorrelation function
  594. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  595. ax.plot(time, moving_autocorr, label='Moving', color='tab:blue')
  596. ax.plot(time, stationary_autocorr, label='Stationary', color='tab:red')
  597. ax.set_xlabel('Time (s)')
  598. ax.set_ylabel('Autocorrelation')
  599. ax.set_title('Autocorrelation function')
  600. ax.legend()
  601. ax.set_xlim(0,2)
  602. plt.show()
  604. # %% [markdown]
  605. # ## Heatmap energy landscape and theta
  607. # %% [markdown]
  608. # ### Average over filters
  610. # %%
  611. moving_stim_psd.shape
  613. # %%
  614. # compute the average spectrum over all filters
  615. mean_moving_psd = np.mean(moving_stim_psd, axis=(0, 1, 2, 3, 4, 5))
  616. mean_stationary_psd = np.mean(stationary_stim_psd, axis=(0, 1, 2, 3, 4, 5))
  618. # plot the mean spectrum
  619. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  620. ax.plot(frequencies[1:-1], mean_stationary_psd[1:-1], label='Stationary', color='tab:gray')
  621. ax.plot(frequencies[1:-1], mean_moving_psd[1:-1], label='Moving', color='tab:orange')
  622. ax.set_xlabel('Frequency (Hz)')
  623. ax.set_ylabel('Power')
  624. ax.set_title('Mean power spectrum')
  625. ax.set_yscale('log')
  626. ax.legend()
  628. plt.savefig(f'../manuscript_figures/fig3_psd_mean_full.pdf', format='pdf', bbox_inches='tight')
  630. # %%
  631. # normalize the spectrum such that the are under the curve is 1
  632. freq_spacing = np.diff(frequencies)[0]
  633. norm_factor = np.sum(mean_stationary_psd*freq_spacing)
  634. mean_stationary_psd /= norm_factor
  635. mean_moving_psd /= norm_factor
  637. # %%
  638. np.cumsum(mean_stationary_psd*freq_spacing)
  640. # %%
  641. np.where(np.cumsum(mean_moving_psd*freq_spacing) <= auc_val)
  643. # %%
  644. # compute the average spectrum over all filters
  645. # plot the mean spectrum
  646. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  647. ax.plot(frequencies[1:-1], mean_moving_psd[1:-1], label='Moving', color='tab:orange')
  648. ax.plot(frequencies[1:-1], mean_stationary_psd[1:-1], label='Stationary', color='tab:gray')
  649. ax.set_xlabel('Frequency (Hz)')
  650. ax.set_ylabel('Power')
  651. ax.set_title('Mean power spectrum')
  652. ax.set_yscale('log')
  653. # ax.set_xscale('log')
  654. ax.legend()
  656. # find where the area under the stationary psd is a certain value
  657. # auc_val = 1
  658. # aoc_idx_stat = np.where(np.cumsum(mean_stationary_psd*freq_spacing) <= auc_val)[0][-1]
  659. # # find the frequency at that index
  660. # freq_val_stat = frequencies[aoc_idx_stat]
  661. # auc_idx_mov = np.where(np.cumsum(mean_moving_psd*freq_spacing) <= auc_val)[0][-1]
  662. # # find the frequency at that index
  663. # freq_val_mov = frequencies[auc_idx_mov]
  665. auc_val = .95
  666. cum_stationary = np.cumsum(mean_stationary_psd * freq_spacing)
  667. cum_moving = np.cumsum(mean_moving_psd * freq_spacing)
  669. # Interpolate to find the frequency where the cumulative sum reaches auc_val
  670. freq_val_stat = np.interp(auc_val, cum_stationary, frequencies)
  671. freq_val_mov = np.interp(auc_val, cum_moving, frequencies)
  673. # # plot filled area under the curve
  674. # ax.fill_between(frequencies[1:-1], mean_stationary_psd[1:-1], where=(frequencies[1:-1] <= freq_val_stat), color='tab:gray', alpha=0.5)
  675. # ax.fill_between(frequencies[1:-1], mean_moving_psd[1:-1], where=(frequencies[1:-1] <= freq_val_mov), color='tab:orange', alpha=0.5)
  677. # Interpolate the PSD at the cutoff frequencies
  678. psd_val_stat = np.interp(freq_val_stat, frequencies, mean_stationary_psd)
  679. psd_val_mov = np.interp(freq_val_mov, frequencies, mean_moving_psd)
  681. # For fill_between, extend the frequency and PSD arrays to include the interpolated cutoff point
  682. def extend_for_fill(frequencies, psd, cutoff_freq, cutoff_psd):
  683. mask = frequencies <= cutoff_freq
  684. # Find the last index before the cutoff
  685. last_idx = np.where(mask)[0][-1]
  686. # Insert the cutoff point after last_idx
  687. new_freqs = np.insert(frequencies[mask], last_idx + 1, cutoff_freq)
  688. new_psd = np.insert(psd[mask], last_idx + 1, cutoff_psd)
  689. return new_freqs, new_psd
  691. freqs_stat_fill, psd_stat_fill = extend_for_fill(frequencies[1:-1], mean_stationary_psd[1:-1], freq_val_stat, psd_val_stat)
  692. freqs_mov_fill, psd_mov_fill = extend_for_fill(frequencies[1:-1], mean_moving_psd[1:-1], freq_val_mov, psd_val_mov)
  694. # plot filled area under the curve using the interpolated cutoff
  695. ax.fill_between(freqs_stat_fill, psd_stat_fill, color='tab:gray', alpha=0.5)
  696. ax.fill_between(freqs_mov_fill, psd_mov_fill, color='tab:orange', alpha=0.5)
  699. plt.savefig(f'../manuscript_figures/fig3_psd_mean_normalize_2Hz_crop_auc_{auc_val}.pdf', format='pdf', bbox_inches='tight')
  701. # %%
  702. psd_stat_fill
  704. # %%
  705. # plot the energy vs. theta
  706. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  707. ax.plot(frequencies[1:], np.cumsum(mean_stationary_psd*freq_spacing)[1:], label='Stationary', color='tab:gray')
  708. ax.plot(frequencies[1:], np.cumsum(mean_moving_psd*freq_spacing)[1:], label='Moving', color='tab:orange')
  709. ax.axhline(y=auc_val, color='k', linestyle='--')
  710. ax.set_xlabel('Frequency (Hz)')
  711. ax.set_ylabel('Energy')
  712. ax.set_title('Energy vs. theta')
  713. ax.set_yscale('log')
  714. ax.set_xscale('log')
  715. ax.legend()
  716. plt.savefig(f'../manuscript_figures/fig3_energy_vs_freq_aoc_{auc_val}.pdf', format='pdf', bbox_inches='tight')
  718. # %%
  719. # repeat the process again for a range of theta values
  720. num_theta_values = 1000
  721. theta_values = np.linspace(0, 15, num_theta_values)
  722. moving_stim_energy = np.zeros(num_theta_values)
  723. stationary_stim_energy = np.zeros(num_theta_values)
  724. for i in tqdm(range(num_theta_values)):
  725. # compute the filter
  726. filter = temporal_filter(frequencies, theta_values[i])
  727. # compute the filtered spectra
  728. filtered_mean_stationary_psd = np.multiply(mean_stationary_psd, filter)
  729. filtered_mean_moving_psd = np.multiply(mean_moving_psd, filter)
  730. # compute the energy
  731. moving_stim_energy[i] = np.sum(filtered_mean_moving_psd*freq_spacing, axis=-1)
  732. stationary_stim_energy[i] = np.sum(filtered_mean_stationary_psd*freq_spacing, axis=-1)
  734. # %%
  735. # plot the energy vs. theta
  736. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  737. ax.plot(theta_values[:], moving_stim_energy[:], label='Moving', color='tab:orange')
  738. ax.plot(theta_values[:], stationary_stim_energy[:], label='Stationary', color='tab:gray')
  739. ax.set_xlabel('Theta')
  740. ax.set_ylabel('Energy')
  741. ax.set_title('Energy vs. theta')
  742. ax.set_yscale('log')
  743. ax.set_xscale('log')
  744. ax.legend()
  746. # %%
  747. relative = True
  749. # for given energy capacities in the stationary case (between 0 and 1)
  750. energy_capacity_stat = np.linspace(0, 1, 101)
  751. # energy_capacity_mov = np.linspace(0, 2, 101)
  752. delta_energy_capacity = np.linspace(0.01, 10, 1000)
  753. # energy_capacity_stat = np.logspace(-1, 0, 100)
  754. # delta_energy_capacity = np.logspace(-1, 1, 101)
  756. # interpolate the energy values to get a smooth curve
  757. stationary_theta_interp = np.interp(energy_capacity_stat, stationary_stim_energy, theta_values)
  759. if relative:
  760. energy_capacity_mov = np.outer(delta_energy_capacity, energy_capacity_stat) + energy_capacity_stat[np.newaxis, :] # relative
  761. delta_energy_capacity *= 100 # convert to percent
  762. else:
  763. energy_capacity_mov = delta_energy_capacity[:, np.newaxis] + energy_capacity_stat[np.newaxis, :] # absolute
  765. moving_theta_interp = np.interp(energy_capacity_mov, moving_stim_energy, theta_values)
  768. # %%
  769. import matplotlib.colors as mcolors
  770. from matplotlib import cm
  772. norm = mcolors.Normalize(vmin=np.min(theta_values), vmax=np.max(theta_values))
  773. two_color_norm = mcolors.TwoSlopeNorm(vmin=-5, vcenter=0, vmax=5)
  775. # example_points = [(50, 50), (75, 75), (95, 100)]
  776. example_points = [(75, 9), (90, 19), (95, 29), (99, 49)]
  778. fig, ax = plt.subplots(1, 3, figsize=(16, 5), width_ratios=[5, 5, 5], sharey=True)
  779. im0 = ax[0].imshow(np.repeat(stationary_theta_interp[np.newaxis, :], len(delta_energy_capacity), axis=0), origin='lower', norm=norm, cmap='gray_r', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  780. ax[0].set_title('Stationary theta')
  781. im1 = ax[1].imshow(moving_theta_interp, origin='lower', norm=norm, cmap='gray_r', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  782. ax[1].set_title('Moving theta')
  783. colors_map = mcolors.LinearSegmentedColormap.from_list('custom_cmap', [(1, 127/255, 14/255), (1, 1, 1), (0.5, 0.5, 0.5)], N=256)
  784. im2 = ax[2].imshow(stationary_theta_interp[np.newaxis, :]-moving_theta_interp, origin='lower', norm=two_color_norm, cmap=colors_map, extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  785. ax[2].set_title('Difference in theta [stationary - moving]')
  786. # ax[2].set_xlim(0.5, 1)
  787. # ax[2].set_ylim(0.5, 1)
  788. # add axis labels
  789. if relative:
  790. ax[0].set_ylabel('Percentage increase in energy capacity')
  791. else:
  792. ax[0].set_ylabel('Increase in energy capacity')
  794. fig.supxlabel('Energy capacity')
  796. # plot the example points in three different colors which are distinct from red, white and blue
  797. colors_map = cm.get_cmap('Blues', len(example_points)+2)
  798. # plot the example points in the first subplot
  799. for i, point in enumerate(example_points):
  800. ax[0].scatter(energy_capacity_stat[point[0]], delta_energy_capacity[point[1]], color=colors_map(i+1), s=100, edgecolor='black', label='Example point ' + str(i+1))
  801. ax[1].scatter(energy_capacity_stat[point[0]], delta_energy_capacity[point[1]], color=colors_map(i+1), s=100, edgecolor='black')
  802. ax[2].scatter(energy_capacity_stat[point[0]], delta_energy_capacity[point[1]], color=colors_map(i+1), s=100, edgecolor='black')
  804. ax[0].set_yscale('log')
  806. # add a colorbar without disrupting the layout of the other subplots
  807. cbar = fig.colorbar(im2, ax=ax[2], orientation='vertical', pad=0.04)
  808. cbar = fig.colorbar(im1, ax=ax[1], orientation='vertical', pad=0.04)
  809. cbar = fig.colorbar(im0, ax=ax[0], orientation='vertical', pad=0.04)
  810. plt.savefig(f'../manuscript_figures/fig3_diff_theta_stat_energy_vs_delta_energy.pdf', format='pdf', bbox_inches='tight')
  812. # %%
  813. stationary_theta_interp[example_points[-1][0]]
  815. # %%
  816. example_filter_mov.shape
  818. # %%
  819. from scipy.fft import irfft
  821. # compute the autocorrelation function for the example points by extracting the theta values in the stationary and moving conditions and then applying the filter to white noise and then computing the autocorrelation function as the inverse fourier transform of the PSD
  822. alpha_arr = [0.25, 0.5, 0.75, 1]
  823. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  824. for i in range(len(example_points)):
  825. # compute the filter
  826. example_theta_mov = moving_theta_interp[example_points[i][1], example_points[i][0]]
  827. example_theta_stat = stationary_theta_interp[example_points[i][0]]
  828. # compute the filter
  829. example_filter_mov = temporal_filter(frequencies, example_theta_mov)
  830. example_filter_stat = temporal_filter(frequencies, example_theta_stat)
  831. # compute the inverse fourier transform
  832. example_moving_autocorr = irfft(example_filter_mov/freq_spacing)
  833. example_stationary_autocorr = irfft(example_filter_stat/freq_spacing)
  834. # normalize the autocorrelation function
  835. example_moving_autocorr /= np.max(example_moving_autocorr)
  836. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  837. # plot the autocorrelation function
  838. time = np.arange(0, len(example_moving_autocorr)) / 30*1000
  839. ax.plot(time, example_moving_autocorr, label='Moving', color='tab:orange', alpha=alpha_arr[i])
  840. ax.plot(time, example_stationary_autocorr, label='Stationary', color='tab:gray', alpha=alpha_arr[i])
  841. ax.set_xlabel('Time (ms)')
  842. ax.set_ylabel('Autocorrelation')
  843. # ax.legend()
  844. ax.set_xlim(0, 500)
  845. plt.savefig(f'../manuscript_figures/fig3_autocorr_example_{i}.pdf', format='pdf', bbox_inches='tight')
  846. plt.show()
  847. plt.close(fig)
  849. # %% [markdown]
  850. # ### Individual filters
  852. # %%
  853. stationary_stim_psd_avg.shape
  855. # %%
  856. # normalize the spectrum such that the are under the curve is 1
  857. norm_factor = np.sum(stationary_stim_psd_avg[:, :, :, :, :]*freq_spacing, axis=-1, keepdims=True)
  858. stationary_stim_psd_avg_norm = stationary_stim_psd_avg[:, :, :, :, :] / norm_factor
  859. moving_stim_psd_avg_norm = moving_stim_psd_avg[:, :, :, :, :] / norm_factor
  861. mat_shape = stationary_stim_psd_avg_norm.shape
  863. # repeat the process again for a range of theta values
  864. num_theta_values = 1000
  865. theta_values = np.linspace(0, 15, num_theta_values)
  866. moving_stim_energy = np.zeros((mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3], num_theta_values))
  867. stationary_stim_energy = np.zeros((mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3], num_theta_values))
  868. for i in tqdm(range(num_theta_values)):
  869. # compute the filter
  870. filter = temporal_filter(frequencies, theta_values[i])
  871. # compute the filtered spectra
  872. filtered_mean_stationary_psd = np.multiply(stationary_stim_psd_avg_norm, filter)
  873. filtered_mean_moving_psd = np.multiply(moving_stim_psd_avg_norm, filter)
  874. # compute the energy
  875. moving_stim_energy[:, :, :, :, i] = np.sum(filtered_mean_moving_psd*freq_spacing, axis=-1)
  876. stationary_stim_energy[:, :, :, :, i] = np.sum(filtered_mean_stationary_psd*freq_spacing, axis=-1)
  878. # %%
  879. # for given energy capacities in the stationary case (between 0 and 1)
  880. relative = True
  882. num_energy_levels = 100
  883. energy_capacity_stat = np.linspace(0, 1, num_energy_levels+1)
  884. delta_energy_capacity = np.linspace(0.01, 10, 10*num_energy_levels)
  885. if relative:
  886. energy_capacity_mov = np.outer(delta_energy_capacity, energy_capacity_stat) + energy_capacity_stat[np.newaxis, :] # relative
  887. delta_energy_capacity *= 100 # convert to percent
  888. else:
  889. energy_capacity_mov = delta_energy_capacity[:, np.newaxis] + energy_capacity_stat[np.newaxis, :]
  890. # energy_capacity_stat = np.logspace(-1, 0, 100)
  891. # delta_energy_capacity = np.logspace(-1, 0, 100)
  893. stationary_theta_interp = np.zeros((num_energy_levels+1, mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3]))
  894. moving_theta_interp = np.zeros((10*num_energy_levels, num_energy_levels+1, mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3]))
  896. # interpolate the energy values to get a smooth curve
  897. for i in range(mat_shape[0]):
  898. for j in range(mat_shape[1]):
  899. for k in range(mat_shape[2]):
  900. for l in range(mat_shape[3]):
  901. # interpolate the energy values to get a smooth curve
  902. stationary_theta_interp[:, i, j, k, l] = np.interp(energy_capacity_stat, stationary_stim_energy[i, j, k, l], theta_values)
  903. for m in range(num_energy_levels):
  904. # interpolate the energy values to get a smooth curve
  905. moving_theta_interp[:, m, i, j, k, l] = np.interp(energy_capacity_mov[:, m], moving_stim_energy[i, j, k, l], theta_values)
  907. # %%
  908. # for each filter, find the cutoff value of the autocorrelation function for different points in the energy landscape
  909. # create a matrix to store the cutoff values for each filter
  910. cutoff_frequency_mov = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3]))
  911. cutoff_frequency_stat = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3]))
  912. autocorr_func_mov = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3], 2*(len(frequencies)-1)))
  913. autocorr_func_stat = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3], 2*(len(frequencies)-1)))
  914. cutoff = 0.5
  915. for m in range(len(example_points)):
  916. for i in range(mat_shape[0]):
  917. for j in range(mat_shape[1]):
  918. for k in range(mat_shape[2]):
  919. for l in range(mat_shape[3]):
  920. # compute the filter
  921. example_theta_mov = moving_theta_interp[example_points[m][1], example_points[m][0], i, j, k, l]
  922. example_theta_stat = stationary_theta_interp[example_points[m][0], i, j, k, l]
  923. # compute the filter
  924. example_filter_mov = temporal_filter(frequencies, example_theta_mov)
  925. example_filter_stat = temporal_filter(frequencies, example_theta_stat)
  926. # compute the filtered spectra
  927. example_filtered_mean_stationary_psd = example_filter_stat # np.multiply(stationary_stim_psd_avg_norm[i, j, k, l], example_filter_stat)
  928. example_filtered_mean_moving_psd = example_filter_mov # np.multiply(moving_stim_psd_avg_norm[i, j, k, l], example_filter_mov)
  929. # compute the inverse fourier transform
  930. example_moving_autocorr = irfft(example_filtered_mean_moving_psd)
  931. example_stationary_autocorr = irfft(example_filtered_mean_stationary_psd)
  932. # normalize the autocorrelation function
  933. example_moving_autocorr /= np.max(example_moving_autocorr)
  934. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  935. # store the autocorrelation function
  936. autocorr_func_mov[m, i, j, k, l, :] = example_moving_autocorr
  937. autocorr_func_stat[m, i, j, k, l, :] = example_stationary_autocorr
  938. # find the cutoff frequency
  939. try:
  940. cutoff_frequency_mov[m, i, j, k, l] = np.where(example_moving_autocorr < cutoff)[0][0]
  941. except IndexError:
  942. cutoff_frequency_mov[m, i, j, k, l] = len(example_moving_autocorr)
  943. try:
  944. cutoff_frequency_stat[m, i, j, k, l] = np.where(example_stationary_autocorr < cutoff)[0][0]
  945. except IndexError:
  946. cutoff_frequency_stat[m, i, j, k, l] = len(example_stationary_autocorr)
  949. # %%
  950. # plot and save all autocorrelation functions
  951. # make a directory called example_autocorrelation_sim in ../manuscript_figures/
  952. import os
  953. if not os.path.exists('../manuscript_figures/example_autocorrelation_sim_-1_1_white_noise_band_pass'):
  954. os.makedirs('../manuscript_figures/example_autocorrelation_sim_-1_1_white_noise_band_pass')
  957. point_idx = 0
  958. position_idx = 4
  959. time = np.arange(0, len(example_moving_autocorr)) / 30*1000
  960. for ori_idx in range(mat_shape[0]):
  961. for phase_idx in range(mat_shape[1]):
  962. for wavelength_idx in range(mat_shape[2]):
  963. fig, ax = plt.subplots(1, 1, figsize=(4, 3))
  964. ax.plot(time, autocorr_func_mov[point_idx, ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:orange')
  965. ax.plot(time, autocorr_func_stat[point_idx, ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Stationary', color='tab:gray')
  966. ax.set_xlim(0, 500)
  967. ax.set_xlabel('Time (ms)')
  968. ax.set_ylabel('Autocorrelation')
  969. ax.set_title(f'Point {point_idx}, Ori {ori_idx}, Phase {phase_idx}, Wavelength {wavelength_idx}, Position {position_idx}')
  970. ax.legend()
  971. ax.set_ylim(-1, 1)
  972. plt.savefig(f'../manuscript_figures/example_autocorrelation_sim_-1_1_white_noise_band_pass/fig3_autocorr_point_{point_idx}_ori_{ori_idx}_phase_{phase_idx}_wavelength_{wavelength_idx}_position_{position_idx}.pdf', format='pdf', bbox_inches='tight')
  973. # close the figure
  974. plt.close(fig)
  976. # %%
  977. # plot the average autocorrelation function for each point
  978. avg_autocorr_func_mov = np.mean(autocorr_func_mov, axis=(1, 2, 3, 4))
  979. avg_autocorr_func_stat = np.mean(autocorr_func_stat, axis=(1, 2, 3, 4))
  980. std_autocorr_func_mov = np.std(autocorr_func_mov, axis=(1, 2, 3, 4))
  981. std_autocorr_func_stat = np.std(autocorr_func_stat, axis=(1, 2, 3, 4))
  982. time = np.arange(0, len(example_moving_autocorr)) / 30*1000
  983. # plot the average autocorrelation function for each point
  984. for point_idx in range(len(example_points)):
  985. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  986. ax.plot(time, avg_autocorr_func_mov[point_idx, :], label='Moving', color='tab:orange')
  987. ax.plot(time, avg_autocorr_func_stat[point_idx, :], label='Stationary', color='tab:gray')
  988. ax.fill_between(time, avg_autocorr_func_mov[point_idx, :] - std_autocorr_func_mov[point_idx, :], avg_autocorr_func_mov[point_idx, :] + std_autocorr_func_mov[point_idx, :], color='tab:orange', alpha=0.2)
  989. ax.fill_between(time, avg_autocorr_func_stat[point_idx, :] - std_autocorr_func_stat[point_idx, :], avg_autocorr_func_stat[point_idx, :] + std_autocorr_func_stat[point_idx, :], color='tab:gray', alpha=0.2)
  990. ax.set_xlim(0, 500)
  991. ax.set_ylim(-1, 1)
  992. plt.savefig(f'../manuscript_figures/fig3_avg_autocorr_point_{point_idx}.pdf', format='pdf', bbox_inches='tight')
  993. plt.show()
  994. plt.close(fig)
  996. # %%
  997. cutoff_frequency_stat /= 30
  998. cutoff_frequency_mov /= 30
  1000. # %%
  1001. cutoff_frequency_stat *= 1000
  1002. cutoff_frequency_mov *= 1000
  1004. # %%
  1005. np.unique(diff_cutoff)
  1007. # %%
  1008. np.linspace(-300, 300, 20)
  1010. # %%
  1011. # take the difference in cutoffs and plot the histogram
  1012. diff_cutoff = cutoff_frequency_mov - cutoff_frequency_stat
  1014. bins = np.linspace(-300, 300, 20)
  1016. idx = 3
  1018. plt.hist(diff_cutoff[idx].flatten(), bins=bins, density=True)
  1019. plt.xlabel('Difference of cutoff (moving - stationary) in ms')
  1020. plt.axvline(np.mean(diff_cutoff[idx]), color='k', linestyle='--')
  1021. plt.savefig('../manuscript_figures/fig3_hist_cutoff_diff_theory.pdf', format='pdf', bbox_inches='tight')
  1023. # %%
  1024. # plot the cutoff frequencies as a density scatter plot with jitter
  1025. point = 3
  1026. jitter_strength = 10
  1028. # Ensure the x and y axes span the same range
  1029. min_val = 0
  1030. max_val = max(np.max(cutoff_frequency_stat[point]), np.max(cutoff_frequency_mov[point]))
  1032. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1033. ax.set_title('Cutoff frequency')
  1034. ax.set_xlabel('Stationary cutoff frequency')
  1035. ax.set_ylabel('Moving cutoff frequency')
  1036. # add jitter to the points
  1037. cutoff_stat_jitter = cutoff_frequency_stat[point, :, :, :, :].flatten() + np.random.normal(0, jitter_strength, cutoff_frequency_stat[point, :, :, :, :].flatten().shape)
  1038. cutoff_mov_jitter = cutoff_frequency_mov[point, :, :, :, :].flatten() + np.random.normal(0, jitter_strength, cutoff_frequency_mov[point, :, :, :, :].flatten().shape)
  1040. # use gaussian kde to plot the density of points
  1041. from scipy.stats import gaussian_kde
  1043. kde = gaussian_kde([cutoff_stat_jitter, cutoff_mov_jitter])
  1045. # plot scatter plot colored by density
  1046. density = kde([cutoff_stat_jitter, cutoff_mov_jitter])
  1047. ax.scatter(cutoff_stat_jitter, cutoff_mov_jitter, c=density, cmap='viridis', s=1)
  1048. # add a colorbar
  1049. cbar = fig.colorbar(ax.collections[0], ax=ax, orientation='vertical')
  1050. cbar.set_label('Density')
  1051. # add a line with slope 1 which goes from 0 to the max value axes
  1052. ax.plot([min_val, max_val], [min_val, max_val], color='black', linestyle='--')
  1053. ax.set_xlim(0, 500)
  1054. ax.set_ylim(0, 500)
  1055. plt.savefig(f'../manuscript_figures/fig3_cutoff_frequency_point_{point}.pdf', format='pdf', bbox_inches='tight')
  1058. # %%
  1059. # plot the cutoff frequencies as a density scatter plot with jitter
  1060. point = 3
  1062. # Ensure the x and y axes span the same range
  1063. min_val = 0
  1064. max_val = max(np.max(cutoff_frequency_stat[point]), np.max(cutoff_frequency_mov[point]))
  1065. frame_length = 1/30*1000
  1066. bins = np.arange(min_val, max_val+frame_length, frame_length)
  1068. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1069. ax.set_title('Cutoff frequency')
  1070. ax.set_xlabel('Stationary cutoff lag (ms)')
  1071. ax.set_ylabel('Moving cutoff lag (ms)')
  1073. ax.hist2d(cutoff_frequency_stat[point, :, :, :, :].flatten(), cutoff_frequency_mov[point, :, :, :, :].flatten(), bins=[bins, bins], cmap='Blues')
  1074. # add a colorbar
  1075. cbar = fig.colorbar(ax.collections[0], ax=ax, orientation='vertical')
  1076. cbar.set_label('Density')
  1077. # add a line with slope 1 which goes from 0 to the max value axes
  1078. ax.plot([min_val, max_val], [min_val, max_val], color='black', linestyle='--')
  1080. plt.savefig(f'../manuscript_figures/fig3_cutoff_frequency_density_hist_point_{point}.pdf', format='pdf', bbox_inches='tight')
  1083. # %%
  1084. # print the fraction above the diagonal
  1085. above_diagonal = np.sum(cutoff_frequency_mov[point] > cutoff_frequency_stat[point]) / (cutoff_frequency_mov[point].flatten().shape[0])
  1086. print(f'Fraction of points above the diagonal: {above_diagonal:.2f}')
  1088. # %%
  1089. import matplotlib.colors as mcolors
  1091. norm = mcolors.LogNorm(vmin=np.min(theta_values), vmax=np.max(theta_values))
  1092. two_color_norm = mcolors.TwoSlopeNorm(vmin=-5, vcenter=0, vmax=5)
  1094. random_draws = 10
  1096. for i in range(random_draws):
  1097. ori_idx = np.random.randint(0, mat_shape[0])
  1098. phase_idx = np.random.randint(0, mat_shape[1])
  1099. wavelength_idx = np.random.randint(0, mat_shape[2])
  1100. position_idx = np.random.randint(0, mat_shape[3])
  1102. stationary_theta_test = stationary_theta_interp[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  1103. moving_theta_test = moving_theta_interp[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]
  1105. fig, ax = plt.subplots(2, 3, figsize=(16, 10))
  1106. im0 = ax[0, 0].imshow(np.repeat(stationary_theta_test[np.newaxis, :], 100, axis=0), origin='lower', norm=norm, cmap='viridis', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  1107. ax[0, 0].set_title('Stationary theta')
  1108. im1 = ax[0, 1].imshow(moving_theta_test, origin='lower', norm=norm, cmap='viridis', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  1109. ax[0, 1].set_title('Moving theta')
  1110. im2 = ax[0, 2].imshow(stationary_theta_test[np.newaxis, :]-moving_theta_test, origin='lower', norm=two_color_norm, cmap='bwr', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  1111. ax[0, 2].set_title('Difference in theta [stationary - moving]')
  1112. # add axis labels
  1113. ax[0, 0].set_ylabel('Percentage increase in energy capacity')
  1114. ax[0, 1].set_xlabel('Energy capacity')
  1115. fig.suptitle(f'Orientation: {orientation_arr[ori_idx]}, Phase: {phase_arr[phase_idx]}, Frequency: {freq_arr[low_spatial_freq_idx][wavelength_idx]}, Position: {position_arr[position_idx]}')
  1117. # add a colorbar without disrupting the layout of the other subplots
  1118. cbar = fig.colorbar(im2, ax=ax[0, 2], orientation='vertical', fraction=0.02, pad=0.04)
  1119. cbar = fig.colorbar(im1, ax=ax[0, 1], orientation='vertical', fraction=0.02, pad=0.04)
  1120. cbar = fig.colorbar(im0, ax=ax[0, 0], orientation='vertical', fraction=0.02, pad=0.04)
  1122. # plot the spectra and the energy vs. theta
  1123. ax[1, 0].plot(frequencies, stationary_stim_psd_avg[ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Stationary', color='tab:red')
  1124. ax[1, 0].plot(frequencies, moving_stim_psd_avg[ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:blue')
  1125. ax[1, 0].set_xlabel('Frequency (Hz)')
  1126. ax[1, 0].set_ylabel('Power')
  1127. ax[1, 0].set_title('Power spectrum')
  1128. ax[1, 0].set_yscale('log')
  1129. ax[1, 1].plot(theta_values, stationary_stim_energy[ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Stationary', color='tab:red')
  1130. ax[1, 1].plot(theta_values, moving_stim_energy[ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:blue')
  1131. ax[1, 1].set_xlabel('Theta')
  1132. ax[1, 1].set_ylabel('Energy')
  1133. ax[1, 1].set_title('Energy vs. theta')
  1134. # make last plot empty
  1135. ax[1, 2].axis('off')
  1137. plt.show()
  1139. # %% [markdown]
  1140. # ### Block trial stimuli
  1142. # %%
  1143. import numpy as np
  1144. # create white noise block structure, have 1 second of white noise with low variance and then 1 second of white noise with high variance
  1145. # the sampling rate is 10 ms.
  1146. num_samples = 100 # number of samples in the white noise block
  1147. low_noise = 0.1
  1148. high_noise = 1
  1149. sampling_rate = 100
  1151. # repeat this 100 times with new sampling so that in the end I have an array of shape (100, 300)
  1152. num_blocks = 100
  1153. white_noise_blocks = np.array([np.concatenate([
  1154. np.random.normal(0, low_noise, 20), # low variance
  1155. np.random.normal(0, high_noise, num_samples), # high variance
  1156. np.random.normal(0, low_noise, 80), # low variance
  1157. ]) for _ in range(num_blocks)])
  1159. # take the fft of the white noise blocks
  1160. white_noise_blocks_fft = np.fft.rfft(white_noise_blocks, axis=-1)
  1161. white_noise_frequencies = np.fft.rfftfreq(white_noise_blocks.shape[1], d=1/sampling_rate)
  1163. # %%
  1164. avg_white_noise_psd = np.mean(np.abs(white_noise_blocks_fft)**2, axis=0)[4:]
  1166. # %%
  1167. import matplotlib.pyplot as plt
  1168. plt.plot(white_noise_frequencies[4:], avg_white_noise_psd)
  1170. # %%
  1171. # for each filter, find the cutoff value of the autocorrelation function for different points in the energy landscape
  1172. # create a matrix to store the cutoff values for each filter
  1173. cutoff_frequency_mov = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3]))
  1174. cutoff_frequency_stat = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3]))
  1175. autocorr_func_mov = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3], len(white_noise_frequencies[4:])//2))
  1176. autocorr_func_stat = np.zeros((len(example_points), mat_shape[0], mat_shape[1], mat_shape[2], mat_shape[3], len(white_noise_frequencies[4:])//2))
  1177. cutoff = 0.5
  1178. for m in range(len(example_points)):
  1179. for i in range(mat_shape[0]):
  1180. for j in range(mat_shape[1]):
  1181. for k in range(mat_shape[2]):
  1182. for l in range(mat_shape[3]):
  1183. # compute the filter
  1184. example_theta_mov = moving_theta_interp[example_points[m][1], example_points[m][0], i, j, k, l]
  1185. example_theta_stat = stationary_theta_interp[example_points[m][0], i, j, k, l]
  1186. # compute the filter
  1187. example_filter_mov = temporal_filter(white_noise_frequencies[4:], example_theta_mov)
  1188. example_filter_stat = temporal_filter(white_noise_frequencies[4:], example_theta_stat)
  1189. # compute the filtered spectra
  1190. example_filtered_mean_stationary_psd = np.multiply(avg_white_noise_psd, example_filter_stat)
  1191. example_filtered_mean_moving_psd = np.multiply(avg_white_noise_psd, example_filter_mov)
  1192. # compute the inverse fourier transform
  1193. example_moving_autocorr = irfft(example_filtered_mean_moving_psd)[:len(example_filtered_mean_moving_psd)//2]
  1194. example_stationary_autocorr = irfft(example_filtered_mean_stationary_psd)[:len(example_filtered_mean_stationary_psd)//2]
  1195. # normalize the autocorrelation function
  1196. example_moving_autocorr /= np.max(example_moving_autocorr)
  1197. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  1198. # store the autocorrelation function
  1199. autocorr_func_mov[m, i, j, k, l, :] = example_moving_autocorr
  1200. autocorr_func_stat[m, i, j, k, l, :] = example_stationary_autocorr
  1201. # # find the cutoff frequency
  1202. try:
  1203. cutoff_frequency_mov[m, i, j, k, l] = np.where(example_moving_autocorr < cutoff)[0][0]
  1204. cutoff_frequency_stat[m, i, j, k, l] = np.where(example_stationary_autocorr < cutoff)[0][0]
  1205. except IndexError:
  1206. # if the cutoff is not found, set it to nan
  1207. cutoff_frequency_mov[m, i, j, k, l] = np.nan
  1208. cutoff_frequency_stat[m, i, j, k, l] = np.nan
  1211. # %%
  1212. # save the autocorrelation and cutoff frequency matrices
  1213. np.savez('../data/fig3_autocorr_cutoff_frequency_white_noise.npz',
  1214. autocorr_func_mov=autocorr_func_mov,
  1215. autocorr_func_stat=autocorr_func_stat,
  1216. cutoff_frequency_mov=cutoff_frequency_mov,
  1217. cutoff_frequency_stat=cutoff_frequency_stat,
  1218. white_noise_frequencies=white_noise_frequencies[4:],
  1219. example_points=example_points,
  1220. frequencies=frequencies[99:],
  1221. theta_values=theta_values,
  1222. stationary_theta_interp=stationary_theta_interp,
  1223. moving_theta_interp=moving_theta_interp)
  1225. # %%
  1226. time
  1228. # %%
  1229. autocorr_func_mov.shape
  1231. # %%
  1232. point_idx = 1
  1233. position_idx = 5
  1234. ori = 0
  1235. phase = 0
  1236. wavelength = 3
  1237. time = np.arange(0, len(white_noise_frequencies[4:])//2*10, 10)
  1238. plt.plot(time, autocorr_func_mov[point_idx, ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:orange')
  1239. plt.plot(time, autocorr_func_stat[point_idx, ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:gray')
  1240. plt.xlim(0, 500)
  1241. plt.axhline(0, color='black', linestyle='--')
  1243. # %%
  1244. # plot and save all autocorrelation functions
  1245. # make a directory called example_autocorrelation_sim in ../manuscript_figures/
  1246. import os
  1247. if not os.path.exists('../manuscript_figures/example_autocorrelation_sim_-1_1_block_stim'):
  1248. os.makedirs('../manuscript_figures/example_autocorrelation_sim_-1_1_block_stim')
  1251. point_idx = 1
  1252. position_idx = 4
  1253. for ori_idx in range(mat_shape[0]):
  1254. for phase_idx in range(mat_shape[1]):
  1255. for wavelength_idx in range(mat_shape[2]):
  1256. fig, ax = plt.subplots(1, 1, figsize=(4, 3))
  1257. ax.plot(time, autocorr_func_mov[point_idx, ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Moving', color='tab:orange')
  1258. ax.plot(time, autocorr_func_stat[point_idx, ori_idx, phase_idx, wavelength_idx, position_idx, :], label='Stationary', color='tab:gray')
  1259. ax.set_xlabel('Time (ms)')
  1260. ax.set_ylabel('Autocorrelation')
  1261. ax.set_title(f'Point {point_idx}, Ori {ori_idx}, Phase {phase_idx}, Wavelength {wavelength_idx}, Position {position_idx}')
  1262. ax.legend()
  1263. ax.set_ylim(-1, 1)
  1264. ax.set_xlim(0, 500)
  1265. plt.savefig(f'../manuscript_figures/example_autocorrelation_sim_-1_1_block_stim/fig3_autocorr_point_{point_idx}_ori_{ori_idx}_phase_{phase_idx}_wavelength_{wavelength_idx}_position_{position_idx}_block_stim.pdf', format='pdf', bbox_inches='tight')
  1266. # close the figure
  1267. plt.close(fig)
  1269. # %%
  1270. # plot the average autocorrelation function for each point
  1271. avg_autocorr_func_mov = np.mean(autocorr_func_mov, axis=(1, 2, 3, 4))
  1272. avg_autocorr_func_stat = np.mean(autocorr_func_stat, axis=(1, 2, 3, 4))
  1273. std_autocorr_func_mov = np.std(autocorr_func_mov, axis=(1, 2, 3, 4))
  1274. std_autocorr_func_stat = np.std(autocorr_func_stat, axis=(1, 2, 3, 4))
  1275. time = np.arange(0, len(avg_autocorr_func_mov[1])*10, 10)
  1276. # plot the average autocorrelation function for each point
  1277. for point_idx in range(len(example_points)):
  1278. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1279. ax.plot(time, avg_autocorr_func_mov[point_idx, :], label='Moving', color='tab:orange')
  1280. ax.plot(time, avg_autocorr_func_stat[point_idx, :], label='Stationary', color='tab:gray')
  1281. ax.fill_between(time, avg_autocorr_func_mov[point_idx, :] - std_autocorr_func_mov[point_idx, :], avg_autocorr_func_mov[point_idx, :] + std_autocorr_func_mov[point_idx, :], color='tab:orange', alpha=0.2)
  1282. ax.fill_between(time, avg_autocorr_func_stat[point_idx, :] - std_autocorr_func_stat[point_idx, :], avg_autocorr_func_stat[point_idx, :] + std_autocorr_func_stat[point_idx, :], color='tab:gray', alpha=0.2)
  1283. # ax.set_xlim(0, 1000)
  1284. ax.set_ylim(-1, 1)
  1285. plt.show()
  1286. plt.savefig(f'../manuscript_figures/fig3_avg_autocorr_point_{point_idx}_block_stimuli.pdf', format='pdf', bbox_inches='tight')
  1287. plt.close(fig)
  1289. # %%
  1290. # plot the cutoff frequencies as a density scatter plot with jitter
  1291. point = 2
  1293. jitter_strength = 10
  1295. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1296. ax.set_title('Cutoff frequency')
  1297. ax.set_xlabel('Stationary cutoff frequency')
  1298. ax.set_ylabel('Moving cutoff frequency')
  1299. # add jitter to the points and remove nan values
  1300. cutoff_frequency_stat_temp = cutoff_frequency_stat[point, :, :, :, :].flatten()
  1301. cutoff_frequency_mov_temp = cutoff_frequency_mov[point, :, :, :, :].flatten()
  1302. # check what percent of values in cutoff_frequency_mov are nan
  1303. print(f'Percent of NaN values in cutoffs: {np.sum(np.isnan(cutoff_frequency_mov_temp))/ cutoff_frequency_mov_temp.size * 100:.2f}%')
  1305. cutoff_frequency_stat_temp = cutoff_frequency_stat_temp[~np.isnan(cutoff_frequency_stat_temp)]*10
  1306. cutoff_frequency_mov_temp = cutoff_frequency_mov_temp[~np.isnan(cutoff_frequency_mov_temp)]*10
  1308. # Ensure the x and y axes span the same range
  1309. min_val = 0
  1310. max_val = max(np.max(cutoff_frequency_stat_temp), np.max(cutoff_frequency_mov_temp))
  1312. cutoff_stat_jitter = cutoff_frequency_stat_temp + np.random.normal(0, jitter_strength, cutoff_frequency_stat_temp.shape)
  1313. cutoff_mov_jitter = cutoff_frequency_mov_temp + np.random.normal(0, jitter_strength, cutoff_frequency_mov_temp.shape)
  1315. # use gaussian kde to plot the density of points
  1316. from scipy.stats import gaussian_kde
  1318. kde = gaussian_kde([cutoff_stat_jitter, cutoff_mov_jitter])
  1320. # plot scatter plot colored by density
  1321. density = kde([cutoff_stat_jitter, cutoff_mov_jitter])
  1322. ax.scatter(cutoff_stat_jitter, cutoff_mov_jitter, c=density, cmap='viridis', s=1)
  1323. # add a colorbar
  1324. cbar = fig.colorbar(ax.collections[0], ax=ax, orientation='vertical')
  1325. cbar.set_label('Density')
  1326. # add a line with slope 1 which goes from 0 to the max value axes
  1327. ax.plot([min_val, max_val], [min_val, max_val], color='black', linestyle='--')
  1328. ax.set_xlim(0, 500)
  1329. ax.set_ylim(0, 500)
  1331. plt.savefig(f'../manuscript_figures/fig3_cutoff_frequency_density_scatter_block_stim_point_{point}.pdf', format='pdf', bbox_inches='tight')
  1334. # %%
  1335. # plot the cutoff frequencies as a density scatter plot with jitter
  1336. point = 2
  1338. # Ensure the x and y axes span the same range
  1339. # add jitter to the points and remove nan values
  1340. cutoff_frequency_stat_temp = cutoff_frequency_stat[point, :, :, :, :].flatten()
  1341. cutoff_frequency_mov_temp = cutoff_frequency_mov[point, :, :, :, :].flatten()
  1342. # check what percent of values in cutoff_frequency_mov are nan
  1343. print(f'Percent of NaN values in cutoffs: {np.sum(np.isnan(cutoff_frequency_mov_temp))/ cutoff_frequency_mov_temp.size * 100:.2f}%')
  1345. cutoff_frequency_stat_temp = cutoff_frequency_stat_temp[~np.isnan(cutoff_frequency_stat_temp)]*10
  1346. cutoff_frequency_mov_temp = cutoff_frequency_mov_temp[~np.isnan(cutoff_frequency_mov_temp)]*10
  1348. # Ensure the x and y axes span the same range
  1349. min_val = 0
  1350. max_val = max(np.max(cutoff_frequency_stat_temp), np.max(cutoff_frequency_mov_temp))
  1351. frame_length = 10
  1352. bins = np.arange(min_val, 510, frame_length)
  1354. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1355. ax.set_title('Cutoff frequency')
  1356. ax.set_xlabel('Stationary cutoff lag (ms)')
  1357. ax.set_ylabel('Moving cutoff lag (ms)')
  1359. ax.hist2d(cutoff_frequency_stat_temp, cutoff_frequency_mov_temp, bins=[bins, bins], cmap='Blues')
  1360. # add a colorbar
  1361. cbar = fig.colorbar(ax.collections[0], ax=ax, orientation='vertical')
  1362. cbar.set_label('Density')
  1363. # add a line with slope 1 which goes from 0 to the max value axes
  1364. ax.plot([min_val, max_val], [min_val, max_val], color='black', linestyle='--')
  1366. plt.savefig(f'../manuscript_figures/fig3_cutoff_frequency_density_scatter_block_stim_point_{point}.pdf', format='pdf', bbox_inches='tight')
  1369. # %% [markdown]
  1370. # ## Normalized Utility
  1372. # %%
  1373. # compute the max info and energy for each filter
  1374. moving_max_info = np.max(moving_stim_information, axis=0)
  1375. stationary_max_info = np.max(stationary_stim_information, axis=0)
  1376. moving_max_energy = np.max(moving_stim_energy, axis=0)
  1377. stationary_max_energy = np.max(stationary_stim_energy, axis=0)
  1378. max_info = np.maximum(moving_max_info, stationary_max_info)
  1379. max_energy = np.maximum(moving_max_energy, stationary_max_energy)
  1380. # normalize
  1381. moving_stim_information_norm = moving_stim_information / max_info
  1382. stationary_stim_information_norm = stationary_stim_information / max_info
  1383. moving_stim_energy_norm = moving_stim_energy / max_energy
  1384. stationary_stim_energy_norm = stationary_stim_energy / max_energy
  1386. # %%
  1387. ori_idx = 0
  1388. phase_idx = 0
  1389. wavelength_idx = 0
  1390. position_idx = 4
  1391. moving_test_informations = moving_stim_information_norm[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  1392. stationary_test_informations = stationary_stim_information_norm[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  1393. moving_test_energy = moving_stim_energy_norm[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  1394. stationary_test_energy = stationary_stim_energy_norm[:, ori_idx, phase_idx, wavelength_idx, position_idx]
  1396. # create a subplot with 1 row and 3 columns. First column plot the information vs. theta, second column plot the energy vs. theta, third column plot the information vs. energy
  1397. fig, axs = plt.subplots(1, 3, figsize=(15, 5))
  1398. axs[0].plot(theta_values, moving_test_informations, 'o-', color='tab:blue', label='Moving test')
  1399. axs[0].plot(theta_values, stationary_test_informations, 'o-', color='tab:red', label='Stationary test')
  1400. axs[0].set_xlabel('Theta')
  1401. axs[0].set_ylabel('Information')
  1402. axs[0].set_title('Information vs. theta')
  1403. axs[0].legend()
  1404. axs[1].plot(theta_values, moving_test_energy, 'o-', color='tab:blue', label='Moving test')
  1405. axs[1].plot(theta_values, stationary_test_energy, 'o-', color='tab:red', label='Stationary test')
  1406. axs[1].set_xlabel('Theta')
  1407. axs[1].set_ylabel('Energy')
  1408. axs[1].set_title('Energy vs. theta')
  1409. axs[1].legend()
  1410. axs[2].plot(moving_test_energy, moving_test_informations, 'o-', color='tab:blue', label='Moving test')
  1411. axs[2].plot(stationary_test_energy, stationary_test_informations, 'o-', color='tab:red', label='Stationary test')
  1412. axs[2].set_xlabel('Energy')
  1413. axs[2].set_ylabel('Information')
  1414. axs[2].set_title('Information vs. energy')
  1415. axs[2].legend()
  1416. plt.tight_layout()
  1417. plt.show()
  1419. # %%
  1420. # compute the utility for a variety of lambdas
  1421. lmbd_arr = np.linspace(0, 100, 100)
  1423. moving_utility = np.zeros((lmbd_arr.shape[0], num_theta_values, len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  1424. stationary_utility = np.zeros((lmbd_arr.shape[0], num_theta_values, len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  1426. for i, lmbd in enumerate(lmbd_arr):
  1427. moving_utility[i] = moving_stim_information_norm - lmbd*moving_stim_energy_norm
  1428. stationary_utility[i] = stationary_stim_information_norm - lmbd*stationary_stim_energy_norm
  1430. # %%
  1431. # compute the optimal value of theta for each filter
  1432. optimal_theta_idx_moving = np.argmax(moving_utility, axis=1)
  1433. optimal_theta_idx_stationary = np.argmax(stationary_utility, axis=1)
  1435. optimal_theta_moving = theta_values[optimal_theta_idx_moving]
  1436. optimal_theta_stationary = theta_values[optimal_theta_idx_stationary]
  1438. # %%
  1439. optimal_theta_moving.shape
  1441. # %%
  1442. # take the difference and then the average
  1443. optimal_theta_diff = optimal_theta_moving - optimal_theta_stationary
  1444. optimal_theta_diff_avg = np.mean(optimal_theta_diff, axis=(1, 2, 4))
  1446. # %%
  1447. optimal_theta_diff.shape
  1449. # %%
  1450. plt.hist(np.min(optimal_theta_diff, axis=0).flatten(), bins=20)
  1452. # %%
  1453. # plot the average difference in optimal theta, color them with viridis
  1455. from matplotlib import cm
  1456. from matplotlib.colors import Normalize
  1458. colors_map = cm.get_cmap('viridis', len(wavelength_arr[low_spatial_freq_idx]))
  1460. fig, ax = plt.subplots(1, 1, figsize=(8, 6))
  1461. for i in range(wavelength_arr[low_spatial_freq_idx].shape[0]):
  1462. ax.plot(lmbd_arr, optimal_theta_diff_avg[:, i], color=colors_map(i))
  1466. # %%
  1467. ori_idx = 0
  1468. phase_idx = 0
  1469. wavelength_idx = 0
  1470. position_idx = 2
  1472. # plot the delta theta moving - stationary for a specific filter
  1473. fig, ax = plt.subplots(1, 1, figsize=(8, 6))
  1475. ax.plot(lmbd_arr, optimal_theta_moving[:, ori_idx, phase_idx, wavelength_idx, position_idx] - optimal_theta_stationary[:, ori_idx, phase_idx, wavelength_idx, position_idx])
  1477. # %%
  1478. import matplotlib.colors as mcolors
  1480. norm = mcolors.TwoSlopeNorm(vmin=np.min([moving_utility, stationary_utility]), vcenter=0, vmax=np.max([moving_utility, stationary_utility]))
  1482. # create a heatmap of the utility for theta vs. lambda for stationary and moving
  1483. fig, ax = plt.subplots(1, 2, figsize=(16, 6), sharey=True)
  1484. ax[1].imshow(moving_utility[:, :, ori_idx, phase_idx, wavelength_idx, position_idx], aspect='auto', cmap='bwr', origin='lower')
  1485. ax[1].set_title('Moving utility')
  1486. ax[0].imshow(stationary_utility[:, :, ori_idx, phase_idx, wavelength_idx, position_idx], aspect='auto', cmap='bwr', origin='lower')
  1487. ax[0].set_title('Stationary utility')
  1488. ax[0].set_xlabel('Theta')
  1489. ax[1].set_xlabel('Theta')
  1490. ax[0].set_ylabel('Lambda')
  1491. # set the tick labels
  1492. ax[0].set_xticks(np.arange(0, num_theta_values, 10))
  1493. ax[0].set_xticklabels(np.round(theta_values[::10], 2))
  1494. ax[1].set_xticks(np.arange(0, num_theta_values, 10))
  1495. ax[1].set_xticklabels(np.round(theta_values[::10], 2))
  1496. ax[0].set_yticks(np.arange(0, len(lmbd_arr), 10))
  1497. ax[0].set_yticklabels(np.round(lmbd_arr[::10], 2))
  1499. # %% [markdown]
  1500. # ## Gaussian filter
  1502. # %%
  1503. # repeat the process again for a range of theta values
  1504. num_theta_values = 50
  1505. array_shape = moving_stim_psd_avg.shape
  1506. theta_values_u = [0] # np.linspace(0, np.max(frequencies), num_theta_values)
  1507. theta_values_s = np.logspace(-1, 1, num_theta_values)
  1508. moving_stim_information = np.zeros((len(theta_values_u), len(theta_values_s), array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  1509. stationary_stim_information = np.zeros((len(theta_values_u), len(theta_values_s), array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  1510. moving_stim_energy = np.zeros((len(theta_values_u), len(theta_values_s), array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  1511. stationary_stim_energy = np.zeros((len(theta_values_u), len(theta_values_s), array_shape[0], array_shape[1], array_shape[2], array_shape[3]))
  1512. for i in tqdm(range(len(theta_values_u))):
  1513. for j in range(len(theta_values_s)):
  1514. # compute the filter
  1515. filter = temporal_filter_gauss(frequencies, theta_values_u[i], theta_values_s[j])
  1516. # compute the filtered spectra
  1517. filtered_stationary_stim_psd_avg = np.multiply(stationary_stim_psd_avg, filter)
  1518. filtered_moving_stim_psd_avg = np.multiply(moving_stim_psd_avg, filter)
  1519. # compute the information
  1520. moving_stim_information[i, j] = np.sum(np.log2(1 + filtered_moving_stim_psd_avg / noise), axis=-1)
  1521. stationary_stim_information[i, j] = np.sum(np.log2(1 + filtered_stationary_stim_psd_avg / noise), axis=-1)
  1522. # compute the energy
  1523. moving_stim_energy[i, j] = np.sum(filtered_moving_stim_psd_avg, axis=-1)
  1524. stationary_stim_energy[i, j] = np.sum(filtered_stationary_stim_psd_avg, axis=-1)
  1527. # %%
  1528. import matplotlib.colors as mcolors
  1530. norm = mcolors.CenteredNorm()
  1532. # plot the information as a heatmap of theta_u and theta_s
  1533. fig, ax = plt.subplots(1, 3, figsize=(24, 6))
  1534. im = ax[0].imshow(moving_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1535. ax[0].set_xlabel('SD')
  1536. ax[0].set_ylabel('Mean')
  1537. ax[0].set_title('Info')
  1539. ax[1].imshow(moving_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1540. ax[1].set_xlabel('SD')
  1541. ax[1].set_ylabel('Mean')
  1542. ax[1].set_title('Energy')
  1544. ax[2].imshow(moving_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]) - 0.1*moving_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1545. ax[2].set_xlabel('SD')
  1546. ax[2].set_ylabel('Mean')
  1549. # %%
  1550. import matplotlib.colors as mcolors
  1552. norm = mcolors.CenteredNorm()
  1554. # plot the information as a heatmap of theta_u and theta_s
  1555. fig, ax = plt.subplots(1, 3, figsize=(24, 6))
  1556. im = ax[0].imshow(stationary_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1557. ax[0].set_xlabel('SD')
  1558. ax[0].set_ylabel('Mean')
  1559. ax[0].set_title('Info')
  1561. ax[1].imshow(stationary_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1562. ax[1].set_xlabel('SD')
  1563. ax[1].set_ylabel('Mean')
  1564. ax[1].set_title('Energy')
  1566. ax[2].imshow(stationary_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_information[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]) - 0.1*stationary_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_energy[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1567. ax[2].set_xlabel('SD')
  1568. ax[2].set_ylabel('Mean')
  1571. # %%
  1572. plt.plot(moving_stim_information[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_information[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]) - 2 * moving_stim_energy[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_energy[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]), 'o-', color='tab:blue', label='Moving test')
  1573. plt.plot(stationary_stim_information[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_information[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]) - 2 * stationary_stim_energy[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]/np.max(moving_stim_energy[0, :, ori_idx, phase_idx, wavelength_idx, position_idx]), 'o-', color='tab:red', label='Stationary test')
  1575. # %%
  1576. # compute the max info and energy for each filter
  1577. moving_max_info = np.max(moving_stim_information[0], axis=0)
  1578. stationary_max_info = np.max(stationary_stim_information[0], axis=0)
  1579. moving_max_energy = np.max(moving_stim_energy[0], axis=0)
  1580. stationary_max_energy = np.max(stationary_stim_energy[0], axis=0)
  1581. max_info = np.maximum(moving_max_info, stationary_max_info)
  1582. max_energy = np.maximum(moving_max_energy, stationary_max_energy)
  1583. # normalize
  1584. moving_stim_information_norm = moving_stim_information[0] / max_info
  1585. stationary_stim_information_norm = stationary_stim_information[0] / max_info
  1586. moving_stim_energy_norm = moving_stim_energy[0] / max_energy
  1587. stationary_stim_energy_norm = stationary_stim_energy[0] / max_energy
  1589. # %%
  1590. moving_stim_information_norm.shape
  1592. # %%
  1593. # compute the utility for a variety of lambdas
  1594. lmbd_arr = np.linspace(0, 2, 50)
  1596. moving_utility = np.zeros((lmbd_arr.shape[0], num_theta_values, len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  1597. stationary_utility = np.zeros((lmbd_arr.shape[0], num_theta_values, len(orientation_arr), len(phase_arr), len(wavelength_arr[low_spatial_freq_idx]), len(position_arr)))
  1599. for i, lmbd in enumerate(lmbd_arr):
  1600. moving_utility[i] = moving_stim_information_norm - lmbd*moving_stim_energy_norm
  1601. stationary_utility[i] = stationary_stim_information_norm - lmbd*stationary_stim_energy_norm
  1603. # %%
  1604. moving_utility.shape
  1606. # %%
  1609. # %%
  1610. import matplotlib.colors as mcolors
  1612. norm = mcolors.TwoSlopeNorm(vmin=np.min([moving_utility, stationary_utility]), vcenter=0, vmax=np.max([moving_utility, stationary_utility]))
  1614. ori_idx = 0
  1615. phase_idx = 0
  1616. wavelength_idx = 10
  1617. position_idx = 4
  1619. # utility as as a heatmaps of lambda and theta
  1621. fig, ax = plt.subplots(1, 2, figsize=(16, 7))
  1622. ax[0].imshow(stationary_utility[:, :, ori_idx, phase_idx, wavelength_idx, position_idx].T, aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1623. ax[0].set_xlabel('Lambda')
  1624. ax[0].set_ylabel('SD')
  1625. im = ax[1].imshow(moving_utility[:, :, ori_idx, phase_idx, wavelength_idx, position_idx].T, aspect='auto', cmap='coolwarm', origin='lower', norm=norm)
  1626. ax[1].set_xlabel('Lambda')
  1627. ax[1].set_ylabel('SD')
  1628. ax[0].set_title('Stationary Utility')
  1629. ax[1].set_title('Moving Utility')
  1630. # for each lambda find the optimal theta and then plot a line overlaid on the heatmap connecting these values
  1631. optimal_theta_stationary_idx = np.argmin(np.abs(stationary_utility[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), axis=1)
  1632. optimal_theta_moving_idx = np.argmin(np.abs(moving_utility[:, :, ori_idx, phase_idx, wavelength_idx, position_idx]), axis=1)
  1633. optimal_theta_stationary = theta_values_s[optimal_theta_stationary_idx]
  1634. optimal_theta_moving = theta_values_s[optimal_theta_moving_idx]
  1635. ax[0].plot(np.arange(len(lmbd_arr)), optimal_theta_stationary_idx, color='k')
  1636. ax[1].plot(np.arange(len(lmbd_arr)), optimal_theta_moving_idx, color='k')
  1638. # %%
  1639. optimal_theta_stationary_idx.shape
  1641. # %%
  1642. plt.plot(optimal_theta_moving - optimal_theta_stationary)
  1644. # %% [markdown]
  1645. # ## Streamlined windowed and demeaned
  1647. # %%
  1648. # parse the data and compute fourier transform for each filter and video
  1649. environments = ['field', 'forest', 'orchard', 'tall_grass', 'pond']
  1651. num_videos = 10
  1652. vid_length = 50*30
  1653. low_spatial_freq_idx = np.arange(0, 31)
  1654. # low_spatial_freq_idx = np.arange(0, 19)
  1655. # high_spatial_freq_idx = np.arange(35, 70)
  1656. # low_spatial_freq_idx = np.arange(35, 70)
  1658. stationary_stim = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), vid_length))
  1659. moving_stim = np.zeros((len(environments), num_videos, len(orientation_arr), len(phase_arr), len(low_spatial_freq_idx), len(position_arr), vid_length))
  1661. # get the responses for each environment
  1662. for i, env_key in enumerate(environments):
  1663. stationary_count = 0
  1664. moving_count = 0
  1665. all_gabor_responses_env = all_gabor_responses[env_key]
  1666. print(all_gabor_responses_env.keys())
  1667. for vid_key in all_gabor_responses_env.keys():
  1668. if 'stationary' in vid_key:
  1669. stationary_stim[i, stationary_count, :, :, :, :, :] = all_gabor_responses_env[vid_key][()][:, :, low_spatial_freq_idx, :, :vid_length]
  1670. stationary_count += 1
  1671. if 'moving' in vid_key and 'free_moving' not in vid_key:
  1672. moving_stim[i, moving_count, :, :, :, :, :] = all_gabor_responses_env[vid_key][()][:, :, low_spatial_freq_idx, :, :vid_length]
  1673. moving_count += 1
  1675. # %%
  1676. from scipy.signal import welch
  1678. def compute_psd_welch(x, fs, window_size, overlap):
  1679. """
  1680. Compute one-sided PSD using Welch's method for a real-valued signal with variable window size and full overlap.
  1682. Parameters:
  1683. - x: 1D real-valued input signal
  1684. - fs: sampling frequency (Hz)
  1685. - window_size: length of each segment (samples)
  1687. Returns:
  1688. - freqs: frequency bins (Hz)
  1689. - psd: power spectral density (power/Hz)
  1690. """
  1691. freqs, psd = welch(
  1692. x,
  1693. fs=fs,
  1694. window='boxcar',
  1695. nperseg=window_size,
  1696. noverlap=overlap, # full overlap
  1697. return_onesided=True,
  1698. scaling='density',
  1699. detrend='constant'
  1700. )
  1701. return freqs, psd
  1703. # %%
  1704. # Compute the spectra for multiple window sizes: 10, 5, 2, and 1 seconds
  1705. sampling_rate = 30
  1706. window_lengths_sec = [5] # [10, 5, 2]
  1707. window_sizes = [int(w * sampling_rate) for w in window_lengths_sec]
  1709. # Prepare arrays to store PSDs for each window size
  1710. stationary_stim_psd_all = []
  1711. moving_stim_psd_all = []
  1712. frequencies_all = []
  1714. for window_size in window_sizes:
  1715. # Pre-allocate arrays for this window size
  1716. psd_shape = (
  1717. len(environments), num_videos, len(orientation_arr), len(phase_arr),
  1718. len(low_spatial_freq_idx), len(position_arr), window_size // 2 + 1
  1719. )
  1720. stationary_stim_psd = np.zeros(psd_shape)
  1721. moving_stim_psd = np.zeros(psd_shape)
  1722. frequencies = None
  1724. for i, env_key in enumerate(environments):
  1725. for j in range(num_videos):
  1726. _, stationary_stim_psd[i, j] = compute_psd_welch(
  1727. stationary_stim[i, j], sampling_rate, window_size, window_size // 2
  1728. )
  1729. frequencies, moving_stim_psd[i, j] = compute_psd_welch(
  1730. moving_stim[i, j], sampling_rate, window_size, window_size // 2
  1731. )
  1732. stationary_stim_psd_all.append(stationary_stim_psd)
  1733. moving_stim_psd_all.append(moving_stim_psd)
  1734. frequencies_all.append(frequencies)
  1736. # %%
  1737. stationary_stim_psd_avg_all = []
  1738. moving_stim_psd_avg_all = []
  1740. for i in range(len(window_sizes)):
  1741. stationary_stim_psd_avg_all.append(np.mean(stationary_stim_psd_all[i], axis=(0, 1)))
  1742. moving_stim_psd_avg_all.append(np.mean(moving_stim_psd_all[i], axis=(0, 1)))
  1744. # %%
  1745. def temporal_filter(temporal_freqs, theta, low_freq_cutoff=2):
  1746. # create the filter
  1747. filter = np.zeros(temporal_freqs.shape)
  1748. # set the filter to 1 for frequencies below theta
  1749. filter[temporal_freqs <= theta] = 1
  1750. # filter[temporal_freqs < low_freq_cutoff] = 1 - (low_freq_cutoff - temporal_freqs[temporal_freqs < low_freq_cutoff])
  1751. # create a linear dropoff for frequencies above theta with a fixed slope
  1752. filter[temporal_freqs > theta] = 1 - (temporal_freqs[temporal_freqs > theta] - theta)
  1753. # ensure that the filter is non-negative
  1754. filter[filter < 0] = 0
  1755. return np.square(filter)
  1757. # %%
  1758. mean_moving_psd_all = []
  1759. mean_stationary_psd_all = []
  1760. sd_moving_psd_all = []
  1761. sd_stationary_psd_all = []
  1763. for win_idx, (moving_stim_psd, stationary_stim_psd, frequencies) in enumerate(zip(moving_stim_psd_all, stationary_stim_psd_all, frequencies_all)):
  1764. mean_moving_psd = np.mean(moving_stim_psd, axis=(0, 1, 2, 3, 4, 5))
  1765. mean_stationary_psd = np.mean(stationary_stim_psd, axis=(0, 1, 2, 3, 4, 5))
  1766. sd_moving_psd = np.std(moving_stim_psd, axis=(0, 1, 2, 3, 4, 5))
  1767. sd_stationary_psd = np.std(stationary_stim_psd, axis=(0, 1, 2, 3, 4, 5))
  1769. mean_moving_psd_all.append(mean_moving_psd)
  1770. mean_stationary_psd_all.append(mean_stationary_psd)
  1771. sd_moving_psd_all.append(sd_moving_psd)
  1772. sd_stationary_psd_all.append(sd_stationary_psd)
  1774. # Plot and save
  1775. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1776. ax.plot(frequencies[1:-1], mean_stationary_psd[1:-1], label='Stationary', color='tab:gray')
  1777. ax.plot(frequencies[1:-1], mean_moving_psd[1:-1], label='Moving', color='tab:orange')
  1778. # Optionally add SD shading:
  1779. # ax.fill_between(frequencies[1:-1], mean_stationary_psd[1:-1] - sd_stationary_psd[1:-1], mean_stationary_psd[1:-1] + sd_stationary_psd[1:-1], color='tab:gray', alpha=0.2)
  1780. # ax.fill_between(frequencies[1:-1], mean_moving_psd[1:-1] - sd_moving_psd[1:-1], mean_moving_psd[1:-1] + sd_moving_psd[1:-1], color='tab:orange', alpha=0.2)
  1781. ax.set_xlabel('Frequency (Hz)')
  1782. ax.set_ylabel('Power')
  1783. ax.set_title(f'Mean power spectrum ({window_lengths_sec[win_idx]}s window)')
  1784. ax.set_yscale('log')
  1785. ax.legend()
  1786. plt.savefig(f'../manuscript_figures/fig3_psd_mean_full_window_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  1787. plt.show()
  1788. plt.close(fig)
  1790. # %%
  1791. # compute the average spectrum over all filters
  1792. # plot the mean spectrum
  1793. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  1794. ax.plot(frequencies[1:-1], mean_moving_psd[1:-1], label='Moving', color='tab:orange')
  1795. ax.plot(frequencies[1:-1], mean_stationary_psd[1:-1], label='Stationary', color='tab:gray')
  1796. ax.set_xlabel('Frequency (Hz)')
  1797. ax.set_ylabel('Power')
  1798. ax.set_yscale('log')
  1799. # ax.set_xscale('log')
  1800. ax.legend()
  1802. # find where the area under the stationary psd is a certain value
  1803. # auc_val = 1
  1804. # aoc_idx_stat = np.where(np.cumsum(mean_stationary_psd*freq_spacing) <= auc_val)[0][-1]
  1805. # # find the frequency at that index
  1806. # freq_val_stat = frequencies[aoc_idx_stat]
  1807. # auc_idx_mov = np.where(np.cumsum(mean_moving_psd*freq_spacing) <= auc_val)[0][-1]
  1808. # # find the frequency at that index
  1809. # freq_val_mov = frequencies[auc_idx_mov]
  1811. auc_val_stat = 1
  1812. auc_val_mov = 2
  1813. cum_stationary = np.cumsum(mean_stationary_psd * freq_spacing)
  1814. cum_moving = np.cumsum(mean_moving_psd * freq_spacing)
  1816. # Interpolate to find the frequency where the cumulative sum reaches auc_val
  1817. freq_val_stat = np.interp(auc_val_stat, cum_stationary, frequencies)
  1818. freq_val_mov = np.interp(auc_val_mov, cum_moving, frequencies)
  1820. # # plot filled area under the curve
  1821. # ax.fill_between(frequencies[1:-1], mean_stationary_psd[1:-1], where=(frequencies[1:-1] <= freq_val_stat), color='tab:gray', alpha=0.5)
  1822. # ax.fill_between(frequencies[1:-1], mean_moving_psd[1:-1], where=(frequencies[1:-1] <= freq_val_mov), color='tab:orange', alpha=0.5)
  1824. # Interpolate the PSD at the cutoff frequencies
  1825. psd_val_stat = np.interp(freq_val_stat, frequencies, mean_stationary_psd)
  1826. psd_val_mov = np.interp(freq_val_mov, frequencies, mean_moving_psd)
  1828. # For fill_between, extend the frequency and PSD arrays to include the interpolated cutoff point
  1829. def extend_for_fill(frequencies, psd, cutoff_freq, cutoff_psd):
  1830. mask = frequencies <= cutoff_freq
  1831. # Find the last index before the cutoff
  1832. last_idx = np.where(mask)[0][-1]
  1833. # Insert the cutoff point after last_idx
  1834. new_freqs = np.insert(frequencies[mask], last_idx + 1, cutoff_freq)
  1835. new_psd = np.insert(psd[mask], last_idx + 1, cutoff_psd)
  1836. return new_freqs, new_psd
  1838. freqs_stat_fill, psd_stat_fill = extend_for_fill(frequencies[1:-1], mean_stationary_psd[1:-1], freq_val_stat, psd_val_stat)
  1839. freqs_mov_fill, psd_mov_fill = extend_for_fill(frequencies[1:-1], mean_moving_psd[1:-1], freq_val_mov, psd_val_mov)
  1841. # plot filled area under the curve using the interpolated cutoff
  1842. ax.fill_between(freqs_stat_fill[:-1], psd_stat_fill[:-1], color='tab:gray', alpha=0.5)
  1843. ax.fill_between(freqs_mov_fill, psd_mov_fill, color='tab:orange', alpha=0.5)
  1844. ax.set_xlim(0, 15)
  1845. ax.set_ylim(1e-3, 4)
  1846. # remove top and right spines
  1847. ax.spines['top'].set_visible(False)
  1848. ax.spines['right'].set_visible(False)
  1849. ax.set_title(f'Mean power spectrum, moving cutoff: {freq_val_mov:.2f} Hz, stationary cutoff: {freq_val_stat:.2f} Hz')
  1852. plt.savefig(f'../manuscript_figures/fig3_psd_mean_normalize_2Hz_crop_auc_stat_{auc_val_stat}_auc_mov_{auc_val_mov}.pdf', format='pdf', bbox_inches='tight')
  1854. # %%
  1855. # Normalize the spectrum for each window so that the area under the curve is 1
  1856. normalized_mean_stationary_psd_all = []
  1857. normalized_mean_moving_psd_all = []
  1859. for win_idx in range(len(window_sizes)):
  1860. frequencies = frequencies_all[win_idx]
  1861. mean_stationary_psd = mean_stationary_psd_all[win_idx]
  1862. mean_moving_psd = mean_moving_psd_all[win_idx]
  1863. freq_spacing = np.diff(frequencies)[0]
  1864. norm_factor = np.sum(mean_stationary_psd * freq_spacing)
  1865. mean_stationary_psd_norm = mean_stationary_psd / norm_factor
  1866. mean_moving_psd_norm = mean_moving_psd / norm_factor
  1867. normalized_mean_stationary_psd_all.append(mean_stationary_psd_norm)
  1868. normalized_mean_moving_psd_all.append(mean_moving_psd_norm)
  1870. # %%
  1871. # Repeat the process for a range of theta values for each window size
  1872. num_theta_values = 1000
  1873. theta_values = np.linspace(0, 15, num_theta_values)
  1875. all_moving_stim_energy = []
  1876. all_stationary_stim_energy = []
  1878. for win_idx in range(len(window_sizes)):
  1879. frequencies = frequencies_all[win_idx]
  1880. mean_stationary_psd = normalized_mean_stationary_psd_all[win_idx]
  1881. mean_moving_psd = normalized_mean_moving_psd_all[win_idx]
  1882. freq_spacing = np.diff(frequencies)[0]
  1883. moving_stim_energy = np.zeros(num_theta_values)
  1884. stationary_stim_energy = np.zeros(num_theta_values)
  1885. for i in tqdm(range(num_theta_values), desc=f"Window {window_lengths_sec[win_idx]}s"):
  1886. # compute the filter
  1887. filter = temporal_filter(frequencies, theta_values[i])
  1888. # compute the filtered spectra
  1889. filtered_mean_stationary_psd = mean_stationary_psd * filter
  1890. filtered_mean_moving_psd = mean_moving_psd * filter
  1891. # compute the energy
  1892. moving_stim_energy[i] = np.sum(filtered_mean_moving_psd * freq_spacing)
  1893. stationary_stim_energy[i] = np.sum(filtered_mean_stationary_psd * freq_spacing)
  1894. all_moving_stim_energy.append(moving_stim_energy)
  1895. all_stationary_stim_energy.append(stationary_stim_energy)
  1897. # %%
  1898. # Plot the energy vs. theta for each window size
  1899. fig, ax = plt.subplots(1, len(window_sizes), figsize=(3*1.5*len(window_sizes), 3*1.5))
  1900. if len(window_sizes) == 1:
  1901. ax = [ax]
  1902. for win_idx in range(len(window_sizes)):
  1903. ax[win_idx].axhline(1, color='k', linestyle='--')
  1904. ax[win_idx].axhline(2, color='tab:orange', linestyle='--')
  1905. ax[win_idx].plot(theta_values, all_moving_stim_energy[win_idx], label='Moving', color='tab:orange')
  1906. ax[win_idx].plot(theta_values, all_stationary_stim_energy[win_idx], label='Stationary', color='tab:gray')
  1907. ax[win_idx].set_xlabel('Theta')
  1908. ax[win_idx].set_title(f'Window: {window_lengths_sec[win_idx]}s')
  1909. ax[win_idx].set_yscale('log')
  1910. ax[win_idx].set_xlim(0, 15)
  1911. ax[win_idx].set_ylim(.3, 4)
  1912. # remove the right and top spines
  1913. ax[win_idx].spines['right'].set_visible(False)
  1914. ax[win_idx].spines['top'].set_visible(False)
  1915. if win_idx == 0:
  1916. ax[win_idx].set_ylabel('Energy')
  1917. ax[win_idx].legend()
  1918. plt.tight_layout()
  1919. plt.savefig(f'../manuscript_figures/fig3_energy_vs_theta_filtered.pdf', format='pdf', bbox_inches='tight')
  1920. plt.show()
  1922. # %%
  1923. freq_val_stat - freq_val_mov
  1925. # %% [markdown]
  1926. # ### Full 2D optimization
  1928. # %%
  1929. relative = True
  1931. # for given energy capacities in the stationary case (between 0 and 1)
  1932. energy_capacity_stat = np.linspace(0, 1, 101)
  1933. # energy_capacity_mov = np.linspace(0, 2, 101)
  1934. delta_energy_capacity = np.linspace(0.01, 10, 1000)
  1935. # energy_capacity_stat = np.logspace(-1, 0, 100)
  1936. # delta_energy_capacity = np.logspace(-1, 1, 101)
  1938. # interpolate the energy values to get a smooth curve
  1939. stationary_theta_interp = np.interp(energy_capacity_stat, stationary_stim_energy, theta_values)
  1941. if relative:
  1942. energy_capacity_mov = np.outer(delta_energy_capacity, energy_capacity_stat) + energy_capacity_stat[np.newaxis, :] # relative
  1943. delta_energy_capacity *= 100 # convert to percent
  1944. else:
  1945. energy_capacity_mov = delta_energy_capacity[:, np.newaxis] + energy_capacity_stat[np.newaxis, :] # absolute
  1947. moving_theta_interp = np.interp(energy_capacity_mov, moving_stim_energy, theta_values)
  1950. # %%
  1951. import matplotlib.colors as mcolors
  1952. from matplotlib import cm
  1954. norm = mcolors.Normalize(vmin=np.min(theta_values), vmax=np.max(theta_values))
  1955. two_color_norm = mcolors.TwoSlopeNorm(vmin=-5, vcenter=0, vmax=5)
  1957. # example_points = [(50, 50), (75, 75), (95, 100)]
  1958. example_points = [(75, 9), (90, 19), (95, 29), (99, 49)]
  1960. fig, ax = plt.subplots(1, 3, figsize=(16, 5), width_ratios=[5, 5, 5], sharey=True)
  1961. im0 = ax[0].imshow(np.repeat(stationary_theta_interp[np.newaxis, :], len(delta_energy_capacity), axis=0), origin='lower', norm=norm, cmap='gray_r', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  1962. ax[0].set_title('Stationary theta')
  1963. im1 = ax[1].imshow(moving_theta_interp, origin='lower', norm=norm, cmap='gray_r', extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  1964. ax[1].set_title('Moving theta')
  1965. colors_map = mcolors.LinearSegmentedColormap.from_list('custom_cmap', [(1, 127/255, 14/255), (1, 1, 1), (0.5, 0.5, 0.5)], N=256)
  1966. im2 = ax[2].imshow(stationary_theta_interp[np.newaxis, :]-moving_theta_interp, origin='lower', norm=two_color_norm, cmap=colors_map, extent=[energy_capacity_stat[0], energy_capacity_stat[-1], delta_energy_capacity[0], delta_energy_capacity[-1]], aspect='auto')
  1967. ax[2].set_title('Difference in theta [stationary - moving]')
  1968. # ax[2].set_xlim(0.5, 1)
  1969. # ax[2].set_ylim(0.5, 1)
  1970. # add axis labels
  1971. if relative:
  1972. ax[0].set_ylabel('Percentage increase in energy capacity')
  1973. else:
  1974. ax[0].set_ylabel('Increase in energy capacity')
  1976. fig.supxlabel('Energy capacity')
  1978. # plot the example points in three different colors which are distinct from red, white and blue
  1979. colors_map = cm.get_cmap('Blues', len(example_points)+2)
  1980. # plot the example points in the first subplot
  1981. for i, point in enumerate(example_points):
  1982. ax[0].scatter(energy_capacity_stat[point[0]], delta_energy_capacity[point[1]], color=colors_map(i+1), s=100, edgecolor='black', label='Example point ' + str(i+1))
  1983. ax[1].scatter(energy_capacity_stat[point[0]], delta_energy_capacity[point[1]], color=colors_map(i+1), s=100, edgecolor='black')
  1984. ax[2].scatter(energy_capacity_stat[point[0]], delta_energy_capacity[point[1]], color=colors_map(i+1), s=100, edgecolor='black')
  1986. ax[0].set_yscale('log')
  1988. # add a colorbar without disrupting the layout of the other subplots
  1989. cbar = fig.colorbar(im2, ax=ax[2], orientation='vertical', pad=0.04)
  1990. cbar = fig.colorbar(im1, ax=ax[1], orientation='vertical', pad=0.04)
  1991. cbar = fig.colorbar(im0, ax=ax[0], orientation='vertical', pad=0.04)
  1992. plt.savefig(f'../manuscript_figures/fig3_diff_theta_stat_energy_vs_delta_energy_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  1994. # %%
  1995. from scipy.fft import irfft
  1997. # compute the autocorrelation function for the example points by extracting the theta values in the stationary and moving conditions and then applying the filter to white noise and then computing the autocorrelation function as the inverse fourier transform of the PSD
  1998. alpha_arr = [0.25, 0.5, 0.75, 1]
  1999. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2000. for i in range(len(example_points)):
  2001. # compute the filter
  2002. example_theta_mov = moving_theta_interp[example_points[i][1], example_points[i][0]]
  2003. example_theta_stat = stationary_theta_interp[example_points[i][0]]
  2004. # compute the filter
  2005. example_filter_mov = temporal_filter(frequencies, example_theta_mov)
  2006. example_filter_stat = temporal_filter(frequencies, example_theta_stat)
  2007. # compute the inverse fourier transform
  2008. example_moving_autocorr = irfft(example_filter_mov/freq_spacing)
  2009. example_stationary_autocorr = irfft(example_filter_stat/freq_spacing)
  2010. # normalize the autocorrelation function
  2011. example_moving_autocorr /= np.max(example_moving_autocorr)
  2012. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  2013. # plot the autocorrelation function
  2014. time = np.arange(0, len(example_moving_autocorr)) / 30*1000
  2015. ax.plot(time, example_moving_autocorr, label='Moving', color='tab:orange', alpha=alpha_arr[i])
  2016. ax.plot(time, example_stationary_autocorr, label='Stationary', color='tab:gray', alpha=alpha_arr[i])
  2017. ax.set_xlabel('Time (ms)')
  2018. ax.set_ylabel('Autocorrelation')
  2019. # ax.legend()
  2020. ax.set_xlim(0, 500)
  2021. plt.savefig(f'../manuscript_figures/fig3_autocorr_example_{i}_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  2022. plt.show()
  2023. plt.close(fig)
  2025. # %% [markdown]
  2026. # ### 1D optimization
  2028. # %%
  2029. # Repeat for all windows, but keep energy_capacity_stat fixed at 1
  2031. relative = True
  2033. energy_capacity_stat = 1.0 # Always use 1 for stationary energy capacity
  2034. delta_energy_capacity = np.linspace(0.01, 3, 300)
  2036. all_stationary_theta_interp = []
  2037. all_moving_theta_interp = []
  2039. for win_idx in range(len(window_sizes)):
  2040. stationary_stim_energy = all_stationary_stim_energy[win_idx]
  2041. moving_stim_energy = all_moving_stim_energy[win_idx]
  2042. # theta_values should be defined as before
  2043. # If not, use: theta_values = np.linspace(0, 15, stationary_stim_energy.shape[0])
  2045. # Interpolate stationary theta for energy_capacity_stat = 1
  2046. stationary_theta_interp = np.interp(energy_capacity_stat, stationary_stim_energy, theta_values)
  2048. if relative:
  2049. energy_capacity_mov = delta_energy_capacity*energy_capacity_stat + energy_capacity_stat # relative to stat
  2050. delta_energy_capacity_scaled = delta_energy_capacity * 100 # percent
  2051. else:
  2052. energy_capacity_mov = delta_energy_capacity + energy_capacity_stat # absolute
  2054. # Interpolate moving theta for each delta
  2055. moving_theta_interp = np.interp(energy_capacity_mov, moving_stim_energy, theta_values)
  2057. all_stationary_theta_interp.append(stationary_theta_interp)
  2058. all_moving_theta_interp.append(moving_theta_interp)
  2060. # %%
  2061. import matplotlib.colors as mcolors
  2062. from matplotlib import cm
  2064. # example_points = [(50, 50), (75, 75), (95, 100)]
  2065. example_points = [49, 99, 149]
  2067. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2068. for win_idx in range(len(window_lengths_sec)):
  2069. ax.plot(delta_energy_capacity*100, all_stationary_theta_interp[win_idx]-all_moving_theta_interp[win_idx], label='Window: ' + str(window_lengths_sec[win_idx]) + 's', color=cm.Reds((win_idx+1)/len(window_lengths_sec)))
  2070. ax.scatter(delta_energy_capacity[example_points]*100, all_stationary_theta_interp[win_idx]-all_moving_theta_interp[win_idx][example_points], color=cm.viridis((np.arange(len(example_points))+1)/len(example_points)), s=50, edgecolor='black')
  2071. ax.set_xlabel('Increase in energy capacity (%)')
  2072. ax.set_ylabel('Theta difference (stationary - moving)')
  2073. ax.legend()
  2075. plt.savefig(f'../manuscript_figures/fig3_diff_theta_fixed_stat_energy_vs_delta_energy_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  2077. # %%
  2078. from scipy.fft import irfft
  2080. # Repeat for each window size
  2081. alpha_arr = [0.25, 0.5, 0.75, 1]
  2082. for win_idx in range(len(window_sizes)):
  2083. frequencies = frequencies_all[win_idx]
  2084. freq_spacing = np.diff(frequencies)[0]
  2085. stationary_theta_interp = all_stationary_theta_interp[win_idx]
  2086. moving_theta_interp = all_moving_theta_interp[win_idx]
  2087. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2088. for i in range(len(example_points)):
  2089. # compute the filter
  2090. example_theta_mov = moving_theta_interp[example_points[i]]
  2091. example_theta_stat = stationary_theta_interp
  2092. # compute the filter
  2093. example_filter_mov = temporal_filter(frequencies, example_theta_mov)
  2094. example_filter_stat = temporal_filter(frequencies, example_theta_stat)
  2095. # set the DC to zero
  2096. example_filter_mov[0] = 0
  2097. example_filter_stat[0] = 0
  2098. # compute the inverse fourier transform
  2099. example_moving_autocorr = irfft(example_filter_mov)
  2100. example_stationary_autocorr = irfft(example_filter_stat)
  2101. # normalize the autocorrelation function
  2102. example_moving_autocorr /= np.max(example_moving_autocorr)
  2103. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  2104. # plot the autocorrelation function
  2105. time = np.arange(0, len(example_moving_autocorr)) / 30*1000
  2106. ax.plot(time, example_moving_autocorr, color=cm.viridis((np.arange(len(example_points))+1)/len(example_points))[i], label=f'$\\Delta E$ {delta_energy_capacity[example_points[i]]*100}')
  2107. ax.set_xlabel('Time (ms)')
  2108. ax.set_ylabel('Autocorrelation')
  2109. ax.set_title(f'Autocorrelation for {window_lengths_sec[win_idx]}s window')
  2110. ax.set_xlim(0, 500)
  2111. ax.plot(time, example_stationary_autocorr, color='tab:gray', label='stationary')
  2112. ax.legend()
  2113. plt.savefig(f'../manuscript_figures/fig3_autocorr_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  2114. plt.show()
  2115. plt.close(fig)
  2117. # %% [markdown]
  2118. # ### Individual filters
  2120. # %%
  2121. # Repeat the normalization and energy computation for all window sizes
  2122. all_moving_stim_energy = []
  2123. all_stationary_stim_energy = []
  2125. num_theta_values = 1000
  2126. theta_values = np.linspace(0, 15, num_theta_values)
  2128. all_stationary_theta_interp = []
  2129. all_moving_theta_interp = []
  2131. num_energy_levels = 100
  2132. energy_capacity_stat = np.linspace(0, 1, num_energy_levels+1)
  2133. delta_energy_capacity = np.linspace(0.01, 10, 10*num_energy_levels)
  2134. relative = True
  2136. for win_idx in range(len(window_sizes)):
  2137. stationary_stim_psd_avg = stationary_stim_psd_avg_all[win_idx]
  2138. moving_stim_psd_avg = moving_stim_psd_avg_all[win_idx]
  2139. # Use the correct frequency array for this window size
  2140. frequencies = frequencies_all[win_idx]
  2141. freq_spacing = np.diff(frequencies)[0]
  2143. # Normalize the spectrum such that the area under the curve is 1
  2144. norm_factor = np.sum(stationary_stim_psd_avg * freq_spacing, axis=-1, keepdims=True)
  2145. stationary_stim_psd_avg_norm = stationary_stim_psd_avg / norm_factor
  2146. moving_stim_psd_avg_norm = moving_stim_psd_avg / norm_factor
  2148. mat_shape = stationary_stim_psd_avg_norm.shape
  2150. # Compute energy for a range of theta values
  2151. moving_stim_energy = np.zeros((*mat_shape[:-1], num_theta_values))
  2152. stationary_stim_energy = np.zeros((*mat_shape[:-1], num_theta_values))
  2153. for i in tqdm(range(num_theta_values), desc=f"Window {window_lengths_sec[win_idx]}s"):
  2154. filter = temporal_filter(frequencies, theta_values[i])
  2155. filtered_mean_stationary_psd = stationary_stim_psd_avg_norm * filter
  2156. filtered_mean_moving_psd = moving_stim_psd_avg_norm * filter
  2157. moving_stim_energy[:, :, :, :, i] = np.sum(filtered_mean_moving_psd * freq_spacing, axis=-1)
  2158. stationary_stim_energy[:, :, :, :, i] = np.sum(filtered_mean_stationary_psd * freq_spacing, axis=-1)
  2159. all_moving_stim_energy.append(moving_stim_energy)
  2160. all_stationary_stim_energy.append(stationary_stim_energy)
  2162. # --- Interpolation block (from your prompt) ---
  2163. if relative:
  2164. energy_capacity_mov = np.outer(delta_energy_capacity, energy_capacity_stat) + energy_capacity_stat[np.newaxis, :]
  2165. delta_energy_capacity_scaled = delta_energy_capacity * 100 # convert to percent
  2166. else:
  2167. energy_capacity_mov = delta_energy_capacity[:, np.newaxis] + energy_capacity_stat[np.newaxis, :]
  2168. stationary_theta_interp = np.zeros((num_energy_levels+1, *mat_shape[:-1]))
  2169. moving_theta_interp = np.zeros((10*num_energy_levels, num_energy_levels+1, *mat_shape[:-1]))
  2171. for i in range(mat_shape[0]):
  2172. for j in range(mat_shape[1]):
  2173. for k in range(mat_shape[2]):
  2174. for l in range(mat_shape[3]):
  2175. stationary_theta_interp[:, i, j, k, l] = np.interp(
  2176. energy_capacity_stat, stationary_stim_energy[i, j, k, l], theta_values
  2177. )
  2178. for m in range(num_energy_levels):
  2179. moving_theta_interp[:, m, i, j, k, l] = np.interp(
  2180. energy_capacity_mov[:, m], moving_stim_energy[i, j, k, l], theta_values
  2181. )
  2182. all_stationary_theta_interp.append(stationary_theta_interp)
  2183. all_moving_theta_interp.append(moving_theta_interp)
  2185. # %%
  2186. # For each window size, compute the cutoff value of the autocorrelation function for different points in the energy landscape
  2188. cutoff = 0.5
  2189. example_points = [(75, 9), (90, 19), (95, 29), (99, 49)] # or use your own
  2191. all_cutoff_frequency_mov = []
  2192. all_cutoff_frequency_stat = []
  2193. all_autocorr_func_mov = []
  2194. all_autocorr_func_stat = []
  2196. for win_idx in range(len(all_stationary_theta_interp)):
  2197. stationary_theta_interp = all_stationary_theta_interp[win_idx]
  2198. moving_theta_interp = all_moving_theta_interp[win_idx]
  2199. frequencies = frequencies_all[win_idx]
  2200. mat_shape = stationary_theta_interp.shape[1:] # (ori, phase, freq, pos)
  2202. cutoff_frequency_mov = np.zeros((len(example_points), *mat_shape))
  2203. cutoff_frequency_stat = np.zeros((len(example_points), *mat_shape))
  2204. autocorr_func_mov = np.zeros((len(example_points), *mat_shape, 2*(len(frequencies)-1)))
  2205. autocorr_func_stat = np.zeros((len(example_points), *mat_shape, 2*(len(frequencies)-1)))
  2207. for m, (stat_idx, mov_idx) in enumerate(example_points):
  2208. for i in range(mat_shape[0]):
  2209. for j in range(mat_shape[1]):
  2210. for k in range(mat_shape[2]):
  2211. for l in range(mat_shape[3]):
  2212. # Get theta values for this filter and energy point
  2213. example_theta_stat = stationary_theta_interp[stat_idx, i, j, k, l]
  2214. example_theta_mov = moving_theta_interp[mov_idx, stat_idx, i, j, k, l]
  2215. # Compute temporal filters
  2216. example_filter_stat = temporal_filter(frequencies, example_theta_stat)
  2217. example_filter_mov = temporal_filter(frequencies, example_theta_mov)
  2218. example_filter_stat[0] = 0
  2219. example_filter_mov[0] = 0
  2220. # Compute autocorrelation (inverse FFT of filter)
  2221. example_stationary_autocorr = irfft(example_filter_stat)
  2222. example_moving_autocorr = irfft(example_filter_mov)
  2223. # Normalize
  2224. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  2225. example_moving_autocorr /= np.max(example_moving_autocorr)
  2226. # Store
  2227. autocorr_func_stat[m, i, j, k, l, :] = example_stationary_autocorr
  2228. autocorr_func_mov[m, i, j, k, l, :] = example_moving_autocorr
  2229. # Find cutoff
  2230. try:
  2231. cutoff_frequency_stat[m, i, j, k, l] = np.where(example_stationary_autocorr < cutoff)[0][0]
  2232. except IndexError:
  2233. cutoff_frequency_stat[m, i, j, k, l] = len(example_stationary_autocorr)
  2234. try:
  2235. cutoff_frequency_mov[m, i, j, k, l] = np.where(example_moving_autocorr < cutoff)[0][0]
  2236. except IndexError:
  2237. cutoff_frequency_mov[m, i, j, k, l] = len(example_moving_autocorr)
  2238. all_cutoff_frequency_mov.append(cutoff_frequency_mov)
  2239. all_cutoff_frequency_stat.append(cutoff_frequency_stat)
  2240. all_autocorr_func_mov.append(autocorr_func_mov)
  2241. all_autocorr_func_stat.append(autocorr_func_stat)
  2243. # Now all_cutoff_frequency_mov, all_cutoff_frequency_stat, all_autocorr_func_mov, all_autocorr_func_stat
  2244. # are lists, one per window size, containing the results for each window.
  2246. # %%
  2247. # Plot the average autocorrelation function for each point and window size
  2248. for win_idx in range(len(all_autocorr_func_mov)):
  2249. avg_autocorr_func_mov = np.mean(all_autocorr_func_mov[win_idx], axis=(1, 2, 3, 4))
  2250. avg_autocorr_func_stat = np.mean(all_autocorr_func_stat[win_idx], axis=(1, 2, 3, 4))
  2251. std_autocorr_func_mov = np.std(all_autocorr_func_mov[win_idx], axis=(1, 2, 3, 4))
  2252. std_autocorr_func_stat = np.std(all_autocorr_func_stat[win_idx], axis=(1, 2, 3, 4))
  2253. # Use the correct time axis for this window size
  2254. time = np.arange(0, avg_autocorr_func_mov.shape[1]) / 30 * 1000
  2255. for point_idx in range(len(example_points)):
  2256. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2257. ax.plot(time, avg_autocorr_func_mov[point_idx, :], label='Moving', color='tab:orange')
  2258. ax.plot(time, avg_autocorr_func_stat[point_idx, :], label='Stationary', color='tab:gray')
  2259. ax.fill_between(
  2260. time,
  2261. avg_autocorr_func_mov[point_idx, :] - std_autocorr_func_mov[point_idx, :],
  2262. avg_autocorr_func_mov[point_idx, :] + std_autocorr_func_mov[point_idx, :],
  2263. color='tab:orange', alpha=0.2
  2264. )
  2265. ax.fill_between(
  2266. time,
  2267. avg_autocorr_func_stat[point_idx, :] - std_autocorr_func_stat[point_idx, :],
  2268. avg_autocorr_func_stat[point_idx, :] + std_autocorr_func_stat[point_idx, :],
  2269. color='tab:gray', alpha=0.2
  2270. )
  2271. ax.set_xlim(0, 500)
  2272. ax.set_ylim(-1, 1)
  2273. ax.set_title(f'Window size: {window_lengths_sec[win_idx]}s, Point {point_idx}')
  2274. plt.savefig(f'../manuscript_figures/fig3_avg_autocorr_point_{point_idx}_window_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  2275. plt.show()
  2276. plt.close(fig)
  2278. # ...existing code...
  2280. # %%
  2281. all_cutoff_frequency_mov = np.array(all_cutoff_frequency_mov)
  2282. all_cutoff_frequency_stat = np.array(all_cutoff_frequency_stat)
  2283. all_cutoff_frequency_stat = all_cutoff_frequency_stat / 30 * 1000
  2284. all_cutoff_frequency_mov = all_cutoff_frequency_mov / 30 * 1000
  2285. diff_cutoff_frequency = all_cutoff_frequency_mov - all_cutoff_frequency_stat
  2287. # %%
  2288. # take the difference in cutoffs and plot the histogram
  2289. bins = np.linspace(-300, 300, 20)
  2291. window_idx = 3
  2292. point_idx = 3
  2294. for window_idx in range(len(window_sizes)):
  2295. for point_idx in range(len(example_points)):
  2296. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2297. ax.hist(diff_cutoff_frequency[window_idx, point_idx].flatten(), bins=bins, density=True)
  2298. ax.set_xlabel('Difference of cutoff (moving - stationary) in ms')
  2299. ax.axvline(np.mean(diff_cutoff_frequency[window_idx, point_idx]), color='k', linestyle='--')
  2300. ax.set_title(f'Window size: {window_sizes[window_idx]/30}s, Point {point_idx}')
  2301. plt.savefig(f'../manuscript_figures/fig3_hist_cutoff_diff_theory_window_{window_sizes[window_idx]}_point_{point_idx}.pdf', format='pdf', bbox_inches='tight')
  2302. plt.show()
  2303. plt.close(fig)
  2305. # %%
  2306. diff_cutoff_frequency.shape
  2308. # %%
  2309. avg_diff_cutoff_frequency = np.mean(diff_cutoff_frequency, axis=(2, 3, 4, 5))
  2310. median_diff_cutoff_frequency = np.median(diff_cutoff_frequency, axis=(2, 3, 4, 5))
  2312. # %%
  2313. norm = mcolors.TwoSlopeNorm(vmin=-100, vcenter=0, vmax=100)
  2314. energy_capacity_stat = np.linspace(0, 1, num_energy_levels+1)
  2315. delta_energy_capacity = np.linspace(0.01, 10, 10*num_energy_levels)
  2317. im = plt.imshow(avg_diff_cutoff_frequency.T, norm=norm, cmap='bwr', aspect='auto')
  2318. plt.xticks(np.arange(len(window_sizes)), np.array(window_sizes)/30)
  2319. plt.xlabel('Window size (s)')
  2320. plt.ylabel('Constraint point')
  2321. plt.yticks(np.arange(len(example_points)), [f'$E_s$ {energy_capacity_stat[example_points[i][0]]:.2f}\n$\\Delta E$ {delta_energy_capacity[example_points[i][1]]:.2f}' for i in range(len(example_points))])
  2322. plt.colorbar(im, label='moving cutoff - stationary cutoff (ms)')
  2323. plt.suptitle('Average difference in cutoff frequency for different energy constraints')
  2324. plt.savefig(f'../manuscript_figures/fig3_avg_diff_cutoff_frequency.pdf', format='pdf', bbox_inches='tight')
  2327. # %%
  2328. norm = mcolors.TwoSlopeNorm(vmin=-100, vcenter=0, vmax=100)
  2329. energy_capacity_stat = np.linspace(0, 1, num_energy_levels+1)
  2330. delta_energy_capacity = np.linspace(0.01, 10, 10*num_energy_levels)
  2332. im = plt.imshow(median_diff_cutoff_frequency.T, norm=norm, cmap='bwr', aspect='auto')
  2333. plt.xticks(np.arange(len(window_sizes)), np.array(window_sizes)/30)
  2334. plt.xlabel('Window size (s)')
  2335. plt.ylabel('Constraint point')
  2336. plt.yticks(np.arange(len(example_points)), [f'$E_s$ {energy_capacity_stat[example_points[i][0]]:.2f}\n$\\Delta E$ {delta_energy_capacity[example_points[i][1]]:.2f}' for i in range(len(example_points))])
  2337. plt.colorbar(im, label='moving cutoff - stationary cutoff (ms)')
  2338. plt.suptitle('Median difference in cutoff frequency for different energy constraints')
  2339. plt.savefig(f'../manuscript_figures/fig3_median_diff_cutoff_frequency.pdf', format='pdf', bbox_inches='tight')
  2342. # %% [markdown]
  2343. # #### 1D optimization
  2345. # %%
  2346. # Repeat the normalization and energy computation for all window sizes
  2347. all_moving_stim_energy = []
  2348. all_stationary_stim_energy = []
  2350. num_theta_values = 1000
  2351. theta_values = np.linspace(0, 15, num_theta_values)
  2353. all_stationary_theta_interp = []
  2354. all_moving_theta_interp = []
  2356. num_energy_levels = 30
  2357. energy_capacity_stat = 1.0
  2358. delta_energy_capacity = np.linspace(0.01, 3, 10*num_energy_levels)
  2359. relative = True
  2361. for win_idx in range(len(window_sizes)):
  2362. stationary_stim_psd_avg = stationary_stim_psd_avg_all[win_idx]
  2363. moving_stim_psd_avg = moving_stim_psd_avg_all[win_idx]
  2364. # Use the correct frequency array for this window size
  2365. frequencies = frequencies_all[win_idx]
  2366. freq_spacing = np.diff(frequencies)[0]
  2368. # Normalize the spectrum such that the area under the curve is 1
  2369. norm_factor = np.sum(stationary_stim_psd_avg * freq_spacing, axis=-1, keepdims=True)
  2370. stationary_stim_psd_avg_norm = stationary_stim_psd_avg / norm_factor
  2371. moving_stim_psd_avg_norm = moving_stim_psd_avg / norm_factor
  2373. mat_shape = stationary_stim_psd_avg_norm.shape
  2375. # Compute energy for a range of theta values
  2376. moving_stim_energy = np.zeros((*mat_shape[:-1], num_theta_values))
  2377. stationary_stim_energy = np.zeros((*mat_shape[:-1], num_theta_values))
  2378. for i in tqdm(range(num_theta_values), desc=f"Window {window_lengths_sec[win_idx]}s"):
  2379. filter = temporal_filter(frequencies, theta_values[i])
  2380. filtered_mean_stationary_psd = stationary_stim_psd_avg_norm * filter
  2381. filtered_mean_moving_psd = moving_stim_psd_avg_norm * filter
  2382. moving_stim_energy[:, :, :, :, i] = np.sum(filtered_mean_moving_psd * freq_spacing, axis=-1)
  2383. stationary_stim_energy[:, :, :, :, i] = np.sum(filtered_mean_stationary_psd * freq_spacing, axis=-1)
  2384. all_moving_stim_energy.append(moving_stim_energy)
  2385. all_stationary_stim_energy.append(stationary_stim_energy)
  2387. # --- Interpolation block (from your prompt) ---
  2388. if relative:
  2389. energy_capacity_mov = delta_energy_capacity*energy_capacity_stat + energy_capacity_stat # relative to stat
  2390. delta_energy_capacity_scaled = delta_energy_capacity * 100 # percent
  2391. else:
  2392. energy_capacity_mov = delta_energy_capacity[:, np.newaxis] + energy_capacity_stat[np.newaxis, :]
  2393. stationary_theta_interp = np.zeros(mat_shape[:-1])
  2394. moving_theta_interp = np.zeros((10*num_energy_levels, *mat_shape[:-1]))
  2396. for i in range(mat_shape[0]):
  2397. for j in range(mat_shape[1]):
  2398. for k in range(mat_shape[2]):
  2399. for l in range(mat_shape[3]):
  2400. stationary_theta_interp[i, j, k, l] = np.interp(
  2401. energy_capacity_stat, stationary_stim_energy[i, j, k, l], theta_values
  2402. )
  2403. moving_theta_interp[:, i, j, k, l] = np.interp(
  2404. energy_capacity_mov, moving_stim_energy[i, j, k, l], theta_values
  2405. )
  2406. all_stationary_theta_interp.append(stationary_theta_interp)
  2407. all_moving_theta_interp.append(moving_theta_interp)
  2409. # %%
  2410. moving_theta_interp[49, 0, -1, :, 4]
  2412. # %%
  2413. from scipy.fft import irfft
  2415. cutoff = 0.5
  2416. example_points = [49, 99, 149] # or your own, for 1D fixed stationary energy
  2418. all_cutoff_frequency_mov = []
  2419. all_cutoff_frequency_stat = []
  2420. all_autocorr_func_mov = []
  2421. all_autocorr_func_stat = []
  2423. for win_idx in range(len(all_stationary_theta_interp)):
  2424. stationary_theta_interp = all_stationary_theta_interp[win_idx] # shape: (filters...)
  2425. moving_theta_interp = all_moving_theta_interp[win_idx] # shape: (delta_energy, filters...)
  2426. frequencies = frequencies_all[win_idx]
  2427. mat_shape = stationary_theta_interp.shape # (ori, phase, freq, pos)
  2429. cutoff_frequency_mov = np.zeros((len(example_points), *mat_shape))
  2430. cutoff_frequency_stat = np.zeros((len(example_points), *mat_shape))
  2431. autocorr_func_mov = np.zeros((len(example_points), *mat_shape, 2*(len(frequencies)-1)))
  2432. autocorr_func_stat = np.zeros((len(example_points), *mat_shape, 2*(len(frequencies)-1)))
  2434. for m, mov_idx in enumerate(example_points):
  2435. for i in range(mat_shape[0]):
  2436. for j in range(mat_shape[1]):
  2437. for k in range(mat_shape[2]):
  2438. for l in range(mat_shape[3]):
  2439. # Get theta values for this filter and energy point
  2440. example_theta_stat = stationary_theta_interp[i, j, k, l]
  2441. example_theta_mov = moving_theta_interp[mov_idx, i, j, k, l]
  2442. # Compute temporal filters
  2443. example_filter_stat = temporal_filter(frequencies, example_theta_stat)
  2444. example_filter_mov = temporal_filter(frequencies, example_theta_mov)
  2445. example_filter_stat[0] = 0
  2446. example_filter_mov[0] = 0
  2447. # Compute autocorrelation (inverse FFT of filter)
  2448. example_stationary_autocorr = irfft(example_filter_stat)
  2449. example_moving_autocorr = irfft(example_filter_mov)
  2450. # Normalize
  2451. example_stationary_autocorr /= np.max(example_stationary_autocorr)
  2452. example_moving_autocorr /= np.max(example_moving_autocorr)
  2453. # Store
  2454. autocorr_func_stat[m, i, j, k, l, :] = example_stationary_autocorr
  2455. autocorr_func_mov[m, i, j, k, l, :] = example_moving_autocorr
  2456. # Find cutoff
  2457. try:
  2458. cutoff_frequency_stat[m, i, j, k, l] = np.where(example_stationary_autocorr < cutoff)[0][0]
  2459. except IndexError:
  2460. cutoff_frequency_stat[m, i, j, k, l] = len(example_stationary_autocorr)
  2461. try:
  2462. cutoff_frequency_mov[m, i, j, k, l] = np.where(example_moving_autocorr < cutoff)[0][0]
  2463. except IndexError:
  2464. cutoff_frequency_mov[m, i, j, k, l] = len(example_moving_autocorr)
  2465. all_cutoff_frequency_mov.append(cutoff_frequency_mov)
  2466. all_cutoff_frequency_stat.append(cutoff_frequency_stat)
  2467. all_autocorr_func_mov.append(autocorr_func_mov)
  2468. all_autocorr_func_stat.append(autocorr_func_stat)
  2470. # %%
  2471. # do a paired t-test on the cutoff frequencies for each point and window size
  2472. from scipy.stats import ttest_rel
  2473. mov_cutoff = all_cutoff_frequency_mov[0, 0]
  2474. stat_cutoff = all_cutoff_frequency_stat[0, 0]
  2476. t_stat, p_value = ttest_rel(mov_cutoff.flatten(), stat_cutoff.flatten())
  2477. print(f"T-statistic: {t_stat}, P-value: {p_value}")
  2479. # %%
  2480. import sys
  2481. print(sys.float_info.min)
  2483. # %%
  2484. # Plot the average autocorrelation function for each point and window size
  2485. for win_idx in range(len(all_autocorr_func_mov)):
  2486. avg_autocorr_func_mov = np.mean(all_autocorr_func_mov[win_idx], axis=(1, 2, 3, 4))
  2487. avg_autocorr_func_stat = np.mean(all_autocorr_func_stat[win_idx], axis=(1, 2, 3, 4))
  2488. std_autocorr_func_mov = np.std(all_autocorr_func_mov[win_idx], axis=(1, 2, 3, 4))
  2489. std_autocorr_func_stat = np.std(all_autocorr_func_stat[win_idx], axis=(1, 2, 3, 4))
  2490. # Use the correct time axis for this window size
  2491. time = np.arange(0, avg_autocorr_func_mov.shape[1]) / 30 * 1000
  2492. for point_idx in range(len(example_points)):
  2493. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2494. ax.plot(time, avg_autocorr_func_mov[point_idx, :], label='Moving', color='tab:orange')
  2495. ax.plot(time, avg_autocorr_func_stat[point_idx, :], label='Stationary', color='tab:gray')
  2496. ax.fill_between(
  2497. time,
  2498. avg_autocorr_func_mov[point_idx, :] - std_autocorr_func_mov[point_idx, :],
  2499. avg_autocorr_func_mov[point_idx, :] + std_autocorr_func_mov[point_idx, :],
  2500. color='tab:orange', alpha=0.2
  2501. )
  2502. ax.fill_between(
  2503. time,
  2504. avg_autocorr_func_stat[point_idx, :] - std_autocorr_func_stat[point_idx, :],
  2505. avg_autocorr_func_stat[point_idx, :] + std_autocorr_func_stat[point_idx, :],
  2506. color='tab:gray', alpha=0.2
  2507. )
  2508. ax.set_xlim(0, 500)
  2509. ax.set_ylim(-1, 1)
  2510. ax.set_title(f'Window size: {window_lengths_sec[win_idx]}s, Point {point_idx}')
  2511. plt.savefig(f'../manuscript_figures/fig3_avg_autocorr_point_{point_idx}_window_{window_lengths_sec[win_idx]}s.pdf', format='pdf', bbox_inches='tight')
  2512. plt.show()
  2513. plt.close(fig)
  2515. # ...existing code...
  2517. # %%
  2518. all_autocorr_func_mov[1].shape
  2520. # %%
  2521. import os
  2523. # Directory to save figures
  2524. save_dir = "../manuscript_figures/autocorr_individual_corrected"
  2525. os.makedirs(save_dir, exist_ok=True)
  2527. # plot each individual autocorrelation function for a given windows size and position, for all example points, orientations, phases, and wavelengths
  2528. window_idx = 1
  2529. position_idx = 4 # choose a position index, e.g., 0 for the first position
  2530. example_points = [49] # or your own, for 1D fixed stationary energy
  2531. frequencies = frequencies_all[window_idx]
  2532. mat_shape = all_autocorr_func_mov[window_idx].shape[1:5] # (ori, phase, freq, pos)
  2534. i = 0
  2535. point_idx = 49
  2536. for ori_idx in range(mat_shape[0]):
  2537. for phase_idx in range(mat_shape[1]):
  2538. for freq_id in range(mat_shape[2]):
  2539. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2540. ax.plot(
  2541. np.arange(0, len(all_autocorr_func_mov[window_idx][i, ori_idx, phase_idx, freq_id, position_idx])) / 30 * 1000,
  2542. all_autocorr_func_mov[window_idx][i, ori_idx, phase_idx, freq_id, position_idx],
  2543. label='Moving', color='tab:orange'
  2544. )
  2545. ax.plot(
  2546. np.arange(0, len(all_autocorr_func_stat[window_idx][i, ori_idx, phase_idx, freq_id, position_idx])) / 30 * 1000,
  2547. all_autocorr_func_stat[window_idx][i, ori_idx, phase_idx, freq_id, position_idx],
  2548. label='Stationary', color='tab:gray'
  2549. )
  2550. ax.set_xlim(0, 500)
  2551. ax.set_ylim(-1, 1)
  2552. ax.set_title(f'Window size: {window_lengths_sec[window_idx]}s\n$\\Delta E$ {delta_energy_capacity[point_idx]*100:.2f}, Ori {orientation_arr[ori_idx]:.2f}, Phase {phase_arr[phase_idx]:.2f}, Freq {freq_arr[freq_id]:.2f} cpd')
  2553. ax.set_xlabel('Time (ms)')
  2554. ax.set_ylabel('Autocorrelation')
  2555. ax.legend()
  2556. # Save figure
  2557. fname = (
  2558. f"autocorr_win{window_lengths_sec[window_idx]}s_point{point_idx}_"
  2559. f"ori{ori_idx}_phase{phase_idx}_freq{freq_id}_pos{position_idx}.pdf"
  2560. )
  2561. plt.savefig(os.path.join(save_dir, fname), bbox_inches='tight', format='pdf')
  2562. plt.close(fig)
  2564. # %%
  2565. all_cutoff_frequency_mov = np.array(all_cutoff_frequency_mov)
  2566. all_cutoff_frequency_stat = np.array(all_cutoff_frequency_stat)
  2567. all_cutoff_frequency_stat = all_cutoff_frequency_stat / 30 * 1000
  2568. all_cutoff_frequency_mov = all_cutoff_frequency_mov / 30 * 1000
  2569. diff_cutoff_frequency = all_cutoff_frequency_mov - all_cutoff_frequency_stat
  2571. # %%
  2572. # take the difference in cutoffs and plot the histogram
  2573. bins = np.linspace(-300, 300, 20)
  2575. window_idx = 3
  2576. point_idx = 3
  2578. for window_idx in range(len(window_sizes)):
  2579. for point_idx in range(len(example_points)):
  2580. fig, ax = plt.subplots(1, 1, figsize=(4*1.5, 3*1.5))
  2581. ax.hist(diff_cutoff_frequency[window_idx, point_idx].flatten(), bins=bins, density=True)
  2582. ax.set_xlabel('Difference of cutoff (moving - stationary) in ms')
  2583. ax.axvline(np.mean(diff_cutoff_frequency[window_idx, point_idx]), color='k', linestyle='--')
  2584. ax.set_title(f'Window size: {window_sizes[window_idx]/30}s, Point {point_idx}')
  2585. plt.savefig(f'../manuscript_figures/fig3_hist_cutoff_diff_theory_window_{window_sizes[window_idx]}_point_{point_idx}.pdf', format='pdf', bbox_inches='tight')
  2586. plt.show()
  2587. plt.close(fig)
  2589. # %%

temporal_filtering_analysis.ipynb at commit 24c076f, no license · at the source

Overview

Authors: Jonathan M. Gant1,2, Wiktor F. Młynarski1,3
ORCID iDs: Jonathan M. Gant, Wiktor F. Młynarski
  1. Faculty of Biology, LMU, Munich, Germany
  2. Graduate School of Systemic Neurosciences, Munich, Germany
  3. Bernstein Center for Computational Neuroscience Munich, Munich, Germany
Institutions: Ludwig-Maximilians-Universität München (Germany); Bernstein Center for Computational Neuroscience Munich (Germany)
Journal: Science advances, volume 12, issue 35, article eaed4172
Dates: received 27 October 2025; accepted 20 July 2026; published online 28 August 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1126/sciadv.aed4172 · PMID 42664357 · PMCID PMC13524049 · OpenAlex W7154616600
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: cellular / molecular (subfield)
Methods: Smoothing, state filtering, decompositions, Statistics, Connectivity, Machine learning, Single-unit activity, calcium imaging
MeSH: Locomotion*, Models, Neurological*, Primates*, Rodentia*, Sensory Receptor Cells*, Animals (* major topic)
Journal subjects: Neuroscience
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Citations: not cited yet (Europe PMC); 64 references in the paper

Abstract

Behavior modulates the activity of sensory systems in multiple ways: from gain changes in individual neurons to changing interactions in neural populations. These effects are not universal; while movement has a strong influence on sensory coding in rodents, its impact on primates is less prominent. The diversity of effects that locomotion exerts on sensory neurons, as well as disparities between species, raises questions about the existence of universal principles that may underlie sensation during behavior. We propose that sensory systems are internally modulated to match systematic changes in stimulus statistics caused by locomotion, to facilitate an accurate and efficient sensory code. We find that model neurons, adapted to stimuli recorded during movement in natural environments, predict and reproduce a broad spectrum of experimental observations in rodents and primates. This simple principle of maintaining coding efficiency across behavioral states reconciles the diversity of ways in which locomotion modulates visual coding in different animal species.

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

Zenodo 20624372

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data, code, and materials availability:”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (31 files), h5py (21 files), Matplotlib (20 files), SciPy (19 files), scikit-learn (7 files), OpenCV (6 files), pandas (2 files), PyTorch (2 files), Numba (1 file), rpy2 (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
32 files
At the source: doi.org/10.5281/zenodo.20624372

mlynarski-group/locomotion-modulation

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 24c076f12556f265925d8408ef9f6bff8131385f, 10 June 2026
Languages: Jupyter (18), Python (13)
Size: 42 files, 31 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, environment (requirements.txt), 18 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (31 files), h5py (21 files), Matplotlib (20 files), SciPy (19 files), scikit-learn (7 files), OpenCV (6 files), pandas (2 files), PyTorch (2 files), Numba (1 file), rpy2 (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
32 files
At the source: github.com/Mlynarski-Group/locomotion-modulation/tree/v1.0.0

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

Datasets cited

Data, code, and materials availability

All data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. The data (natural videos) used in the analysis are available at https://doi.org/10.12751/g-node.14m4gq. The code is publicly available at https://doi.org/10.5281/zenodo.20624372. No new materials were generated in this study.

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, 6 MeSH terms, 59 references.

Cite

This paper

Gant, J. M., & Młynarski, W. F. (2026). Locomotion optimizes sensory representations through a computational principle shared by rodents and primates. Science advances, 12(35), eaed4172. https://doi.org/10.1126/sciadv.aed4172

BibTeX

@article{gant2026locomotion,
author = {Gant, Jonathan M. and Młynarski, Wiktor F.},
title = {{Locomotion optimizes sensory representations through a computational principle shared by rodents and primates}},
journal = {Science advances},
year = {2026},
month = aug,
volume = {12},
number = {35},
pages = {eaed4172},
publisher = {American Association for the Advancement of Science},
issn = {2375-2548},
doi = {10.1126/sciadv.aed4172},
url = {https://doi.org/10.1126/sciadv.aed4172},
pmid = {42664357},
pmcid = {PMC13524049}
}

RIS

TY - JOUR
AU - Gant, Jonathan M.
AU - Młynarski, Wiktor F.
TI - Locomotion optimizes sensory representations through a computational principle shared by rodents and primates
T2 - Science advances
J2 - Sci Adv
PY - 2026
DA - 2026/08/28
VL - 12
IS - 35
SP - eaed4172
SN - 2375-2548
PB - American Association for the Advancement of Science
DO - 10.1126/sciadv.aed4172
UR - https://doi.org/10.1126/sciadv.aed4172
LA - en
ER -

CSL-JSON

{
"id": "10.1126/sciadv.aed4172",
"type": "article-journal",
"title": "Locomotion optimizes sensory representations through a computational principle shared by rodents and primates",
"container-title": "Science advances",
"author": [
{
"family": "Gant",
"given": "Jonathan M."
},
{
"family": "Młynarski",
"given": "Wiktor F."
}
],
"container-title-short": "Sci Adv",
"volume": "12",
"issue": "35",
"page": "eaed4172",
"DOI": "10.1126/sciadv.aed4172",
"PMID": "42664357",
"PMCID": "PMC13524049",
"ISSN": "2375-2548",
"publisher": "American Association for the Advancement of Science",
"URL": "https://doi.org/10.1126/sciadv.aed4172",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
28
]
]
}
}

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.1371/journal.pcbi.1014123 [code]
Subunit-specific behavioral modulation of sensory tuning in the visual cortex.
Journal: PLoS computational biology
In common: pandas, Matplotlib, NumPy, 10 references, author Wiktor F. Młynarski
[2] doi:10.1038/s41467-026-75347-4 [code]
Sleep reveals dynamics integrating and segregating movement and stimulus representations in V1.
Journal: Nature communications
In common: Numba, PyTorch, scikit-learn, 4 other tools, 7 references
[3] doi:10.1038/s41467-026-71667-7
Behavioural states control binocular vision through input-specific mechanisms.
Journal: Nature communications
In common: 9 references
[4] doi:10.1016/j.isci.2026.116055 [code]
Mapping the transcriptional diversity of calcium signaling in the mouse and human brain.
Journal: iScience
In common: rpy2, Numba, h5py, 7 other tools
[5] doi:10.1038/s41586-026-10348-3 [code]
An enteric neuron ionotropic receptor regulates salt stress resistance.
Journal: Nature
In common: Numba, OpenCV, h5py, 7 other tools, cellular / molecular
[6] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: rpy2, Numba, OpenCV, 6 other tools, cellular / molecular
[7] doi:10.7554/elife.109717 [code]
Retrosplenial cortex enables context-dependent goal-directed sensorimotor transformation.
Journal: eLife
In common: Numba, OpenCV, h5py, 6 other tools, 1 reference
[8] doi:10.1038/s41593-026-02267-3 [code]
Spatial proteomic analysis in human Alzheimer's disease brains enables identification of microenvironment-dependent microglial cell states.
Journal: Nature neuroscience
In common: rpy2, Numba, h5py, 6 other tools, cellular / molecular
[9] doi:10.1016/j.isci.2026.116825 [code]
Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: Numba, OpenCV, h5py, 7 other tools
[10] doi:10.1038/s41467-026-72057-9 [code]
Sex-specific behavioral feedback modulates sensorimotor processing and drives flexible social behavior.
Journal: Nature communications
In common: Numba, OpenCV, h5py, 7 other tools

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.