import math import copy def spike_times(time,vrec,V_min_peak=-20,V_max_valley=0): valley_reached = 1 sptime = [] for j in range(1,len(time)-1): if valley_reached and vrec[j] >= V_min_peak and vrec[j] > vrec[j-1] and vrec[j] >= vrec[j+1]: valley_reached = 0 sptime.append(time[j]) elif valley_reached==False and vrec[j] <= V_max_valley: valley_reached = 1 return sptime #membpotderivs(time,vrec): Given the membrane potentials (vrec) at time points time[0],time[1],...,time[N-1], #return the derivatives at time points time[1],time[2],...,time[N-3] def membpotderivs(time,vrec): N = len(time) tdiff = [x-y for x,y in zip(time[1:N-1],time[0:N-2])] vdiff = [x-y for x,y in zip(vrec[1:N-1],vrec[0:N-2])] mderiv = [x/y for x,y in zip(vdiff,tdiff)] return [0.5*(x+y) for x,y in zip(mderiv[1:N-2],mderiv[0:N-3])] #limitcyclescaledv(v1,dv1,v2,dv2): Give the coefficient for memb. pot. derivative that one has to use in order to make #the difference on the derivative axis as significant as the difference on the memb. pot. axis def limitcyclescaledv(v1,dv1,v2,dv2): maxv = max(max(v1),max(v2)) minv = min(min(v1),min(v2)) maxdv = max(max(dv1),max(dv2)) mindv = min(min(dv1),min(dv2)) return 1.0*(maxv-minv)/(maxdv-mindv) def limitcyclediff(v1,dv1,v2,dv2,dvcoeff=0.1): N1 = len(v1) N2 = len(v2) dv1 = [dvcoeff*x for x in dv1] dv2 = [dvcoeff*x for x in dv2] Dmin = N1*[0] for i in range(0,N1): Dmin[i] = math.sqrt(min([(x-v1[i])**2+(y-dv1[i])**2 for x,y in zip(v2,dv2)])) vdiff = [x-y for x,y in zip(v1[1:N1],v1[0:N1-1])] dvdiff = [x-y for x,y in zip(dv1[1:N1],dv1[0:N1-1])] h = [math.sqrt(x**2+y**2) for x,y in zip(vdiff,dvdiff)] #use the trapezoid rule for integration: Dminmean = [(x+y)/2.0 for x,y in zip(Dmin[1:N1],Dmin[0:N1-1])] #print "hsum="+str(sum(h)) return sum([x*y for x,y in zip(Dminmean,h)]) def interpolate(tref,vref,tint): #Assumes that the trefs come sorted! vint = len(tint)*[0.0] addedOne = False #print tref #print tint #if tref[len(tref)-1] == tint[len(tint)-1]: # tref.append(tref[len(tref)-1]+0.0001) # vref.append(vref[len(tref)-1]) # addedOne = True if tref[0] > tint[0] or tref[len(tref)-1] < tint[len(tint)-1]: print "Extrapolation needed!" return len(tint)*[-1] indvrecnow = 0 for j in range(0,len(tint)): while tref[indvrecnow+1] <= tint[j]: indvrecnow = indvrecnow + 1 if indvrecnow == len(tref)-1: # It must be the last index if this happens vint[j:len(tint)] = [vref[indvrecnow]]*(len(tint)-j) return vint vint[j] = vref[indvrecnow] + 1.0*(tint[j]-tref[indvrecnow])/(tref[indvrecnow+1]-tref[indvrecnow])*(vref[indvrecnow+1]-vref[indvrecnow]) return vint def interpolate_extrapolate_constant(tref,vref,tint): #Assumes that the trefs come sorted! vint = len(tint)*[0.0] addedOne = False #print tref #print tint #if tref[len(tref)-1] == tint[len(tint)-1]: # tref.append(tref[len(tref)-1]+0.0001) # vref.append(vref[len(tref)-1]) # addedOne = True indvrecnow = 0 for j in range(0,len(tint)): if tint[j] < tref[0]: vint[j] = vref[0] continue if indvrecnow >= len(tref) - 1: vint[j] = vref[-1] while tref[indvrecnow+1] <= tint[j]: indvrecnow = indvrecnow + 1 if indvrecnow == len(tref)-1: # It must be the last index if this happens vint[j:len(tint)] = [vref[indvrecnow]]*(len(tint)-j) return vint vint[j] = vref[indvrecnow] + 1.0*(tint[j]-tref[indvrecnow])/(tref[indvrecnow+1]-tref[indvrecnow])*(vref[indvrecnow+1]-vref[indvrecnow]) return vint #kronecker product of list A and list B def kron(A,B): C = [] if type(B[0]) is int or type(B[0]) is float: for i in range(0,len(A)): for j in range(0,len(B)): print "asdf" print B[j] C.append(A[i]*B[j]) elif type(B[0][0]) is int or type(B[0][0]) is float: for i in range(0,len(A)): for j in range(0,len(B)): C.append([x*A[i] for x in B[j]]) return C def cumprod(A): B = len(A)*[0]; B[0]=A[0] for j in range(1,len(A)): B[j] = B[j-1]*A[j] return B def printlistlen(A): #TODO: recursive method might work out but needs some thought... #toCheck = A #lens = [len(x) for x in A] #levels = [0 for x in A] #while type(toCheck) is list and len(toCheck) > 0: # while type(toCheck[0]) is list and len(toCheck[0]) > 0: # toCheck.append(toCheck[0][0]) # lens.append[len # toCheck[0].pop(0) nan = -1 if type(A) is list: B = copy.deepcopy(A) listFound0 = False for i0 in range(0,len(B)): if type(B[i0]) is list: listFound0 = True listFound1 = False for i1 in range(0,len(B[i0])): if type(B[i0][i1]) is list: listFound1 = True listFound2 = False for i2 in range(0,len(B[i0][i1])): if type(B[i0][i1][i2]) is list: listFound2 = True listFound3 = False for i3 in range(0,len(B[i0][i1][i2])): if type(B[i0][i1][i2][i3]) is list: listFound3 = True listFound4 = False for i4 in range(0,len(B[i0][i1][i2][i3])): if type(B[i0][i1][i2][i3][i4]) is list: listFound4 = True listFound5 = False for i5 in range(0,len(B[i0][i1][i2][i3][i4])): if type(B[i0][i1][i2][i3][i4][i5]) is list: listFound5 = True listFound6 = False for i6 in range(0,len(B[i0][i1][i2][i3][i4][i5])): if type(B[i0][i1][i2][i3][i4][i5][i6]) is list: listFound6 = True B[i0][i1][i2][i3][i4][i5][i6] = len(B[i0][i1][i2][i3][i4][i5][i6]) else: B[i0][i1][i2][i3][i4][i5][i6] = nan if not listFound6: B[i0][i1][i2][i3][i4][i5] = len(B[i0][i1][i2][i3][i4][i5]) else: B[i0][i1][i2][i3][i4][i5] = nan if not listFound5: B[i0][i1][i2][i3][i4] = len(B[i0][i1][i2][i3][i4]) else: B[i0][i1][i2][i3][i4] = nan if not listFound4: B[i0][i1][i2][i3] = len(B[i0][i1][i2][i3]) else: B[i0][i1][i2][i3] = nan if not listFound3: B[i0][i1][i2] = len(B[i0][i1][i2]) else: B[i0][i1][i2] = nan if not listFound2: B[i0][i1] = len(B[i0][i1]) else: B[i0][i1] = nan if not listFound1: B[i0] = len(B[i0]) else: B[i0] = nan if not listFound0: B = len(B) else: B = nan print B