Merge branch 'master' of https://whale.am28.uni-tuebingen.de/git/jgrewe/gp_neurobio

code/base_chirps.py

@@ -8,7 +8,7 @@ from IPython import embed
 
 
 data_dir = "../data"
-dataset = "2018-11-13-ah-invivo-1"
+dataset = "2018-11-13-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", "2018-11-13-ac-invivo-1", "2018-11-13-ad-invivo-1", "2018-11-13-ah-invivo-1", "2018-11-13-ai-invivo-1", "2018-11-13-aj-invivo-1", "2018-11-13-ak-invivo-1", "2018-11-13-al-invivo-1", "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", "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"]
 
 
@@ -41,7 +41,9 @@ for i in sort_df:
 chirp_spikes = read_chirp_spikes(os.path.join(data_dir, dataset))
 df_map = map_keys(chirp_spikes)
 sort_df = sorted(df_map.keys())
-dct_phase = plot_std_chirp(sort_df, df_map, chirp_spikes, chirp_mods)
+example = [-50, 200, 400]
+
+dct_phase = plot_std_chirp(example, df_map, chirp_spikes, chirp_mods)
 plt.show()
 plt.close('all')
 

code/base_eod.py

@@ -14,10 +14,16 @@ time,eod = read_baseline_eod(os.path.join(data_dir, dataset))
 zeit = np.asarray(time)
 
 
-plt.plot(zeit[0:1000], eod[0:1000])
+
-plt.title('A.lepto EOD', fontsize = 18)#Plottitelk
+inch_factor = 2.54
-plt.xlabel('time [ms]', fontsize = 16)#Achsentitel
+fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
-plt.ylabel('amplitude[mv]', fontsize = 16)#Achsentitel
+
-plt.xticks(fontsize = 14)
+plt.plot(zeit[0:1000], eod[0:1000], color = 'darkblue')
-plt.yticks(fontsize = 14)
+plt.title('A.lepto EOD', fontsize = 24)#Plottitel
+plt.xlabel('time [ms]', fontsize = 22)#Achsentitel
+plt.ylabel('amplitude[mv]', fontsize = 22)#Achsentitel
+plt.tick_params(axis='both', which='major', labelsize = 22)
+ax.spines['right'].set_visible(False)
+ax.spines['top'].set_visible(False)
+fig.tight_layout()
 plt.show()

code/base_spikes.py

