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


def fish_indirect(delta_t, amplitude, size, beat, phase_zero, sampling, show_figure = False):
    for a in range(len(amplitude)):
        for s in range(len(size)):
            for p in range(len(phase_zero)):
                for b in range(len(beat)):

                    # sigma calculation
                    sigma = delta_t /math.sqrt((2*math.log(10))) # width of the chirp

                    # time calculation
                    #if np.abs(1 / beat[b]) * 10 < 6 * delta_t:
                    time = np.arange(-3 * delta_t, 3 * delta_t, 1 / sampling)
                    #else:
                    #    time = np.arange(-np.abs(1/beat[b]) * 3, abs(1/beat[b]) * 3, 1/ sampling_frequency)

                    # frequency calculation
                    chirp = size[s] * np.exp((-(time**2)/(2*sigma**2)))
                    frequency = beat[b] + chirp

                    # phase calculation
                    I1 = []
                    for i in range(len(time)):
                        y1 = ((np.pi**0.5)/2)*math.erf(time[i]/sigma) - ((np.pi**0.5)/2)*math.erf(-np.inf)
                        I1.append(y1)
                    phase = 2 * np.pi * beat[b] * time + 2 * np.pi * size[s] * sigma * np.array(I1) + phase_zero[p]

                    # am calculation
                    am = amplitude[a] * np.cos(phase)
                    # plot
                    plt.plot(time*1000, am, label='Amplitude Modulation', color='red', linewidth=2)
                    plt.axvline(x=-7.5, color='black', linewidth=1, linestyle='-')
                    plt.axvline(x=7.5, color='black', linewidth=1, linestyle='-')
                    plt.xlabel('time [ms]')
                    plt.ylabel('AM [mV]')
                    plt.title('beat(%5.2f' % (beat[b]) + 'Hz) ' + 'size(%5.2f' % (size[s]) + 'Hz) ' + 'phase(%5.2f' % (phase_zero[p]) + ') ' + 'amplitude(%1.2f' % (amplitude[a]) + ')')
                    #plt.show()
                    plt.savefig(fname='../results/' + 'simulations/' +  'amplitude(%1.2f' % (amplitude[a]) +')_size(%5.2f' % (size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')_beat(%5.2f' % (beat[b]) + ')_AM_fish%5.2f' % sampling + '.png')
                    if show_figure:
                        plt.show()
                    else:
                        plt.close()

def stimulus_indirect(delta_t,amplitude ,size, beat, phase_zero, show_figure = False):
    for a in range(len(amplitude)):
        for s in range(len(size)):
            for p in range(len(phase_zero)):
                for b in range(len(beat)):
                    # sigma calculation
                    sigma = delta_t /math.sqrt((2*math.log(10))) # width of the chirp

                    # time assignment
                    if np.abs(1 / beat[b]) * 2 < 3 * delta_t:
                        time = np.arange(-1.5 * delta_t, 1.5 * delta_t, 1 / 40000)
                    else:
                        time = np.arange(-np.abs(1/beat[b]) * 1, abs(1/beat[b]) * 1, 1/40000)

                    # frequency calculation
                    chirp = size[s] * np.exp((-(time**2)/(2*sigma**2)))
                    frequency = beat[b] + chirp

                    # phase calculation
                    I1 = []
                    for i in range(len(time)):
                        y1 = ((np.pi**0.5)/2)*math.erf(time[i]/sigma) - ((np.pi**0.5)/2)*math.erf(-np.inf)
                        I1.append(y1)
                    phase = 2 * np.pi * beat[b] * time + 2 * np.pi * size[s] * sigma * np.array(I1) + phase_zero[p]

                    # am calculation
                    am = amplitude[a] * np.cos(phase)
                    plt.plot(time*1000, am, label='Amplitude Modulation', color='red', linewidth=2)
                    plt.axvline(x=-7.5, color='black', linewidth=1, linestyle='-')
                    plt.axvline(x=7.5, color='black', linewidth=1, linestyle='-')
                    plt.xlabel('time [ms]')
                    plt.ylabel('AM [mV]')
                    plt.title('beat(%5.2f' % (beat[b]) + 'Hz) ' + 'size(%5.2f' % (size[s]) + 'Hz) '+ 'phase(%5.2f' % (phase_zero[p]) + ') '+ 'amplitude(%1.2f' % (amplitude[a]) + ')')
                    plt.savefig(fname='../results/' + 'simulations/'  +  'amplitude(%1.2f' % (amplitude[a]) +')_size(%5.2f' % (size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')_beat(%5.2f' % (beat[b]) + ')_AM_stimulus.png')
                    if show_figure:
                        plt.show()
                    else:
                        plt.close()

def create_distance(conv_chirp,conv_beat,a):
    dm_wo_max = np.sqrt(np.mean(np.transpose(
        (np.transpose(conv_chirp) - np.mean(conv_chirp, a) - np.transpose(conv_beat) + np.mean(conv_beat, a)) ** 2), a))
    range_chirp = np.max(conv_chirp,a)-np.min(conv_chirp,a)/np.sqrt(2)
    range_beat = np.max(conv_beat,a)-np.min(conv_beat,a)/np.sqrt(2)
    dm_metz = dm_wo_max/np.max(np.transpose([range_chirp,range_beat]),a)
    dm_wo_mean2 = np.sqrt(np.mean((conv_chirp - conv_beat) ** 2,a)) / np.max(np.transpose([np.max(conv_chirp,a)-np.min(conv_chirp,a)/np.sqrt(2),np.max(conv_beat,a)-np.min(conv_beat,a)/np.sqrt(2)]),a)
    #dm_wo_mean2_max = np.sqrt(np.mean((conv_chirp - conv_beat) ** 2,a))
    dm_wo_mean2_max  = mean_euc_d(conv_beat, conv_chirp, a)
    d_euclidean = euc_d(conv_beat, conv_chirp, a)
    #dst = distance.euclidean(conv_beat, conv_chirp)
    return dm_metz, dm_wo_mean2,dm_wo_max,dm_wo_mean2_max,range_chirp,range_beat, d_euclidean

def snippets2(peaks_interp, maxima_interp, am_indirect, am, d, period_shift, period_distance, middle_am, eod_convolved_up, middle_conv, delta_t, sampling, deviation, time, eod_conv_downsampled):
    conv_chirp = eod_convolved_up[
                 middle_conv - int(0.5 * period_distance * (sampling)):middle_conv + int(0.5 * period_distance * (sampling))]

    peaks_interp_chirp = peaks_interp[
                 3*deviation[d]+middle_conv - int(0.5 * period_distance * (sampling)):3 * deviation[d] + middle_conv + int(0.5 * period_distance * (sampling))]
    maxima_interp_chirp = peaks_interp[
                 3*deviation[d]+middle_conv - int(0.5 * period_distance * (sampling)):3 * deviation[d] + middle_conv + int(0.5 * period_distance * (sampling))]
    am_indirect_chirp = am_indirect[
                 middle_conv - int(0.5 * period_distance * (sampling)):middle_conv + int(0.5 * period_distance * (sampling))]

    time_chirp = time[middle_conv - int(0.5 * period_distance * (sampling)):middle_conv + int(
        0.5 * period_distance * (sampling))]

    am_fish_chirp = am[middle_am - int(0.5 * period_distance * (eod_fr)):middle_am + int(
        0.5 * period_distance * (eod_fr))]
    conv_chirp_corr = eod_conv_downsampled[middle_am - int(0.5 * period_distance * (eod_fr)):middle_am + int(
        0.5 * period_distance * (eod_fr))]

    shift = int(period_shift * (sampling) * 2)
    shift_fish = int(period_shift * (eod_fr) * 2)

    am_fish_beat = am[middle_am - int(0.5 * period_distance * (eod_fr)) - shift_fish:middle_am + int(
        0.5 * period_distance * (eod_fr)) - shift_fish]

    conv_beat = eod_convolved_up[middle_conv - int(0.5 * period_distance * (sampling)) - shift:middle_conv + int(
        0.5 * period_distance * (sampling)) - shift]
    peaks_interp_beat = peaks_interp[3 * deviation[d] + middle_conv - int(0.5 * period_distance * (sampling)) - shift:3 * deviation[d] + middle_conv + int(
        0.5 * period_distance * (sampling)) - shift]
    maxima_interp_beat = maxima_interp[3 * deviation[d] + middle_conv - int(0.5 * period_distance * (sampling)) - shift:3 * deviation[d] + middle_conv + int(
        0.5 * period_distance * (sampling)) - shift]
    am_indirect_beat = am_indirect[middle_conv - int(0.5 * period_distance * (sampling)) - shift:middle_conv + int(
        0.5 * period_distance * (sampling)) - shift]

    time_beat = time[middle_conv - int(0.5 * period_distance * (sampling)) - shift: middle_conv + int(
        0.5 * period_distance * (sampling)) - shift]
    return am_fish_beat,maxima_interp_beat,peaks_interp_beat,maxima_interp_chirp, peaks_interp_chirp,am_indirect_beat,am_indirect_chirp,time_beat,conv_beat,conv_chirp_corr,am_fish_chirp,time_chirp,conv_chirp


def snippets(am_indirect,am, d,period, middle_am,eod_convolved_up,middle_conv,delta_t,sampling,deviation,time,eod_conv_downsampled):
    if period <= delta_t*2:
        conv_chirp = eod_convolved_up[
                     middle_conv - int(delta_t * (sampling)):middle_conv + int( delta_t * (sampling))]
        time_chirp = time[middle_conv - int( delta_t * (sampling)):middle_conv + int(
             delta_t * (sampling))]
        am_indirect_chirp = am_indirect[
                     middle_conv - int(delta_t * (sampling)):middle_conv + int(delta_t * (sampling))]
        am_fish_chirp = am[middle_am - int(delta_t * (eod_fr)):middle_am + int( delta_t * (eod_fr))]
        conv_chirp_corr = eod_conv_downsampled[middle_am - int( delta_t * (eod_fr)):middle_am + int(
             delta_t * (eod_fr))]

        if period == 0:
            shift = int((delta_t)*3* sampling)
        else:
            shift = int((5 * period) * sampling)*2
        #if period == 0:
        #    shift = int(((delta_t + 3 * period) * sampling))
        #else:
        #    shift = int(((np.ceil(delta_t / period) * (period) + 3 * period) * sampling))
        conv_beat = eod_convolved_up[middle_conv - int(delta_t * (sampling)) - shift:middle_conv + int(
             delta_t * (sampling)) - shift]
        am_indirect_beat = am_indirect[middle_conv - int( delta_t * (sampling)) - shift:middle_conv + int(
             delta_t * (sampling)) - shift]
        time_beat = time[
                    middle_conv - int( delta_t * (sampling)) - shift:middle_conv + int(
                         delta_t * (sampling)) - shift]

    else:
        conv_chirp = eod_convolved_up[
                     middle_conv - int( period * (sampling)):middle_conv + int( period * (sampling))]
        am_indirect_chirp = am_indirect[
                     middle_conv - int(period * (sampling)):middle_conv + int( period * (sampling))]

        time_chirp = time[middle_conv - int( period * (sampling)):middle_conv + int(
             period * (sampling))]

        am_fish_chirp = am[middle_am - int( period * (eod_fr)):middle_am + int(
            period * (eod_fr))]
        conv_chirp_corr = eod_conv_downsampled[middle_am - int( period * (eod_fr)):middle_am + int(
             period * (eod_fr))]

        shift = int(period * (sampling) * 2)
        conv_beat = eod_convolved_up[middle_conv - int( period * (sampling)) - shift:middle_conv + int(
             period * (sampling)) - shift]
        am_indirect_beat = am_indirect[middle_conv - int(period * (sampling)) - shift:middle_conv + int(
             period * (sampling)) - shift]

        time_beat = time[middle_conv - int( period * (sampling)) - shift: middle_conv + int(
             period * (sampling)) - shift]
    return am_indirect_beat,am_indirect_chirp,time_beat,conv_beat,conv_chirp_corr,am_fish_chirp,time_chirp,conv_chirp



def global_maxima_array(period_fish_e,period_fish_r,eod_rectified_up):
    period_length = np.max([period_fish_e, np.ones(len(period_fish_e))*len(period_fish_r)],0)
    split_windows = np.linspace(period_length, len(eod_rectified_up), len(eod_rectified_up)/period_length)
    splits = np.split(eod_rectified_up, split_windows)
    steps = np.arange(0, len(eod_rectified_up), len(splits[0]))
    maxima_values = np.max(splits[0:-1], 1)
    maxima_index = np.argmax(splits[0:-1], 1)
    maxima_index = maxima_index + steps[0:-1]
    maxima_interp = np.interp(np.arange(0, len(eod_rectified_up), 1), maxima_index, maxima_values)
    return maxima_values,maxima_index, maxima_interp

def global_maxima(sampling, eod_fr, period_fish_e,period_fish_r,eod_rectified_up):
    #period_length = max(len(period_fish_e), len(period_fish_r))
    period_length = len(period_fish_r)
    period_length = int(np.round((1/eod_fr)*sampling))
    #embed()
    if period_length >len(eod_rectified_up):
        maxima_values = np.max(eod_rectified_up)
        maxima_index = np.argmax(eod_rectified_up)
        maxima_interp = [maxima_values]*len(eod_rectified_up)
    else:
        split_windows = np.arange(period_length, len(eod_rectified_up), period_length)
        splits = np.split(eod_rectified_up, split_windows)
        steps = np.arange(0, len(eod_rectified_up), len(splits[0]))
        maxima_values = np.max(splits[0:-1], 1)
        maxima_index = np.argmax(splits[0:-1], 1)
        maxima_index = maxima_index + steps[0:-1]
        maxima_interp = np.interp(np.arange(0, len(eod_rectified_up), 1), maxima_index, maxima_values)
    return maxima_values,maxima_index, maxima_interp



def snippets3(skretch,peaks_interp, maxima_interp, am_indirect,am_fish, d,period, middle_am,eod_convolved_up,middle_conv,delta_t,sampling,deviation,time,eod_conv_downsampled):
    #period_distance = 0.1
    period_dist = []
    if period <= delta_t*skretch:
        if period == 0:
            period_shift = delta_t*skretch
            period_distance =  delta_t*skretch
        else:
            period_shift = np.ceil(delta_t*skretch/period)*period
            period_distance = np.ceil(delta_t*skretch/period)*period
        am_fish_beat,maxima_interp_beat,peaks_interp_beat,maxima_interp_chirp, peaks_interp_chirp,am_indirect_beat,am_indirect_chirp,time_beat,conv_beat,conv_chirp_corr,am_fish_chirp,time_chirp,conv_chirp = \
            snippets2(peaks_interp, maxima_interp, am_indirect, am_fish, d, period_shift,period_distance, middle_am, eod_convolved_up, middle_conv,
            delta_t, sampling, deviation, time, eod_conv_downsampled)
    else:
        period_shift = period
        period_distance = period
        am_fish_beat,maxima_interp_beat,peaks_interp_beat,maxima_interp_chirp, peaks_interp_chirp,am_indirect_beat,am_indirect_chirp,time_beat,conv_beat,conv_chirp_corr,am_fish_chirp,time_chirp,conv_chirp = \
            snippets2(
            peaks_interp, maxima_interp, am_indirect, am_fish, d, period_shift,period_distance, middle_am, eod_convolved_up, middle_conv,
            delta_t, sampling, deviation, time, eod_conv_downsampled)
    period_dist.append(period_distance)
    return  period_dist, am_fish_beat,maxima_interp_beat,peaks_interp_beat,maxima_interp_chirp, peaks_interp_chirp,am_indirect_beat,am_indirect_chirp,time_beat,conv_beat,conv_chirp_corr,am_fish_chirp,time_chirp,conv_chirp

def conv(eod_fr, sampling, cut,deviation, eod_rectified_up, eod_rectified_down):
    if deviation* 5 % 2:
        points = deviation * 5
    else:
        points = deviation * 5 - 1
    gaussian = signal.gaussian(points, std=deviation, sym=True)
    gaussian_normalised = (gaussian * 2) / np.sum(gaussian)
    length_convolved = int(len(gaussian_normalised) / 2)
    eod_convolved_up = np.convolve(gaussian_normalised, eod_rectified_up)
    eod_convolved_up = eod_convolved_up[length_convolved + cut:-length_convolved - cut]
    eod_convolved_down = np.convolve(gaussian_normalised, eod_rectified_down)
    eod_convolved_down = eod_convolved_down[length_convolved + cut:-length_convolved - cut]
    middle_conv = int(len(eod_convolved_up) / 2)
    eod_conv_downsampled = eod_convolved_up[0:-1:int(np.round(sampling / eod_fr))]

    return middle_conv, eod_convolved_down, eod_convolved_up,eod_conv_downsampled

def find_periods(a_fe, time, eod_fr,a_fr,eod_fe,e):
    time_fish_r = time * 2 * np.pi * eod_fr
    eod_fish_r = a_fr * np.sin(time_fish_r)
    period_fish_r = time_fish_r[(time_fish_r <= np.mean(time_fish_r)+2 * np.pi) & (time_fish_r > np.mean(time_fish_r))]
    time_fish_e = time * 2 * np.pi * eod_fe[e]
    eod_fish_e = a_fe * np.sin(time_fish_e)
    period_fish_e = time_fish_e[(time_fish_e <= np.mean(time_fish_e)+ 2 * np.pi) & (time_fish_e > np.mean(time_fish_e))]
    return eod_fish_e,eod_fish_r,period_fish_r,period_fish_e

def rectify(eod_fish_r,eod_fe_chirp):
    eod_rec_up = eod_fish_r + eod_fe_chirp
    eod_rectified_down = eod_fish_r + eod_fe_chirp
    eod_rec_up[eod_rec_up < 0] = 0  # rectify
    eod_rectified_down[eod_rectified_down > 0] = 0  # rectify
    return eod_rectified_down, eod_rec_up

def integrate_chirp(a_fe,time,beat,phase_zero,size, sigma):
    I = ((np.pi ** 0.5) / 2) * sp.special.erf(time / sigma) - ((np.pi ** 0.5) / 2) * sp.special.erf(-np.inf)
    phase = time * 2 * np.pi *  beat+ 2 * np.pi * size * sigma * I + phase_zero
    eod_fe_chirp = a_fe * np.sin(phase)
    return eod_fe_chirp

def find_lm(eod_rec_up):
    x = signal.find_peaks(eod_rec_up)
    index_peaks = x[0]
    value_peaks = eod_rec_up[index_peaks]
    peaks_interp = np.interp(np.arange(0, len(eod_rec_up), 1), index_peaks, value_peaks)
    return index_peaks, value_peaks, peaks_interp

def find_times(left_c,right_c, sampling,deviation_s):
    for_conv = 5 * deviation_s
    time = np.arange(int(np.round(left_c))-1000, int(np.round(right_c))+1000, 1)
    time = time[(time >left_c) &(time < right_c)]
    time = time/sampling
    #time = np.arange(-for_conv+left_c,for_conv+right_c, 1 / sampling)
    cut = int(np.ceil(for_conv*sampling))
    time_cut = time[cut:-cut]
    return time, time_cut, cut


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 == True:
        eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
    else:
        eod_fe_chirp = eod_fish_e

    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
    period_length = len(period_fish_r)
    #embed()
    sampled = np.abs(eod_overlayed_chirp[0:len(eod_overlayed_chirp):period_length])
    time_new = time[0:len(eod_overlayed_chirp):period_length]
    maxima_values, maxima_index, maxima_interp = global_maxima(sampling, eod_fr, 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_new,sampled,  cut, 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 same_lengths(maxima_interp_b, maxima_interp_c, time_b, time_c, am_fish_c, am_fish_b, conv_c, conv_b, am_indirect_b,
                 am_indirect_c, peaks_interp_b, peaks_interp_c):
    diff_s = len(time_b) - len(time_c)
    diff_f = len(am_fish_c) - len(am_fish_b)
    if len(time_b) > len(time_c):
        time_b = time_b[0:len(time_c)]
        conv_b = conv_b[0:len(conv_c)]
        am_indirect_b = am_indirect_b[0:len(am_indirect_c)]
        peaks_interp_b = peaks_interp_b[0:len(peaks_interp_c)]
        maxima_interp_b = maxima_interp_b[0:len(maxima_interp_c)]
    if len(time_c) > len(time_b):
        time_c = time_c[0:len(time_b)]
        conv_c = conv_c[0:len(conv_b)]
        am_indirect_c = am_indirect_c[0:len(am_indirect_b)]
        peaks_interp_c = peaks_interp_c[0:len(peaks_interp_b)]
        maxima_interp_c = maxima_interp_c[0:len(maxima_interp_b)]
    if len(am_fish_b) > len(am_fish_c):
        am_fish_b = am_fish_b[0:len(am_fish_c)]
    if len(am_fish_c) > len(am_fish_b):
        am_fish_c = am_fish_c[0:len(am_fish_b)]
    return diff_s, diff_f, maxima_interp_b, maxima_interp_c, time_b, time_c, am_fish_c, am_fish_b, conv_c, conv_b, am_indirect_b, am_indirect_c, peaks_interp_b, peaks_interp_c

def power_func(restrict = [],a_fr = 1, a_fe = 0.2, eod_fr = 500, eod_fe = [500], win = 'w2', deviation_s = [0.0015], sigma = 50, sampling = 100000, deviation_ms = [1.5], beat_corr = [0], size = [120],phase_zero =  [0], delta_t = 0.014,deviation_dp = [180], show_figure = True, plot_dist = False, save = False,bef_c = -2,aft_c = -1):
    print('plot_eod')
    results = [[]]*10
    if len(restrict) >0:
        results = [[]]*len(restrict)
    for d in range(len(deviation_s)):
        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')
        for s in range(len(size)):
            for p in range(len(phase_zero)):
                for e in range(len(eod_fe)):
                    #embed()
                    left_b = [bef_c * sampling] * len(beat_corr)
                    right_b = [aft_c * sampling] * len(beat_corr)
                    time_fish,sampled, cut, cubed, 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_dp,
                                                                                                beat_corr, chirp = False)
                    thresholding = eod_overlayed_chirp*1
                    thresholding[eod_overlayed_chirp<0] = 0
                    squaring = eod_overlayed_chirp**2
                    fifth = (thresholding)**5
                    #embed()
                    rectified = np.abs(eod_overlayed_chirp)#'df ds','corr','corr ds',
                    name = ['df','conv','thresholding','rectified','squaring','cubed','fifth','global max','local max', 'samped threshold']#am_df_b,am_df_ds_b, am_corr_b, am_corr_ds_b,
                    var = [am_df_b,conv_b, thresholding, rectified, squaring,cubed ,fifth, maxima_interp_b,peaks_interp_b, sampled]
                    samp = [sampling,sampling,sampling, sampling, sampling,sampling,sampling,sampling, sampling, eod_fr]#eod_fr,sampling,eod_fr,
                    #['Hilbert', 'Thresholded', 'Rectified', 'Squaring', 'Threshold cubed']
                    #embed()
                    if len(restrict) == 1:
                        var = [np.array(var)[restrict[0]]]
                        name = [np.array(name)[restrict[0]]]
                        samp = [np.array(samp)[restrict[0]]]
                    elif len(restrict)>1:
                        var = np.array(var)[restrict]
                        name = np.array(name)[restrict]
                        samp = np.array(samp)[restrict]
                    #embed()
                    for i in range(len(results)):
                        if samp[i]>1000:
                            nfft = int(4096 * samp[0] / 40000 * 2)
                        else:
                            nfft = 4096
                        results[i],pp,f = power_parameters(results[i], var[i], samp[i],nfft,name[i],eod_fr)
                results = pd.DataFrame(results)
                if save:
                    results = pd.DataFrame(results)
                    try:
                        results.to_pickle('../data/power_simulation.pkl')
                        np.save('../data/Ramona/power_simulation.npy', results)
                    except:
                        a = 0
    return results

def power_parameters(results, var,samp,nfft,name, eod_fr,max_all = 'max_all',max = 'max',integral = 'amp',plot = False,title = '',freq_half = 'f',freq_whole = 'ff'):
    pp, f = ml.psd(var - np.mean(var), Fs=samp, NFFT=nfft,
                   noverlap=nfft / 2)  #
    # embed()
    # print(nfft * 5 < len(var[i]):


    if plot == True:
        plt.figure()
        plt.subplot(1, 2, 1)
        plt.plot(var)
        plt.title(name)
        plt.subplot(1, 2, 2)
        plt.plot(f, pp)
        # plt.xlim([0,2000])
        plt.show()
    # print(results)
    freq_step = np.abs(f[0] - f[1])
    if type(results) != dict:
        #embed()
        results = {}
        results['type'+title] = name
        # embed()
        results['result_frequency_whole'+title] = list([f[np.argmax(pp)]])
        results['result_frequency'+title] = list([f[np.argmax(pp[f < 0.5 * eod_fr])]])
        results['result_amplitude'+title] = list([np.sqrt(si.trapz(pp, f, np.abs(f[1] - f[0])))])
        results['result_amplitude_max'+title] = list([np.sqrt(np.max(pp[f < 0.5 * eod_fr]) * freq_step)])
        results['result_amplitude_max_whole'+title] = list([np.sqrt(np.max(pp) * freq_step)])
    else:
        #embed()
        results['result_frequency'+title].extend(list([f[np.argmax(pp[f < 0.5 * eod_fr])]]))
        results['result_frequency_whole'+title].extend(list([f[np.argmax(pp)]]))
        # embed()
        results['result_amplitude'+title].extend(list([np.sqrt(si.trapz(pp, f, np.abs(f[1] - f[0])))]))
        results['result_amplitude_max'+title].extend(list([np.sqrt(np.max(pp[f < 0.5 * eod_fr]) * freq_step)]))
        results['result_amplitude_max_whole'+title].extend(list([np.sqrt(np.max(pp) * freq_step)]))
    return results,pp,f

def distances_func(disc, 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):
    distance_metz_all = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_all = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_max_all = only_one_dict_str(d_ms, eod_fe)
    distance_wo_max_all = only_one_dict_str(d_ms, eod_fe)
    range_chirp_all = only_one_dict_str(d_ms, eod_fe)
    range_beat_all = only_one_dict_str(d_ms, eod_fe)
    corr_all = only_one_dict_str(d_ms, eod_fe)
    distance_metz_all_corr = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_all_in = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_max_all_in = only_one_dict_str(d_ms, eod_fe)
    distance_wo_max_all_in = only_one_dict_str(d_ms, eod_fe)
    range_chirp_all_in = only_one_dict_str(d_ms, eod_fe)
    range_beat_all_in = only_one_dict_str(d_ms, eod_fe)
    distance_metz_all_gm = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_all_gm = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_max_all_gm = only_one_dict_str(d_ms, eod_fe)
    distance_wo_max_all_gm = only_one_dict_str(d_ms, eod_fe)
    distance_metz_all_lm = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_all_lm = only_one_dict_str(d_ms, eod_fe)
    distance_wo_mean_max_all_lm = only_one_dict_str(d_ms, eod_fe)
    distance_wo_max_all_lm = only_one_dict_str(d_ms, eod_fe)
    range_chirp_all_lm = only_one_dict_str(d_ms, eod_fe)

    distance_metz_all_df_ds,distance_metz_all_corr_ds,range_beat_all_lm,range_chirp_all_gm ,range_beat_all_gm,distance_metz_all_df,distance_wo_mean_all_df,distance_wo_mean_max_all_df,distance_wo_max_all_df,range_chirp_all_df ,euc_all_in,euc_all_df,euc_all_gm,euc_all_lm,euc_all,range_beat_all_df   = (only_one_dict_str(d_ms, eod_fe) for i in range(16))
    period_distance_c,period_distance_b,period_shift = ({} for i in range(3))
    for d in range(len(deviation)):
        left_c, right_c, left_b, right_b, period_distance_c[d_ms[d]], period_distance_b[d_ms[d]], period_shift[d_ms[d]], 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)):
                    time_fish,sampled, cut, cube_b, 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_b = 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_fish,sampled,cut, cube_c, time_c, conv_c, am_corr_c, peaks_interp_c, maxima_interp_c,am_corr_ds_c, am_df_ds_c, am_df_c,eod_overlayed_chirp_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)


                    #print(len(time_c))
                    #print(len(time_b))
                    #diff_s[e], diff_f[e],maxima_interp_b, maxima_interp_c, time_b, time_c, am_df_c, am_df_b, conv_c,conv_b, am_corr_b, am_corr_c, peaks_interp_b, peaks_interp_c = same_lengths(
                    #    maxima_interp_b, maxima_interp_c, time_b, time_c, am_df_c, am_df_b,conv_c, conv_b, am_corr_b,
                    #    am_corr_c, peaks_interp_b, peaks_interp_c)
                    #time_b[e], conv_b[e], am_corr_b[e], peaks_interp_b[e], maxima_interp_b[e], am_df_b[e] = snip(left_c-shift[e], right_c-shift[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[e], conv_c[e], am_corr_c[e], peaks_interp_c[e], maxima_interp_c[e], am_df_c[e] = snip(left_c, right_c, e,
                    #                                                                            sampling, deviation_s,
                    #                                                                            d, eod_fr, a_fr, eod_fe,
                    # #                                                                           phase_zero, p, size, s,
                    #                                                                            sigma, a_fe, deviation,
                    #                                                                            beat_corr)

                    distance_metz_all_df_ds[d_ms[d]][e], _,_,_,_,_,_  = create_distance(np.array(am_df_ds_c), np.array(am_df_ds_b), None)
                    distance_metz_all_corr_ds[d_ms[d]][e], _,_,_,_,_,_ = create_distance(np.array(am_corr_ds_c), np.array(am_corr_ds_b), None)

                    distance_metz_all[d_ms[d]][e], distance_wo_mean_all[d_ms[d]][e],distance_wo_max_all[d_ms[d]][e],distance_wo_mean_max_all[d_ms[d]][e],range_chirp_all[d_ms[d]][e],range_beat_all[d_ms[d]][e],euc_all[d_ms[d]][e] = create_distance(np.array(conv_c),np.array(conv_b),None)
                    distance_metz_all_corr[d_ms[d]][e], distance_wo_mean_all_in[d_ms[d]][e],distance_wo_max_all_in[d_ms[d]][e],distance_wo_mean_max_all_in[d_ms[d]][e],range_chirp_all_in[d_ms[d]][e],range_beat_all_in[d_ms[d]][e],euc_all_in[d_ms[d]][e] = create_distance(np.array(am_corr_c),np.array(am_corr_b),None)
                    distance_metz_all_gm[d_ms[d]][e], distance_wo_mean_all_gm[d_ms[d]][e],distance_wo_max_all_gm[d_ms[d]][e],distance_wo_mean_max_all_gm[d_ms[d]][e],range_chirp_all_gm[d_ms[d]][e],range_beat_all_gm[d_ms[d]][e],euc_all_gm[d_ms[d]][e] = create_distance(np.array(maxima_interp_b),np.array(maxima_interp_c),None)
                    distance_metz_all_lm[d_ms[d]][e], distance_wo_mean_all_lm[d_ms[d]][e], distance_wo_max_all_lm[d_ms[d]][e], distance_wo_mean_max_all_lm[d_ms[d]][e], range_chirp_all_lm[d_ms[d]][e], range_beat_all_lm[d_ms[d]][e],euc_all_lm[d_ms[d]][e] = create_distance(
                        np.array(peaks_interp_c), np.array(peaks_interp_b),None)
                    distance_metz_all_df[d_ms[d]][e], distance_wo_mean_all_df[d_ms[d]][e], distance_wo_max_all_df[d_ms[d]][e], distance_wo_mean_max_all_df[d_ms[d]][e], range_chirp_all_df[d_ms[d]][e], range_beat_all_df[d_ms[d]][e],euc_all_df[d_ms[d]][e] = create_distance(
                        np.array(am_df_c), np.array(am_df_b),None)

                    #distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, range_chirp_all, range_beat_all, euc_all = create_distance(
                    #    np.array(conv_c), np.array(conv_b), None)
                    #distance_metz_all_corr, distance_wo_mean_all_in, distance_wo_max_all_in, distance_wo_mean_max_all_in, range_chirp_all_in, range_beat_all_in, euc_all_in = create_distance(
                    #    np.array(am_corr_c), np.array(am_corr_b), None)
                    #distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_max_all_gm, distance_wo_mean_max_all_gm, range_chirp_all_gm, range_beat_all_gm, euc_all_gm = create_distance(
                    #    np.array(maxima_interp_b), np.array(maxima_interp_c), None)
                    #distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_max_all_lm, distance_wo_mean_max_all_lm, range_chirp_all_lm, range_beat_all_lm, euc_all_lm = create_distance(
                    #    np.array(peaks_interp_c), np.array(peaks_interp_b), None)
                    #distance_metz_all_df, distance_wo_mean_all_df, distance_wo_max_all_df, distance_wo_mean_max_all_df, range_chirp_all_df, range_beat_all_df, euc_all_df = create_distance(
                    #    np.array(am_df_b), np.array(am_df_c), None)


            # correlation
            #corr = np.corrcoef(am_chirp, conv_chirp_corr)
            #corr_all.append(corr[1,0])
                if save:
                    np.save('../data/' + 'sim_dist3_' + str(
                            d_ms[d]) + 'ms' + 'beats(%1.2f' % (beats[0]) + '_%1.2f' % (beats[-1]) + '_%1.2f' % (abs(beats[1] - beats[0])) + ')' + save_n + '.npy',[distance_metz_all_df_ds,distance_metz_all_corr_ds,period_shift,period_distance_c, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_df, range_beat_all_df, range_chirp_all_df, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm,
                                         beat_corr, distance_metz_all_df, distance_wo_mean_all_df,distance_wo_mean_max_all_df, distance_wo_max_all_df, distance_metz_all_lm,
                                         distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm,
                                         distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm,
                                         distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling,
                                         d_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_corr, distance_wo_mean_all_in, eod_fe, eod_fr,distance_metz_all, distance_wo_mean_all, distance_wo_max_all,distance_wo_mean_max_all])
                    with open('../data/' + 'sim_dist3_' + str(
                            d_ms[d]) + 'ms' + 'beats(%1.2f' % (beats[0]) + '_%1.2f' % (beats[-1]) + '_%1.2f' % (abs(beats[1] - beats[0])) + ')' + save_n + '.pkl', 'wb') as f:  # Python 3: open(..., 'wb')
                            pickle.dump([distance_metz_all_df_ds,distance_metz_all_corr_ds,period_shift,period_distance_c, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_df, range_beat_all_df, range_chirp_all_df, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm,
                                        beat_corr, distance_metz_all_df, distance_wo_mean_all_df,
                                         distance_wo_mean_max_all_df, distance_wo_max_all_df, distance_metz_all_lm,
                                         distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm,
                                         distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm,
                                         distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling,
                                         d_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in,
                                         distance_metz_all_corr, distance_wo_mean_all_in, eod_fe, eod_fr,
                                         distance_metz_all, distance_wo_mean_all, distance_wo_max_all,
                                         distance_wo_mean_max_all], f)
    return  distance_metz_all_df_ds,distance_metz_all_corr_ds,period_shift,period_distance_c,period_distance_b, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_df, range_beat_all_df, range_chirp_all_df, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm \
               , beat_corr, distance_metz_all_df, distance_wo_mean_all_df, distance_wo_mean_max_all_df, distance_wo_max_all_df, \
           distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, \
           distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, \
           distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling, \
           d_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, \
           distance_metz_all_corr, distance_wo_mean_all_in, eod_fe, eod_fr, \
           distance_metz_all, distance_wo_mean_all, distance_wo_max_all, \
           distance_wo_mean_max_all

