Activity-mediated_accumulation_potassium_induces_switch_firing_pattern/functions/f_post_simulation_analysis.py

import numpy as np
from scipy.signal import find_peaks_cwt
from scipy.signal import argrelextrema
from scipy.optimize import curve_fit

################important define!
global precision_convergence_ss

precision_convergence_ss=0.00001

def tic():
    #Homemade version of matlab tic and toc functions
    import time
    global startTime_for_tictoc
    startTime_for_tictoc = time.time()

def toc():
    import time
    if 'startTime_for_tictoc' in globals():
        print( "Elapsed time is " + str(time.time() - startTime_for_tictoc) + " seconds.")
    else:
        print( "Toc: start time not set")

def detect_peaks(signal, threshold=0.5, spacing=3):#https://github.com/MonsieurV/py-findpeaks/blob/master/tests/libs/tony_beltramelli_detect_peaks.py
    limit=None
    data=signal
    len = data.size
    x = np.zeros(len+2*spacing)
    x[:spacing] = data[0]-1.e-6
    x[-spacing:] = data[-1]-1.e-6
    x[spacing:spacing+len] = data
    peak_candidate = np.zeros(len)
    peak_candidate[:] = True
    for s in range(spacing):
        start = spacing - s - 1
        h_b = x[start : start + len]  # before
        start = spacing
        h_c = x[start : start + len]  # central
        start = spacing + s + 1
        h_a = x[start : start + len]  # after
        peak_candidate = np.logical_and(peak_candidate, np.logical_and(h_c > h_b, h_c > h_a))

    peak_indexes = np.argwhere(peak_candidate)
    peak_indexes = peak_indexes.reshape(peak_indexes.size)
    if limit is not None:
        peak_indexes = ind[data[ind] > limit]
    return peak_indexes

def firing_rate(t,V):
    fr=np.zeros(size(t))
    t_sp=np.zeros(size(t))
    d=[]
    sp=[]
    i_n=[]
    d=detect_peaks(V)
    if d==[] or len(d)<2:
        fr=np.zeros(size(t))
        t_sp=np.zeros(size(t))
    else:
        sp=V[d]>thresh
        if np.count_nonzero(sp)>2:
            i=np.transpose(sp)*d
            i_n=i[np.nonzero(i)]
            t_sp=t[i_n]
            t_sp=t_sp.reshape(size(t_sp),1)
            fr=1/(np.diff(t_sp, n=1, axis=0)/1000)
            t_sp=t_sp[1:]
        else:
            fr=np.zeros(size(t))
            t_sp=t

    return t_sp,fr,i_n

def extract_spike_ampl(V, i_nn,compress=[]):
    c=0
    sp_V_max=[]
    sp_V_min=[]
    sp_ampl=[]
    if compress !=[]:
        i_n=[int(i/compress) for i in i_nn]
    else:
        i_n=i_nn
    d_i_ni=int(i_n[0]/2)

    for i in i_n:
        if c<len(i_n)-1:
            d_i_nd=int((i_n[c+1]-i)/2)
        else:
            d_i_nd=int((len(V)-i)/2)
        sp_V_max.append(max(V[i-d_i_ni:i+d_i_nd]))
        sp_V_min.append(min(V[i-d_i_ni:i+d_i_nd]))
        sp_ampl.append(max(V[i-d_i_ni:i+d_i_nd])-min(V[i-d_i_ni:i+d_i_nd]))
        c+=1
        d_i_ni=int((i_n[c-1]-i)/2)

    return sp_ampl,sp_V_max,sp_V_min 



def doub_expon_fun_2(t,tau1,tau2,a,b,c):
    if b > 0 and a > 0 and c>0 and tau1>0 and tau2>0:## Constraint to have all values positive
        return a * np.exp(-t/tau1 ) + b * np.exp(-t/tau2 )+c
    else:
        return 1e10
    
def doub_expon_fit_2(t,x,p=[],std_error=False):
    x=np.array(x)
    t=np.array(t)
    perr=[]
    im=np.argmax(x)
    popt=[]
    tad=[]
    tpeak=[]
    adapt=False
    if x[-1]<x[im]:
        adapt=True
    if adapt:
        tad=t[im:]-t[im]
        frad=x[im:]
        frad02=x[im:]
        frad01=x[im:int(len(x[im:])/10)]
        half02=(max(frad02)-min(frad02))/2
        half01=(max(frad01)-min(frad01))/2
        halfid01=(np.abs(frad01-(max(frad01)-half01))).argmin()
        halfid02=(np.abs(frad02-(max(frad02)-half02))).argmin()
        ss_data01=frad[-int(len(frad01)/5):]
        ss_data02=frad[-int(len(frad02)/5):]
        ss_data01=ss_data01[abs(ss_data01 - np.mean(ss_data01)) < 2 * np.std(ss_data01)]
        ss_data02=ss_data02[abs(ss_data02 - np.mean(ss_data02)) < 2 * np.std(ss_data02)]
        if p!=[]:
            popt,pcov=curve_fit(doub_expon_fun_2,tad.flatten(),frad.flatten(), p0=p)
        else:
            popt,pcov=curve_fit(doub_expon_fun_2,tad.flatten(),frad.flatten(), p0=[tad[halfid01],tad[halfid02]*10,x[im]-mean(ss_data01[-int(len(ss_data01)/5):]),(x[im]-mean(ss_data02[-int(len(ss_data02)/5):])),mean(ss_data02[-int(len(ss_data02)/5):])])

        tpeak=t[im]
        perr=np.sqrt(np.diag(pcov))
    if std_error:
        return tad, popt,tpeak, perr
    else:
        return tad, popt,tpeak