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

code/NixFrame.py

@@ -0,0 +1,188 @@
+import nixio as nix
+from IPython import embed
+import numpy as np
+import os
+import pandas as pd
+import pickle
+
+def DataFrame(nixfile, savefile=False, saveto='./'):
+    '''
+    opens a nix file, extracts the data and converts it to a pandas.DataFrame
+
+    :param nixfile (string): path and name of .nix file
+    :param savefile (string): if not False, the dataframe will be saved as <savefile>.pickle
+    :param saveto (string): path to save the files in NOT IMPLEMENTED YET
+    :return dataframe (pandas.DataFrame): pandas.DataFrame with available nix data
+    '''
+
+    block = nix.File.open(nixfile,'r').blocks[0]
+
+    data_arrays = block.data_arrays
+    names = [data_arrays[i].name for i in range(len(data_arrays))]
+    shapes = [x.shape for x in data_arrays]
+    data_names = np.array([[x,i] for i,x in enumerate(names) if (shapes[i][0] >= 0.999*shapes[0][0])])
+    data_traces = np.array([data_arrays[name][:] for name,idx in data_names])
+    time = data_arrays[1].dimensions[0].axis(data_arrays[1].shape[0])
+    dt = time[1]-time[0]
+
+    block_metadata = {}
+    block_metadata[block.id] = getMetadataDict(block.metadata)
+
+    tag = block.tags
+    tag_metadata = {}
+    tag_id_times = {}
+    for i in range(len(tag)):
+        meta = tag[i].metadata
+        tag_metadata[meta.id] = getMetadataDict(meta)
+        tag_id_times[meta.id] = [tag[i].position[0], tag[i].position[0]+tag[i].extent[0]]
+
+    data = []
+    stim_num = -1
+    protocol_idcs = np.where([' onset times' in name for name in names])[0]
+    for i in range(len(protocol_idcs)):
+        # print(names[int(protocol_idcs[i])].split(' onset times')[0])
+        protocol = names[protocol_idcs[i]].split(' onset times')[0]
+
+        #skip certain protocols
+        if 'VC=' in protocol:
+            # print('skip this protocol')
+            continue
+
+        #number of meta data entries
+        if i == len(protocol_idcs)-1:
+            meta_len = len(names) - protocol_idcs[i]
+        else:
+            meta_len = protocol_idcs[i+1] - protocol_idcs[i]
+
+        #get new line for every sweep and save the data, make a pn subtraction if necessary
+        if any([protocol + '_pn' == string for string in names[protocol_idcs[i]:protocol_idcs[i]+meta_len]]):
+            pn = data_arrays[protocol + '_pn'][0]
+            sweeps = np.arange(np.abs(pn),len(data_arrays[int(protocol_idcs[i])][:]),(np.abs(pn)+1), dtype=int)
+        else:
+            pn = np.nan
+            sweeps = np.arange(len(data_arrays[int(protocol_idcs[i])][:]), dtype=int)
+
+        for sweep in sweeps:
+            stim_num +=1
+            data.append({})
+
+            # save protocol names
+            split_vec = protocol.split('-')
+            if len(split_vec)>2:
+                prot_name = split_vec[0]
+                prot_num = int(split_vec[-1])
+                for j in range(len(split_vec)-2):
+                    prot_name += '-' + split_vec[j+1]
+            else:
+                prot_name = split_vec[0]
+                prot_num = split_vec[-1]
+            data[stim_num]['protocol'] = prot_name
+            data[stim_num]['protocol_number'] = prot_num
+
+            #save id
+            data[stim_num]['id'] = data_arrays[int(protocol_idcs[i])].id
+
+            #save rest of stored data
+            for idx in range(meta_len):
+                j = int(protocol_idcs[i] + idx)
+                if (' durations' in names[j]) or (' onset times' in names[j]):
+                    continue
+                if len(data_arrays[j][sweep]) == 1:
+                    data[stim_num][names[j].split(protocol + '_')[-1]] = data_arrays[j][sweep][0]
+                else:
+                    data[stim_num][names[j].split(protocol+'_')[-1]] = data_arrays[j][sweep]
+            data[stim_num]['samplingrate'] = 1/dt
+
+            #save data arrays
+            onset = data_arrays[protocol + ' onset times'][sweep]
+            dur = data_arrays[protocol + ' durations'][sweep]
+            t0 = int(onset/dt)
+            t1 = int((onset+dur)/dt+1)
+            data[stim_num]['onset time'] = onset
+            data[stim_num]['duration'] = dur
+
+            for name,idx in data_names:
+                data[stim_num][name] = data_traces[int(idx)][t0:t1]
+
+            for j in np.arange(int(idx)+1,protocol_idcs[0]):
+                bool_vec = (data_arrays[names[j]][:]>=onset) & (data_arrays[names[j]][:]<=onset+dur)
+                data[stim_num][names[j]] = np.array(data_arrays[names[j]])[bool_vec]
+
+            data[stim_num]['time'] = time[t0:t1] - data[stim_num]['onset time']
+
+            #pn-subtraction (if necessary)
+            '''
+            change the location of the pn (its already in the metadata, you dont need it as option
+            '''
+            if pn != np.nan and np.abs(pn)>0:
+                pn_curr = np.zeros(len(data[stim_num][name]))
+                idx = np.where(data_names[:,0] == 'Current-1')[0][0]
+                for j in range(int(np.abs(pn))):
+                    onset = data_arrays[protocol + ' onset times'][sweep-j-1]
+                    t0 = int(onset / dt)
+                    t1 = int(onset/dt + len(data[stim_num]['Current-1']))
+                    pn_curr += data_traces[int(idx),t0:t1]
+
+                data[stim_num]['Current-2'] = data[stim_num]['Current-1'] - pn/np.abs(pn)*pn_curr #- data[stim_num][name][0] - pn_curr[0]
+
+
+            '''
+            this one saves the complete metadata in EVERY line
+            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!THINK OF SOMETHING BETTER!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+            '''
+
+            tag_id = None
+            for key in tag_id_times.keys():
+                if (data[stim_num]['onset time'] >= tag_id_times[key][0]) and (data[stim_num]['onset time'] <= tag_id_times[key][1]):
+                    tag_id = key
+            # # save metadata
+            data[stim_num]['block_meta'] = block_metadata[list(block_metadata.keys())[0]]
+            data[stim_num]['tag_meta'] = tag_metadata[tag_id]
+
+            # add block id
+            data[stim_num]['block_id'] = list(block_metadata.keys())[0]
+            data[stim_num]['tag_id'] = tag_id
+
+    data = pd.DataFrame(data)
+    if savefile != False:
+        if savefile == True:
+            savefile = nixfile.split('/')[-1].split('.nix')[0]
+
+        with open(savefile + '_dataframe.pickle', 'wb') as f:
+            pickle.dump(data, f, -1)  # create pickle-files, using the highest pickle-protocol
+    # embed()
+    return data
+
+def NixToFrame(folder):
+    '''
+    searches subfolders of folder to convert .nix files to a pandas dataframe and saves them in the folder
+
+    :param folder: path to folder that contains subfolders of year-month-day-aa style that contain .nix files
+    '''
+    if folder[-1] != '/':
+        folder = folder + '/'
+
+    dirlist = os.listdir(folder)
+    for dir in dirlist:
+        if os.path.isdir(folder + dir):
+            for file in os.listdir(folder+dir):
+
+                if '.nix' in file:
+                    print(file)
+                    DataFrame(folder+dir+'/'+file, True, folder)
+
+def load_data(filename):
+    with open(filename, 'rb') as f:
+        data = pickle.load(f)  # load data with pickle
+    return data
+
+def getMetadataDict(metadata):
+    def unpackMetadata(sec):
+        metadata = dict()
+        metadata = {prop.name: sec[prop.name] for prop in sec.props}
+        if hasattr(sec, 'sections') and len(sec.sections) > 0:
+            metadata.update({subsec.name: unpackMetadata(subsec) for subsec in sec.sections})
+        return metadata
+
+    return unpackMetadata(metadata)
+

