added py files I that are not primary plot versions here

This commit is contained in:
saschuta 2020-11-17 10:49:04 +01:00
parent 343945ed06
commit d758ff8722
10 changed files with 4043 additions and 0 deletions

209
differentcells_filter.py Normal file
View File

@ -0,0 +1,209 @@
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 axes import label_axes, labelaxes_params
from myfunctions import auto_rows
from myfunctions import default_settings
from myfunctions import remove_tick_marks
from myfunctions import remove_tick_ymarks
import matplotlib.gridspec as gridspec
from rotated import ps_df
def plot_single_cells(ax, data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']):
colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
# labelaxes_params(xoffs=-3, yoffs=0, labels='A', font=dict(fontweight='bold'))
#baseline = pd.read_pickle('data_baseline.pkl')
#data_beat = pd.read_pickle('data_beat.pkl')
data_all = pd.read_pickle('beat_results_smoothed.pkl')
#data = np.unique(data_all['dataset'])[0]
#data = np.unique(data_all['dataset'])
end = ['original', '005','05', '2' ]
end = ['original','005']
y_sum = [[]]*len(data)*len(end)
counter = 0
for dd,set in enumerate(data):
for ee, e in enumerate(end):
d = data_all[data_all['dataset'] == set]
#x = d['delta_f'] / d['eodf'] + 1
#embed()
#y = d['result_frequency_' + e]
y = d['result_amplitude_max_' + e]
#y2 = d['result_amplitude_max_' + e]
y_sum[counter] = np.nanmax(y)
counter += 1
#print(np.nanmax(y))
#embed()
lim = np.max(y_sum)
hd = 0.3
ws = 0.35
rows = 1
cols = 3
#embed()
main_grid = gridspec.GridSpec(1, 2,hspace=0.4, width_ratios = [1,7])
filter_grid = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=main_grid[0], hspace=0.4)
upper_filter_g = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=filter_grid[0], wspace=ws, hspace=0.4) # ,
grid1 = gridspec.GridSpecFromSubplotSpec(rows,cols,subplot_spec=upper_filter_g[0], wspace=ws, hspace=0.4)#,
wish_df = 150
fc = 'lightgrey'
ec = 'grey'
sampling_rate = 40000
data_beat = pd.read_pickle('data_beat.pkl')
df, p, f, db = ps_df(data_beat, d=set, wish_df=wish_df, window='no', sampling_rate=sampling_rate)
ax = {}
ax_nr = 0
colors = ['brown']
#embed()
ax[ax_nr] = plt.subplot(grid1[0])
ax[ax_nr].spines['right'].set_visible(False)
ax[ax_nr].spines['left'].set_visible(False)
ax[ax_nr].spines['top'].set_visible(False)
ax[ax_nr].spines['bottom'].set_visible(False)
ax[ax_nr].set_yticks([])
ax[ax_nr].set_xticks([])
ax[ax_nr].plot(p, f, color=colors[0])
eodf = d['eodf'].iloc[0]
mult = 1.1
ax[ax_nr].fill_between([np.min(p), np.max(p)], [eodf / 2, eodf / 2], color=fc, edgecolor=ec)
ax[ax_nr].set_ylim([0, eodf *0.6])
ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
upper_filter_g = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=filter_grid[1], wspace=ws, hspace=0.4) # ,
grid1 = gridspec.GridSpecFromSubplotSpec(rows,cols,subplot_spec=upper_filter_g[0], wspace=ws, hspace=0.4)#,
i = 0
sigma = 0.0005
df, p, f, db = ps_df(data_beat, d=set, wish_df=wish_df, window=sigma * sampling_rate,
sampling_rate=sampling_rate)
ax_nr = 1
sigmaf = 1 / (2 * np.pi * sigma)
gauss = np.exp(-(f ** 2 / (2 * sigmaf ** 2)))
stepsize = np.abs(f[0]-f[1])
wide = 2
scale = 1
prev_height = np.max((p[int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale)
now_height = np.max((p[int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) *scale)
ax[ax_nr] = plt.subplot(grid1[1])
ax[ax_nr].spines['right'].set_visible(False)
ax[ax_nr].spines['left'].set_visible(False)
ax[ax_nr].spines['top'].set_visible(False)
ax[ax_nr].set_yticks([])
ax[ax_nr].set_xticks([])
ax[ax_nr].spines['bottom'].set_visible(False)
ax[ax_nr].fill_between(max(p) * gauss ** 2, f, facecolor=fc, edgecolor=ec)
ax[ax_nr].plot([prev_height, now_height+440],[np.abs(df), np.abs(df)], color = 'black')
ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
grid = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=main_grid[1], hspace=0.4)
#grid = gridspec.GridSpecFromSubplotSpec(2,1,subplot_spec=grid[1], hspace = 0.4)
grid1 = gridspec.GridSpecFromSubplotSpec(rows,cols,subplot_spec=grid[0], wspace=ws, hspace=0.4)#,
end = ['original']
color_modul = 'steelblue'
color_mpf = 'red'
y_sum1 = plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,title = True,xlabel = False, label =False,remove =True)
grid2 = gridspec.GridSpecFromSubplotSpec(rows, cols, subplot_spec=grid[1], wspace=ws, hspace=0.4)#,
end = ['05']
#y_sum = [[]] * len(data)
color_modul = 'steelblue'
color_mpf = 'red'
y_sum2 = plot_single(lim, data, end, data_all, grid2, color_mpf, color_modul,title = False, label = True,remove =False)
#for dd, d in enumerate(data):
# embed()
#embed()
#ax[1].set_ylim([0, np.nanmax([y_sum1,y_sum2])])
#ax[1, dd].set_ylim([0, 350])
plt.subplots_adjust(wspace = 0.4,left = 0.17, right = 0.96,bottom = 0.2)
def plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,label = True, xlabel = True,remove =True, title = True):
y_sum = [[]] * len(data)
ax = {}
for dd,set in enumerate(data):
for ee, e in enumerate(end):
if title == True:
plt.title('Cell '+str(dd))
d = data_all[data_all['dataset'] == set]
x = d['delta_f'] / d['eodf'] + 1
#embed()
y = d['result_frequency_' + e]
y2 = d['result_amplitude_max_' + e]
y_sum[dd] = np.nanmax(y)
ff = d['delta_f'] / d['eodf'] + 1
fe = d['beat_corr']
#fig.suptitle(set)
grid2 = gridspec.GridSpecFromSubplotSpec(2,1, subplot_spec=grid1[dd], wspace=0.02, hspace=0.2) # ,
ax[0] = plt.subplot(grid2[0])
ax[0].plot(ff, fe, color='grey', linestyle='--')
ax[0].plot(x, y, color=color_mpf)
ax[set+e+str(1)] = plt.subplot(grid2[1])
if (set == data[0]) and (label == True):
ax[set+e+str(1)].set_ylabel('Modulation ')
ax[0].set_ylabel('MPF [EODf]')
if (set == data[0]) and (xlabel == True):
ax[set + e + str(1)].set_xlabel('EOD multiples')
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
ax[set+e+str(1)].plot(x, y2, color=color_modul)
# ax[1,0].set_ylabel('modulation depth [Hz]')
#ax[2, ee].plot(x, y2, color=colors[0])
#ax[2, 0].set_ylabel(' modulation depth [Hz]')
# ax[1,ee].annotate("", xy=(0.53, 16.83), xytext=(0.53, 17.33), arrowprops=dict(arrowstyle="->"))
# ax[1,ee].annotate("", xy=(1.51, 16.83), xytext=(1.51, 17.33), arrowprops=dict(arrowstyle="->"))
#ax[1, 0].set_xlabel('stimulus frequency [EODf]')
#ax[1, 2].set_xlabel('stimulus frequency [EODf]')
#ax[2, 3].set_xlabel('stimulus frequency [EODf]')
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[set+e+str(1)].spines['right'].set_visible(False)
ax[set+e+str(1)].spines['top'].set_visible(False)
#ax[2, ee].spines['right'].set_visible(False)
#ax[2, ee].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
#plt.tight_layout()
# fig.label_axes()
ax[set+e+str(1)].set_ylim([0, lim])
#ax[0].set_ylim([0, 240])
#ax[set+e+str(1)] = remove_tick_marks(ax[0])
ax[0] = remove_tick_marks(ax[0])
if set != data[0]:
ax[0] = remove_tick_ymarks(ax[0])
ax[set + e + str(1)] = remove_tick_ymarks(ax[set+e+str(1)] )
if remove == True:
ax[set+e+str(1)] = remove_tick_marks(ax[set+e+str(1)])
#ax[0].set_ylim([0, lim])
return y_sum
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']
data = ['2019-10-21-aa-invivo-1', '2019-10-21-au-invivo-1', '2019-10-28-aj-invivo-1']
default_settings(data,intermediate_width = 6.7,intermediate_length = 5)
#fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True)
ax = {}
plot_single_cells(ax, data = data)
#fig.savefig()
plt.savefig('differentcells_filter.pdf')
plt.savefig('../highbeats_pdf/differentcells_filter.pdf')
# plt.subplots_adjust(left = 0.25)
plt.show()
#plt.close()

191
differentcells_trans.py Normal file
View File

@ -0,0 +1,191 @@
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 axes import label_axes, labelaxes_params
from myfunctions import auto_rows
from myfunctions import default_settings
from myfunctions import remove_tick_marks
import matplotlib.gridspec as gridspec
import string
def plot_single_cells(ax, data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']):
colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B']
# labelaxes_params(xoffs=-3, yoffs=0, labels='A', font=dict(fontweight='bold'))
#baseline = pd.read_pickle('data_baseline.pkl')
#data_beat = pd.read_pickle('data_beat.pkl')
data_all = pd.read_pickle('beat_results_smoothed.pkl')
#data = np.unique(data_all['dataset'])[0]
#data = np.unique(data_all['dataset'])
end = ['original', '005','05', '2' ]
end = ['original','005']
y_sum = [[]]*len(data)*len(end)
counter = 0
for dd,set in enumerate(data):
for ee, e in enumerate(end):
d = data_all[data_all['dataset'] == set]
#x = d['delta_f'] / d['eodf'] + 1
#embed()
#y = d['result_frequency_' + e]
y = d['result_amplitude_max_' + e]
#y2 = d['result_amplitude_max_' + e]
y_sum[counter] = np.nanmax(y)
counter += 1
#print(np.nanmax(y))
#embed()
lim = np.max(y_sum)
#embed()
grid = gridspec.GridSpec(1,3,width_ratios = [0.2,4,4],)
grid0 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=grid[0], wspace=0.02, hspace=0.1) # ,
label1 = plt.subplot(grid0[0])
label2 = plt.subplot(grid0[1])
label3 = plt.subplot(grid0[2])
labels = [label1,label2,label3]
titels = ['Cell 1','Cell 2','Cell 3']# < 0.5 EODf,with 0.5 ms wide
fs_big = 11
weight = 'bold'
for ll,l in enumerate(labels):
#embed()
l.spines['right'].set_visible(False)
l.spines['left'].set_visible(False)
l.spines['top'].set_visible(False)
l.spines['bottom'].set_visible(False)
l.set_yticks([])
l.set_xticks([])
l.set_ylabel(titels[ll],labelpad = 15, fontsize = fs_big, fontweight=weight, rotation = 0)
hd = 0.3
grid1 = gridspec.GridSpecFromSubplotSpec(3,1,subplot_spec=grid[1], wspace=0.02, hspace=0.4)#,
end = ['original']
color_modul = 'steelblue'
color_mpf = 'red'
y_sum1,moduls = plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,arrows = False,nrs = [0,1,2], title = 'Binary spike trains',xlabel = True,yticks = False)
grid2 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=grid[2], wspace=0.02, hspace=0.4)#,
end = ['05']
end = ['2']
#y_sum = [[]] * len(data)
color_modul = 'steelblue'
color_mpf = 'red'
y_sum2, moduls = plot_single(lim, data, end, data_all, grid2, color_mpf, color_modul,arrows = True, mods = moduls, nrs = [3,4,5], title = '2 ms Gaussian',xlabel = False)
#for dd, d in enumerate(data):
# embed()
#embed()
#ax[1].set_ylim([0, np.nanmax([y_sum1,y_sum2])])
#ax[1, dd].set_ylim([0, 350])
plt.subplots_adjust(wspace = 0.4,left = 0.1, right = 0.96,bottom = 0.1)
def plot_single(lim, data, end, data_all, grid1, color_mpf, color_modul,mods = [], arrows = True, nr_size = 12, nrs = [0,1,2], xlabel = True,title = '',yticks =True):
y_sum = [[]] * len(data)
ax = {}
moduls = [[]]*len(data)
for dd,set in enumerate(data):
for ee, e in enumerate(end):
d = data_all[data_all['dataset'] == set]
x = d['delta_f'] / d['eodf'] + 1
#embed()
y = d['result_frequency_' + e]
y2 = d['result_amplitude_max_' + e]
y_sum[dd] = np.nanmax(y)
ff = d['delta_f'] / d['eodf'] + 1
fe = d['beat_corr']
#fig.suptitle(set)
grid2 = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=grid1[dd], wspace=0.02, hspace=0.2) # ,
ax[0] = plt.subplot(grid2[0])
ax[0].plot(ff, fe, color='grey', linestyle='--')
ax[0].plot(x, y, color=color_mpf)
ax[0].text(-0.1, 1.1, string.ascii_uppercase[nrs[dd]], transform=ax[0].transAxes,
size=nr_size, weight='bold')
if dd == 0:
plt.title(title)
ax[set+e+str(1)] = plt.subplot(grid2[1])
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
#if (e == 2) and xlabel == True:
ax[set+e+str(1)].plot(x, y2, color=color_modul,zorder = 2)
if arrows == True:
plt.fill_between(x, y2,mods[dd], color = 'gainsboro', edgecolor= 'grey',zorder = 1)
array = [0.65, 1.65, 2.65]
small_arrows = False
if small_arrows == True:
for a in range(len(array)):
#embed()
pos = np.argmin(np.abs(np.array(x)-array[a]))
x_present = np.array(x)[pos]
y2 = np.array(y2)
#embed()
pos_change = 2
nr = 25
#embed()
if (np.array(mods[dd])[pos]-y2[pos])>nr:
plt.plot([x_present, x_present],
[y2[pos] + nr, np.max(np.array(mods[dd])[pos - pos_change:pos + pos_change])],
color='black')
plt.scatter([x_present],[y2[pos]+nr], marker = 'v',s = 10, color='black')
moduls[dd] = y2
# ax[1,0].set_ylabel('modulation depth [Hz]')
#ax[2, ee].plot(x, y2, color=colors[0])
#ax[2, 0].set_ylabel(' modulation depth [Hz]')
# ax[1,ee].annotate("", xy=(0.53, 16.83), xytext=(0.53, 17.33), arrowprops=dict(arrowstyle="->"))
# ax[1,ee].annotate("", xy=(1.51, 16.83), xytext=(1.51, 17.33), arrowprops=dict(arrowstyle="->"))
#ax[1, 0].set_xlabel('stimulus frequency [EODf]')
if (dd == 2) and xlabel == True:
ax[set+e+str(1)].set_xlabel('stimulus frequency [EODf]')
ax[0].set_ylabel('MPF [EODf]')
ax[set + e + str(1)].set_ylabel('Modulation ')
#ax[1, 2].set_xlabel('stimulus frequency [EODf]')
#ax[2, 3].set_xlabel('stimulus frequency [EODf]')
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[set+e+str(1)].spines['right'].set_visible(False)
ax[set+e+str(1)].spines['top'].set_visible(False)
#ax[2, ee].spines['right'].set_visible(False)
#ax[2, ee].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
ax[set+e+str(1)].set_xlim([0, 5])
#plt.tight_layout()
# fig.label_axes()
ax[set+e+str(1)].set_ylim([0, lim])
#ax[0].set_ylim([0, 240])
ax[0] = remove_tick_marks(ax[0])
if set != data[-1]:
ax[set+e+str(1)] = remove_tick_marks(ax[set+e+str(1)])
if yticks == True:
remove_tick_ymarks(ax[set+e+str(1)])
remove_tick_ymarks(ax[0])
#ax[0].set_ylim([0, lim])
return y_sum, moduls
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1']
data = ['2019-10-21-aa-invivo-1', '2019-10-21-au-invivo-1', '2019-10-28-aj-invivo-1']
data = ['2019-09-23-ad-invivo-1', '2019-10-21-au-invivo-1', '2019-10-28-aj-invivo-1']
default_settings(data,intermediate_width = 6.7,intermediate_length = 7.5)
#fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True)
ax = {}
plot_single_cells(ax, data = data)
#fig.savefig()
plt.savefig('differentcells_trans.pdf')
plt.savefig('../highbeats_pdf/differentcells_trans.pdf')
# plt.subplots_adjust(left = 0.25)
plt.show()
#plt.close()

428
examples.py Normal file
View File

@ -0,0 +1,428 @@
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 axes import label_axes, labelaxes_params
from myfunctions import auto_rows
#from differentcells import default_settings
#from differentcells import plot_single_cells
import matplotlib.gridspec as gridspec
from functionssimulation import single_stim
import math
from functionssimulation import find_times
from functionssimulation import rectify
from functionssimulation import global_maxima
from functionssimulation import integrate_chirp
from functionssimulation import find_periods
from myfunctions import default_settings
from axes import labelaxes_params,label_axes
from mpl_toolkits.axes_grid1 import host_subplot
import mpl_toolkits.axisartist as AA
import string
def plot_single_cells(ax,colors = ['#BA2D22', '#F47F17', '#AAB71B', '#3673A4', '#53379B'], data = ['2019-10-21-aa-invivo-1','2019-11-18-af-invivo-1','2019-10-28-aj-invivo-1'], var = '05'):
# labelaxes_params(xoffs=-3, yoffs=0, labels='A', font=dict(fontweight='bold'))
data_all = pd.read_pickle('beat_results_smoothed.pkl')
end = ['original', '005','05', '2' ]
end = [var]
y_sum = [[]] * len(data)
axis = {}
for dd,set in enumerate(data):
for ee, e in enumerate(end):
d = data_all[data_all['dataset'] == set]
eod = d['eodf'].iloc[0]
x = d['delta_f'] / d['eodf'] + 1
xx = d['delta_f']
y = d['result_frequency_' + e]
y2 = d['result_amplitude_max_' + e]
y_sum[dd] = np.nanmax(y)
axis[1] = plt.subplot(ax[0])
axis[1].plot(x, y, zorder = 1,color=colors[0])
axis[1].set_ylabel('AF [Hz]')
axis[1].set_xlim([0, 4])
labels = [item.get_text() for item in axis[1].get_xticklabels()]
empty_string_labels = [''] * len(labels)
axis[1].set_xticklabels(empty_string_labels)
axis[2] = host_subplot(ax[1], axes_class=AA.Axes)
#axis[2] = plt.subplot(ax[1])
#host = host_subplot(ax[1], axes_class=AA.Axes)
#host.spines['right'].set_visible(False)
#host.spines['top'].set_visible(False)
#axis[2] = host.twiny()
axis[2].plot(xx, y2, label="Beats [Hz]", zorder = 2,color=colors[0])
axis[2].set_ylabel('Modulation ')
axis[1].spines['right'].set_visible(False)
axis[1].spines['top'].set_visible(False)
axis[2].spines['right'].set_visible(False)
axis[2].spines['top'].set_visible(False)
axis[1].set_xlim([0, np.max(x)])
axis[2].set_xlim([-eod, np.max(xx)])
nr_size = 10
axis[2].text(-0.02, 1.1, string.ascii_uppercase[4],
transform=axis[2].transAxes,
size=nr_size, weight='bold')
axis[1].text(-0.02, 1.1, string.ascii_uppercase[3],
transform=axis[1].transAxes,
size=nr_size, weight='bold')
axis[3] = axis[2].twiny()
axis[3].set_xlabel('EOD multiples')
offset = -40
new_fixed_axis = axis[3].get_grid_helper().new_fixed_axis
axis[3].axis["bottom"] = new_fixed_axis(loc="bottom", axes=axis[3],
offset=(0,offset))
axis[3].spines['right'].set_visible(False)
axis[3].spines['top'].set_visible(False)
axis[3].axis["bottom"].toggle(all=True)
axis[2].set_xlabel("Difference frequency [Hz]")
#par2.set_xlim([np.min(xx), np.max(xx)])
axis[3].set_xlim([0, np.max(x)])
#p1, = host.plot([0, 1, 2], [0, 1, 2], label="Density")
#p2, = par1.plot([0, 1, 2], [0, 3, 2], label="Temperature")
p3, = axis[3].plot(x, y2,color = 'grey',zorder = 1)
#embed()
axis[2].set_xticks(np.arange(-eod,np.max(xx),eod/2))
#ax['corr'].set_yticks(np.arange(eod_fe[0] - eod_fr, eod_fe[-1] - eod_fr, eod_fr / 2))
axis[2].set_ylim([0, np.nanmax(y_sum)])
plt.subplots_adjust(wspace = 0.4,left = 0.17, right = 0.96,bottom = 0.2)
return axis,y2
def plot_beat_corr(ax,lower,beat_corr_col = 'red',df_col = 'pink',ax_nr = 3,multiple = 3):
eod_fr = 500
eod_fe = np.arange(0, eod_fr * multiple, 5)
beats = eod_fe - eod_fr
beat_corr = eod_fe % eod_fr
beat_corr[beat_corr > eod_fr / 2] = eod_fr - beat_corr[beat_corr > eod_fr / 2]
#gs0 = gridspec.GridSpec(3, 1, height_ratios=[4, 1, 1], hspace=0.7)
#plt.figure(figsize=(4.5, 6))
style = 'dotted'
color_v = 'black'
color_b = 'silver'
# plt.subplot(3,1,1)
ax['corr'] = plt.subplot(lower[ax_nr])
np.max(beats) / eod_fr
ax['corr'].set_xticks(np.arange((eod_fe[0]-eod_fr)/eod_fr+1, (eod_fe[-1]-eod_fr)/eod_fr+1,(eod_fr/2)/eod_fr+1))
ax['corr'].set_yticks(np.arange((eod_fe[0]-eod_fr)/eod_fr+1, (eod_fe[-1]-eod_fr)/eod_fr+1,(eod_fr/2)/eod_fr+1))
ax['corr'].set_xticks(np.arange(0,10,0.5))
ax['corr'].set_yticks(np.arange(0, 10, 0.5))
# plt.axvline(x = -250, Linestyle = style,color = color_v)
# plt.axvline(x = 250, Linestyle = style,color = color_v)
# plt.axvline(x = 750, Linestyle = style,color = color_v)
# plt.axvline(x = 1500, Linestyle = style)
# plt.subplot(3,1,2)
plt.xlabel('Beats [Hz]')
plt.ylabel('Difference frequency [Hz]')
#plt.subplot(gs0[1])
if beat_corr_col != 'no':
plt.plot(beats/eod_fr+1, beat_corr/(eod_fr+1), color=beat_corr_col, alpha = 0.7)
plt.ylim([0,np.max(beat_corr/(eod_fr+1))*1.4])
plt.xlim([(beats/eod_fr+1)[0],(beats/eod_fr+1)[-1]])
if df_col != 'no':
plt.plot(beats/eod_fr+1, np.abs(beats)/(eod_fr+1), color=df_col, alpha = 0.7)
#plt.axvline(x=-250, Linestyle=style, color=color_v)
#plt.axvline(x=250, Linestyle=style, color=color_v)
#plt.axvline(x=750, Linestyle=style, color=color_v)
plt.xlabel('EOD adjusted beat [Hz]')
ax['corr'] .spines['right'].set_visible(False)
ax['corr'] .spines['top'].set_visible(False)
ax['corr'] .spines['left'].set_visible(True)
ax['corr'] .spines['bottom'].set_visible(True)
# plt.axvline(x = 1250, Linestyle = style,color = color_v)
# plt.axvline(x = 1500, Linestyle = style,color = color_v)
mult = np.array(beats) / eod_fr + 1
# plt.subplot(3,1,3)
plt.xlabel('EOD multiples')
plt.ylabel('EOD adj. beat [Hz]', fontsize = 10)
plt.grid()
#plt.subplot(gs0[2])
#plt.plot(mult, beat_corr, color=color_b)
# plt.axvline(x = 0, Linestyle = style)
#plt.axvline(x=0.5, Linestyle=style, color=color_v)
# plt.axvline(x = 1, Linestyle = style)
#plt.axvline(x=1.5, Linestyle=style, color=color_v)
#plt.axvline(x=2.5, Linestyle=style, color=color_v)
#plt.xlabel('EOD multiples')
#plt.ylabel('EOD adj. beat [Hz]', fontsize = 10)
return ax
def try_resort_automatically():
diffs = np.diff(dfs)
fast_sampling = dfs[np.concatenate([np.array([True]),diffs <21])]
second_derivative = np.diff(np.diff(fast_sampling))
first_index = np.concatenate([np.array([False]),second_derivative <0])
second_index = np.concatenate([second_derivative > 0,np.array([False])])
remaining = fast_sampling[np.concatenate([np.array([True]),second_derivative == 0, np.array([True])])]
first = np.arange(0,len(first_index),1)[first_index]
second = np.arange(0, len(second_index), 1)[second_index]-1
residual = []
indeces = []
for i in range(len(first)):
index = np.arange(first[i],second[i],2)
index2 = np.arange(first[i], second[i], 1)
indeces.append(index2)
residual.append(fast_sampling[index])#first[i]:second[i]:2
indeces = np.concatenate(indeces)
remaining = fast_sampling[~indeces]
residual = np.concatenate(residual)
new_dfs = np.sort(np.concatenate([residual, remaining]))
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1']
data = ['2019-09-23-ad-invivo-1']
labelaxes_params(xoffs=1, yoffs=0, labels='A', font=dict(fontweight='bold'))
labelaxes_params(xoffs=-6, yoffs=1, labels='A', font=dict(fontweight='bold'))
default_settings(data,intermediate_width = 6.29,intermediate_length = 7.5, ts = 6, ls = 8, fs = 9)
fig = plt.figure()
#fig, ax = plt.subplots(nrows=2, ncols=3, sharex=True)
ax = {}
#ax = plt.subplot(grid[2])
data_all = pd.read_pickle('data_beat.pkl')
d = data_all[data_all['dataset'] == data[0]]
eod = d['eodf'].iloc[0]
dfs = np.unique(d['df'])
#embed()
grid = gridspec.GridSpec(2, 4, wspace=0.0, height_ratios=[6, 2], width_ratios=[1,1,0.3,3], hspace=0.2)
low_nr = 60
from_middle = 45 #20
example_df = [low_nr- eod,eod / 2 - from_middle - eod,eod - low_nr - eod,low_nr,eod / 2 - from_middle, low_nr + eod]
#example_df = [1, eod / 2 - 20 - eod, eod - low_nr - eod, low_nr, eod / 2 - 20, low_nr + eod]
rows = len(example_df)
cols = 1
power_raster = gridspec.GridSpecFromSubplotSpec(rows, cols,
subplot_spec=grid[0, 0],wspace = 0.05, hspace=0.3)
max_p = [[]]*len(example_df)
for i in range(len(example_df)):
power = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[1,1.7],hspace = 0.2, wspace = 0.2, subplot_spec = power_raster[i])
first = ['crimson', 'lightcoral', 'darkviolet']
second = ['hotpink', 'deeppink', 'mediumvioletred']
third = ['khaki', 'yellow', 'gold']
third = ['orange', 'orangered', 'darkred']
fourth = ['DarkGreen', 'LimeGreen', 'YellowGreen']
fith = ['SkyBlue', 'DeepSkyBlue', 'Blue']
colors = np.concatenate([fourth, third, first])
ax_nr = 0
ax['scatter_small'+str(i)] = plt.subplot(power[ax_nr])
eod_fr = eod
eod_fe = [example_df[i] + eod]
e = 0
factor = 200
sampling = 500 * factor
minus_bef = -250
plus_bef = -200
#minus_bef = -2100
#plus_bef = -100
f_max, lims, _ = single_stim(ax, [colors[i]], 1, 1, eod_fr, eod_fe, e, power,delta_t = 0.001, add = 'no',minus_bef =minus_bef, plus_bef = plus_bef,sampling = sampling,
col_basic = 'silver',col_hline = 'no',labels = False,a_fr=1, ax_nr=ax_nr , phase_zero=[0], shift_phase=0,df_col = 'no',beat_corr_col='no', size=[120], a_fe=0.8)
ax['between'] = plt.subplot(grid[0, 2])
ax['between'].spines['right'].set_visible(False)
ax['between'].spines['top'].set_visible(False)
ax['between'].spines['left'].set_visible(False)
ax['between'].spines['bottom'].set_visible(False)
ax['between'].set_ylim([np.min(dfs), np.max(dfs)])
ax['between'].set_xlim([-0.5,30])
ax['between'].set_xticks([])
ax['between'].set_yticks([])
ax['between'].set_ylim(ax['between'].get_ylim()[::-1])
nr_size = 10
ax['scatter'] = plt.subplot(grid[0,1])
ax['scatter'].spines['right'].set_visible(False)
ax['scatter'].spines['top'].set_visible(False)
counter = 0
new_dfs = np.concatenate([dfs[0:25], dfs[25:40:2], dfs[40:53:2], dfs[54:-1]])
for i in range(len(new_dfs)):
spikes = d[d['df'] == new_dfs[i]]['spike_times']
counter += 1
ll = 0.1
ul = 0.3
transformed_spikes = spikes.iloc[0]-spikes.iloc[0][0]
used_spikes = transformed_spikes[transformed_spikes>ll]
used_spikes = used_spikes[used_spikes<ul]*1000
ax['scatter'].scatter(used_spikes,np.ones(len(used_spikes))*new_dfs[i],s = 0.2,color = 'silver')
#plt.gca().invert_yaxis()
#ax = plt.gca()
ax['scatter'].set_ylim([np.min(dfs),np.max(dfs)])
ax['scatter'].set_xlim([ll*1000,ul*1000])
ax['scatter'].set_ylabel('Difference frequency [Hz]')
ax['scatter'].set_xlabel('Time [ms]')
ax['scatter'].set_ylim(ax['scatter'].get_ylim()[::-1])
ax['scatter'].text(-0.1, 1.025, string.ascii_uppercase[0], transform=ax['scatter'].transAxes,
size= nr_size, weight='bold')
#embed()
axis = gridspec.GridSpecFromSubplotSpec(2, 1,
subplot_spec=grid[1,:], wspace=0, hspace=0.5)
x = d['df'] / d['eodf'] + 1
main_color = 'darkgrey'
var = '05'
var = 'original'
axis,y2 = plot_single_cells(axis,colors = [main_color], data = data,var = var)
new_dfs = np.concatenate([dfs[0:25], dfs[25:40:2], dfs[40:53:2], dfs[54:-1]])
low_nr = 60
low = [low_nr- eod, low_nr, low_nr + eod, low_nr + eod * 2, low_nr + eod * 3]
high_nr = eod / 2 - 20
high = [high_nr - eod, high_nr, high_nr + eod, high_nr + eod * 2, high_nr + eod * 3]
high_nr = eod - low_nr
low_mirrowed = [high_nr - eod, high_nr, high_nr + eod, high_nr + eod * 2, high_nr + eod * 3]
first = ['crimson','lightcoral','darkviolet']
second = ['hotpink','deeppink','mediumvioletred']
third = ['khaki','yellow','gold']
third = ['orange','orangered','darkred']
fourth = ['DarkGreen','LimeGreen','YellowGreen']
fith = ['SkyBlue','DeepSkyBlue','Blue']
colors = np.concatenate([fourth, third,first ])
example_df = np.concatenate([first, high, low_mirrowed])
#embed()
new = np.transpose([low, high, low_mirrowed])
example_df = np.concatenate([new[0], new[1]])#new[2]new[2]
rows = len(example_df)
cols = 1
power_raster = gridspec.GridSpecFromSubplotSpec(rows, cols,
subplot_spec=grid[0, 3],wspace = 0.05, hspace=0.3)
#plt.tight_layout(power_raster)
#embed()
low_nr = 60
from_middle = 45 #20
example_df = [low_nr- eod,eod / 2 - from_middle - eod,eod - low_nr - eod,low_nr,eod / 2 - from_middle, low_nr + eod]
#example_df = [1, eod / 2 - 20 - eod, eod - low_nr - eod, low_nr, eod / 2 - 20, low_nr + eod]
max_p = [[]]*len(example_df)
for i in range(len(example_df)):
power = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[1,1.7],hspace = 0.2, wspace = 0.2, subplot_spec = power_raster[i])
ax_nr = 0
ax['scatter_small'+str(i)] = plt.subplot(power[ax_nr])
eod_fr = eod
eod_fe = [example_df[i] + eod]
e = 0
factor = 200
sampling = 500 * factor
minus_bef = -250
plus_bef = -200
#minus_bef = -2100
#plus_bef = -100
f_max, lims, _ = single_stim(ax, [colors[i]], 1, 1, eod_fr, eod_fe, e, power,delta_t = 0.001, add = 'no',minus_bef =minus_bef, plus_bef = plus_bef,sampling = sampling,
col_basic = 'silver',col_hline = 'no',labels = False,a_fr=1, ax_nr=ax_nr , phase_zero=[0], shift_phase=0,df_col = 'no',beat_corr_col='no', size=[120], a_fe=0.8)
ax['between'].scatter(0.12,example_df[i],zorder=2, s=25,marker = '<',color=colors[i])
ax['between'].scatter(0.12,example_df[i], zorder=2, s=25,marker = '<',color=colors[i])
ll = np.abs(plus_bef)
ul = np.abs(minus_bef)
df = new_dfs[np.argmin(np.abs(new_dfs - example_df[i]))]
spikes = d[d['df'] == df]['spike_times']
tranformed_spikes = spikes.iloc[0]*1000-spikes.iloc[0][0]*1000
used_spikes = tranformed_spikes[tranformed_spikes>ll]
used_spikes = used_spikes[used_spikes<ul]
used_spikes = used_spikes-used_spikes[0]
ax['scatter_small'+str(i)].scatter((used_spikes),np.ones(len(used_spikes))*-2,zorder=2,s = 2,marker = '|',color = 'black')#color = colors[i]
ax['scatter_small'+str(i)].set_ylim([-2.5,2.2])
ax['power'+str(i)] = plt.subplot(power[1])
nfft = 4096 #
sampling_rate = 40000
# embed()
pp = [[]]*len(spikes)
for s in range(len(spikes)):
new_spikes = list(map(int, (spikes.iloc[s] - spikes.iloc[s][0]) * sampling_rate))
array = np.zeros(new_spikes[-1] + 2)
array[new_spikes] = 1
array = array*sampling_rate
if var == '05':
window05 = 0.0005 * sampling_rate
array = gaussian_filter(array, sigma=window05) * sampling_rate
pp[s], f = ml.psd(array - np.mean(array), Fs=sampling_rate, NFFT=nfft, noverlap=nfft / 2)
#embed()
p = np.mean(pp, axis = 0)
#embed()
diff = d['eodf'].iloc[0] * (df / d['eodf'].iloc[0] - int(df / d['eodf'].iloc[0]))
if diff > d['eodf'].iloc[0] * 0.5:
diff = diff - d['eodf'].iloc[0] * 0.5
plt.plot(f, p, zorder=1 ,color=main_color)
max_p[i] = np.max(p)
ax['power'+str(i)].scatter(f[np.argmax(p[f < 0.5 * eod])],max(p[f < 0.5 * eod]),zorder=2,color = colors[i], s = 25)
ax['power' + str(i)].scatter(f[f == f[np.argmin(np.abs(f-eod))]], p[f == f[np.argmin(np.abs(f-eod))]]*0.90, zorder=2,
s=25, color = 'darkgrey',edgecolor = 'black')
ax['power' + str(i)].axvline(x = eod/2, color = 'black', linestyle = 'dashed', lw = 0.5)
plt.xlim([-40, 1600])
axis[3].scatter(example_df[i]/(eod)+1, np.sqrt(np.max(p[f < 0.5 * eod])*np.abs(f[0]-f[1])),zorder=3, s=20,marker = 'o',color=colors[i])
axis[1].scatter(example_df[i]/(eod)+1,f[np.argmax(p[f < 0.5 * eod])],zorder=2, s=20,marker = 'o',color=colors[i])
if i != rows-1:
#ax['power'+str(i)].set_xticks([])
#ax['scatter_small'].set_xticks([])
labels = [item.get_text() for item in ax['scatter_small'+str(i)].get_xticklabels()]
empty_string_labels = [''] * len(labels)
ax['scatter_small'+str(i)].set_xticklabels(empty_string_labels)
labels = [item.get_text() for item in ax['power'+str(i)].get_xticklabels()]
empty_string_labels = [''] * len(labels)
ax['power'+str(i)].set_xticklabels(empty_string_labels)
else:
ax['power'+str(i)].set_xlabel('Frequency [Hz]')
ax['scatter_small'+str(i)].set_xlabel('Time [ms]')
ax['power' + str(i)].set_yticks([])
ax['power'+str(i)].spines['left'].set_visible(False)
ax['scatter_small'+str(i)].spines['left'].set_visible(False)
ax['scatter_small'+str(i)].set_yticks([])
ax['power'+str(i)].spines['right'].set_visible(False)
ax['power'+str(i)].spines['top'].set_visible(False)
ax['scatter_small'+str(i)].spines['right'].set_visible(False)
ax['scatter_small'+str(i)].spines['top'].set_visible(False)
for i in range(len(example_df)):
ax['power'+str(i)].set_ylim([0,np.max(max_p)])
ax['power'+str(0)].text(-0.1, 1.1, string.ascii_uppercase[2], transform=ax['power'+str(0)].transAxes,
size= nr_size, weight='bold')
ax['scatter_small'+str(0)].text(-0.1, 1.1, string.ascii_uppercase[1], transform=ax['scatter_small'+str(0)].transAxes,
size= nr_size, weight='bold')
plt.subplots_adjust(left = 0.11, bottom = 0.18, top = 0.94)
#fig.label_axes()
#fig.label_axes()
#embed()
#grid.format(
# xlabel='xlabel', ylabel='ylabel', suptitle=titles[mode],
# abc=True, abcloc='ul',
# grid=False, xticks=25, yticks=5)
plt.savefig('singlecellexample5.pdf')
plt.savefig('../highbeats_pdf/singlecellexample5.pdf')
# plt.subplots_adjust(left = 0.25)
plt.show()
#plt.close()

