wuhu

code/plot_eodform_spikehist.py

@@ -9,8 +9,8 @@ from IPython import embed
 inch_factor = 2.54
 sampling_rate = 40000
 data_dir = '../data'
-dataset = '2018-11-09-ad-invivo-1'
+#dataset = '2018-11-09-ad-invivo-1'
-#dataset = '2018-11-13-aa-invivo-1'
+dataset = '2018-11-14-ad-invivo-1'
 
 # read eod and time of baseline
 time, eod = read_baseline_eod(os.path.join(data_dir, dataset))

code/response_beat.py

@@ -8,29 +8,33 @@ from IPython import embed
 # define data path and important parameters
 data_dir = "../data"
 sampling_rate = 40 #kHz
-cut_window = 40
+cut_window = 100
 cut_range = np.arange(-cut_window * sampling_rate, 0, 1)
 window = 1
 
+'''
 # norm: -150, 150, 300  aa, #ac, aj??
-data = ["2018-11-13-aa-invivo-1"]#, "2018-11-13-ad-invivo-1", "2018-11-13-ah-invivo-1", "2018-11-13-ai-invivo-1",
+data = ["2018-11-13-aa-invivo-1", "2018-11-13-ad-invivo-1", "2018-11-13-ah-invivo-1", "2018-11-13-ai-invivo-1",
-        #"2018-11-13-ak-invivo-1", "2018-11-13-al-invivo-1"]
+        "2018-11-13-ak-invivo-1", "2018-11-13-al-invivo-1"]
 
-'''
 # norm: -50
 data = ["2018-11-20-aa-invivo-1", "2018-11-20-ab-invivo-1", "2018-11-20-ac-invivo-1","2018-11-20-ad-invivo-1", 
         "2018-11-20-ae-invivo-1", "2018-11-20-af-invivo-1", "2018-11-20-ag-invivo-1", "2018-11-20-ah-invivo-1",
         "2018-11-20-ai-invivo-1"]
-
-data = ["2018-11-14-aa-invivo-1", "2018-11-14-ac-invivo-1", "2018-11-14-ad-invivo-1", "2018-11-14-af-invivo-1",
-        "2018-11-14-ag-invivo-1", "2018-11-14-ah-invivo-1", "2018-11-14-ai-invivo-1", "2018-11-14-ak-invivo-1",
-        "2018-11-14-al-invivo-1", "2018-11-14-am-invivo-1", "2018-11-14-an-invivo-1"]
 '''
+
+data = ["2018-11-14-ad-invivo-1", "2018-11-14-ak-invivo-1", "2018-11-14-am-invivo-1"]
+
+#data = ["2018-11-14-ac-invivo-1", "2018-11-14-ad-invivo-1", "2018-11-14-ak-invivo-1"]
 #data = ["2018-11-09-ad-invivo-1", "2018-11-14-af-invivo-1"]
 
+#data = ["2018-11-20-ad-invivo-1", "2018-11-13-ad-invivo-1"]
+#data = ["2018-11-09-ad-invivo-1"]
+
 rates = {}
 
 for dataset in data:
+    rates[dataset] = {}
     print(dataset)
     # read baseline spikes
     base_spikes = read_baseline_spikes(os.path.join(data_dir, dataset))
@@ -66,21 +70,48 @@ for dataset in data:
                 # also save as binary, 0 no spike, 1 spike
                 binary_spikes = np.isin(cut_range, spikes_idx) * 1
                 smoothed_data = smooth(binary_spikes, window, 1 / sampling_rate)
-                train = smoothed_data[window:beat_window+window]
+                #train = smoothed_data[window:beat_window+window]
-                norm_train = train*1000/spikerate
+                #norm_train = train*1000/spikerate
-                rep_rates.append(np.std(norm_train))#/spikerate)
+                #df_rate = np.std(norm_train)
+                #rates[dataset][df] = [df_rate]
+                rep_rates.append(np.std(smoothed_data))#/spikerate)
+                '''
+                if df in rates[dataset].keys():
+                    rates[dataset][df].append(np.std(norm_train))
+                else:
+                    rates[dataset][df] = [np.std(norm_train)]
+                '''
                 break
+            #break
         df_rate = np.mean(rep_rates)
+        #df_rate = rep_rates
+        rates[dataset][df] = df_rate
+
         #embed()
         #exit()
+        '''
         if df in rates.keys():
-            rates[df].append(df_rate)
+            rates[dataset][df].append(df_rate)
         else:
-            rates[df] = [df_rate]
+            rates[dataset][df] = [df_rate]
+        '''
+
+colors = ['royalblue', 'red', 'green', 'violet', 'orange', 'black', 'gray']
+
+fig, ax = plt.subplots()
+for i, cell in enumerate(rates.keys()):
+    for j, df in enumerate(sorted(rates[cell].keys())):
+        ax.plot(df, rates[cell][df], 'o', color=colors[i])
+#ax.legend(sorted(rates.keys()), loc='upper left', bbox_to_anchor=(1.04, 1))
+fig.tight_layout()
+plt.show()
 
+'''
 fig, ax = plt.subplots()
-for i, k in enumerate(sorted(rates.keys())):
+for i, cell in enumerate(rates.keys()):
-    ax.plot(np.ones(len(rates[k]))*k, rates[k], 'o')
+    for j, df in enumerate(sorted(rates[cell].keys())):
+        ax.plot(np.ones(len(rates[cell][df]))*df, rates[cell][df], 'o', color=colors[i])
 #ax.legend(sorted(rates.keys()), loc='upper left', bbox_to_anchor=(1.04, 1))
 fig.tight_layout()
 plt.show()
+'''

code/spikes_analysis.py

@@ -8,7 +8,8 @@ from IPython import embed
 # define sampling rate and data path
 sampling_rate = 40 #kHz
 data_dir = "../data"
-dataset = "2018-11-13-ah-invivo-1"
+dataset = "2018-11-13-al-invivo-1"
+#dataset = "2018-11-09-ad-invivo-1"
 
 '''
 data = ["2018-11-09-ad-invivo-1", "2018-11-09-ae-invivo-1", "2018-11-09-ag-invivo-1", "2018-11-13-aa-invivo-1",
@@ -90,14 +91,14 @@ for deltaf in df_map.keys():
 
 
 # make dictionaries for csi and beat
-csi_trains = {}
+#csi_trains = {}
 csi_rates = {}
-beat = {}
+#beat = {}
 # for plotting and calculating iterate over delta f and phases
 for df in df_phase_time.keys():
-    csi_trains[df] = []
+    #csi_trains[df] = {}
-    csi_rates[df] = []
+    csi_rates[df] = {}
-    beat[df] = []
+    #beat[df] = []
     beat_duration = int(abs(1/df*1000)*sampling_rate) #steps
     beat_window = 0
     # beat window is at most 20 ms long, multiples of beat_duration
@@ -123,9 +124,12 @@ for df in df_phase_time.keys():
 
         std_chirp = np.std(np.mean(train_chirp, axis=0))
         std_beat = np.std(np.mean(train_beat, axis=0))
-        beat[df].append(std_beat)
+        #beat[df].append(std_beat)
         csi_spikerate = (std_chirp - std_beat) / (std_chirp + std_beat)
 
+        csi_rates[df][phase] = np.mean(csi_spikerate)
+
+        '''
         rcs = []
         rbs = []
         for i, train in enumerate(train_chirp):
@@ -145,9 +149,7 @@ for df in df_phase_time.keys():
 
         # add the csi to the dictionaries with the correct df and phase
         csi_trains[df].append(csi_train)
-        csi_rates[df].append(np.mean(csi_spikerate))
 
-        '''
         # plot
         plot_trials = df_phase_time[df][phase]
         plot_trials_binary = np.mean(df_phase_binary[df][phase], axis=0)
@@ -170,37 +172,39 @@ for df in df_phase_time.keys():
         plt.show()
         '''
 
+
+colors = ['k', 'k', 'k',
+          'k', 'k', 'k',
+          'k', 'k', 'k',
+          'k', 'k', 'firebrick']
+
+sizes = [12, 12, 12,
+         12, 12, 12,
+         12, 12, 12,
+         12, 12, 18]
+
 upper_limit = np.max(sorted(csi_rates.keys()))+30
 lower_limit = np.min(sorted(csi_rates.keys()))-30
 
 fig, ax = plt.subplots()
-for i, k in enumerate(sorted(csi_rates.keys())):
-    ax.scatter(np.ones(len(csi_rates[k]))*k, csi_rates[k], s=20)
-    #ax.plot(i, np.mean(csi_rates[k]), 'o', markersize=15)
-#ax.legend(sorted(csi_rates.keys()), loc='upper left', bbox_to_anchor=(1.04, 1))
 ax.plot([lower_limit, upper_limit], np.zeros(2), 'silver', linewidth=2, linestyle='--')
-#ax.set_xticklabels(sorted(csi_rates.keys()))
+for i, df in enumerate(sorted(csi_rates.keys())):
+    for j, phase in enumerate(sorted(csi_rates[df].keys())):
+        ax.plot(df, csi_rates[df][phase], 'o', color=colors[j], ms=sizes[j])
 fig.tight_layout()
 plt.show()
 
 '''
 fig, ax = plt.subplots()
-for i, k in enumerate(sorted(csi_trains.keys())):
+for i, k in enumerate(sorted(beat.keys())):
-    ax.plot(np.ones(len(csi_trains[k]))*i, csi_trains[k], 'o')
+    ax.plot(np.ones(len(beat[k]))*k, beat[k], 'o')
-    #ax.plot(i, np.mean(csi_trains[k]), 'o', markersize=15)
-ax.legend(sorted(csi_trains.keys()), loc='upper left', bbox_to_anchor=(1.04, 1))
-ax.plot(np.arange(-1, len(csi_trains.keys())+1), np.zeros(len(csi_trains.keys())+2), 'silver', linewidth=2, linestyle='--')
-#ax.set_xticklabels(sorted(csi_trains.keys()))
 fig.tight_layout()
 plt.show()
-'''
 
-'''
 fig, ax = plt.subplots()
-for i, k in enumerate(sorted(beat.keys())):
+for i, k in enumerate(sorted(csi_trains.keys())):
-    ax.plot(np.ones(len(beat[k]))*i, beat[k], 'o')
+    ax.plot(np.ones(len(csi_trains[k]))*i, csi_trains[k], 'o')
-ax.legend(sorted(beat.keys()), loc='upper left', bbox_to_anchor=(1.04, 1))
+ax.plot(np.arange(-1, len(csi_trains.keys())+1), np.zeros(len(csi_trains.keys())+2), 'silver', linewidth=2, linestyle='--')
-#ax.set_xticklabels(sorted(csi_trains.keys()))
 fig.tight_layout()
 plt.show()
 '''

code/spikes_beat.py

@@ -0,0 +1,62 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from read_chirp_data import *
+from read_baseline_data import *
+from utility import *
+from IPython import embed
+
+# define data path and important parameters
+data_dir = "../data"
+sampling_rate = 40 #kHz
+cut_window = 100
+cut_range = np.arange(-cut_window * sampling_rate, 0, 1)
+window = 1
+#dataset = "2018-11-13-ad-invivo-1"
+#dataset = "2018-11-13-aj-invivo-1"
+#dataset = "2018-11-13-ak-invivo-1" #al
+#dataset = "2018-11-14-ad-invivo-1"
+dataset = "2018-11-20-af-invivo-1"
+
+base_spikes = read_baseline_spikes(os.path.join(data_dir, dataset))
+base_spikes = base_spikes[1000:2000]
+spikerate = len(base_spikes) / base_spikes[-1]
+print(spikerate)
+
+# read spikes during chirp stimulation
+spikes = read_chirp_spikes(os.path.join(data_dir, dataset))
+df_map = map_keys(spikes)
+
+rates = {}
+# iterate over df
+for deltaf in df_map.keys():
+    rates[deltaf] = {}
+    beat_duration = int(abs(1 / deltaf) * 1000)
+    beat_window = 0
+    while beat_window + beat_duration <= cut_window/2:
+        beat_window = beat_window + beat_duration
+    for x, repetition in enumerate(df_map[deltaf]):
+        for phase in spikes[repetition]:
+            # get spikes some ms before the chirp first chirp
+            spikes_to_cut = np.asarray(spikes[repetition][phase])
+            spikes_cut = spikes_to_cut[(spikes_to_cut > -cut_window) & (spikes_to_cut < 0)]
+            spikes_idx = np.round(spikes_cut * sampling_rate)
+            # also save as binary, 0 no spike, 1 spike
+            binary_spikes = np.isin(cut_range, spikes_idx) * 1
+            smoothed_data = smooth(binary_spikes, window, 1 / sampling_rate)
+            #train = smoothed_data[window*sampling_rate:beat_window*sampling_rate+window*sampling_rate]
+            modulation = np.std(smoothed_data)
+            rates[deltaf][x] = modulation
+            break
+
+fig, ax = plt.subplots()
+for i, df in enumerate(sorted(rates.keys())):
+    for j, rep in enumerate(rates[df].keys()):
+        if j == 15:
+            farbe = 'royalblue'
+            gro = 18
+        else:
+            farbe = 'k'
+            gro = 12
+        ax.plot(df, rates[df][rep], marker='o', color=farbe, ms=gro)
+fig.tight_layout()
+plt.show()

code/spikes_chirp.py

@@ -0,0 +1,102 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from read_chirp_data import *
+from utility import *
+from IPython import embed
+
+# define sampling rate and data path
+sampling_rate = 40 #kHz
+data_dir = "../data"
+dataset = "2018-11-13-al-invivo-1"
+
+# parameters for binning, smoothing and plotting
+cut_window = 20
+chirp_duration = 14 #ms
+neuronal_delay = 5 #ms
+chirp_start = int((-chirp_duration / 2 + neuronal_delay + cut_window * 2) * sampling_rate) #index
+chirp_end = int((chirp_duration / 2 + neuronal_delay + cut_window * 2) * sampling_rate) #index
+number_bins = 12
+window = 1 #ms
+time_axis = np.arange(-cut_window*2, cut_window*2, 1/sampling_rate) #steps
+spike_bins = np.arange(-cut_window*2, cut_window*2) #ms
+
+colors = ['k', 'k', 'k',
+          'k', 'k', 'k',
+          'k', 'k', 'k',
+          'k', 'k', 'firebrick']
+
+sizes = [12, 12, 12,
+         12, 12, 12,
+         12, 12, 12,
+         12, 12, 18]
+
+# differentiate between phases
+phase_vec = np.arange(0, 1 + 1 / number_bins, 1 / number_bins)
+cut_range = np.arange(-cut_window*2*sampling_rate, cut_window*2*sampling_rate, 1)
+
+df_phase_binary = {}
+
+spikes = read_chirp_spikes(os.path.join(data_dir, dataset))
+df_map = map_keys(spikes)
+
+for deltaf in df_map.keys():
+    df_phase_binary[deltaf] = {}
+    for rep in df_map[deltaf]:
+        chirp_size = int(rep[-1].strip('Hz'))
+        if chirp_size == 150:
+            continue
+        for phase in spikes[rep]:
+            for idx in np.arange(number_bins):
+                # check the phase
+                if phase[1] > phase_vec[idx] and phase[1] < phase_vec[idx+1]:
+
+                    # get spikes between 40 ms before and after the chirp
+                    spikes_to_cut = np.asarray(spikes[rep][phase])
+                    spikes_cut = spikes_to_cut[(spikes_to_cut > -cut_window*2) & (spikes_to_cut < cut_window*2)]
+                    spikes_idx = np.round(spikes_cut*sampling_rate)
+                    # save as binary, 0 no spike, 1 spike
+                    binary_spikes = np.isin(cut_range, spikes_idx)*1
+
+                    # add the spikes to the dictionary with the correct df and phase
+                    if idx in df_phase_binary[deltaf].keys():
+                        df_phase_binary[deltaf][idx] = np.vstack((df_phase_binary[deltaf][idx], binary_spikes))
+                    else:
+                        df_phase_binary[deltaf][idx] = binary_spikes
+
+csi_rates = {}
+
+for df in df_phase_binary.keys():
+    csi_rates[df] = {}
+    beat_duration = int(abs(1/df*1000)*sampling_rate) #steps
+    beat_window = 0
+    # beat window is at most 20 ms long, multiples of beat_duration
+    while beat_window+beat_duration <= cut_window*sampling_rate:
+        beat_window = beat_window+beat_duration
+    for phase in df_phase_binary[df].keys():
+        # csi calculation
+        trials_binary = df_phase_binary[df][phase]
+
+        train_chirp = []
+        train_beat = []
+        for i, trial in enumerate(trials_binary):
+            smoothed_trial = smooth(trial, window, 1/sampling_rate)
+            train_chirp.append(smoothed_trial[chirp_start:chirp_end])
+            train_beat.append(smoothed_trial[chirp_start-beat_window:chirp_start])
+
+        std_chirp = np.std(np.mean(train_chirp, axis=0))
+        std_beat = np.std(np.mean(train_beat, axis=0))
+        csi_spikerate = (std_chirp - std_beat) / (std_chirp + std_beat)
+
+        csi_rates[df][phase] = np.mean(csi_spikerate)
+
+
+upper_limit = np.max(sorted(csi_rates.keys()))+30
+lower_limit = np.min(sorted(csi_rates.keys()))-30
+
+fig, ax = plt.subplots()
+ax.plot([lower_limit, upper_limit], np.zeros(2), 'silver', linewidth=2, linestyle='--')
+for i, df in enumerate(sorted(csi_rates.keys())):
+    for j, phase in enumerate(sorted(csi_rates[df].keys())):
+        ax.plot(df, csi_rates[df][phase], 'o', color=colors[j], ms=sizes[j])
+fig.tight_layout()
+plt.show()