From d758ff8722080a6132b30c3780da83186238a02f Mon Sep 17 00:00:00 2001
From: saschuta <56550328+saschuta@users.noreply.github.com>
Date: Tue, 17 Nov 2020 10:49:04 +0100
Subject: [PATCH] added py files I that are not primary plot versions here
---
differentcells_filter.py | 209 +++++++++++
differentcells_trans.py | 191 ++++++++++
examples.py | 428 ++++++++++++++++++++++
localmaxima.py | 185 ++++++++++
rotated_singlethree.py | 522 ++++++++++++++++++++++++++
rotatedps_single.py | 521 ++++++++++++++++++++++++++
rotatedps_singleall.py | 744 ++++++++++++++++++++++++++++++++++++++
rotatedps_singleamodul.py | 722 ++++++++++++++++++++++++++++++++++++
rotatedps_singlesingle.py | 521 ++++++++++++++++++++++++++
singlecellexample5.pdf | Bin 160758 -> 145493 bytes
10 files changed, 4043 insertions(+)
create mode 100644 differentcells_filter.py
create mode 100644 differentcells_trans.py
create mode 100644 examples.py
create mode 100644 localmaxima.py
create mode 100644 rotated_singlethree.py
create mode 100644 rotatedps_single.py
create mode 100644 rotatedps_singleall.py
create mode 100644 rotatedps_singleamodul.py
create mode 100644 rotatedps_singlesingle.py
diff --git a/differentcells_filter.py b/differentcells_filter.py
new file mode 100644
index 0000000..aeb12e8
--- /dev/null
+++ b/differentcells_filter.py
@@ -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()
diff --git a/differentcells_trans.py b/differentcells_trans.py
new file mode 100644
index 0000000..2c3c186
--- /dev/null
+++ b/differentcells_trans.py
@@ -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()
diff --git a/examples.py b/examples.py
new file mode 100644
index 0000000..59f6cea
--- /dev/null
+++ b/examples.py
@@ -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
ll]
+ used_spikes = used_spikes[used_spikes 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()
diff --git a/localmaxima.py b/localmaxima.py
new file mode 100644
index 0000000..4b0fcdd
--- /dev/null
+++ b/localmaxima.py
@@ -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)
\ No newline at end of file
diff --git a/rotated_singlethree.py b/rotated_singlethree.py
new file mode 100644
index 0000000..afee2e1
--- /dev/null
+++ b/rotated_singlethree.py
@@ -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])
+ #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) 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])
+ #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) 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],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])
+ #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) 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],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])
+ #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) 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])
+ #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)%Ve|x7Or;|I7ZrCHLPfkN^Ik|N7VefAOz>`w#LL{vRRw
z_kaB>Z?_iK_E(JcR$Gu?!X>2bmeN+XaQyTCa{TkM)Z5xedHntT_xAz*`(MoV&Z2F%
zB7d|k>sISu@;|rd&*ESI_y7La|Mb^?`|IESLlo)$wcKinTMnVFMW!8O+Tk$m;m7|z
z{`Nop>%WVtxwoyao>a?i-BMe^Dk8}rFhRVf
zrNj{0x@46=IHcP0)>^992xN#FQxWy1)#6Y?j?0#oXmLobw%FDCAtA#}WK2OM4aUz~
zsbv!jSBpS?mP={kODaLSt+A|C=9rQtqi3rvr)3Q_XI+MJf6MFAHW?vWny#_R99Lms
zaM;8U;&LgPifnNWL1Z*#ShCC`x279YIE0edOyhwf#kgXyOcO&j6OclpISB
zEZtTIQ(TRjf9|h&6`wWr^xCJ>t+ZteoxrqJT-PGWBU?j=u}jFJ1>IlrZ3#7{?$f5G
zUpthf`WRB~$=!^ryw#E@ef-CN{l~w)t7s(AcKr#BvXomB>uo8E34f#xp!+OIZ8>TV
z2`HA4fUE}uMbds7(uGA(tRcurEg>(DKu$amgyj~Oe^R!d`^_=Ls&@iwj8Yf*)eFqU
zS-3x6vILoE$ut$^yvrQSIlNY+1VolDp-S3xy}^{h^~hg)!|u}NEL}bAVP)Eis>@tp
ziXv+hU5g}^kiJ9i1lOnedw+2`c=U0Jd?Ku<|+
zi}6S5=+z%2Al)MYTGtw(S_G7s2ZHssq#$)(D1yO(9S8};Iqplic6pZMVMsvLV<3=4
zo>(ER$1*&rNPm_dfom#n4q0+`vk6Yjff}w^e}dEQW~%wImTo$M7^)U&(2t94l##0V
zEWz*Ht%g;SNQy<*!(m6jG8H#<)+|g@{5CcC4*S|6B{e)G2$r6H?{L15Zfc;u5k0PF
zw_#7gQlx1s&zNsWaIe$!Ve3SXn{tyh-b+)jg&)h~-Yk4;Qm23}x3!v=5Z2@kVH!*P
ze>L^(K$F9-lb~*^FKkw14d42PU6E!b-qMz(zNFOBkJQP~%hPyBkXwpXtQQReZm0)>
zEDG8B=FFjAp5vfJ8WK>20s{F-tbKp4sRT(XR4fv0FNQ$>uzqebpgjZ7QxtpVD1l;e
z;7~OPTvJ(U;>HqXjGGK0>v7Xvs7}9Ie~(>8JAjWFk?c?nX?40u_>vdEsI!N5}nP#b7?Lev$hmx6Cq`xYN~+>;B-*X$0*F9
zYn36yFHH`v9y}SIpVUlSNR3{GDm@!Ez8+rZpMz7ILdGqocKCiVErJ{So~CO{f6|wl
zK706%<2tn&R~Hlc_jM~t+-dn~+{z+ymsV?+J^x6ZoqRI&bM4bLG9c?C0cjqGD4Q6+
z>(Jo99*6{#PEZEwBtV9eL=i(;otwp(X)dWpS}dlJt|mPKiQ?9GM>d@x_D*+gVacmC
z&1dg=@#n6hv;24)b_AvsO$%AHf2DX*97@CchuC6RhECz_wsbTpwv?2fdh)p}gA-R$
z8n#u0IBa}OE0U5b21UrZ^%8126g9V2>{|qUNWcmTYnr?Vg#^@^PUoLR7E^ti*L`;1
zrPB3q2#f2P+zSgg#rcW)^$?Qc7HOf&v84iQwvtL}qN(3k6Vs!p@4KAje=H%n#m?W@
zjr4+$!(A6e#A*^tDX+D+|36X}Z{9$_r)dh9H_Hga*KHhvc5BOMoTcP#?Y(3%IIxo<
z0eKJ(Asd}>rXpa$xNTJm}OI-#nSkap6B68Zs2B#CVI;E
z&(=Z2D&nKDl^kMxT*$0m#;TyXweM^7SXK$NZ+XZeR;jYbBkL#wxyL{dBq7xhhs6nS
z;Jr>9g3Jp;v>;it7dUI0^dTQ+nPn|o*JIK|2DD*|ilpUm)dO|0X6bM^GfPTXI89l;
zGhrrd?RztskQnv
z&4J1fa-iZctp`_zZq3rqEL@erp9MAPvsuvP@CHp%?>L*weHo}VL>V5Nqt&zp_1MUP
zR+iHhi54Ss^+f9ikOOs}h8|XngP)44$k1MRl#_~2e>jY!z_|94K33@Qu3-Jno~D0`
zy`vq7R7rGOORi;Y^+)LDxnKj;NAa)qxYR@r*iAMs3V=b~Z)ocr4mGY?BK&$1>_0dQ
ziUOPkeWbQ#p%Z8>3NlEq5JNTzPRt?{nM6GX8p=Si46|UQ43vSsz0PS;E5s-yB|S9v
zj%;8$e+-pT6p~<1(A>1MfQniqsdfT=h!eBen#9D=MWd>5od=VMwJ-avsL&4Wq07oT
zcpXySRP4(XQg??J4(rT-S3@WORj*lpb)(Lep(p(5kz(H96@1{Imhy8Ud+8!}K+ht43g?CuaPe?9QR6lbw-W0`tz)0Y&&*0Z$9fxg7d
zWd?0KZ`#~$`l8j4m18XqS402%x*b_#r$<@MTThxX=CFj^ex%+ECw<1jX1IA#2{boH
z4pkb-o(`hLf!6fMAxR1B4UtZuS#R@F4shUerof?pZd?aYlnvaRwhe$lY!Le-7hO@x
zelPrM+NKs<@maoMIW9
zw;hmxR?|p8-)TVxE+4b*(6W?|4hNzOsg+T~9kR&*8=`4{+P%TF<(?cYTl(oG#1NA;
zNo}D!D*`tYn4f4xH|GtmKmBPv4FYY8e;M6{%J8haJTsU+3vIDv$#6aYtWA9WxUYd6
zXdN<}J@})B<#DFV(G~zYOyA58mQNr2#Y``{zkFLivq+~NO6l!m?`M`~i&OjD?^q6>*Rs7%xScnw7f7HTJ
zc#Sb_6D~wF2m=j$U=W5f!$GJUrcqW|Shl`AHb-GlC2JJJvRX1H_hAUW>ABEu1r0(A
zR2qaz<459=@+pM^RbUE*w$=P|?B6o&12uOw^})>wNlLq5ABhE5>x=+w27^5pT9(l$
z9SybLTrrZ^sbKKuvVRP*Lym#Yf98WBr2EO)V9$l{*yi1Iulr^|QQ;?T&?=DzAuQet
zv9Ol~-%D-(d1UOxES*Qt(#47#*g1I6FNdg6x?LYn@!dpr4G%G4-JZw#Z8e5
zQ=!#9QE;>DV*}4xoL23M4EMDTl_0C1~e?%Ti>UEJxi3hzb
zLp*rj4|x<>Se0SR1yNdQ&mVb6Qb5;#{X+t$rLTubp(q$A
zQpu8?<1zhukeSG0f2kt&aG#9CV6V0W%_aqJ5(9;jTTx-0kETbFLL${_$P_My=)A_1
zd^D}GRM#Bh@uQ22Y)gU`DXhcQ_fMEsmR=H`Lx7^BybwT1b3cGReEpeMDK}N~#^>Ex
zlys9OA-{;H6M}Lb0gBk40~w%EO+n3QaRohLd)}ExLASQFfBr$;piuY1EH5=}eLRE
zp-Nk>WjP)$`Dfxu>6cN`v~6uW-XENILO(_M#%Jhxl5~;E(T*pS&r14zJyokc$=%n6
zve4|{9Fo9ee=d7ZmB)!~r9)0EhYJHpp%ibRP`&{Jg|dW09vtOB5oLjwJF0v(d2r*4
zJlI))Liu6|6tX&(L?Kf#DAZDujUYDV-~*^77}^zqGn#FbY^%~
zr{T6Pi#UMj)Ewn`-ofN?!V5VuPJ5w}evyYLWZPcpe?$372^49Qs7psYdeRsjyT72Z
z{q14j$4reWr|stBwn)#s4X0gDW4ty*jqP2xnWgwrA@@p4>h__4lcX%cHwXDI0}EQ
z01D;Y?yLhT%Ob-p*);9*FpQxN?P@(a&>GXae;x0;nGR&mMcE@jNx7Dx;kXOI$fv8D
zBq<-5K}q?f%Cw6q`z~lX?zt=%GiIw@GbGFHhX8CJb3{N_{Ovcy(9G9fI1QSdqy@WG^qCrW-$W$%leG6^o`8A{p9T2M@+Jil*Xez(Eh^;
zf4;)I(LkZBXn;cbhzk_Y6<#zcxC*#9oEJu6hg;oj1K5aLKS$Qj&D*1J>{xOb6dH;o
z@=#h49EEogTo^|hg_i=8qiFEiGzvEwI0}a-jzZw}GzuM=#X-1X#z7t_zSx8K3iC(t
zr?QW7Sz)pSk${WakQafXixAf1^AaQnf44Iz2d@KhwgljQL8orWDMb+|lI#2AN95de9}81TqLx&DK4u+ySGxDlD3EmsMo;MNUoOHcoj
zq`Y*p{M4Qc=hpeza?6CaarbS+vn$U
z!cYC#`Z>t=Q+g8p6MlLkLhQ|CToYX9ee^?)i<9|X@Umt%DkMsJlt^WqwzC-K(KeF84
zL-9@@`-ioDViYdE-)fiBIoY0z(t3SBA16j}ee6X3@FjhNC%;9GcUF(@f6pK9@Y8qn
z4wvQlwl;FP^;B@?5yBaH
zD5s>5f_K=Ff>)gLv-K%>Jqfc-X<63+sFx(A4db^@Rz)Op@mozfY?bnHf1r5oup6_@
zvfxayt;5w;l!Lw$Aps|fe`$&nw&yQJNWn)iV6veVVR1Z*j)L+hocWM^2Za;XrzTr-
zhbTUNMYb$!D&32P;b@l#gdtNDc%qLIV2AN^(-c
zBW2J>4W!^JVb{el)x}>MNDMk0Nn#N81<65Wgv>?>4Jin(mc*cM&Ljq9?^!Xl$IH#h
z+M`gm!hqysdNc}Uf6b9X!@fcaMoztue*D?e3Xf(>M|{lLnuAcuq(RnIs@<`*n;*oV
z2_4^=x7+i6wn@QIau?*FPWQrWrA3YF5nz12IGspQ^jYmf-Jm7z*}Cy%9YpE39c6LZ
zXVVY0_V1&u^wawFfteo>&5v%LBZZQVTqr)I;MfUzMk}Tie{2At=vTc`U)nc|;#m-*
zISMCuM8KmPI_6vriYN)K-(DXQtBch~l5$2f&_t^_7o}Z7vvBIYWEb;?+@m0WTI$j~
z>GU1Vlgw@SR~t-0m5hl7p*#0v6*}sK4ytqyE*3I{65V*Y
zU^uR1(+;Y%f3TB*ayVdSCyxCdi7Cek%dJ~k!X2>6vK%1@{X+reCl1=uBRS|c1<66j
zw#Y8D!%21_BzhW!*9DN!L+G!qh(5sjU)CM;{e}Jj_K#iiy4~{qx_QU?`H>nx37-%u
ze|#mnxddSg>BtiCpj+)op~N48&T$MQ4=&3V$Y-!9f0QEN0vi$!`ba@M2m^#Pwgs82
z$YRb#Y3q>3o=Z%vKP_i#o!=sTnM~GKQr!>l
zJzmI?f824+$C__Ax*$S4xYkYcP%8Vp;pmRzq~WwrPiR^xv>~)3fZSZH@y6zSW%F%j
z1j?877G>#kDfcz&CcQ6CXK5Sx!4e@!Jx2|jdQW4_jvoh#G_Zw=xo1a?3a
zTn$kUZayq%Vc&yn%e_Mmf=V;BZtH^!6~;a%y?H-%W1sJBQbh3w4iA&{LXzTbA)t7o
z&lzGUuY;4E{BdjEk6l5PtpgW4mL%my$Ax6QkoV~{5alR(AC)8(P22^wW5L&{FX=mS
ze-wcvBE_>2F}1ENvgUJ?k&(ZxRq63k(jyZ}b0Gb>m_f?#`gr^$x_EaVq2Q9|qG4X(
z@eVoI(9=AWh6hke`B8G@p*SsKp%OSe@!-8B;GwiLh*3*r0$C8QGO^&i=mL9Qa1YL{
z^QPqvU|tT+i_@c=(1$;Y9}F&!ylh)bf5SF4I9_~fn+aSzEwH*uAeM9uMms2%~~PpMWem?cqg{J9u_M8V$Ve*(qP
zT=*+9!qM?&1n}Us6HXT2G`uiZ&S*%%+bYv@m`U#m7+)|<8=)GK|;72lOwvSH58bv)^
zGVi9pPRJTybNot7)_wx6zsT2T$$@w)`4dj~IQUsjiv6RSsgC`I
z)#|&LL7`T6Lr{AV$J5VfR5=`io=$!Lyo2%gA&Z9**6NU5@54yb-I+uHuK
zXIn$-r&6Vqk!>B$&!Zgsf8&8>BJ}2R%8^}-Dc@w#k?FPV3ts4_TmJHkJh-=;k0|b#
zkO!^CEtQ%MbRuK>Nz5zxCs@z}Sn(?%Z2Nv%%l9wk<*!`ot-n^dY+G9QqdOI*&O1y1
zWSC+t+ptckdPY0!ssn!D6GYgcP~!qe;rtDwC~FOCl- hy~{npHdN@aN*{L#f>ma>`li`5Z-$f>ZY3
z$~>Jt{JHQ$NwcF!e?PSJtUnxto=!pUNFQ77w_IM74rK*9&%>9)chh&Ipm&7r8^Yx+
zJdE`5Z7tK%4~H#3h!fH>m0QV6TO_r`2ijX~kOv6n1}T)&y$#gInG<=iI{^hfM2>^R
zr8mCAh9wTd=^1!DF))6?=!eRe;dh`q=w;SwyPwpzUyM>
zRE-RrE${4E4pgp>l4OCn#9_<(A8{Z8g|bu{e6xIu1J&!2zGb#8$1-eoq5nK~!-8!V
zu4m{!oWG}4s_g8cN>jUkbU9qN_-9ngn>)yZ8zP)7f90NcAU?r}v!&-P=51RR$=PG0
ziXOzR?4DKPR4C357a2+B8vIP($l`72!&a{+%eeZaD$ORW9G1!Qk8(bETYEtlS7mU(c?OSfg
zYGpwwnsZS$Mp0wD5WrD5rA_Kodwh~c;a(Rhct1nYe^?kV-Og)_Qwfd|wrxG4@Xq?N
z_M+njb#PyXN{`FhhIf50MBHAj%u(8!YaH<)EJ+PHru1QJ-5KSekIRy(n5w-033z4-6p;8i!|oP9HS&n
zC;qpTe@ZtRlpUaxrSi*KqSD_-mS0_^QCgxqkVM#)So5PLN(k&o3ML9tA5w5dO90u~
zR~w}r7jdQ>)uk|baOy-kxM0yJxL1q>ydohnD7q42RJK;qafOd8p77+hoD@1lLH!?!
zAroEhci~LSaw;JT>UVII67yDuPlWRx#>s|8f9X5t%Mm#I!dxi%(kK$iaYz&w;>n?q
zQQESF;n3Brq}tXUQnc+SIT1qaVSk#AvJrG6j%ep(ErRk6faO$x5>An3A{bsf$UUKr&K3X>kUhT%iO
z{3uG9KnmK4GBswKR_a+95Wt8B=K-YPfB2i!n6`IEw(ZwG%PSryB=y6^^|4T-I(NL1
z=8j^^ao)#FKYh0Lw#LdL$!n-JpU@arMbsEQwDgfQ=JI&9cZNO=K4xc!zsJY+cew1?
z`ugxG{JGHgU((nbJNAc3@uvLYYAi0b)Ns6X@d=G7D;XCSoHW*QuJyRde8Gc!e+onE
z-`4Rg9XLGvEJJc`qT#JnQR3qmQ5vuyi}XWM5!WMTZ;Gyokb+y0=sF})a6v-~F1;|9
zw$&8U(L%_At&9}xu}HywFe#{-1z;*2IjZC%%;&ruC(Po{g*O4F=iEKerre#m$Io~dR_HcbLOQ#ZF`^OvR
z^43fUw*9tVCf`KF5`NY~(&Wj4qdR{GjCzQTV4)P`2z$!8*F|q83ic+8f5MbQ@uq=d
zy_Hltp!ETUh5!r8Jn&Ck|Tsgeeef*T@>!kEf>e^#Wgi{nS9!hOj6WPR`OZJ+zM
znd9B-3$}MX{Vm?%d%CL+<&!EJKjmD;?|<`_b2X-PJM;SBy_-o2vHsd>8g6IJ&sM4L
z^VWYyG`RXu&KH0_It1m?e}+$U(_ZN7VNx^u6JO?2dZKk;UfE^Cb7mN$ZJB_N={h_
z1Ix{zhn`6+*qdCOSfXHmn;(U9$-Eps56(wGkvpi>a7_~_lzN4ue{hte@jJNwHU{9l
z^&(_$zU551%IGkpUbt*J26vL`~$
zRVmkiFLjRd$>SHP;eXa00T!nSivsBFH&S#HiZEig+^gl
za3MJ7C4Dv#$>iZ`>>ZNk(%Le9Mw*sGITit6p=*^|6_3vZj1;jalaTL^=Gsc-2t*g+
z6{`i`Y%7(Mngw|M3m(24d5ftISAj2>-*v+ij@m+K1HQ@ye@Erjw-x**g|bCBKU+T>
z_~iZNplha|p3teEVthGV3S}iTZ87ZOr_Fh#`*Jv=wA)stf(*bo
z^B$%Yt9i@0Uh1>(L&2x;N8$YLBFLlC_p9ccACso>e|@F2^Euy*
zgEmrpTd$OEOFxFy!M)pGIc>M?{+vn@hk1zYA+#cavM_UjVy6dDn&yF`-1^%+ZHyZi
ziev{3fEXVIW1*c4yt#T@RTBW{K&fCqV;PntrrGrnr_8->;Inf{mC#z|F
zIc5_VT$i@Ja8h;U;lj&b;Nee}6CMIY^T&>XBH!faVp-9-ZF?tZms^Ut#O3%v#Jr^b
zC?^7+cT=u=scn}%Q~FN$o0#Mt7X38dOF&zVf0#To%HNaqeK+M!2i&x_O$Iqaa!re&
z>;}yHpzp`$W75o}9TzJ7!n>3cN%J1c>yDF|^xJx+zPs7Ko!+mN+xkC5`)V!Qg9ZT=5a!$gxBB#_g)a=
ze;pd_5Y^+!@p+XJjLo#pl}rR0jS>lKUqjeDe3f#SVj=3!sThhP|w??`xPOx#Y)k
zU%I55tnP9Ov2A<9BPs^bc
zZ!@9
zOd$>5rlvM|
z^PofVcX{ECwQb~51tEVr-jwF?6FJ{aE+;~PD>sHhPPOI%VRu>`)RIjaEy?0TT8D9)
z6?p6ioBn|CTVe!zCfsn%C*xncf2agPZc9nW7pE?a;&ud<4}M=yR{QbGpB>Pg@~ws%
zVy>w?Fr~7K*oze0duJDqcYYMkcr*&tyICGu2;(UHbhxRHVp8BLAp;KX5<6<}aZCQf
zRD6ZeGgu~%YEn>?lfoWvkNwFM7hKF$*rUg=Rqi;cZ?;noMpc9;Fu&sm?M|LPq`ItW!FC_WF(hDy3_b7~3>K7`@YjwP-usCJIEwf8n`4=jB2Mj}+P~*z<#M>V-k*c|v5NxZ@;+n!Q~x*O5YX
zX%}*uDA*rApmTQ+ewCI~%4KFw3>id}p6ig#rU&7(D;F<F_I5b-Y+@AkeV;Q|_jt
z(`S;4h8=F`n1k3;aVTPcX1w3ZwLOqm?L@JIZyMTAlly$g(mu$1e=u>$osU>YP<7ex
zfi&rao7%Ej>M*9Fwo23LP`VgYIz0mYi8Q&ty=*)B@6+V-$zwC@y(FmiV`vq64fP`h
z(9e@}K`?(QBKj(hUVe^n>kZaCTa
zRNaWRY{PN0`BKa+>U`PgN3sWO=?kUU%V2M*#D;TxS;r-Uh%%P;(_i+ED})^G(eZe7
zFDFpWQeqgMG{`Q*FY#anyAZ#~Ln$*JhBy*1NZCjcZeogD?l-_K9(;6##9&`Rb?l|x
zpituzdF-Dse{6N|HW{jP|M+=TrCP>wv8;)LOBv3Udsr7$*>u753w`mP{;Uh~W2NYV
z{nEwO7v-j)N~1ifbhW#5;joDZw`48|IjMAC`zcEZ9f-ktis>iH4d#oiq$)RN~S
zFOL*KdU#of#^Wl|g(5}@l6;IF$&w_jlhGpiKB^M5^h9e!BHC(3|qDtFc;Ene6dF
zrdfTcI}6kK_C$s`m$Sb&S$#KGHX^M%FD9NZeZ82a-OdqSl_
zx~2QhgzFwgOsnlSP$;{cR$ee;sYRK|fBLeD`W5gU#^3XTom+U(_mxoW9U(I;4~^w6
zPVL7_I>_Hc>=n$&oi!kK&8{t;q`$&3Z{W7(vOsaqWw`K6^0
zY~!mP^}I^GsP&7>iKsnqSv%rfe@{PA7+fBIlWCj${a~gSv*dqYM2RO>wq^9;W7v?B@LKR10
zkG($1^+A-49|N8D@=VmK3p~WEQe%cHgJN*SL>FAEEQ^K+M>#kL7!J0We-zJ}YD%iC
zyl|Ag4w*Vts_uNSywi_=;l>?>str&RE@n54!hOHt$D1{XQc`h{y~D8mn0M;u{g@t=
zBR_6i%{^L7h57O1@lboi>zUHf`oY6&Kjt0-N8ug=N8z1T7arVnF_S{s>2jYyxR8+d
zWua;P63dqR9@OPg&Ug^_f6DVIiXIfXkfJUZ^zUn_N?6wA2azwweE*u(N)|0Y#BSr%
zWksIDfS8$GpjhOgctUf{h&*^@jyyP*AOJKy>_M~bpYzIap~6M1MIQU3u!sK8g=VJN
zQj0k`S#@EB6ztV5Vm2wP=B=&6ho29)#|h)R4&+g24{FeC!*=gde-EGATr_lYk1Q#z
z^(Bbu9@OiiIyhm?8{b#h-{T%!c|S2*-lCkJE%#kk_ORW9zmpo-BmBZW+~ncg-cK*Q
z0zrLobw_=bxGu#G7OnHgxpKoaKUv>7_`v_P`zW8}(DbV~ddD^}3z+gf*3ySrS|)$#9i{nTJoe~xeKx@0Momv#@ABwbg(dfw0Mw9@v|SeqMNK#}5-@9)YtokHE2
z2Odf(vhawb#Y576TI;aohIq1E5|IbHf(tZ_JeEzI<@gQ-sg$D}D3otd$b)wetx*h>
zs;Tu8JUHuFbFq226uHU*ZM?>U;L7w`J2|14HMikY7dj#Be^#Uav{hp62zE|B$W7=x
z5B@rYJoawKb}`;Xm^a+FxR{%x{E#~_+jq#}@<=(}FMS_wj=31%rX!odplmT(9%N=@
z<;34S5AH&d$0DN)&lp_rWIw_9{^8r!_n`9;Wy|5uQ_*2t`n4@dk6mv@}f2AI`0dOw7EsH#sA|5oL56#Cu
zC+OKMC3&uNIFRG@airiE6-f?iz2>vSGmjY>(-vK%)7sxcK7Jo$PzUq}xATN4;+n
z9xg@@5AHPPtxtm_|uSqS!f=;myn&6_z0M>(p4
z;}6xr`xQ6}ml_;}%N35oQI6{1Y=i2c!=m&0;B(|ip)7nXCp3-!7Jcx#7@d%gy<7-N
z?uh1t()U8{0p^2}D?I9hD;kbMpYTW@e?{tdJYI&q;DxT%8Ia>bRJ!_5x+WZDO?6v`
z)uG{qwodtQwV4XdQTV1AMdg~7mWDw%qz_)Fd_o_7Y4D@+Q5sAZ@2CvihoUy<
zVDH&N`EA`yzbC)tM}Cnx-B$4w$v8Pap7||+HL8?6ph&H*%kiim6098-z`1_$t@I<&`4
zNgdoG;3%B(kb*yFqB?l%4&~rdgrjh>nUwUDKKyVwuMcw+rFEMQmlI*@NO}ktQ|rpM
z?7Uq#fm`8X`nn%3kI9D1?(uhJe~%OD@WX{W*?C9av$~=5H>D3ZjrrMB>AC0Ahaa|v
zBY$6?v5051dS4ynkZx^DIj>9GU%xF`pV(pJ3MqIW0V!C014R}dB2r&k7mJ7DEI0}~
z3-aIqG-i7+BT{hPoSdzMM#~#eM#xQZ=XYytakqhoAyV>N^wnq6x`YtSHVcZjVq4A
z&Vm$d0;FJPF%W#KF&7Qi4o7)zd?$-S*>beB9sPT-V39wJK49vaR*!g@PP9#IYnCNe2%|w4}Z|I4Mp~RylbCh>`q4y*XKf
zhwAl~I!qo)^i3O1+PW6UQUpAnmqnT>sb>~{vL_73RTU}t>mE#2+3cG4QeG0zOX^0g
zXE<5#)?C8?um1_N_=5U|bA51WO*+(9f^yG=6kK*cF%{pkf4t>}ljUZ4ezxb6#~ZSK
zFj9kL{hV`qsSjndbUvPT%lZ1aT)}M1eOn8$be)Q~-&?A@wA;EL&e(eh14r?o9H3A>
zIw1wS6Hq8i8W&n{UJfOQD}IAb9+eCWn5}a@&@itK&NVKG
z@4OtIAFDeA7QVMh;d{E$$RfZlf0lE>PcLZQ)kjxnSci{Y2*JB75>i6;i(v5l^i7fP
z6>+UueROwk@Rl${?7df4n0H`oZQ0y_-eit)LQ?L|T?|DUg*zRf
z*pKBN9L_~qL>~Jc;b>AQc>w3aH^D5}WYhXel)E^LGzteX+LkjG>SMnVW6nicXBzrg
zZAp2Je>klV<(8l2>*j13KAm{bl~)(Q_kycaEgA@-th(1V_c~&-SVBwjj!0$-%KZZh
zy5CeW;z>>5R*o_?s+=O6Q-}6UnPoP6DsORPHxVjVNKhSo(R4bj`8kXEWO^da@#)}6
zN!9J6dAo3*F+W?sl=|t#O=A3RLit8OS`kEwf8f%7UWdy|5hJY92Bfg<-3&qo66}2e
z@H*vYB?92F84-4n&4@hsuAvJ^OFa1IoC~f?`k>GM79OoV6<-~c7G2Omq)>{L3*?9t
zImNW*;iEc&V=~|ns)G-);9R&C&3nh+q7NncqNH3kkq38>E|@Iw;5Gwi%PAg#mH9Xg
ze_D_{=8}-3U
ze9Uw7z!&)kWAg#Uv*-SefB4Cgdw&=7nDoKD?0g;M9Uv>ynHoEh>3$viXEL3#-gZIe
z&Az~t<3ujh3*nb_RsDic&7uB$zQ{}Jf9U6ujfC#;xYCO7cVwd#(GkhtkOwb(kb?6f
zQP9GNB&D?^O3Lq5BacNg{1T54+9U^cck`0+k_07fN%GZt?Cy|)j}ej_6n{iP{R~M;
zQBFLlQ=QisU1i6~a`Q%p!_8z)!}t_
z(}T>YeT%Qt5jFC4def35rHz3LmUHLA3p7($-1*&>BG9nTIC1$~JiDLH)S
z4ey#&Nh9;lQG7{#Io?6Rt<3yv&-3_=<-SE9=l$aye)>;)@;mhLe};7Y$=QBi
zVU`ly9!F&JvUE_I2c-Q&Y!{G_C|G?1g}NR=*g85!h&+^Kq6^|jlCl{mDVDn(w{kdn
zfh>4uou)#sc+gb1?i#q6*;J|G_h^=u(tYE`REN^!<6QRV^K3~8YjRdnHS6JAR!Jk`
zVTT1R&V_ana4Ni#z^PFCf9Rr_Ua-pk?Z);n`szkg;r+$S9`>8LJm1YR?Y2pZTH5jT
zrVD0y$6VaN-_o|VZD`sp3On-Yu}L$Dnb!aPT-3@tuXR7&`laI?lKQ?*Uqn8lQ#TsE
z#gqBH1M=im*4o<j7(;4R5I-=$N9cuYG+?BR=-yLf;Ea`!of4?Gp9#sJ#;M0gB
zh4$h%4Z_+39{d(J@!*xSi>r#GD7^}f!g&Bk;ng>eLhm!q>qGhegXq3=Yl0-@*g;7j
z9|5enXbA77*3G#n&F-W=;9dneAbsGy;DFrKPdlK}#u^H*wrg<7G~kr0Jqo?yJFgF=
z{ClXYWEN@s?vITbf0^T;Qrj$Hm1LP>F~x8b*om~@J1CupopVrp9eg(JsQc~6Bg7EX
z3!|TiRPH`LksrBloD70LDk)dWr^tghYz!XBm&2e?mt4q$6DRWE
z3S!O0u$-C)kOwcy=6P_|LmpIevzsW81vgX3gZIFI2fqtMJh*5g4|+dtHf&E(5CR^B
zKHZk3kP{!he_y3C^`V^RLLStIAPZhHpm%U@InU!Ak?-?(N7VW}{HXOKpO>Y`gSVZH
z$ZvMka_csqJ>1CWWfVr--k6!5&b{4y;Pb}-Dw%o5NEFJfwJYF7Hi|dbwU&`@`e*mWa22?s8Q%4?r6$0hRCI3i$
zSq`=90v^wo*`qlZ!j&aCxH6+oIoBW$?sleis(jP|IaWylThD00l!J2(>XbXk59m}$
zSU4A6N#b0%3k4p0#hvEDVUBYVbu9g1K-;$8qEp_?ML*?@pXr!WHl}edd}@5!*H6Tp
zpT~HGfA2zDBM<&kk8Vn@y_s%mP^gEGQKi&D&Swwa(nF1Mw=)@WXXI$HE$x+qJ$Kkm
zZ77W!YLvHsP@}w*wgS!$mUlR*uTnqoC5@(TI74Gv+x5?=)Qz_%x;jku%i6Be-L%HX
zS}AS)iIV$e$s8UJHLNQrm>U3q;^n6cF7gI&2k-Y_42_-gb>A-BxX;VMJpl5c*FG-j!;f3uhyqmLZQaTehEGwz+uc*QlHl254g4Q5e#5u7deL|+QDoR%#-uE?Y<+|uK9%DQ(e71W0qHb~G)
ze?_O?h1MzJEu_8$G*|!!;Zy=Va?7!doA1Cw`Mzq663tPh8g1+GRUjILmR~@i+;E3M
zqO5yT7|t=GpzQX*gWtlP=W#-%-ei@ue+h#qY148r!zya@St_+Xj0-RDDAJ;C+W?V(
z>fn5Zqwu#GR61`+e?aM+(B?f%)SU*zmJ30P4~ost)is;alj9(~p@t5|UrW)!Ht8^e
zo#f7Fw5+kT7jrq|V9L<})CaE(kOwc$aV}gk(Ze3k@2z&gIbnTxzuL6i37szDf2A+<
zaVlte{fV>XHSD~P(PMEm3U4pZM=d?#Ge28DYWWPze2U^_-=wRnJ5^tb3cWSSQtQ1X
zPh9}IIxj)8h^VkQ87KO5d)Lfs&QTFIF
zBvosbHHtx@7WZkhh(gmcehjs+f2j?9sGrth6q#uZVfYw2Ykbcv8l3;E1NojVW?5zB
z9`Bo}77tpWpuRZsVz6*0j8mc2E)wv987LSukj%zQdutRs;JCQY+wHf`aaf<&ruWRI
zWOp2d-f+DTOng8Mt-_mF)0x!C0Hi=$zfT$f@`W7dOm-@h+WDGGJC0H$n!4r8P=D}k
z?+>!u-cHxNEr?2|hx2Bm$yfTJ&9%eHEi5NY%Nqjo4#XRW--t_Hx<1nS^DSQ1Aval1
zw_K|%Ssqd&{3OLE;&rlm5doKL$mSXV#h@C31Y9pIF;<--$y}waUXGx#^P?zB1Q=xv
zIk)jd4vxZ!!~(}QDHwK`l)h~DH-F;}3ssbZ>#`-Ot#xkeEUnwdKd*Ij6n??ooUQHZ
zlzMTpd>jjv&iUQ&^k77bnpjbeeIM0y!V>dV#?vA5PUu9wC_PK#W~<{wSf9u)d|V54
zzVF8ULZf~fS_Ist5CtW)ks{y1CcWcg5@-QT6y!|v
zg7Qbhk>Zef4T^)?+mcdXkbk6<21rusg-E5LEqOpbrctQ#oA(Z1QcfkL(~#;q1lI*G
z#mqV{C}Tj7O67;kc5H;_ozM>_H&`U!$wZGkA=87LFf3&}azR=Ik-lo~9|R9uZS!Vx
zeUv(Blp=*K*WnY#yt-U5{TxVJDWt6;F+E_a6A!A}7m}bWM}ONR4S!#)FNC9?;@r7#
z?I%f>TI%rflc^6=j<=MJ{Y!u$$6<|RmknR1)W=-UakfsMIqnyQ_JmITBv-Rc_lPl9
zzY-FA>2D$AKG)QD!6e;mDW<&E_;^O74v?bdmSw@SjENNioGBe40k;@P!F~@E{Oy1!
zXn{x+^l5-5D@&$rcz@y5n(V>vr`5skhJtcnOj4}3MJmX;9j|oE%W;B&6X*O~=wx;Ej@{O|<+lKA_R^)X_
zsoy~Y+GDujrd$J%f~yrwHuXDt=agaYKA{i
z7#tteHuWgK)Q%W*yt6v5^b_&q8_rYMiGcG9`uRj?d#N9BZbG)y@nt<#M^M=MXZZW)
zGYRqATIwbhNlMZ~>;aT?yDe=Me;b;`=mz;b1@A#2d@6Yp60ij<7AA$_2^OfjsSYJO
z&r@)AMzbtUGJpIr!Q)h@&L9Ct&NRhaH0ttD^1DT0Yt-lA`;TuGaq%S_(GRUquodRU
z$~%
*IQcI>Exv4Ck}lCw0Em{=xW&4|~A}
z-B9pd(AUO0DEth?i;L?y%%IR-h=XSuedw*Xjd|xNi#DCJy!C6Vnx2=DQ
z3uVoZcUh_nJX--cG?pamyCWvZNTHoy9VymZSc^DpJQm)OM?dKxVLqVzqo{nNgC;{!
z4s%I4w!Y!L7%VgdPNd+6Sx5|aS{EmpC^+fO^WdBYbK$4P=jGsoLl^Y%+Hgzn&wCgj
zj{K4synn@Hq2=3?eS_hY{EQTo;^!^PO_zl_Y($)zz^A98d=;MdL1_(B49Amt$b;~L
zrrqbnTxn%69rg^bgggkuk>uo^R2+nWesBwA{UbBYt6{ai`tU>I_a_2lh6;}>Id!<1rCjy452w9gnEXytDFr?5xi;+Sa
z_fm!vcSykxgpU-lBgThyOspf1?Us7R%BJmQlz#YdIFETDcko8WG!IXwr^4FLVegoW
z8-Gr|8~d5z4Ke$M^m8Ka-9&jJ$@!_&_a9f%HzZ2sGX^FPdZm1x2ftB|fzOZ_^&3H)
z`}1MavWK4sRa33ET;j5{T-L{<7S)1_{o@H9Vki-?Ga)7E6mr2@95|K*dB#9I$RiX^7OQ}W!&hvo?d*r*gAbj@_0jTctZkG&e41&_J4Bk
z4Yz(Ftb-<{TJAt_LcHev5G;a0vIP*7lB3H7r8Gw|3Am#jD9Y`5ZP^^e7gR|DE*$Or
zR5(6iDp~So9NleNP;Qq#VK&}wnYM~Qiqgrf>)5=zXr(4aG3DU&4{{V)45#t>%bej(nT6X=
z6l05F8K9e{#WIU)IEwaO(&{^?rRn)#bMB
z4`A#^V_gtPqTqV!g07l8$}dKV9lvwF5^&Z1cxc*%Hcp=#4`Lxbe
zbK1pDWJ1^Yj9+_U6o1~`#8Glu*71<{^|_cD<41N(H?@uLq|?Le2h*YdeI1cZb-O>-
zR{lzByVY9y9h+7kh|1N)?yU%(Nl@d&gG=7#n$ZPvG$pkL;kN-VOqSkCx-i+t8wAE6
zr{?Mo@_dbJGkI{ubD?O72WMxTEIm?0v*pm4Hp>aqzGJqWT7Pl2T#H;_)OnpMx9=uL
ziIRAv2-ps!t%zOd=Xq_Mu)VL1#{*S|eVu9B^W#$U>Tss=ru@UT+r^Q+K3U7OYO{m0
zrDvb!o%F6@U};^&-+$I5ODplqAE_I{$fhmU^`QTeAcz};9ENrjS+HtG3N3dd4{myi2fegO
zJm|I=&gJ=LZ!zaWV2~sSB{Y(QLL7PABS+@58Iqpi!KvngPJO%32N%;W=3Mw=3EG7&
zjFMR>k(lZ*40lSUyJ-!IQ!4e{gAnma9WqE?cx783Zhy#e8$xnWgFthkch-pq)`$yU
zO1%#9=$G$8KIV5)2Q6lafs;p>Ce+3iaqf~9JddsJ8Hb@HvB
z#JSh1$A9VQ=3@Pg6snb33X0yEzjy)t^ObfFEns2{Yjo}*tjP-=+W5=hAXSMU|p-tgX{9!wXWS0q-dcP
zB59plKh)TszwLP9nmKQ|9SX#pi_&wUUFsIs;eR;|;=#S2Wfw#1+RFVSJo1{C+B?XL
zfvw#qo;93r{4+V~gjrly_j&M>8)n!VlB(_DyUp?kwyT;flva4&wC8J_tuU{~F_0Cp
z*gJCD_PpXp9x0b)cvJAK4}M{{BlYd)JXc!dEv{{;Ygo3&(c~6gj|B>)1RxwBI}`HY
zbbkdr`1lC&Sj!sf@wJhWK>JWeVi39t@!)iYJfdXx;XoZuMQN^4>6~hClzj`zn2Wkk
zF*(Z8Un>)>9svN6hb%@~8_uVaqy!O2lG4jgB`*u9mGi!Zy8)06jZ^(4d@@@0z
zg|e3hTT*V!<}<$=`&&r&mrkurcgFhGT)wYaS!AZAqY-b?t@PivEysSl>LE7b(SLFQ
zM~Z#_EbA0%H;BNLTn6Ul@F}=(S|F|_4_WQj&7O+dR#^a*=0=UD?xPggp`d9!f!o`04!_H)p^$+V+|>ncV|liV=8%l-b+amp@kaS0{L
zdTZzqJ^bc*AM|J9=MKM+x&g-@OHUf*%Nms~NLj=~7jb~UQmMBsZvCCkSRN01Hy6hh
zDfa05gMT^C<8Pv1UqMV-dUc{%XrI;M@oZ(H#Y5RmL>`Nbl8@(hkOiMHL4O=h!e1wK
za1ub%@>0`9>Leb#q_NCm8_pdi_6w-n@R`X4VWm;HNVt&2rd^uFhWBOuc{J)w1OcRr8`
zeTp2V!1MzzW|mSmvoIRget+0?*aDc>hi}`zB^plfcta68QKG2r`ci1*MO1aap_rr<
zDIHu!@^?wGJ;X++U)C>BJgY#qD7fPUg6+OPJA4ZSiV*JLx5tEqg0c~jf_qY&3hxo%
zRCsS}dM?W*>yaUEPESS2@Q6Oi2ZT@zdW7YI!XAq++q6oaz)`p(ntwM7A9z8t#9UL~
zowr?Vf0Dp-vAyzUI~Js`Z~Xvh(kav>9@>So^R$1Qh)Lc`c0s;LV_QvYJ>qiAhs$#m
zb1u&ZDNKziof%p`OFy$5(Nak3O$?ick5woscg(1D?uald_ZNpuicQ?OEW=^c`FN2g
zA^vQOU$;JY|2PWw<$p%J6fgECD!9EKV&ok-GqEP$|M`1&wq?{M$C1qDY
zNqMtypwO@vNDi*OhSm=&qiG&Wr-5^!XNhPOdRy~?fB2H}mVY%$T4iOu4yj_^KWOW7
zda_DBLP5EF1I6=eatY02iAl$QG{n6N(c-F9*+T#dW#`7q4YnLkZqT43Kaqh>xCbvZCo|Bg
zZ0cr@zbCLgUssPA<0!mxoQ{SQxsrC<=X0f7Iy?u`o~u}@L+dBNjdEy4=zW@gB=bi
zlz$ccq#OqbP!Kw
zH9eX`Q=#f^s>7ymB%!Id)RiFh--6{v`IL07;4iApJ2td
zi#c(vcCqsrXbY_!k35jl_v_uZ(&AAnmba(;itiHal4T)QlWcmXFmYU@6A7u6PcrOagiLnGI#MiO$t*EM(HK3
zbIt?9?HeY^)L4gK4r66Teeg1J+O{X0kgGf5%yKh~b9tnJV9r+Al$*ATk`)cK$!O?Q
zt!HRk{&tA^Aarfg2lWv1rgbTl1uv?D_sr*$1FsESFv$x}$bTKvywV>5WxJT4=uTvH
z0&t&CLCPnBX;)`Z)OouoYsbm*qOSXDHC1Y>!&jI~q2z(-
zxhVU@sMA#fe1Axbz+#oZzYb7i@Ar-~UaAxUq#*DPQW(K=k-~Tv0eKLz@r8)vgf;IT
zt`jOLZQg)b?spvwI|2c1Z#jDJK&EC}1tOdp0yHI(V?O94w|
z#TDw$LxCyv8>w5$d0pD$NxBUI<2d_T;CDr*(B39Q9!e!aIXLT~9JEF2A`Y69wFlv}
z2t4>9FXZuj$SawWav`zsNKLaSDb2--Bz^Gay9=>5Z+xX(xhS8rc2QG_HCt0--fYzs
z8jeEuFMnwiPQIX1ex}_DsMXe~AJ}x91)Y!p{D}p<@I!)fSA~MoJHTcJvR&OD#SeY%
zMCS9H*A-
zlj24H(H_~n;2yl1K=+`>g2{3{l*al#Ia_r{8t1~P6|GN?9+IRniZc&~=O(qTUQ--hq*#`;)ZXz_YZS%;CEIeRIw`4!i+W+Uk5W@jPnNQ}taaI9oF{OOZOUm0w$bUJ)+dc3VCn=N3c)6f(bm7FTfy&vhKC?VNjgP4o;>!@P9Bz
z$y@0FKF8baRxUK+rWm{fIPE|ug26Z3uT}bvU^rDwc`Nc{
z8vC|>M%E#sQ@9GAHdPRHMwD;9@
z?z{BA2nIJvy(Pi8q5l>~=C;a`G?_fODt^L~oxqMaiKZ6YaCII%ncA7jtbgts|CD?1
znrq%o-TI>(ngjqqV@@;yHa-6#TVOMW_PlEtZ+@M05rutQoDhE)zPk2H&ivaL2b7Ztp~
z8pV+GgjGt<-s49N7g2InrGKZ~_zkw_`akc*&SYhjB#qnfe)0tmd&gY<7KNP;7K1`O
ztS$X&ORT0F`#coz6c^QW*7&I}S5qF(cw3{G0jD&>sE;j%Wn7U>CZJPp%0N?=f7p)f
zkHWpoteZZ+?VE1eaR!%dXm+kg+y~&4zNC61jw!ZXRcRr^A8
z8ND}fDaV~qts}}Qe>;)*y%;-X!8goDYrikW@;mA-e*<`sd|T(z@28g*f6z};?ASZ`
z@oSem_@+YeMx~}932PZB_Ul_qGf1+;QVFe$g%QV5Y;}Rem48+r3C
zRkv@9+p^H!`IRUmEDnycDmEN$YuY*Cg|Kf)8z)?0O|dP*)uSs$e2=*8jxS=nf4(iT
zmHXQjmb09TGk>LH_<(O`(SPgLtqxF#ze@fuWH}K;OaHcX{FHV{bxRWOzTY(Yh3e=|
z25T;MD2({KXnhK6mk`$#$}NXA#FHiDcga%e4d+9V?r}dTzO0q%f{1#}GU^REn~m9`7wa!`V_
z{LL0f)>{+5NCTb~k3+{oh=hUp=iNm^0iT|;9=51{+GS?V(G@as%T(f0UW%e6JS0(y
z`MmnICTBJj53u!?EE-Z5sb!PgwhlmCSA6?~>U>u?qg-Xf54k8ir50}N3;lb!rQK2s
zvI?_xdw(t~({T~DGvO$5geaQ?7c|GK02jmfN<6!wv|G$IwtB?Dyb?rb97^d9KTwI@
zMZ&0?Y$L*Ve2{!ums3JN)KwqHrBhh@@!uZ_+Ag<+uzIkki*5-d+KmMw(JBrM)W4ms
zJ>yveosdPl5J<9U?&kt!0!g~{RzE^j0+O{KiGLcd4U)M0+#jn>>6Ki#RUoOirAmfa
zj#meHAiD8wlUU-48Qaz&LJ#M3h=frY*f9xEjcr5-;8$kkI)hqgfr3vSk)faG%B@J#
zzjwJP!bHLs&{!>WD}connmm4ga6y2{|M~GQ7Dd@f9!OSMSO*E|mkCH~KTwwa$wJFS
z+JCEtcDUO5KjFE275IYBX*C~#XQk-dKV@WW9z*G+TqxgpNj3N__Q)kmaJK!~NmEyvyyPAp(hBYX*{b+twt`JEkDy
zqNLqf5W0pk;ei!7G
zd|?ml`og15_=4wj8vZqrJl|olWH~&d7Ou;pw`{W{+5=`h>Bnoo*_u(7#-p$Lju)Sq
zM3P62eGkka(eB3B&R~&f=p|QzL~CC#NO<8rAN1ZZ`bNTqyrArF0T->=8c3=Pw12d9
zI8;8LUG!DWd`8uF5-u1~`9eLPqqLHbfwqJE+p0Sut}9Etg?^X0n-K^4w$35R_?r}m
zdmN_F?>)-gq-TAeJV97N#o`u8F-f%GcToX}MET%X##mXT$4lFI(CdOIoI>WAJdvSo
zttiQCpeVOpv&o*1w$Q?HlzrafCx1LkD=!qMR0;&tRJL&jl6>2?MY{0gU7brtRB|FB
z|Lz-8FAS~p1TJW%tru&cQosCxQjv86I-z$6}#tBJJq#iezwY{%d
zjp%LaK(gL;NtPEA@Bq8}2EVQCh_@ourTvtt6qc^%CKYmv>n^v8nMWj=mR;OMAkj}(
zAs4oli!{h2N*8%STIEEOC9TCpZ|V>yRM{B;S+pxMF4io|qHM7v3HM2sw^&T1)+~#V
z^c!Qp(j48tyuwARk@JzoCyi>lC&uc$-UtILMwj6YmF+JIi(Wn_#NazteyxVi_5olCx1(S41E1L68@Js
zxb&=B%~JiA=ONMon=ajIC}A6*V1Yyr6(mubeG606BKh+~eS9H{Qe#~JR3>2ut&q(6
zJ+JXoo{LO7J&b-))rBZz5+zZ6MEhE(AlR<53yC&-+=D#s(12dnEafr1RKL;?SWapA
zt(ZBBc7IME{X)4~s(-?PnuOW7Pdq;%t;0dZDVHTDz4AoB%Q(_AB>1d(At#*hZKWci
zZr9kBZQP%np0idQF<5uhgBp95)Ak==nirz@
zcJ^HEoIm-Ft1TEW)Ly^dDQn8JX8yMJC5rCXT7Rz5i%UssEPq~p2i4!@_Aob@BGr)A
zF;WoV?@{@ZMWU}9kqbd)F&B1E7k!V$VJP}}F1iOd2sh#*s_?8V=d9@4NS1SJ3&$65fF$(iG9_xqd$BIYD$AmDp)T$y4@8il)2=JNu!*~@
zLv}WTHM}74jDM(~3%cNQC$yp*D4FGQB8pDQ!suXEqDX7;=ktPgX7)k?@;p~Ri!+~5
zmHyquA-$*FaFBXFU&~ziY&|iPOlvRj7EPY
zjQDAsEB^9?@V7!*zil?JUF{?jTAChln2*N8ON+KF+85IKOs7-=7ctm2(Tk{8wzyFL
zN>^aH?thEK;K`!jfw8jgiEQ|Si=U*F%TDwDuQieAx%>95trmVo{8An9Q>s}czOA>2
z9eQ+mJW-RDa4RXwa;l%9{SiAXZXyakir+}0-fKcGN;r`$T>gj)uPl%YKlfnEvYIne
zB96H5$^vDHEw*$Aur^3Eq;?wUd4npN;wrZvNPk?81X3MWNFw2qLj!U6lcl*>n6kVW
zh&u<gCo4MY!o(Ttk3h+#Y|Lwebtvp%8SUzqV+9PyS9iA&X*#^I_&xr5$3
zyyUuEiZA%(^O+Yb)oHKmu)A=JV(P`9$hR0b=|zWgtyjk3+IAs#D<2Su3q9w3MgJ!f
zkbj#xPB`5O=hp76Qg-i0EtdhYfHcP!7Q&apr>1^k_~-RWaXm@&DF!jdQg^u#@Qho<
zkwkmk7fAGQ0xnu=bphva9Iem6aU=&V^?qK(8b@n|a6*b_p<3u+YMzji@lag8eLO9T
zo*{
zLp$xde<0P9Z6dAx_>5lm-r@xprJqGJ?!EO3Gb*jn1%FYR__Kxa+uDL&T0?qtFn@8Y
z;um>oi=?FIhJuf-7NYN=5p>8OE0BwpqJ1nxelinGM)B_;*~V!27tbdnsq6?(?>8|WiN78+*1{xRY}ZkT$otcz1BcU0=7tA
z7~UtaA|SuWKqoS@WyUWn`Xb1-^*!mkQ^u{x$}DepFs%I*xPb6>d!j$68Gm8ib_G>~
zMBleTkVt;^1_#p01CdQ5Ff^Z5t!kf?c##Cst1(40i*N
zWQ}oMmg8%hw(5dGa!7g>7rl=`66KSU4`Da$QisPnnn9vHG)a0<@XBdM$|4|eT?{CZotp0fmffKi%
z#BaGv>2ergeLw
z7`>=^U)C0s<(B)ara$RPQl@m;_nGxI?L%&W3fsgT?r$O=bcmpoA%BM0k^Q1cqQLkD
z1Y%Q?q7Ou9l58f4awsJLY?-va`EC}Ai(d0lFUk_dhq;
z%XWWO0;Bc&3VJti?hC%Y&hJTN?d~OsWS6^)3X(oiI@STMJ%9KJTS&s!L_^VC>sW~q
zrb{VGfP}YfWnJrVJ{G4H*$B|T^dJ`{iz64MqeTPojUMv4y%yLudTxJMBr?$2
zkxa*KX*TRdyGx0W)7Y-R*vUg{s300clG0UQhZ-XvkO7n}%?k!lc5CKkk;SI?|A3r@
zBx_ki9Z(o3>YJ=6prgXH%V$JzK8!wpdOlxR^B3{9`+wE0zhZ0s$ZP4y;q7_OZ#6>u
z;1etOgpw7JizqN`3
z_4=Gom`GnY`yg{Tkk*$0iL%Lu6H+qBytim2e_E71J-kqgr#C&vgCi(F0W^@(&$^JE
zG#@^Y?|)<2(R`GT`AEVj4Hq0v%d0**AIc&tl9EQO+38R{;gshrq)4(k6D5heuwXQ0WpVB!4wB<;A+zNmfF56fLNjaEp@3xnzMd=+{pQ82KOMf9zHnS6g@Y`zWRMl}*f;=5
zlz)D#ZsTchE95OM;a1jly}#LExl6jW)MDGlgMv6|e(!WTL#fC#<#Dt=x1
zqgeOEF76=3t-qTb`bX=B+&)HAKmaD?%zrX+Q5q{CdFED&OOSS6`pIno7#c~kI7O=~
z>}vYRY&aurR};lmKC1d?jMD?@sR&7wtq&woJQPWkRde(6459(nHHSFlGeBL_J|v?k
z%AWlNQIykI$VJ(0_faaRm8-7_P!^@Q`QVZ0m#weJ*5U1IG{Ck>qm#xC-G{v+*MC*k
z`96>zPFmSMy`WsrM0$&Y5v_d5M+ih|sS50I&1X8T7lky&sO^-+l!XgA6u*v5%2ifk
zmYuCRpC|zOs^mh7dgQ|J!d<)>=?gmKqFgumMl`^F*gn(}ezu-si+%U!0SCCt?ZYbt
zz$@yE0x(#QvO1|j3
zEVeA#ipy}9Xbq(A!}_3(vnEm@K#7U6;x%q)(TfUMY*X%@2%(kUw=V^nGI-Z
zoTdEbb^P9R!H5JjPg+qng?vQ2>1=r_6}J9PLjOGTB3e(#;>}23W_=?@U$6L09I4#9
z(EEt_QS$efks^YN?dOI5_$QJm#u`ba5ciFO@#tvufwm%o9ufQ_G2x@g$VA!LBobbo
z&_E1hfof4&QXJ^Hua-N@FxucI*dq_2rZCUsCgd*ZN
zDMcSK)q=G1<2w(;MG0;oQ5xh3kVeWzIC4?$`XLFQ$BKZiB!yV~Z^*46NvQp^>2NKR
z=EJKg97jv@IAuNoNPhzHM>`zb-oTb7%sKCR;wbs|aI*0!;8
zKwvO>K9DBFIEtiyKGuaK%B^Hlu9D_`1imYV_B={RHu2o9mKT%eEqQZB)}-}1i<91$
zTrX6|ht82itJaowS6f`hMVJkGy%q9dFfzq@OB)PJ9b%F^wh`6Me8C5uT;
zPIdTxej)uUU(iTeeie;0r!3cINzS>Y7mFtno!-QMDRf>fYkg$XSWAb(sQ2S}`!he$
zM`((GiXfQ7Rj?ewz#Of#t~4F#686j<8lhbRj%LSK*~`PBvHOY5Q*Y`Flprw
zFlm=>F(HYTet&^PYjY6Miy=>@XLQ0w>J=p#NLisrBSly<()S#Es6R4*QY}ziC5xkt
z_?svVq^zeF}y)_xT5^^EvS1P!G0OE@F`%yY$z%96{Pw8cez+>Zw0
zLy0tylId|Ee%z7_>VFp@rjoWZtjekT@sv8cppQg3+ukW2A7iD(|
z#m!ONU%1%$UF+jf|9m%TIv6|hd~A8V)l9yiEZDvgJFoQp!FIoNo}-MIr0eWQPGWrS
z?IS1>f5Ac$t>z%k5JMzK#2Hd{%K?5#d8~@KJc9D*6fVpYQYg^%(-YD?3rUocg8*gx
zlsa+Q1AoERzqW|CNNl}mi%1`R?8>D1`z*<#e2TC<%9cd|0KZ~JB@-c{7hi}%?|$}(
zwk%taRquGPoMJ>-BKoM06eG$8>1@2JU@%CMqWIJJQauf%if4>#Pu9guYa@JJvAfc3?D@RpMM-6Gb)=g$c2G)$afhhdnC!X?r+Pc
zaQU)w-32L2sy{7cti5>U+~N|Yjef|D=stXXFd`P0Vsb>ri}@7;C4u;0UdTn+SM!ne
zaGd8+NH6
zmVaibN(v8UjhX%AV9n@-GbsB27yx`wkhm!O;HVd6KMic8fUjJT<#{>AGNTIlL#k6&
z%{Xag$uaMk3UK|E%%FUqoXxm;gNQhnd^``iC^s@`M#|;~2JE^mcP76o5N>KdG3gTl
zn{P`^v906!=R#`UT?+l>%L;x89rGis;eXzGOEcP`?(h46^|?UzvBe_7Y%y>wa#5ZT
zKrZ}3Nd%B5trnmtbbC1hBvL=RwB(n6ReM&a%S8|i;93}q`cEX
z11W0ufwo8%3x0tD;kJlXOKiYOOFF_-SXkDkV&gi
z6gQq~hB^0WTHE-A*GE)CS+oQBs9fdz*t~L|?qZ3nK3j_eY2^Y3Qoc6BK!4g7Gz%E7
zngeN9QqTa(CM=>7Dwz{O*H-Z<=y`plxhok+Y5#mgHF7Dwvq_d#PP3tQm2>^`+SQ8j
zY&JX*^(Qj(pKJf$FOV(G@E
z7ROxT^N^p9(G&r%==u37nSUhOArK3Uu9_qYj1vx|Eu3e$=zRfTdTHB*x^zOW7hpuTws3qqs_w=)(b$Gbq8011Wjtf~8Mrmml9D
zE=r9=ODl)`5J*uuplVFIS&pe6?U)r;-wOB97Li2j%zO+=8b~Q!vwxP>8~NE8>DySd
z$wPf~`AVpHzN>G`UkEkd0NTP2BYiGQL*B#NT!ZTf(X({t8yuaCe-Zpo-$(@v*11?VpP3O5a;WGWC>
zJA^%JBmG-F45WYn%*Ud(L2AMHS-zmYH#r?>r#wn4188XFz7P7PlH@IMUywMt@SJkJ
zK$Dgrq`E@A_mLK-)uQ4%&Cl6Q?)n)AA1@HOESt>vMTO8mA%Dy(?nEkdb*nvP7{3t;
zu*6%CT6a&{qQZ5#r3XYn-QVVkK0+LlXr7l9sO*
zAmW{Jxyw>qwz()xx(|wnQ`Qm)a#7Nq57vkSg|KV^PLv}TW%-ViR$v4Wp-gEH5q?s+
zTx8A2(u;mZ-hWb^ZALAd&2q`Nu&rY1p_n5R&dxqo9Whbz9nMGDu|Om=26
z8=V&dO*tAhsX6_a{e^M-AS!Ly;If1demFgjzS^_A-=2_m%?~xF+=Q^yZ`TmF%;*R)
z$VJ)TLV(2m1rIyU+)1CvEjn({*-3v{+tBd@Q%7yNqklX|ft7l`iLPae`y)UJ;k>xw
zBt(-c{gC94<)^u5p`8GqOgXTh5EG1-<}?o5YE~4r_7M}#&nWHlVso|!;?4rODCZkc
z7G+V5Bu}_kOD$@DGCgT^u}coMpZZU>>iB-NX~h?Y=BDh*q`evW9;GHCR7;g&UdQhO
z#D&{7GJiBT`Q$GA-ErD4PsIFTsbKk~S}Z=IT<+SYb*o}|o9XApK%TUxGL&Z^6$bN)
zYoAaHw};b-Ln1Net&9)C&_FvR#S3Tg20GL;d!W$w4acon)BSDLDX^*OD-j8kxRZ2a
zeZNS=DWyL-UhbX4{*q8a=p}YgADk%x#FH4a=6|lAVslX(6yT=#T|lJJV#R>9YPPh>
znHb`tq(ngB+LslF&yyGCtnBVhbJ0WF0;;oAmwFezEU)yUeq~1OKA&N?NYn>&P!_F)
zoSn4#zK^LF+W=Z2He5Hpa0aEzL8P?(ttZwhvH5kloEkK`d2WgCV}#CK~*Eqq*mTDd3URXrDlT$C+v${l=()zr%iGoDBRQ|V
zYO*N{Z{5wBNIgxCBnnjKY(T5`PtiikodBGWa*7kVD9@yt6S6I&poUMU4u8FG!$g;b
zW*a+EA2gDph(Dl{cKNuX%SMq)vQ;Qjv&()oF9B$h
za*{0}u$5A*fs3~IKzuU3*@**bD^Uvu${I+|5GbzF!GXBSxm*-i`w-*9WBq^Zy-Tky
z$#R}I-@jrF?EyFQei@3UM1RA8VVFGdMDRdp9ZD3sDT@>>*uUQKL_}6))mLk;y>>Un
z1GsTOVs-5NsA}UkzvdNj`)S4TNT+%4t&D8%{p#1u(4U@3Fn`-Dem11i#cC=th7$|M`+a9SdKl$PlSg0nvk!Sb?(ei#5BBq&
z1nS}b_D&`Bkb?Fe!&86vMGvK;KK(_{ccvc_%;h`t*h3KJ{>=|+=st-LcKUn#W$!&j
z?3b9L6HdMRln2lIvrWEZ{GSb_`%_lP{@sT>7(pL%!+mxRi+@az3~RsCt$TUMgRk{5
zGueHD9*99-*6q5be}8D7tgZXd#yQ(n#oy=W!63LV%GX@Dj~VDby=vUTJse8+p?d}Jy?pe8WpV$aS90wa&G5=7`=Ws!LVEYUJOrlR>+bub
z`xG#`)6YFv1)m+=X9Mxlqq?C^YmVU-Z#7dcNtr
z&(QlUWVz`#>cp15e^Y6}v`bwVt?%=++7QRaxez=XkZ*@PohWDX+umJxg3H6xe
zFn@B=z+Nh3Q8rnjIk#1~C*R3~9yG+AJm_Jw+@~)EX(4v2D7;>;PVajZY3~o*!;$|n7u9DU@@y1*%+U4n
zl!r^FJE2vLV_q3XUKYZGO8l5x?*5b~mw)#MhW|l5+~zq9^lW`)R2)sP_6r0^fB?Z=
zgIjP2?he6Sg1aof$RrTlb#ZqM?h@RCy98NW7x&A(-{0@woUT4SGc|pxYf7HykQZ{3
zfN{TsFEmIZEC2_!!%5N^Ms#qfmoqln89~HF8j$>wX=;Sz1;d6
zXRVG2w5>9UU5u)e8MNIdF8U~ueIPJ?fDvW;t1xwSMMVeRbrcI0zxOj@64B^-xWd8p
z4$m-L#p!qjrpU9jHaI~W7h-15jWEgBo^kc2j=w=1mu(-Ou^ji!fq|TW6rr&mW!)<#siDnmCM7Tcf+O&2)Ga
z4%|}rl2UYJu|xLTGuzNjx5)YzScv@WQL()suDZ1H%-OG_4JT{+iziAea_plO-GzY1
zl~LGdEqa(cwl`Uu<7M*VOW=$ituXyJH!VWKOX6l!VE)8EwO(dKP>fxjJ5_?^Isj
zLog5x9x80c-V!oL1YYtXfr(wY2V|#@BD>PQ2cr^KXJ%N_5LyhaKYO0f*v`8uaU)eD
z;R-@tNLnK&cBPKQgS+k*MH(xR6U&-a=|~*FKz+FEAxhbvGJe
zr1{qNpS`A@i+ju$z(Aw?m
zwb3p{(|p33y4v^V^1!+%)ychW`WMe}KHW2&)-?;Gkq{^U!b@)(7!iukQe}0`7C9d^
znY=JWL=N11j$2Z!NeRhs_APv=WWH*}S|1tQ5rrzJU4eN)4xpma)AA>?)_h2(U>HhE4e&V(gNaU47Jzpi^#M#1|{9twqd_T-;Opu*cthAPQ}~S;%hq1
zFsbh2xy2X4gzHq?F30yKZ$-vPn1<6M8yS(0Sz*z*D!jCTBD+z;bv*;=i^?Ndn#NSA`muZms{auOh#qs`B1#z
zQ~KuIDBZvQ!~mBtU2`;7|9FBFx5uM*_)Z=5d$n9B7JhK)M|`^Hc)_6SAG}hzWdr}~
z@k#AR#@`c$!*Gdtj39&m_h$b?yE|p_Zj2HJ0aHr0^ovdKW?lJX&FaRR0nI4=%X(H^
zkVd}v-$w2a;iuh7tRWwi{d?1Y9bx`XEP}LmzO8M4Lsl^8cVnb)+mc5^B2EQ_A1;<3
z*MwmH+D2iyMxFPY;FG!sSdF93(h&4mBN?d8*IC_v{<8Zn+3fau7I?J4#=_O!{A=6^
zJO#Bsbot#!h_pZ7cZE4!o?rlvgOvom=bs78!zYC>o3>j
zZ%Kx(F<4UD=czSF)5(@omy?F)cilCWYQvd4%5qkFWlKnD;zuYk+WwdIVylsDU#o$Y_wc7t(X=_*(Q0Q>R2P75K
z3o=p9B9$q)qfnQgrHaolq)JC##`!on0*nBe;?A_CEo7f?WaN}hj2WmD{nS!;HVC(BI%Pf+Bku4)77jrL$q%D8$NX1
z(@q%}Ln9BUcNkm`1f0$X^!d>h_-yis4?MfkY$B(8mC^wG+&5p?0_V#*m^0Xj$Bsn;crN9`lmmnUPX
ztGXwV$Fg(}ynD#)wYJ^r=Wqeg$pNh~d7F@F*vmGA4JeoCyKp*PGc@d)?TlWIeD}HQ
z?7?qJf;N<<{Eif1`n7lF+uHKewQQ<~-F&BMV*e
zN?ex{a|F|#!s3+g&$h4Fc$FQ_Zi)vEo}e<@SD9zck=R^=QQcSWDLZ@_jGfVHQ6YD?
zVh{L)pIsk(N=!A5pnJW}?e+IZUEsL2tLl^)1(D!1>?G5pmXGr;q`Sv+n=8#P3O;%)jfH%3R0Rb>vT5n0h+Ob3QVNn
z2o7jFKgTg5uf`b7_sxWT7N%2tG;TLS?bvKQ;Pb0AW*x*k2t$kH6VRqRMP2P(y_@Cr
z#>XcGaSZPrn_jr;DIsL*v?lsQZ8;9x+%(_^aGyn3gflhh?C6_+R*2ckjPQzm_v^+t
zB1NTc;n!CB1yhy&U#e5%korlAfw3pU2w7loVVz$>YO%A3-gvg!A^#4F=q%)fW2lf-
zdiqxu&qdDEqhGmyN#koa_RL+{_rAF5`X&e%o7^lwr;`ryt$;16RWFR`*=a~>ttpER
zHZqdVo!HiX$Ra;}DjZzq|5EuZ0JAhO*uGJ%^jG2^?cPAjbr9cg6_2Oto
z@mDyX0uU+ZJsG=7zZ&U64de{cd0k@GHLiJBz*5}Zxm?6E+#L?u!nV-stOdx4h{bmw
z?C9f6W-ST}bF=x>x~g`#*j@C=?N{#RVo-KU#HCUUWUOwg)gK56pvpz>{yWR=GB-_v
zVXzpRFo0HzUtTn%zKuWJQrsItuI^axRs)`W9OSQf)lDFNMtfhwzs!AK~)ZX90
zZAs3{Z?9BWD}96Oxu2Mr=Le>WN3rbl{#bn9pI7!mr4}5
zpZT)L&f#@6mzP;{VOH24fkLdLwJ4a`_#V@4W2vgJ@Tg)at}5fCDT!832L1Srai}~9
z3v2N`g*oJgW{6Ija1msP3d=Qa0p<(qNR>EL3aZte??Cb23nzNkdaf@_XD4(5t1{~@
zgGw&ZC|B#F=~1XEfrP?3UX>)`U!_b|Z>7i4e?M`zY(LC1ANh}|ZW(q)7Wk2hD0DrI
z@3E~`c`Rd%kUON>m@a*{OpMiU;o0bGA~(^7=TvOMKZZT3@O|h<7h^N>uN`(0k>ePs
zwLw5Ne;sD+dXyHkL5#ny&5X1$(K}H55QUZOXQXt*E%wL{R2?j0w}E
zKfNS4py)z+=%@F;a~
z=`g-Cp$;%x$$V0()`r-0Ro`pB?t%m^nL2)*Lz}`Xu{$5A_O>fPJJOo=rA-;ls^tcL
z6IFB@lt6M<9=($y+kKt%{bcy>ceih@pHjnryvwk&-`#q(XaEf2WDE-}-
zA`bXjA*F&Q*a=*^3X|v+#o^zM%RiF}3Q(kJ8K^Ic*CUiyQzZ|gA1vHxQvDK_*4W~K
zFFvP33MQOVf3(*xT~z2CdYKKs)o9T6!9e+&^cQGxOU}c|Esb*apQl|7`G@1w>n`k0
zyjX_}k>m4e8(+9LTp1It49K1s>W-@#lzwm;YS8ZNlh;eD8hL)1&a-rLk6#)y`!Zyvc)Q?~-{pCo
zQu>@f!S9#$vSFM1%CM%O!-3B4_guErKL8}EtNK<%z>6|UL~^$tpt$fd{KuHJpjkr!
zv4XfviOUvZe4jV`ell&?X1Ei=xSH*2JKjOw7QW`=cewqGsxsFRGg!43H_9en0h#an
zWiHn5_M`mj)jVS{OANodI=T
z{+q(vDo5MAODIkkqfV}DW~eRySa)-S!My<=mRKzrx$$O6`Jpw|1Vy;0$6TG|K!
zs%Wx&=^YN+aqo%!^3=Zlvt{t4Vl-RM!yp{@gtQ`C$MEf3GMQC-oZp+od-dU#LzIWc
z0#!%9+UY}x?|!L|sZSl@d@-Asfg4l89+uko+~B>jxSe-a~IpDok`8pVy;$^3CWM
zlmpGIEuxQ*9j8Z~4R!|qz3CO{fAy+2q3q<^?6<+Qp;ErxR4ePBZpUXhz%VL|WH?)8
zOTcf_cv6dzMDw|5eWnvgh#O%y!|Q1|v8d~AMf(IRqKbugl(*2P<2}grVTDPk8fN+K
z*LM)NA*~Adm5n|lsLZuu3|56+!{^={j67LK<#L#{2&W7&{z>-XrN~9MNVY2ul!6Xi
zBk&$r@Ugt7)3P`-T$uI2U!
zP7>UfeTkiSv49j)AUQYWtAvK@C%UrpD;ymxTiFHCSmP{G{-2{t%$f0*iatFZ<@C_1
z=TTGYrgR*
zEVQ%9XB?gO{RT9RS$)A%E4s0A#-bqquOT*UnIZcf8o>LsxmLqzx~D9x`fCq4@)i_K
z)Ff??yj1+`)lz2Wm^R3Y3D{)6z-0g9mea`OciJu5l~Baw7ctYxeIQ&B-<5&Wyo=wH8qWK{z8RMr|3J>
zQN+s?Zov3Q4OU^nC1$>&K1`{wZT{42}8_2q&a*WuQd2
z(X?%}Whtgm*El)@jfXWejF-KFc-F+Xt;jqFSA19PPKqqV)Y!lCoq(4zfzs`^5}8y8
zjlyt_UI&rKYE%>nv&9w_QiSQo=fQO5CI071mm9t_kaz2qPznh)huv}b&2;+YDJTGi
z8}QDP!(A|Djx0H`*z$3x2HB83by&fNcWUGHLWvJB9x}WbO_^M?PEeiBSh`%I0p9KVjKsl)tav7>STEV*uUUjJ^z#OFW
zQEfESLHyV_TBpm+PeJ5I_kP%nHqMqkhjO?q^?6=5lyMBZ4EC)
z)zw+65ljhbkrH_o8$O)*i^eoBG#F-~lbNbPGLwfVczHNhu;r%cz3W4Z`}KF|g#(N;
zK_2O6jojB8k~xmAy+!gjZEW={KyUan`H!FuL(YmxqbRhWO?Kf)k;+9J)t#OSHQ^bN
zLSFLosisn|7EMb5Yibo=ik5%pd%sCmSHhttgt}JYxNASN|nO>#mA3J&%XuBd(n2Tr8BOXf)aLrJ05Vy|9-|*Ww;;*##
z5L|f-)At6)I?bsra*kYO;618!*^bgn+=s7MG0KFg?+@Z7l^CV>vE(cIs2BUqC94K6
zeM>LYarj`7aq|tbfj#mYFpw>uZW_7yTL3J6VwY8DrAs>qPB^Q*vC1P1M2ibRUrx17
zu8($7$7kj~{=knoK5b*R2-ixNbKo43i?^Niqe#o&E;*%awU^imsYj&Q7)03E(UVSi)Xf6x)I;#Z9y*Mg#S&wp*}t8D6d|=6t@g
zj-O86piSWGd)DIg{_B*G9UiX0Dss6(MfZ(Uz~Zk{uUxEB5#40+&;VBA!9{x$XPi=N
zY8bjL@7qn6eD3Ih)|wl|#?XOqVHADZ^m_6R&!gSsVLVRD@YFOPTK#uqEZ`ZA)d|JM
zAu?9#w5Kt1?3R*~0-3F>VLXfJmcpqWR{r$9cbWSF%j#w`92*m8O+y8CdreO~WK`su
zr}qx3ym?ly4)+e9y@&Tp^)P7${>)UKZSfTPtN29d?Ui`}UsqFaXhH1~YkMsiOWAO!
zS$p+pf3{b+Nzn@etusATLa6!QgYlm$SV&pRdeRQr9qye?$_pFP)M{Be#no3FLc=1#
zmz*lB|Eec+6j|xi6RzjD0{-hS@NoT`oq)^qm4$n>T;OugT7~q)=Otsv0-9vw4YGPM
zFhLz%u8V92r~y`Inf0_~wKu7tDxNHAaNOBIaer?L@~Qhb)V4-zRXob%da
zJO`B^u2kodB{2GLR2GT%K;RC;J7iNI?L$JhE?vI;$KQxuGJ*avNY7*+AB2dme7N>+
zttZc>ayIhC{{I5t|4Ck?^2L~(u0JWqz-q
z8U+5tm9HT|1m9j;)~U(tbCU-h8mM9AN1RpL$H)B>53Zg0(|3YZ9>nip_D#7Y3zH7b
zgLRoSY0q!9Un~347WO5$@QgH}_|NP!)@)q?Hf}z^#W!v`6gdy|rNmLhc=S1h*VxpX?QaDCy5en{d}R^DYThUeM-c>TEa$A4t&Y2uTJXuq677?;Eau1gGi
zE5(Kod;jYZ@N)Zlb=mb;4cunGzMP~Z0I9&E@ztc2hgtEvWTB~+X`?!(Ex<$^AK%y`m;`@;d{)LUljR*7a3v}pLJ^kG?
z{{1~^D%h@umu1i8wf|B;lI7!dl*zl6r##>1`kvAzA#IVuwGRbf5S+ZVPhTH^$I&L^
zmxI062hZ1qzrfSrnd<(E^oF}J(5KAj@oFi0I_`OKjoil2^14OU^=w|p$R@G`UC2I0
z+I>lgc)Cd#pKpIMhZqC3>6woztZ)CTJg8T&1VjJlVaG_LB6G+_Csd5i7yoT2JRGRv^h!7S9Uo-oe8NyKq^Jz0`W9Bg*xqlb6ME2f{_CN@p6tXj#`LO4b$p&TAwSx1-Jf1{i3Odb&7^>N{5*43nm>^Q;s
z#v{G@vaK=IPY)4ezwHb#%8%hHGWJKJq&%cD2i<@eynTv2>
zoZH+8s>!t`Efhqq|8(QtI?XxX=kJ+7Wk5ign{G?j^HmCZz~V
zq1auJf9&BLLj>SX$G{(ipPXEIBIgGZzBHF>X@4=-;BNz~Rc9F24>!hi$@*p7$dFOD
zcD;~1ROat7J7m|7@6F7C=w93t5&YaAz8(@uihyU-V#9u7HLL^P0*Sg
z4Tr34j-b^>LYR|Dm&e=_&~!TBb6duUa8_F3a#JmL`2oWplLnws4iaq&E_|0L{S9Kw
zv`61s-zrh-xbU9Rdz4?oh1wKX>kt3oo?!^X4J|lkI72&T=<2;+yWi9o{ieP3H|^y#
z`5P2LkFR(__HIuupXB62k
z!!pLr=c8t)@PaFYfvs4hQjME=I_}|@oxp*e<*0;9^Rp|lm|9?+p*^ynRT#<9=a{|E
z3WGJF#@#c$^<<{e#$OMoOL@~DKX(pULx&t#!V-}Tl9JvzhW+JbtWfVvkf%}7ZBzeN
zS{;vvI47sm|6w_Q8J@TWcUL!8_)#;ww_s|I2^Wu5+cr(5R%FT6Y~1Svc_A;TtW2Ac
zIWB=T{V@myHF=ar;N#5&Ux7AbMdh{h8?S?ReopI^8Fj@nxFb&A!MIuBEffIV+dlzR
zgk`{IgTib+)lIsqozCeO#n%@zS{x=!c(ja3&4e#!Y-#k=NLa4*!@xRJ#-(NfY&LRo
zZm7{58)W4|I?~<=KzO$xFaFl%3`sPvx@UJGa^92YKsZtN$19nyUSDG+O^|?_R+n
z$US4|F+uh9N2MjZ+wp%vX748QfQ3O`3?&MdJ41Kn^>3Gqr@B>$N(jsga-0!8F}$Cy
zS(mlbu07MuYXtHvP(*mAGF(x*V;UFOJ4I-B7QgMD2?hH)OyMvd
zt|e~5!*wtJq|DgRe4VJI!G(N}k1D0#BSjry&D{4Jz$kaUDC=fxa20eGaH7XIRGF<|
zNf5j>N92n2n8F67g)6A+##lK>tdJAG;TOv}a(#mbvCGyroKe3UWj*;km*7?bo2@O(
zvn$g-IqqR7(>gXmfU{=MFp@HyrV2^1Os34jm1j=k9r-WC0FRaAXh`hMFZ3+qFPl9p
zlx{C+pPy!C&GwEJR+pFsR%~7w;{yaf3|?wh7|!QGKiE^b`L
z&kGr&+-5{hyvD9|Z9rlm;p=ndHJ4scq9+R&wwThD766R|C?_ERbQ~sOL;Qi}7F2@{
z=V(T${4Ui3XEYN8B%h`!v!dN({B%PzlXb}Q5@=xJLF+Vh3Z11H&ENDujY>t&I~v7c?wdTI<6Geu~tV+*yyv}2RBqZ}{2<)Q&dl_7JogVsD|QjrxU
ze6v^-mn2b$
zw-KG_Es8s2v2Ceev#(`Z6bYjJ&i`PQ`vnJ4#zHo=Z@s4hU6
zYS1XHufkv(9IiBW!9RfI>OKtRFpFJWaZ2vWmX!jEFT?xOTG|AeoRy)V>BQz`1SYw)
z5f)YxWPRd-qdhL?s|Je
zpc-II_&;R_)@)L(%Pe*0qz%{KirPE1>T{u@gW+knH`uX~THAS>2$H}{!FiBJr(eHd
zX}LGDiDm^qU)cd@%|VNbt4FXjW;1qN%#)uCe-6PgvS?y{6s*XvFJRKGcXlB(5+5iV
z$!OJ(qLuQwexiaXo7PUwWqCM`yqvtE$P=LR|8_>w3U9$s{
zuPwJGx2@Us4ge>()czqfpXGrdoN6jXp&zJocEzcvrRA$QxOUI)jN!Pq?S8A~GZ3cU
z*G>xAk{*}T<>EdiMD%!qFTpvLIX3vn;FHO!B~WE_iMujNiC@faO{TDb66G;Q4FT51
z_!WKl+Z2yPA9xwq?3q#-e<&E$poUbLpPRgNh${1$xb@DpDAVqks^#EIuK#ngXALPO
zv9i0qYfJ$xA2OP&){XR_A2
zq(AiC1WV4T%>N4=!PHGQQ1+NZEdz>nzq}{3_}4Ej#+H*{A)Foo{amx1>IO`S_zg?N
zaH$U~wVFLxBv&Fu5@{ffUoafS^xId%)7m~q9DhDM{r)(ji
zHsx|6S2qGp`I+F~G1u(oFJM!PwYvt5$0Y
z_Co_RZp=@h6k29V-*q^4iwOJ76_HdwPy`eo{s&U(*Qsx{f}7N(Vv?(5I7q;yz`jrP
zBOJGnM?x~w;7518x*4{glmuEoXiajt%3o)-#RjS$I@rhN1_#UdXcXGqL#ILv{k
zgKTjI>P6NoNLd@59$}&C;#DqarA|8}W$58r!Z&T@QA6;qi1>H$)HEZ0rWnYe?()MU
zcT4->)$$=paH<4jrv#ul|Fua6(#`-lu%qnjE%YK+DbAS%Bppkryug+4aEDoRgtNEr
z{E|%{q5GF~2&y6ZvuNYPui=aQ-Lzh-!p~ZtI12kL73re;+YYOaL+XcnzjL3P0Valo3b^v6U}8EhTA
z8BsuF!M*kgqX;k^mZMR#xmu$YSDWE_P{SbmhO6k0;}R
zE*>7fMEpc3YFpJxAUXtHt49go^Z+{;t=
z|LiD*$ttY}{Y)BcTCRNlH$4p%UC=l%Ym;JoC8+s=;{%+GSUe=UcK=Eqqe#qOF^-^_
zLJ3^rudVPPH8ptoZfQFHoou|E^-u%!w%IY?4EK3+>~+<--f>~w>FOXKl
za**{AanJ%CLl*9y>h41kQuAD#bk*zwoh7z;qxAKvB*dtRwMyZn9D5|?>(6p|?sk?q
z5%{To=rwG5p?Ih7!8^^((Kf2`tZ$kd%ZO)$^QR)O1a}uWNiSR1x(iL$$m7*2@?QrZvd`>I``nMQzPw1Sf3@)F!
zdrOA(FUfE`=1j{~8mQ>+_e8|19(9`DQs@t~x<0yg4#UvcLp;1#8r&YNe3#iF`18^X
zb`d5(iv$$s_y>pGlhmq$0RtTNONXnPa`z+PXyrlOtgmfHA0A#-#X$EvDFdpdv`%Qx
z8Iv*glh%e1Xin=IPtvJSV8td?Y#Bm|e({1yGFM6a@K;n=GkN-EXjzIzpp1EBSSJbU
zZcd?%QzC^6Wjnt+jmvyyxx@MsP+}WW{viWM^Q>ct_|`~~m*VeeP3w52FS%xC&>@{&
zHo&_j=IvIXbNQla#-OQ!9J!7-$OaDUoKUWwN&kk2r-jXL^hcFUM^WT1YPJk7rUxIM
zNZbiGA{aiP4(+o&k70-H`om%O59^eF?}>{o3_(kEM1sXNe7@yEkpgocaqPo(VZ$sy
zm?0$>PmDX+*Im2q!mbTeUVl#^U7drtwa`EHDCc+R$a|QZ;~>`a^Po`9HSWX$cN%SU
z|I|LOSZ@Wsi^#yXbHz8{9JFm1;_p;D6b?zjO+wIL;Jt;LgG$w*`jG~*yYoNk_A-x&
z@2^UCF_6kOD!u>O(W(A@d7d
zZU~)7juBG~2`hv}YM^)oTfk=6XU|vdKGH~Ee~0u-MK=4f^3m>Rcv`xQ=zt8GUXpdi
zRjT-kqL8qs-vcxdHe(n3dtp~(j**@KTS@;R%~Mu#78=fs{?Iy|?~L?7M3Mi!M&>T9
z{gj<5JSZcB8C88~$U=Z!H=~48S4R^1W1zO7IJ!~k1&Mej=Mtq$h2`$Ib?u_
zXaNHzwv=g&_lv8hL#_F>1ZK&PYD#%Pb6;h>$s-d~1=3ZBjh&TqV
z?K+7^V=E*4L^Z|+H-ZpBqB8Y5u%q^cKQw;P2!CO9FG8G_Oe;a$Atq9y0Fqh>sdvzDcR}R2ucS4?y9VUe&A4jh$4;t
zD`!F-aZ!-MdBuX2VXnrEhuM_CO#B+d{!)()VkY~wQOMEZ^vMt-_;SAJ@f}08W+>Bi
z2x3Rk6mn{F%=avBZD|f;4KFHJ#7T>`mCDu;3-sqk3Jl&Q)Li8XPKV0*vIu=t+w{@d
zqPQXqp|}YXPb&(bR-C8&6Vex#H6dSCbsN2QXLKai?%CbUnXM*9a!MRe^gwb(L_gI8
zV@UW{1MKL|Ns0=p#`{-LQQpwtPGaZgMS0`JuZcobknk4}+b3bCiTgvD?eb>hr(E1jZ
z+Kg93dgH=_2a7!Um$bf=fsN%$#GPz^U$l-w%V*6Oh)GQ271miEm`(3<0lY(T*V%jK
ziyR?LEmeH=Irtl^?MM%P-_PUd#mAp=WJ=Y+PXmPP4VW3E>uo*&6VIpWzQuD!IVlix!5^*mxeXFWFQr+f{cI2-{P~2!KU%B~!lnyedUiQHDNsiziirt&tiKxYp@r+zC{JivYU&nJ%
zUu&R~w-%{!o+|9!Jmi8A!+I@Cy~I~!!;haM{w4x43gAup3DY)@>;FiMzk=acg`7Md?`e{f
ze|~+b^x!OG1!L#Q=h+#^#5`4PJ;BI~Z&rBkYZEksXZdA-T}cvcOHp+?`FR40)x=tV
zca?wXpzbi)SZcipanKaTXPBX
z(YRaP$xb@q*X1^!