floss dance

code/spikes_analysis.py

@@ -1,5 +1,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
+import scipy.stats as ss
 from read_chirp_data import *
 from utility import *
 from IPython import embed
@@ -9,14 +10,15 @@ sampling_rate = 40 #kHz
 data_dir = "../data"
 dataset = "2018-11-09-ad-invivo-1"
 # parameters for binning, smoothing and plotting
+cut_window = 60
 chirp_size = 14 #ms
 neuronal_delay = 5 #ms
-chirp_start = int((-chirp_size/2+neuronal_delay+50)*sampling_rate)
-chirp_end = int((chirp_size/2+neuronal_delay+51)*sampling_rate)
+chirp_start = int((-chirp_size/2+neuronal_delay+cut_window)*sampling_rate)
+chirp_end = int((chirp_size/2+neuronal_delay+cut_window+1)*sampling_rate)
 num_bin = 12
 window = 1 #ms
-time_axis = np.arange(-50, 50, 1/sampling_rate) #steps
-spike_bins = np.arange(-50, 51) #ms
+time_axis = np.arange(-cut_window, cut_window, 1/sampling_rate) #steps
+spike_bins = np.arange(-cut_window, cut_window+1) #ms
 
 # read data from files
 spikes = read_chirp_spikes(os.path.join(data_dir, dataset))
@@ -24,17 +26,11 @@ eod = read_chirp_eod(os.path.join(data_dir, dataset))
 chirp_times = read_chirp_times(os.path.join(data_dir, dataset))
 
 # make a delta f map for the quite more complicated keys
-df_map = {}
-for k in spikes.keys():
-    df = k[1]
-    if df in df_map.keys():
-        df_map[df].append(k)
-    else:
-        df_map[df] = [k]
+df_map = map_keys(spikes)
 
 # differentiate between phases
 phase_vec = np.arange(0, 1+1/num_bin, 1/num_bin)
-cut_range = np.arange(-50*sampling_rate, 50*sampling_rate, 1)
+cut_range = np.arange(-cut_window*sampling_rate, cut_window*sampling_rate, 1)
 
 # make dictionaries for spiketimes
 df_phase_time = {}
@@ -52,7 +48,7 @@ for deltaf in df_map.keys():
 
                     # get spikes between 50 ms befor and after the chirp
                     spikes_to_cut = np.asarray(spikes[rep][phase])
-                    spikes_cut = spikes_to_cut[(spikes_to_cut > -50) & (spikes_to_cut < 50)]
+                    spikes_cut = spikes_to_cut[(spikes_to_cut > -cut_window) & (spikes_to_cut < cut_window)]
                     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
@@ -65,28 +61,47 @@ for deltaf in df_map.keys():
                         df_phase_time[deltaf][idx] = [spikes_cut]
                         df_phase_binary[deltaf][idx] = binary_spikes
 
+
 # for plotting and calculating iterate over delta f and phases
 for df in df_phase_time.keys():
+    beat_duration = int(abs(1/df*1000)) #ms
+    beat_window = 0
+    while beat_window+beat_duration <= cut_window:
+        beat_window = beat_window+beat_duration
     for phase in df_phase_time[df].keys():
 
         trials_binary = df_phase_binary[df][phase]
 
-        sr_chirp = np.zeros(len(trials_binary))
-        sr_beat = np.zeros(len(trials_binary))
+        train_chirp = []
+        train_beat = []
+        spikerate_chirp = np.zeros(len(trials_binary))
+        spikerate_beat = np.zeros(len(trials_binary))
         for i, trial in enumerate(trials_binary):
             smoothed_trial = smooth(trial, window, 1/sampling_rate)
-            sr_chirp[i] = np.mean(smoothed_trial[chirp_start:chirp_end])
-            sr_beat[i] = np.mean(smoothed_trial[0:chirp_start])
-
-        for rate_chirp in sr_chirp:
-            for rate_beat in sr_beat:
-                r = np.corrcoef(rate_chirp, rate_beat)
-                print(r)
-
-        embed()
-        exit()
-
-        #csi = (spikerate_chirp-spikerate_befor)/(spikerate_chirp+spikerate_befor)
+            train_chirp.append(smoothed_trial[chirp_start:chirp_end])
+            train_beat.append(smoothed_trial[chirp_start-beat_window:chirp_start])
+            spikerate_chirp[i] = np.mean(smoothed_trial[chirp_start:chirp_end])
+            spikerate_beat[i] = np.mean(smoothed_trial[chirp_start-beat_window:chirp_start])
+
+        rcs = []
+        rbs = []
+        for i, train in enumerate(train_chirp):
+            for j, train2 in enumerate(train_chirp):
+                if np.array_equal(train, train2):
+                    continue
+                else:
+                    rc, _ = ss.pearsonr(train, train2)
+                    rb, _ = ss.pearsonr(train_beat[i], train_beat[j])
+                    rcs.append(rc)
+                    rbs.append(rb)
+
+        r_train_chirp = np.mean(rcs)
+        r_train_beat = np.mean(rbs)
+
+        csi_train = (r_train_chirp - r_train_beat) / (r_train_chirp + r_train_beat)
+        print(csi_train)
+        csi_rate = (np.std(spikerate_chirp) - np.std(spikerate_beat)) / (np.std(spikerate_chirp) + np.std(spikerate_beat))
+        print(csi_rate)
 
         # plot
         #plot_trials = df_phase_time[df][phase]