123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140 |
- 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
|