Processing experimental data for flashed gratings

exp/Stimulus.py

@@ -0,0 +1,43 @@
+import numpy as np
+import pandas as pd
+
+
+
+data_folder = "../data/"
+stim_file = "Stiminfo_PVCre_2021_0012_s06_e14.csv"
+# Stimulus DataFrame
+stim = pd.read_csv(data_folder+stim_file)
+num_trial = len(stim)
+stim_val = {}
+trial_stim_id = {}
+stim_val['pair'], trial_stim_id['pair'] = np.unique(stim[['grat_orientation', 'grat_phase']], return_inverse=True, axis=0) # 50 trials per orientation-phase pair
+stim_val['orientation'], trial_stim_id['orientation'] = np.unique(stim['grat_orientation'], return_inverse=True) # 1000 trials per orientation
+stim_val['phase'], trial_stim_id['phase'] = np.unique(stim['grat_phase'], return_inverse=True) # 1000 trials per phase
+key_list = ['pair', 'orientation', 'phase']
+stim_id_trial = {}
+num_stim = {}
+for key in key_list:
+    stim_id_trial[key] = [np.where(trial_stim_id[key] == i)[0] for i in range(len(stim_val[key]))]
+    num_stim[key] = len(stim_val[key])
+
+# stim_id_trial[i] returns an array of trial indices where the stimulus has the value of the ith stimulus.
+def pair_id(i):
+    return (int(i//num_stim['orientation']), i%num_stim['orientation'])
+
+trial_pair_id = np.array([pair_id(i) for i in trial_stim_id['pair']], dtype=object)
+pair_val = stim_val['pair'].reshape(num_stim['orientation'], num_stim['phase'],2)
+
+# for each (orientation_id = i, phase_id = j) find the trial indices
+pair_trial_id = np.ndarray((num_stim['orientation'], num_stim['phase']), dtype=object)
+for i in range(num_stim['pair']):
+    pair_trial_id[pair_id(i)] = stim_id_trial['pair'][i]
+
+# s.close()
+
+stim_data = {'stim_val':stim_val, 'trial_stim_id':trial_stim_id, 'key_list':key_list, 'num_trial':num_trial,
+             'trial_pair_id':trial_pair_id, 'pair_val':pair_val, 'pair_trial_id':pair_trial_id, 'stim_id_trial':stim_id_trial, 'num_stim':num_stim}
+
+import pickle
+file = open('stim_data.pkl', 'wb')
+pickle.dump(stim_data, file)
+file.close()

exp/load.py

@@ -0,0 +1,44 @@
+import pandas as pd
+import numpy as np
+import pickle as pkl 
+
+data_folder = "../data/"
+
+stim_file = "Stiminfo_PVCre_2021_0012_s06_e14.csv"
+stim = pd.read_csv(data_folder+stim_file)
+
+spike_times_file = "Spiketimes_PVCre_2021_0012_s06_e14.npy"
+spike_times = np.load(data_folder+spike_times_file, allow_pickle=True)
+active = [len(spike_times[i]) > 0 for i in range(len(spike_times))]
+spike_times = spike_times[np.where(active)]
+
+num_unit = len(spike_times)
+num_trial = len(stim)
+
+# sort by firing rate
+
+num_spike = list(map(len, spike_times))
+# num_spike = np.array([len(spike_times[i]) for i in range(len(spike_times))])
+spike_times = spike_times[np.argsort(num_spike)[::-1]]
+
+files = ['stim_data.pkl', 'spike_data.pkl', 'corr_data.pkl']
+for file_name in files:
+    print(file_name)
+    file = open(file_name, 'rb')
+    data = pkl.load(file)
+    print(data.keys())
+    for key, value in data.items():
+        globals()[key] = value
+
+    file.close()
+
+def pkl_save(filename:str, object:object):
+    file = open(filename+'.pkl', 'wb')
+    pkl.dump(object, file, pkl.HIGHEST_PROTOCOL)
+    file.close()
+
+def pkl_load(filename:str):
+    file = open(filename+'.pkl', 'rb')
+    a = pkl.load(file)
+    file.close()
+    return a

model/sym_model.py

@@ -1,9 +1,9 @@
 from __future__ import annotations
 import typing
 from dataclasses import dataclass, field
+from itertools import chain
 from functools import partial
 
-import matplotlib as mpl
 import sympy as sp
 import numpy as np
 
@@ -143,20 +143,20 @@ class Population:
                 )
 
     @classmethod
-    def random(cls, n: int, 
+    def sample(cls, n_orient: int, n_phases: int,
                phi_dist: Distribution = phi_dist_uni,
                theta_dist: Distribution = theta_dist_uni,
                sigma_dist: Distribution = defaults[sigma],
                k_val: float = defaults[k],
                xy_dist: Distribution = get_uniform_dist(-5, 5)):
         return cls(cells=[
-            Cell(phi_val=phi_val, theta_val=theta_val, sigma_val=sigma_val, x0_val=x0_val, y0_val=y0_val, k_val=k_val) 
+            Cell(phi_val=phi_val, theta_val=theta_val, sigma_val=sigma_val, x0_val=x0_val * 0, y0_val=y0_val * 0, k_val=k_val) 
             for phi_val, theta_val, sigma_val, x0_val, y0_val in zip(
-                sample_distribution(phi_dist, n),
+                list(sample_distribution(phi_dist, n_phases)) * n_orient,
-                sample_distribution(theta_dist, n),
+                list(chain(*[[x] * n_phases for x in sample_distribution(theta_dist, n_orient)])),
-                sample_distribution(sigma_dist, n),
+                sample_distribution(sigma_dist, n_phases * n_orient),
-                sample_distribution(xy_dist, n),
+                sample_distribution(xy_dist, n_phases * n_orient),
-                sample_distribution(xy_dist, n))
+                sample_distribution(xy_dist, n_phases * n_orient))
         ])
 
     def get_response(self, phi_deg: float, theta_deg: float, coef: float = 4, use_sigmoid: bool = True) -> np.ndarray: