import matplotlib
#matplotlib.use('Qt5Agg')
import numpy as np
import math
import scipy
import scipy.signal as sig
import os
import pickle as pkl
from scipy import constants as cst
from scipy import io
from scipy import signal
from PyQt5 import QtGui, QtWidgets
from numpy import linalg as lin
import pyqtgraph as pg
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import pandas as pd
import h5py
from timeit import default_timer as timer
import platform
from scipy.optimize import curve_fit
#from sys import platform
import pyabf
from matplotlib.ticker import EngFormatter
from pprint import pprint

# if not 'verbose' in globals():
#     verbose = False
#
# verboseprint = print if verbose else lambda *a, **k: None
#

def GetKClConductivity(Conc, Temp):
    p = pkl.load(open('KCl_ConductivityValues.p', 'rb'))
    if Conc==1.0:
        Conc = np.uint(1)
    return np.polyval(p[str(Conc)], Temp)

def GetTempFromKClConductivity(Conc, Cond):
    p = pkl.load(open('KCl_ConductivityValues.p', 'rb'))
    if Conc==1.0:
        Conc = np.uint(1)
    return (Cond-p[str(Conc)][1])/p[str(Conc)][0]

def RecursiveLowPassFast(signal, coeff, samplerate):
    padlen = np.uint64(samplerate)
    prepadded = np.ones(padlen)*np.mean(signal[0:1000])
    signaltofilter = np.concatenate((prepadded, signal))

    mltemp = scipy.signal.lfilter([1 - coeff['a'], 0], [1, -coeff['a']], signaltofilter)
    vltemp = scipy.signal.lfilter([1 - coeff['a'], 0], [1, -coeff['a']], np.square(signaltofilter - mltemp))

    ml = np.delete(mltemp, np.arange(padlen))
    vl = np.delete(vltemp, np.arange(padlen))

    sl = ml - coeff['S'] * np.sqrt(vl)
    Ni = len(signal)
    points = np.array(np.where(signal <= sl)[0])
    to_pop = np.array([])
    for i in range(1, len(points)):
        if points[i] - points[i - 1] == 1:
            to_pop = np.append(to_pop, i)
    points = np.unique(np.delete(points, to_pop))
    RoughEventLocations = []
    NumberOfEvents = 0

    for i in points:
        if NumberOfEvents is not 0:
            if i >= RoughEventLocations[NumberOfEvents-1][0] and i <= RoughEventLocations[NumberOfEvents-1][1]:
                continue
        NumberOfEvents += 1
        start = i
        El = ml[i] + coeff['E'] * np.sqrt(vl[i])
        Mm = ml[i]
        Vv = vl[i]
        duration = 0
        endp = start
        if (endp + 1) < len(signal):
            while signal[endp + 1] < El and endp < (Ni - 2):# and duration < coeff['eventlengthLimit']*samplerate:
                duration += 1
                endp += 1
        if duration >= coeff['eventlengthLimit'] * samplerate or endp > (Ni - 10): # or duration <= coeff['minEventLength'] * samplerate:
            NumberOfEvents -= 1
            continue
        else:
            k = start
            while signal[k] < Mm and k > 1:
                k -= 1
            start = k - 1
            k2 = i + 1
            #while signal[k2] > Mm:
            #    k2 -= 1
            #endp = k2
            if start < 0:
                start = 0
            RoughEventLocations.append((start, endp, ml[start], vl[start]))

    return np.array(RoughEventLocations)#, ml, vl, sl

def RecursiveLowPassFastUp(signal, coeff, samplerate):
    ml = scipy.signal.lfilter([1 - coeff['a'], 0], [1, -coeff['a']], signal)
    vl = scipy.signal.lfilter([1 - coeff['a'], 0], [1, -coeff['a']], np.square(signal - ml))
    sl = ml + coeff['S'] * np.sqrt(vl)
    Ni = len(signal)
    points = np.array(np.where(signal>=sl)[0])
    to_pop=np.array([])
    for i in range(1,len(points)):
        if points[i] - points[i - 1] == 1:
            to_pop=np.append(to_pop, i)
    points = np.delete(points, to_pop)

    points =np.delete(points, np.array(np.where(points == 0)[0]))

    RoughEventLocations = []
    NumberOfEvents=0
    for i in points:
        if NumberOfEvents is not 0:
            if i >= RoughEventLocations[NumberOfEvents-1][0] and i <= RoughEventLocations[NumberOfEvents-1][1]:
                continue
        NumberOfEvents += 1
        start = i
        El = ml[i] + coeff['E'] * np.sqrt(vl[i])
        Mm = ml[i]
        duration = 0
        while signal[i + 1] > El and i < (Ni - 2) and duration < coeff['eventlengthLimit']*samplerate:
            duration += 1
            i += 1
        if duration >= coeff['eventlengthLimit']*samplerate or i > (Ni - 10):
            NumberOfEvents -= 1
        else:
            k = start
            while signal[k] > Mm and k > 2:
                k -= 1
            start = k - 1
            k2 = i + 1
            while signal[k2] > Mm:
                k2 -= 1
            endp = k2
            RoughEventLocations.append((start, endp, ml[start], vl[start]))

    return np.array(RoughEventLocations)

def ImportABF(datafilename):
    # abf = pyabf.ABF(datafilename)
    # # abf.info()  # shows what is available
    # output={'type': 'Clampfit', 'graphene': 0, 'samplerate': abf.pointsPerSec, 'i1': -20000./65536 * abf.dataY, 'v1': abf.dataC, 'filename': datafilename}
    # return output
    abf = pyabf.ABF(datafilename)
    #abf.info()  # shows what is available
    #output={'type': 'Clampfit', 'graphene': 0, 'samplerate': abf.pointsPerSec, 'i1': -20000./65536 * abf.dataY, 'v1': abf.dataC, 'filename': datafilename}
    output = {'type': 'Clampfit', 'graphene': 0, 'samplerate': abf.dataRate,
              'i1': abf.data[0] * 1e-12,
              'v1': abf.sweepC,
              'filename': datafilename}
    return output

def ImportAxopatchData(datafilename):
    x=np.fromfile(datafilename, np.dtype('>f4'))
    f=open(datafilename, 'rb')
    graphene=0
    for i in range(0, 10):
        a=str(f.readline())
        print(a)
        if 'Acquisition' in a or 'Sample Rate' in a:
            samplerate=int(''.join(i for i in a if i.isdigit()))/1000
        if 'FEMTO preamp Bandwidth' in a:
            femtoLP=int(''.join(i for i in a if i.isdigit()))
        if 'I_Graphene' in a:
            graphene=1
            print('This File Has a Graphene Channel!')
    end = len(x)
    if graphene:
        #pore current
        i1 = x[250:end-3:4]
        #graphene current
        i2 = x[251:end-2:4]
        #pore voltage
        v1 = x[252:end-1:4]
        #graphene voltage
        v2 = x[253:end:4]
        print('The femto was set to : {} Hz, if this value was correctly entered in the LabView!'.format(str(femtoLP)))
        output={'FemtoLowPass': femtoLP, 'type': 'Axopatch', 'graphene': 1, 'samplerate': samplerate, 'i1': i1, 'v1': v1, 'i2': i2, 'v2': v2, 'filename': datafilename}
    else:
        i1 = np.array(x[250:end-1:2])
        v1 = np.array(x[251:end:2])
        output={'type': 'Axopatch', 'graphene': 0, 'samplerate': samplerate, 'i1': i1, 'v1': v1, 'filename': datafilename}
    return output

def ImportChimeraRaw(datafilename):
    matfile=io.loadmat(str(os.path.splitext(datafilename)[0]))
    #buffersize=matfile['DisplayBuffer']
    data = np.fromfile(datafilename, np.dtype('<u2'))
    samplerate = np.float64(matfile['ADCSAMPLERATE'])
    TIAgain = np.int32(matfile['SETUP_TIAgain'])
    preADCgain = np.float64(matfile['SETUP_preADCgain'])
    currentoffset = np.float64(matfile['SETUP_pAoffset'])
    ADCvref = np.float64(matfile['SETUP_ADCVREF'])
    ADCbits = np.int32(matfile['SETUP_ADCBITS'])
    if 'blockLength' in matfile:
        blockLength=np.int32(matfile['blockLength'])
    else:
        blockLength=1048576

    closedloop_gain = TIAgain * preADCgain
    bitmask = (2 ** 16 - 1) - (2 ** (16 - ADCbits) - 1)
    data = -ADCvref + (2 * ADCvref) * (data & bitmask) / 2 ** 16
    data = (data / closedloop_gain + currentoffset)
    data.shape = [data.shape[1], ]
    output = {'matfilename': str(os.path.splitext(datafilename)[0]),'i1raw': data, 'v1': np.float64(matfile['SETUP_biasvoltage']), 'samplerateRaw': np.int64(samplerate), 'type': 'ChimeraRaw', 'filename': datafilename, 'graphene': 0, 'blockLength':blockLength}
    return output

