import nixio as nix import os from IPython import embed #from utility import * import matplotlib.pyplot as plt import numpy as np import pandas as pd import matplotlib.mlab as ml import scipy.integrate as si from scipy.ndimage import gaussian_filter from IPython import embed from myfunctions import * from myfunctions import auto_rows from functionssimulation import default_settings import matplotlib.gridspec as gridspec from myfunctions import remove_tick_marks import string from myfunctions import load_cell def ps_df(data, d = '2019-09-23-ad-invivo-1', wish_df = 310, window = 'no',sampling_rate = 40000): #nfft = 4096 #trial_cut = 0.1 #freq_step = sampling_rate / nfft #d = load_cell(data[0], fname='singlecellexample5_2', big_file='beat_results_smoothed') data_cell = data[data['dataset'] == d]# dfs = np.unique(data_cell['df']) df_here = dfs[np.argmin(np.abs(dfs - wish_df))] dfs310 = data_cell[data_cell['df'] == df_here] #pp = [[]]*len(dfs310) pp = [] ppp = [] trial_cut = 0.1 for i in range(len(dfs310)): duration = dfs310.iloc[i]['durations'] #cut_vec = np.arange(0, duration, trial_cut) cut_vec = np.arange(0, duration, trial_cut) #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)] #for j, cut in enumerate(cut_vec): # # print(j) # spike_times = dfs310.iloc[i]['spike_times'] # spikes = spike_times - spike_times[0] # spikes_cut = spikes[(spikes > cut) & (spikes < cut_vec[j + 1])] # if cut == cut_vec[-2]: # #counter_cut += 1 # break # if len(spikes_cut) < 10: # #counter_spikes += 1 # break # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1) # spikes_idx = np.round((spikes_cut - trial_cut * j) * sampling_rate) # for spike in spikes_idx: # spikes_mat[int(spike)] = 1# # # #spikes_mat = np.zeros(int(spikes[-1]* sampling_rate + 5)) # #spikes_idx = np.round((spikes) * sampling_rate) # #for spike in spikes_idx: # # spikes_mat[int(spike)] = 1 # spikes_mat = spikes_mat * sampling_rate # if type(window) != str: # spikes_mat = gaussian_filter(spikes_mat, sigma=window) # # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate # # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate # else: # smoothened = spikes_mat * 1 # nfft = 4096 # p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2) # pp.append(p) spike_times = dfs310.iloc[i]['spike_times'] if len(spike_times) < 3: counter_spikes += 1 break spikes = spike_times - spike_times[0] spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)] if len(spikes_cut) < 3: counter_cut += 1 break spikes_mat = np.zeros(int(spikes[-1] * sampling_rate + 5)) spikes_idx = np.round((spikes) * sampling_rate) for spike in spikes_idx: spikes_mat[int(spike)] = 1 spikes_mat = spikes_mat * sampling_rate if type(window) != str: spikes_mat = gaussian_filter(spikes_mat, sigma=window) # smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05) * sampling_rate # smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2) * sampling_rate else: spikes_mat = spikes_mat*1 nfft = 4096 p, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2) ppp.append(p) #spike_times = data_cell.iloc[i]['spike_times']# #if len(spike_times) < 3: # counter_spikes += 1 # break #spikes = spike_times - spike_times[0] # cut trial into snippets of 100 ms #cut_vec = np.arange(0, duration, trial_cut) #spikes_cut = spikes[(spikes > 0.05) & (spikes < 0.95)] #if len(spikes_cut) < 3: # counter_cut += 1 # break #spikes_new = spikes_cut - spikes_cut[0] #spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2) # spikes_mat = np.zeros(int(trial_cut * sampling_rate) + 1) #spikes_idx = np.round((spikes_new) * sampling_rate) #for spike in spikes_idx: # spikes_mat[int(spike)] = 1 #spikes_mat = spikes_mat * sampling_rate #nfft = 4096 #p, f = ml.psd(smoothened - np.mean(smoothened), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2) #ppp.append(p) #p_mean = np.mean(pp,axis = 0) p_mean2 = np.mean(ppp, axis=0) #ref = (np.max(p_mean2)) # db = 10 * np.log10(p_mean2 / np.max(p_mean2)) #ref = (np.max(p_mean2)) #db2 = 10 * np.log10(p_mean2 / ref) #embed() return df_here,p_mean2,f,db def plot_example_ps(grid,lw_hline = 0.8,line_style = (0,(1,10)), hs = 0.5,hr = [1,1], multiples = 1.1,colors = ['brown'],fc = 'lightgrey',line_col = 'black',input = ['2019-10-21-aa-invivo-1'],sigma = [0.00005,0.00025,0.0005, 0.002],wish_df = 150, color_eod = 'orange',color_stim = 'red', color_df = 'green'): sampling_rate = 40000 #colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'] plt.rcParams['lines.linewidth'] = 1.5 plt.rcParams['lines.markersize'] = 6 #data = pd.read_pickle('../pictures_highbeats/data_beat.pkl') #iter = np.unique(data['dataset']) iter = ['2019-05-07-by-invivo-1'] iter = ['2019-09-23-ad-invivo-1'] iter = input for cell in iter: data = pd.read_pickle('data_beat.pkl') beat_results = pd.read_pickle('beat_results_smoothed.pkl') #embed() eodf = int(beat_results[beat_results['dataset'] == cell]['eodf'].iloc[0]) df = [[]] * (len(sigma) + 1) p = [[]] * (len(sigma) + 1) f = [[]] * (len(sigma) + 1) db = [[]] * (len(sigma) + 1) sigmaf = [[]] * (len(sigma) + 1) gauss = [[]] * (len(sigma) + 1) df[0], p[0], f[0], db[0] = ps_df(data, d=cell, wish_df= wish_df, window='no', sampling_rate=sampling_rate) for i in range(len(sigma)): df[1+i], p[1+i], f[1+i], db[1+i] = ps_df(data, d=cell, wish_df= wish_df, window = sigma[i]*sampling_rate,sampling_rate = sampling_rate) sigmaf[i + 1] = 1 / (2 * np.pi * sigma[i]) gauss[i + 1] = np.exp(-(f[1+i] ** 2 / (2 * sigmaf[i + 1] ** 2))) db = 'no' stepsize = f[0][1] - f[0][0] if db == 'db': p = db # fig.suptitle(d, labelpad = 25) #print(d) ax = {} ec = 'grey' #fc = 'moccasin' #ec = 'wheat' scale = 1 #ax = plot_whole_ps(f, ax, grid, colors, eodf, stepsize, p, df, scale = scale, ax_nr = 0,nr=0, filter='whole' ,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec) #ax[0].legend( loc=(0,1), # ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1), #ax[0] = remove_tick_marks(ax[0]) plots = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=grid[0], height_ratios = hr,wspace=0.4, hspace=hs) ax = plot_whole_ps(f,ax, plots, colors, eodf, stepsize, p, df,mult = multiples, scale = scale, ax_nr = 0,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec) ax[0].legend( loc=(0,1), ncol=3, mode="expand", borderaxespad=0.)#bbox_to_anchor=(0.4, 1, 0.6, .1), ax[0] = remove_tick_marks(ax[0]) nr_size = 10 ax[0].text(-0.1, 1.1, string.ascii_uppercase[0], transform=ax[0].transAxes, size=nr_size, weight='bold') #ax[0].set_ylim([0, 2000]) wide = 2 #embed() for i in range(len(sigma)): plots = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=grid[i+1],height_ratios = hr, wspace=0.4, hspace=hs) ax[i+2] = plt.subplot( plots[0]) plot_filter(ax, i+2, f[1+i], p,i+1, colors, gauss[1+i], eodf, stepsize, wide, df[1+i],scale = scale,color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec) ax[i+2].set_ylim([0, eodf*multiples]) ax[2].text(-0.1, 1.1, string.ascii_uppercase[1], transform=ax[2].transAxes, size=nr_size, weight='bold') ax[3].text(-0.1, 1.1, string.ascii_uppercase[2], transform=ax[3].transAxes, size=nr_size, weight='bold') ax[2] = remove_tick_marks(ax[2]) #embed() #if db == 'db': # ax[0].set_ylim([np.min([p]),0])#p[0][,p[1][0:2000],p[2][0:2000],p[3][0:2000] #else: # ax[0].set_ylim([ 0,np.max([p])]) ax[int(len(df))].set_ylabel('Frequency [Hz]') # ax[1].set_ylabel(r'power [Hz$^2$/Hz]') #ax[0].ticklabel_format(axis='y', style='sci', scilimits=[0, 0]) #print(df[3]) for i in [0,2,3]: ax[i].spines['right'].set_visible(False) ax[i].spines['top'].set_visible(False) cols = grid.ncols rows = grid.nrows ax[int(len(df))].set_xlabel(' psd [Hz²/Hz]') #ax[2].set_ylabel('Hz²/Hz') #ax[3].set_ylabel('Hz²/Hz') #ax[0].set_ylabel('Hz²/Hz') for i in [0,2,3]: ax[i].axhline(y = eodf/2, color = line_col, linestyle = line_style, linewidth = lw_hline) plt.tight_layout() #embed() #fig.label_axes() def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, mult = 1.05,ax_nr = 0,nr = 0, filter = 'original', scale = 1, color_eod = 'orange',color_stim = 'red', color_df = 'green',fc = 'lightgrey', ec = 'grey',): ax[ax_nr] = plt.subplot(grid[ax_nr]) if filter == 'whole': #ax[nr].set_facecolor('lightgrey') ax[ax_nr].plot(p[nr], f[nr], color=colors[0]) ax[ax_nr].fill_between([np.min(p), np.max(p)], [f[0][-1],f[0][-1]], color=fc,edgecolor=ec) ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0], color=color_df, marker='o', linestyle='None', label='Df') ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o', linestyle='None', label='stimulus') ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod, marker='o', linestyle='None', label='EODf') # = '+str(int(eodf))+' Hz') elif filter == 'original': #ax[nr].fill_between([eodf] * len(p[nr]), p[nr], color='lightgrey') #ax[nr].fill_between([max(p[0])]*len(f[nr]),f[nr], color = 'lightgrey') #embed() ax[ax_nr].plot(p[nr], f[nr], color=colors[0]) #ax[ax_nr].plot(p[nr][f[nr] eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0]) #embed() ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0], color=color_df, marker='o', linestyle='None', label='Df') ax[ax_nr].plot(p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, color=color_stim, marker='o', linestyle='None', label='stimulus') ax[ax_nr].plot(np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, eodf - 1, color=color_eod, marker='o', linestyle='None', label='EODf') # = '+str(int(eodf))+' Hz') ax[ax_nr].fill_between([np.min(p),np.max(p)], [eodf/2,eodf/2], color=fc,edgecolor=ec) #ax[ax_nr].plot(np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale, df[0], # color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black' #ax[ax_nr].plot(0, df[nr] + eodf, color=color_stim, marker='o', # linestyle='None', # label='stimulus',zorder = 2)#,edgecolors = 'black' #ax[ax_nr].plot(0, eodf - 1, color=color_eod, # marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz') #plt.plot([np.min(p),np.max(p)],[eodf,eodf], color = 'red') #embed() #ax[nr].plot([0]*5) #ax[nr].plot([1000]*5) # ax[0].fill_between( [max(p[0])]*len(f[1]),f[0], facecolor='lightgrey', edgecolor='grey') ax[ax_nr].set_ylim([0, eodf * mult]) ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1]) return ax def plot_filter(ax, ax_nr, f, p4,array_nr, colors, gauss3, eodf, stepsize, wide, df, fc = 'lightgrey', scale = 1, ec = 'grey',color_eod = 'orange',color_stim = 'red', color_df = 'green'): ax[ax_nr].plot( p4[array_nr],f, color=colors[0]) prev_height = np.max((p4[0][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale) now_height = np.max((p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale) ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black') ax[ax_nr].scatter( now_height+440, np.abs(df), marker = '>', color='black', zorder = 2) #embed() ax[ax_nr].fill_between(max(p4[0]) * gauss3 ** 2,f, facecolor=fc, edgecolor=ec) ax[ax_nr].plot(np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale, eodf, color=color_eod, marker='o', linestyle='None') ax[ax_nr].plot( np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,abs(df), color=color_df, marker='o', linestyle='None') ax[ax_nr].plot( np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,df + eodf, color=color_stim, marker='o', linestyle='None') #p[nr][int((df[nr] + eodf) / stepsize) + 1], df[nr] + eodf, #np.max(p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale, df + eodf #embed() ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1]) return ax def plot_amp(ax, mean1, dev,name = 'amp',nr = 1): np.unique(mean1['type']) all_means = mean1[mean1['type'] == name +' mean'] original = all_means[all_means['dev'] == 'original'] #m005 = all_means[all_means['dev'] == '005'] m05 = all_means[all_means['dev'] == '05'] m2 = all_means[all_means['dev'] == '2'] # fig, ax = plt.subplots(nrows=4, ncols = 3, sharex=True) versions = [original, m05, m2] #m005, for i in range(len(versions)): keys = [k for k in versions[i]][2::] try: data = np.array(versions[i][keys])[0] except: break axis = np.arange(0, len(data), 1) axis_new = axis * 1 similarity = [keys, data] sim = np.argsort(similarity[0]) # similarity[sim] all_means = mean1[mean1['type'] == name+' std'] std = all_means[all_means['dev'] == dev[i]] std = np.array(std[keys])[0] #ax[1, 1].set_ylabel('Modulation depth') #ax[nr,i].set_title(dev[i] + ' ms') all_means = mean1[mean1['type'] == name+' 95'] std95 = all_means[all_means['dev'] == dev[i]] std95 = np.array(std95[keys])[0] all_means = mean1[mean1['type'] == name+' 05'] std05 = all_means[all_means['dev'] == dev[i]] std05 = np.array(std05[keys])[0] ax[nr,i].fill_between(np.array(keys)[sim], list(std95[sim]), list(std05[sim]), color='gainsboro') ax[nr,i].fill_between(np.array(keys)[sim], list(data[sim] + std[sim]), list(data[sim] - std[sim]), color='darkgrey') # ax[i].plot(data_tob.ff, data_tob.fe, color='grey', linestyle='--', label='AMf') ax[nr,i].plot(np.array(keys)[sim], data[sim], color='black') # ax[0].plot(data1.x, data1.freq20, color=colors[1], label='20 %') #embed() return ax def create_beat_corr(hz_range, eod_fr): beat_corr = hz_range%eod_fr beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2] return beat_corr def plot_mean_cells( grid,lw_hline = 0.8, hs = 0.5,line_style = (0,(1,10)), multiples = 1.1,data = ['2019-10-21-aa-invivo-1'],hr = [1,1],line_col = 'black',lw = 0.5, sigma = ['original','05','2'],colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], wish_df = 150, color_eod = 'black',color_df = 'orange', size = 17, color_modul = ['steelblue']): #mean1 = pd.read_pickle('mean.pkl') data_all = pd.read_pickle('beat_results_smoothed.pkl') d = data_all[data_all['dataset'] == data[0]] #embed() inch_factor = 2.54 plt.rcParams['lines.markersize'] = 8 plt.rcParams['legend.loc'] = 'upper right' plt.rcParams["legend.frameon"] = False # load data for plot # data1 = pd.read_csv('ma_allcells_unsmoothed.csv') # data2 = pd.read_csv('ma_allcells_05.csv') # data3 = pd.read_csv('ma_allcells_2.csv') # data_tob = pd.read_csv('ma_toblerone.csv') # smothed = df_beat[df_beat['dev'] == 'original'] # data1 = smothed[smothed['type'] == 'amp'] x = np.arange(0, 2550, 50) corr = create_beat_corr(x, np.array([500] * len(x))) #np.unique(mean1['type']) #all_means = mean1[mean1['type'] == 'max mean'] #versions = [[]]*len(dev) #for i in range(len(dev)): version =[[]]*len(sigma) version2 = [[]] * len(sigma) dev = [[]] * len(sigma) limits = [[]]*len(sigma) minimum = [[]] * len(sigma) y_max = [[]] * len(sigma) y_min = [[]] * len(sigma) ax ={} for i, e in enumerate(sigma): y2 = d['result_amplitude_max_' + e] y_max[i] = np.max(y2) y_min[i] = np.min(y2) for i,e in enumerate(sigma): dev[i] = sigma[i] plots = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=grid[i], height_ratios = hr,wspace=0.4, hspace=hs) d = data_all[data_all['dataset'] == data[0]] x = d['delta_f'] / d['eodf'] + 1 #embed() data = ['2019-10-21-aa-invivo-1'] #end = ['original', '005', '05', '2'] y = d['result_frequency_' + e] #embed() y2 = d['result_amplitude_max_' + e] #y_sum[i] = np.nanmax(y) ff = d['delta_f'] / d['eodf'] + 1 fe = d['beat_corr'] #fig.suptitle(set) ax[0] = plt.subplot(plots[0]) #if e != sigma[-1]: ax[0] = remove_tick_marks(ax[0]) tob_col = 'black' tob_lw = 0.8 ax[0].plot(ff, fe, color=tob_col, linestyle='--', linewidth=tob_lw, zorder=1) ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw) #embed() eod = d['eodf'].iloc[0] ax[0].axhline(y=eod / 2, color=line_col, linestyle=line_style,linewidth = lw_hline) if np.max(y)