Python pylab.find函数代码示例

def classify_epoch(epoch,rate):
    """
    This function returns a sleep stage classification (integers: 1 for NREM
    stage 1, 2 for NREM stage 2, and 3 for NREM stage 3/4) given an epoch of 
    EEG and a sampling rate.
    """
    Pxx, freqs = m.psd(epoch,NFFT=256,Fs=rate)
    nPxx = Pxx/float(sum(Pxx))
    delta_f = plt.find((0<freqs) & (freqs <=3))
 #   delta_power = sum(Pxx[delta_f])
    ndelta_power = sum(nPxx[delta_f])
    spindles_f = plt.find((11 <= freqs) & (freqs <= 15))
 #   spindle_power = sum(Pxx[spindles_f])
    nspindle_power = sum(nPxx[spindles_f])
  #  ratio = spindle_power/delta_power
    
    if (ndelta_power > 0.8): #suggests stage 3
        stage = 3
    else:  
        if (nspindle_power > 0.03): #suggests stage 2
            stage = 2
        else: 
            stage = 1
    
    return stage

def compute_firing_rate_per_motion( direction_of_motions, spk_times, trials, time_bin ):
    
   avg_firing_rate = np.zeros( len( direction_of_motions ) )
   
   for d in range( len( direction_of_motions ) ):
       
       # extract the trial times for each direction
       time_indx = plt.find( direction_of_motions[d] == trials[:,0] ) 
       times_motion_started = trials[ time_indx, 1 ]
       
       total_spikes_per_trial = 0
       
       # for each time of the trial where the movement started,
       # count the number of spikes that occurred in the time window
       for t in range( len(times_motion_started) ):
           spike_indx = plt.find( (spk_times >= times_motion_started[t] - time_bin) & (spk_times <= times_motion_started[t] + time_bin)  )
           total_spikes_per_trial += len( spike_indx )
       
       # compute the average firing rate
       avg_rate = ( (1.0*total_spikes_per_trial)  / (1.0*len(times_motion_started))) / time_bin
       avg_rate /= 2
       
       # append the firing rate to 
       avg_firing_rate[d] = avg_rate
          
   #append the direction d and the firing rate to the initially empty array
   dir_rates = np.column_stack( (direction_of_motions, avg_firing_rate)  )
   
   return dir_rates

def get_features(data_epoch,rate):
    pxxs = []
    for eeg_epoch in data_epoch:
        pxx, freqs = m.psd(eeg_epoch, Fs=rate)
    pxxs.append(sum(pxx[plt.find((10 <= freqs) & (freqs <= 15))]))
    pxxs.append(max(pxx[plt.find((10 <= freqs) & (freqs <= 15))]))
    pxxs.append(np.mean(pxx[plt.find((10 <= freqs) & (freqs <= 15))]))
    return pxxs    

def detector_tester(APTimes, actualTimes):
    """
    returns percentTrueSpikes (% correct detected) and falseSpikeRate
    (extra APs per second of data)
    compares actual spikes times with detected spike times
    This only works if we give you the answers!
    """
    
    JITTER = 0.0025 #2 ms of jitter allowed
    
    #first match the two sets of spike times. Anything within JITTER_MS
    #is considered a match (but only one per time frame!)
    
    #order the lists
    detected = np.sort(APTimes)
    actual = np.sort(actualTimes)
    
    #remove spikes with the same times (these are false APs)
    temp = np.append(detected, -1)
    detected = detected[plt.find(plt.diff(temp) != 0)]
 
    #find matching action potentials and mark as matched (trueDetects)
    trueDetects = [];
    for sp in actual:
        z = plt.find((detected >= sp-JITTER) & (detected <= sp+JITTER))
        if len(z)>0:
            for i in z:
                zz = plt.find(trueDetects == detected[i])
                if len(zz) == 0:
                    trueDetects = np.append(trueDetects, detected[i])
                    break;
    percentTrueSpikes = 100.0*len(trueDetects)/len(actualTimes)
    
    
    #everything else is a false alarm
    totalTime = (actual[len(actual)-1]-actual[0])
    falseSpikeRate = (len(APTimes) - len(actualTimes))/totalTime
    
    # Added this for auto-evaluation based on criteria 
    pct_spike_eval = "PASS" if percentTrueSpikes > 90.0 else "FAIL"
    false_spike_eval = "PASS" if falseSpikeRate < 2.5 else "FAIL"
    overall_result = "FAIL" if pct_spike_eval == "FAIL" or false_spike_eval == "FAIL" else "PASS"    
    
    print 'Action Potential Detector Performance performance: '
    print '     Correct number of action potentials = %d' % len(actualTimes)
    print '     %s: Percent True Spikes = %f' % (pct_spike_eval, percentTrueSpikes)
    print '     %s: False Spike Rate = %f spikes/s' % (false_spike_eval, falseSpikeRate)
    print ''
    print 'Overall Evaluation: %s' % overall_result
    print ''
    return {'Percent True Spikes':percentTrueSpikes, 'False Spike Rate':falseSpikeRate}

def scoresStats(scoresList):
    scoresList = array(scoresList)
    nbPlayers = scoresList.shape[1]
    print "Player | wins loses (ties) meanWinMargin meanLossMargin"
    for player in xrange(nbPlayers):
        otherPlayer = (player+1) % nbPlayers
        winning = find(scoresList[:,player] > scoresList[:,otherPlayer])
        losing  = find(scoresList[:,player] < scoresList[:,otherPlayer])
        wins,loses = len(winning),len(losing)
        ties = scoresList.shape[0]-wins-loses
        meanWinningMargin = mean(sum(scoresList[winning],1))
        meanLosingMargin  = mean(sum(scoresList[losing],1))
        print "     %c | %d %d (%d) +%f -%f" % \
            (chr(65+player), wins, loses, ties, \
             meanWinningMargin, meanLosingMargin)

def plot_waveforms(time,voltage,APTimes,titlestr):
    """
    plot_waveforms takes four arguments - the recording time array, the voltage
    array, the time of the detected action potentials, and the title of your
    plot.  The function creates a labeled plot showing the waveforms for each
    detected action potential
    """

    plt.figure()

    ## Your Code Here
    for n in range(0,APTimes.size):
        ind = time[(time>(APTimes[n] - 0.003)) & (time<(APTimes[n] + 0.003))]
        sp = np.zeros(ind.size)
    
        for i in range(0,ind.size):
            sp[i] = plt.find(time == ind[i])
        sp=sp.astype(np.int64)
        x = np.linspace(-3.0e-3, 3.0e-3, sp.size)
        plt.plot(x,voltage[sp])
        plt.hold(True)
        
    plt.xlabel('Time (s)')
    plt.ylabel('Voltages (s)')
    plt.title(titlestr)    
    plt.show()

def analys_trans(fname, maxSr=15, srHstBins=10):
    srHstBinSize = maxSr/srHstBins

    (time, sr) = np.load(fname)
#    sr = gs_filter(sr, 5)

    srHst = np.histogram(sr, bins=srHstBins, range=(0, maxSr))[0]

    srMidpntIdx = len(srHst)/2
    firstMaxIdx = np.argmax(srHst[:srMidpntIdx])
    secondMaxIdx = srMidpntIdx + np.argmax(srHst[srMidpntIdx:])
    srMidpnt = np.argmin(srHst[firstMaxIdx:secondMaxIdx])*srHstBinSize
#    srMidpnt = (firstMaxIdx + secondMaxIdx)/2
    print(srMidpnt)
    srMidpnt = 50.

    thr = np.array(sr >= srMidpnt, dtype='int32')
    df = np.diff(thr)
    indices_up = pl.find(df == 1)
    indices_down = pl.find(df == -1)

    TimesUp = []
    ActMeanUp = []
    TimesDown = []
    ActMeanDown = []
    Periods = []
    # all time in sec
    if sr[0] > srMidpnt:
        for up, down in zip(indices_up, indices_down):
            TimesDown.append(up - down)
            ActMeanDown.append(np.mean(sr[down:up]))
        for up, down in zip(indices_up, indices_down[1:]):
            TimesUp.append(down - up)
            ActMeanUp.append(np.mean(sr[up:down]))
        Periods = np.diff(indices_down)
    elif sr[0] <= srMidpnt:
        for up, down in zip(indices_up, indices_down):
            TimesUp.append(down - up)
            ActMeanUp.append(np.mean(sr[up:down]))
        for up, down in zip(indices_up[1:], indices_down):
            TimesDown.append(up - down)
            ActMeanDown.append(np.mean(sr[down:up]))
        Periods = np.diff(indices_up)
    else:
        raise Exception("Unexpected type of sr")
    return np.array(Periods), np.array(TimesDown), np.array(TimesUp), np.array(ActMeanDown), np.array(ActMeanUp)

def detector_tester(APTimes, actualTimes):
    """
    returns percentTrueSpikes (% correct detected) and falseSpikeRate
    (extra APs per second of data)
    compares actual spikes times with detected spike times
    This only works if we give you the answers!
    """
    
    #print APTimes
    #print actualTimes
    JITTER = 0.0025 #2 ms of jitter allowed
    
    #first match the two sets of spike times. Anything within JITTER_MS
    #is considered a match (but only one per time frame!)
    
    #order the lists
    detected = np.sort(APTimes)
    actual = np.sort(actualTimes)
    
    #remove spikes with the same times (these are false APs)
    temp = np.append(detected, -1)
    detected = detected[plt.find(plt.diff(temp) != 0)]
 
    #find matching action potentials and mark as matched (trueDetects)
    trueDetects = [];
    for sp in actual:
        z = plt.find((detected >= sp-JITTER) & (detected <= sp+JITTER))
        if len(z)>0:
            for i in z:
                zz = plt.find(trueDetects == detected[i])
                if len(zz) == 0:
                    trueDetects = np.append(trueDetects, detected[i])
                    break;
    percentTrueSpikes = 100.0*len(trueDetects)/len(actualTimes)
    
    #everything else is a false alarm
    totalTime = (actual[len(actual)-1]-actual[0])
    falseSpikeRate = (len(APTimes) - len(actualTimes))/totalTime
   
    
    print 'Action Potential Detector Performance performance: '
    print '     Correct number of action potentials = ' + str(len(actualTimes))
    print '     Percent True Spikes = ' + str(percentTrueSpikes)
    print '     False Spike Rate = ' + str(falseSpikeRate) + ' spikes/s'
    print 
    return {'Percent True Spikes':percentTrueSpikes, 'False Spike Rate':falseSpikeRate}

def time_spent_in_stages(stages):
    i=0
    timespent_array=[]
    for i in range(8):
        indices=plt.find(stages==i)
        j=len(indices)*30
        timespent_array.append(j)
    return timespent_array

def plot_waveforms(time,voltage,APTimes,titlestr):
    """
    plot_waveforms takes four arguments - the recording time array, the voltage
    array, the time of the detected action potentials, and the title of your
    plot.  The function creates a labeled plot showing the waveforms for each
    detected action potential
    """
    # calculate sample rate (frames per ms)

    sample_rate = float(len(time)-1)/(1000*(time[-1]-time[0]))
    plt.figure()
    # cycle for each detected spike
    
    for i in range(0, len(APTime)): 
        # find the index for the i-th detected spike time in the time, 
        # using plt.find
        idx_spike = plt.find(time == APTime[i])
        #please note that the use of "==" as an condition could be dangerous,
        #since it requires that the the detectedSpikeTimesMS[i] to be exact the same with one element in the timesMS, """
        #to avoild this, we can instead use  
        idx_spike = min(plt.find(abs(time - APTime[i])<0.00001))
        #to allow 0.01ms jittering
        #I use min() here in case it returns several indices
               
        # find the start and end index for plotting
        idx_start = idx_spike - int(3*sample_rate)
        idx_end   = idx_spike + int(3*sample_rate)

     
        #remeber that sample_rate may not be an interger, 
        #so we have to use int() to convert it into an interger
        
        #note if the index could potentially be negative or bigger
        #than the length of the array at this point, an error will occur
        #we have to make sure that this does not happen 
        if (idx_start>=0) & (idx_end < len(time)-1):
            # plot the waveform for the i-th spike
            xx = np.linspace(-3,3,sample_rate*6)
            yy = voltage[range(idx_start, idx_end)]
            plt.plot(xx, yy, 'b',hold=True,)
    
    # add axis labels and title
    plt.xlabel('Time (ms)')
    plt.ylabel('Voltage (uV)')  
    plt.title(titlestr)
    plt.show()

def subject34_stage4_list():
    """
    list of indices for stage 4s
    """
    stage4_indices = [[[], []], [[], []]]
    stages, stagelengths = pull_subject34_stages()
    for sidx in range(2):
        for bridx in range(2):
            stage4_indices[sidx][bridx] = plt.find(stages[sidx, bridx] == 4)
    gc.collect()
    return np.array(stage4_indices), stagelengths

def preferred_direction(fit_curve):
    """
    The function takes a 2-dimensional array with the x-values of the fit curve
    in the first column and the y-values of the fit curve in the second.  
    It returns the preferred direction of the neuron (in degrees).
    """
    
    # find the motion that is associated with the highest firing rate
    prefered_value = max( fit_curve[:,1] )
    preferred_direction = plt.find(  fit_curve[:,1] == prefered_value )
    
    return fit_curve[preferred_direction[0],0]

def good_AP_finder2(time,voltage):
    """
    This function takes the following input:
        time - vector where each element is a time in seconds
        voltage - vector where each element is a voltage at a different time
        
        We are assuming thaht he two vectors are in correspondance (meaning
        that at a given index, the time in one corresponds to the voltage in
        the other). The vectors must be the same size or the code
        won't run
    
    This function returns the following output:
        APTime - all the times where a spike (action potential) was detected
    """
    
    #Let's make sure the input looks at least reasonable
    if (len(voltage) != len(time)):
        print "Can't run - the vectors aren't the same length!"
        APTime = []
        return APTime
     
    
    # Pick a threshold. You can eyeball it by looking at the plot, or you can 
    # write code to find it.  Code would be better, but isn't 100% necessary.
    thrd = -80
    
    # find all the indices whose corresponding voltage is lower than the threshold
    detectedAPIndex = plt.find(voltage < thrd)
    # note that now several neighboring indices could correspond to the same spike
    
    
    # we only want the first index for one spike
    # so we will throw away several frames following the first one
    
    # calculate difference of the picked neiboring indices
    diff_detectedAPIndex= plt.diff(detectedAPIndex)
    
    # if diff_detectedAPIndex>1, we know that it's a new spike
    # note that diff omits the first element, which is a spike, so I insert the first one
    detectedAPIndex_select = np.insert(diff_detectedAPIndex>1, 0, True)
    # detectedAPIndex_select is a boolean array with the same length of detectedAPIndex
    
    # we selecte the indices that correspond to the begginning frame of each spikes
    detectedAPIndex = detectedAPIndex[detectedAPIndex_select]
    
    # find the time im ms based on the indices
    APTime =list(time[i] for i in detectedAPIndex)
    
    return APTime

def plotIvsR(ix, iy, fx, fy, cx, cy):
    RF = sqrt((fx - cx)**2 + (fy - cy)**2)
    mD2 = []
    mI2 = []
    for i in range(len(RF)):
        r = 5
        D = sqrt((ix - fx[i])**2 + (iy - fy[i])**2)
        I = find( D < r )
        for j in I:
            mD2.append(D[j])
            mI2.append(i)
    #plot(mI, mD, 'ro', mI2, mD2, 'bo')
    plot(mI2, mD2, 'bo')
    xlabel('Index star number')
    ylabel('Match distance')

def bin_spikes(trials, spk_times, time_bin = 0.08):
    """
    bin_spikes takes the trials array (with directions and times) and the spk_times
    array with spike times and returns the average firing rate for each of the
    eight directions of motion, as calculated within a time_bin before and after
    the trial time (time_bin should be given in seconds).  For example,
    time_bin = .1 will count the spikes from 100ms before to 100ms after the 
    trial began.
    
    dir_rates should be an 8x2 array with the first column containing the directions
    (in degrees from 0-360) and the second column containing the average firing rate
    for each direction
    """
    directions = np.unique(trials[:,0])
    dir_rates = np.zeros( (8, 2) ) # intialize an 8x2 array of zeros 
    for direction in directions:
        nmbr_trials = float(len(plt.find(trials[:,0]==direction)))
        trial_times = [ t[1] for t in trials if t[0] == direction ]
        nmbr_spks = 0.0
        for time in trial_times:
            nmbr_spks += len(np.where(np.logical_and(spk_times>time-time_bin,spk_times<time+time_bin))[0])
        firing_rate = nmbr_spks/nmbr_trials/(2*time_bin)
        dir_rates[plt.find(directions==direction)] = [direction, firing_rate] 
    return dir_rates