def plot_dist_interactiv3(bef_c, aft_c,d, sampling, size, phase_zero, sigma, beat_corr, delta_t, a_fr, a_fe, deviation_dp, deviation_s, d_ms, eod_fe, eod_fr, share = False, show = False):
    fig = plt.figure(figsize = [13,6])
    plt.title('convolved (red), indirect (magenta), loc max (green), glo max (black),sampled f (orange), period (purple)')
    gs0 = gridspec.GridSpec(3, 2)
    win = ['w1', 'w2', 'w3', 'w4']
    ax = {}
    gs00 = gridspec.GridSpecFromSubplotSpec(7, 2, subplot_spec=gs0[3])
    ax[0] = plt.subplot(gs00[4])
    p = 0
    s = 0
    for g in range(4):
        period_shift, period_distance_c, period_distance_b, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds, range_beat_all_ds, range_chirp_all_ds, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm, beat_corr_org, distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling, d_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, \
        distance_wo_mean_max_all = distances_func(bef_c, aft_c, win[g], deviation_s, sigma, sampling, d_ms, beat_corr, size,
                                                  phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation_dp,
                                                  show_figure=False, plot_dist=False, save=True)



        gs00 = gridspec.GridSpecFromSubplotSpec(7, 2, subplot_spec=gs0[g])
        #plt.title('w1')
        #fig = frame_subplots(6.5, 3, 'Beats [Hz]', '', 'Convolved: EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (d_ms[d]) + ') '+win)
        #two_ylabels(1,2,1,fig, 'distance (not normalised)', 'distance (normalised)')
        subplots = 14
        #ax = create_subpl(0,fig, subplots, 7, 2, share2 = share)
        plots = [distance_metz_all[d_ms[d]], euc_all[d_ms[d]],distance_metz_all_in[d_ms[d]],euc_all_in[d_ms[d]],distance_metz_all_lm[d_ms[d]], euc_all_lm[d_ms[d]], distance_metz_all_gm[d_ms[d]],euc_all_gm[d_ms[d]],  distance_metz_all_ds[d_ms[d]], euc_all_ds[d_ms[d]],period_distance_c[d_ms[d]],period_distance_b[d_ms[d]],period_shift[d_ms[d]],period_shift[d_ms[d]]]
        labels = ['conv mean','conv sum', 'indirect mean', 'indirect sum', 'lm mean', 'lm sum', 'gm mean', 'gm sum','ds mean', 'ds sum','period lengths','','period shift','']
        colors = ['red','red', 'magenta','magenta','green','green','black','black','orange','orange','darkblue','darkblue','lime','lime']
        I = {}


        for i in range(subplots):
            if share:
                ax[i] = plt.subplot(gs00[i],sharex=ax[0])
            else:
                ax[i] = plt.subplot(gs00[i])
            I[i] = ax[i].plot(eod_fe - eod_fr, plots[i], label=labels[i], color=colors[i], linestyle='-')
            ax[i].set_yticks([])
        axvline_everywhere(ax, subplots, eod_fe, eod_fr)
        ax[0].set_title(win[g] +' Mean ')
        ax[1].set_title(win[g]+' Sum (Euclidean) ')
        plots = [distance_wo_mean_all[d_ms[d]],distance_wo_mean_all_in[d_ms[d]],distance_wo_mean_all_lm[d_ms[d]],distance_wo_mean_all_gm[d_ms[d]],distance_wo_mean_all_ds[d_ms[d]]]
        for i in np.arange(0,subplots-4,2):
            ax[i].plot(eod_fe - eod_fr, plots[int(i/2)], color=colors[i], linestyle='--')

        plots1 = [distance_wo_max_all[d_ms[d]],distance_wo_max_all_in[d_ms[d]],distance_wo_max_all_lm[d_ms[d]],distance_wo_max_all_gm[d_ms[d]],distance_wo_max_all_ds[d_ms[d]]]
        plots2 = [distance_wo_mean_max_all[d_ms[d]], distance_wo_mean_max_all_in[d_ms[d]], distance_wo_mean_max_all_lm[d_ms[d]], distance_wo_mean_max_all_gm[d_ms[d]],
                  distance_wo_mean_max_all_ds[d_ms[d]]]
        for i in np.arange(0, subplots - 4, 2):
            ax[20+i] = ax[i].twinx()
            ax[20 +i].set_yticks([])
            ax[20+i].plot(eod_fe - eod_fr, plots1[int(i/2)],color=colors[i],linestyle='dashed')
            ax[20+i].plot(eod_fe - eod_fr, plots2[int(i/2)],linestyle='-.',
                  color=colors[i])
        plt.subplots_adjust(hspace = 0.15, wspace = 0.15)
        #fig.legend(I,labels=labels)
        #fig.legend([I[1],I[3],I[5],I[7],I[9],I[11]],labels = ['convolved','indirect','lm','gm','ds','period'])
        if win[g] == 'w1':
            ax[0].set_ylabel('convolved',color = 'red', rotation = 360, labelpad = 40)
            ax[2].set_ylabel('perfect', color='magenta', rotation = 360, labelpad = 40)
            ax[4].set_ylabel('local max', color='green', rotation = 360, labelpad = 40)
            ax[6].set_ylabel('global max', color='black', rotation = 360, labelpad = 40)
            ax[8].set_ylabel('sampled f', color='orange', rotation = 360, labelpad = 40)
            ax[10].set_ylabel('period dist', color='darkblue', rotation = 360, labelpad = 40)
            ax[12].set_ylabel('period shift', color='lime', rotation=360, labelpad=40)
            ax[8].set_xticks([])
            ax[10].set_xticks([])
        if win[g] == 'w3':
            ax[0].set_ylabel('w1: pd and ps, mult', rotation=360, labelpad=60)
            ax[2].set_ylabel('w2: same pd',  rotation=360, labelpad=60)
            ax[4].set_ylabel('w3: pd constant high f',  rotation=360, labelpad=60)
            ax[6].set_ylabel('w4, same pd and ps', rotation=360, labelpad=60)
        if win[g] == 'w2':
            ax[8].set_xticks([])
            ax[10].set_xticks([])

        #ax = plt.gca()
        #leg = ax.get_legend()
        #leg.legendHandles[0].set_color('red')
        #leg.legendHandles[1].set_color('darkblue')

        def onclick(event):
            e = np.argmin(np.abs((eod_fe - eod_fr) - event.xdata))
            onclick_effect(beat_corr, deviation_dp,delta_t,sampling, deviation_s,d,a_fe, eod_fr, a_fr, eod_fe, e,phase_zero,p, size,s, sigma)
            plt.show()

        if share == False:
            cid = fig.canvas.mpl_connect('button_press_event', onclick)
        #plt.savefig(
        #    fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
        #        size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (
        #              eod_fr) + ')_beat_%1.2f' % (
        #              beats[0]) + '_%1.2f' % (abs(beats[1] - beats[0])) +'_%1.2f' % (beats[-1]) + ')_deviation(%1.2f' % (
        #              d_ms[d]) + '_distances_all' + 'several' + '.png')

    g = 4
    gs00 = gridspec.GridSpecFromSubplotSpec(4, 2, subplot_spec=gs0[g])
    with open('C:/Users/alexi/Desktop/Masterarbeit/labrotation/data/consp4_1.5msbeats(-500.00_1995.00_5.00).pkl',
              'rb') as f:  # Python 3: open(..., 'rb')
        eod_fe, beat_corr, beats, integral_ds, integral_lm, integral_gm, integral_ind, integral_conv, max_values_gm, max_positions_gm, max_values_lm, max_positions_lm, max_values_conv, max_positions_conv, max_values_ind, max_positions_ind, max_values_ds, max_positions_ds = pickle.load(
            f)
    d = 1
    #two_ylabels(1,2,1,fig, 'Integral beneath the chirp Consp', '')
    #two_ylabels(1, 2, 2, fig, '', 'Conspicousness [range 0 to 1]')
    subplots = 8
    #ax = create_subpl(0.3,fig, subplots, 4, 2, share2 = share)
    plots = [integral_conv, max_values_conv,integral_lm,max_values_lm, integral_gm,  max_values_gm, integral_ind, max_values_ind]
    labels = ['Integral conv','Max Values conv','Integral lm','Max Values lm','Integral gm','Max Values gm','Integral ind','Max Values ind']
    colors = ['red','red', 'green','green','black','black','magenta','magenta','orange','orange','darkblue','darkblue']
    for i in range(subplots):
        if share:
            ax[i] = plt.subplot(gs00[i], sharex=ax[0])
        else:
            ax[i] = plt.subplot(gs00[i])
        ax[i].plot(eod_fe - eod_fr, plots[i], label=labels[i], color=colors[i], linestyle='-')
    axvline_everywhere(ax,subplots, eod_fe, eod_fr)
    #ax[0].set_ylabel('convolved',color = 'red', rotation)
    #ax[2].set_ylabel('loc max', color='green')
    #ax[4].set_ylabel('glo max', color='black')
    #ax[6].set_ylabel('perfect beat corr', color='magenta')
    ax[0].set_title('Integral')
    ax[1].set_title('Max Value')
    for i in np.arange(0,8,2):
        ax[i].set_ylim([0,int(20/3) * 3])
    for i in np.arange(1, 8, 2):
        ax[i].set_ylim([0, 1.1])
    #plt.legend(loc='upper left',bbox_to_anchor=(0, 2))
    period = period_func(beat_corr, sampling, 'stim')
    def onclick2(event):

        e = np.argmin(np.abs((eod_fe-eod_fr) - event.xdata))
        onclick_effect(beat_corr, deviation_dp,delta_t, sampling, deviation_s, d, a_fe, eod_fr, a_fr, eod_fe, e, phase_zero, p, size, s, sigma)
        rang = [0.5*sampling]
        time_fish,sampled, cut, cube_c, time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp = snip([-rang[0] * 0.5]*len(eod_fe), [rang[0] * 0.5]*len(eod_fe), e,e,sampling, deviation_s,
                                                                                   d, eod_fr, a_fr, eod_fe,
                                                                                   phase_zero, 0, size, 0,
                                                                                   sigma, a_fe, deviation,
                                                                                   beat_corr)
        B_conv, B_indirect,B_downsampled,B_lm,B_gm,time = snippets4(deviation_s,a_fr, sigma,period,time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp, eod_fe, beat_corr, a_fe, eod_fr, e, size, 0, phase_zero, 0, deviation, d, rang[0], delta_t, sampling, plot = False,show = False)

        #fig = plt.figure(figsize = [7,6])
        gaussian = gauss(10,sampling,0.02,1)
        pl = True
        sh = False
        if event.inaxes == ax[4] or event.inaxes == ax[5]:
            max_position, max_value, I,ConSpi2 = single_consp(1,eod_fe[e],eod_fr,'black',gaussian,B_gm, gm_chirp, sampling,plot = pl,show = sh)
        if event.inaxes == ax[2] or event.inaxes == ax[3]:
            max_position, max_value, I,ConSpi2  = single_consp(1,eod_fe[e],eod_fr,'green',gaussian, B_lm, lm_chirp, sampling,plot = pl,show = sh)
        if event.inaxes == ax[0] or event.inaxes == ax[1]:
            max_position, max_value, I,ConSpi2  = single_consp(1,eod_fe[e],eod_fr,'red',gaussian, B_conv,  conv_chirp, sampling,plot = pl,show = sh)
        if event.inaxes == ax[6] or event.inaxes == ax[7]:
            max_position, max_value, I,ConSpi2 = single_consp(1,eod_fe[e],eod_fr,'magenta',gaussian,B_indirect, am_indirect_chirp, sampling,plot = pl,show = sh)
        plt.xlabel('Time [ms]')
        plt.show()

    if share == False:
        cid = fig.canvas.mpl_connect('button_press_event', onclick)

    if show:
        plt.show()



def plot_dist_interactiv4(bef_c, aft_c,disc,d, sampling, size, phase_zero, sigma, beat_corr, delta_t, a_fr, a_fe, deviation_dp, deviation_s, d_ms, eod_fe, eod_fr, share = False, show = False):
    fig = plt.figure(figsize = [12,5])
    plt.title('convolved (red), indirect (magenta), loc max (green), glo max (black),sampled f (orange), period (purple)')
    gs0 = gridspec.GridSpec(2, 2)
    win = ['w1', 'w2', 'w3', 'w4']
    ax = {}
    gs00 = gridspec.GridSpecFromSubplotSpec(7, 2, subplot_spec=gs0[3])
    ax[0] = plt.subplot(gs00[4])
    p = 0
    s = 0
    for g in range(4):
        period_shift, period_distance_c, period_distance_b, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds, range_beat_all_ds, range_chirp_all_ds, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm, beat_corr_org, distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling, d_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, \
        distance_wo_mean_max_all = distances_func(disc, bef_c, aft_c, win[g], deviation_s, sigma, sampling, d_ms, beat_corr, size,
                                                  phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation_dp,
                                                  show_figure=False, plot_dist=False, save=True)



        gs00 = gridspec.GridSpecFromSubplotSpec(7, 2, subplot_spec=gs0[g])
        #plt.title('w1')
        #fig = frame_subplots(6.5, 3, 'Beats [Hz]', '', 'Convolved: EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (d_ms[d]) + ') '+win)
        #two_ylabels(1,2,1,fig, 'distance (not normalised)', 'distance (normalised)')
        subplots = 14
        #ax = create_subpl(0,fig, subplots, 7, 2, share2 = share)
        plots = [distance_metz_all[d_ms[d]], euc_all[d_ms[d]],distance_metz_all_in[d_ms[d]],euc_all_in[d_ms[d]],distance_metz_all_lm[d_ms[d]], euc_all_lm[d_ms[d]], distance_metz_all_gm[d_ms[d]],euc_all_gm[d_ms[d]],  distance_metz_all_ds[d_ms[d]], euc_all_ds[d_ms[d]],period_distance_c[d_ms[d]],period_distance_b[d_ms[d]],period_shift[d_ms[d]],period_shift[d_ms[d]]]
        labels = ['conv mean','conv sum', 'indirect mean', 'indirect sum', 'lm mean', 'lm sum', 'gm mean', 'gm sum','ds mean', 'ds sum','period lengths','','period shift','']
        colors = ['red','red', 'magenta','magenta','green','green','black','black','orange','orange','darkblue','darkblue','lime','lime']
        I = {}


        for i in range(subplots):
            if share:
                ax[i] = plt.subplot(gs00[i],sharex=ax[0])
            else:
                ax[i] = plt.subplot(gs00[i])
            I[i] = ax[i].plot(eod_fe - eod_fr, plots[i], label=labels[i], color=colors[i], linestyle='-')
            ax[i].set_yticks([])
        axvline_everywhere(ax, subplots, eod_fe, eod_fr)
        ax[0].set_title(win[g] +' Mean ')
        ax[1].set_title(win[g] +' Sum (Euclidean) ')
        plots = [distance_wo_mean_all[d_ms[d]],distance_wo_mean_all_in[d_ms[d]],distance_wo_mean_all_lm[d_ms[d]],distance_wo_mean_all_gm[d_ms[d]],distance_wo_mean_all_ds[d_ms[d]]]
        for i in np.arange(0,subplots-4,2):
            ax[i].plot(eod_fe - eod_fr, plots[int(i/2)], color=colors[i], linestyle='--')

        plots1 = [distance_wo_max_all[d_ms[d]],distance_wo_max_all_in[d_ms[d]],distance_wo_max_all_lm[d_ms[d]],distance_wo_max_all_gm[d_ms[d]],distance_wo_max_all_ds[d_ms[d]]]
        plots2 = [distance_wo_mean_max_all[d_ms[d]], distance_wo_mean_max_all_in[d_ms[d]], distance_wo_mean_max_all_lm[d_ms[d]], distance_wo_mean_max_all_gm[d_ms[d]],
                  distance_wo_mean_max_all_ds[d_ms[d]]]
        for i in np.arange(0, subplots - 4, 2):
            ax[20+i] = ax[i].twinx()
            ax[20 +i].set_yticks([])
            ax[20+i].plot(eod_fe - eod_fr, plots1[int(i/2)],color=colors[i],linestyle='dashed')
            ax[20+i].plot(eod_fe - eod_fr, plots2[int(i/2)],linestyle='-.',
                  color=colors[i])
        plt.subplots_adjust(hspace = 0.15, wspace = 0.15)
        #fig.legend(I,labels=labels)
        #fig.legend([I[1],I[3],I[5],I[7],I[9],I[11]],labels = ['convolved','indirect','lm','gm','ds','period'])
        if win[g] == 'w1':
            ax[0].set_ylabel('convolved',color = 'red', rotation = 360, labelpad = 40)
            ax[2].set_ylabel('perfect', color='magenta', rotation = 360, labelpad = 40)
            ax[4].set_ylabel('local max', color='green', rotation = 360, labelpad = 40)
            ax[6].set_ylabel('global max', color='black', rotation = 360, labelpad = 40)
            ax[8].set_ylabel('sampled f', color='orange', rotation = 360, labelpad = 40)
            ax[10].set_ylabel('period dist', color='darkblue', rotation = 360, labelpad = 40)
            ax[12].set_ylabel('period shift', color='lime', rotation=360, labelpad=40)
            ax[8].set_xticks([])
            ax[10].set_xticks([])
        if win[g] == 'w3':
            ax[0].set_ylabel('w1: pd and ps, mult', rotation=360, labelpad=60)
            ax[2].set_ylabel('w2: same pd',  rotation=360, labelpad=60)
            ax[4].set_ylabel('w3: pd constant high f',  rotation=360, labelpad=60)
            ax[6].set_ylabel('w4, same pd and ps', rotation=360, labelpad=60)
        if win[g] == 'w2':
            ax[8].set_xticks([])
            ax[10].set_xticks([])

        #ax = plt.gca()
        #leg = ax.get_legend()
        #leg.legendHandles[0].set_color('red')
        #leg.legendHandles[1].set_color('darkblue')

        def onclick(event):
            e = np.argmin(np.abs((eod_fe - eod_fr) - event.xdata))
            onclick_effect(beat_corr, deviation_dp,delta_t,sampling, deviation_s,d,a_fe, eod_fr, a_fr, eod_fe, e,phase_zero,p, size,s, sigma)
            plt.show()

        if share == False:
            cid = fig.canvas.mpl_connect('button_press_event', onclick)
        #plt.savefig(
        #    fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
        #        size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (
        #              eod_fr) + ')_beat_%1.2f' % (
        #              beats[0]) + '_%1.2f' % (abs(beats[1] - beats[0])) +'_%1.2f' % (beats[-1]) + ')_deviation(%1.2f' % (
        #              d_ms[d]) + '_distances_all' + 'several' + '.png')
    if show:
        plt.show()


def plot_dist_interactiv2_more(beat_corr, deviation_dp, delta_t, deviation_s, a_fe, a_fr, sigma,disc,win,period_shift,d,d_ms, period_distance, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds, beat_corr_org, distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s,  size, amplitude, phase_zero, beats, sampling, deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, share = False, show = False):
    fig = frame_subplots(12, 5, 'Beats [Hz]', '', 'Convolved: EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (d_ms[d]) + ') '+win)
    two_ylabels(110,1,2,1,fig, 'distance (not normalised)', 'distance (normalised)')
    subplots = 14
    ax = create_subpl(0,fig, subplots, 7, 2, share2 = share)
    plots = [distance_metz_all[d_ms[d]], euc_all[d_ms[d]],distance_metz_all_in[d_ms[d]],euc_all_in[d_ms[d]],distance_metz_all_lm[d_ms[d]], euc_all_lm[d_ms[d]], distance_metz_all_gm[d_ms[d]],euc_all_gm[d_ms[d]],  distance_metz_all_ds[d_ms[d]], euc_all_ds[d_ms[d]],period_distance[d_ms[d]],period_distance[d_ms[d]],period_shift[d_ms[d]],period_shift[d_ms[d]]]
    labels = ['conv mean','conv sum', 'indirect mean', 'indirect sum', 'lm mean', 'lm sum', 'gm mean', 'gm sum','ds mean', 'ds sum','period lengths','','period shift','']
    colors = ['red','red', 'magenta','magenta','green','green','black','black','orange','orange','darkblue','darkblue','lime','lime']
    I = {}
    for i in range(subplots):
        I[i] = ax[i].plot(eod_fe - eod_fr, plots[i], label=labels[i], color=colors[i], linestyle='-')
    axvline_everywhere(ax, subplots, eod_fe, eod_fr)
    ax[0].set_title('Mean (Mean instead sum)')
    ax[1].set_title('Sum (Euclidean distance)')

    plots = [distance_wo_mean_all[d_ms[d]],distance_wo_mean_all_in[d_ms[d]],distance_wo_mean_all_lm[d_ms[d]],distance_wo_mean_all_gm[d_ms[d]],distance_wo_mean_all_ds[d_ms[d]]]
    for i in np.arange(0,subplots-4,2):
        ax[i].plot(eod_fe - eod_fr, plots[int(i/2)], color=colors[i], linestyle='--')

    plots1 = [distance_wo_max_all[d_ms[d]],distance_wo_max_all_in[d_ms[d]],distance_wo_max_all_lm[d_ms[d]],distance_wo_max_all_gm[d_ms[d]],distance_wo_max_all_ds[d_ms[d]]]
    plots2 = [distance_wo_mean_max_all[d_ms[d]], distance_wo_mean_max_all_in[d_ms[d]], distance_wo_mean_max_all_lm[d_ms[d]], distance_wo_mean_max_all_gm[d_ms[d]],
              distance_wo_mean_max_all_ds[d_ms[d]]]
    for i in np.arange(0, subplots - 4, 2):
        ax[20+i] = ax[i].twinx()
        ax[20+i].plot(eod_fe - eod_fr, plots1[int(i/2)],color=colors[i],linestyle='dashed')
        ax[20+i].plot(eod_fe - eod_fr, plots2[int(i/2)],linestyle='-.',
              color=colors[i])
    #fig.legend(I,labels=labels)
    #fig.legend([I[1],I[3],I[5],I[7],I[9],I[11]],labels = ['convolved','indirect','lm','gm','ds','period'])
    ax[0].set_ylabel('convolved',color = 'red',rotation = 360, labelpad = 40)
    ax[2].set_ylabel('indirect', color='magenta',rotation = 360, labelpad = 40)
    ax[4].set_ylabel('loc max', color='green',rotation = 360, labelpad = 40)
    ax[6].set_ylabel('glo max', color='black',rotation = 360, labelpad = 40)
    ax[8].set_ylabel('sampled f', color='orange',rotation = 360, labelpad = 40)
    ax[10].set_ylabel('period distance', color='darkblue',rotation = 360, labelpad = 40)
    ax[12].set_ylabel('period shift', color='lime', rotation=360, labelpad=40)
    #ax = plt.gca()
    #leg = ax.get_legend()
    #leg.legendHandles[0].set_color('red')
    #leg.legendHandles[1].set_color('darkblue')


    def onclick(event):
        e = np.argmin(np.abs((eod_fe - eod_fr) - event.xdata))
        onclick_effect(beat_corr, deviation_dp,delta_t,sampling, deviation_s,d,a_fe, eod_fr, a_fr, eod_fe, e,phase_zero,p, size,s, sigma)
        plt.show()

    if share == False:
        cid = fig.canvas.mpl_connect('button_press_event', onclick)

    plt.savefig(
        fname= disc+'Masterarbeit/labrotation/results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
            size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (
                  eod_fr) + ')_beat_%1.2f' % (
                  beats[0]) + '_%1.2f' % (abs(beats[1] - beats[0])) +'_%1.2f' % (beats[-1]) + ')_deviation(%1.2f' % (
                  d_ms[d]) + '_distances_all' + 'several' + '.png')

    if show:
        plt.show()

def plot_dist_interactiv2(df_ds,corr_ds,beat_corr, deviation_dp, delta_t, deviation_s, a_fe, a_fr, sigma, disc, win, period_shift, d, d_ms, period_distance, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds, beat_corr_org, distance_metz_all_df, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s, size, amplitude, phase_zero, beats, sampling, deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_corr, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, share = False, show = False):
    #embed()
    fig = frame_subplots(7, 7, 'Beats [Hz]', '', 'Convolved: EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (d_ms[d]) + ') '+win)

    two_ylabels(130,1,2,1,fig, 'distance (normalised)', 'distance (normalised)')
    plots = [distance_metz_all[d_ms[d]],distance_metz_all_df[d_ms[d]], df_ds[d_ms[d]],distance_metz_all_corr[d_ms[d]],corr_ds[d_ms[d]],
             distance_metz_all_gm[d_ms[d]], distance_metz_all_lm[d_ms[d]],  period_distance[d_ms[d]],
             period_shift[d_ms[d]]]
    subplots = len(plots)

    ax = create_subpl(0,fig, subplots, subplots, 1, share2 = share)
    #plots = [distance_metz_all[d_ms[d]], euc_all[d_ms[d]], distance_metz_all_corr[d_ms[d]], euc_all_in[d_ms[d]], distance_metz_all_lm[d_ms[d]], euc_all_lm[d_ms[d]], distance_metz_all_gm[d_ms[d]], euc_all_gm[d_ms[d]], distance_metz_all_df[d_ms[d]], euc_all_ds[d_ms[d]], period_distance[d_ms[d]], period_distance[d_ms[d]], period_shift[d_ms[d]], period_shift[d_ms[d]]]
    labels = ['conv mean', 'df','df ds','corr mean','corr ds', 'gm mean','lm mean', 'period lengths','']
    colors = ['black','magenta','pink','orange','moccasin','red','green','darkblue','lime']
    I = {}


    for i in range(subplots):
        I[i] = ax[i].plot(eod_fe - eod_fr, plots[i], label=labels[i], color=colors[i], linestyle='-')
    axvline_everywhere(ax, subplots, eod_fe, eod_fr)
    ax[0].set_title('Mean (Mean instead sum)')
    #ax[1].set_title('Sum (Euclidean distance)')

    #plots = [distance_wo_mean_all[d_ms[d]],distance_wo_mean_all_in[d_ms[d]],distance_wo_mean_all_lm[d_ms[d]],distance_wo_mean_all_gm[d_ms[d]],distance_wo_mean_all_ds[d_ms[d]]]
    #for i in np.arange(0,subplots,1):
    #    ax[i].plot(eod_fe - eod_fr, plots[int(i/2)], color=colors[i], linestyle='--')

    #plots1 = [distance_wo_max_all[d_ms[d]],distance_wo_max_all_in[d_ms[d]],distance_wo_max_all_lm[d_ms[d]],distance_wo_max_all_gm[d_ms[d]],distance_wo_max_all_ds[d_ms[d]]]
    #plots2 = [distance_wo_mean_max_all[d_ms[d]], distance_wo_mean_max_all_in[d_ms[d]], distance_wo_mean_max_all_lm[d_ms[d]], distance_wo_mean_max_all_gm[d_ms[d]],
    #          distance_wo_mean_max_all_ds[d_ms[d]]]
    #for i in np.arange(0, subplots - 4, 2):
    #    ax[20+i] = ax[i].twinx()
    #    ax[20+i].plot(eod_fe - eod_fr, plots1[int(i/2)],color=colors[i],linestyle='dashed')
    #    ax[20+i].plot(eod_fe - eod_fr, plots2[int(i/2)],linestyle='-.',
    #          color=colors[i])
    #fig.legend(I,labels=labels)
    #fig.legend([I[1],I[3],I[5],I[7],I[9],I[11]],labels = ['convolved','indirect','lm','gm','ds','period'])
    ax[0].set_ylabel('convolved', color='black', rotation=360, labelpad=40)

    ax[1].set_ylabel('DF', color='magenta', rotation=360, labelpad=40)
    ax[2].set_ylabel('DF downsampled', color='pink', rotation=360, labelpad=40)

    ax[3].set_ylabel('beat corr', color='orange',rotation = 360, labelpad = 40)
    ax[4].set_ylabel('beat corr downsampled', color='moccasin', rotation=360, labelpad=40)

    ax[5].set_ylabel('global max', color='red', rotation=360, labelpad=40)
    ax[6].set_ylabel('local max', color='green', rotation=360, labelpad=40)

    ax[7].set_ylabel('segments length', color='darkblue',rotation = 360, labelpad = 40)
    ax[8].set_ylabel('timewindow', color='lime', rotation=360, labelpad=40)
    #ax = plt.gca()
    #leg = ax.get_legend()
    #leg.legendHandles[0].set_color('red')
    #leg.legendHandles[1].set_color('darkblue')
    plt.subplots_adjust(left = 0.3)

    def onclick(event):
        e = np.argmin(np.abs((eod_fe - eod_fr) - event.xdata))
        onclick_effect(beat_corr, deviation_dp,delta_t,sampling, deviation_s,d,a_fe, eod_fr, a_fr, eod_fe, e,phase_zero,p, size,s, sigma)
        plt.show()

    if share == False:
        cid = fig.canvas.mpl_connect('button_press_event', onclick)

    plt.savefig(
        fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
            size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (
                  eod_fr) + ')_beat_%1.2f' % (
                  beats[0]) + '_%1.2f' % (abs(beats[1] - beats[0])) +'_%1.2f' % (beats[-1]) + ')_deviation(%1.2f' % (
                  d_ms[d]) + '_distances_all' + 'several' + '.png')

    if show:
        plt.show()

def onclick_effect(beat_corr, deviation_dp, delta_t,sampling, deviation_s,d,a_fe, eod_fr, a_fr, eod_fe, e,phase_zero,p, size,s, sigma):
    #embed()
    left_c = -200 * delta_t*sampling
    right_c = 200 * delta_t*sampling
    time, time_cut, cut = find_times(left_c, right_c, 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)
    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(sampling, eod_fr, 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_dp[d], eod_rec_up,
                                                                         eod_rec_down)  # convolve
    am_corr = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p], size[s],
                                 sigma)  # indire# ct am calculation
    am_df = integrate_chirp(a_fe, time_cut, eod_fe[e]-eod_fr, phase_zero[p], size[s],
                                 sigma)  # indire# ct am calculation
    _, time_fish, cut_f = find_times(left_c, right_c, eod_fr, deviation_s[d])  # downsampled through fish EOD
    am_df_ds = integrate_chirp(a_fe, time_fish, eod_fe[e]-eod_fr, phase_zero[p], size[s],
                                 sigma)  # indire# ct am calculation
    am_corr_ds = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
    middle_am = int(len(am_corr_ds) / 2)
    plt.figure(figsize=[6.5, 5])

    #conv
    plt.plot(time_cut * 1000, eod_conv_up + 3.3, color='black',label='convolved', linewidth=2)

    plt.plot(time_cut * 1000, am_df + 3.4, color='magenta', label='DF', linewidth=2)
    plt.plot(
        time_fish * 1000,
        am_df_ds +2.8,
        color='pink',
        label='DF downsampled', linewidth=2)

    plt.plot(time_cut * 1000, am_corr + 2.2, color='orange', label='Beat corr', linewidth=2)
    plt.plot(
        time_fish * 1000,
        am_corr_ds + 1.6,
        color='moccasin',
        label='Beat corr downsampled', linewidth=2)


    plt.plot(time_cut * 1000, eod_conv_down - 1.3, color='black',  linewidth=2)

    #signal
    plt.plot(time_cut * 1000, eod_overlayed_chirp,
             label='EOD both fish',
             color='silver', linewidth=0.5)


    #AM maxima
    plt.scatter((index_peaks - 0.5 * len(eod_rec_up[cut:-cut])) / (sampling / 1000), value_peaks,
                color='green', s=10, linewidth=2)
    plt.plot((index_peaks - 0.5 * len(eod_rec_up[cut:-cut])) / (sampling / 1000),
             value_peaks,
             color='green', label='all maxima', linewidth=2)

    plt.scatter((maxima_index - 0.5 * len(eod_rec_up[cut:-cut])) / (sampling / 1000), maxima_values,
                color='red', s=10)
    plt.plot((maxima_index - 0.5 * len(eod_rec_up[cut:-cut])) / (sampling / 1000), maxima_values,
             color='red', label='maxima on lower frequency', linewidth=2)






    # beat corr

    # df

    plt.subplots_adjust(left = 0.4,hspace = 0.8)
    #embed()
    plt.title('Mult' + str(int(((eod_fe[e] - eod_fr)/eod_fe[e]+1)*100)/100))
    plt.legend(loc=(-0.75,0.6), ncol=1)
    plt.xlim([-200, 200])


def plot_cons_interactiv2(deviation_dp, deviation_s, a_fe, a_fr, sigma, delta_t,type,disc,integral_ind,p, s, d, size, amplitude, phase_zero, beats, sampling, deviation,deviation_ms, eod_fe, eod_fr, beat_corr, max_values_ds, max_positions_ds, integral_ds, max_values_lm, max_positions_lm, integral_lm, max_values_gm, max_positions_gm, integral_gm, max_values_ind, max_positions_ind, integral_corr, max_values_conv, max_positions_conv, integral_conv, share = False):
    fig = frame_subplots(13, 6, 'Beats [Hz]', '', 'Convolved: EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (deviation_ms[d]) + ')')
    two_ylabels(60,1,2,1,fig, 'Integral beneath the chirp Consp', '')
    two_ylabels(60,1, 2, 2, fig, '', 'Conspicousness [range 0 to 1]')
    subplots = 8
    ax = create_subpl(0.3,fig, subplots, 4, 2, share2 = share)
    plots = [integral_conv, max_values_conv,integral_lm,max_values_lm, integral_gm,  max_values_gm, integral_ind, max_values_ind]
    labels = ['Integral conv','Max Values conv','Integral lm','Max Values lm','Integral gm','Max Values gm','Integral ind','Max Values ind']
    colors = ['red','red', 'green','green','black','black','magenta','magenta','orange','orange','darkblue','darkblue']
    for i in range(subplots):
        ax[i].plot(eod_fe - eod_fr, plots[i], label=labels[i], color=colors[i], linestyle='-')
    axvline_everywhere(ax,subplots, eod_fe, eod_fr)
    ax[0].set_ylabel('convolved',color = 'red')
    ax[2].set_ylabel('loc max', color='green')
    ax[4].set_ylabel('glo max', color='black')
    ax[6].set_ylabel('perfect beat corr', color='magenta')
    ax[0].set_title('Integral')
    ax[1].set_title('Max Value')
    for i in np.arange(0,8,2):
        ax[i].set_ylim([0,int(20/3) * 3])
    for i in np.arange(1, 8, 2):
        ax[i].set_ylim([0, 1.1])
    #plt.legend(loc='upper left',bbox_to_anchor=(0, 2))
    period = period_func(beat_corr, sampling,'stim')
    def onclick(event):

        e = np.argmin(np.abs((eod_fe-eod_fr) - event.xdata))
        onclick_effect(beat_corr, deviation_dp,delta_t, sampling, deviation_s, d, a_fe, eod_fr, a_fr, eod_fe, e, phase_zero, p, size, s, sigma)
        rang = [0.5*sampling]
        time_fish,sampled, cut, cube_c, time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp = snip([-rang[0] * 0.5]*len(eod_fe), [rang[0] * 0.5]*len(eod_fe), e,e,sampling, deviation_s,
                                                                                   d, eod_fr, a_fr, eod_fe,
                                                                                   phase_zero, 0, size, 0,
                                                                                   sigma, a_fe, deviation,
                                                                                   beat_corr)
        B_conv, B_indirect,B_downsampled,B_lm,B_gm,time = snippets4(deviation_s,a_fr, sigma,period,time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp, eod_fe, beat_corr, a_fe, eod_fr, e, size, 0, phase_zero, 0, deviation, d, rang[0], delta_t, sampling, plot = False,show = False)

        #fig = plt.figure(figsize = [7,6])
        gaussian = gauss(10,sampling,0.02,1)
        pl = True
        sh = False
        if event.inaxes == ax[4] or event.inaxes == ax[5]:
            max_position, max_value, I,ConSpi2 = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'black',gaussian,B_gm, gm_chirp, sampling,'stim',plot = pl,show = sh)
        if event.inaxes == ax[2] or event.inaxes == ax[3]:
            max_position, max_value, I,ConSpi2  = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'green',gaussian, B_lm, lm_chirp, sampling,'stim',plot = pl,show = sh)
        if event.inaxes == ax[0] or event.inaxes == ax[1]:
            max_position, max_value, I,ConSpi2  = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'red',gaussian, B_conv,  conv_chirp, sampling,'stim',plot = pl,show = sh)
        if event.inaxes == ax[6] or event.inaxes == ax[7]:
            max_position, max_value, I,ConSpi2 = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'magenta',gaussian,B_indirect, am_indirect_chirp, sampling,'stim',plot = pl,show = sh)
        plt.xlabel('Time [ms]')
        plt.show()

    if share == False:
        cid = fig.canvas.mpl_connect('button_press_event', onclick)

    plt.savefig(

        fname='../results/' + 'simulations/' +str(type[0]) +'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
            size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (
                  eod_fr) + ')_beat_%1.2f' % (
                  beats[0]) + '_%1.2f' % (abs(beats[1] - beats[0])) +'_%1.2f' % (beats[-1])+ ')_deviation(%1.2f' % (
                  deviation_ms[d])  + '_consp' + 'several' + '.png')
    plt.show()

def plot_dist_interactiv(distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, maxima_index_all, maxima_all, am_indirect_all, am_fish_all, time_fish_all, eod_convolved_down_all, eod_convolved_up_all, value_peaks_all, p, s, d, size, amplitude, phase_zero, beats, sampling, eod_rectified_up_all, index_peaks_all, time, eod_overlayed_chirp_all, deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, interact2 = False):
    fig = plt.figure(figsize=[13, 6])
    #plt.axes(frameon=False)
    ax = fig.add_subplot(1, 1, 1)
    plt.xticks([])
    plt.yticks([])
    ax.set_xticks([])
    ax.set_yticks([])
    ax0 = ax.twinx()
    ax0.set_xticks([])
    ax0.set_yticks([])
    ax1 = fig.add_subplot(5, 1, 1)
    ax2 = fig.add_subplot(5, 1, 4)
    ax3 = fig.add_subplot(5, 1, 3)
    ax4 = fig.add_subplot(5, 1, 2)
    ax5 = fig.add_subplot(5, 1, 5)
    ax1.plot(eod_fe - eod_fr, distance_metz_all, label='Metzen', color='red',linestyle=':')
    #ax1.set_xticks([])
    ax1.plot(eod_fe - eod_fr, distance_wo_mean_all, label='Metzen ohne Mean', color='red',linestyle='-')
    ax11 = ax1.twinx()
    ax11.plot(eod_fe - eod_fr, distance_wo_max_all, label='Metzen ohne Normierung', color='red',linestyle='--')
    ax11.plot(eod_fe - eod_fr, distance_wo_mean_max_all, label='Metzen ohne Mean und Normierung',linestyle='-.',
              color='red')
    ax2.plot(eod_fe - eod_fr, distance_metz_all_in, label='Metzen',linestyle=':',color ='magenta')
    #ax2.set_xticks([])
    ax2.plot(eod_fe - eod_fr, distance_wo_mean_all_in, label='Metzen ohne Mean',linestyle='-',color ='magenta')
    ax22 = ax2.twinx()
    ax22.plot(eod_fe - eod_fr, distance_wo_max_all_in, label='Metzen ohne Normierung',linestyle='--',color ='magenta')
    ax22.plot(eod_fe - eod_fr, distance_wo_mean_max_all_in, label='Metzen ohne Mean und Normierung',linestyle='-.',color ='magenta')
    ax3.plot(eod_fe - eod_fr, distance_metz_all_lm, label='Metzen', color='green',linestyle=':')
    ax3.plot(eod_fe - eod_fr, distance_wo_mean_all_lm, label='Metzen ohne Mean', color='green',linestyle='-')
    ax33 = ax3.twinx()
    ax33.plot(eod_fe - eod_fr, distance_wo_max_all_lm, label='Metzen ohne Normierung', color='green',linestyle='--')
    ax33.plot(eod_fe - eod_fr, distance_wo_mean_max_all_lm, label='Metzen ohne Mean und Normierung',linestyle='-.',
              color='green')
    ax4.plot(eod_fe - eod_fr, distance_metz_all_gm, label='Metzen', color='black',linestyle=':')
    ax4.plot(eod_fe - eod_fr, distance_wo_mean_all_gm, label='Metzen ohne Mean', color='black',linestyle='-')
    ax44 = ax4.twinx()
    ax44.plot(eod_fe - eod_fr, distance_wo_max_all_gm, label='Metzen ohne Normierung', color='black',linestyle='--')
    ax44.plot(eod_fe - eod_fr, distance_wo_mean_max_all_gm, label='Metzen ohne Mean und Normierung',
              color='black',linestyle='-.')
    ax5.plot(eod_fe - eod_fr, distance_metz_all_ds, label='Metzen', color='orange',linestyle=':')
    ax5.plot(eod_fe - eod_fr, distance_wo_mean_all_ds, label='Metzen ohne Mean', color='orange',linestyle='-')
    ax55 = ax5.twinx()
    ax55.plot(eod_fe - eod_fr, distance_wo_max_all_ds, label='Metzen ohne Normierung', color='orange',linestyle='--')
    ax55.plot(eod_fe - eod_fr, distance_wo_mean_max_all_ds, label='Metzen ohne Mean und Normierung',
              color='orange',linestyle='-.')
    ax1.axvline(x=-eod_fr, color='black', linewidth=1, linestyle='-')
    for i in range(int(np.max(eod_fe) / np.max(eod_fr))):
        ax1.axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')
        ax2.axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')
        ax3.axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')
        ax4.axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')
        ax5.axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')

    #ax2.axvline(x=-eod_fr, color='black', linewidth=1, linestyle='-')
    #for i in range(int(np.max(eod_fe) / np.max(eod_fr))):

    # ax1.legend(loc = 'upper left',bbox_to_anchor=(1.15, 1.6))
    # ax11.legend(loc='upper right', bbox_to_anchor=(1.15, 1.6))
    ax1.legend(loc='upper left',bbox_to_anchor=(0, 2))
    ax11.legend(loc='upper right',bbox_to_anchor=(1.0, 2))
    #ax2.legend(loc='upper left')
    #ax22.legend(loc='upper right',bbox_to_anchor=(1, 1))
    #ax3.legend(loc='upper left')
    #ax33.legend(loc='upper right',bbox_to_anchor=(1, 1))
    #ax4.legend(loc='upper left')
    #ax44.legend(loc='upper right',bbox_to_anchor=(1, 1))

    ax1.set_title('Convolved: EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (deviation_ms[d]) + ')')
    ax2.set_title('Beat corr')
    ax3.set_title('Lobal maxima')
    ax4.set_title('Global maxima')
    ax5.set_title('Down sampled')
    #ax1.set_xlabel('beat [hz]')
    #ax1.set_ylabel('distance')
    #ax2.set_xlabel('beat [hz]')
    #ax2.set_ylabel('distance')
    ax5.set_xlabel('beat [hz]')
    #ax3.set_ylabel('distance')
    ax.set_ylabel('distance (not normalised)', labelpad = 40)
    ax0.set_ylabel('distance (normalised)', rotation = 270,labelpad = 40)

    plt.subplots_adjust(hspace = 0.65)

    def onclick(event):
        i = np.argmin(np.abs((eod_fe-eod_fr) - event.xdata))
        plt.figure(figsize = [7,6])
        plt.plot(time * 1000, eod_overlayed_chirp_all[i][3 * deviation_dp[d]:-3 * deviation_dp[d]],
                 label='EOD both fish',
                 color='blue', linewidth=1)
        #ax3.title('beat:' + str(beat) + 'Hz ' + 'rf:' + str(eod_fr) + ' ef:' + str(eod_fe[e]))
        plt.scatter((index_peaks_all[i] - 0.5 * len(eod_rectified_up_all[i])) / (sampling / 1000), value_peaks_all[i],
                    color='green', s=10, linewidth=2)
        plt.plot((index_peaks_all[i] - 0.5 * len(eod_rectified_up_all[i])) / (sampling / 1000),
                 value_peaks_all[i],
                 color='green', label='all maxima', linewidth=2)
        plt.plot(time * 1000, eod_convolved_up_all[i], color='red', linewidth=2)
        plt.plot(time * 1000, eod_convolved_down_all[i], color='red', label='convolved', linewidth=2)

        plt.plot(
            time_fish_all[i][int(3 * deviation_dp[d] / factor):int(-3 * deviation_dp[d] / factor)] * 1000,
            am_fish_all[i][int(3 * deviation_dp[d] / factor):int(-3 * deviation_dp[d] / factor)] - 0.3,
            color='orange',
            label='indirect am -  downgesampled', linewidth=2)
        plt.scatter((maxima_index_all[i] - 0.5 * len(eod_rectified_up_all[i])) / (sampling / 1000), maxima_all[i],
                   color='black', s=10)
        plt.plot((maxima_index_all[i] - 0.5 * len(eod_rectified_up_all[i])) / (sampling / 1000), maxima_all[i],
                color='black', label='all maxima', linewidth=2)
        plt.plot(time * 1000, am_indirect_all[i] + 0.3, color='magenta', label='beat corr', linewidth=2)
        plt.title('beat'+str(eod_fe[i]-eod_fr))
        plt.legend(loc = 'lower center', ncol = 3, bbox_to_anchor = (0.5,1.05))
        plt.xlim([-400,400])
        plt.show()

    if interact2:
        cid = fig.canvas.mpl_connect('button_press_event', onclick)

    plt.savefig(
        fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
            size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (
                  eod_fr) + ')_beat_%1.2f' % (
                  beats[0]) + '_%1.2f' % (abs(beats[1] - beats[0])) +'_%1.2f' % (beats[-1])+ ')_deviation(%1.2f' % (
                  deviation_ms[d])  + '_distances_all' + 'several' + '.png')
    plt.show()

def plot_eod_phase(beat_corr, minus_bef ,plus_bef,row,col,a_fe, delta_t, a_fr, size,beat_corr_org,s,pp, amplitude, phase_zero, deviation_ms, eod_fe,deviation, eod_fr, sampling,factor, fs = [6.2992, 4], shift_phase = 0):
    beats = eod_fe - eod_fr
    #fig = plt.figure(figsize=[6.2992/3, 4])
    fig = plt.figure(figsize=fs)
    ax = {}
    counter = -1
    ax[0] = plt.subplot(row, col, 1)  # int(np.ceil(np.sqrt(len(eod_fe))))
    for p,p_nr in zip(pp,range(len(pp))):
        d = 0
        for e in range(len(eod_fe)):
            counter += 1
            #time, time_cut = find_times(time_range, sampling, deviation[d], 1)
            left_c = minus_bef * delta_t * sampling
            right_c = plus_bef * delta_t * sampling
            time, time_cut,cut = find_times(left_c, right_c, sampling, deviation[d]/(1000*sampling))
            #embed()
            time_fish_both = time * 2 * np.pi * (eod_fr-eod_fe[e])
            eod_fish_both = 0.05 * np.sin(time_fish_both)
            eod_fish_e, eod_fish_r, period_fish_r, period_fish_e = find_periods(a_fe,time, eod_fr, a_fr, eod_fe, e)
            eod_fish_both = integrate_chirp(a_fe, time, eod_fe[e]-eod_fr, phase_zero[p]+shift_phase, size[s], sigma)
            eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
            eod_overlayed_chirp = eod_fish_r + eod_fe_chirp
            eod_rectified_down, eod_recitified_up = rectify(eod_fish_r, eod_fe_chirp)  # rectify
            maxima_values, maxima_index, maxima_interp = global_maxima(sampling, eod_fr, period_fish_e, period_fish_r,
                                                                       eod_recitified_up)  # global maxima
            index_peaks, value_peaks, peaks_interp = find_lm(eod_recitified_up)  # local maxima
            middle_conv, eod_convolved_down, eod_convolved_up,eod_conv_downsampled = conv(eod_fr,sampling, cut,deviation[d], eod_recitified_up, eod_rectified_down)  # convolve
            left_c = -200 * delta_t * sampling
            right_c = 200 * delta_t * sampling
            _, time_fish,_ = find_times(left_c,right_c, eod_fr, deviation[d])  # downsampled through fish EOD
            am_fish = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
            middle_am = int(len(am_fish) / 2)
            print(beat_corr[e])
            am_corr = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p]+shift_phase, size[s],
                                          sigma)  # indirect am calculation
            if plus_bef < 0:
                green_true = False
            else:
                green_true = True
                plt.axvline(x=-7.5, color='black', linestyle='dotted', linewidth=1)
                plt.axvline(x=7.5, color='black', linestyle='dotted', linewidth=1)
            ax[counter] = pl_eods(row, col,eod_fish_both,cut, beats,beat_corr, peaks_interp, maxima_interp, maxima_index, maxima_values,fig, ax[0],counter, e, time_cut ,am_corr,eod_fe,eod_overlayed_chirp,deviation,d,eod_fr,sampling,value_peaks, time_fish,am_fish,factor,eod_convolved_down,index_peaks,eod_convolved_up,eod_recitified_up,green_true = green_true,add = 'simple')
            ax[counter].title.set_text('Beat:' + str(eod_fe[e] - eod_fr) + 'Hz'+' phase'+str(int(phase_zero[p]*100)/100))
            #embed()
            xticks = 'off'
            yticks = 'off'
            plot_pos = col * row - col + 1
            ax[counter].set_xlabel('Time [ms]', labelpad=5)
            ax[counter].set_ylabel('[mv]', labelpad=5)

            #if counter == 2:
            #    ax[counter].set_xlabel('Time [ms]', labelpad=5)
            #    ax[counter].set_ylabel('[mv]', labelpad=5)
            #else:
            #    if xticks == 'off':
            #        ax[counter].set_xticks([])
            #    if yticks == 'off':
            #        ax[counter].set_yticks([])
            #lower_left_label(e+1, col, row, 'Time [ms]', '[mv]',xticks = 'off',yticks = 'off',)
        #plt.xlabel('Time [ms]')
    plt.legend(loc=(-1, -1.3),ncol = 2)
    plt.subplots_adjust(wspace=.2, hspace=.7,bottom = 0.3)
    #embed()
    plt.savefig(
        fname='../results/' + 'thesis/ '+'phasepic'+str(p)+str(eod_fe[e])+'.png')
    plt.show()


def plot_eod(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):
    beats = eod_fe - eod_fr
    for d in range(len(deviation)):
        fig = plt.figure(figsize=fs)
        ax = {}

        #shift_phase = 0

        ax[d] = fig.add_subplot(row,col, 1)#int(np.ceil(np.sqrt(len(eod_fe))))
        for e in range(len(eod_fe)):
            #time, time_cut = find_times(time_range, sampling, deviation[d], 1)
            left_c = minus_bef  * sampling
            right_c = plus_bef  * sampling
            time, time_cut,cut = find_times(left_c, right_c, sampling, deviation[d]/(1000*sampling))
            #embed()
            time_fish_both = time * 2 * np.pi * (eod_fr-eod_fe[e])
            eod_fish_both = 0.05 * np.sin(time_fish_both)
            eod_fish_e, eod_fish_r, period_fish_r, period_fish_e = find_periods(a_fe, time, eod_fr, a_fr, eod_fe, e)
            time_test = period_fish_e * 2 * np.pi * eod_fe[e]
            eod_test = a_fe * np.sin(time_test)
            #embed()
            eod_fish_both = integrate_chirp(a_fe, time, eod_fe[e]-eod_fr, phase_zero[p]+shift_phase, size[s], sigma)
            eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
            eod_overlayed_chirp = eod_fish_r + eod_fe_chirp
            eod_rectified_down, eod_recitified_up = rectify(eod_fish_r, eod_fe_chirp)  # rectify


            maxima_values, maxima_index, maxima_interp = global_maxima(sampling, eod_fr, period_fish_e, period_fish_r,
                                                                       eod_recitified_up)  # global maxima
            index_peaks, value_peaks, peaks_interp = find_lm(eod_recitified_up)  # local maxima
            middle_conv, eod_convolved_down, eod_convolved_up,eod_conv_downsampled = conv(eod_fr,sampling, cut,deviation[d], eod_recitified_up, eod_rectified_down)  # convolve
            left_c = -200 * delta_t * sampling
            right_c = 200 * delta_t * sampling
            _, time_fish,_ = find_times(left_c,right_c, eod_fr, deviation[d])  # downsampled through fish EOD
            am_fish = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
            middle_am = int(len(am_fish) / 2)
            print(beat_corr[e])
            am_corr = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p]+shift_phase, size[s],
                                          sigma)  # indirect am calculation
            if plus_bef < 0:
                green_true = False
            else:
                green_true = True
                plt.axvline(x=-7.5, color='black', linestyle='dotted', linewidth=1)
                plt.axvline(x=7.5, color='black', linestyle='dotted', linewidth=1)
            ax[e] = pl_eods(row, col,eod_fish_both,cut, beats,beat_corr, peaks_interp, maxima_interp, maxima_index, maxima_values,fig, ax[d],e, e, time_cut ,am_corr,eod_fe,eod_overlayed_chirp,deviation,d,eod_fr,sampling,value_peaks, time_fish,am_fish,factor,eod_convolved_down,index_peaks,eod_convolved_up,eod_recitified_up,add = False,green_true = green_true)
            #embed()
            xticks = 'off'
            yticks = 'off'
            plot_pos = col * row - col + 1
            #if e+1 == plot_pos:
            #    ax[e].set_xlabel('Time [ms]', labelpad=5)
            if e+1 in np.arange(row*col-col+1,row*col+1,1):
                ax[e].set_xlabel('Time [ms]', labelpad=5)
            if e + 1 in np.arange(1,row*col+1,col):
                ax[e].set_ylabel('[mv]', labelpad=5)
            ax[e].set_yticks([])
            #else:
            #    if xticks == 'off':
            #        ax[e].set_xticks([])
            #    if yticks == 'off':
            #        ax[e].set_yticks([])
            #lower_left_label(e+1, col, row, 'Time [ms]', '[mv]',xticks = 'off',yticks = 'off',)
        #plt.xlabel('Time [ms]')
        plt.legend(loc=(-1, -0.7),ncol = 2)
        plt.subplots_adjust(wspace=.2, hspace=.4,bottom = 0.20)
        #embed()
        #scalebar = ScaleBar(50,'ms')
        #plt.gca().add_artist(scalebar)
        #embed()
        # plt.xlim([-200,200])
        plt.savefig(
            fname='../highbeats_pdf/cell_simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
                size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')_deviation(%1.2f' % (
                      deviation_ms[d]) + 'beats(%1.2f' % (beats[0]) + '_%1.2f' % (beats[-1]) + '_%1.2f' % (
                      abs(beats[1] - beats[0])) + ')_eods_several_range'+str(minus_bef)+'_'+str(plus_bef)+  '.png')
        plt.show()
        embed()


def lower_left_label(number, ncol, nrow, x, y,xticks = 'off',yticks = 'off',labelpad = 5):
    plot_pos = ncol * nrow - ncol+1
    if number ==  plot_pos:
        plt.xlabel(x,labelpad = labelpad)
        plt.ylabel(y,labelpad = labelpad)
    else:
        if xticks == 'off':
            plt.xticks([])
        if yticks == 'off':
            plt.yticks([])

def plot_am(a_fe, delta_t, a_fr, size,beat_corr_org,s,p, amplitude, phase_zero, deviation_ms, eod_fe,deviation, eod_fr, sampling,factor):
    beats = eod_fe - eod_fr
    for d in range(len(deviation)):
        fig = plt.figure(figsize=[12, 6])
        ax = {}
        ax = fig.add_subplot(int(np.ceil(np.sqrt(len(eod_fe)))), int(np.ceil(np.sqrt(len(eod_fe)))), 1)
        for e in range(len(eod_fe)):
            time, time_cut = find_times(time_range, sampling, deviation[d], 1)
            eod_fish_e, eod_fish_r, period_fish_r, period_fish_e = find_periods(a_fe, time, eod_fr, a_fr, eod_fe, e)
            eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
            eod_overlayed_chirp = eod_fish_r + eod_fe_chirp
            eod_rectified_down, eod_rec_up = rectify(eod_fish_r, eod_fe_chirp)  # rectify
            maxima_values, maxima_index, maxima_interp = global_maxima(sampling, eod_fr, period_fish_e, period_fish_r,
                                                                       eod_rec_up)  # global maxima
            index_peaks, value_peaks, peaks_interp = find_lm(eod_rec_up)  # local maxima
            middle_conv, eod_conv_down, eod_conv_up,eod_conv_downsampled = conv(eod_fr,sampling, deviation[d], eod_rec_up, eod_rectified_down)  # convolve
            _, time_fish = find_times(time_range, eod_fr, deviation[d], 50)  # downsampled through fish EOD
            am_fish = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
            middle_am = int(len(am_fish) / 2)
            am_indirect = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p], size[s],
                                          sigma)  # indirect am calculation
            period_dist,period_dist,am_beat, maxima_interp_beat, peaks_interp_beat, maxima_interp_chirp, peaks_interp_chirp, am_indirect_beat, \
            time_fish,sampled,cube, am_indirect_chirp, time_beat, conv_beat, conv_chirp_corr, am_chirp, time_chirp, conv_chirp = snippets3(
                1.2, peaks_interp, maxima_interp, am_indirect, am_fish, d, period, middle_am, eod_convolved_up,
                middle_conv,
                delta_t, sampling, deviation, time, eod_conv_downsampled)

            sav, ax = plot_amp(eod_fr,time, ax, e, fig, 'indirect', eod_fe, time_chirp, time_beat, am_indirect_chirp,
                               am_indirect_beat, am_indirect)
        plt.xlabel('Time [ms]')
        plt.legend(loc=(1, 0))
        plt.subplots_adjust(wspace=.2, hspace=.7)
        # plt.xlim([-200,200])
        plt.savefig(
            fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
                size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')_deviation(%1.2f' % (
                      deviation_ms[d]) + 'beats(%1.2f' % (beats[0]) + '_%1.2f' % (beats[-1]) + '_%1.2f' % (
                      abs(beats[1] - beats[0])) + ')_am.png')
        plt.show()

def pl_eods(row,col,eod_fish_both, cut, beats, beat_corr, peaks_interp, maxima_interp, maxima_index, maxima, fig, ax, i, e, time, am_corr, eod_fe, eod_overlayed_chirp, deviation, d, eod_fr, sampling, value_peaks, time_fish, am_fish, factor, eod_convolved_down, index_peaks, eod_convolved_up, eod_rectified_up, add = False, lw = 1, add1 = False,share = False,green_true = True):

    if share == True:
        ax = fig.add_subplot(row, col, i + 1, sharex=ax,
                             sharey=ax)
    else:
        ax = fig.add_subplot(row, col, i + 1)
    #embed()
    #ax = plt.subplot(1,1,1)
    ax.plot(time * 1000, eod_overlayed_chirp[cut:-cut],
            label='EOD both fish',
            color='silver', linewidth = lw)
    #ax.plot(time * 1000, eod_overlayed_chirp[3 * deviation[d]:-3 * deviation[d]],
    #        label='EOD both fish',
    #        color='blue', linewidth=2)
    if add == True:
        ax.title.set_text('Beat:' + str(eod_fe[e]-eod_fr) + 'Hz ' + 'rf:' + str(eod_fr) + ' ef:' + str(eod_fe[e]))
    elif add == 'simple':
        ax.title.set_text('Beat:' + str(eod_fe[e] - eod_fr) + 'Hz')
    else:
        ax.title.set_text('Beat:' + str(eod_fe[e]-eod_fr) + 'Hz, Mult:'+str(int(((eod_fe[e]-eod_fr)/eod_fr+1)*100)/100))
    if eod_fe[e]-eod_fr <0:
        green_true = True
        black_true = False
    else:
        green_true = False
        black_true = True
    if green_true == True:
        #ax.scatter((index_peaks - 0.5 * len(eod_rectified_up[cut:-cut])) / (sampling / 1000), value_peaks,
        #           color='green', s=10)
        ax.plot(time * 1000, peaks_interp[cut:-cut], color='red', label = 'all maxima',linewidth = lw)
        #embed()
    #ax.plot((index_peaks - 0.5 * len(eod_rectified_up)) / (sampling / 1000), value_peaks,
    #        color='green', label='all maxima')
    #plt.axvline(x=-7.5, color='black', linestyle='dotted', linewidth=1)
    #plt.axvline(x=7.5, color='black', linestyle='dotted', linewidth=1)
    if add == True:
        ax.plot(time * 1000, eod_convolved_up, color='red',linewidth = lw)
        ax.plot(time * 1000, eod_convolved_down, color='red', label='convolved',linewidth = lw)
    if add1 == True:
        ax.plot(time_fish[int(3 * deviation[d] / factor):int(-3 * deviation[d] / factor)] * 1000,
                am_fish[int(3 * deviation[d] / factor):int(-3 * deviation[d] / factor)]+ 0.4, color='purple',
                label='indirect am -  downgesampled',linewidth = lw)


    ax.plot(time * 1000, am_corr +1.55, color='orange', label='EOD adjusted beat', linewidth = lw)
    if add1 == True:
        ax.scatter((maxima_index - 0.5 * len(eod_rectified_up)) / (sampling / 1000), maxima,
                   color='red', s=10)
    if black_true == True:
        ax.plot(time*1000, maxima_interp[cut:-cut], color='red', label= 'AM',linewidth = lw)#[int(3 * deviation[d]):int(-3 * deviation[d])]
    ax.plot(time*1000,eod_fish_both[cut:-cut]+ 2.10,color='magenta',label= 'Difference frequency', linewidth = lw)

    #ax.set_xlim([-(2.5*1/beat_corr)*1000,(2.5*1/beat_corr)*1000])
    #ax.plot((maxima_index  - 0.5 * len(eod_rectified_up)) / (sampling / 1000), maxima,
    #        color='black', label='one maxima pro lower frequency')
    #ax.set_ylabel('Time [ms]')
    #ax.set_ylabel('[mV]')
    return ax

def plot_amp(eod_fr,time ,ax, e, fig, sav, eod_fe,time_chirp,time_beat,conv_chirp,conv_beat, eod_convolved_up):
    ax = fig.add_subplot(int(np.ceil(np.sqrt(len(eod_fe)))), int(np.ceil(np.sqrt(len(eod_fe)))), e + 1,sharex = ax,sharey = ax)
    ax.plot(time * 1000, eod_convolved_up,color='red', linewidth=2,label = 'beat')
    ax.plot(time_chirp * 1000, conv_chirp, color='blue', linewidth=2, label = 'chirp')
    ax.plot(time_beat * 1000, conv_beat, color='green', linewidth=2,label = 'beat')
    ax.set_title('beat' + str(eod_fe[e]-eod_fr) + 'Hz')
    #ax.set_ylim([0.2, 1])
    return sav,ax


def plot_all_distances(d,s,p,show_figure,beats,amplitude,size,phase_zero,corr_all,deviation_ms,eod_fe,eod_fr,distance_metz_all,distance_wo_mean_all,distance_wo_max_all,distance_wo_mean_max_all,range_chirp_all,range_beat_all, version):
    # plt.title('deviation(%1.2f' % (deviation[d]) + ')_beat(_%1.2f' % (beats[0]) + '_%1.2f' % (
    # beats[-1]) + '_%1.2f' % (abs(beats[1] - beats[0])) + ')recieverEOD(_%1.2f' % (eod_fr) + ')')
    fig = plt.figure(figsize=[10, 10])
    ax1 = fig.add_subplot(3, 1, 1)
    ax2 = fig.add_subplot(3, 1, 2, sharex = ax1)
    ax3 = fig.add_subplot(3, 1, 3, sharex = ax1)
    ax1.plot(eod_fe - eod_fr, distance_metz_all, label='Metzen', color='orange')
    ax1.plot(eod_fe - eod_fr, distance_wo_mean_all, label='Metzen ohne Mean', color='blue')
    ax11 = ax1.twinx()
    ax11.plot(eod_fe - eod_fr, distance_wo_max_all, label='Metzen ohne Normierung', color='red')
    ax11.plot(eod_fe - eod_fr, distance_wo_mean_max_all, label='Metzen ohne Mean und Normierung', color='green')
    ax22 = ax2.twinx()
    ax2.plot(eod_fe - eod_fr, range_chirp_all, label='range chirps',linestyle='-')
    ax2.plot(eod_fe - eod_fr, range_beat_all, label='range beat',linestyle='--')
    ax1.axvline(x=-eod_fr, color='black', linewidth=1, linestyle='-')
    for i in range(int(np.max(eod_fe) / np.max(eod_fr))):
        ax1.axvline(x=eod_fr * i, color='black', linewidth=1, linestyle='-')
    # ax1.legend(loc = 'upper left',bbox_to_anchor=(1.15, 1.6))
    # ax11.legend(loc='upper right', bbox_to_anchor=(1.15, 1.6))
    ax1.legend(loc='upper left')
    ax11.legend(loc='upper right')

    ax1.set_title('EOD reciever ' + str(eod_fr) + '_deviation(%1.2f' % (deviation_ms[d]) + ')')
    ax2.legend()
    ax1.set_xlabel('beat [hz]')
    ax1.set_ylabel('distance')
    ax2.set_xlabel('Amplitude')
    ax2.set_ylabel('distance')
    ax3.plot(eod_fe - eod_fr, corr_all)
    ax3.set_xlabel('beat [hz]')
    ax3.set_ylabel('Correlation')

    plt.savefig(
        fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
        size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')recieverEOD(_%1.2f' % (eod_fr) + '_%1.2f' % (
              beats[-1]) + '_%1.2f' % (abs(beats[1] - beats[0])) + ')_deviation(%1.2f' % (
              deviation_ms[d]) + ')_beat(_%1.2f' % (beats[0]) + '_distances_all'+version+'.png')
    if show_figure:
        plt.show()
    else:
        plt.close()



def plot_distances(eod_fe,p, s,deviation_ms,d, distance,size,phase_zero,amplitude, show_figure = False):
    for e in range(len(eod_fe)):
        time, time_cut = find_times(time_range, sampling, deviation_dp[d], 1)
        eod_fish_e, eod_fish_r, period_fish_r, period_fish_e = find_periods(a_fe,time, eod_fr, a_fr, eod_fe, e)
        eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
        eod_overlayed_chirp = eod_fish_r + eod_fe_chirp
        eod_rectified_down, eod_rec_up = rectify(eod_fish_r, eod_fe_chirp)  # rectify
        maxima_values, maxima_index, maxima_interp = global_maxima(sampling, eod_fr, period_fish_e, period_fish_r,
                                                                   eod_rec_up)  # global maxima
        index_peaks, value_peaks, peaks_interp = find_lm(eod_rec_up)  # local maxima
        middle_conv, eod_conv_down, eod_conv_up,eod_conv_downsampled = conv(eod_fr,sampling, deviation_dp[d], eod_rec_up, eod_rectified_down)  # convolve
        _, time_fish = find_times(time_range, eod_fr, deviation_dp[d], 50)  # downsampled through fish EOD
        am_fish = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
        middle_am = int(len(am_fish) / 2)
        am_indirect = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p], size[s],
                                      sigma)  # indirect am calculation
        period_dist,am_beat, maxima_interp_beat, peaks_interp_beat, maxima_interp_chirp, peaks_interp_chirp, am_indirect_beat, \
        cube,am_indirect_chirp, time_beat, conv_beat, conv_chirp_corr, am_chirp, time_chirp, conv_chirp = snippets3(
            1.2, peaks_interp, maxima_interp, am_indirect, am_fish, d, period, middle_am, eod_convolved_up,
            middle_conv,
            delta_t, sampling, deviation_dp, time, eod_conv_downsampled)
        plt.plot(distance)
        plt.plot(conv_chirp)
        plt.plot(conv_beat)
        if show_figure:
            plt.show()
        else:
            plt.close()
        plt.savefig(
            fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
                size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')_deviation(%1.2f' % (deviation_ms[d])+')_beat(%5.2f' % (
                      beat) + ')_absolute_distance.png')

def plot_snippets(time_chirp, time_beat, time_range,p, sigma, s,delta_t,period,beat,sampling, deviation,d, time, eod_convolved_up,conv_chirp,conv_beat,size,phase_zero,amplitude,show_figure = False):
    plt.plot(time* 1000, eod_convolved_up, label='Amplitude Modulation',color='red', linewidth=2)
    plt.plot(time_chirp * 1000, conv_chirp, label='Amplitude Modulation', color='blue', linewidth=2)
    plt.plot(time_beat * 1000, conv_beat, label='Amplitude Modulation', color='green', linewidth=2)
    if show_figure:
        plt.show()
    else:
        plt.close()
    plt.savefig(
        fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
            size[s]) + ')_phase(%1.2f' % (phase_zero[p]) +  ')_deviation(%1.2f' % (deviation_ms[d])+'_beat(%5.2f' % (
                  beat) + ')_distances_snippets_single.png')




def plot_subplots(eod_fe,sides,p,s,deviation,d, deviation_ms,sampling,size,phase_zero,amplitude, show_figure = False, chirp = True):
    # frequency calculation
    for e in range(len(eod_fe)):
        time, time_cut = find_times(time_range, sampling, deviation[d], 1)
        eod_fish_e, eod_fish_r, period_fish_r, period_fish_e = find_periods(a_fe,time, eod_fr, a_fr, eod_fe, e)
        if chirp == True:
            eod_fe_chirp = integrate_chirp(a_fe, time, eod_fe[e], phase_zero[p], size[s], sigma)
        else:
            eod_fe_chirp = eod_fish_e
        eod_overlayed_chirp = eod_fish_r + eod_fe_chirp
        eod_rectified_down, eod_rec_up = rectify(eod_fish_r, eod_fe_chirp)  # rectify
        maxima_values, maxima_index, maxima_interp = global_maxima(sampling, eod_fr, period_fish_e, period_fish_r,
                                                                   eod_rec_up)  # global maxima
        index_peaks, value_peaks, peaks_interp = find_lm(eod_rec_up)  # local maxima
        middle_conv, eod_conv_down, eod_conv_up,eod_conv_downsampled = conv(eod_fr,sampling, deviation[d], eod_rec_up, eod_rectified_down)  # convolve
        _, time_fish = find_times(time_range, eod_fr, deviation[d], 50)  # downsampled through fish EOD
        am_fish = integrate_chirp(a_fe, time_fish, beat_corr[e], phase_zero[p], size[s], sigma)
        middle_am = int(len(am_fish) / 2)
        am_indirect = integrate_chirp(a_fe, time_cut, beat_corr[e], phase_zero[p], size[s],
                                      sigma)  # indirect am calculation

        chirp = size[s] * np.exp((-(time ** 2) / (2 * sigma ** 2)))
        frequency = eod_fr - eod_fe[e] + chirp

        fig = plt.figure(figsize=[10, 6])
        ax = {}
        ax[0] = fig.add_subplot(3, 2, 1)

        for i in range(6):
            ax[i] = fig.add_subplot(3, 2, i + 1, sharex = ax[0])
            ax[i].axvline(x=-sides, color='black', linewidth=1, linestyle='-')
            ax[i].axvline(x=sides, color='black', linewidth=1, linestyle='-')
        ax[0].plot(time * 1000, eod_fish_r[3 * deviation[d]:-3 * deviation[d]])
        ax[0].title.set_text('EOD reciever fish')
        ax[2].plot(time * 1000, eod_fe_chirp[3 * deviation[d]:-3 * deviation[d]])
        ax[2].title.set_text('EOD and chirp of emmiting fish')
        ax[4].plot(time * 1000, frequency)
        ax[4].title.set_text('Frequency')
        ax[1].plot(time* 1000, eod_overlayed_chirp[3 * deviation[d]:-3 * deviation[d]])
        # ax[3].plot(time[3 * deviation[d]:-3 * deviation[d]] * 1000,eod_rectified[3 * deviation[d]:-3 * deviation[d]])

        ax[1].title.set_text('EOD both fish + chirp')
        (ax[1]).plot(time * 1000, eod_convolved_up, color='red')
        (ax[1]).plot(time * 1000, eod_convolved_down, color='red')
        corr_up = np.array([eod_convolved_up - np.min(eod_convolved_up) - np.mean(eod_convolved_up)])
        # calculate AM directly
        time_sam = np.arange(-100 * delta_t, 100 * delta_t, 1 / 500*50)
        I1 = []
        for i in range(len(time_sam)):
            y1 = ((np.pi ** 0.5) / 2) * math.erf(time_sam[i] / sigma) - ((np.pi ** 0.5) / 2) * math.erf(-np.inf)
            I1.append(y1)
        phase = 2 * np.pi * beat * time_sam + 2 * np.pi * size[s] * sigma * np.array(I1) + phase_zero[p]

        # am calculation
        am = amplitude * np.cos(phase)
        (ax[3]).plot(time_sam[3 * deviation[d]:-3 * deviation[d]] * 1000, am[
                                                                      3 * deviation[d]: - 3 * deviation[d]], color='green')

        #ax33 = ax[3].twinx()
        #ax33.twinx().plot(time * 1000, eod_convolved_up, color='red')

        corr_down = np.array([eod_convolved_down - np.min(eod_convolved_down) - np.mean(eod_convolved_down)])



        # (ax[5]).plot(time * 1000, eod_convolved_down[
        #                                                            length_convolved + 3 * deviation[
        #                                                                d]:-length_convolved - 3 *deviation[d]]-np.min(eod_convolved_down)-np.mean(eod_convolved_down), color = 'red')

        #time = np.arange(-5 * delta_t, 5 * delta_t, 1 / sampling)
        # calculate AM directly
        time_fish = np.arange(-100 * delta_t, 100 * delta_t, 1 / 500)
        I1 = []
        for i in range(len(time_fish)):
            y1 = ((np.pi ** 0.5) / 2) * math.erf(time_fish[i] / sigma) - ((np.pi ** 0.5) / 2) * math.erf(-np.inf)
            I1.append(y1)
        phase = 2 * np.pi * beat * time_fish + 2 * np.pi * size[s] * sigma * np.array(I1) + phase_zero[p]

        # am calculation
        am = amplitude * np.cos(phase)
        #FIXME
        #(ax[5]).twinx().plot(time * 1000, eod_convolved_down,
        #                     color='red')
        (ax[5]).plot(time_fish[int(np.round((3 * deviation[d]) / sampling)):int(
            np.round((-3 * deviation[d]) / sampling))] * 1000, am[int(
            np.round((3 * deviation[d]) / sampling)):int(np.round((-3 * deviation[d]) / sampling))],
                     color='green')
        ax[1].set_xlim([-60,60])
        ax[2].set_xlim([-60, 60])
        ax[3].set_xlim([-60, 60])
        ax[4].set_xlim([-60, 60])
        ax[5].set_xlim([-60, 60])
        # eod_convolved2 = (eod_convolved[length_convolved + int(3.395 * deviation[d]): len(
        #    eod_convolved) - length_convolved - int(3.395 * deviation[d])])
        plt.subplots_adjust(wspace=.2, hspace=.7)
        plt.savefig(fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (
        size[s]) + ')_phase(%1.2f' % (phase_zero[p]) +  ')_deviation(%1.2f' % (deviation_ms[d])+'_beat(%5.2f' % (beat) + ')_subplots.png')

        if show_figure:
            plt.show()
        else:
            plt.close()



def stimulus_harmonics(size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation, show_figure = False):
    sampling = 40000
    time = np.arange(-5 * delta_t, 5 * delta_t, 1 / sampling)
    for s in range(len(size)):
        for p in range(len(phase_zero)):
            for d in range(len(deviation)):
                save_div = deviation[d] / (sampling / 1000)
                if deviation[d] * 5 % 2:
                    points = deviation[d] * 5
                else:
                    points = deviation[d] * 5 -1
                gaussian = signal.gaussian(points, std=deviation[d], sym=True)
                gaussian_normalised = (gaussian * 100) / np.sum(gaussian)
                length_convolved = int(len(gaussian_normalised) / 2)
                #cut_convolved = np.arange(length_convolved + (int(3.395 * deviation[d])),
                #                          (len(time) - length_convolved - int(3.395 * deviation[d])), 1)
                fig = plt.figure(figsize=[5, 4])
                ax = {}
                # plot only the maxima
                for i in range(1):
                    ax[i] = fig.add_subplot(1, 1, i + 1)
                    ax[i].axvline(x=-7.5, color='black', linewidth=1, linestyle='-')
                    ax[i].axvline(x=7.5, color='black', linewidth=1, linestyle='-')
                for e in range(len(eod_fe)):
                    eod_fish_r = a_fr*np.sin(time * 2 * np.pi * eod_fr)
                    sigma = delta_t / math.sqrt((2 * math.log(10)))  # width of the chirp
                    I1 = []
                    for i in range(len(time)):
                        y1 = ((np.pi ** 0.5) / 2) * math.erf(time[i] / sigma) - ((np.pi ** 0.5) / 2) * math.erf(-np.inf)
                        I1.append(y1)
                    phase1 = time * 2 * np.pi * eod_fe[e] + 2 * np.pi * size[s] * sigma * np.array(I1) + phase_zero[p]
                    eod_fe_chirp = a_fe*np.sin(phase1)
                    eod_overlayed_chirp = eod_fish_r + eod_fe_chirp
                    eod_rectified = eod_fish_r + eod_fe_chirp
                    amplitude = a_fe
                    beat = eod_fe[e] - eod_fr
                    # rectify and gaussian
                    eod_rectified[np.where(eod_overlayed_chirp < 0)] = 0  # rectify
                    eod_convolved = np.convolve(gaussian_normalised, eod_rectified)
                    # deviation = 60  # 40 = (1/1000)/(1 / 40000): 40 datapoints for 1 milisecond: 40*5 = 200
                    # frequency calculation
                    chirp = size[s] * np.exp((-(time**2)/(2*sigma**2)))
                    frequency = beat + chirp
                    #ax[0].plot(time[(length_convolved + (int(3.395 * deviation[d]))):(len(time) - length_convolved - int(3.395 * deviation[d]))]* 1000, eod_convolved[(length_convolved + (int(3.395 * deviation[d]))):(len(time) - length_convolved - int(3.395 * deviation[d]))],label=str(beat))
                    ax[0].plot(time[3*deviation[d]:-3*deviation[d]] * 1000, eod_convolved[length_convolved+3*deviation[d]:-length_convolved-3*deviation[d]], label=str(beat))
                    ax[0].title.set_text('deviation('+str(save_div)+')')
                    ax[0].legend(loc = 'lower left')
                plt.subplots_adjust(wspace=.2, hspace=.7)
                plt.savefig(fname='../results/' + 'simulations/' + 'amplitude(%1.2f' % (amplitude) + ')_size(%5.2f' % (size[s]) + ')_phase(%1.2f' % (phase_zero[p]) + ')_beat(%5.2f' % (beat) + ')_AM_stimulus_illustration3.png')
                if show_figure:
                    plt.show()
                else:
                    plt.close()

def snippets4(deviation_s,a_fr, sigma, period, time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp, eod_fe, beat_corr, a_fe,eod_fr, e,size, ss, phase_zero, p, deviation, d, window, delta_t, sampling, plot = False,show = False):
    shifts_b = 100
    B_conv = [[]]*shifts_b
    B_indirect = [[]]*shifts_b
    B_downsampled = [[]]*shifts_b
    B_lm = [[]]*shifts_b
    time = [[]] * shifts_b
    B_gm = [[]] * shifts_b
    ran = np.linspace(0, 3* period[e], shifts_b)
    shift = np.round(window * 2 + ran)
    for s in range(shifts_b):
        time_fish,sampled,cut, time[s],  B_conv[s] , B_indirect[s], B_lm[s], B_gm[s] , B_downsampled[s] = snip(-window*0.5-shift, window*0.5-shift, e,s,
                                                                                         sampling, deviation_s,
                                                                                         d, eod_fr, a_fr, eod_fe,
                                                                                         phase_zero, p, size, ss,
                                                                                         sigma, a_fe, deviation,
                                                                                         beat_corr)
    if plot:
        fig = frame_subplots(6, 6, 'Time', 'AM', str(eod_fe[e]-eod_fr)+'Hz')
        ax = create_subpl(0.2, fig, 5, 3, 2, share2=True)
        ax[0].set_title('convovled')
        ax[0].plot(np.transpose(B_conv[0:100:20]))
        ax[1].set_title('global maximum')
        ax[1].plot(np.transpose(B_gm[0:100:20]))
        ax[2].set_title('local maximum')
        ax[2].plot(np.transpose(B_lm[0:100:20]))
        ax[3].set_title('downsampled')
        ax[3].plot(np.transpose(B_downsampled[0:100:20]))
        ax[4].set_title('indirect')
        ax[4].plot(np.transpose(B_indirect[0:100:20]))
    if show:
        plt.show()
    return B_conv, B_indirect,B_downsampled,B_lm,B_gm,time




def conspicousy(disc,deviation_s,sigma, type,eod_fe, beats, deviation_ms, d, sampling, beat_corr, a_fe, delta_t, eod_fr, a_fr, size, phase_zero, deviation):
    max_values_gm = [[]]*len(eod_fe)
    max_positions_gm = [[]]*len(eod_fe)
    max_values_lm = [[]]*len(eod_fe)
    max_positions_lm = [[]]*len(eod_fe)
    max_values_conv = [[]]*len(eod_fe)
    max_positions_conv = [[]]*len(eod_fe)
    max_values_ind = [[]]*len(eod_fe)
    max_positions_ind = [[]]*len(eod_fe)
    max_values_ds = [[]]*len(eod_fe)
    max_positions_ds = [[]]*len(eod_fe)
    integral_ds = [[]]*len(eod_fe)
    integral_lm= [[]]*len(eod_fe)
    integral_gm= [[]]*len(eod_fe)
    integral_ind= [[]]*len(eod_fe)
    integral_conv= [[]]*len(eod_fe)
    rang = [0.4*sampling]
    max = [[]]*len(eod_fe)
    min = [[]] * len(eod_fe)
    period = period_func(beat_corr,sampling, 'stim')
    for e in range(len(eod_fe)):
        time_fish,sampled,cut, time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp = snip([-rang[0] * 0.5]*len(eod_fe), [rang[0] * 0.5]*len(eod_fe), e,e,sampling, deviation_s,
                                                                                   d, eod_fr, a_fr, eod_fe,
                                                                                   phase_zero, 0, size, 0,
                                                                                   sigma, a_fe, deviation,
                                                                                   beat_corr)

        B_conv, B_indirect,B_downsampled,B_lm,B_gm, time = snippets4(deviation_s,a_fr, sigma,period,time_c, conv_chirp, am_indirect_chirp, lm_chirp, gm_chirp, am_chirp, eod_fe, beat_corr, a_fe, eod_fr, e, size, 0, phase_zero, 0, deviation, d, rang[0], delta_t, sampling, plot = False,show = False)



        gaussian = gauss(10,sampling,0.02,1)
        max[e] = np.max(np.concatenate(time))
        min[e] = np.min(np.concatenate(time))
        pl = False
        sh = False
        max_positions_gm[e], max_values_gm[e], integral_gm[e],ConSpi2 = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'black',gaussian,B_gm, gm_chirp, sampling,'stim',plot = pl, show = sh)
        max_positions_lm[e], max_values_lm[e], integral_lm[e],ConSpi2 = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'green',gaussian,B_lm, lm_chirp, sampling,'stim',plot = pl,show = sh)
        max_positions_conv[e], max_values_conv[e], integral_conv[e],ConSpi2 = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'red',gaussian,B_conv,  conv_chirp, sampling,'stim',plot = pl,show = sh)
        max_positions_ind[e], max_values_ind[e], integral_ind[e],ConSpi2 = single_consp(type,0,1,eod_fe[e]-eod_fr,eod_fr,'magenta',gaussian,B_indirect, am_indirect_chirp,  sampling,'stim',plot = pl,show = sh)
        #max_positions_ds[e], max_values_ds[e], integral_ds[e],ConSpi2 = single_consp(gaussian,B_downsampled, am_chirp, eod_fr)

    with open(disc +'Masterarbeit/labrotation/data/' + 'consp4_' +str(type[0])+ str(
                deviation_ms[d]) + 'ms' + 'beats(%1.2f' % (beats[0]) + '_%1.2f' % (beats[-1]) + '_%1.2f' % (
                  abs(beats[1] - beats[0])) + ')' + '.pkl', 'wb') as f:  # Python 3: open(..., 'wb')
                pickle.dump(
                [eod_fe,beat_corr,beats,integral_ds, integral_lm, integral_gm, integral_ind, integral_conv, max_values_gm, max_positions_gm, max_values_lm, max_positions_lm, max_values_conv, max_positions_conv, max_values_ind, max_positions_ind, max_values_ds, max_positions_ds], f)

    return eod_fe,beat_corr,beats,integral_ds,  integral_lm,  integral_gm, integral_ind,  integral_conv, max_values_gm,max_positions_gm,max_values_lm,max_positions_lm,max_values_conv,max_positions_conv,max_values_ind,max_positions_ind,max_values_ds,max_positions_ds

def print_len(time):
    for i in range(len(time)):
        print(len(time[i]))






def test_periods():
    plt.plot(period)
    plt.plot(period <= delta_t * skretch)
    plt.plot(period_shift)
    plt.plot(np.ceil(delta_t * skretch / period) * period)
    plt.plot(np.ceil(delta_t * skretch / period))


    a = np.ones(100)
    b = np.ones(100)*2
    C1 = []
    C2 = []
    C3 = []
    C4 = []
    euc = [[]]*len(a)
    mean = [[]] * len(a)
    for c in range(len(a)):
        g = np.sqrt(np.mean((a[0:c]-np.mean(a[0:c]) - b[0:c] -np.mean(b[0:c]) ) ** 2))
        C1.append(g)
        d = np.sqrt(np.mean((a[0:c]  - b[0:c]) ** 2))
        C2.append(d)
        k = np.sqrt(np.mean((a[0:c]-np.mean(a[0:c]) - b[0:c] -np.mean(b[0:c]) )** 2))
        #C3.append(k)
        o = np.mean(np.sqrt(((a[0:c]-np.mean(a[0:c]) - b[0:c] -np.mean(b[0:c]) ) ** 2)))
        k = np.sqrt(np.sum((a[0:c]-np.mean(a[0:c]) - b[0:c] -np.mean(b[0:c]) )** 2)/len(a[0:c]))
        C3.append(k)
        o = np.sum(np.sqrt(((a[0:c]-np.mean(a[0:c]) - b[0:c] -np.mean(b[0:c]) ) ** 2)))/len(a[0:c])
        C4.append(o)
        mean[c] = dm_wo_mean2_max = mean_euc_d(a[0:c],  b[0:c], None)
        euc[c] = d_euclidean = euc_d(a[0:c],  b[0:c], None)



def find_beats(start,end,step,eod_fr):
    eod_fe = np.arange(start, end, step)
    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]
    return  eod_fe, beat_corr, beats

def find_dev(x, sampling):
    deviation_ms = np.array(x)
    deviation_s = deviation_ms/1000
    deviation_dp = sampling*deviation_s
    deviation_dp = list(map(int, deviation_dp))
    return deviation_ms, deviation_s, deviation_dp


def options(run, win,type):
    disc = 'D:/'
    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
    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

    x = [0.5, 1.5, 2.5]
    deviation_ms, deviation_s, deviation_dp = find_dev(x, sampling)
    time_range = 200 * delta_t
    if run == 'run_int':
        start = 0
        end = 2500
        step = 10
        eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
        eod_fe,beat_corr,beats, integral_ds,  integral_lm,  integral_gm, integral_ind,  integral_conv, max_values_gm,max_positions_gm,max_values_lm,max_positions_lm,max_values_conv,max_positions_conv,max_values_ind,max_positions_ind,max_values_ds,max_positions_ds = conspicousy(disc,deviation_s, sigma, type,eod_fe,beats, deviation_ms, 1, sampling, beat_corr, a_fe, delta_t, eod_fr, a_fr, size, phase_zero, deviation_dp)
        d = 1
        plot_cons_interactiv2(deviation_dp, deviation_s, a_fe, a_fr, sigma, delta_t,type,disc,integral_ind,0, 0, d, size, a_fe, phase_zero, beats, sampling,deviation_dp, deviation_ms, eod_fe, eod_fr, beat_corr, max_values_ds, max_positions_ds, integral_ds, max_values_lm, max_positions_lm, integral_lm, max_values_gm, max_positions_gm, integral_gm, max_values_ind, max_positions_ind, integral_ind, max_values_conv, max_positions_conv, integral_conv, share = True)
    if run == 'show_int':
        with open(disc +'Masterarbeit/labrotation/data/consp4_1.5msbeats(-495.00_1995.00_10.00).pkl', 'rb') as f:  # Python 3: open(..., 'rb')
            eod_fe,beat_corr,beats, integral_ds, integral_lm, integral_gm, integral_ind, integral_conv, max_values_gm, max_positions_gm, max_values_lm, max_positions_lm, max_values_conv, max_positions_conv, max_values_ind, max_positions_ind, max_values_ds, max_positions_ds   = pickle.load(f)
        d = 1
        plot_cons_interactiv2(deviation_dp, deviation_s, a_fe, a_fr, sigma, delta_t,type,disc,integral_ind,0, 0, d, size, a_fe, phase_zero, beats, sampling,deviation_dp, deviation_ms, eod_fe, eod_fr, beat_corr, max_values_ds, max_positions_ds, integral_ds, max_values_lm, max_positions_lm, integral_lm, max_values_gm, max_positions_gm, integral_gm, max_values_ind, max_positions_ind, integral_ind, max_values_conv, max_positions_conv, integral_conv, share = False)
    if run == 'run_dist':
        start = 5
        end = 2500
        step = 10
        eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
        period_shift,period_distance_c,period_distance_b,euc_all_in,euc_all_gm,euc_all_lm,euc_all,euc_all_ds,range_beat_all_ds, range_chirp_all_ds, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm,beat_corr_org, distance_metz_all_ds, distance_wo_mean_all_ds,distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm,distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm,distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm,distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling,deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in,distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr,distance_metz_all, distance_wo_mean_all, distance_wo_max_all,\
        distance_wo_mean_max_all = distances_func(disc, bef_c, aft_c, win, deviation_s, sigma, sampling, deviation_ms, beat_corr, size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation_dp, show_figure = True, plot_dist = False, save = True)
        #plot_dist_interactiv2(period_shift,0,deviation_ms,period_distance_c,euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds,beat_corr_org,distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s,size, amplitude, phase_zero, beats, sampling, deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, share = False,show = False)
        plot_dist_interactiv2(beat_corr, deviation_dp, delta_t, deviation_s,  a_fe, a_fr, sigma,disc,win,period_shift,1,deviation_ms,period_distance_c,euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds,beat_corr_org,distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s,size, amplitude, phase_zero, beats, sampling, deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, share = False,show = True)
    if run == 'show_dist':

        with open('../data/sim_dist3_1.5msbeats(-500.00_1995.00_5.00)'+win+'.pkl',
                  'rb') as f:  # Python 3: open(..., 'rb')
            period_distance, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds, range_beat_all_ds, range_chirp_all_ds, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm, beat_corr_org, distance_metz_all_ds, distance_wo_mean_all_ds, \
            distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, \
            distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, \
            distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, \
            distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling, \
            deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, \
            distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, \
            distance_metz_all, distance_wo_mean_all, distance_wo_max_all, \
            distance_wo_mean_max_all = pickle.load(f)
        plot_dist_interactiv2(period_distance, euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds, beat_corr_org,
                              distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds,
                              distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm,
                              distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm,
                              distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s, d,
                              size, amplitude, phase_zero, beats, sampling, deviation_ms, distance_wo_mean_max_all_in,
                              distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr,
                              distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all,
                              share=True)
    if run == 'run_all_dist':
        start = 5
        end = 2500
        step = 10
        eod_fe, beat_corr, beats = find_beats(start, end, step, eod_fr)
        plot_dist_interactiv4(bef_c, aft_c,disc,1, sampling, size, phase_zero, sigma, beat_corr, delta_t, a_fr, a_fe, deviation_dp,
                                  deviation_s, deviation_ms, eod_fe, eod_fr, share=True, show=True)