def Reshape1DTo2D(inputarray, buffersize):
    npieces = int(len(inputarray)/buffersize)
    voltages = np.array([], dtype=np.float64)
    currents = np.array([], dtype=np.float64)

    for i in range(1, npieces+1):
        if i % 2 == 1:
            currents = np.append(currents, inputarray[(i-1)*buffersize:i*buffersize-1], axis=0)
            #print('Length Currents: {}'.format(len(currents)))
        else:
            voltages = np.append(voltages, inputarray[(i-1)*buffersize:i*buffersize-1], axis=0)
            #print('Length Voltages: {}'.format(len(voltages)))

    v1 = np.ones((len(voltages)), dtype=np.float64)
    i1 = np.ones((len(currents)), dtype=np.float64)
    v1[:]=voltages
    i1[:]=currents

    out = {'v1': v1, 'i1': i1}
    print('Currents:' + str(v1.shape))
    print('Voltages:' + str(i1.shape))
    return out

def ImportChimeraData(datafilename):
    matfile = io.loadmat(str(os.path.splitext(datafilename)[0]))
    samplerate = matfile['ADCSAMPLERATE']
    if samplerate<4e6:
        data = np.fromfile(datafilename, np.dtype('float64'))
        buffersize = int(matfile['DisplayBuffer'])
        out = Reshape1DTo2D(data, buffersize)
        output = {'i1': out['i1'], 'v1': out['v1'], 'samplerate':float(samplerate), 'type': 'ChimeraNotRaw', 'filename': datafilename,'graphene': 0}
    else:
        output = ImportChimeraRaw(datafilename)
    return output

def OpenFile(filename = '', ChimeraLowPass = 10e3, approxImpulseResponse=False, Split=False):
    if ChimeraLowPass==None:
        ChimeraLowPass=10e3
    if filename == '':
        datafilename = QtGui.QFileDialog.getOpenFileName()
        datafilename=datafilename[0]
        print(datafilename)
    else:
        datafilename=filename
    print('Loading file... ' +filename)
    if datafilename[-3::] == 'dat':
        isdat = 1
        output = ImportAxopatchData(datafilename)
    elif datafilename[-3::] == 'log':
        isdat = 0
        output = ImportChimeraData(datafilename)
        if output['type'] is 'ChimeraRaw':  # Lowpass and downsample
            print('length: ' + str(len(output['i1raw'])))
            Wn = round(2 * ChimeraLowPass / output['samplerateRaw'], 4)  # [0,1] nyquist frequency
            b, a = signal.bessel(4, Wn, btype='low', analog=False)  # 4-th order digital filter
            if approxImpulseResponse:
                z, p, k = signal.tf2zpk(b, a)
                eps = 1e-9
                r = np.max(np.abs(p))
                approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r)))
                Filt_sig=(signal.filtfilt(b, a, output['i1raw'], method='gust', irlen=approx_impulse_len))
            else:
                Filt_sig=(signal.filtfilt(b, a, output['i1raw'], method='gust'))
            ds_factor = np.int16(np.ceil(output['samplerateRaw'] / (4 * ChimeraLowPass)))

            howmanytimes = np.floor(ds_factor/10)
            output['i1'] = Filt_sig

            for i in np.arange(howmanytimes):
                output['i1'] = scipy.signal.decimate(Filt_sig, 10)
            output['i1'] = scipy.signal.decimate(Filt_sig, int(np.mod(ds_factor, 10)))

            output['samplerate'] = output['samplerateRaw'] / ds_factor
            output['v1'] = output['v1']*np.ones(len(output['i1']))
            print('Samplerate after filtering:' + str(output['samplerate']))
            print('new length: ' + str(len(output['i1'])))

            if Split:
                splitNr = len(output['i1raw'])/output['blockLength']
                assert(splitNr.is_integer())
                rawSplit = np.array_split(output['i1raw'], splitNr)
                Filt_sigSplit = []
                for raw in rawSplit:
                    if approxImpulseResponse:
                        z, p, k = signal.tf2zpk(b, a)
                        eps = 1e-9
                        r = np.max(np.abs(p))
                        approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r)))
                        signalSplit=signal.filtfilt(b, a, raw, method='gust', irlen=approx_impulse_len)
                    else:
                        signalSplit=signal.filtfilt(b, a, raw, method = 'gust')
                    Filt_sigSplit.append(scipy.signal.resample(signalSplit, int(len(raw) / ds_factor)))
                output['i1Cut'] = Filt_sigSplit

        if max(abs(output['i1'])) > 1e-6:
            print('converting to SI units')
            output['i1']=1e-9*output['i1']
            output['v1']=1e-3*output['v1']

    elif datafilename[-3::] == 'abf':
        output = ImportABF(datafilename)
    st = os.stat(datafilename)
    if platform.system() == 'Darwin':
        print('Platform is ' + platform.system())
        output['TimeFileWritten'] = st.st_birthtime
        output['TimeFileLastModified'] = st.st_mtime
        output['ExperimentDuration'] = st.st_mtime - st.st_birthtime
    elif platform.system() == 'Windows':
        print('Platform is WinShit')
        output['TimeFileWritten'] = st.st_ctime
        output['TimeFileLastModified'] = st.st_mtime
        output['ExperimentDuration'] = st.st_mtime - st.st_ctime
    else:
        print('Platform is ' + platform.system() +
            ', might not get accurate results.')
        try:
            output['TimeFileWritten'] = st.st_ctime
            output['TimeFileLastModified'] = st.st_mtime
            output['ExperimentDuration'] = st.st_mtime - st.st_ctime
        except:
            raise Exception('Platform not detected')

    return output

def MakeIVData(output, approach = 'mean', delay = 0.7, UseVoltageRange=0):
    (ChangePoints, sweepedChannel) = CutDataIntoVoltageSegments(output)
    if ChangePoints is 0:
        return 0

    if output['graphene']:
        currents = ['i1', 'i2']
    else:
        currents = ['i1']

    Values = output[sweepedChannel][ChangePoints]
    Values = np.append(Values, output[sweepedChannel][::-1][0])
    delayinpoints = np.int64(delay * output['samplerate'])
    if delayinpoints>ChangePoints[0]:
        raise ValueError("Delay is longer than Changepoint")

    #   Store All Data
    All={}
    for current in currents:
        Item = {}
        ItemExp = {}
        l = len(Values)
        Item['Voltage'] = np.zeros(l)
        Item['StartPoint'] = np.zeros(l,dtype=np.uint64)
        Item['EndPoint'] = np.zeros(l, dtype=np.uint64)
        Item['Mean'] = np.zeros(l)
        Item['STD'] = np.zeros(l)
        Item['SweepedChannel'] = sweepedChannel
        ItemExp['SweepedChannel'] = sweepedChannel
        Item['YorkFitValues'] = {}
        ItemExp['YorkFitValues'] = {}
        ### MEAN METHOD###
        # First item
        Item['Voltage'][0] = Values[0]
        trace = output[current][0 + delayinpoints:ChangePoints[0]]
        Item['StartPoint'][0] = 0 + delayinpoints
        Item['EndPoint'][0] = ChangePoints[0]
        Item['Mean'][0] = np.mean(trace)
        isProblematic = math.isclose(Values[0], 0, rel_tol=1e-3) or math.isclose(Values[0], np.min(Values), rel_tol=1e-3) or math.isclose(Values[0], np.max(Values), rel_tol=1e-3)
        if sweepedChannel is 'v2' and isProblematic:
            print('In Problematic If for Values: {}, index {}'.format(Values[0],0))
            Item['STD'][0] = np.std(trace[-np.int64(output['samplerate']):])
        else:
            Item['STD'][0] = np.std(trace)
        for i in range(1, len(Values) - 1):
            trace = output[current][ChangePoints[i - 1] + delayinpoints : ChangePoints[i]]
            Item['Voltage'][i] = Values[i]
            Item['StartPoint'][i] = ChangePoints[i - 1]+ delayinpoints
            Item['EndPoint'][i] = ChangePoints[i]
            Item['Mean'][i] = np.mean(trace)
            isProblematic = math.isclose(Values[i], 0, rel_tol=1e-3) or math.isclose(Values[i], np.min(Values),
                                                                                     rel_tol=1e-3) or math.isclose(Values[i], np.max(Values), rel_tol=1e-3)
            if sweepedChannel is 'v2' and isProblematic:
                print('In Problematic If for Values: {}, index {}'.format(Values[i], i))
                Item['STD'][i] = np.std(trace[-np.int64(output['samplerate']):])
            else:
                Item['STD'][i] = np.std(trace)
            print('{}, {},{}'.format(i, ChangePoints[i - 1] + delayinpoints, ChangePoints[i]))
        # Last
        if 1:
            trace = output[current][ChangePoints[len(ChangePoints) - 1] + delayinpoints : len(output[current]) - 1]
            Item['Voltage'][-1:] = Values[len(Values) - 1]
            Item['StartPoint'][-1:] = ChangePoints[len(ChangePoints) - 1]+delayinpoints
            Item['EndPoint'][-1:] = len(output[current]) - 1
            Item['Mean'][-1:] = np.mean(trace)
            isProblematic = math.isclose(Values[-1:], 0, rel_tol=1e-3) or math.isclose(Values[-1:], np.min(Values),
                                                                                     rel_tol=1e-3) or math.isclose(
                Values[-1:], np.max(Values), rel_tol=1e-3)
            if sweepedChannel is 'v2' and isProblematic:
                print('In Problematic If for Values: {}, index {}'.format(Values[-1:], i))
                Item['STD'][-1:] = np.std(trace[-np.int64(output['samplerate']):])
            else:
                Item['STD'][-1:] = np.std(trace)

        ## GET RECTIFICATION FACTOR
        parts = {'pos': Item['Voltage'] > 0, 'neg': Item['Voltage'] < 0}
        a = {}
        b = {}
        sigma_a = {}
        sigma_b = {}
        b_save = {}
        x_values = {}

        for part in parts:
            (a[part], b[part], sigma_a[part], sigma_b[part], b_save[part]) = YorkFit(
                Item['Voltage'][parts[part]].flatten(), Item['Mean'][parts[part]].flatten(),
                1e-12 * np.ones(len(Item['Voltage'][parts[part]].flatten())),
                Item['STD'][parts[part]].flatten())
        factor = b['neg'] / b['pos']
        Item['Rectification'] = factor
        ItemExp['Rectification'] = factor

        ###EXPONENTIAL METHOD###
        if approach is not 'mean':
            l = 1000
            ItemExp['StartPoint'] = np.zeros(l, dtype=np.uint64)
            ItemExp['EndPoint'] = np.zeros(l, dtype=np.uint64)
            baselinestd = np.std(output[current][0:np.int64(1 * output['samplerate'])])
            movingstd = pd.rolling_std(output[current], 10)

            # Extract The Desired Parts
            ind = (np.abs(movingstd - baselinestd) < 1e-9) & (np.abs(output[current]) < np.max(output[current]) / 2)

            # How Many Parts?
            restart = True
            counter = 0
            for i, value in enumerate(ind):
                if value and restart:
                    ItemExp['StartPoint'][counter] = i
                    print('Start: {}\t'.format(i))
                    restart = False
                elif not value and not restart:
                    if ((i - 1) - ItemExp['StartPoint'][counter]) > 1 * output['samplerate']:
                        ItemExp['EndPoint'][counter] = i - 1
                        print('End: {}'.format(i - 1))
                        restart = True
                        counter += 1
                    else:
                        restart = True

            if not restart:
                counter += 1
                ItemExp['EndPoint'][counter] = len(ind)

            ItemExp['StartPoint'] = ItemExp['StartPoint'][np.where(ItemExp['StartPoint'])]
            ItemExp['EndPoint'] = ItemExp['EndPoint'][np.where(ItemExp['EndPoint'])]
            NumberOfPieces = len(ItemExp['StartPoint'])
            ItemExp['Voltage'] = np.zeros(NumberOfPieces)
            ItemExp['ExpValues'] = np.zeros((3, NumberOfPieces))
            ItemExp['STD'] = np.zeros(NumberOfPieces)
            ItemExp['Mean'] = np.zeros(NumberOfPieces)

            # Make It Piecewise
            for i in np.arange(0, NumberOfPieces):
                trace = output[current][ItemExp['StartPoint'][i]:ItemExp['EndPoint'][i]]
                popt, pcov = MakeExponentialFit(np.arange(len(trace)) / output['samplerate'], trace)
                ItemExp['Voltage'][i] = output[sweepedChannel][ItemExp['StartPoint'][i]]
                if popt[0]:
                    ItemExp['ExpValues'][:, i] = popt
                    ItemExp['Mean'][i] = popt[2]
                    try:
                        ItemExp['STD'][i] = np.sqrt(np.diag(pcov))[2]
                    except RuntimeWarning:
                        ItemExp['STD'][i] = np.std(trace)/np.sqrt(len(trace))
                        print('ExpFit: STD of voltage {} failed calculating...'.format(ItemExp['Voltage'][i]))
                else:
                    print('Exponential Fit on for ' + current + ' failed at V=' + str(ItemExp['Voltage'][0]))
                    ItemExp['Mean'][i] = np.mean(trace)
                    ItemExp['STD'][i] = np.std(trace)

                ## FIT THE EXP CALCULATED VALUES ###
                (a, b, sigma_a, sigma_b, b_save) = YorkFit(ItemExp['Voltage'], ItemExp['Mean'],
                                                           1e-12 * np.ones(len(ItemExp['Voltage'])), ItemExp['STD'])
                x_fit = np.linspace(min(Item['Voltage']), max(Item['Voltage']), 1000)
                y_fit = scipy.polyval([b, a], x_fit)
                ItemExp['YorkFitValues'] = {'x_fit': x_fit, 'y_fit': y_fit, 'Yintercept': a, 'Slope': b,
                                            'Sigma_Yintercept': sigma_a,
                                            'Sigma_Slope': sigma_b, 'Parameter': b_save}
                All[current + '_Exp'] = ItemExp

        ##Remove NaNs
        nans = np.isnan(Item['Mean'][:])
        nans = np.logical_not(nans)
        Item['Voltage'] = Item['Voltage'][nans]
        Item['StartPoint'] = Item['StartPoint'][nans]
        Item['EndPoint'] = Item['EndPoint'][nans]
        Item['Mean'] = Item['Mean'][nans]
        Item['STD'] = Item['STD'][nans]

        ## Restrict to set voltage Range
        if UseVoltageRange is not 0:
            print('Voltages Cut')
            ind = np.argwhere((Item['Voltage'] >= UseVoltageRange[0]) & (Item['Voltage'] <= UseVoltageRange[1]))
            Item['Voltage'] = Item['Voltage'][ind].flatten()
            Item['StartPoint'] = Item['StartPoint'][ind].flatten()
            Item['EndPoint'] = Item['EndPoint'][ind].flatten()
            Item['Mean'] = Item['Mean'][ind].flatten()
            Item['STD'] = Item['STD'][ind].flatten()

        ## FIT THE MEAN CALCULATED VALUES ###
        (a, b, sigma_a, sigma_b, b_save) = YorkFit(Item['Voltage'].flatten(), Item['Mean'].flatten(), 1e-9 * np.ones(len(Item['Voltage'])), Item['STD'].flatten())
        x_fit = np.linspace(min(Item['Voltage']), max(Item['Voltage']), 1000)
        y_fit = scipy.polyval([b, a], x_fit)
        Item['YorkFitValues'] = {'x_fit': x_fit, 'y_fit': y_fit, 'Yintercept': a, 'Slope': b, 'Sigma_Yintercept': sigma_a,
                         'Sigma_Slope': sigma_b, 'Parameter': b_save}

        All[current] = Item
        All['Currents'] = currents

    return All

def ExpFunc(x, a, b, c):
    return a * np.exp(-b * x) + c