@@ -12,25 +12,28 @@ data_base = ("2018-11-09-ab-invivo-1", "2018-11-09-ad-invivo-1", "2018-11-13-aa-
 data_chirps = ("2018-11-09-ad-invivo-1", "2018-11-09-ae-invivo-1", "2018-11-09-ag-invivo-1", "2018-11-13-aa-invivo-1", "2018-11-13-ac-invivo-1", "2018-11-13-ad-invivo-1", "2018-11-13-ah-invivo-1", "2018-11-13-ai-invivo-1", "2018-11-13-aj-invivo-1", "2018-11-13-ak-invivo-1", "2018-11-13-al-invivo-1", "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", "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") 
 
 dataset = "2018-11-13-ad-invivo-1"
-
+inch_factor = 2.54
 
 #for dataset in data_base:
 
 spike_times = read_baseline_spikes(os.path.join(data_dir, dataset))	
 spike_iv = np.diff(spike_times)
 x = np.arange(0.001, 0.01, 0.0001)
-plt.hist(spike_iv,x)
+
+fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
+plt.hist(spike_iv,x, color = 'darkblue')
 
 mu = np.mean(spike_iv)
 sigma = np.std(spike_iv)
 cv = sigma/mu
 
-plt.title('A.lepto ISI Histogramm', fontsize = 18)
+plt.title('A.lepto ISI Histogramm', fontsize = 24)
-plt.xlabel('duration ISI[ms]', fontsize = 16)
+plt.xlabel('duration ISI[ms]', fontsize = 22)
-plt.ylabel('number of ISI', fontsize = 16)
+plt.ylabel('number of ISI', fontsize = 22)
-
+plt.tick_params(axis='both', which='major', labelsize = 22)
-plt.xticks(fontsize = 14)
+ax.spines['right'].set_visible(False)
-plt.yticks(fontsize = 14)
+ax.spines['top'].set_visible(False)
+plt.tight_layout()
 plt.show()
 
 
@@ -47,6 +50,7 @@ sort_df = sorted(df_map.keys())
 dct_rate, over_r = spike_rates(sort_df, df_map, chirp_spikes)
 
 plt.figure()
+fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 ls_mean = plot_df_spikes(sort_df, dct_rate)
 plt.show()
 
@@ -54,14 +58,17 @@ plt.show()
 
 	#mittlere Feuerrate einer Frequenz auf Frequenz:
 
-plt.figure()
+fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
-plt.plot(np.arange(0,len(ls_mean),1),ls_mean)
+plt.plot(np.arange(0,len(ls_mean),1),ls_mean, color = 'darkblue')
 plt.scatter(np.arange(0,len(ls_mean),1), np.ones(len(ls_mean))*over_r, color = 'green')
-plt.title('Mean firing rate of a cell for a range of frequency differences', fontsize = 18)
+plt.title('Mean firing rate of a cell for a range of frequency differences', fontsize = 24)
 plt.xticks(np.arange(1,len(sort_df),1), (sort_df))
-plt.xlabel('Range of frequency differences [Hz]', fontsize = 16)
+plt.xlabel('Range of frequency differences [Hz]', fontsize = 22)
-plt.ylabel('Mean firing rate of the cell', fontsize = 16)
+plt.ylabel('Mean firing rate of the cell', fontsize = 22)
-plt.tick_params(axis='both', which='major', labelsize = 14)
+plt.tick_params(axis='both', which='major', labelsize = 18)
+ax.spines['right'].set_visible(False)
+ax.spines['top'].set_visible(False)
+plt.tight_layout()
 plt.show()
 
 
@@ -70,12 +77,15 @@ plt.show()
 	#wie viel Prozent der Anfangsrate macht die Adaption von Zellen aus?
 
 adapt = adaptation_df(sort_df, dct_rate)
-plt.figure()
+fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 plt.boxplot(adapt)
-plt.title('Adaptation of cell firing rate during a trial', fontsize = 18)
+plt.title('Adaptation of cell firing rate during a trial', fontsize = 24)
-plt.xlabel('Cell', fontsize = 16)
+plt.xlabel('Cell', fontsize = 22)
-plt.ylabel('Adaptation size [Hz]', fontsize = 16)
+plt.ylabel('Adaptation size [Hz]', fontsize = 22)
-plt.tick_params(axis='both', which='major', labelsize = 14)
+plt.tick_params(axis='both', which='major', labelsize = 18)
+ax.spines['right'].set_visible(False)
+ax.spines['top'].set_visible(False)
+plt.tight_layout()
 plt.show()	
 
 

code/func_chirp.py

@@ -2,8 +2,10 @@ from read_baseline_data import *
 from read_chirp_data import *
 from utility import *
 import matplotlib.pyplot as plt
+import math
 import numpy as np
 
+inch_factor = 2.54
 
 def chirp_eod_plot(df_map, eod, times):
 	#die äußere Schleife geht für alle Keys durch und somit durch alle dfs
@@ -11,7 +13,7 @@ def chirp_eod_plot(df_map, eod, times):
 
 	for i in df_map.keys():
 		freq = list(df_map[i]) 
-		fig,axs = plt.subplots(2, 2, sharex = True, sharey = True)
+		fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 
 		for idx, k in enumerate(freq):
 			ct = times[k]
@@ -19,22 +21,28 @@ def chirp_eod_plot(df_map, eod, times):
 			zeit = e1[0]
 			eods = e1[1]
 	
-			if idx <= 3:
+			if idx <= 1:
-				axs[0, 0].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+				ax.plot(zeit, eods, color= 'darkblue')
-				axs[0, 0].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
+				ax.scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
-			elif 4<= idx <= 7:
+			#elif 4<= idx <= 7:
-				axs[0, 1].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+			#	axs[0, 1].plot(zeit, eods, color= 'blue', linewidth = 0.25)
-				axs[0, 1].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
+			#	axs[0, 1].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
-			elif 8<= idx <= 11:
+			#elif 8<= idx <= 11:
-				axs[1, 0].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+			#	axs[1, 0].plot(zeit, eods, color= 'blue', linewidth = 0.25)
-				axs[1, 0].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
+			#	axs[1, 0].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
-			else: 
+			else: 	
-				axs[1, 1].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+				continue
-				axs[1, 1].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
+				#axs[1, 1].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+				#axs[1, 1].scatter(np.asarray(ct), np.ones(len(ct))*np.mean(eods), color = 'green', s= 22)
+
+		ax.set_ylabel('Amplitude [mV]', fontsize = 22)
+		ax.set_xlabel('Time [ms]', fontsize = 22)		
+		ax.tick_params(axis='both', which='major', labelsize = 18)
+		ax.spines['right'].set_visible(False)
+		ax.spines['top'].set_visible(False)
+		fig.suptitle('EOD for chirps', fontsize = 24)
+		fig.tight_layout()
 
-	axs[0,1].set_ylabel('Amplitude [mV]')
-	axs[1,0].set_xlabel('Time [ms]')			
-	fig.suptitle('EOD for chirps', fontsize = 16)
 
 
 
@@ -69,7 +77,7 @@ def cut_chirps(freq, eod, times):
 
 def plot_std_chirp(sort_df, df_map, chirp_spikes, chirp_mods):
 		
-	plt.figure()	
+	fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 	dct_phase = {}
 	num_bin = 12
 	phase_vec = np.arange(0, 1+1/num_bin, 1/num_bin)
@@ -82,10 +90,13 @@ def plot_std_chirp(sort_df, df_map, chirp_spikes, chirp_mods):
 				dct_phase[i].append(phase[1])
 	
 	for i in sort_df:
-		plt.scatter(dct_phase[i], chirp_mods[i], label = i)
+		norm = np.asarray(dct_phase[i]) *2*math.pi
-	plt.title('Change of std depending on the phase where the chirp occured')	
+		plt.scatter(norm, chirp_mods[i], label = i, s = 22)
-	plt.xlabel('Phase')
+	plt.title('Change of std depending on the phase where the chirp occured', fontsize = 24)	
-	plt.ylabel('Standard deviation of the amplitude modulation')
+	plt.xlabel('Phase', fontsize = 22)
+	plt.ylabel('Standard deviation of the amplitude modulation', fontsize = 22)
+	plt.xticks([0, math.pi/2, math.pi, math.pi*1.5, math.pi*2], ('0', '$\pi$ /2', '$\pi$', '1.5 $\pi$', '2$\pi$'))
+	plt.tick_params(axis='both', which='major', labelsize = 18)
 	plt.legend()
 	return(dct_phase)
 

code/func_spike.py

@@ -43,6 +43,8 @@ def spike_rates(sort_df, df_map, chirp_spikes):
 
 def plot_df_spikes(sort_df, dct_rate):
 #gibt die Feuerrate gegen die Frequenz aufgetragen
+	inch_factor = 2.54
+	fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 	ls_mean = []
 	for h in sort_df:	
 		mean = np.mean(dct_rate[h])
@@ -50,10 +52,13 @@ def plot_df_spikes(sort_df, dct_rate):
 		plt.plot(np.arange(0,len(dct_rate[h]),1),dct_rate[h], label = h)
 
 	plt.legend()
-	plt.title('Firing rate of the cell for all trials, sorted by df', fontsize = 18)
+	plt.title('Firing rate of the cell for all trials, sorted by df', fontsize = 24)
-	plt.xlabel('# of trials', fontsize = 16)
+	plt.xlabel('# of trials', fontsize = 22)
-	plt.ylabel('Instant firing rate of the cell', fontsize = 16)
+	plt.ylabel('Instant firing rate of the cell', fontsize = 22)
-	plt.tick_params(axis='both', which='major', labelsize = 14)
+	plt.tick_params(axis='both', which='major', labelsize = 18)
+	ax.spines['right'].set_visible(False)
+	ax.spines['top'].set_visible(False)
+	plt.tight_layout()
 	return(ls_mean)
 
 

code/order_eff.py

@@ -12,6 +12,7 @@ data_chirps = ("2018-11-09-ad-invivo-1", "2018-11-09-ae-invivo-1", "2018-11-09-a
 
 
 
+inch_factor = 2.54
 data_rate_dict = {}
 for dataset in data_chirps:
 	
@@ -31,14 +32,17 @@ for dataset in data_chirps:
 		rate = len(spikes)/ 1.2
 		data_rate_dict[dataset].append(rate)		
 
-
+fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 for dataset in data_rate_dict:
 	plt.plot(data_rate_dict[dataset])
-plt.title('Test for sequence effects', fontsize = 20)
-plt.xlabel('Number of stimulus presentations', fontsize = 18)
-plt.ylabel('Firing rates of cells', fontsize = 18)
-plt.tick_params(axis='both', which='major', labelsize = 16)
 
+plt.title('Test for sequence effects', fontsize = 24)
+plt.xlabel('Number of stimulus presentations', fontsize = 22)
+plt.ylabel('Firing rates of cells', fontsize = 22)
+plt.tick_params(axis='both', which='major', labelsize = 22)
+ax.spines['right'].set_visible(False)
+ax.spines['top'].set_visible(False)
+fig.tight_layout()
 plt.show()
 
 

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/plot_spikesduringbaselineactivity.py

@@ -17,6 +17,7 @@ interspikeintervals = np.diff(spikes)*1000
 
 fig, ax = plt.subplots(figsize=(20/inch_factor, 10/inch_factor))
 plt.hist(interspikeintervals, bins=np.arange(0, np.max(interspikeintervals), 0.0001), color='darkblue')
+#Titel fehlt!!
 plt.xlabel("time [ms]", fontsize = 22)
 plt.xticks(fontsize = 18)
 plt.ylabel("Number of \n Interspikeinterval", fontsize = 22)

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()