code/analysis_rs.py

@@ -1,60 +1,94 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from read_baseline_data import *
+from NixFrame import *
+from utility import *
 from IPython import embed
 
+# plot and data values
 inch_factor = 2.54
 data_dir = '../data'
 dataset = '2018-11-09-ad-invivo-1'
+
+# read eod and time of baseline
 time, eod = read_baseline_eod(os.path.join(data_dir, dataset))
 
-fig = plt.figure(figsize=(12/inch_factor, 8/inch_factor))
+fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
-ax = fig.add_subplot(111)
 ax.plot(time[:1000], eod[:1000])
 ax.set_xlabel('time [ms]', fontsize=12)
 ax.set_ylabel('voltage [mV]', fontsize=12)
-plt.xticks(fontsize = 8)
+plt.xticks(fontsize=8)
-plt.yticks(fontsize = 8)
+plt.yticks(fontsize=8)
 fig.tight_layout()
-plt.savefig('eod.pdf')
+#plt.savefig('eod.pdf')
-
+plt.show()
-#interspikeintervalhistogram, windowsize = 1 ms
-#plt.hist
-#coefficient of variation
-#embed()
-#exit()
 
+# read spikes during baseline activity
 spikes = read_baseline_spikes(os.path.join(data_dir, dataset))
+# calculate interpike intervals and plot them
 interspikeintervals = np.diff(spikes)
