From a7217fb4c651b4796a018505c616594b9ce05a8e Mon Sep 17 00:00:00 2001 From: Ramona Date: Thu, 29 Nov 2018 17:06:17 +0100 Subject: [PATCH] wuhu --- code/plot_eodform_spikehist.py | 4 +- code/response_beat.py | 61 +++++++++++++++----- code/spikes_analysis.py | 56 +++++++++--------- code/spikes_beat.py | 62 ++++++++++++++++++++ code/spikes_chirp.py | 102 +++++++++++++++++++++++++++++++++ 5 files changed, 242 insertions(+), 43 deletions(-) create mode 100644 code/spikes_beat.py create mode 100644 code/spikes_chirp.py diff --git a/code/plot_eodform_spikehist.py b/code/plot_eodform_spikehist.py index 8452399..9512ca6 100644 --- a/code/plot_eodform_spikehist.py +++ b/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-13-aa-invivo-1' +#dataset = '2018-11-09-ad-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)) diff --git a/code/response_beat.py b/code/response_beat.py index 2d28c12..e899e6b 100644 --- a/code/response_beat.py +++ b/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", - #"2018-11-13-ak-invivo-1", "2018-11-13-al-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"] -''' # 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] - norm_train = train*1000/spikerate - rep_rates.append(np.std(norm_train))#/spikerate) + #train = smoothed_data[window:beat_window+window] + #norm_train = train*1000/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())): - ax.plot(np.ones(len(rates[k]))*k, rates[k], 'o') +for i, cell in enumerate(rates.keys()): + 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() +''' \ No newline at end of file diff --git a/code/spikes_analysis.py b/code/spikes_analysis.py index ea8648e..53750d8 100644 --- a/code/spikes_analysis.py +++ b/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_rates[df] = [] - beat[df] = [] + #csi_trains[df] = {} + csi_rates[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())): - ax.plot(np.ones(len(csi_trains[k]))*i, csi_trains[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())) +for i, k in enumerate(sorted(beat.keys())): + ax.plot(np.ones(len(beat[k]))*k, beat[k], 'o') fig.tight_layout() plt.show() -''' -''' fig, ax = plt.subplots() -for i, k in enumerate(sorted(beat.keys())): - ax.plot(np.ones(len(beat[k]))*i, beat[k], 'o') -ax.legend(sorted(beat.keys()), loc='upper left', bbox_to_anchor=(1.04, 1)) -#ax.set_xticklabels(sorted(csi_trains.keys())) +for i, k in enumerate(sorted(csi_trains.keys())): + ax.plot(np.ones(len(csi_trains[k]))*i, csi_trains[k], 'o') +ax.plot(np.arange(-1, len(csi_trains.keys())+1), np.zeros(len(csi_trains.keys())+2), 'silver', linewidth=2, linestyle='--') fig.tight_layout() plt.show() ''' \ No newline at end of file diff --git a/code/spikes_beat.py b/code/spikes_beat.py new file mode 100644 index 0000000..c9f124c --- /dev/null +++ b/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() diff --git a/code/spikes_chirp.py b/code/spikes_chirp.py new file mode 100644 index 0000000..fcdda70 --- /dev/null +++ b/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()