highbeats_pdf/simulationexamples.py
2020-12-01 12:00:50 +01:00

432 lines
21 KiB
Python

from plot_eod_chirp import find_lm,find_periods,integrate_chirp,global_maxima,rectify,find_dev,power_func,power_parameters
from plot_eod_chirp import rectify, find_beats,find_dev,find_times,global_maxima, find_lm, conv
from functionssimulation import single_stim,pl_eods
from matplotlib.colors import LinearSegmentedColormap
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 myfunctions import default_settings
import os
import string
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, chirp = True):
time, time_cut, cut = find_times(left_c[g], right_c[g], sampling, deviation_s[d])
eod_fish_e, eod_fish_r, period_fish_r, period_fish_e = find_periods(a_fe, time, eod_fr, a_fr, eod_fe, e)
#embed()
if chirp == False:
eod_fe_chirp = eod_fish_e
else:
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]
threshold_cube = (eod_rec_up)**3
#embed()
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 threshold_cube , 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_func2(bef_c = -2, aft_c = -1, a_fr = 1, factor = 200, sampling = 10000, phase_zero = np.arange(0, 2 * np.pi, 2 * np.pi / 10),new = True, eod_fe = [50], beat_corr = [50], beats = [50],eod_fr = 500,size = [120] ,delta_t = 0.014,a_fe = 0.2, win = 'w2',deviation = [150], save = False):
print('simulation')
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
# maximal frequency excursion during chirp / 60 or 100 here
# eod fish reciever
# amplitude fish reciever
#sampling = 500
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)
if (not os.path.exists('power_simulation.pkl')) or (new == True):
results = [[]] * 1
for d in range(len(deviation)):
bef_c = bef_c
aft_c = aft_c
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)):
#embed()
left_b = [-0.3*sampling]*len(eod_fe)
right_b = [-0.1 * sampling]*len(eod_fe)
threshold_cube, 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
if sampling>1000:
#nfft = int((4096 * sampling/ 10000) * 2)
nfft = int(4096 * sampling / 40000 * 2)
else:
nfft = 4096
#embed()
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]
name = ['global max']
var = [maxima_interp_b]
samp = [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(samp)):
#print(i)
plot = False
pp, f = ml.psd(var[i] - np.mean(var[i]), Fs=samp[i], NFFT=nfft,
noverlap=nfft / 2)
#plot = True
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()
#eod_fr2 = eod_fr*3
eod_fr2 = eod_fr
try:
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_fr2])]])
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_fr2])*np.abs(f[1]-f[0]))])
else:
results[i]['f'].extend(list([f[np.argmax(pp[f < 0.5 *eod_fr2])]]))
#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_fr2]) * np.abs(f[1] - f[0]))]))
except:
embed()
if save:
#embed()
results = pd.DataFrame(results)
results.to_pickle('power_simulation.pkl')
np.save('power_simulation.npy', results)
else:
results = pd.read_pickle('power_simulation.pkl')
return results
def plot_power(a_fe, a_fr,beats, beat_corr, eod_fe, sampling, ax, eod_fr,nr_size = 11, left = 0.15,nrs = [6,7], bef_c = -2, aft_c = -1):
plt.rcParams['figure.figsize'] = (3, 3)
plt.rcParams["legend.frameon"] = False
colors = ['black', 'magenta', 'pink', 'orange', 'moccasin', 'red', 'green', 'silver']
colors = ['red','pink']
#eod_fr = 500
#start = 5
#end = 2500
#step = 5
#eod_fe, beat_corr, beats = find_beats(start, end, step, eod_fr)
#embed()
#beats = eod_fe - eod_fr
factor = 200
#sampling = eod_fr * factor,
#(a_fr
# = 1, a_fe = 0.2, eod_fr = eod_fr, eod_fe = eod_fe, win = 'w2', deviation_s = deviation_s, sigma = sigma, sampling = 100000, deviation_ms = deviation_ms, beat_corr = beat_corr, size =[120], phase_zero =[phase_zero[0]], delta_t = delta_t, deviation_dp = deviation_dp, show_figure = True, plot_dist = False, save = False, bef_c = -2, aft_c = -1)
results = power_func2(bef_c = bef_c, aft_c = aft_c, a_fe = a_fe, a_fr = a_fr, sampling = sampling, phase_zero = [0], eod_fe = eod_fe, beat_corr = beat_corr, eod_fr = eod_fr, save = True)
#results = [results.iloc[5]]
results = [results.iloc[0]]
#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])
#embed()
ax[0].plot(beats / eod_fr + 1, np.array(results[i]['f']) / eod_fr + 1, color=colors[i])
ax[0].text(-left, 1.1, string.ascii_uppercase[nrs[0]], transform=ax[0].transAxes,
size=nr_size, weight='bold')
# 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[1].text(-0.1, 1.1, string.ascii_uppercase[nrs[1]], transform=ax[1].transAxes,
size=nr_size, weight='bold')
#remove_tick_marks(ax[1])
remove_tick_marks(ax[0])
#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('[Hz]')
ax[1].set_ylabel('Modulation')
#ax[0, 2].set_title('Modulation depth (same scale)')
for i in [0,1]:
ax[1].set_xlabel('EOD multiples')
ax[i].spines['right'].set_visible(False)
ax[i].spines['top'].set_visible(False)
ax[i].set_xlim([0,5])
#plt.subplots_adjust(bottom = 0.2)
#plt.savefig('localmaxima.pdf')
#plt.savefig('../highbeats_pdf/localmaxima.pdf')
#plt.show()
def plot_eod(beat_corr_upper, beats_upper,eod_fe_upper,beat_corr, minus_bef ,plus_bef,row,col,a_fe, delta_t, a_fr, size,beat_corr_org,s,p, amplitude, phase_zero, deviation_ms, eod_fe,deviation, eod_fr, sampling,factor,fs = [6.2992, 4], shift_phase = 0, bef_c = -2, aft_c = -1):
beats = eod_fe - eod_fr
nr_size = 11
plt.rcParams['axes.titlesize'] = 10
plt.rcParams['axes.labelsize'] = 10
plt.rcParams['font.size'] = 8
plt.rcParams["legend.frameon"] = False
beat_corr_col = 'goldenrod'
df_col = 'steelblue'
interval = np.hstack([np.linspace(0.25, 0.97)])
colors = plt.cm.Blues(interval)
colors = plt.cm.Reds(interval)
cmap = LinearSegmentedColormap.from_list('name', colors)
colors = cmap(np.linspace(0, 1, len(eod_fe)))
d = 0
fig = plt.figure(figsize=fs)#hspace = 0.55
outer = gridspec.GridSpec(5, 1, top = 0.9, hspace = 0, wspace = 0, height_ratios = [2,1.2,3.2,1.2,2])
#fig, ax = plt.subplots(nrows=row, ncols=col, figsize=fs)
ax = {}
#shift_phase = 0
#lower = gridspec.GridSpecFromSubplotSpec(1, 1, hspace = 0.63, subplot_spec=outer[0])
ax['fill_up'] = plt.subplot(outer[1])
ax['fill_up'].spines['right'].set_visible(False)
ax['fill_up'].spines['top'].set_visible(False)
ax['fill_up'].spines['left'].set_visible(False)
ax['fill_up'].spines['bottom'].set_visible(False)
ax['fill_lower'] = plt.subplot(outer[3])
ax['fill_lower'].spines['right'].set_visible(False)
ax['fill_lower'].spines['top'].set_visible(False)
ax['fill_lower'].spines['left'].set_visible(False)
ax['fill_lower'].spines['bottom'].set_visible(False)
left = 0.15
grid0 = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=outer[2], height_ratios =[2,1,1], wspace=0.02, hspace=0.2) # ,
ax['upper'] = plt.subplot(grid0[0])
#start =
#embed()
remove_tick_marks(ax['upper'])
#embed()
ax['upper'].plot((eod_fe_upper - eod_fr)/eod_fr+1, (beat_corr_upper- eod_fr)/eod_fr+1, color = beat_corr_col)
ax['upper'].set_xlabel('EOD multiples')
ax['upper'].set_ylabel('[Hz]')
ax['upper'].plot((eod_fe_upper - eod_fr) / eod_fr + 1,np.abs(eod_fe_upper - eod_fr) / (eod_fr), color=df_col, alpha=0.7,zorder = 1)
#ax['upper'].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))
#ax['upper'].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))
ax['upper'].set_xticks(
np.arange((eod_fe_upper[0] - eod_fr) / eod_fr + 1, (eod_fe_upper[-1] - eod_fr) / eod_fr + 1, (eod_fr / 2) / eod_fr))
ax['upper'].set_yticks(
np.arange((eod_fe_upper[0] - eod_fr) / eod_fr + 1, (eod_fe_upper[-1] - eod_fr) / eod_fr + 1, (eod_fr / 2) / eod_fr))
nrs = [5]
ax['upper'].text(-left, 1.1, string.ascii_uppercase[nrs[0]], transform=ax['upper'].transAxes,
size=nr_size, weight='bold')
plt.grid()
colors = np.concatenate([['hotpink']*col,['red']*col,['darkred']*col])
plt.xlim([(np.min(eod_fe_upper)- eod_fr)/eod_fr+1, (np.max(eod_fe_upper)- eod_fr)/eod_fr+1])
#ax[d] = fig.add_subplot(row,col, 1)#int(np.ceil(np.sqrt(len(eod_fe))))
#sampling_rate = sampling
plt.ylim([-0.2,1])
ax[0] = plt.subplot(grid0[1])
ax[1] = plt.subplot(grid0[2])
plot_power(a_fe, a_fr, beats_upper, beat_corr_upper, eod_fe_upper, sampling, ax, eod_fr,left = left, nr_size = nr_size,bef_c = bef_c, aft_c = aft_c)
upper = gridspec.GridSpecFromSubplotSpec(int(row / 2), col, hspace=0.5, subplot_spec=outer[0])
nfft = 4096
nfft = int(4096 * sampling / 10000 * 2)
full = eod_fe
colors_full = colors
#embed()
eod_fe = full[5:10]
colors_plot = colors_full[5:10]
colors = colors_full[5:10]
colors = ['black'] * len(eod_fe)
#colors = ['black']*len(eod_fe)
nrs = [0, 1, 2, 3, 4, 5]
for e in range(len(eod_fe)):
f_max,lims,ax = single_stim(ax,colors_plot, int(row/2), col, eod_fr, eod_fe, e, upper,s = s, p = p, d = d, df_col = df_col, factor = factor, beat_corr_col = beat_corr_col, nfft = nfft, minus_bef = minus_bef, delta_t = delta_t, sampling = sampling,
deviation = deviation,plus_bef = plus_bef, a_fr = a_fr, phase_zero = phase_zero, shift_phase = shift_phase, size = size, a_fe = a_fe)
ax[e].text(-left, 1.1, string.ascii_uppercase[nrs[e]], transform=ax[e].transAxes,
size=nr_size, weight='bold')
x = (eod_fe[e] - eod_fr)/eod_fr
ax['upper'].scatter(x, (f_max)/eod_fr, color=colors_full[e], s=19,zorder = 2)
ax['fill_up'].scatter(x, 1, marker = 'v',color=colors[e], s=19,zorder = 2)
ax['fill_up'].plot([x,x], [3,1], color=colors[e], zorder = 2)
axis = (eod_fe_upper - eod_fr) / eod_fr
#embed()
ax['fill_up'].set_xlim([axis[0],axis[-1]])
ax['fill_up'].set_ylim([0.72,6])
ax['fill_up'].set_yticks([])
lower = gridspec.GridSpecFromSubplotSpec(int(row / 2), col, hspace=0.5, subplot_spec=outer[4])
eod_fe = full[0:5]
colors_plot = colors_full[0:5]
colors_plot = colors_full[5:10]
colors = colors_full[5:10]
colors = ['black'] * len(eod_fe)
nrs = [8,9, 10, 11, 12, 13, 14]
for e in range(len(eod_fe)):
f_max,lims,ax = single_stim(ax,colors_plot, int(row/2), col, eod_fr, eod_fe, e, lower, s = s, p = p, d = d, df_col = df_col, factor = factor, beat_corr_col = beat_corr_col, nfft = nfft, minus_bef = minus_bef, delta_t = delta_t, sampling = sampling,
deviation = deviation,plus_bef = plus_bef, a_fr = a_fr, phase_zero = phase_zero, shift_phase = shift_phase, size = size, a_fe = a_fe)
ax[e].text(-left, 1.1, string.ascii_uppercase[nrs[e]], transform=ax[e].transAxes,
size=nr_size, weight='bold')
ax['upper'].scatter((eod_fe[e] - eod_fr)/eod_fr+1, (f_max)/eod_fr, color=colors_full[e], s=19, zorder = 2)
#ax['upper'].scatter((eod_fe[e] - eod_fr) / eod_fr, -0.17, marker = '^', color=colors[e], s=19)
x = (eod_fe[e] - eod_fr) / eod_fr
ax['fill_lower'].scatter(x, 4.5, marker = '^',color=colors[e], s=19, zorder = 2)
ax['fill_lower'].plot([x,x], [1.5,4.5], color=colors[e], zorder = 2)
axis = (eod_fe_upper - eod_fr) / eod_fr
#embed()
ax['fill_lower'].set_xlim([axis[0],axis[-1]])
ax['fill_lower'].set_ylim([0.5,7.3])
ax['fill_lower'].set_yticks([])
#for e in range(int(len(eod_fe)/2)):
# f_max,lims = single_stim(factor, beats, beat_corr, plus_bef,ax,s,df_col, beat_corr_col, colors, lower, row, col, fig, nfft, p, minus_bef, delta_t, sampling,
# deviation, d, eod_fr, eod_fe, e, a_fr, phase_zero, shift_phase, size, sigma, a_fe)
# ax['upper'].scatter((eod_fe[e] - eod_fr)/eod_fr, (f_max)/eod_fr, color=colors[e], s=19)
#plt.xlabel('Time [ms]')
plt.legend(loc=(-5, -0.7),ncol = 4)
plt.subplots_adjust(bottom = 0.1, top = 0.8,left = 0.13)
#embed()
#scalebar = ScaleBar(50,'ms')
#plt.gca().add_artist(scalebar)
#embed()
# plt.xlim([-200,200])
#plt.show()
return fig
if __name__ == "__main__":
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)
a_fr = 1 # amplitude fish reciever
amplitude = a_fe = 0.5 # amplitude fish emitter
factor = 200
eod_fr =500 # eod fish reciever
eod_fr = 537
eod_fr = 734
eod_fr = 500
initial = np.array([60, 310])
stim_matrix = np.transpose([initial, initial + eod_fr, initial + eod_fr*2, initial + eod_fr*3, initial + eod_fr*4])
eod_fe = np.concatenate(stim_matrix)
sampling = eod_fr * factor
sampling = 112683 #122300 #500
sampling_fish = 1000
start = 5
step = 10
eod_fe_upper = np.arange(start,np.int(np.max(eod_fe)/eod_fr)*eod_fr+eod_fr,step)
beats_upper = eod_fe_upper - eod_fr
beat_corr_upper = eod_fe_upper % eod_fr
beat_corr_upper[beat_corr_upper > eod_fr / 2] = eod_fr - beat_corr_upper[beat_corr_upper > eod_fr / 2]
# start = 0
# end = 2500
# step = 10
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)
#embed()
pp = [1, 3, 6, 9] # [2,7]
minus_bef = -20
plus_bef = -10
#minus_bef = -10
#plus_bef = -10
bef_c = -2
aft_c = -1
# embed()
# eod_fe = np.array([initial,60+500,310+500,1010,60+1000,310+1000,1010+500,60+1000+500,310+1000+500])
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]
row, col = np.shape(stim_matrix)
fig = plot_eod(beat_corr_upper, beats_upper,eod_fe_upper,beat_corr, minus_bef, plus_bef, row, col, a_fe, delta_t, a_fr, size, beat_corr, 0, 0, a_fe, phase_zero,
deviation_ms, eod_fe, deviation_dp, eod_fr, sampling, factor, fs=[6.28, 6],
shift_phase=(2 * np.pi) / 4, bef_c = bef_c, aft_c = aft_c) # 6.2992,6.69
#fig.savefig()
plt.savefig(
fname= 'simulationexamples.pdf')
plt.savefig('../highbeats_pdf/simulationexamples.pdf')
plt.show()