-fig = plt.figure()
+
+fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
 plt.hist(interspikeintervals, bins=np.arange(0, np.max(interspikeintervals), 0.0001))
 plt.show()
 
+# calculate coefficient of variation
 mu = np.mean(interspikeintervals)
 sigma = np.std(interspikeintervals)
 cv = sigma/mu
 print(cv)
 
-# calculate zero crossings of the eod
+# calculate eod times and indices by zero crossings
-# plot mean of eod circles
+threshold = 0
-# plot std of eod circles
-# plot psth into the same plot
-# calculate vector strength
-
-threshold = 0;
 shift_eod = np.roll(eod, 1)
 eod_times = time[(eod >= threshold) & (shift_eod < threshold)]
 sampling_rate = 40000.0
 eod_idx = eod_times*sampling_rate
 
-fig = plt.figure()
+# align eods and spikes to eods
-for i, idx in enumerate(eod_idx):
+max_cut = int(np.max(np.diff(eod_idx)))
-	#embed()
+eod_cuts = np.zeros([len(eod_idx)-1, max_cut])
-	#exit()
+spike_times = []
-	plt.plot(time[int(idx):int(eod_idx[i+1])], eod[int(idx):int(eod_idx[i+1])])
+eod_durations = []
+
+for i, idx in enumerate(eod_idx[:-1]):
+	eod_cut = eod[int(idx):int(eod_idx[i+1])]
+	eod_cuts[i, :len(eod_cut)] = eod_cut
+	eod_cuts[i, len(eod_cut):] = np.nan
+	time_cut = time[int(idx):int(eod_idx[i+1])]
+	spike_cut = spikes[(spikes > time_cut[0]) & (spikes < time_cut[-1])]
+	spike_time = spike_cut - time_cut[0]
+	if len(spike_time) > 0:
+		spike_times.append(spike_time[:][0]*1000)
+		eod_durations.append(len(eod_cut)/sampling_rate*1000)
 
+# calculate vector strength
+vs = vector_strength(spike_times, eod_durations)
 
+# determine means and stds of eod for plot
+# determine time axis
+mu_eod = np.nanmean(eod_cuts, axis=0)
+std_eod = np.nanstd(eod_cuts, axis=0)*3
+time_axis = np.arange(max_cut)/sampling_rate*1000
 
+# plot eod form and spike histogram
+fig, ax1 = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
+ax1.hist(spike_times, color='crimson')
+ax1.set_xlabel('time [ms]', fontsize=12)
+ax1.set_ylabel('number', fontsize=12)
+ax1.tick_params(axis='y', labelcolor='crimson')
+plt.yticks(fontsize=8)
+ax1.spines['top'].set_visible(False)
 
+ax2 = ax1.twinx()
+ax2.fill_between(time_axis, mu_eod+std_eod, mu_eod-std_eod, color='dodgerblue', alpha=0.5)
+ax2.plot(time_axis, mu_eod, color='black', lw=2)
+ax2.set_ylabel('voltage [mV]', fontsize=12)
+ax2.tick_params(axis='y', labelcolor='dodgerblue')
+
+plt.xticks(fontsize=8)
+plt.yticks(fontsize=8)
+fig.tight_layout()
 plt.show()
 
+#NixToFrame(data_dir)

code/base_chirps.py