def YorkFit(X, Y, sigma_X, sigma_Y, r=0):
    N_itermax=10 #maximum number of interations
    tol=1e-15 #relative tolerance to stop at
    N = len(X)
    temp = np.matrix([X, np.ones(N)])
    #make initial guess at b using linear squares

    tmp = np.matrix(Y)*lin.pinv(temp)
    b_lse = np.array(tmp)[0][0]
    #a_lse=tmp(2);
    b = b_lse #initial guess
    omega_X = np.true_divide(1,np.power(sigma_X,2))
    omega_Y = np.true_divide(1, np.power(sigma_Y,2))
    alpha=np.sqrt(omega_X*omega_Y)
    b_save = np.zeros(N_itermax+1) #vector to save b iterations in
    b_save[0]=b

    for i in np.arange(N_itermax):
        W=omega_X*omega_Y/(omega_X+b*b*omega_Y-2*b*r*alpha)

        X_bar=np.sum(W*X)/np.sum(W)
        Y_bar=np.sum(W*Y)/np.sum(W)

        U=X-X_bar
        V=Y-Y_bar

        beta=W*(U/omega_Y+b*V/omega_X-(b*U+V)*r/alpha)

        b=sum(W*beta*V)/sum(W*beta*U)
        b_save[i+1]=b
        if np.abs((b_save[i+1]-b_save[i])/b_save[i+1]) < tol:
            break

    a=Y_bar-b*X_bar
    x=X_bar+beta
    y=Y_bar+b*beta
    x_bar=sum(W*x)/sum(W)
    y_bar=sum(W*y)/sum(W)
    u=x-x_bar
    #%v=y-y_bar
    sigma_b=np.sqrt(1/sum(W*u*u))
    sigma_a=np.sqrt(1./sum(W)+x_bar*x_bar*sigma_b*sigma_b)
    return (a, b, sigma_a, sigma_b, b_save)

def DoubleExpFunc(x, Y0, percentFast, plateau, Kfast, Kslow):
    SpanFast = (Y0-plateau) * percentFast * 0.01
    SpanSlow = (Y0-plateau) *(100-percentFast) * 0.01
    return plateau + ((Y0-plateau) * percentFast * 0.01)*np.exp(-Kfast*x) + ((Y0-plateau) *(100-percentFast) * 0.01)*np.exp(-Kslow*x)

def MakeExponentialFit(xdata, ydata):
    try:
        popt, pcov = curve_fit(ExpFunc, xdata, ydata, method='lm', xtol = 1e-20, ftol = 1e-12, gtol = 1e-20)
        #popt, pcov = curve_fit(DoubleExpFunc, xdata, ydata, p0 = (ydata[0], 50, ydata[-1], 1 / 15, 1 / 60))
        return (popt, pcov)
    except RuntimeError:
        popt = (0,0,0)
        pcov = 0
        return (popt, pcov)

def CutDataIntoVoltageSegments(output):
    sweepedChannel = ''
    if output['type'] == 'ChimeraNotRaw' or (output['type'] == 'Axopatch' and not output['graphene']):
        ChangePoints = np.where(np.diff(output['v1']))[0]
        sweepedChannel = 'v1'
        if len(ChangePoints) is 0:
            print('No voltage sweeps in this file:\n{}'.format(output['filename']))
            return (0,0)
    elif (output['type'] == 'Axopatch' and output['graphene']):
        ChangePoints = np.where(np.diff(output['v1']))[0]
        sweepedChannel = 'v1'
        if len(ChangePoints) is 0:
            ChangePoints = np.where(np.diff(output['v2']))[0]
            if len(ChangePoints) is 0:
                print('No voltage sweeps in this file')
                return (0,0)
            else:
                sweepedChannel = 'v2'
    print('Cutting into Segments...\n{} change points detected in channel {}...'.format(len(ChangePoints), sweepedChannel))
    return (ChangePoints, sweepedChannel)

def CalculatePoreSize(G, L, s):
    #https://doi.org/10.1088/0957-4484/22/31/315101
    return (G+np.sqrt(G*(G+16*L*s/np.pi)))/(2*s)

def CalculateCapillarySize(G, D=0.3e-3, t=3.3e-3, s=10.5):
    #DOI: 10.1021/nn405029j
    #s=conductance (10.5 S/m at 1M KCl)
    #t=taper length (3.3 mm)
    #D=diameter nanocapillary at the shaft (0.3 mm)
    return G*(4*t/np.pi+0.5*D)/(s*D-0.5*G)

def CalculateShrunkenCapillarySize(G, D=0.3e-3, t=3.3e-3, s=10.5, ts=500e-9, Ds=500e-9):
    #DOI: 10.1021/nn405029j
    #s=conductance (10.5 S/m at 1M KCl)
    #t=taper length (3.3 mm)
    #D=diameter nanocapillary at the shaft (0.3 mm)
    #ts=taper length of the shrunken part (543nm fit)
    #Ds=diameter nanocapillary at the shrunken part (514nm fit)
    return G * D * (8 * ts / np.pi + Ds)/((2 * s * D * Ds) - (G * (8 * t / np.pi + Ds)))

def CalculateConductance(L, s, d):
    return s*1/(4*L/(np.pi*d**2)+1/d)

def GetSurfaceChargeFromLeeEquation(s, c, D, G, L, A, B, version=1):
    lDu = 1/(8*np.pi*B*D*s)*(
        version * np.sqrt((4*np.pi*A*D**2*s+np.pi*B*D**2*s - 4*B*G*L - 8*np.pi*D*G)**2
                 - 16*np.pi*B*D*s*(np.pi*A*D**3*s - 4*A*D*G*L - 2*np.pi*D**2*G))
        - 4*np.pi*A*D**2*s - np.pi*B*D**2*s + 4*B*G*L + 8*np.pi*D*G
    )
    e = cst.elementary_charge
    Na = cst.Avogadro
    return lDu*(2 * c * Na * 1e3 * e)
def ConductivityFromConductanceModel(L, d, G):
    return G * (4 * L / (np.pi * d ** 2) + 1 / d)


def RefinedEventDetection(out, AnalysisResults, signals, limit):
    for sig in signals:
        #Choose correct reference
        if sig is 'i1_Up':
            sig1 = 'i1'
        elif sig is 'i2_Up':
            sig1 = 'i2'
        else:
            sig1 = sig
        if sig is 'i1' or 'i1_Up':
            volt='v1'
        elif sig is 'i2' or 'i2_Up':
            volt='v2'

        if len(AnalysisResults[sig]['RoughEventLocations']) > 1:
            startpoints = np.uint64(AnalysisResults[sig]['RoughEventLocations'][:, 0])
            endpoints = np.uint64(AnalysisResults[sig]['RoughEventLocations'][:, 1])
            localBaseline = AnalysisResults[sig]['RoughEventLocations'][:, 2]
            localVariance = AnalysisResults[sig]['RoughEventLocations'][:, 3]

            CusumBaseline=500
            numberofevents = len(startpoints)
            AnalysisResults[sig]['StartPoints'] = startpoints
            AnalysisResults[sig]['EndPoints'] = endpoints
            AnalysisResults[sig]['LocalBaseline'] = localBaseline
            AnalysisResults[sig]['LocalVariance'] = localVariance
            AnalysisResults[sig]['NumberOfEvents'] = len(startpoints)

            deli = np.zeros(numberofevents)
            dwell = np.zeros(numberofevents)
            fitvalue = np.zeros(numberofevents)
            AllEvents = []

            #####FIX THE UP AND DOWN STORY!! ALL THE PROBLEMS ARE IN HERE...
            for i in range(numberofevents):
                length = endpoints[i] - startpoints[i]
                if out[sig1][startpoints[i]] < 0 and (sig == 'i1_Up' or sig == 'i2_Up'):
                    if length <= limit and length > 3:
                        # Impulsion Fit to minimal value
                        fitvalue[i] = np.max(out[sig1][startpoints[i]+np.uint(1):endpoints[i]-np.uint(1)])
                    elif length > limit:
                        fitvalue[i] = np.mean(out[sig1][startpoints[i]+np.uint(5):endpoints[i]-np.uint(5)])
                    else:
                        fitvalue[i] = np.max(out[sig1][startpoints[i]:endpoints[i]])
                    deli[i] = -(localBaseline[i] - fitvalue[i])

                elif out[sig1][startpoints[i]] < 0 and (sig == 'i1' or sig == 'i2'):
                    if length <= limit and length > 3:
                        # Impulsion Fit to minimal value
                        fitvalue[i] = np.min(out[sig1][startpoints[i]+np.uint(1):endpoints[i]-np.uint(1)])
                    elif length > limit:
                        fitvalue[i] = np.mean(out[sig1][startpoints[i]+np.uint(5):endpoints[i]-np.uint(5)])
                    else:
                        fitvalue[i] = np.min(out[sig1][startpoints[i]:endpoints[i]])
                    deli[i] = -(localBaseline[i] - fitvalue[i])

                elif out[sig1][startpoints[i]] > 0 and (sig == 'i1_Up' or sig == 'i2_Up'):
                    if length <= limit and length > 3:
                        # Impulsion Fit to minimal value
                        fitvalue[i] = np.max(out[sig1][startpoints[i]+np.uint(1):endpoints[i]-np.uint(1)])
                    elif length > limit:
                        fitvalue[i] = np.mean(out[sig1][startpoints[i]+np.uint(5):endpoints[i]-np.uint(5)])
                    else:
                        fitvalue[i] = np.max(out[sig1][startpoints[i]:endpoints[i]])
                    deli[i] = (localBaseline[i] - fitvalue[i])

                elif out[sig1][startpoints[i]] > 0 and (sig == 'i1' or sig == 'i2'):
                    if length <= limit and length > 3:
                        # Impulsion Fit to minimal value
                        fitvalue[i] = np.min(out[sig1][startpoints[i]+np.uint(1):endpoints[i]-np.uint(1)])
                    elif length > limit:
                        fitvalue[i] = np.mean(out[sig1][startpoints[i]+np.uint(5):endpoints[i]-np.uint(5)])
                    else:
                        fitvalue[i] = np.min(out[sig1][startpoints[i]:endpoints[i]])
                    deli[i] = (localBaseline[i] - fitvalue[i])
                else:
                    print('Strange event that has to come from a neighboring universe...Computer will self-destruct in 3s!')
                dwell[i] = (endpoints[i] - startpoints[i]) / out['samplerate']
                AllEvents.append(out[sig1][startpoints[i]:endpoints[i]])
            frac = deli / localBaseline
            dt = np.array(0)
            dt = np.append(dt, np.diff(startpoints) / out['samplerate'])
            numberofevents = len(dt)
            #AnalysisResults[sig]['CusumFits'] = AllFits
            AnalysisResults[sig]['FractionalCurrentDrop'] = frac
            AnalysisResults[sig]['DeltaI'] = deli
            AnalysisResults[sig]['DwellTime'] = dwell
            AnalysisResults[sig]['Frequency'] = dt
            AnalysisResults[sig]['LocalVoltage'] = out[volt][startpoints]
            AnalysisResults[sig]['AllEvents'] = AllEvents
            AnalysisResults[sig]['FitLevel'] = fitvalue
        else:
            AnalysisResults[sig]['NumberOfEvents'] = 0

    return AnalysisResults

def CorrelateTheTwoChannels(AnalysisResults, DelayLimit, Ch1 = 'i1', Ch2 = 'i2'):
    if len(AnalysisResults[Ch1]['RoughEventLocations']) is not 0:
        i1StartP = np.int64(AnalysisResults[Ch1]['StartPoints'][:])
    else:
        i1StartP = []
    if len(AnalysisResults[Ch2]['RoughEventLocations']) is not 0:
        i2StartP = np.int64(AnalysisResults[Ch2]['StartPoints'][:])
    else:
        i2StartP = []

    # Common Events, # Take Longer
    CommonEventsi1Index = np.array([], dtype=np.uint64)
    CommonEventsi2Index = np.array([], dtype=np.uint64)

    for k in i1StartP:
        val = i2StartP[(i2StartP > k - DelayLimit) & (i2StartP < k + DelayLimit)]
        if len(val) == 1:
            CommonEventsi2Index = np.append(CommonEventsi2Index, np.uint64(np.where(i2StartP == val)[0][0]))
            CommonEventsi1Index = np.append(CommonEventsi1Index, np.uint64(np.where(i1StartP == k)[0][0]))
        if len(val) > 1:
            diff = np.absolute(val-k)
            CommonEventsi2Index = np.append(CommonEventsi2Index, np.uint64(np.where(i2StartP == val[np.argmin(diff)]))[0][0])
            CommonEventsi1Index = np.append(CommonEventsi1Index, np.uint64(np.where(i1StartP == k)[0][0]))
    # Only i1
    Onlyi1Indexes = np.delete(range(len(i1StartP)), CommonEventsi1Index)
    #Onlyi1Indexes=[]
    # Only i2
    Onlyi2Indexes = np.delete(range(len(i2StartP)), CommonEventsi2Index)
    #Onlyi2Indexes=[]

    CommonIndexes={}
    CommonIndexes[Ch1]=CommonEventsi1Index
    CommonIndexes[Ch2]=CommonEventsi2Index
    OnlyIndexes={}
    OnlyIndexes[Ch1] = Onlyi1Indexes
    OnlyIndexes[Ch2] = Onlyi2Indexes
    return (CommonIndexes, OnlyIndexes)

def PlotEvent(fig1, t1, i1, t2 = [], i2 = [], fit1 = np.array([]), fit2 = np.array([]), channel = 'i1'):
    if len(t2)==0:
        ax1 = fig1.add_subplot(111)
        ax1.plot(t1, i1*1e9, 'b')
        if len(fit1) is not 0:
            ax1.plot(t1, fit1*1e9, 'y')
        ax1.set_ylabel(channel + ' Current [nA]')
        ax1.set_xlabel(channel + ' Time [s]')
        ax1.ticklabel_format(useOffset=False)
        ax1.ticklabel_format(useOffset=False)
        return ax1
    else:
        ax1 = fig1.add_subplot(211)
        ax2 = fig1.add_subplot(212, sharex=ax1)
        ax1.plot(t1, i1*1e9, 'b')
        if len(fit1) is not 0:
            ax1.plot(t1, fit1*1e9, 'y')
        ax2.plot(t2, i2*1e9, 'r')
        if len(fit2) is not 0:
            ax2.plot(t2, fit2*1e9, 'y')
        ax1.set_ylabel('Ionic Current [nA]')
        #ax1.set_xticklabels([])
        ax2.set_ylabel('Transverse Current [nA]')
        ax2.set_xlabel('Time [s]')
        ax2.ticklabel_format(useOffset=False)
        ax2.ticklabel_format(useOffset=False)
        ax1.ticklabel_format(useOffset=False)
        ax1.ticklabel_format(useOffset=False)
        return ax1, ax2

def SaveAllPlots(CommonIndexes, OnlyIndexes, AnalysisResults, directory, out, buffer, withFit = 1):
    if len(CommonIndexes['i1']) is not 0:
        # Plot All Common Events
        pp = PdfPages(directory + '_SavedEventsCommon.pdf')
        ind1 = np.uint64(CommonIndexes['i1'])
        ind2 = np.uint64(CommonIndexes['i2'])

        t = np.arange(0, len(out['i1']))
        t = t / out['samplerate'] * 1e3
        count=1
        for eventnumber in range(len(ind1)):
            parttoplot = np.arange(AnalysisResults['i1']['StartPoints'][ind1[eventnumber]] - buffer,
                                   AnalysisResults['i1']['EndPoints'][ind1[eventnumber]] + buffer, 1, dtype=np.uint64)
            parttoplot2 = np.arange(AnalysisResults['i2']['StartPoints'][ind2[eventnumber]] - buffer,
                                    AnalysisResults['i2']['EndPoints'][ind2[eventnumber]] + buffer, 1, dtype=np.uint64)

            fit1 = np.concatenate([np.ones(buffer) * AnalysisResults['i1']['LocalBaseline'][ind1[eventnumber]],
                                   np.ones(AnalysisResults['i1']['EndPoints'][ind1[eventnumber]] - AnalysisResults['i1']['StartPoints'][
                                       ind1[eventnumber]]) * (
                                       AnalysisResults['i1']['LocalBaseline'][ind1[eventnumber]] - AnalysisResults['i1']['DeltaI'][ind1[eventnumber]]),
                                   np.ones(buffer) * AnalysisResults['i1']['LocalBaseline'][ind1[eventnumber]]])

            fit2 = np.concatenate([np.ones(buffer) * AnalysisResults['i2']['LocalBaseline'][ind2[eventnumber]],
                                   np.ones(AnalysisResults['i2']['EndPoints'][ind2[eventnumber]] - AnalysisResults['i2']['StartPoints'][
                                       ind2[eventnumber]]) * (
                                       AnalysisResults['i2']['LocalBaseline'][ind2[eventnumber]] - AnalysisResults['i2']['DeltaI'][ind2[eventnumber]]),
                                   np.ones(buffer) * AnalysisResults['i2']['LocalBaseline'][ind2[eventnumber]]])
            if withFit:
                fig = PlotEvent(t[parttoplot], out['i1'][parttoplot], t[parttoplot2], out['i2'][parttoplot2],
                                   fit1=fit1, fit2=fit2)
            else:
                fig = PlotEvent(t[parttoplot], out['i1'][parttoplot], t[parttoplot2], out['i2'][parttoplot2])

            if not np.mod(eventnumber+1,200):
                pp.close()
                pp = PdfPages(directory + '_SavedEventsCommon_' + str(count) + '.pdf')
                count+=1
            pp.savefig(fig)
            print('{} out of {} saved!'.format(str(eventnumber), str(len(ind1))))
            print('Length i1: {}, Fit i1: {}'.format(len(out['i1'][parttoplot]), len(fit1)))
            print('Length i2: {}, Fit i2: {}'.format(len(out['i2'][parttoplot2]), len(fit2)))
            fig.clear()
            plt.close(fig)
        pp.close()

    if len(OnlyIndexes['i1']) is not 0:
        # Plot All i1
        pp = PdfPages(directory + '_SavedEventsOnlyi1.pdf')
        ind1 = np.uint64(OnlyIndexes['i1'])

        t = np.arange(0, len(out['i1']))
        t = t / out['samplerate'] * 1e3
        count=1
        for eventnumber in range(len(ind1)):
            parttoplot = np.arange(AnalysisResults['i1']['StartPoints'][ind1[eventnumber]] - buffer,
                                   AnalysisResults['i1']['EndPoints'][ind1[eventnumber]] + buffer, 1, dtype=np.uint64)

            fit1 = np.concatenate([np.ones(buffer) * AnalysisResults['i1']['LocalBaseline'][ind1[eventnumber]],
                                   np.ones(AnalysisResults['i1']['EndPoints'][ind1[eventnumber]] - AnalysisResults['i1']['StartPoints'][
                                       ind1[eventnumber]]) * (
                                       AnalysisResults['i1']['LocalBaseline'][ind1[eventnumber]] - AnalysisResults['i1']['DeltaI'][ind1[eventnumber]]),
                                   np.ones(buffer) * AnalysisResults['i1']['LocalBaseline'][ind1[eventnumber]]])

            fig = PlotEvent(t[parttoplot], out['i1'][parttoplot], t[parttoplot], out['i2'][parttoplot], fit1=fit1)
            if not np.mod(eventnumber+1,200):
                pp.close()
                pp = PdfPages(directory + '_SavedEventsCommon_' + str(count) + '.pdf')
                count+=1
            pp.savefig(fig)
            print('{} out of {} saved!'.format(str(eventnumber), str(len(ind1))))
            fig.clear()
            plt.close(fig)
        pp.close()

    if len(OnlyIndexes['i2']) is not 0:
        # Plot All i2
        pp = PdfPages(directory + '_SavedEventsOnlyi2.pdf')
        ind1 = np.uint64(OnlyIndexes['i2'])

        t = np.arange(0, len(out['i2']))
        t = t / out['samplerate'] * 1e3
        count=1
        for eventnumber in range(len(ind1)):
            parttoplot = np.arange(AnalysisResults['i2']['StartPoints'][ind1[eventnumber]] - buffer,
                                   AnalysisResults['i2']['EndPoints'][ind1[eventnumber]] + buffer, 1, dtype=np.uint64)

            fit1 = np.concatenate([np.ones(buffer) * AnalysisResults['i2']['LocalBaseline'][ind1[eventnumber]],
                                   np.ones(AnalysisResults['i2']['EndPoints'][ind1[eventnumber]] - AnalysisResults['i2']['StartPoints'][
                                       ind1[eventnumber]]) * (
                                       AnalysisResults['i2']['LocalBaseline'][ind1[eventnumber]] - AnalysisResults['i2']['DeltaI'][ind1[eventnumber]]),
                                   np.ones(buffer) * AnalysisResults['i2']['LocalBaseline'][ind1[eventnumber]]])

            fig = PlotEvent(t[parttoplot], out['i1'][parttoplot], t[parttoplot], out['i2'][parttoplot], fit2=fit1)
            if not np.mod(eventnumber+1,200):
                pp.close()
                pp = PdfPages(directory + '_SavedEventsCommon_' + str(count) + '.pdf')
                count+=1
            pp.savefig(fig)
            print('{} out of {} saved!'.format(str(eventnumber), str(len(ind1))))
            fig.clear()
            plt.close(fig)
        pp.close()

    # Derivative
    if len(CommonIndexes['i1']) is not 0:
        # Plot All i1
        pp = PdfPages(directory + '_i1vsderivi2.pdf')
        ind1 = np.uint64(CommonIndexes['i1'])
        ind2 = np.uint64(CommonIndexes['i2'])

        t = np.arange(0, len(out['i1']))
        t = t / out['samplerate'] * 1e3
        count=1
        for eventnumber in range(len(ind1)):
            parttoplot = np.arange(AnalysisResults['i1']['StartPoints'][ind1[eventnumber]] - buffer,
                                   AnalysisResults['i1']['EndPoints'][ind1[eventnumber]] + buffer, 1, dtype=np.uint64)
            parttoplot2 = np.arange(AnalysisResults['i2']['StartPoints'][ind2[eventnumber]] - buffer,
                                    AnalysisResults['i2']['EndPoints'][ind2[eventnumber]] + buffer, 1, dtype=np.uint64)

            fig = PlotEvent(t[parttoplot], out['i1'][parttoplot], t[parttoplot2][:-1],
                               np.diff(out['i2'][parttoplot2]))

            if not np.mod(eventnumber+1,200):
                pp.close()
                pp = PdfPages(directory + '_SavedEventsCommon_' + str(count) + '.pdf')
                count+=1
            pp.savefig(fig)
            print('{} out of {} saved!'.format(str(eventnumber), str(len(ind1))))
            fig.clear()
            plt.close(fig)
        pp.close()

def PlotRecursiveLPResults(RoughEventLocations, inp, directory, buffer, channel='i2'):
    pp = PdfPages(directory + '_' + channel + '_DetectedEventsFromLPFilter.pdf')
    a=1
    for i in RoughEventLocations['RoughEventLocations']:
        startp = np.uint64(i[0]-buffer*inp['samplerate'])
        endp = np.uint64(i[1]+buffer*inp['samplerate'])
        t = np.arange(startp, endp)
        t = t / inp['samplerate'] * 1e3
        fig = PlotEvent(t, inp[channel][startp:endp], channel=channel)
        pp.savefig(fig)
        print('{} out of {} saved!'.format(str(a), str(len(RoughEventLocations['RoughEventLocations']))))
        a+=1
        fig.clear()
        plt.close(fig)
    pp.close()

def SaveAllAxopatchEvents(AnalysisResults, directory, out, buffer, withFit = 1):
    # Plot All Common Events
    pp = PdfPages(directory + '_SavedEventsAxopatch.pdf')
    t = np.arange(0, len(out['i1']))
    t = t / out['samplerate'] * 1e3

    for eventnumber in range(AnalysisResults['i1']['NumberOfEvents']):
        parttoplot = np.arange(AnalysisResults['i1']['StartPoints'][eventnumber] - buffer,
                               AnalysisResults['i1']['EndPoints'][eventnumber] + buffer, 1, dtype=np.uint64)

        fit1 = np.concatenate([np.ones(buffer) * AnalysisResults['i1']['LocalBaseline'][eventnumber],
                               np.ones(AnalysisResults['i1']['EndPoints'][eventnumber] -
                                       AnalysisResults['i1']['StartPoints'][
                                           eventnumber]) * (
                                   AnalysisResults['i1']['LocalBaseline'][eventnumber] -
                                   AnalysisResults['i1']['DeltaI'][eventnumber]),
                               np.ones(buffer) * AnalysisResults['i1']['LocalBaseline'][eventnumber]])

        if withFit:
            fig = PlotEvent(t[parttoplot], out['i1'][parttoplot], fit1=fit1)
        else:
            fig = PlotEvent(t[parttoplot], out['i1'][parttoplot])
        try:
            pp.savefig(fig)
        except:
            print('Problem at {} !'.format(str(eventnumber)))

        print('{} out of {} saved!'.format(str(eventnumber), str(AnalysisResults['i1']['NumberOfEvents'])))
        #print('Length i1: {}, Fit i1: {}'.format(len(out['i1'][parttoplot]), len(fit1)))
        fig.clear()
        plt.close(fig)
    pp.close()

def MakeIV(filenames, directory, conductance=10, title=''):
    if len(conductance) is 1:
        conductance=np.ones(len(filenames))*conductance
    for idx,filename in enumerate(filenames):
        print(filename)

        # Make Dir to save images
        output = OpenFile(filename)
        if not os.path.exists(directory):
            os.makedirs(directory)

        AllData = MakeIVData(output, delay=2)
        if AllData == 0:
            print('!!!! No Sweep in: ' + filename)
            continue

        # Plot Considered Part
        # figExtracteParts = plt.figure(1)
        # ax1 = figExtracteParts.add_subplot(211)
        # ax2 = figExtracteParts.add_subplot(212, sharex=ax1)
        # (ax1, ax2) = uf.PlotExtractedPart(output, AllData, current = 'i1', unit=1e9, axis = ax1, axis2=ax2)
        # plt.show()
        # figExtracteParts.savefig(directory + os.sep + 'PlotExtracted_' + str(os.path.split(filename)[1])[:-4] + '.eps')
        # figExtracteParts.savefig(directory + os.sep + 'PlotExtracted_' + str(os.path.split(filename)[1])[:-4] + '.png', dpi=150)


        # Plot IV
        if output['graphene']:
            figIV2 = plt.figure(3)
            ax2IV = figIV2.add_subplot(111)
            ax2IV = PlotIV(output, AllData, current='i2', unit=1e9, axis=ax2IV, WithFit=1, title=title[idx])
            figIV2.tight_layout()
            figIV2.savefig(directory + os.sep + str(os.path.split(filename)[1]) + '_IV_i2.png', dpi=300)
            figIV2.savefig(directory + os.sep + str(os.path.split(filename)[1]) + '_IV_i2.eps')
            figIV2.clear()

        figIV = plt.figure(2)
        ax1IV = figIV.add_subplot(111)
        ax1IV = PlotIV(output, AllData, current='i1', unit=1e9, axis=ax1IV, WithFit=1, PoreSize=[conductance[idx], 20e-9], title=title[idx])
        figIV.tight_layout()

        # Save Figures
        figIV.savefig(directory + os.sep + str(os.path.split(filename)[1]) + 'IV_i1.png', dpi=300)
        figIV.savefig(directory + os.sep + str(os.path.split(filename)[1]) + 'IV_i1.eps')
        figIV.clear()

def TwoChannelAnalysis(filenames, outdir, UpwardsOn=0):

    for filename in filenames:
        print(filename)
        inp = OpenFile(filename)
        file = os.sep + str(os.path.split(filename)[1][:-4])

        if not os.path.exists(outdir):
            os.makedirs(outdir)

        # Low Pass Event Detection
        AnalysisResults = {}

        if inp['graphene']:
            chan = ['i1', 'i2']  # For looping over the channels
        else:
            chan = ['i1']

        for sig in chan:
            AnalysisResults[sig] = {}
            AnalysisResults[sig]['RoughEventLocations'] = RecursiveLowPassFast(inp[sig], pm.coefficients[sig],
                                                                                 inp['samplerate'])
            if UpwardsOn:  # Upwards detection can be turned on or off
                AnalysisResults[sig + '_Up'] = {}
                AnalysisResults[sig + '_Up']['RoughEventLocations'] = RecursiveLowPassFastUp(inp[sig],
                                                                                               pm.coefficients[sig],
                                                                                               inp['samplerate'])

        # f.PlotRecursiveLPResults(AnalysisResults['i2'], inp, directory, 1e-3, channel='i2')
        # f.PlotRecursiveLPResults(AnalysisResults['i2_Up'], inp, directory+'UP_', 1e-3, channel='i2')

        # Refine the Rough Event Detection done by the LP filter and Add event infos
        AnalysisResults = RefinedEventDetection(inp, AnalysisResults, signals=chan,
                                                  limit=pm.MinimalFittingLimit * inp['samplerate'])

        # Correlate the two channels
        if inp['graphene']:
            (CommonIndexes, OnlyIndexes) = f.CorrelateTheTwoChannels(AnalysisResults, 10e-3 * inp['samplerate'])
            print('\n\nAnalysis Done...\nThere are {} common events\n{} Events on i1 only\n{} Events on i2 only'.format(
                len(CommonIndexes['i1']), len(OnlyIndexes['i1']), len(OnlyIndexes['i2'])))
            # Plot The Events
            f.SaveAllPlots(CommonIndexes, OnlyIndexes, AnalysisResults, directory, inp, pm.PlotBuffer)
        else:
            # Plot The Events
            f.SaveAllAxopatchEvents(AnalysisResults, directory, inp, pm.PlotBuffer)

def SaveToHDF5(inp_file, AnalysisResults, coefficients, outdir):
    file = str(os.path.split(inp_file['filename'])[1][:-4])
    f = h5py.File(outdir + file + '_OriginalDB.hdf5', "w")
    general = f.create_group("General")
    general.create_dataset('FileName', data=inp_file['filename'])
    general.create_dataset('Samplerate', data=inp_file['samplerate'])
    general.create_dataset('Machine', data=inp_file['type'])
    general.create_dataset('TransverseRecorded', data=inp_file['graphene'])
    general.create_dataset('TimeFileWritten', data=inp_file['TimeFileWritten'])
    general.create_dataset('TimeFileLastModified', data=inp_file['TimeFileLastModified'])
    general.create_dataset('ExperimentDuration', data=inp_file['ExperimentDuration'])

    segmentation_LP = f.create_group("LowPassSegmentation")
    for k,l in AnalysisResults.items():
        set1 = segmentation_LP.create_group(k)
        lpset1 = set1.create_group('LowPassSettings')
        for o, p in coefficients[k].items():
             lpset1.create_dataset(o, data=p)
        for m, l in AnalysisResults[k].items():
            if m is 'AllEvents':
                eventgroup = set1.create_group(m)
                for i, val in enumerate(l):
                    eventgroup.create_dataset('{:09d}'.format(i), data=val)
            elif m is 'Cusum':
                eventgroup = set1.create_group(m)
                for i1, val1 in enumerate(AnalysisResults[k]['Cusum']):
                    cusevent = eventgroup.create_group('{:09d}'.format(i1))
                    cusevent.create_dataset('NumberLevels', data=np.uint64(len(AnalysisResults[k]['Cusum'][i1]['levels'])))
                    if len(AnalysisResults[k]['Cusum'][i1]['levels']):
                        cusevent.create_dataset('up', data=AnalysisResults[k]['Cusum'][i1]['up'])
                        cusevent.create_dataset('down', data=AnalysisResults[k]['Cusum'][i1]['down'])
                        cusevent.create_dataset('both', data=AnalysisResults[k]['Cusum'][i1]['both'])
                        cusevent.create_dataset('fit', data=AnalysisResults[k]['Cusum'][i1]['fit'])
                        # 0: level number, 1: current, 2: length, 3: std
                        cusevent.create_dataset('levels_current', data=AnalysisResults[k]['Cusum'][i1]['levels'][1])
                        cusevent.create_dataset('levels_length', data=AnalysisResults[k]['Cusum'][i1]['levels'][2])
                        cusevent.create_dataset('levels_std', data=AnalysisResults[k]['Cusum'][i1]['levels'][3])
            else:
                set1.create_dataset(m, data=l)

def PlotExtractedPart(output, AllData, current = 'i1', unit=1e9, axis = '', axis2 = '', exp='_Exp'):
    time = np.arange(0, len(output[current])) / output['samplerate']
    axis.plot(time, output[current] * unit, 'b', label=str(os.path.split(output['filename'])[1])[:-4])
    axis.set_ylabel('Current ' + current + ' [nA]')
    axis.set_title('Time Trace')
    for i in range(0, len(AllData[current+exp]['StartPoint'])):
        x = time[AllData[current+exp]['StartPoint'][i]:AllData[current+exp]['EndPoint'][i]]
        axis.plot(x,output[current][AllData[current+exp]['StartPoint'][i]:AllData[current+exp]['EndPoint'][i]] * unit, 'r')
        axis.plot(x, unit*ExpFunc(x-np.min(x), AllData[current+exp]['ExpValues'][0][i], AllData[current+exp]['ExpValues'][1][i], AllData[current+exp]['ExpValues'][2][i]),'y')
    axis2.plot(time, output[AllData[current+exp]['SweepedChannel']], 'b', label=str(os.path.split(output['filename'])[1])[:-4])
    axis2.set_ylabel('Voltage ' + AllData[current+exp]['SweepedChannel'] + ' [V]')
    axis2.set_xlabel('Time')
    return (axis, axis2)