185
localmaxima.py Normal file
View File

@ -0,0 +1,185 @@
import matplotlib.pyplot as plt
import numpy as np
from IPython import embed
import matplotlib as matplotlib
import math
import scipy.integrate as integrate
from scipy import signal
from scipy.interpolate import interp1d
from scipy.interpolate import CubicSpline
import scipy as sp
import pickle
from scipy.spatial import distance
from myfunctions import *
import time
from matplotlib import gridspec
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.mlab as ml
import scipy.integrate as si
import pandas as pd
from functionssimulation import find_times
from functionssimulation import find_periods
from functionssimulation import integrate_chirp
from functionssimulation import rectify, find_beats,find_dev
from functionssimulation import global_maxima, find_lm, conv
def snip(left_c,right_c,e,g,sampling, deviation_s,d,eod_fr, a_fr, eod_fe,phase_zero,p, size,s, sigma,a_fe,deviation,beat_corr):
time, time_cut, cut = find_times(left_c[g], right_c[g], sampling, deviation_s[d])
eod_fish_r, period_fish_r, period_fish_e = find_periods(time, eod_fr, a_fr, eod_fe, e)
#embed()
eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
eod_rec_down, eod_rec_up = rectify(eod_fish_r, eod_fe_chirp) # rectify
eod_overlayed_chirp = (eod_fish_r + eod_fe_chirp)[cut:-cut]
maxima_values, maxima_index, maxima_interp = global_maxima(period_fish_e, period_fish_r,
eod_rec_up[cut:-cut]) # global maxima
index_peaks, value_peaks, peaks_interp = find_lm(eod_rec_up[cut:-cut]) # local maxima
middle_conv, eod_conv_down, eod_conv_up, eod_conv_downsampled = conv(eod_fr,sampling, cut, deviation[d], eod_rec_up,
eod_rec_down) # convolve
eod_fish_both = integrate_chirp(a_fe, time, eod_fe[e] - eod_fr, phase_zero[p], size[s], sigma)
am_corr_full = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p], size[s],
sigma) # indirect am calculation
_, time_fish, cut_f = find_times(left_c[g], right_c[g], eod_fr, deviation_s[d]) # downsampled through fish EOD
am_corr_ds = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
am_df_ds = integrate_chirp(a_fe, time_fish, eod_fe[e] - eod_fr, phase_zero[p], size[s],
sigma) # indirect am calculation
return time_cut, eod_conv_up, am_corr_full, peaks_interp, maxima_interp, am_corr_ds,am_df_ds,eod_fish_both,eod_overlayed_chirp
def power_func(bef_c, aft_c, win, deviation_s, sigma, sampling, d_ms, beat_corr, size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation, show_figure = False, plot_dist = False, save = False):
results = [[]]*7
for d in range(len(deviation)):
bef_c = 0.3
aft_c = -0.1
left_c, right_c, left_b, right_b, period_distance_c, period_distance_b, _, period, to_cut,exclude,consp_needed,deli,interval = period_calc(
[float('Inf')]*len(beat_corr), 50, win, deviation_s[d], sampling, beat_corr, bef_c, aft_c, 'stim')
save_n = win
for s in range(len(size)):
for p in range(len(phase_zero)):
beats = eod_fe - eod_fr
for e in range(len(eod_fe)):
left_b = [-0.3*sampling]*len(beat_corr)
right_b = [-0.1 * sampling]*len(beat_corr)
time_b, conv_b, am_corr_b, peaks_interp_b, maxima_interp_b,am_corr_ds_b, am_df_ds_b, am_df_b,eod_overlayed_chirp = snip(left_b, right_b,e,e,
sampling, deviation_s,
d, eod_fr, a_fr, eod_fe,
phase_zero, p, size, s,
sigma, a_fe, deviation,
beat_corr)
#time_c, conv_c, am_corr_c, peaks_interp_c, maxima_interp_c,am_corr_ds_c, am_df_ds_c, am_df_c = snip(left_c, right_c,e,e,
# sampling, deviation_s,
# d, eod_fr, a_fr, eod_fe,
# phase_zero, p, size, s,
# sigma, a_fe, deviation,
# beat_corr)
#embed()
nfft = 4096
name = ['conv','df','df ds','corr','corr ds','global max','local max']
var = [conv_b,am_df_b,am_df_ds_b, am_corr_b, am_corr_ds_b,maxima_interp_b,peaks_interp_b ]
samp = [sampling,sampling,eod_fr,sampling,eod_fr,sampling,sampling]
#pp, f = ml.psd(eod_overlayed_chirp - np.mean(eod_overlayed_chirp), Fs=sampling, NFFT=nfft, noverlap=nfft / 2)
for i in range(len(results)):
plot = False
pp, f = ml.psd(var[i] - np.mean(var[i]), Fs=samp[i], NFFT=nfft,
noverlap=nfft / 2)
if plot == True:
plt.figure()
plt.subplot(1,2,1)
plt.plot(var[i])
plt.title(name[i])
plt.subplot(1, 2, 2)
plt.plot(f,pp)
#plt.xlim([0,2000])
plt.show()
#print(results)
#embed()
if type(results[i]) != dict:
results[i] = {}
results[i]['type'] = name[i]
#embed()
results[i]['f'] = list([f[np.argmax(pp[f < 0.5 * eod_fr])]])
results[i]['amp'] = list([np.sqrt(si.trapz(pp, f, np.abs(f[1]-f[0])))])
results[i]['max'] = list([np.sqrt(np.max(pp[f < 0.5 * eod_fr])*np.abs(f[1]-f[0]))])
else:
results[i]['f'].extend(list([f[np.argmax(pp[f < 0.5 *eod_fr])]]))
#embed()
results[i]['amp'].extend(list([np.sqrt(si.trapz(pp, f, np.abs(f[1]-f[0])))]))
results[i]['max'].extend(list([np.sqrt(np.max(pp[f < 0.5 * eod_fr]) * np.abs(f[1] - f[0]))]))
#if save:
# results = pd.DataFrame(results)
# results.to_pickle('../data/power_simulation.pkl')
# np.save('../data/Ramona/power_simulation.npy', results)
return results
def plot_power(results):
plt.rcParams['figure.figsize'] = (3, 3)
plt.rcParams["legend.frameon"] = False
colors = ['black', 'magenta', 'pink', 'orange', 'moccasin', 'red', 'green', 'silver']
colors = ['red','pink']
results = [results[5]]
fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True)
all_max = [[]] * len(results)
#embed()
for i in range(len(results)):
#embed()
#ax[0].set_ylabel(results[i]['type'], rotation=0, labelpad=40, color=colors[i])
ax[0].plot(beats / eod_fr + 1, np.array(results[i]['f']) / eod_fr + 1, color=colors[i])
# plt.title(results['type'][i])
ax[1].plot(beats / eod_fr + 1, np.array(results[i]['amp']), color=colors[0])
ax[1].plot(beats / eod_fr + 1, np.array(results[i]['max']), color=colors[1])
#ax[2].plot(beats / eod_fr + 1, np.array(results[i]['amp']), color=colors[i])
all_max[i] = np.max(np.array(results[i]['amp']))
#for i in range(len(results)):
# ax[2].set_ylim([0, np.max(all_max)])
plt.subplots_adjust(left=0.25)
#ii, jj = np.shape(ax)
ax[0].set_ylabel('MPF')
ax[1].set_ylabel('Modulation depth')
#ax[0, 2].set_title('Modulation depth (same scale)')
for i in range(len(ax)):
ax[1].set_xlabel('EOD multiples')
ax[i].spines['right'].set_visible(False)
ax[i].spines['top'].set_visible(False)
plt.subplots_adjust(bottom = 0.2)
plt.savefig('localmaxima.pdf')
plt.savefig('../highbeats_pdf/localmaxima.pdf')
plt.show()
delta_t = 0.014 # ms
interest_interval = delta_t * 1.2
bef_c = interest_interval / 2
aft_c = interest_interval / 2
sigma = delta_t / math.sqrt((2 * math.log(10))) # width of the chirp
size = [120] # maximal frequency excursion during chirp / 60 or 100 here
phase_zero = [0] # phase when the chirp occured (vary later) / zero at the peak o a beat cycle
phase_zero = np.arange(0,2*np.pi,2*np.pi/10)
eod_fr = 500 # eod fish reciever
a_fr = 1 # amplitude fish reciever
amplitude = a_fe = 0.2 # amplitude fish emitter
factor = 200
sampling = eod_fr * factor
sampling_fish = 500
#start = 0
#end = 2500
#step = 10
start = 510
end = 3500
step = 500
win = 'w2'
d = 1
x = [ 1.5, 2.5,0.5,]
x = [ 1.5]
time_range = 200 * delta_t
deviation_ms, deviation_s, deviation_dp = find_dev(x, sampling)
start = 5
end = 2500
step = 25
eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
results = power_func( bef_c, aft_c, 'w2', deviation_s, sigma, sampling, deviation_ms, beat_corr, size, [phase_zero[0]], delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation_dp, show_figure = True, plot_dist = False, save = True)
plot_power(results)

