diff --git a/code/NixFrame.py b/code/NixFrame.py new file mode 100644 index 0000000..5fd5743 --- /dev/null +++ b/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 .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) + diff --git a/code/analysis_rs.py b/code/analysis_rs.py index f7973f1..950b572 100644 --- a/code/analysis_rs.py +++ b/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)) -ax = fig.add_subplot(111) +fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor)) 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.yticks(fontsize = 8) +plt.xticks(fontsize=8) +plt.yticks(fontsize=8) fig.tight_layout() -plt.savefig('eod.pdf') - -#interspikeintervalhistogram, windowsize = 1 ms -#plt.hist -#coefficient of variation -#embed() -#exit() +#plt.savefig('eod.pdf') +plt.show() +# 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 -# plot mean of eod circles -# plot std of eod circles -# plot psth into the same plot -# calculate vector strength - -threshold = 0; +# calculate eod times and indices by zero crossings +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() -for i, idx in enumerate(eod_idx): - #embed() - #exit() - plt.plot(time[int(idx):int(eod_idx[i+1])], eod[int(idx):int(eod_idx[i+1])]) +# align eods and spikes to eods +max_cut = int(np.max(np.diff(eod_idx))) +eod_cuts = np.zeros([len(eod_idx)-1, max_cut]) +spike_times = [] +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) \ No newline at end of file diff --git a/code/base_chirps.py b/code/base_chirps.py new file mode 100644 index 0000000..dfdfe37 --- /dev/null +++ b/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 diff --git a/code/base_eod.py b/code/base_eod.py new file mode 100644 index 0000000..888084a --- /dev/null +++ b/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() diff --git a/code/base_spikes.py b/code/base_spikes.py new file mode 100644 index 0000000..7ab065f --- /dev/null +++ b/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() + + diff --git a/code/code_cw.py b/code/code_cw.py deleted file mode 100644 index 23f2af2..0000000 --- a/code/code_cw.py +++ /dev/null @@ -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() - - diff --git a/code/eod_sim_rs.py b/code/eod_sim_rs.py deleted file mode 100644 index e69de29..0000000 diff --git a/code/liste.py b/code/liste.py new file mode 100644 index 0000000..7a51eda --- /dev/null +++ b/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 diff --git a/code/read_baseline_data.py b/code/read_baseline_data.py index 6b3cd29..9e951a9 100644 --- a/code/read_baseline_data.py +++ b/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]))] diff --git a/code/read_chirp_data.py b/code/read_chirp_data.py index 4283c48..c3afe76 100644 --- a/code/read_chirp_data.py +++ b/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 diff --git a/code/spikes_analysis.py b/code/spikes_analysis.py new file mode 100644 index 0000000..5ae3da4 --- /dev/null +++ b/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]) \ No newline at end of file diff --git a/code/utility.py b/code/utility.py index 86f2407..68102d3 100644 --- a/code/utility.py +++ b/code/utility.py @@ -1,8 +1,8 @@ import numpy as np -def zero_crossing(eod,time): - threshold = 0; +def zero_crossing(eod, time): + 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) - alphas = spike_times/ eod_durations - cs = (1/len(spike_times))*np.sum(np.cos(alphas))^2 - sn = (1/len(spike_times))*np.sum(np.sin(alphas))^2 - vs = np.sprt(cs+sn) +def vector_strength(spike_times, eod_durations): + n = len(spike_times) + phase_times = np.zeros(len(spike_times)) + for i, idx in enumerate(spike_times): + 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 + diff --git a/code/vector_phase.py b/code/vector_phase.py new file mode 100644 index 0000000..c902e6e --- /dev/null +++ b/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)