import numpy as np import matplotlib.pyplot as plt from IPython import embed from scipy import signal import time from scipy.stats import alpha from scipy.stats import norm from random import random import os import pandas as pd def load_cell(set, fname = 'singlecellexample5', big_file = 'data_beat', new = False): if (not os.path.exists(fname +'.pkl')) or (new == True): data_all = pd.read_pickle(big_file +'.pkl') d = data_all[data_all['dataset'] == set] results = pd.DataFrame(d) results.to_pickle(fname +'.pkl') #np.save('singlecellexample5.npy', results) else: d = pd.read_pickle(fname +'.pkl') return d def default_settings(data,lw = 1,intermediate_width = 2.6*3,intermediate_length = 3, ts = 10, ls = 10, fs = 10): inch_factor = 2.54 plt.rcParams['figure.figsize'] = (intermediate_width, intermediate_length) plt.rcParams['font.size'] = fs plt.rcParams['axes.titlesize'] = ts plt.rcParams['axes.labelsize'] = ls plt.rcParams['lines.linewidth'] = lw plt.rcParams['lines.markersize'] = 8 plt.rcParams['legend.loc'] = 'upper right' plt.rcParams["legend.frameon"] = False def remove_tick_marks(ax): labels = [item.get_text() for item in ax.get_xticklabels()] empty_string_labels = [''] * len(labels) ax.set_xticklabels(empty_string_labels) return ax def remove_tick_ymarks(ax): labels = [item.get_text() for item in ax.get_yticklabels()] empty_string_labels = [''] * len(labels) ax.set_yticklabels(empty_string_labels) return ax def euc_d(conv_beat,conv_chirp,a): d_euclidean = np.sqrt(np.nansum((conv_beat - conv_chirp) ** 2, a)) return d_euclidean def mean_euc_d(conv_beat,conv_chirp,a): mean = np.sqrt(np.nanmean((conv_beat - conv_chirp) ** 2,a)) return mean def ten_per_to_std(width): deviation = width / 4.29 return deviation def gauss(mult,sampling,width,x,conv_fact = 100): # width in s deviation = ten_per_to_std(width) deviation = int(sampling * deviation) gaussian_n = just_gauss(conv_fact,mult,deviation,x) return gaussian_n def just_gauss(conv_fact,mult, deviation,x, func = 'alpha',test = False): points,to_cut = points_cut(deviation, mult,conv_fact) if func == 'gaussian': gaussian = signal.gaussian(points, std=deviation, sym=True) elif func == 'exp': gaussian = signal.exponential(points, 0, deviation, False) if test == True: check_alpha() check_exponentials() elif func == 'exp_self': x_axis = np.arange(0,points, 1) gaussian = np.exp(-x_axis / deviation) elif func == 'alpha': x_axis = np.arange(0, points, 1) gaussian = x_axis* np.exp(-2.45 * x_axis / deviation) gaussian_n = (gaussian*x) / np.sum(gaussian) return gaussian_n def points_cut(deviation, mult,conv_fact): if deviation * 10 % 2: points = deviation * mult else: points = deviation * mult - 1 to_cut = (int(points / 2))/conv_fact return points,to_cut def check_exponentials(): devs = [1, 5, 25, 50, 75, 100, 150, 250, 350, 500, 700] for i in range(len(devs)): if deviation * 10 % 2: points = devs[i] * mult else: points = devs[i] * mult - 1 gaussian = signal.exponential(points, 0, devs[i], False) plt.plot(gaussian) plt.show() def check_alpha(): conv_fact = 100 aa = [2, 3, 4, 5, 6, 7, 10] for i in range(len(aa)): a = aa[i] # x = np.arange(alpha.ppf(0.01, a), # alpha.ppf(0.99, a), 1/100000) x = np.linspace(alpha.ppf(0.01, a), alpha.ppf(0.99, a), 1000) # alp = alpha.pdf(x, a)/(np.sum(alpha.pdf(x, a))) # print(np.argmax(alp)/100) alp = alpha.pdf(x, a) print(alpha.ppf(0.90, a)) # plt.scatter(alpha.ppf(0.90, a),alpha.pdf(0.90, a)) plt.scatter(alpha.pdf(0.90, a), alpha.ppf(0.90, a)) # print(len(alp)) plt.plot(np.arange(0, 1000, 1) / conv_fact, alp, label=str(a)) plt.xlabel('ms') #alpha_func = 1/(x**2*norm.cdf(a)*np.sqrt(2*np.pi))*np.exp(-0.5*(a-1/x)*2) plt.legend() plt.show() rv = alpha(a) x = np.linspace(0,10,points) plt.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf') gaussian = signal.gaussian(points, std=deviation, sym=True) def shift_gauss(step_ms,sampling,C,gaussian): G_all = [] step = int(step_ms*sampling/1000) for i in np.arange(0, len(C) - len(gaussian)+1, step): # shift by 3 ms = 300 data points G = np.zeros(len(C)) G[0 + i:len(gaussian) + i] = gaussian G_all.append(G) return G_all, step_ms def shift_chirp(step_ms,sampling,C,gaussian): step = int(step_ms * sampling / 1000) a = np.arange(0, len(C) - len(gaussian)+1, step) G_all = [[]]*len(a) for i in a: # shift by 3 ms = 300 data points G_all[int(i/step)]= C[0 + i:len(gaussian) + i] #G_all.append(chirp) G_all = np.asarray(G_all) return G_all, step_ms def shift_beat(step_ms,sampling,C,gaussian): G_all = [] step = int(step_ms*sampling/1000) for i in np.arange(0, len(C[0]) - len(gaussian)+1, step): # shift by 3 ms = 300 data points chirp= C[:,0 + i:len(gaussian) + i] G_all.append(chirp) G_all = np.asarray(G_all) return G_all, step_ms def consp_create(type,G_all, C, B): nom = np.dot(G_all, np.transpose(C * B)) if type == 'corr_p': denom = np.transpose((np.transpose((np.sqrt(np.dot(G_all, np.transpose(np.array(B) * np.array(B))))))* np.transpose(((np.sqrt(np.nansum(G_all * C * C, axis=1))))))) if type == 'corr_a': denom = np.transpose((np.transpose( (0.5 * np.dot(G_all, np.transpose(np.array(B) * np.array(B))))) + np.transpose( ((0.5 * np.nansum(G_all * C * C, axis=1)))))) index = nom/denom #if denom.all() == 0 and nom.all() == 0: # S_shift2 = np.ones(len(G_all)) #else: S_shift2 = np.nanmax(index, axis=1) ConSpi2 = 1 - np.abs(np.asarray(S_shift2)) all = 1-index return ConSpi2,all def consp_max(chir, step_ms, ConSpi2, left_c): if chir == 'stim': side = 17 / step_ms I = ConSpi2[int(0.5 * len(ConSpi2)-side):int(0.5 * len(ConSpi2)+side)]#FIXME check if the middle is really the middle if chir == 'chirp': I = ConSpi2 Integral = np.sum(I) max_position = np.argmax(ConSpi2) -left_c max_value = np.nanmax(I) return Integral, max_position, max_value def remove_mean(chirp, beat): C = chirp - np.mean(chirp) B = np.transpose(np.transpose(beat) - np.mean(beat, axis = 1)) return C, B def single_dist(left_c, beat_window, eod_fe, eod_fr, col, filter, beat, chirp, sampling, chir, plot = False, show = False,stages = 'two',test = False): C, B = remove_mean(chirp, beat) step_ms = 1 if stages == 'two': C_all, step_ms = shift_chirp(step_ms, sampling, C, filter) B_all, step_ms = shift_beat(step_ms, sampling, B, filter) a, b, c = np.shape(B_all) distances = [[]] * a for i in range(a): distances[i] = mean_euc_d(B_all[i], C_all[i], 1) optimal_phase = np.nanmin(distances, axis=1) max_value = np.nanmax(optimal_phase) else: C_all = C B_all = B a, b = np.shape(B_all) distances = [[]] * a for i in range(a): distances[i] = mean_euc_d(B_all[i], C_all, None) optimal_phase = np.nanmin(distances) optimal_phase_position = np.nanargmin(distances) max_value = optimal_phase position = optimal_phase_position if test == True: for i in range(len(B_all)): plt.title(str(i)+' '+str(distances[i])) #plt.plot(B_all[0], color = 'red') plt.plot(B_all[int(len(B_all)/2)], color='red') plt.plot(B_all[i], color='blue') plt.plot(C_all,color = 'green') plt.show() if plot: plot_cons(filter, left_c, beat_window, optimal_phase, eod_fe, col, sampling, B, C, 1, chir) if show: plt.show() return max_value,distances,position def single_consp(type,left_c, beat_window, eod_fe, eod_fr, col, filter, beat, chirp, sampling, chir, plot = False, show = False): C, B = remove_mean(chirp, beat) step_ms = 1 G_all, step_ms = shift_gauss(step_ms, sampling, C, filter) ConSpi2,all = consp_create(type,G_all, C, B) Integral, max_position, max_value = consp_max(chir, step_ms, ConSpi2, left_c) if plot: plot_cons(filter, left_c, beat_window, ConSpi2, eod_fe, col, sampling, B, C, 0, chir) if show: plt.show() return max_position, max_value, Integral,ConSpi2,all def plot_cons(gaussian, left_c,beat_window,ConSpi2,beat,col, sampling, B, C, side, var): fig = plt.figure(figsize=[4, 3]) ax = create_subpl(0.4, fig, 3, 3, 1, share2=True) if var == 'chirp': ax[0].plot((np.arange(0, len(ConSpi2), 1) - 0.5 * len(ConSpi2))+len(B[0])/(100*2) -left_c, ConSpi2, color=col) ax[0].axvline(x = len(B[0])/(100*2)-3.5) ax[0].plot(np.linspace(-left_c+(len(gaussian)/100)/2, len(B[0]) / 100 - left_c-(len(gaussian)/100)/2, len(ConSpi2)) , ConSpi2, color='red') if var == 'stim': ax[0].plot((np.arange(0, len(ConSpi2), 1) - 0.5 * len(ConSpi2)), ConSpi2, color=col) ax[0].axvline(x=beat_window, color = 'red') ax[0].axvline(x=-beat_window, color = 'red') ax[1].axvline(x=-beat_window-13.5, color = 'red') ax[1].axvline(x=+beat_window+13.5, color = 'red') ax[1].axvline(x=0, color='red') ax[0].set_title(beat) ax[0].set_ylabel('Conspicousy [0-1]') if var == 'stim': ax[0].axvline(x=int(side), color=col) ax[0].axvline(x=-int(side), color=col) if var == 'chirp': ax[0].axvline(x=int(side/2 +len(B[0]) / (100 * 2)-left_c), color=col) ax[0].axvline(x=int(-side/2 +len(B[0]) / (100 * 2)-left_c), color=col) ax[0].set_ylim([0, 1]) #if chirp == 'chirp': # B = np.transpose(B) if var == 'chirp': ax[2].plot(np.linspace(-left_c,len(B[0])/100 -left_c,len(B[0])), np.transpose(B[0:int(len(B)):1])) if var == 'stim': ax[2].plot(np.arange(-len(np.transpose(B[0:int(len(B)):1])) * 0.5 * 1000 / sampling, len(np.transpose(B[0:int(len(B)):1])) * 0.5 * 1000 / sampling, 1000 / sampling), np.transpose(B[0:int(len(B)):1])) ax[2].set_title('10 shifted Beats') ax[2].set_ylabel('AM') if var == 'chirp': ax[1].plot(np.linspace(-left_c,len(B[0])/100 -left_c,len(B[0])), C, color=col) if var == 'stim': ax[1].plot(np.arange(-len(np.transpose(B[0:int(len(B)):1])) * 0.5 * 1000 / sampling, len(np.transpose(B[0:int(len(B)):1])) * 0.5 * 1000 / sampling, 1000 / sampling), C, color=col) ax[1].set_ylabel('AM') ax[1].set_title('Chirp time window') def mean_euc_d3(conv_beat,conv_chirp,a): mean = np.sqrt((np.nansum((conv_beat - conv_chirp)/len(conv_beat) ** 2,a))) return mean def mean_euc_d2(conv_beat,conv_chirp,a): mean = np.sqrt(np.nanmean((conv_beat - conv_chirp) ** 2,a))/len(conv_beat) return mean def frame_subplots(w, h, xlabel, ylabel,title): fig = plt.figure(figsize=[w, h]) plt.axes(frameon=False) plt.xticks([]) plt.yticks([]) plt.xlabel(xlabel,labelpad=30).set_visible(True) plt.ylabel(ylabel, labelpad=50).set_visible(True) ttl = plt.title(title) ttl.set_position([.5, 1.05]) return fig def create_subpl(h,fig, subplots,nrow, ncol, share2 = False): ax = {} if share2: ax[0] = fig.add_subplot(nrow, ncol, 1) for i in range(subplots): if share2: ax[i] = fig.add_subplot(nrow,ncol, i+1,sharex=ax[0]) else: ax[i] = fig.add_subplot(nrow,ncol, i+1) plt.subplots_adjust(hspace=h) return ax def windows_periods2(bef_c_t, aft_c_t,beat_window, min_bw, win, period, time_transform, period_shift, period_distance, perc_bef, perc_aft): length = ((min_bw / time_transform) / (1 + perc_bef))* time_transform #embed() if length < bef_c_t+aft_c_t: print('length in win2 is smaller as the chirp') elif length > (bef_c_t+aft_c_t)*2.5: length = (bef_c_t+aft_c_t)*2.5 for p in range(len(period)): if win == 'w2': #embed() length_exp = 0.020 period_shift[p] = np.ceil(((length_exp+0.012)* time_transform)/period[p])*period[p] #0.12 is the length of the deviations i guess #period_shift[p] = np.ceil((length_exp + 0.04 * time_transform) / period[p]) * period[ # p] # 0.12 is the length of the deviations i guess #embed() #period_shift[p] = 2*np.ceil((length_exp * time_transform) / period[p]) * period[p] period_distance[p] = length_exp * time_transform if (period_shift[p]+period_distance[p] *perc_bef)>beat_window[p]: print('period'+str(p)+' uncool 1') period_shift[p] = np.ceil((length_exp * time_transform) / period[p]) * period[p] if (period_shift[p]+period_distance[p] *perc_bef) > beat_window[p]: print('period' + str(p) + ' uncool 2') period_shift[p] = 0.020 * time_transform if win == 'w4': period_shift[p] = length period_distance[p] = length period_left = period_distance *perc_bef period_right = period_distance * perc_aft #embed() return period_shift,period_left, period_right def auto_rows(data): ncol = int(np.sqrt(len(data))) nrow = int(np.sqrt(len(data))) if ncol * nrow < len(data): ncol = int(np.sqrt(len(data))) nrow = int(np.sqrt(len(data)) + 1) if ncol * nrow < len(data): ncol = int(np.sqrt(len(data)) + 1) nrow = int(np.sqrt(len(data)) + 1) return nrow, ncol def period_calc(beat_window, m, win, deviation_s, time_transform, beat_corr, bef_c, aft_c, chir, ver = '',conv_fact = 100,shuffled = ''): period, bef_c_t, aft_c_t, deviation_t, period_shift, period_distance, perc_bef, perc_aft,nr_of_cycles,interval = pre_window(beat_corr, time_transform, deviation_s, bef_c, aft_c, chir) #if win == 'w1'or win == 'w3': period_shift, period_left, period_right,exclude,consp_needed = windows_periods(ver,m, time_transform, win, beat_window,perc_bef, perc_aft, period_shift, period_distance, period, interval, nr_of_cycles, bef_c_t, aft_c_t,shuffled = shuffled) #if win == 'w2' or win == 'w4': # period_shift, period_left, period_right = windows_periods2(bef_c_t, aft_c_t,beat_window, m, win, period, time_transform, period_shift, period_distance, perc_bef, perc_aft) #if win == 'w3': # period_shift,period_left, period_right = windows_periods3(beat_window,perc_bef, perc_aft, period_shift, period, interval, nr_of_cycles, bef_c_t, aft_c_t) period_distance_c, period_distance_b, left_b, right_c, left_c, right_b, to_cut = post_window(conv_fact,period_left, period_right, deviation_t, period_shift) #embed() #left_c = left_c+right_c+np.abs(right_b) #right_c = right_c+right_c+np.abs(right_b) #left_b_consp = left_b-consp_needed dels = {} dels['left_c'] = left_c dels['right_c'] = right_c dels['left_b'] = left_b dels['right_b'] = right_b dels['consp_needed'] = consp_needed dels['to_cut'] = to_cut return left_c, right_c, left_b, right_b,period_distance_c, period_distance_b,period_shift,period,to_cut,exclude,consp_needed,dels,interval def choose_randomly(beat_window,interval,l_all = 400,r_all = 800): #np.array(beat_window) * 2 - np.abs(bef_c_t) - np.abs(aft_c_t) #left_c = np.zeros(len(beat_window)) #for i in range(len(beat_window)): ra = random() # print(ra) # left_c[i] = ra * (np.array(beat_window[i]) * 2 - interval) - beat_window[i] left_c = ra * (r_all + np.abs(l_all) - interval) - np.abs(l_all) right_c = left_c + interval return left_c, right_c def windows_periods3(beat_window,perc_bef, perc_aft, period_shift, period, interval, nr_of_cycles, bef_c_t, aft_c_t): period_left = np.zeros(len(period)) period_right = np.zeros(len(period)) for p in range(len(period)): # & period[p] * perc_bef+ if (period[p] < interval) or ((period[p]+ period[p] * perc_bef)>beat_window[p]): period_shift[p] = (period[p] * nr_of_cycles[p]) period_left[p] = bef_c_t period_right[p] = aft_c_t else: period_left[p] = period[p] * perc_bef period_right[p] = period[p] * perc_aft period_shift[p] = period[p] if (period[p]+ bef_c_t)>beat_window[p]: period_shift[p] = (bef_c_t+aft_c_t) period_left[p] = bef_c_t period_right[p] = aft_c_t return period_shift,period_left, period_right def windows_periods(ver, min_bw, time_transform, win, beat_window,perc_bef, perc_aft,period_shift, period_distance, period, interval, nr_of_cycles, bef_c_t, aft_c_t, make_interval_shorter = False,closest_to_2m = True,shuffled = ''): exclude = np.zeros(len(period)) exclude = [[]]*len(period) consp_needed = [[]] * len(period) interval = abs(aft_c_t) + abs(bef_c_t) nr_of_cycles = ((interval) / period).astype(int) + 1 # nr of cycles in the length of the chirp puffer = abs(bef_c_t) -4 if win == 'w1' or win == 'w2' or win == 'w3': for p in range(len(period)): if 'consp' in ver: consp_needed[p] = period[p] * 0.5 consp_puffered = consp_needed[p] - puffer if consp_puffered <0: consp_puffered = 0 else: consp_needed[p] = 0 consp_puffered = 0 mult_interval = period[p] * nr_of_cycles[p] shorter_bw_aft = interval * perc_aft - 13.5 if period[p] < interval: if period[p] == 0: period_shift[p] = (interval) * 2 period_distance[p] = (interval) exclude[p] = True else: if win == 'w1': period_distance[p] = mult_interval elif win == 'w3': period_distance[p] = interval elif win == 'w2': period_distance[p] = interval if closest_to_2m == True: if ((mult_interval*2) <= beat_window[p]+shorter_bw_aft-period_distance[p]* perc_bef -consp_needed[p]) and (mult_interval * 2 >= period_distance[p] +consp_puffered): period_shift[p] = mult_interval * 2 exclude[p] = True else: period_shift[p] = (((beat_window[p]+shorter_bw_aft- period_distance[p] * perc_bef-consp_needed[p]) / period[p]).astype(int)) * period[p] exclude[p] = True if period_shift[p] < mult_interval+consp_puffered: if make_interval_shorter == False: if ver == 'consp': #embed() exclude[p] = False #period_shift[p] = [] elif ver == 'dist': exclude[p] = False #period_shift[p] = [] elif ver == 'conspdist': consp_needed[p] = 0 if ((mult_interval * 2) <= beat_window[p]+shorter_bw_aft - period_distance[p] * perc_bef - consp_needed[p]) and (mult_interval * 2 >= period_distance[p] + consp_needed): period_shift[p] = mult_interval * 2 exclude[p] = True else: period_shift[p] = (((beat_window[p]+shorter_bw_aft - period_distance[p] * perc_bef - consp_needed[p]) / period[p]).astype(int)) * period[p] if period_shift[p] < mult_interval + consp_needed[p]: exclude[p] = False #period_shift[p] = [] if make_interval_shorter == True: short_interval = 5 #bef_c_t, aft_c_t, perc_bef, perc_aft = percentage(bef_c, aft_c, time_transform) # FIXME a period distance adapt for w2 (and also for the backwards calculation # FIXME in beat window adapt 0.8 times period plus minus # FIXME look that it doesnt has to be 2 periods in long periods shorter_bw_aft = (interval)*perc_aft-13.5 period_distance[p], period_shift[p] = exclude_interval(consp_needed[p],shorter_bw_aft, p,interval,short_interval,period, win, period_distance, period_shift, beat_window, perc_bef) exclude[p] = True if period_shift[p] <= mult_interval+consp_needed[p]: short_interval = 7.5 #bef_c_t, aft_c_t, perc_bef, perc_aft = percentage(bef_c, aft_c, time_transform) shorter_bw_aft = (interval) * perc_aft - 13.5 period_distance[p], period_shift[p] = exclude_interval(consp_needed[p],shorter_bw_aft, p,interval, short_interval, period, win, period_distance, period_shift, beat_window, perc_bef) exclude[p] = True if period_shift[p] < mult_interval+consp_needed[p]: exclude[p] = False #period_shift[p] = [] print('exclude') else: exclude[p] = True elif closest_to_2m == False: if (mult_interval * 2) <= beat_window[p] - period_distance[p] * perc_bef - consp_needed[p]: period_shift[p] = mult_interval * 2 exclude[p] = True if ((period[p] * nr_of_cycles[p] + np.round((12.5/period[p]))*period[p]+ period[p] * nr_of_cycles[p] * perc_bef)>beat_window[p]): period_shift[p] = period[p] * nr_of_cycles[p] else: period_shift[p] = period[p] * nr_of_cycles[p] + np.round((12.5/period[p]))*period[p] #elif ((period[p] + period[p] * perc_bef) > beat_window[p]+shorter_bw_aft): # #embed() # if (win == 'w3') or (win == 'w1'): # exclude[p] = False # print('should be excluded, not at least a period!!') # else: # period_distance[p] = interval # if ((period[p] + period_distance[p] * perc_bef) > beat_window[p] + shorter_bw_aft): # exclude[p] = False # print('should be excluded, not at least a period!!') # else: # exclude[p] = True # if ((bef_c_t + aft_c_t) * 2 + bef_c_t) > beat_window[p]: # period_shift[p] = (interval) # period_distance[p] = (interval) # else: # period_shift[p] = (interval) * 2 # period_distance[p] = (interval) else: if win == 'w1': period_distance[p] = period[p] elif win == 'w3': period_distance[p] = period[p] elif win == 'w2': period_distance[p] = interval if ((period[p]*2 + period_distance[p] * perc_bef) <= beat_window[p]-consp_needed[p]+shorter_bw_aft) and (period[p]*2 >= period_distance[p] + consp_puffered): exclude[p] = True period_shift[p] = period[p]*2 elif((period[p] + period_distance[p] * perc_bef) <= beat_window[p]-consp_needed[p]+shorter_bw_aft) and (period[p] >= period_distance[p] + consp_puffered): exclude[p] = True period_shift[p] = period[p] else: if ver == 'conspdist': consp_needed[p] = 0 if ((period[p] * 2 + period_distance[p] * perc_bef) <= beat_window[ p] - consp_needed[p] + shorter_bw_aft) and ( period[p] * 2 >= period_distance[p] +consp_puffered): exclude[p] = True period_shift[p] = period[p] * 2 elif ((period[p] + period_distance[p] * perc_bef) <= beat_window[ p] - consp_needed[p] + shorter_bw_aft) and (period[p] >= period_distance[p] +consp_puffered): exclude[p] = True period_shift[p] = period[p] else: exclude[p] = False #period_shift[p] = [] else: exclude[p] = False #period_shift[p] = [] #elif ((period[p] + period_distance[p] * perc_bef) < beat_window[p]+shorter_bw_aft): # if ver == 'consp': # exclude[p] = True # else: # exclude[p] = False # period_shift[p] = period[p] #if period_shift[p]beat_window[p]: # print('error') #if win == 'w2': # period_distance[p] = interval #embed() #break #if win == 'w2' or win == 'w4': # for p in range(len(period)): # if win == 'w2': # #embed() # length_exp = 0.020 # period_shift[p] = np.ceil(((length_exp+0.012)* time_transform)/period[p])*period[p] #0.12 is the length of the deviations i guess # #period_shift[p] = np.ceil((length_exp + 0.04 * time_transform) / period[p]) * period[ # # p] # 0.12 is the length of the deviations i guess # #embed() # #period_shift[p] = 2*np.ceil((length_exp * time_transform) / period[p]) * period[p] # period_distance[p] = length_exp * time_transform # if (period_shift[p]+period_distance[p] *perc_bef)>beat_window[p]: # print('period'+str(p)+' uncool 1') # period_shift[p] = np.ceil((length_exp * time_transform) / period[p]) * period[p] # if (period_shift[p]+period_distance[p] *perc_bef) > beat_window[p]: # print('period' + str(p) + ' uncool 2') # period_shift[p] = 0.020 * time_transform elif win == 'w4': for p in range(len(period)): if interval*perc_bef+interval*2<=beat_window[p]: period_shift[p] = interval*2 period_distance[p] = interval exclude[p] = True # FIXME here das mit den consp noch reinmachen elif interval*perc_bef+interval<=beat_window[p]: period_shift[p] = interval period_distance[p] = interval exclude[p] = True else: exclude[p] = False #period_shift[p] = [] #length = ((min_bw / time_transform) / (1 + perc_bef))* time_transform #embed() #exclude[p] = False #if length < bef_c_t+aft_c_t: # exclude[p] = True # print('length in win2 is smaller as the chirp') #elif length > (bef_c_t+aft_c_t)*2.5: # exclude[p] = False # length = (bef_c_t+aft_c_t)*2.5 #else: # exclude[p] = False #period_shift[p] = length*2 #period_distance[p] = length period_left = period_distance *perc_bef period_right = period_distance * perc_aft #embed() return period_shift, period_left, period_right,exclude,consp_needed def exclude_interval(consp_needed,bw, p, interval, short_interval, period, win, period_distance, period_shift, beat_window, perc_bef): nr_of_cycles = ((interval - short_interval) / period).astype(int) + 1 interval = interval - short_interval if win == 'w1': period_distance[p] = period[p] * nr_of_cycles[p] elif win == 'w3': period_distance[p] = interval elif win == 'w2': period_distance[p] = interval period_shift[p] = (((beat_window[p]+bw - period[p] * nr_of_cycles[p] * perc_bef -consp_needed) / period[ p]).astype(int)) * period[p] return period_distance[p], period_shift[p] def post_window(conv_fact,period_left, period_right, deviation_t, period_shift): points, to_cut = points_cut(deviation_t*conv_fact,10,conv_fact) #to_cut = 5 * deviation_t #embed() left_c = np.round((-period_left- to_cut)*conv_fact)/conv_fact right_c = np.round((period_right+ to_cut)*conv_fact)/conv_fact left_b = np.round((-period_left - to_cut -period_shift)*conv_fact)/conv_fact right_b = np.round((period_right+ to_cut -period_shift)*conv_fact)/conv_fact period_distance_c = abs(right_c - left_c)-2*to_cut period_distance_b = abs(right_b - left_b)-2*to_cut return period_distance_c, period_distance_b, left_b, right_c, left_c, right_b,to_cut def one_dict_str(beats, spike_phases): AUCI = {} for b in range(len(beats)): AUCI[beats[b]] = [[]] * len(spike_phases[beats[b]]) return AUCI def one_array_str(beats, spike_phases): AUCI = [[]]*len(beats) for b in range(len(beats)): AUCI[b] = [[]] * len(spike_phases[beats[b]]) return AUCI def only_one_dict_str(beats, spike_phases): AUCI = {} for b in range(len(beats)): AUCI[beats[b]] = [[]] * len(spike_phases) return AUCI def three_dict_str(deviation_list, data): AUCI = {} #embed() for dev in range(len(deviation_list)): AUCI[deviation_list[dev]] = {} for d in range(len(data)): AUCI[deviation_list[dev]][data[d]] = {} return AUCI def four_dict_str(deviation_list, data, beats, bin_array): AUCI = {} for dev in range(len(deviation_list)): AUCI[deviation_list[dev]] = {} for d in range(len(data)): AUCI[deviation_list[dev]][data[d]] = [[]]*len(beats[d]) for b in range(len(beats[d])): AUCI[deviation_list[dev]][data[d]][b] = [[]]*len(bin_array) for p in range(len(bin_array)): AUCI[deviation_list[dev]][data[d]][b][p] = [] return AUCI def pre_window(beat_corr, time_transform,deviation_s,bef_c,aft_c,chir): period = period_func(beat_corr, time_transform, chir) deviation_t = deviation_s * time_transform period_shift = np.zeros(len(period)) period_distance = np.zeros(len(period)) exclude = np.zeros(len(period)) bef_c_t, aft_c_t, perc_bef, perc_aft = percentage(bef_c, aft_c, time_transform) nr_of_cycles = ((bef_c_t + aft_c_t) / period).astype(int) + 1 # nr of cycles in the length of the chirp interval = abs(aft_c_t) + abs(bef_c_t) #embed() return period, bef_c_t, aft_c_t, deviation_t, period_shift, period_distance, perc_bef, perc_aft,nr_of_cycles,interval def percentage(bef_c, aft_c, time_transform): bef_c_t = bef_c * time_transform aft_c_t = aft_c * time_transform perc_bef = (np.abs(bef_c_t) / (np.abs(bef_c_t) + np.abs(aft_c_t))) perc_aft = (np.abs(aft_c_t) / (np.abs(bef_c_t) + np.abs(aft_c_t))) return bef_c_t, aft_c_t, perc_bef, perc_aft def period_func(beat_corr,time_transform, chir): if chir == 'chir': period = np.abs(1 / beat_corr) * time_transform if chir == 'stim': period = (np.round(np.abs(1 / beat_corr) * time_transform)).astype(int) period[period == float('inf')] = np.max(period[period != np.inf]) return period def two_ylabels(y_pad,x,y,z,fig,left, right): ax1 = fig.add_subplot(x, y, z) ax1.set_xticks([]) ax1.set_yticks([]) ax0 = ax1.twinx() ax0.set_xticks([]) ax0.set_yticks([]) ax1.set_ylabel(left, labelpad = y_pad) ax0.set_ylabel(right, rotation = 270,labelpad = 40) def del_ticks(data, ncol, ax): for d in range(len(data) - ncol): ax[d].set_xticks([]) a = np.arange(0, len(data), 1) deleted = np.delete(a, (np.arange(0, len(data), ncol))) for d in deleted: ax[d].set_yticks([]) def axvline_everywhere(ax,subplots,eod_fe,eod_fr): for e in range(subplots): for i in range(int(np.max(eod_fe) / np.max(eod_fr))): ax[e].axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')