522
rotated_singlethree.py Normal file
View File

@ -0,0 +1,522 @@
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
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
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,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])
ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,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] = remove_tick_marks(ax[0])
#ax[0].set_ylim([0, 2000])
wide = 2
#embed()
nr = 1
for i in range(len(sigma)):
ax[i+nr] = plt.subplot(grid[i+nr])
plot_filter(ax, i+nr, 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+nr].set_ylim([0, eodf*1.5])
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))-1].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 range(len(df)):
ax[i].spines['right'].set_visible(False)
ax[i].spines['top'].set_visible(False)
cols = grid.ncols
rows = grid.nrows
ax[int(len(df))-1].set_xlabel(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(len(df)):
ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, 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')
ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
#embed()
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 * 1.5])
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')
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,data = ['2019-10-21-aa-invivo-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
half_page_width = 7.9 / inch_factor
intermediate_width = 12 / inch_factor
whole_page_width = 16 * 2 / inch_factor
small_length = 6 / inch_factor
intermediate_length = 12 * 1.5 / inch_factor
max_length = 25 / inch_factor
whole_page_width = 6.7
intermediate_length = 3.7
#plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
plt.rcParams['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
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( 1,2,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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])
ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
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='dashed')
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
#embed()
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
#ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
#ax[0].set_ylabel('MPF [EODf]')
#ax[0].set_ylabel('Modulation ')
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
ax[1] = plt.subplot(plots[1])
if e != sigma[-1]:
ax[1] = remove_tick_marks(ax[1])
else:
ax[1].set_ylabel('Modulation ')
ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
if e == 'original':
whole_height = height
whole_height = np.max(y2)
if (e != 'whole') and (e!= 'original'):
ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
color = 'black')
#ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
#if i == len(sigma)-1:
# ax[1].set_xlabel('stimulus frequency [EODf]')
ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
ax[1].spines['top'].set_visible(False)
ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].spines['left'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
#ax[1].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
#embed()
ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
#fig.tight_layout()
# fig.label_axes()
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1']
#fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
trans = False
sigma = [0.0005, 0.002] # 0.00005,0.00025,
if trans == True:
col = 1
row = 2
col_small = len(sigma)+2
row_small = 1
l = 6
t = 'horizontal'
wd = [1]
hd = [1,2.5]
left = 1
right = 0
else:
col = 2
row = 1
row_small = len(sigma)+1
col_small = 1
t = 'vertical'
l = 9
wd = [2, 4]
hd = [1]
left = 0
right = 1
default_settings(data, intermediate_width=7.4, intermediate_length=9, ts=6, ls=10, fs=10)
grid = gridspec.GridSpec(row, col, right = 0.95,wspace=0.02,top = 0.95,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
hs = 0.1
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid[0,0], wspace=0.15, hspace=hs)
wish_df = 240#324
color_eod = 'darkgreen'#'orange'
color_stim = 'navy'
color_df = 'orange'#'green'
colors = ['brown']
colors_mpf = ['red']
colors_mod = ['steelblue']
line_col = 'black'
plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],fc = 'lightgrey',line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
#plt.show()
#embed()
#fig.savefig()
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid[left,right], wspace=0.25, hspace=hs)
#embed()
plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025','whole',
plt.savefig('rotatedps_singlethree.pdf')
plt.savefig('../highbeats_pdf/rotatedps_singlethree.pdf')
#plt.savefig('.pdf')
plt.show()
# plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
# plt.show()
# plt.close()
# embed()
# plot_single_tublerones() # original beat_activity