@@ -0,0 +1,75 @@
+from read_chirp_data import *
+#import nix_helpers as nh 
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed
+
+
+
+data_dir = "../data"
+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", "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")
+
+
+
+#for dataset in data:
+eod = read_chirp_eod(os.path.join(data_dir, dataset))
+times = read_chirp_times(os.path.join(data_dir, dataset))
+
+
+
+df_map = {} #Keys werden nach df sortiert ausgegeben
+for k in eod.keys():
+	df = k[1]
+	ch = k[3] 
+	if df in df_map.keys():
+		df_map[df].append(k)
+	else:
+		df_map[df] = [k]
+
+print(ch) #die Chirphöhe wird ausgegeben, um zu bestimmen, ob Chirps oder Chirps large benutzt wurde
+
+
+
+#die äußere Schleife geht für alle Keys durch und somit durch alle dfs
+#die innnere Schleife bildet die 16 Wiederholungen einer Frequenz in 4 Subplots ab
+for idx in df_map.keys():
+	freq = list(df_map[idx]) 
+	fig,axs = plt.subplots(2, 2, sharex = True, sharey = True)
+
+	for idx, k in enumerate(freq):
+		ct = times[k]
+		e1 = eod[k]
+		zeit = e1[0]
+		eods = e1[1]
+	
+		if idx <= 3:
+			axs[0, 0].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+			axs[0, 0].scatter(np.asarray(ct), np.ones(len(ct))*3, color = 'green', s= 22)
+		elif 4<= idx <= 7:
+			axs[0, 1].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+			axs[0, 1].scatter(np.asarray(ct), np.ones(len(ct))*3, color = 'green', s= 22)
+		elif 8<= idx <= 11:
+			axs[1, 0].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+			axs[1, 0].scatter(np.asarray(ct), np.ones(len(ct))*3, color = 'green', s= 22)
+		else: 
+			axs[1, 1].plot(zeit, eods, color= 'blue', linewidth = 0.25)
+			axs[1, 1].scatter(np.asarray(ct), np.ones(len(ct))*3, color = 'green', s= 22)
+
+
+fig.suptitle('EOD for chirps', fontsize = 16)
+plt.show()
+
+
+
+#Problem: axs hat keine label-Funktion, also müsste axes nochmal definiert werden. Momentan erscheint Schrift nur auf einem der Subplots
+
+#ax = plt.gca()
+#ax.set_ylabel('Time [ms]')
+#ax.set_xlabel('Amplitude [mV]')
+#ax.label_outer()
+
+
+
+
+#next Step: relative Amplitudenmodulation berechnen, Max und Min der Amplitude bestimmen, EOD und Chirps zuordnen, Unterschied berechnen

code/base_eod.py

@@ -0,0 +1,21 @@
+from read_baseline_data import *
+#import nix_helpers as nh 
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed #Funktionen importieren
+
+
+data_dir = "../data"
+dataset = "2018-11-09-aa-invivo-1"
+#data = ("2018-11-09-aa-invivo-1", "2018-11-09-ab-invivo-1", "2018-11-09-ac-invivo-1", "2018-11-09-ad-invivo-1", "2018-11-13-aa-invivo-1", "2018-11-13-ab-invivo-1", "2018-11-13-ad-invivo-1", "2018-11-09-af-invivo-1", "2018-11-09-ag-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-ab-invivo-1", "2018-11-14-ac-invivo-1", "2018-11-14-ad-invivo-1", "2018-11-14-ae-invivo-1", "2018-11-14-af-invivo-1", "2018-11-14-ag-invivo-1", "2018-11-14-aa-invivo-1", "2018-11-14-aj-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")
+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')#Plottitelk
+plt.xlabel('Time [ms]', fontsize = 12)#Achsentitel
+plt.ylabel('Amplitude[mv]', fontsize = 12)#Achsentitel
+plt.xticks(fontsize = 12)
+plt.yticks(fontsize = 12)
+plt.show()

code/base_spikes.py

@@ -0,0 +1,34 @@
+from read_baseline_data import *
+#import nix_helpers as nh 
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed #Funktionen importieren
+
+data_dir = "../data"
+dataset = "2018-11-09-aa-invivo-1"
+#data = ("2018-11-09-aa-invivo-1", "2018-11-09-ab-invivo-1", "2018-11-09-ac-invivo-1", "2018-11-09-ad-invivo-1", "2018-11-13-aa-invivo-1", "2018-11-13-ab-invivo-1", "2018-11-13-ad-invivo-1", "2018-11-09-af-invivo-1", "2018-11-09-ag-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-ab-invivo-1", "2018-11-14-ac-invivo-1", "2018-11-14-ad-invivo-1", "2018-11-14-ae-invivo-1", "2018-11-14-af-invivo-1", "2018-11-14-ag-invivo-1", "2018-11-14-aa-invivo-1", "2018-11-14-aj-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")
+spike_times = read_baseline_spikes(os.path.join(data_dir, dataset))
+
+
+#spike_frequency = len(spike_times) / spike_times[-1]
+#inst_frequency = 1. / np.diff(spike_times)
+spike_rate = np.diff(spike_times)
+
+
+x = np.arange(0.001, 0.01, 0.0001)
+plt.hist(spike_rate,x)
+
+mu = np.mean(spike_rate)
+sigma = np.std(spike_rate)
+cv = sigma/mu
+print(cv)
+
+plt.title('A.lepto ISI Histogramm', fontsize = 14)
+plt.xlabel('duration ISI[ms]', fontsize = 12)
+plt.ylabel('number of ISI', fontsize = 12)
+
+plt.xticks(fontsize = 12)
+plt.yticks(fontsize = 12)
+plt.show()
+
+

code/code_cw.py

@@ -1,33 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-
-freq  = 800
-freq2 = 820
-dt = 0.00001
-x = np.arange(0.0, 1.0, dt)
-eod = np.sin(x * 2 * np.pi * freq) + np.sin(x * 2 * np.pi * freq * 2) * 0.5
-eod2 = np.sin(x * 2 * np.pi * freq2) + np.sin(x * 2 * np.pi * freq2 * 2) * 0.5
-  
-fig = plt.figure(figsize=(5., 7.5))
-ax= fig.add_subplot(311)
-ax.plot(x, eod, color="darkgreen", linewidth = 1.0)
-ax.set_xlim(0.0, 0.1)
-ax.set_ylabel("voltage [mV]")
-
-
-ax= fig.add_subplot(312)
-ax.plot(x, eod2, color="crimson", linewidth = 1.0)
-ax.set_xlim(0.0, 0.1)
-ax.set_ylabel("voltage [mV]")
-
-ax= fig.add_subplot(313)
-ax.plot(x, eod + eod2 * 0.05, color="lightblue", linewidth = 1.0)
-ax.set_xlim(0.0, 0.1)
-ax.set_xlabel("time [s]")
-ax.set_ylabel("voltage [mV]")
-
-plt.tight_layout()
-plt.savefig("eods.pdf")
-plt.show()
-
-

code/liste.py

@@ -0,0 +1,43 @@
+# 9.11.18
+
+aa: quality: poor, depth: -1341, base
+ab: quality: poor, depth: -1341, base
+ac: quality: good, depth: -1341, base
+ad: quality: good, depth: -200, base, chirps
+ae: quality: good, depth: -200, chirps
+af: quality: good, depth: -200
+ag: no info.dat, chirps
+
+
+# 13.11.18
+
+aa: good, -30 µm, maybe no reaction, base, chirps
+ab: good, -309 µm, base
+ac: poor, -309 µm, chirps
+ad: fair, -360 µm, base, chirps
+ae: fair, -350 µm
+af: good, -440 µm, bursting, base
+ag: fair, -174 µm, base
+ah: good, -209 µm, base, chirps, FI, SAM
+ai: good, -66.9 µm, base, chirps, SAM
+aj: good, -132 µm, base, chirps
+ak: good, -284 µm, base, chirps
+al: good, -286 µm, base, chirps, SAM
+
+
+# 14.11.18
+
+aa: good, -184 µm, base, chirps, FI, SAM, noise
+ab: fair, -279 µm, no reaction, base
+ac: fair, -60 µm, base, chirps
+ad: good, -357 µm, base, chirps
+ae: fair, -357 µm, base
+af: fair, -527 µm, base, (chirps)
+ag: fair, -533 µm, base, chirps
+ah: poor, -505 µm, chirps
+ai: good, -500 µm, still same cell 3x, chirps, FI, noise
+aj: poor, -314 µm, no modulation, base
+ak: good, -140 µm, base, chirps, FI, SAM, noise
+al: good, -280 µm, base, chirps, SAM
+am: good, -320 µm, base, chirps, FI, SAM, noise
+an: good, -434 µm, base, chirps, FI, SAM, noise

code/read_baseline_data.py

@@ -7,7 +7,13 @@ def read_baseline_eod(dataset):
     base = dataset.split(os.path.sep)[-1] + ".nix"
     nix_file = nix.File.open(os.path.join(dataset, base), nix.FileMode.ReadOnly)
     b = nix_file.blocks[0]
-    t = b.tags["BaselineActivity_1"]
+    if 'BaselineActivity_1' in b.tags:
+        t = b.tags["BaselineActivity_1"]
+    elif "BaselineActivity_2" in b.tags:
+        t = b.tags["BaselineActivity_2"]
+    else:
+        f.close()
+        return [],[]
     eod_da = b.data_arrays["LocalEOD-1"]
     eod = t.retrieve_data("LocalEOD-1")[:]
     time = np.asarray(eod_da.dimensions[0].axis(len(eod)))
@@ -19,7 +25,13 @@ def read_baseline_spikes(dataset):
     base = dataset.split(os.path.sep)[-1] + ".nix"
     nix_file = nix.File.open(os.path.join(dataset, base), nix.FileMode.ReadOnly)
     b = nix_file.blocks[0]
-    t = b.tags["BaselineActivity_1"]
+    if 'BaselineActivity_1' in b.tags:
+        t = b.tags["BaselineActivity_1"]
+    elif "BaselineActivity_2" in b.tags:
+        t = b.tags["BaselineActivity_2"]
+    else:
+        f.close()
+        return [],[]
     spikes_da = b.data_arrays["Spikes-1"]
     spike_times = spikes_da[:spikes_da.shape[0]-5000]
     baseline_spikes = spike_times[(spike_times > t.position[0]) & (spike_times < (t.position[0] + t.extent[0]))]

code/read_chirp_data.py

@@ -2,7 +2,7 @@ import numpy as np
 import os
 
 
-def load_chirp_spikes(dataset):
+def read_chirp_spikes(dataset):
     spikes_file = os.path.join(dataset, "chirpspikess1.dat")
     if not os.path.exists(spikes_file):
         print("found no chirps!")
@@ -15,11 +15,11 @@ def load_chirp_spikes(dataset):
         if "index" in l and "chirp" not in l:
             index = int(l.split(":")[-1])
         if "deltaf" in l and "true" not in l:
-           df = l.split(":")[-1]
+           df = l.split(":")[-1].strip()
         if "contrast" in l and "true" not in l:
-            contrast = l.split(":")[-1]
+            contrast = l.split(":")[-1].strip()
         if "chirpsize" in l:
-            cs = l.split(":")[-1]
+            cs = l.split(":")[-1].strip()
         if "#Key" in l:
             spikes[(index, df, contrast, cs)] = {}
         if "chirp index" in l:
@@ -32,7 +32,7 @@ def load_chirp_spikes(dataset):
     return spikes
 
 
-def load_chirp_eod(dataset):
+def read_chirp_eod(dataset):
     eod_file = os.path.join(dataset, "chirpeodampls.dat")
     if not os.path.exists(eod_file):
         print("found no chirpeodampls.dat!")
@@ -45,11 +45,11 @@ def load_chirp_eod(dataset):
         if "index" in l and "chirp" not in l:
             index = int(l.split(":")[-1])
         if "deltaf" in l and "true" not in l:
-           df = l.split(":")[-1]
+           df = l.split(":")[-1].strip()
         if "contrast" in l and "true" not in l:
-            contrast = l.split(":")[-1]
+            contrast = l.split(":")[-1].strip()
         if "chirpsize" in l:
-            cs = l.split(":")[-1]
+            cs = l.split(":")[-1].strip()
         if "#Key" in l:
             chirp_eod[(index, df, contrast, cs)] = ([], [])
         if len(l.strip()) != 0 and "#" not in l:
@@ -60,7 +60,7 @@ def load_chirp_eod(dataset):
     return chirp_eod
 
 
-def load_chirp_times(dataset):
+def read_chirp_times(dataset):
     chirp_times_file = os.path.join(dataset, "chirpss.dat")
     if not os.path.exists(chirp_times_file):
         print("found no chirpss.dat!")
@@ -73,15 +73,15 @@ def load_chirp_times(dataset):
         if "index" in l and "chirp" not in l:
             index = int(l.split(":")[-1])
         if "deltaf" in l and "true" not in l:
-            df = l.split(":")[-1]
+            df = l.split(":")[-1].strip()
         if "contrast" in l and "true" not in l:
-            contrast = l.split(":")[-1]
+            contrast = l.split(":")[-1].strip()
         if "chirpsize" in l:
-            cs = l.split(":")[-1]
+            cs = l.split(":")[-1].strip()
         if "#Key" in l:
             chirp_times[(index, df, contrast, cs)] = []
         if len(l.strip()) != 0 and "#" not in l:
-            chirp_times[(index, df, contrast, cs)].append(float(l.split()[1]))
+            chirp_times[(index, df, contrast, cs)].append(float(l.split()[1]) * 1000.)
     return chirp_times
 
 

code/spikes_analysis.py

@@ -0,0 +1,51 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from read_chirp_data import *
+from utility import *
+from IPython import embed
+
+data_dir = "../data"
+dataset = "2018-11-09-ad-invivo-1"
+
+spikes = read_chirp_spikes(os.path.join(data_dir, dataset))
+eod = read_chirp_eod(os.path.join(data_dir, dataset))
+times = read_chirp_times(os.path.join(data_dir, dataset))
+
+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]
+
+# make phases together, 12 phases
+spikes_mat = {}
+for deltaf in df_map.keys():
+    for rep in df_map[deltaf]:
+        for phase in spikes[rep]:
+            #print(phase)
+            spikes_one_chirp = spikes[rep][phase]
+            if deltaf == '-50Hz' and phase == (9, 0.54):
+                spikes_mat[deltaf, rep, phase] = spikes_one_chirp
+
+plot_spikes = spikes[(0, '-50Hz', '20%', '100Hz')][(0, 0.789)]
+
+mu = 1
+sigma = 1
+time_gauss = np.arange(-4, 4, 1)
+gauss = gaussian(time_gauss, mu, sigma)
+# spikes during time vec (00010000001)?
+smoothed_spikes = np.convolve(plot_spikes, gauss, 'same')
+window = np.mean(np.diff(plot_spikes))
+time_vec = np.arange(plot_spikes[0], plot_spikes[-1]+window, window)
+
+fig, ax = plt.subplots()
+ax.scatter(plot_spikes, np.ones(len(plot_spikes))*10, marker='|', color='k')
+ax.plot(time_vec, smoothed_spikes)
+plt.show()
+
+#embed()
+#exit()
+#hist_data = plt.hist(plot_spikes, bins=np.arange(-200, 400, 20))
+#ax.plot(hist_data[1][:-1], hist_data[0])

code/utility.py

@@ -1,8 +1,8 @@
 import numpy as np
 
 
-def zero_crossing(eod,time):
+def zero_crossing(eod, time):
-	threshold = 0;
+	threshold = 0
 	shift_eod = np.roll(eod, 1)
 	eod_times = time[(eod >= threshold) & (shift_eod < threshold)]
 	sampling_rate = 40000.0
@@ -10,9 +10,16 @@ def zero_crossing(eod,time):
 	return eod_idx
 
 
-def vector_strength(spike_times, eod_durations)
+def vector_strength(spike_times, eod_durations):
-	alphas = spike_times/ eod_durations
+	n = len(spike_times)
-	cs = (1/len(spike_times))*np.sum(np.cos(alphas))^2
+	phase_times = np.zeros(len(spike_times))
-	sn = (1/len(spike_times))*np.sum(np.sin(alphas))^2
+	for i, idx in enumerate(spike_times):
-	vs = np.sprt(cs+sn)
+		phase_times[i] = (spike_times[i] / eod_durations[i]) * 2 * np.pi
+	vs = np.sqrt((1/n*sum(np.cos(phase_times)))**2 + (1/n*sum(np.sin(phase_times)))**2)
 	return vs
+
+
+def gaussian(x, mu, sig):
+	y = np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
+	return y
+

code/vector_phase.py

@@ -0,0 +1,25 @@
+from read_baseline_data import *
+from utility import *
+#import nix_helpers as nh 
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed #Funktionen importieren
+
+
+#Zeitpunkte einer EOD über Zero-crossings finden, die in einer Steigung liegen
+data_dir = "../data"
+dataset = "2018-11-09-ad-invivo-1"
+time,eod = read_baseline_eod(os.path.join(data_dir, dataset))
+spike_times = read_baseline_spikes(os.path.join(data_dir, dataset))
+print(len(spike_times))
+
+eod_times = zero_crossing(eod,time)
+eod_durations = np.diff(eod_times)
+print(len(spike_times))
+print(len(eod_durations))
+
+#for st in spike_times:
+	#et = eod_times[eod_times < st]
+	#dt = st - et
+	
+#vs = vector_strength(spike_times, eod_durations)