def xcorr(x, y, k, normalize=True):
    n = x.shape[0]

    # initialize the output array
    out = np.empty((2 * k) + 1, dtype=np.double)
    lags = np.arange(-k, k + 1)

    # pre-compute E(x), E(y)
    mu_x = x.mean()
    mu_y = y.mean()

    # loop over lags
    for ii, lag in enumerate(lags):

        # use slice indexing to get 'shifted' views of the two input signals
        if lag < 0:
            xi = x[:lag]
            yi = y[-lag:]
        elif lag > 0:
            xi = x[:-lag]
            yi = y[lag:]
        else:
            xi = x
            yi = y

        # x - mu_x; y - mu_y
        xdiff = xi - mu_x
        ydiff = yi - mu_y

        # E[(x - mu_x) * (y - mu_y)]
        out[ii] = xdiff.dot(ydiff) / n

        # NB: xdiff.dot(ydiff) == (xdiff * ydiff).sum()

    if normalize:
        # E[(x - mu_x) * (y - mu_y)] / (sigma_x * sigma_y)
        out /=  np.std(x) * np.std(y)

    return lags, out


def CUSUM(input, delta, h):
    Nd = 0
    kd = []
    krmv = []
    k0 = 0
    k = 1
    l = len(input)
    m = np.zeros(l)
    m[k0] = input[k0]
    v = np.zeros(l)
    sp = np.zeros(l)
    Sp = np.zeros(l)
    gp = np.zeros(l)
    sn = np.zeros(l)
    Sn = np.zeros(l)
    gn = np.zeros(l)

    while k < l:
        m[k] = np.mean(input[k0:k + 1])
        v[k] = np.var(input[k0:k + 1])

        sp[k] = delta / v[k] * (input[k] - m[k] - delta / 2)
        sn[k] = -delta / v[k] * (input[k] - m[k] + delta / 2)

        Sp[k] = Sp[k - 1] + sp[k]
        Sn[k] = Sn[k - 1] + sn[k]

        gp[k] = np.max([gp[k - 1] + sp[k], 0])
        gn[k] = np.max([gn[k - 1] + sn[k], 0])

        if gp[k] > h or gn[k] > h:
            kd.append(k)
            if gp[k] > h:
                kmin = np.argmin(Sp[k0:k + 1])
                krmv.append(kmin + k0 - 1)
            else:
                kmin=np.argmin(Sn[k0:k+1])
                krmv.append(kmin + k0 - 1)

            #Re-initialize
            k0 = k
            m[k0] = input[k0]
            v[k0] = sp[k0] = Sp[k0] = gp[k0] = sn[k0] = Sn[k0] = gn[k0] = 0

            Nd = Nd + 1
        k += 1

    if Nd == 0:
        mc = np.mean(input) * np.ones(k)
    elif Nd == 1:
        mc = np.append(m[krmv[0]] * np.ones(krmv[0]), m[k] * np.ones(k - krmv[0]))
    else:
        mc = m[krmv[0]] * np.ones(krmv[0])
        for ii in range(1, Nd):
            mc=np.append(mc , m[krmv[ii]] * np.ones(krmv[ii] - krmv[ii - 1]))
        mc=np.append(mc , m[k-1] * np.ones(k - krmv[Nd-1]))
    return (mc, kd, krmv)

def SetCusum2ARL0(deltax,sigmax,Arl0_2=1000,thresholdlevels=[1000,0.1,10]):
    h0min=thresholdlevels[0];
    h0max=thresholdlevels[1]; #detection threshold interval

    ARL0=2*ARL0_2; #for two-sided algo
    # Optimization on h "hopt"
    change=0;
    mu=deltax*(change-deltax/2)/sigmax**2;
    sigma=abs(deltax)/sigmax;
    mun=mu/sigma;
    def f(h):
        (exp(-2*mun*(h/sigma+1.166))-1+2*mun*(h/sigma+1.166))/(2*mun^2)-ARL0

    if(f(h0min)*f(h0max)<0):
        print('test')
        #hopt=fzero(f,[h0min h0max])
    elif(f(h0min) <0 and f(h0max) <0):
        hopt=h0max
    elif(f(h0min) >0 and f(h0max) >0):
        hopt=h0min

    hbook=sigmax*hopt/deltax;
    return hbook, hopt


def event_detection(RawSignal, CusumParameters, RoughEventLocations, cutoff):
    for i in range(len(RoughEventLocations)):
        CusumReferencedEndPoint = RoughEventLocations[i][1] + CusumParameters['BaselineLength'];
        CusumReferencedStartPoint = RoughEventLocations[i][0] - CusumParameters['BaselineLength'];

        if CusumReferencedEndPoint > len(RawSignal):
            CusumReferencedEndPoint = len(RawSignal)
        if CusumReferencedStartPoint < 0:
            CusumReferencedStartPoint = 0

        if RoughEventLocations[i][2] < CusumParameters['ImpulsionLimit'] * CusumParameters['SamplingFrequency']:
            trace = RawSignal[CusumReferencedStartPoint:CusumReferencedEndPoint]
             #[mc,krmv]=ImpulsionFitting(trace,BaselineLength, cutoff, SamplingFrequency);
             #EventType='Impulsion';
             #krmv=[BaselineLength, BaselineLength+RoughEventLocations(i,3)];
             #mc=[ones(1,BaselineLength+1)*mean(trace(1:BaselineLength)), ones(1,krmv(2)-krmv(1)-1)*min(trace), ones(1,BaselineLength-2)*mean(trace(krmv(2):end))];
        else:
            mc, kd, krmv = cusum(RawSignal[CusumReferencedStartPoint:CusumReferencedEndPoint],CusumParameters['delta'],CusumParameters['h'])
            EventType='Standard Event';

        Startpoint = CusumReferencedStartPoint + krmv[0]
        EndPoint = CusumReferencedStartPoint + krm[-1]

        BaselinePart = RawSignal[CusumReferencedStartPoint:CusumReferencedStartPoint]


def CondPoreSizeTick(x, pos):
    formCond = EngFormatter(unit='S')
    formPoreSize = EngFormatter(unit='m')
    outstring = '{}, {}'.format(formCond.format_eng(x), formPoreSize.format_eng(CalculatePoreSize(x, 1e-9, 10)))
    return outstring

def EredoxBefore14062018():
    Eredox = np.array([31.7, 82.9, 135, 185], dtype=float)
    Eredox = Eredox * 1e-3
    return (Eredox-Eredox[0])

def DebyeHuckelParam(c, T = 20):
    ## Valid Only for same valency and c1 = c2
    epsilon_r = 87.740-0.40008*T+9.398e-4*T**2-1.410e-6*T**3
    print('Dielectric of Water is {}'.format(epsilon_r))
    n = c*1e3*cst.Avogadro
    formCond = EngFormatter(unit='m')
    k = np.sqrt((cst.e**2*2*n)/(cst.k*(T+273.15)*cst.epsilon_0*epsilon_r))
    return [k, 1/k]#, '{}'.format(formCond.format_eng(1/k))]

def ElectrostaticPotential(T, sigma, Tau=5e18, pK=7.9, pH=11):
    return cst.k * T / cst.e * (np.log(-sigma / (cst.e * Tau + sigma)) + (pK - pH) * np.log(10))

def GetTempRedoxPotential(Cdiff=1e-3/1e-1, T=293.15):
    return cst.R * T / cst.physical_constants['Faraday constant'][0] * np.log(Cdiff)

def GetRedoxError(cmax, cmin, ErrorPercent = 0.1, T = 293.15):
    cmaxErr = cmax*ErrorPercent
    cminErr = cmin*ErrorPercent
    return cst.R * T / cst.physical_constants['Faraday constant'][0] * np.sqrt((1/cmax**2 * cmaxErr**2 + 1/cmin**2 * cminErr**2))

##Result: Negative ion divided by positive ion
#Goldman–Hodgkin–Katz voltage equation
def GetIonSelectivityWithoutPotential(c_trans, c_cis, Erev, T):
    n_trans = c_trans#*1e3#*cst.Avogadro
    n_cis = c_cis#*1e3#*cst.Avogadro
    A = Erev * cst.physical_constants['Faraday constant'][0] / (cst.R*T)
    return -(n_trans-n_cis*np.exp(-A))*(1-np.exp(A))/((n_trans-n_cis*np.exp(A))*(1-np.exp(-A)))

def GetIonSelectivityWithPotential(c_trans, c_cis, Erev, T, phi):
    sel1 = GetIonSelectivityWithoutPotential(c_trans, c_cis, Erev, T)
    return sel1*np.exp((-2*phi*cst.physical_constants['Faraday constant'][0])/(cst.R * T))

def Selectivity(ConcGradient, V):
    return V * cst.physical_constants['Faraday constant'][0]/((cst.R*293.15) * np.log(ConcGradient))

def GetRedox(cmax, cmin, T = 295.15):
    return cst.R * T / cst.physical_constants['Faraday constant'][0] * np.log(cmax/cmin)
