swhlab/analysis/ap.py
def detectSweep(self,sweep=0):
        """perform AP detection on current sweep."""

        if self.APs is False: # indicates detection never happened
            self.APs=[] # now indicates detection occured

        # delete every AP from this sweep from the existing array
        for i,ap in enumerate(self.APs):
            if ap["sweep"]==sweep:
                self.APs[i]=None
        if self.APs.count(None):
            self.log.debug("deleting %d existing APs from memory",self.APs.count(None))
            while None in self.APs:
                self.APs.remove(None)
        self.log.debug("initiating AP detection (%d already in memory)",len(self.APs))

        self.abf.derivative=True
        self.abf.setsweep(sweep)

        # detect potential AP (Is) by a dV/dT threshold crossing
        Is = cm.where_cross(self.abf.sweepD,self.detect_over)
        self.log.debug("initial AP detection: %d APs"%len(Is))

        # eliminate APs where dV/dT doesn't cross below -10 V/S within 2 ms
        for i,I in enumerate(Is):
            if np.min(self.abf.sweepD[I:I+2*self.abf.pointsPerMs])>-10:
                Is[i]=0
        Is=Is[np.nonzero(Is)]
        self.log.debug("after lower threshold checking: %d APs"%len(Is))

        # walk 1ms backwards and find point of +10 V/S threshold crossing
        for i,I in enumerate(Is):
            stepBack=0
            while(self.abf.sweepD[I-stepBack])>10 and stepBack/self.abf.pointsPerMs<1: #2ms max
                stepBack+=1
            Is[i]-=stepBack

        # analyze each AP
        sweepAPs=[]
        for i,I in enumerate(Is):
            try:
                timeInSweep=I/self.abf.pointsPerSec
                if timeInSweep<self.detect_time1 or timeInSweep>self.detect_time2:
                    continue # skip because it's not within the marks
                ap={} # create the AP entry
                ap["sweep"]=sweep # number of the sweep containing this AP
                ap["I"]=I # index sweep point of start of AP (10 mV/ms threshold crossing)
                ap["Tsweep"]=I/self.abf.pointsPerSec # time in the sweep of index crossing (sec)
                ap["T"]=ap["Tsweep"]+self.abf.sweepInterval*sweep # time in the experiment
                ap["Vthreshold"]=self.abf.sweepY[I] # threshold at rate of -10mV/ms

                # determine how many points from the start dV/dt goes below -10 (from a 5ms chunk)
                chunk=self.abf.sweepD[I:I+5*self.abf.pointsPerMs] # give it 5ms to cross once
                I_toNegTen=np.where(chunk<-10)[0][0]
                chunk=self.abf.sweepD[I+I_toNegTen:I+I_toNegTen+10*self.abf.pointsPerMs] # give it 30ms to cross back
                if not max(chunk)>-10:
                    self.log.debug("skipping unreal AP at T=%f"%ap["T"])
                    self.log.error("^^^ can you confirm this is legit?")
                    continue # probably a pre-AP "bump" to be ignored
                I_recover=np.where(chunk>-10)[0][0]+I_toNegTen+I # point where trace returns to above -10 V/S
                ap["dVfastIs"]=[I,I_recover] # span of the fast component of the dV/dt trace
                ap["dVfastMS"]=(I_recover-I)/self.abf.pointsPerMs # time (in ms) of this fast AP component

                # determine derivative min/max from a 2ms chunk which we expect to capture the fast AP
                chunk=self.abf.sweepD[ap["dVfastIs"][0]:ap["dVfastIs"][1]]
                ap["dVmax"]=np.max(chunk)
                ap["dVmaxI"]=np.where(chunk==ap["dVmax"])[0][0]+I
                ap["dVmin"]=np.min(chunk)
                ap["dVminI"]=np.where(chunk==ap["dVmin"])[0][0]+I
                if ap["dVmax"]<10 or ap["dVmin"]>-10:
                    self.log.debug("throwing out AP with low dV/dt to be an AP")
                    self.log.error("^^^ can you confirm this is legit?")
                    continue

                # before determining AP shape stats, see where trace recovers to threshold
                chunkSize=self.abf.pointsPerMs*10 #AP shape may be 10ms
                if len(Is)-1>i and Is[i+1]<(I+chunkSize): # if slow AP runs into next AP
                    chunkSize=Is[i+1]-I # chop it down
                if chunkSize<(self.abf.pointsPerMs*2):
                    continue # next AP is so soon, it's >500 Hz. Can't be real.
                ap["VslowIs"]=[I,I+chunkSize] # time range of slow AP dynamics
                chunk=self.abf.sweepY[I:I+chunkSize]

                # determine AP peak and minimum
                ap["Vmax"]=np.max(chunk)
                ap["VmaxI"]=np.where(chunk==ap["Vmax"])[0][0]+I
                chunkForMin=np.copy(chunk) # so we can destroy it
                chunkForMin[:ap["VmaxI"]-I]=np.inf # minimum won't be before peak now
                ap["Vmin"]=np.min(chunkForMin) # supposedly the minimum is the AHP
                ap["VminI"]=np.where(chunkForMin==ap["Vmin"])[0][0]+I
                if ap["VminI"]<ap["VmaxI"]:
                    self.log.error("-------------------------------")
                    self.log.error("how is the AHP before the peak?") #TODO: start chunk at the peak
                    self.log.error("-------------------------------")
                #print((I+len(chunk))-ap["VminI"],len(chunk))
                if (len(chunk))-((I+len(chunk))-ap["VminI"])<10:
                    self.log.error("-------------------------------")
                    self.log.error("HP too close for comfort!")
                    self.log.error("-------------------------------")

                ap["msRiseTime"]=(ap["VmaxI"]-I)/self.abf.pointsPerMs # time from threshold to peak
                ap["msFallTime"]=(ap["VminI"]-ap["VmaxI"])/self.abf.pointsPerMs # time from peak to nadir

                # determine halfwidth
                ap["Vhalf"]=np.average([ap["Vmax"],ap["Vthreshold"]]) # half way from threshold to peak
                ap["VhalfI1"]=cm.where_cross(chunk,ap["Vhalf"])[0]+I # time it's first crossed
                ap["VhalfI2"]=cm.where_cross(-chunk,-ap["Vhalf"])[1]+I # time it's second crossed
                ap["msHalfwidth"]=(ap["VhalfI2"]-ap["VhalfI1"])/self.abf.pointsPerMs # time between crossings

                # AP error checking goes here
                # TODO:

                # if we got this far, add the AP to the list
                sweepAPs.extend([ap])
            except Exception as e:
                self.log.error("crashed analyzing AP %d of %d",i,len(Is))
                self.log.error(cm.exceptionToString(e))
                #cm.pause()
                #cm.waitFor(30)
                #self.log.error("EXCEPTION!:\n%s"%str(sys.exc_info()))

        self.log.debug("finished analyzing sweep. Found %d APs",len(sweepAPs))
        self.APs.extend(sweepAPs)
        self.abf.derivative=False
def detectSweep(self,sweep=0):
        """perform AP detection on current sweep."""

        if self.APs is False: # indicates detection never happened
            self.APs=[] # now indicates detection occured

        # delete every AP from this sweep from the existing array
        for i,ap in enumerate(self.APs):
            if ap["sweep"]==sweep:
                self.APs[i]=None
        if self.APs.count(None):
            self.log.debug("deleting %d existing APs from memory",self.APs.count(None))
            while None in self.APs:
                self.APs.remove(None)
        self.log.debug("initiating AP detection (%d already in memory)",len(self.APs))

        self.abf.derivative=True
        self.abf.setsweep(sweep)

        # detect potential AP (Is) by a dV/dT threshold crossing
        Is = cm.where_cross(self.abf.sweepD,self.detect_over)
        self.log.debug("initial AP detection: %d APs"%len(Is))

        # eliminate APs where dV/dT doesn't cross below -10 V/S within 2 ms
        for i,I in enumerate(Is):
            if np.min(self.abf.sweepD[I:I+2*self.abf.pointsPerMs])>-10:
                Is[i]=0
        Is=Is[np.nonzero(Is)]
        self.log.debug("after lower threshold checking: %d APs"%len(Is))

        # walk 1ms backwards and find point of +10 V/S threshold crossing
        for i,I in enumerate(Is):
            stepBack=0
            while(self.abf.sweepD[I-stepBack])>10 and stepBack/self.abf.pointsPerMs<1: #2ms max
                stepBack+=1
            Is[i]-=stepBack

        # analyze each AP
        sweepAPs=[]
        for i,I in enumerate(Is):
            try:
                timeInSweep=I/self.abf.pointsPerSec
                if timeInSweep<self.detect_time1 or timeInSweep>self.detect_time2:
                    continue # skip because it's not within the marks
                ap={} # create the AP entry
                ap["sweep"]=sweep # number of the sweep containing this AP
                ap["I"]=I # index sweep point of start of AP (10 mV/ms threshold crossing)
                ap["Tsweep"]=I/self.abf.pointsPerSec # time in the sweep of index crossing (sec)
                ap["T"]=ap["Tsweep"]+self.abf.sweepInterval*sweep # time in the experiment
                ap["Vthreshold"]=self.abf.sweepY[I] # threshold at rate of -10mV/ms

                # determine how many points from the start dV/dt goes below -10 (from a 5ms chunk)
                chunk=self.abf.sweepD[I:I+5*self.abf.pointsPerMs] # give it 5ms to cross once
                I_toNegTen=np.where(chunk<-10)[0][0]
                chunk=self.abf.sweepD[I+I_toNegTen:I+I_toNegTen+10*self.abf.pointsPerMs] # give it 30ms to cross back
                if not max(chunk)>-10:
                    self.log.debug("skipping unreal AP at T=%f"%ap["T"])
                    self.log.error("^^^ can you confirm this is legit?")
                    continue # probably a pre-AP "bump" to be ignored
                I_recover=np.where(chunk>-10)[0][0]+I_toNegTen+I # point where trace returns to above -10 V/S
                ap["dVfastIs"]=[I,I_recover] # span of the fast component of the dV/dt trace
                ap["dVfastMS"]=(I_recover-I)/self.abf.pointsPerMs # time (in ms) of this fast AP component

                # determine derivative min/max from a 2ms chunk which we expect to capture the fast AP
                chunk=self.abf.sweepD[ap["dVfastIs"][0]:ap["dVfastIs"][1]]
                ap["dVmax"]=np.max(chunk)
                ap["dVmaxI"]=np.where(chunk==ap["dVmax"])[0][0]+I
                ap["dVmin"]=np.min(chunk)
                ap["dVminI"]=np.where(chunk==ap["dVmin"])[0][0]+I
                if ap["dVmax"]<10 or ap["dVmin"]>-10:
                    self.log.debug("throwing out AP with low dV/dt to be an AP")
                    self.log.error("^^^ can you confirm this is legit?")
                    continue

                # before determining AP shape stats, see where trace recovers to threshold
                chunkSize=self.abf.pointsPerMs*10 #AP shape may be 10ms
                if len(Is)-1>i and Is[i+1]<(I+chunkSize): # if slow AP runs into next AP
                    chunkSize=Is[i+1]-I # chop it down
                if chunkSize<(self.abf.pointsPerMs*2):
                    continue # next AP is so soon, it's >500 Hz. Can't be real.
                ap["VslowIs"]=[I,I+chunkSize] # time range of slow AP dynamics
                chunk=self.abf.sweepY[I:I+chunkSize]

                # determine AP peak and minimum
                ap["Vmax"]=np.max(chunk)
                ap["VmaxI"]=np.where(chunk==ap["Vmax"])[0][0]+I
                chunkForMin=np.copy(chunk) # so we can destroy it
                chunkForMin[:ap["VmaxI"]-I]=np.inf # minimum won't be before peak now
                ap["Vmin"]=np.min(chunkForMin) # supposedly the minimum is the AHP
                ap["VminI"]=np.where(chunkForMin==ap["Vmin"])[0][0]+I
                if ap["VminI"]<ap["VmaxI"]:
                    self.log.error("-------------------------------")
                    self.log.error("how is the AHP before the peak?") #TODO: start chunk at the peak
                    self.log.error("-------------------------------")
                #print((I+len(chunk))-ap["VminI"],len(chunk))
                if (len(chunk))-((I+len(chunk))-ap["VminI"])<10:
                    self.log.error("-------------------------------")
                    self.log.error("HP too close for comfort!")
                    self.log.error("-------------------------------")

                ap["msRiseTime"]=(ap["VmaxI"]-I)/self.abf.pointsPerMs # time from threshold to peak
                ap["msFallTime"]=(ap["VminI"]-ap["VmaxI"])/self.abf.pointsPerMs # time from peak to nadir

                # determine halfwidth
                ap["Vhalf"]=np.average([ap["Vmax"],ap["Vthreshold"]]) # half way from threshold to peak
                ap["VhalfI1"]=cm.where_cross(chunk,ap["Vhalf"])[0]+I # time it's first crossed
                ap["VhalfI2"]=cm.where_cross(-chunk,-ap["Vhalf"])[1]+I # time it's second crossed
                ap["msHalfwidth"]=(ap["VhalfI2"]-ap["VhalfI1"])/self.abf.pointsPerMs # time between crossings

                # AP error checking goes here
                # TODO:

                # if we got this far, add the AP to the list
                sweepAPs.extend([ap])
            except Exception as e:
                self.log.error("crashed analyzing AP %d of %d",i,len(Is))
                self.log.error(cm.exceptionToString(e))
                #cm.pause()
                #cm.waitFor(30)
                #self.log.error("EXCEPTION!:\n%s"%str(sys.exc_info()))

        self.log.debug("finished analyzing sweep. Found %d APs",len(sweepAPs))
        self.APs.extend(sweepAPs)
        self.abf.derivative=False
def detectSweep ( self , sweep = 0 ) : if self . APs is False : # indicates detection never happened self . APs = [ ] # now indicates detection occured # delete every AP from this sweep from the existing array for i , ap in enumerate ( self . APs ) : if ap [ "sweep" ] == sweep : self . APs [ i ] = None if self . APs . count ( None ) : self . log . debug ( "deleting %d existing APs from memory" , self . APs . count ( None ) ) while None in self . APs : self . APs . remove ( None ) self . log . debug ( "initiating AP detection (%d already in memory)" , len ( self . APs ) ) self . abf . derivative = True self . abf . setsweep ( sweep ) # detect potential AP (Is) by a dV/dT threshold crossing Is = cm . where_cross ( self . abf . sweepD , self . detect_over ) self . log . debug ( "initial AP detection: %d APs" % len ( Is ) ) # eliminate APs where dV/dT doesn't cross below -10 V/S within 2 ms for i , I in enumerate ( Is ) : if np . min ( self . abf . sweepD [ I : I + 2 * self . abf . pointsPerMs ] ) > - 10 : Is [ i ] = 0 Is = Is [ np . nonzero ( Is ) ] self . log . debug ( "after lower threshold checking: %d APs" % len ( Is ) ) # walk 1ms backwards and find point of +10 V/S threshold crossing for i , I in enumerate ( Is ) : stepBack = 0 while ( self . abf . sweepD [ I - stepBack ] ) > 10 and stepBack / self . abf . pointsPerMs < 1 : #2ms max stepBack += 1 Is [ i ] -= stepBack # analyze each AP sweepAPs = [ ] for i , I in enumerate ( Is ) : try : timeInSweep = I / self . abf . pointsPerSec if timeInSweep < self . detect_time1 or timeInSweep > self . detect_time2 : continue # skip because it's not within the marks ap = { } # create the AP entry ap [ "sweep" ] = sweep # number of the sweep containing this AP ap [ "I" ] = I # index sweep point of start of AP (10 mV/ms threshold crossing) ap [ "Tsweep" ] = I / self . abf . pointsPerSec # time in the sweep of index crossing (sec) ap [ "T" ] = ap [ "Tsweep" ] + self . abf . sweepInterval * sweep # time in the experiment ap [ "Vthreshold" ] = self . abf . sweepY [ I ] # threshold at rate of -10mV/ms # determine how many points from the start dV/dt goes below -10 (from a 5ms chunk) chunk = self . abf . sweepD [ I : I + 5 * self . abf . pointsPerMs ] # give it 5ms to cross once I_toNegTen = np . where ( chunk < - 10 ) [ 0 ] [ 0 ] chunk = self . abf . sweepD [ I + I_toNegTen : I + I_toNegTen + 10 * self . abf . pointsPerMs ] # give it 30ms to cross back if not max ( chunk ) > - 10 : self . log . debug ( "skipping unreal AP at T=%f" % ap [ "T" ] ) self . log . error ( "^^^ can you confirm this is legit?" ) continue # probably a pre-AP "bump" to be ignored I_recover = np . where ( chunk > - 10 ) [ 0 ] [ 0 ] + I_toNegTen + I # point where trace returns to above -10 V/S ap [ "dVfastIs" ] = [ I , I_recover ] # span of the fast component of the dV/dt trace ap [ "dVfastMS" ] = ( I_recover - I ) / self . abf . pointsPerMs # time (in ms) of this fast AP component # determine derivative min/max from a 2ms chunk which we expect to capture the fast AP chunk = self . abf . sweepD [ ap [ "dVfastIs" ] [ 0 ] : ap [ "dVfastIs" ] [ 1 ] ] ap [ "dVmax" ] = np . max ( chunk ) ap [ "dVmaxI" ] = np . where ( chunk == ap [ "dVmax" ] ) [ 0 ] [ 0 ] + I ap [ "dVmin" ] = np . min ( chunk ) ap [ "dVminI" ] = np . where ( chunk == ap [ "dVmin" ] ) [ 0 ] [ 0 ] + I if ap [ "dVmax" ] < 10 or ap [ "dVmin" ] > - 10 : self . log . debug ( "throwing out AP with low dV/dt to be an AP" ) self . log . error ( "^^^ can you confirm this is legit?" ) continue # before determining AP shape stats, see where trace recovers to threshold chunkSize = self . abf . pointsPerMs * 10 #AP shape may be 10ms if len ( Is ) - 1 > i and Is [ i + 1 ] < ( I + chunkSize ) : # if slow AP runs into next AP chunkSize = Is [ i + 1 ] - I # chop it down if chunkSize < ( self . abf . pointsPerMs * 2 ) : continue # next AP is so soon, it's >500 Hz. Can't be real. ap [ "VslowIs" ] = [ I , I + chunkSize ] # time range of slow AP dynamics chunk = self . abf . sweepY [ I : I + chunkSize ] # determine AP peak and minimum ap [ "Vmax" ] = np . max ( chunk ) ap [ "VmaxI" ] = np . where ( chunk == ap [ "Vmax" ] ) [ 0 ] [ 0 ] + I chunkForMin = np . copy ( chunk ) # so we can destroy it chunkForMin [ : ap [ "VmaxI" ] - I ] = np . inf # minimum won't be before peak now ap [ "Vmin" ] = np . min ( chunkForMin ) # supposedly the minimum is the AHP ap [ "VminI" ] = np . where ( chunkForMin == ap [ "Vmin" ] ) [ 0 ] [ 0 ] + I if ap [ "VminI" ] < ap [ "VmaxI" ] : self . log . error ( "-------------------------------" ) self . log . error ( "how is the AHP before the peak?" ) #TODO: start chunk at the peak self . log . error ( "-------------------------------" ) #print((I+len(chunk))-ap["VminI"],len(chunk)) if ( len ( chunk ) ) - ( ( I + len ( chunk ) ) - ap [ "VminI" ] ) < 10 : self . log . error ( "-------------------------------" ) self . log . error ( "HP too close for comfort!" ) self . log . error ( "-------------------------------" ) ap [ "msRiseTime" ] = ( ap [ "VmaxI" ] - I ) / self . abf . pointsPerMs # time from threshold to peak ap [ "msFallTime" ] = ( ap [ "VminI" ] - ap [ "VmaxI" ] ) / self . abf . pointsPerMs # time from peak to nadir # determine halfwidth ap [ "Vhalf" ] = np . average ( [ ap [ "Vmax" ] , ap [ "Vthreshold" ] ] ) # half way from threshold to peak ap [ "VhalfI1" ] = cm . where_cross ( chunk , ap [ "Vhalf" ] ) [ 0 ] + I # time it's first crossed ap [ "VhalfI2" ] = cm . where_cross ( - chunk , - ap [ "Vhalf" ] ) [ 1 ] + I # time it's second crossed ap [ "msHalfwidth" ] = ( ap [ "VhalfI2" ] - ap [ "VhalfI1" ] ) / self . abf . pointsPerMs # time between crossings # AP error checking goes here # TODO: # if we got this far, add the AP to the list sweepAPs . extend ( [ ap ] ) except Exception as e : self . log . error ( "crashed analyzing AP %d of %d" , i , len ( Is ) ) self . log . error ( cm . exceptionToString ( e ) ) #cm.pause() #cm.waitFor(30) #self.log.error("EXCEPTION!:\n%s"%str(sys.exc_info())) self . log . debug ( "finished analyzing sweep. Found %d APs" , len ( sweepAPs ) ) self . APs . extend ( sweepAPs ) self . abf . derivative = False
perform AP detection on current sweep.
perform AP detection on current sweep .