521
rotatedps_single.py Normal file
View File

@ -0,0 +1,521 @@
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
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
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,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])
ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale, ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
ax[1] = remove_tick_marks(ax[1])
#ax[0].set_ylim([0, 2000])
wide = 2
#embed()
for i in range(len(sigma)):
ax[i+2] = plt.subplot(grid[i+2])
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*1.5])
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))-1].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 range(len(df)+1):
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(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(len(df)+1):
ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, 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')
ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
#embed()
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 * 1.5])
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')
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,data = ['2019-10-21-aa-invivo-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
half_page_width = 7.9 / inch_factor
intermediate_width = 12 / inch_factor
whole_page_width = 16 * 2 / inch_factor
small_length = 6 / inch_factor
intermediate_length = 12 * 1.5 / inch_factor
max_length = 25 / inch_factor
whole_page_width = 6.7
intermediate_length = 3.7
#plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
plt.rcParams['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
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( 1,2,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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])
ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
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='dashed')
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
#embed()
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
#ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
#ax[0].set_ylabel('MPF [EODf]')
#ax[0].set_ylabel('Modulation ')
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
ax[1] = plt.subplot(plots[1])
if e != sigma[-1]:
ax[1] = remove_tick_marks(ax[1])
else:
ax[1].set_ylabel('Modulation ')
ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
if e == 'original':
whole_height = height
whole_height = np.max(y2)
if (e != 'whole') and (e!= 'original'):
ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
color = 'black')
#ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
#if i == len(sigma)-1:
# ax[1].set_xlabel('stimulus frequency [EODf]')
ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
ax[1].spines['top'].set_visible(False)
ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].spines['left'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
#ax[1].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
#embed()
ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
#fig.tight_layout()
# fig.label_axes()
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1']
#fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
trans = False
sigma = [0.0005, 0.002] # 0.00005,0.00025,
if trans == True:
col = 1
row = 2
col_small = len(sigma)+2
row_small = 1
l = 6
t = 'horizontal'
wd = [1]
hd = [1,2.5]
left = 1
right = 0
else:
col = 2
row = 1
row_small = len(sigma)+2
col_small = 1
t = 'vertical'
l = 9
wd = [2, 4]
hd = [1]
left = 0
right = 1
default_settings(data, intermediate_width=7.4, intermediate_length=9, ts=6, ls=10, fs=10)
grid = gridspec.GridSpec(row, col, right = 0.95,wspace=0.02,top = 0.97,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
hs = 0.1
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid[0,0], wspace=0.15, hspace=hs)
wish_df = 240#324
color_eod = 'darkgreen'#'orange'
color_stim = 'navy'
color_df = 'orange'#'green'
colors = ['brown']
colors_mpf = ['red']
colors_mod = ['steelblue']
line_col = 'black'
plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],fc = 'lightgrey',line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
#plt.show()
#embed()
#fig.savefig()
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid[left,right], wspace=0.25, hspace=hs)
#embed()
plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['whole','original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025',
plt.savefig('rotatedps_single.pdf')
plt.savefig('../highbeats_pdf/rotatedps_single.pdf')
#plt.savefig('.pdf')
plt.show()
# plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
# plt.show()
# plt.close()
# embed()
# plot_single_tublerones() # original beat_activity

744
rotatedps_singleall.py Normal file
View File

@ -0,0 +1,744 @@
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
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
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)
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_trans(grid,colors = ['brown'],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
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['lines.markersize'] = 6
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
ax = {}
fc = 'lightgrey'
ec = 'grey'
#fc = 'moccasin'
#ec = 'wheat'
scale = 1
ax = plot_whole_ps_trans(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].set_xlim([0, 1000])
ax[0] = remove_tick_marks(ax[0])
ax = plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale, ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
ax[1] = remove_tick_marks(ax[1])
ax[1].set_xlim([0, 1000])
wide = 2
#embed()
for i in range(len(sigma)):
ax[i+2] = plt.subplot(grid[i+2])
plot_filter_trans(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*1.5])
ax[i + 2].set_xlim([0, 1000])
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))-1].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 range(len(df)+1):
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(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(1,len(df)+1):
ax[i].axvline(x = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_example_ps(grid,colors = ['brown'],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 = {}
fc = 'lightgrey'
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])
ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale, ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
ax[1] = remove_tick_marks(ax[1])
#ax[0].set_ylim([0, 2000])
wide = 2
#embed()
for i in range(len(sigma)):
ax[i+2] = plt.subplot(grid[i+2])
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*1.5])
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))-1].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 range(len(df)+1):
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(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(1,len(df)+1):
ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df, 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( f[nr],p[nr], color=colors[0])
ax[ax_nr].fill_between( [f[0][-1],f[0][-1]],[np.min(p), np.max(p)], color=fc,edgecolor=ec)
ax[ax_nr].plot(df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
color=color_df, marker='o', linestyle='None', label='Df')
ax[ax_nr].plot(df[nr] + eodf,p[nr][int((df[nr] + eodf) / stepsize) + 1], color=color_stim, marker='o',
linestyle='None',
label='stimulus')
ax[ax_nr].plot(eodf - 1,np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, 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')
ax[ax_nr].plot(f[nr][f[nr]<eodf/2],p[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(f[nr][f[nr] > eodf / 2],np.zeros(len(f[nr][f[nr] > eodf / 2])), color=colors[0])
#embed()
ax[ax_nr].fill_between( [eodf/2,eodf/2],[np.min(p),np.max(p)], color=fc,edgecolor=ec)
ax[ax_nr].plot( df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
ax[ax_nr].plot(df[nr] + eodf,0, color=color_stim, marker='o',
linestyle='None',
label='stimulus',zorder = 2)#,edgecolors = 'black'
ax[ax_nr].plot(eodf - 1,0, color=color_eod,
marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
#ax[ax_nr].set_ylim([0, eodf * 1.5])
#ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
return ax
def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, 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')
ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
#embed()
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')
ax[ax_nr].set_ylim([0, eodf * 1.5])
ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
return ax
def plot_filter_trans(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(f, p4[array_nr],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([np.abs(df), np.abs(df)],[prev_height, now_height+440], color = 'black')
ax[ax_nr].scatter( np.abs(df), now_height+440, marker = 'v', color='black', zorder = 2)
#embed()
ax[ax_nr].fill_between(f, max(p4[0]) * gauss3 ** 2, facecolor=fc, edgecolor=ec)
ax[ax_nr].plot( eodf, np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale,color=color_eod, marker='o',
linestyle='None')
ax[ax_nr].plot( abs(df),np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,
color=color_df, marker='o', linestyle='None')
#ax[ax_nr].plot(df + eodf,
# np.max(df + eodf,p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,
# color=color_stim, marker='o', linestyle='None')
#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')
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_modul( grid,data = ['2019-10-21-aa-invivo-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['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['lines.markersize'] = 8
plt.rcParams['legend.loc'] = 'upper right'
plt.rcParams["legend.frameon"] = False
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( 1,1,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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']
ax[0] = plt.subplot(plots[0])
eod = d['eodf'].iloc[0]
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
ax[0] = plt.subplot(plots[0])
if e != sigma[-1]:
ax[0] = remove_tick_marks(ax[0])
ax[0].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
if e == 'original':
whole_height = height
whole_height = np.max(y2)
if (e != 'whole') and (e!= 'original'):
ax[0].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
ax[0].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
color = 'black')
ax[0].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
if i == len(sigma)-1:
ax[0].set_xlabel('stimulus frequency [EODf]')
ax[0].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
#ax[1].spines['top'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].spines['left'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
#ax[0].set_xlim([0, 5])
#embed()
#ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-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
half_page_width = 7.9 / inch_factor
intermediate_width = 12 / inch_factor
whole_page_width = 16 * 2 / inch_factor
small_length = 6 / inch_factor
intermediate_length = 12 * 1.5 / inch_factor
max_length = 25 / inch_factor
whole_page_width = 6.7
intermediate_length = 3.7
#plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
plt.rcParams['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
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( 1,1,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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])
if e != 'whole':
ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
#embed()
eod = d['eodf'].iloc[0]
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
#embed()
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
#ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
#ax[0].set_ylabel('MPF [EODf]')
#ax[0].set_ylabel('Modulation ')
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
#ax[1] = plt.subplot(plots[1])
#if e != sigma[-1]:
# ax[1] = remove_tick_marks(ax[1])
#ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
#height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
#if e == 'original':
# whole_height = height
# whole_height = np.max(y2)
#if (e != 'whole') and (e!= 'original'):
#ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
#ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
#color = 'black')
#ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
#if i == len(sigma)-1:
# ax[1].set_xlabel('stimulus frequency [EODf]')
#ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
#ax[1].spines['top'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].spines['left'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
#ax[1].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
#embed()
ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
#fig.tight_layout()
# fig.label_axes()
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1']
#fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
trans = False
sigma = [0.0005, 0.002] # 0.00005,0.00025,
if trans == True:
col = 1
row = 2
col_small = len(sigma)+2
row_small = 1
l = 6
t = 'horizontal'
wd = [1]
hd = [1,2.5]
left = 1
right = 0
else:
col = 2
row = 1
row_small = len(sigma)+2
col_small = 1
t = 'vertical'
l = 9
wd = [3, 4]
hd = [1]
left = 0
right = 1
default_settings(data, intermediate_width=6.7, intermediate_length=9, ts=6, ls=10, fs=10)
grid1 = gridspec.GridSpec(1, 2, left=0.15, right=0.95, wspace=0.25, height_ratios=hd, width_ratios=[1,1],
hspace=0.2) # ,
grid2 = gridspec.GridSpecFromSubplotSpec(row, col, subplot_spec=grid1[0], wspace=0.02,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
hs = 0.1
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid2[0,0], wspace=0.15, hspace=hs)
wish_df = 240#324
color_eod = 'darkgreen'#'orange'
color_stim = 'navy'
color_df = 'orange'#'green'
colors = ['brown']
colors_mpf = ['brown']
colors_mod = ['steelblue']
line_col = 'black'
plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
#plt.show()
#embed()
#fig.savefig()
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid2[left,right], wspace=0.25, hspace=hs)
#embed()
plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['whole','original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025',
grid3 = gridspec.GridSpecFromSubplotSpec(row, col,subplot_spec=grid1[1], wspace=0.02, height_ratios=hd,
width_ratios=wd, hspace=0.2) # ,
hs = 0.1
axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
subplot_spec=grid3[0, 0], wspace=0.15, hspace=hs)
wish_df = 240 # 324
color_eod = 'darkgreen' # 'orange'
color_stim = 'navy'
color_df = 'orange' # 'green'
colors = ['brown']
colors_mpf = ['brown']
colors_mod = ['steelblue']
line_col = 'black'
plot_example_ps_trans(axis, input=['2019-10-21-aa-invivo-1'], line_col=line_col, colors=colors, sigma=sigma,
wish_df=wish_df, color_eod=color_eod, color_stim=color_stim, color_df=color_df)
# plt.show()
# embed()
# fig.savefig()
axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
subplot_spec=grid3[left, right], wspace=0.25, hspace=hs)
# embed()
plot_mean_cells_modul(axis, data=['2019-10-21-aa-invivo-1'], line_col=line_col, lw=1.23, size=22,
sigma=['whole', 'original', '05', '2'], colors=colors_mpf, wish_df=wish_df, color_eod=color_eod,
color_df=color_df, color_modul=colors_mod) # '005','025',
plt.savefig('rotatedps_singleall.pdf')
plt.savefig('../highbeats_pdf/rotatedps_singleall.pdf')
#plt.savefig('.pdf')
plt.show()
# plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
# plt.show()
# plt.close()
# embed()
# plot_single_tublerones() # original beat_activity

722
rotatedps_singleamodul.py Normal file
View File

@ -0,0 +1,722 @@
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
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
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)
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_trans(grid,colors = ['brown'],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
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['lines.markersize'] = 6
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
ax = {}
fc = 'lightgrey'
ec = 'grey'
#fc = 'moccasin'
#ec = 'wheat'
scale = 1
ax = plot_whole_ps_trans(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].set_xlim([0, 1000])
ax[0] = remove_tick_marks(ax[0])
ax = plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale, ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
ax[1] = remove_tick_marks(ax[1])
ax[1].set_xlim([0, 1000])
wide = 2
#embed()
for i in range(len(sigma)):
ax[i+2] = plt.subplot(grid[i+2])
plot_filter_trans(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*1.5])
ax[i + 2].set_xlim([0, 1000])
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_xlabel('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 range(len(df)+1):
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_ylabel(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(1,len(df)+1):
ax[i].axvline(x = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_example_ps(grid,colors = ['brown'],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 = {}
fc = 'lightgrey'
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])
ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale, ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
ax[1] = remove_tick_marks(ax[1])
#ax[0].set_ylim([0, 2000])
wide = 2
#embed()
for i in range(len(sigma)):
ax[i+2] = plt.subplot(grid[i+2])
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*1.5])
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))-1].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 range(len(df)+1):
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(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(1,len(df)+1):
ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_whole_ps_trans(f,ax,grid, colors, eodf, stepsize, p, df, 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( f[nr],p[nr], color=colors[0])
ax[ax_nr].fill_between( [f[0][-1],f[0][-1]],[np.min(p), np.max(p)], color=fc,edgecolor=ec)
ax[ax_nr].plot(df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
color=color_df, marker='o', linestyle='None', label='Df')
ax[ax_nr].plot(df[nr] + eodf,p[nr][int((df[nr] + eodf) / stepsize) + 1], color=color_stim, marker='o',
linestyle='None',
label='stimulus')
ax[ax_nr].plot(eodf - 1,np.max(p[nr][int(eodf / stepsize) - 5:int(eodf / stepsize) + 5]) * scale, 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')
ax[ax_nr].plot(f[nr][f[nr]<eodf/2],p[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(f[nr][f[nr] > eodf / 2],np.zeros(len(f[nr][f[nr] > eodf / 2])), color=colors[0])
#embed()
ax[ax_nr].fill_between( [eodf/2,eodf/2],[np.min(p),np.max(p)], color=fc,edgecolor=ec)
ax[ax_nr].plot( df[0],np.max(p[nr][int(abs(df[nr]) / stepsize) - 5:int(abs(df[nr]) / stepsize) + 5]) * scale,
color=color_df, marker='o',zorder = 2, linestyle='None', label='Df')#edgecolors = 'black'
ax[ax_nr].plot(df[nr] + eodf,0, color=color_stim, marker='o',
linestyle='None',
label='stimulus',zorder = 2)#,edgecolors = 'black'
ax[ax_nr].plot(eodf - 1,0, color=color_eod,
marker='o', linestyle='None', label='EODf',zorder = 2) #edgecolors = 'black', # = '+str(int(eodf))+' Hz')
#ax[ax_nr].set_ylim([0, eodf * 1.5])
#ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
return ax
def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, 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')
ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
#embed()
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')
ax[ax_nr].set_ylim([0, eodf * 1.5])
ax[ax_nr].set_xlim(ax[ax_nr].get_xlim()[::-1])
return ax
def plot_filter_trans(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(f, p4[array_nr],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([np.abs(df), np.abs(df)],[prev_height, now_height+440], color = 'black')
ax[ax_nr].scatter( np.abs(df), now_height+440, marker = 'v', color='black', zorder = 2)
#embed()
ax[ax_nr].fill_between(f, max(p4[0]) * gauss3 ** 2, facecolor=fc, edgecolor=ec)
ax[ax_nr].plot( eodf, np.max(p4[array_nr][int(eodf / stepsize) - wide:int(eodf / stepsize) + wide]) * scale,color=color_eod, marker='o',
linestyle='None')
ax[ax_nr].plot( abs(df),np.max(p4[array_nr][int(abs(df) / stepsize) - wide:int(abs(df) / stepsize) + wide]) * scale,
color=color_df, marker='o', linestyle='None')
#ax[ax_nr].plot(df + eodf,
# np.max(df + eodf,p4[array_nr][int((df + eodf) / stepsize) - wide:int((df + eodf) / stepsize) + wide]) * scale,
# color=color_stim, marker='o', linestyle='None')
#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')
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_modul( grid,data = ['2019-10-21-aa-invivo-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['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['lines.markersize'] = 8
plt.rcParams['legend.loc'] = 'upper right'
plt.rcParams["legend.frameon"] = False
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( 1,1,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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']
ax[0] = plt.subplot(plots[0])
eod = d['eodf'].iloc[0]
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
ax[0] = plt.subplot(plots[0])
if e != sigma[-1]:
ax[0] = remove_tick_marks(ax[0])
ax[0].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
if e == 'original':
whole_height = height
whole_height = np.max(y2)
if (e != 'whole') and (e!= 'original'):
ax[0].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
ax[0].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
color = 'black')
ax[0].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
if i == len(sigma)-1:
ax[0].set_xlabel('stimulus frequency [EODf]')
ax[0].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
#ax[1].spines['top'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
#ax[0].spines['left'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
#ax[0].set_xlim([0, 5])
#embed()
#ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
def plot_mean_cells( grid,data = ['2019-10-21-aa-invivo-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
half_page_width = 7.9 / inch_factor
intermediate_width = 12 / inch_factor
whole_page_width = 16 * 2 / inch_factor
small_length = 6 / inch_factor
intermediate_length = 12 * 1.5 / inch_factor
max_length = 25 / inch_factor
whole_page_width = 6.7
intermediate_length = 3.7
#plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
plt.rcParams['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
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( 1,1,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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])
if e != 'whole':
ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
#embed()
eod = d['eodf'].iloc[0]
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
#embed()
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
#ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
#ax[0].set_ylabel('MPF [EODf]')
#ax[0].set_ylabel('Modulation ')
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
#ax[1] = plt.subplot(plots[1])
#if e != sigma[-1]:
# ax[1] = remove_tick_marks(ax[1])
#ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
#height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
#if e == 'original':
# whole_height = height
# whole_height = np.max(y2)
#if (e != 'whole') and (e!= 'original'):
#ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
#ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
#color = 'black')
#ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
#if i == len(sigma)-1:
# ax[1].set_xlabel('stimulus frequency [EODf]')
#ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
#ax[1].spines['top'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].spines['left'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
#ax[1].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
#embed()
ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
#fig.tight_layout()
# fig.label_axes()
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1']
#fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
trans = False
sigma = [0.0005, 0.002] # 0.00005,0.00025,
if trans == True:
col = 1
row = 2
col_small = len(sigma)+2
row_small = 1
l = 6
t = 'horizontal'
wd = [1]
hd = [1,2.5]
left = 1
right = 0
else:
col = 2
row = 1
row_small = len(sigma)+2
col_small = 1
t = 'vertical'
l = 9
wd = [3, 4]
hd = [1]
left = 0
right = 1
default_settings(data, intermediate_width=6.7, intermediate_length=9, ts=6, ls=10, fs=10)
grid1 = gridspec.GridSpec(1, 1, left=0.15, right=0.95, wspace=0.25, height_ratios=hd,
hspace=0.2) # ,
grid3 = gridspec.GridSpecFromSubplotSpec(row, col,subplot_spec=grid1[0], wspace=0.35, height_ratios=hd,
width_ratios=wd, hspace=0.2) # ,
hs = 0.1
axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
subplot_spec=grid3[0, 0], wspace=0.15, hspace=hs)
wish_df = 240 # 324
color_eod = 'darkgreen' # 'orange'
color_stim = 'navy'
color_df = 'orange' # 'green'
colors = ['brown']
colors_mpf = ['brown']
colors_mod = ['brown']
line_col = 'black'
plot_example_ps_trans(axis, input=['2019-10-21-aa-invivo-1'], line_col=line_col, colors=colors, sigma=sigma,
wish_df=wish_df, color_eod=color_eod, color_stim=color_stim, color_df=color_df)
# plt.show()
# embed()
# fig.savefig()
axis = gridspec.GridSpecFromSubplotSpec(row_small, col_small,
subplot_spec=grid3[left, right], wspace=0.25, hspace=hs)
# embed()
plot_mean_cells_modul(axis, data=['2019-10-21-aa-invivo-1'], line_col=line_col, lw=1.23, size=22,
sigma=['whole', 'original', '05', '2'], colors=colors_mpf, wish_df=wish_df, color_eod=color_eod,
color_df=color_df, color_modul=colors_mod) # '005','025',
plt.savefig('rotatedps_singleamodul.pdf')
plt.savefig('../highbeats_pdf/rotatedps_singleamodul.pdf')
#plt.savefig('.pdf')
plt.show()
# plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
# plt.show()
# plt.close()
# embed()
# plot_single_tublerones() # original beat_activity

521
rotatedps_singlesingle.py Normal file
View File

@ -0,0 +1,521 @@
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
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
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,colors = ['brown'],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 = {}
fc = 'lightgrey'
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])
ax = plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df,scale = scale, ax_nr = 1,nr = 0, filter = 'original',color_eod = color_eod,color_stim = color_stim , color_df = color_df,fc = fc, ec = ec)
ax[1] = remove_tick_marks(ax[1])
#ax[0].set_ylim([0, 2000])
wide = 2
#embed()
for i in range(len(sigma)):
ax[i+2] = plt.subplot(grid[i+2])
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*1.5])
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))-1].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 range(len(df)+1):
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(' power spectral density [Hz²/Hz]')
#ax[2].set_ylabel('Hz²/Hz')
#ax[3].set_ylabel('Hz²/Hz')
#ax[0].set_ylabel('Hz²/Hz')
for i in range(1,len(df)+1):
ax[i].axhline(y = eodf/2, color = line_col, linestyle = 'dashed')
plt.tight_layout()
#embed()
#fig.label_axes()
def plot_whole_ps(f,ax,grid, colors, eodf, stepsize, p, df, 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')
ax[ax_nr].plot(p[nr][f[nr]<eodf/2], f[nr][f[nr]<eodf/2], color=colors[0])
ax[ax_nr].plot(np.zeros(len(f[nr][f[nr] > eodf / 2])), f[nr][f[nr] > eodf / 2], color=colors[0])
#embed()
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 * 1.5])
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')
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,data = ['2019-10-21-aa-invivo-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
half_page_width = 7.9 / inch_factor
intermediate_width = 12 / inch_factor
whole_page_width = 16 * 2 / inch_factor
small_length = 6 / inch_factor
intermediate_length = 12 * 1.5 / inch_factor
max_length = 25 / inch_factor
whole_page_width = 6.7
intermediate_length = 3.7
#plt.rcParams['figure.figsize'] = (whole_page_width, intermediate_length)
plt.rcParams['font.size'] = 11
plt.rcParams['axes.titlesize'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['lines.linewidth'] = 1.5
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( 1,1,
subplot_spec=grid[i], wspace=0.4, hspace=0.5)
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])
if e != 'whole':
ax[0].plot(ff, fe, color='grey', zorder = 1, linestyle='--', linewidth = lw)
ax[0].axhline(y=eod / 2, color=line_col, linestyle='dashed')
ax[0].plot(x, y, color=colors[0], zorder = 2,linewidth = lw)
#embed()
eod = d['eodf'].iloc[0]
if np.max(y)<d['eodf'].iloc[0]*0.6:
color_chosen = color_df
else:
color_chosen = color_eod
#embed()
x_scatter = x.iloc[np.argmin(np.abs(np.array(x) - (wish_df / eod + 1)))]
#ax[0].scatter(x_scatter, y.iloc[np.argmin(np.abs(np.array(x) - (wish_df/eod+1)))], zorder=3, color=color_chosen, s = size)
#ax[0].set_ylabel('MPF [EODf]')
#ax[0].set_ylabel('Modulation ')
#ax[0, dd].set_title(e + ' ms')
ax[0].set_xlim([0, 4])
#ax[1] = plt.subplot(plots[1])
#if e != sigma[-1]:
# ax[1] = remove_tick_marks(ax[1])
#ax[1].plot(x, y2, color=color_modul[0],zorder = 1, linewidth = lw)
#height = y2.iloc[np.argmin(np.abs(np.array(x)-(wish_df/eod+1)))]
#if e == 'original':
# whole_height = height
# whole_height = np.max(y2)
#if (e != 'whole') and (e!= 'original'):
#ax[1].scatter(x_scatter, height+70, zorder=2, marker = 'v', color='black', s=size)
#embed()
#ax[1].plot([x_scatter,x_scatter], [whole_height,height +70], zorder=3,
#color = 'black')
#ax[1].scatter(x_scatter,height,zorder = 2, color = color_chosen, s =size)
#y_all[i] = np.max(y2)
#if i == len(sigma)-1:
# ax[1].set_xlabel('stimulus frequency [EODf]')
#ax[1].set_ylim([np.min(y_min)*0.8, np.max(y_max)*1.2])
#ax[1].spines['top'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
ax[0].spines['right'].set_visible(False)
ax[0].spines['top'].set_visible(False)
ax[0].spines['left'].set_visible(False)
#ax[1].spines['right'].set_visible(False)
#ax[1].spines['top'].set_visible(False)
ax[0].set_xlim([0, 5])
#embed()
ax[0].set_ylim([0, d['eodf'].iloc[0]*1.5])
ax[0].set_yticks([])
# fig.tight_layout()
# fig.label_axes()
print(sigma[i])
#embed()
plt.subplots_adjust(bottom = 0.13)
#fig.tight_layout()
# fig.label_axes()
if __name__ == "__main__":
data = ['2019-10-21-aa-invivo-1']
#fig, ax = plt.subplots(nrows=5, sharex=True, sharey=True)
trans = False
sigma = [0.0005, 0.002] # 0.00005,0.00025,
if trans == True:
col = 1
row = 2
col_small = len(sigma)+2
row_small = 1
l = 6
t = 'horizontal'
wd = [1]
hd = [1,2.5]
left = 1
right = 0
else:
col = 2
row = 1
row_small = len(sigma)+2
col_small = 1
t = 'vertical'
l = 9
wd = [3, 4]
hd = [1]
left = 0
right = 1
default_settings(data, intermediate_width=6.7, intermediate_length=9, ts=6, ls=10, fs=10)
grid = gridspec.GridSpec(row, col, left = 0.15,right = 0.95,wspace=0.02,height_ratios = hd, width_ratios=wd, hspace=0.2)#,
hs = 0.1
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid[0,0], wspace=0.15, hspace=hs)
wish_df = 240#324
color_eod = 'darkgreen'#'orange'
color_stim = 'navy'
color_df = 'orange'#'green'
colors = ['brown']
colors_mpf = ['brown']
colors_mod = ['steelblue']
line_col = 'black'
plot_example_ps(axis,input = ['2019-10-21-aa-invivo-1'],line_col = line_col,colors = colors,sigma = sigma,wish_df = wish_df,color_eod = color_eod,color_stim = color_stim , color_df = color_df)
#plt.show()
#embed()
#fig.savefig()
axis = gridspec.GridSpecFromSubplotSpec(row_small,col_small,
subplot_spec=grid[left,right], wspace=0.25, hspace=hs)
#embed()
plot_mean_cells(axis, data = ['2019-10-21-aa-invivo-1'],line_col = line_col,lw = 1.23,size = 22, sigma = ['whole','original','05','2'],colors = colors_mpf, wish_df = wish_df,color_eod = color_eod,color_df = color_df, color_modul = colors_mod )#'005','025',
plt.savefig('rotatedps_singlesingle.pdf')
plt.savefig('../highbeats_pdf/rotatedps_singlesingle.pdf')
#plt.savefig('.pdf')
plt.show()
# plt.savefig('../results/Ramona/ma_powerspecs_negative_df' + d + '.pdf')
# plt.show()
# plt.close()
# embed()
# plot_single_tublerones() # original beat_activity

Binary file not shown.