def plot_beat_corr(eod_fr):
    eod_fe = np.arange(1, 1500, 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)
    plt.subplot(gs0[0])
    np.max(beats) / eod_fr
    plt.plot(beats, beats, color=color_b)

    # 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])
    plt.plot(beats, beat_corr, color=color_b)
    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]')

    # 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('Beats [Hz]')
    plt.ylabel('EOD adj. beat [Hz]', fontsize = 10)
    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)
    plt.tight_layout()
    plt.subplots_adjust(left = 0.2,top = 0.97)
    plt.savefig(fname='../results/thesis/' + 'beat_corr_illustration')
    plt.show()
    embed()

def plot_chirp(delta_t, sampling, deviation_ms, sigma, eod_fr, eod_fe,minus_bef = -3,plus_bef = +3):
    plt.figure(figsize=[3, 2])
    left_c = minus_bef * delta_t * sampling
    right_c = plus_bef * delta_t * sampling
    time, time_cut, _ = find_times(left_c, right_c, sampling, deviation_ms[0] / (1000 * sampling))
    chirp = 60 * np.exp((-(time ** 2) / (2 * sigma ** 2)))
    frequency = eod_fr - eod_fe[0] + chirp
    plt.plot(time * 1000, frequency, color='black')
    plt.ylabel('Frequency [Hz]')
    plt.xlabel('Time [ms]')
    plt.subplots_adjust(left=0.3, bottom=0.3)
    # plt.tight_layout()
    plt.savefig(fname='../results/thesis/chirpfigure')
    plt.show()

def plot_power(beats,  results, 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,bef_c = 1.1,aft_c =-0.1):
    plt.rcParams['figure.figsize'] = (12, 7)
    plt.rcParams["legend.frameon"] = False
    colors = ['black','blue','purple','magenta','pink','orange', 'brown', 'red', 'green','lime']
    fig, ax = plt.subplots(nrows=len(results), ncols=3, sharex=True)
    all_max = [[]] * len(results)
    #embed()
    for i in range(len(results)):
        ax[i, 0].set_ylabel(results['type'][i], rotation=0, labelpad=40, color=colors[i])
        ax[i, 0].plot(beats / eod_fr + 1, np.array(results['result_frequency'][i]) / eod_fr + 1, color=colors[i])
        # plt.title(results['type'][i])
        ax[i, 1].plot(beats / eod_fr + 1, np.array(results['result_amplitude_max'][i]), color=colors[i])
        ax[i, 2].plot(beats / eod_fr + 1, np.array(results['result_amplitude_max'][i]), color=colors[i])
        ax[i,2].set_ylim([1,1.6])
        all_max[i] = np.max(np.array(results['result_amplitude_max'][i]))
    for i in range(len(results)):
        ax[i, 2].set_ylim([0, np.max(all_max)])
    plt.subplots_adjust(left=0.25)
    ii, jj = np.shape(ax)
    ax[0, 0].set_title('Most popular frequency')
    ax[0, 1].set_title('Modulation depth')
    ax[0, 2].set_title('Modulation depth (same scale)')
    for i in range(ii):

        for j in range(jj):
            ax[-1, j].set_xlabel('EOD multiples')
            ax[i, j].spines['right'].set_visible(False)
            ax[i, j].spines['top'].set_visible(False)

    def onclick(event):
        #embed()
        eod_fe = [(event.xdata-1)*eod_fr]
        nfft = 4096
        e = 0
        #sampling = 121000
        left_b = [bef_c * sampling] * len(beat_corr)
        right_b = [aft_c * sampling] * len(beat_corr)
        p = 0
        s = 0
        d = 0
        time_fish,sampled,cut, cubed, 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, chirp=False)
        #cubed, 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 = onclick_func(e,nfft,beats,  results,        disc, win, deviation_s, sigma,  d_ms, beat_corr, size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation, show_figure = show_figure, plot_dist = plot_dist, save = save,bef_c = bef_c,aft_c =aft_c, sampling = sampling)

        nfft = 4096

        results = [[]] * 1
        name = ['cubed']
        var = [cubed]
        var = [maxima_interp_b]
        samp = [sampling]
        nfft = int(4096 * samp[0] / 10000 * 2)
        i = 0
        pp, f = ml.psd(var[i] - np.mean(var[i]), Fs=samp[i], NFFT=nfft,
                       noverlap=nfft / 2)
        plt.figure()
        plt.subplot(1,2,1)
        plt.plot(f, pp)
        plt.xlim([0,2000])
        #plt.subplot(1,3,2)
        #plt.plot(time_b, cubed)
        plt.subplot(1,2,2)
        plt.plot(time_b, maxima_interp_b)
        plt.show()
    cid = fig.canvas.mpl_connect('button_press_event', onclick)
    plt.savefig('../highbeats_pdf/simulate_power')
    plt.show()

def onclick_func(e,nfft, beats,  results,        disc, win, deviation_s, sigma,  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,bef_c = 1.1,aft_c =-0.1,sampling = 100000):
    left_b = [bef_c * sampling] * len(beat_corr)
    right_b = [aft_c * sampling] * len(beat_corr)
    p = 0
    s = 0
    time_fish,sampled,cut, cubed, 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, chirp=False)

    # 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()
    # plt.plot(cubed)
    # embed()


    name = ['conv', 'df', 'df ds', 'corr', 'corr ds', 'global max', 'local max', 'cubed']
    var = [conv_b, am_df_b, am_df_ds_b, am_corr_b, am_corr_ds_b, maxima_interp_b, peaks_interp_b, cubed]
    samp = [sampling, sampling, eod_fr, sampling, eod_fr, sampling, sampling, sampling]
    # pp, f = ml.psd(eod_overlayed_chirp - np.mean(eod_overlayed_chirp), Fs=sampling, NFFT=nfft, noverlap=nfft / 2)

    return pp, f, cubed, cubed, 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

def test_different_resolutions():
    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
    eod_fr = 537# eod fish reciever


    a_fr = 1  # amplitude fish reciever
    amplitude = a_fe = 0.2  # amplitude fish emitter
    factor = 200
    sampling = eod_fr * factor
    sampling = 100000
    sampling = 121730

    sampling_fish = 500
    #start = 0
    #end = 2500
    #step = 10
    start = 510
    end = 3500
    end = eod_fr*5
    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 = 3500
    step = 250
    eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
    bef_c = 1.1
    aft_c = 1
    bef_c = 0.02
    aft_c = 0.01
    bef_c = -2
    aft_c = -1



    samp = [10000,10000,40000,11000,11234,100000, 110000,1100000]
    fig = plt.figure()
    ax = {}
    nr = 4
    ax['second'+str(-1)] = fig.add_subplot(len(samp), nr, 0 * nr + 2)
    ax['third'+str(-1)] = fig.add_subplot(len(samp), nr, 0 * nr + 3)
    ax['forth' + str(-1)] = fig.add_subplot(len(samp), nr, 0 * nr + 4)
    pp_max = []
    for i in range(len(samp)):
        eod_fe = [eod_fr*0.75]
        eod_fe = [0]
        add = 1
        nfft = 4096
        e = 0
        # sampling = 121000
        #embed()
        left_b = [bef_c * samp[i]] * len(beat_corr)
        right_b = [aft_c * samp[i]] * len(beat_corr)
        p = 0
        s = 0
        d = 0

        time_fish,sampled, cut, cubed, 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,
            samp[i], deviation_s,
            d, eod_fr, a_fr, eod_fe,
            phase_zero, p, size, s,
            sigma, a_fe, deviation_dp,
            beat_corr, chirp=False)
        # cubed, 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 = onclick_func(e,nfft,beats,  results,        disc, win, deviation_s, sigma,  d_ms, beat_corr, size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation, show_figure = show_figure, plot_dist = plot_dist, save = save,bef_c = bef_c,aft_c =aft_c, sampling = sampling)
        if single == True:
            eod_fe = [0]
            add = 1
            time, time_cut, cut = find_times(left_c[e], right_c[e], 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)
            eod_overlayed_chirp = eod_fish_r+1
        nfft = int(4096*samp[i]/samp[0]*2*2)
        results = [[]] * 1
        name = ['cubed']
        var = [cubed]
        var = cubed
        time1 = time_fish
        time2 = time_b
        time1 = time_b
        #samp = [sampling]
        #i = 0
        pp, f = ml.psd(var - np.mean(var), Fs=samp[i], NFFT=nfft,
                       noverlap=nfft / 2)
        #plt.figure()

        #plt.subplot(len(samp),3,i*3+1,sharex=True, sharey=True)
        #ax[i] = fig.add_subplot(len(samp),3,i*3+1, sharex=ax[0])
        fig.add_subplot(len(samp), 4, i * 4 + 1)
        plt.plot(f, pp, label = samp[i])
        #if i == 0:
        #    pp_max = (np.max(pp))
        plt.title(str(len(pp)))
        #plt.ylim([0,pp_max])
        plt.xlim([0, 2000])
        # plt.subplot(len(samp), 3, i*3+2, color = 'red',sharex=True, sharey=True)

        ax['second'+str(i)] = fig.add_subplot(len(samp), nr, i * nr + 2,sharex=ax['second'+str(i-1)],sharey=ax['second'+str(i-1)])
        #embed()
        plt.plot(time1, var[cut:-cut], color='red')#[cut:-cut]
        # plt.subplot(len(samp), 3, i * 3 + 3,sharex=True, sharey=True)
        #ax['third'+str(i)] = fig.add_subplot(1, 1, 1)
        plt.xlim([time_b[0], -1.9])
        ax['third'+str(i)] = fig.add_subplot(len(samp), nr, i * nr + 3,sharey=ax['third'+str(i-1)], sharex=ax['third'+str(i-1)])#sharey=ax['third'+str(i)], sharex=ax['third'+str(i)]
        plt.plot(time2, maxima_interp_b,color = 'red')
        pp, f = ml.psd(eod_overlayed_chirp  - np.mean(eod_overlayed_chirp ), Fs=samp[i], NFFT=nfft,
                       noverlap=nfft / 2)
        plt.plot(time2, eod_overlayed_chirp )
        plt.xlim([time_b[0],-1.9])
        ax['forth' + str(i)] = fig.add_subplot(len(samp), nr, i * nr + 4, sharex=ax['forth' + str(i - 1)],
                                                sharey=ax['forth' + str(i - 1)])

        plt.plot(f, pp)
        plt.xlim([0, 2000])
        # plt.subplot(1,3,2)
        # plt.plot(time_b, cubed)
        #plt.subplot(1, 2, 2)
        #plt.plot(time_b, maxima_interp_b)
    plt.legend()
    plt.show()

def singlemodulation(beat_corr, deviation_s, a_fe, a_fr):
    bef_c = 0
    aft_c = 1
    sampling = 100000
    i = 0
    fig = plt.figure()
    ax = {}
    pp_max = []
    add = 1
    nfft = 4096
    e = 0
    left_c = [bef_c * sampling] * len(beat_corr)
    right_c = [aft_c * sampling] * len(beat_corr)
    p = 0
    s = 0
    d = 0
    eod_fe = [0]
    eod_fr = 200
    time, time_cut, cut = find_times(left_c[e], right_c[e], 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)
    eod_overlayed_chirp = eod_fish_r+1
    nfft = int(4096*2*2)
    ax[1] = fig.add_subplot(1,2,1)#sharey=ax['third'+str(i)], sharex=ax['third'+str(i)]
    plt.plot(time, eod_overlayed_chirp )
    plt.xlim([0,0.05])
    ax[2] = fig.add_subplot(1,2,2)
    pp, f = ml.psd(eod_overlayed_chirp, Fs=sampling, NFFT=nfft,
                   noverlap=nfft / 2)
    plt.plot(f, pp)
    plt.xlim([-50, 2000])
    plt.legend()
    plt.savefig('../highbeats_pdf/psd_singlemodulation.pdf')
    plt.show()


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
    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
    eod_fr = 537# eod fish reciever
    factor = 200
    sampling_fish = 500
    size = [120]

    d = 1
    x = [ 1.5, 2.5,0.5,]
    x = [ 1.5]
    sampling = 100000
    time_range = 200 * delta_t
    deviation_ms, deviation_s, deviation_dp = find_dev(x, sampling)


    #sampling = 750
    #sampling = 522300

    bef_c = -2
    aft_c = -1
    deviation_ms, deviation_s, deviation_dp = find_dev(x, sampling)
    #minus_bef = bef_c/(delta*sampling)#minus_bef = -30
    #plus_bef = aft_c/(delta*sampling)
    minus_bef = bef_c
    plus_bef = aft_c
    #plus_bef = -10
    eod_fe = np.array([510, 560,850,1010])
    eod_fr = 750
    eod_fe = np.array([eod_fr + 10, eod_fr + 100,eod_fr + eod_fr/2,eod_fr + eod_fr/2+100,eod_fr*3 + 10,eod_fr*4 + 10])
    a_fe = 0.2
    a_fr = 1
    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]
    plot_eod(beat_corr, minus_bef, plus_bef,3,2,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.2992, 5],shift_phase = (2 * np.pi) / 4)

    embed()

    start = 5
    end = 3500
    step = 250
    eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
    singlemodulation(beat_corr, deviation_s, a_fe, a_fr)


    embed()
    delta_t = 1
    results = power_func(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)
    plot_power(beats, results,'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,bef_c = bef_c,aft_c =aft_c)

    embed()


    start = 0
    end = 2500
    step = 25
    eod_fe, beat_corr, beats = find_beats(start,end,step,eod_fr)
    distance_metz_all_df_ds,distance_metz_all_corr_ds,period_shift,period_distance_c,period_distance_b,euc_all_in,euc_all_gm,euc_all_lm,euc_all,euc_all_ds,range_beat_all_ds, range_chirp_all_ds, range_chirp_all, range_beat_all, range_chirp_all_in, range_beat_all_in, range_beat_all_gm, range_chirp_all_gm, range_beat_all_lm, range_chirp_all_lm,beat_corr_org, distance_metz_all_ds, distance_wo_mean_all_ds,distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm,distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm,distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm,distance_wo_max_all_gm, p, s, d, size, phase_zero, beats, sampling,deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in,distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr,distance_metz_all, distance_wo_mean_all, distance_wo_max_all,\
    distance_wo_mean_max_all = distances_func(disc, bef_c, aft_c, 'w2', deviation_s, sigma, sampling, deviation_ms, beat_corr, size, phase_zero, delta_t, a_fr, a_fe, eod_fr, eod_fe, deviation_dp, show_figure = True, plot_dist = False, save = True)
    plot_dist_interactiv2(distance_metz_all_df_ds,distance_metz_all_corr_ds,beat_corr, deviation_dp, delta_t, deviation_s, a_fe, a_fr, sigma,disc,win,period_shift,0,deviation_ms,period_distance_c,euc_all_in, euc_all_gm, euc_all_lm, euc_all, euc_all_ds,beat_corr_org,distance_metz_all_ds, distance_wo_mean_all_ds, distance_wo_mean_max_all_ds, distance_wo_max_all_ds, distance_metz_all_lm, distance_wo_mean_all_lm, distance_wo_mean_max_all_lm, distance_wo_max_all_lm, distance_metz_all_gm, distance_wo_mean_all_gm, distance_wo_mean_max_all_gm, distance_wo_max_all_gm, p, s,size, amplitude, phase_zero, beats, sampling, deviation_ms, distance_wo_mean_max_all_in, distance_wo_max_all_in, distance_metz_all_in, distance_wo_mean_all_in, eod_fe, eod_fr, distance_metz_all, distance_wo_mean_all, distance_wo_max_all, distance_wo_mean_max_all, share = False, show = True)
    embed()














    minus_bef = -30
    plus_bef = -10
    eod_fe = np.array([510, 560,850,1010])
    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]
    plot_eod(beat_corr, minus_bef, plus_bef,2,2,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.2992, 5],shift_phase = (2 * np.pi) / 4)


    minus_bef = -4
    plus_bef = +4
    eod_fe = np.array([510, 580,750,1010,1250,1510])
    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]
    plot_eod(beat_corr, minus_bef, plus_bef,3,2,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.2992, 7.5],shift_phase = (2 * np.pi) / 4)



    minus_bef = -4
    plus_bef = +4
    eod_fe = np.array([510])
    eod_fe = np.array([510,750])
    eod_fe = np.array([580])#[510]
    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]
    pp = [1,3,6,9]#[2,7]
    #embed()
    plot_eod_phase(beat_corr, minus_bef, plus_bef,2,2,a_fe, delta_t, a_fr, size, beat_corr, 0, pp, a_fe,
                   phase_zero, deviation_ms, eod_fe, deviation_dp, eod_fr, sampling, factor,fs = [6.2992, 5],shift_phase = (2 * np.pi) / 4)

    embed()



