def streaks(rats):
    
    """ Okay, I want to find consecutive hits during uncued blocks 
        
        So, what I'll do first is grab the 
        """
    
    d_list = []
    for rat in rats.itervalues():
        trials = rat.sessions['trials']
        uncued = [ find(trial['block'] == 'uncued') for trial in trials ]
        cued = [ find(trial['block'] == 'cued') for trial in trials ]
        
        un_blk_lens = [ np.diff(session)-1 for session in cued ]
        un_blk_lens = np.array([ sess[find(sess != 0)] for sess in uncued_lens ])
        
        un_hits = [ trial['hits'][uncued[ii]] for ii, trial in enumerate(trials) ]
        un_hits = [ find(hts == 1) for hts in un_hits ]
        
        # Grab the first block of uncued trials
        un_blk_lens
        
        
        d_list.append(diffs)
    
    return d_list
    

def getdrift_raw(filename,id3,interval,datetime_wanted):
    
  # range_time is a number,unit by one day.  datetime_wanted format is num
  d=ml.load(filename)
  lat1=d[:,8]
  lon1=d[:,7]
  idd=d[:,0]
  year=[]
  for n in range(len(idd)):
      year.append(str(idd[n])[0:2])
  h=d[:,4]
  day=d[:,3]
  month=d[:,2]
  time1=[]
  for i in range(len(idd)):
      time1.append(date2num(datetime.datetime.strptime(str(int(h[i]))+' '+str(int(day[i]))+' '+str(int(month[i]))+' '+str(int(year[i])), "%H %d %m %y")))


  idg1=list(ml.find(idd==id3))
  idg2=list(ml.find(np.array(time1)<=datetime_wanted+interval/24))
  "'0.25' means the usual Interval, It can be changed base on different drift data "
  idg3=list(ml.find(np.array(time1)>=datetime_wanted-0.1))
  idg23=list(set(idg2).intersection(set(idg3)))
  # find which data we need
  idg=list(set(idg23).intersection(set(idg1)))
  print 'the length of drifter data is  '+str(len(idg)),str(len(set(idg)))+'   . if same, no duplicate'
  lat,lon,time=[],[],[]
  
  for x in range(len(idg)):
      lat.append(round(lat1[idg[x]],4))
      lon.append(round(lon1[idg[x]],4))
      time.append(round(time1[idg[x]],4))
  # time is num
  return lat,lon,time

    def Extract_Range(self, begin, end):
        #Extract a date range and use it to make a new time series.
        
        try: 
            datetime.datetime.now() > begin
        except TypeError:
            begin = parser.parse(begin)
            
        try: 
            datetime.datetime.now() > end
        except TypeError:
            end = parser.parse(end)
        
        #Find which rows of timestamps lie in this range.
        upper_rows = mlab.find( self.timestamps >= begin )
        lower_rows = mlab.find( self.timestamps <= end )
        rows = filter( lambda x: x in upper_rows, lower_rows )

        #Now create a time series with the data in this range.
        range = {}
        range['headers'] = self.headers
        range['timestamps'] = self.timestamps[rows]
        range['data'] = self.data[rows,:]
        
        return Series(series=range)

def FWHM(X,Y):
    half_max = (max(Y)+min(Y)) / 2
    d = np.sign(half_max - np.array(Y[0:-1])) - np.sign(half_max - np.array(Y[1:]))
    #find the left and right most indexes
    left_idx = find(d > 0)[0]
    right_idx = find(d < 0)[-1]
    return X[right_idx], X[left_idx], half_max #return the difference (full width)

	def add_trials_from_file(self, filename):
		if self.span_type=='period' and filename:
			bci_stream=FileReader.bcistream(filename)
			sig,states=bci_stream.decode(nsamp='all')
			sig,chan_labels=bci_stream.spatialfilteredsig(sig)
			erpwin = [int(bci_stream.msec2samples(ww)) for ww in bci_stream.params['ERPWindow']]
			x_vec = np.arange(bci_stream.params['ERPWindow'][0],bci_stream.params['ERPWindow'][1],1000/bci_stream.samplingfreq_hz,dtype=float)
			trigchan = bci_stream.params['TriggerInputChan']
			trigchan_ix = find(trigchan[0] in chan_labels)
			trigthresh = bci_stream.params['TriggerThreshold']
			trigdetect = find(np.diff(np.asmatrix(sig[trigchan_ix,:]>trigthresh,dtype='int16'))>0)+1			
			intensity_detail_name = 'dat_TMS_powerA' if self.detail_values.has_key('dat_TMS_powerA') else 'dat_Nerve_stim_output'
			#Get approximate data segments for each trial
			trig_ix = find(np.diff(states['Trigger'])>0)+1
			for i in np.arange(len(trigdetect)):
				ix = trigdetect[i]
				dat = sig[:,ix+erpwin[0]:ix+erpwin[1]]
				self.trials.append(Datum(subject_id=self.subject_id\
                                            , datum_type_id=self.datum_type_id\
                                            , span_type='trial'\
                                            , parent_datum_id=self.datum_id\
                                            , IsGood=1, Number=0))
				my_trial=self.trials[-1]
				my_trial.detail_values[intensity_detail_name]=str(states['StimulatorIntensity'][0,trig_ix[i]])
				if int(bci_stream.params['ExperimentType']) == 1:#SICI intensity
					my_trial.detail_values['dat_TMS_powerB']=str(bci_stream.params['StimIntensityB'])#TODO: Use the state.
					my_trial.detail_values['dat_TMS_ISI']=str(bci_stream.params['PulseInterval'])
				my_trial.store={'x_vec':x_vec, 'data':dat, 'channel_labels': chan_labels}
			Session.commit()

def reject_outliers(x,percent=10,side='both'):
    '''
    Reject outliers from data based on percentiles.

    Parameters
    ----------
    x : ndarary
        1D numeric array of data values
    percent : number
        percent between 0 and 100 to remove
    side : str
        'left' 'right' or 'both'. Default is 'both'. Remove extreme
        values from the left / right / both sides of the data
        distribution. If both, the percent is halved and removed
        from both the left and the right
    Returns
    -------
    ndarray
        Values with outliers removed
    kept
        Indecies of values kept
    removed
        Indecies of values removed
    '''
    N = len(x)
    remove = outliers(x,percent,side)
    to_remove = find(remove==True)
    to_keep   = find(remove==False)
    return x[to_keep], to_keep, to_remove

def make_batches(data, batch_size=100):
    """
    Split the data into minibatches of size batch_size
    
    This procedure generates subsamples ids for batches by only considering 
    features that are active in the minibatch
    
    Args:
        data - the data to be split into minibatches (must be rank 2)
        batch_size - the size of the minibatches 
        
    Returns:
        batches - a list: [(batch, subsample_ids) for batch in minibatchs]
    """
    n = data.shape[0]
    perm = random.permutation(range(n))
    i = 0
    batches = []
    while i < n:
        batch = perm[i:i+batch_size]
        i += batch_size
        batches.append(data[batch])
    try:
        ids = [find((b.sum(0) != 0).A.flatten()) for b in batches]
    except AttributeError:
        ids = [find((b.sum(0) != 0).flatten()) for b in batches]
    batches = [(b[:,i].toarray(), i) for b,i in zip(batches, ids)]
    return batches

def dct_cut_downsampled(data,cutoff,spacing=0.4):
    '''
    like dctCut but also lowers the sampling rate, creating a compact
    representation from which the whole downsampled data could be
    recovered
    '''
    h,w    = np.shape(data)[:2]
    f1     = fftfreq(h*2,spacing)
    f2     = fftfreq(w*2,spacing)
    wl     = 1./abs(reshape(f1,(2*h,1))+1j*reshape(f2,(1,2*w)))
    mask   = int32(wl>=cutoff)
    mirror = reflect2D(data)
    ff     = fft2(mirror,axes=(0,1))
    cut    = (ff.T*mask.T).T # weird shape broadcasting constraints
    empty_cols = find(all(mask==0,0))
    empty_rows = find(all(mask==0,1))
    delete_col = len(empty_cols)/2 #idiv important here
    delete_row = len(empty_rows)/2 #idiv important here
    keep_cols  = w-delete_col
    keep_rows  = h-delete_row
    col_mask = np.zeros(w*2)
    col_mask[:keep_cols] =1
    col_mask[-keep_cols:]=1
    col_mask = col_mask==1
    row_mask = np.zeros(h*2)
    row_mask[:keep_rows] =1
    row_mask[-keep_rows:]=1
    row_mask = row_mask==1
    cut = cut[row_mask,...][:,col_mask,...]
    w,h = keep_cols,keep_rows
    result = ifft2(cut,axes=(0,1))[:h,:w,...]
    return result

    def cb_active_plot(self,start_ms,stop_ms,line_color='b'):
        """
        Plot timing information of time spent in the callback. This is similar
        to what a logic analyzer provides when probing an interrupt.

        cb_active_plot( start_ms,stop_ms,line_color='b')
        """
        # Find bounding k values that contain the [start_ms,stop_ms]
        k_min_idx = mlab.find(np.array(self.DSP_tic)*1000 < start_ms)
        if len(k_min_idx) < 1:
            k_min = 0
        else:
            k_min = k_min_idx[-1]
        k_max_idx = mlab.find(np.array(self.DSP_tic)*1000 > stop_ms)
        if len(k_min_idx) < 1:
            k_max= len(self.DSP_tic)
        else:
            k_max = k_max_idx[0]
        for k in range(k_min,k_max):
            if k == 0:
                plt.plot([0,self.DSP_tic[k]*1000,self.DSP_tic[k]*1000,
                         self.DSP_toc[k]*1000,self.DSP_toc[k]*1000],
                        [0,0,1,1,0],'b')
            else:
                plt.plot([self.DSP_toc[k-1]*1000,self.DSP_tic[k]*1000,self.DSP_tic[k]*1000,
                          self.DSP_toc[k]*1000,self.DSP_toc[k]*1000],[0,0,1,1,0],'b')
        plt.plot([self.DSP_toc[k_max-1]*1000,stop_ms],[0,0],'b')
        
        plt.xlim([start_ms,stop_ms])
        plt.title(r'Time Spent in the callback')
        plt.ylabel(r'Timing')
        plt.xlabel(r'Time (ms)')
        plt.grid();

def __get_level_index_from_chan_slot__(slotnum, channel):
    """ given slot and channel, returns corresponding the level|keep column index """
    global slot
    global __boardTypes__
    # 1. determine board type from slot
    for boardname in slot.keys():
        if slotnum in slot[boardname]:
            # 1a. for driver channels, multiply the chan number by 2
            # to get the UniqueStateArr index
            if boardname == 'drvr':
                channel *= 2
            # 1b. check that channel is valid for board type.
            if channel >= __chan_per_board__[boardname]:
                print "*** INVALID channel (%d) specified for slot (%d,%s) ***"%(channel,slotnum,boardname)
                return -1
            # 2. determine index in global slot for board type (0 for the first board of that type, etc)
            indx_slot = mlab.find(np.array(slot[boardname]) == slotnum)[0]
            # 3. calculate the base index (boardname index * slot index)
            signalPartitions = np.cumsum([ 0,
                                           __chan_per_board__['drvr'] * len(slot['drvr']), ## !driver-speed-keep
                                           __chan_per_board__['lvds'] * len(slot['lvds']),
                                           __chan_per_board__['adc']  * len(slot['adc']),
                                           __chan_per_board__['back'] * len(slot['back']),
                                           __chan_per_board__['hvbd'] * len(slot['hvbd']) ])
            indx_LVL_boardname = mlab.find(np.array(__boardTypes__) == boardname)[0]
            indx_base = signalPartitions[indx_LVL_boardname] + indx_slot * __chan_per_board__[boardname]
            # 4. add the channel offset
            return (indx_base + channel)

def sudokuChecker(matrix):
    '''
    sudokuChecker takes a 9x9 sudoku and return True if it is a valid sudoku,
    i.e. it meets all the rules.
    
    Parameters
    ----------   
    matrix: array_like
            `matrix` is a numpy.array that contains a 9x9 sudoku.
    
    Returns
    -------
    valid : bool
        If any of the numbers in `matrix` does not follow the 3 sudoku's rules, 
        `out` will be False. Otherwise will be True.
    '''
    
    valid = True
    
    if len( np.hstack((find(matrix <1), find(matrix > 9))) ) != 0:
        valid = False
        return valid
        
    for i in range(0,9):
        for j in range(0,9):
            
            k = matrix[i,j]
            
            if len(find(matrix[i,:] == k)) > 1 or \
               len(find(matrix[:,j] == k)) > 1 or \
               len(find(mySquare(matrix,i,j) == k)) > 1:
                   
                   valid = False
                   
        return valid

def relabel(complex_objs, old_tagging, majority=True):
    '''flatten+label(majority or consistent)'''
    val_map={}
    for i,obj in enumerate(complex_objs):
        tag= old_tagging[i]
        for item in obj: #each entity
            if not val_map.has_key(item):
                val_map[item]={}
            if not val_map[item].has_key(tag):
                val_map[item][tag]=0
            val_map[item][tag]+=1
    blorf= [[a] for a,counts in val_map.items() if len(counts.values())==1 or #it is not the case that there is full equality for all tags, so it's 1-1-1 or 2-2-2
        sum(counts.values()) % len(counts.values()) != 0]
    items= array(blorf, dtype=object)
    
    tags=[]
    for i,item in enumerate(items):
        label_counts=val_map[item[0]]
        #label counts is mapping of tag->apperances
        if majority:
            tags.append(label_counts.keys()[argmax(label_counts.values())])
        else:
            vals=array(label_counts.values())
            inds= find(vals>0)
            if len(inds) == 1: #exactly one class
                tags.append(array(label_counts.keys())[inds])
            else:
                tags.append(-1)
    tags=array(tags)
    if majority:
        return items,tags
    idxs=find(tags>=0)
    return items[idxs], tags[idxs]

def extr(x):
    """Extract the indices of the extrema and zero crossings.
    :param x: input signal
    :type x: array-like
    :return: indices of minima, maxima and zero crossings.
    :rtype: tuple
    """
    m = x.shape[0]

    x1 = x[:m - 1]
    x2 = x[1:m]
    indzer = find(x1 * x2 < 0)
    if np.any(x == 0):
        iz = find(x == 0)
        indz = []
        if np.any(np.diff(iz) == 1):
            zer = x == 0
            dz = np.diff([0, zer, 0])
            debz = find(dz == 1)
            finz = find(dz == -1) - 1
            indz = np.round((debz + finz) / 2)
        else:
            indz = iz
        indzer = np.sort(np.hstack([indzer, indz]))

    indmax = argrelmax(x)[0]
    indmin = argrelmin(x)[0]

    return indmin, indmax, indzer

	def find_spinning(mag,gyros):
		""" return the indicies in the magnetometer data when
		the gyro indicates it is spinning on the z axis """

		import numpy
		import scipy.signal
		from matplotlib.mlab import find

		threshold = 40
		spinning = scipy.signal.medfilt(abs(gyros['z'][:,0]),kernel_size=5) > threshold

		# make sure to always find end elements
		spinning = numpy.concatenate((numpy.array([False]),spinning,numpy.array([False])))
		start = find(spinning[1:] & ~spinning[0:-1])
		stop = find(~spinning[1:] & spinning[0:-1])-1

		tstart = gyros['time'][start]
		tstop = gyros['time'][stop]

		idx = numpy.zeros((0),dtype=int)
		for i in arange(tstart.size):

			i1 = abs(mag['time']-tstart[i]).argmin()
			i2 = abs(mag['time']-tstop[i]).argmin()
			
			idx = numpy.concatenate((idx,arange(i1,i2,dtype=int)))

		return idx

def split_and_subtree(query_chosen, recursive_step_obj):
    query_results=array([query_chosen(x) for x in recursive_step_obj.objects])
    pos_inds=find(query_results==1)
    neg_inds=find(query_results!=1)
    recursive_step_obj.left_son= TreeRecursiveSRLStep(recursive_step_obj.objects[neg_inds], recursive_step_obj.tagging[neg_inds], recursive_step_obj.relations, recursive_step_obj.transforms, recursive_step_obj.n, recursive_step_obj.MAX_DEPTH, recursive_step_obj.SPLIT_THRESH, recursive_step_obj.cond)
    recursive_step_obj.right_son=TreeRecursiveSRLStep(recursive_step_obj.objects[pos_inds], recursive_step_obj.tagging[pos_inds], recursive_step_obj.relations, recursive_step_obj.transforms, recursive_step_obj.n, recursive_step_obj.MAX_DEPTH, recursive_step_obj.SPLIT_THRESH, recursive_step_obj.cond)
    return query_chosen,recursive_step_obj.left_son,recursive_step_obj.right_son

def durations():
    
    import os
    import re
    from matplotlib.mlab import find
    
    save_file = 'durations.pkl'
    save_dir = '/home/mat/Dropbox/Working-memory/Metadata/'
    
    # I don't like regex...
    datadir = '/media/hippocampus/NDAQ'
    filelist = os.listdir(datadir)
    ML = find(['ML' in filename for filename in filelist])
    data = [ filelist[ind] for ind in ML ]
    # Now I have a list of the filenames for all of my data
    
    # But, I only want the ns5 data, so do the same thing again.
    ns5 = find(['.ns5' in filename for filename in data])
    data = [ data[ind] for ind in ns5 ]
    
    # Let's check if we've already done this for some data so we
    # don't do it twice.
    
    check_files = os.listdir(save_dir)
    if save_file in check_files:
        with open(save_file,'r') as f:
            info_list = pkl.load(f)
        for file in data:
            reg=re.search('datafile_ML_(\w+)_(\w+)_001.ns5', file)
            rat = reg.group(1)
            date = reg.group(2)
            
            for info in info_list:
                if (info['rat'] == rat) & (info['date'] == date):
                    info_list.remove(info)
    else:
        info_list = []

    for file in data:
        filepath = os.path.join(datadir,file)
        print "Loading %s" % file
        loader = ns5.loader(filepath)
        loader.load_header()
        print "File %s loaded" % file
        samples = loader.header.n_samples
        sample_rate = 30000
        
        duration = samples/float(sample_rate)
        print 'Duration is %s' % duration
        
        reg=re.search('datafile_ML_(\w+)_(\w+)_001.ns5', file)
        rat = reg.group(1)
        date = reg.group(2)
        
        info = {'rat':rat, 'date':date, 'duration':duration}
        info_list.append(info)
    
    with open(''.join([save_dir,save_file]), 'w') as f:
        pkl.dump(info_list, f)

def get_valid_trials(session, area=None):
    if dowarn():
        print("TRIAL IS 1 INDEXED FOR MATLAB COMPATIBILITY")
    if area is None:
        trials = set(find(metaloadvariable(session, "M1", "eventsByTrial")[:, 0]) + 1)
        trials &= set(find(metaloadvariable(session, "PMv", "eventsByTrial")[:, 0]) + 1)
        trials &= set(find(metaloadvariable(session, "PMd", "eventsByTrial")[:, 0]) + 1)
        return np.array(list(trials))
    else:
        return np.array(find(1 == metaloadvariable(session, area, "eventsByTrial")[:, 0]))

def tight(I):
	B  = I>0
	if I.ndim==2:
		ix = B.any(1);   ix0,ix1 = min(find(ix)) , max(find(ix))+1
		iy = B.any(0);   iy0,iy1 = min(find(iy)) , max(find(iy))+1
		return I[ix0:ix1,iy0:iy1]
	elif I.ndim==3:
		ix = B.max(2).any(1);   ix0,ix1 = min(find(ix)) , max(find(ix))+1
		iy = B.max(2).any(0);   iy0,iy1 = min(find(iy)) , max(find(iy))+1
		iz = B.max(0).any(0);   iz0,iz1 = min(find(iz)) , max(find(iz))+1
		return I[ix0:ix1,iy0:iy1,iz0:iz1]

def fortgaugeread (datafile="fort.gauge",setgaugefile="setgauges.data"):
    """
    Read data from fort.gauge files output by GeoClaw.

    Reads the gauge data and returns a list with mgauge elements.
    Each element of the list is a dictionary containing the data for a particular gauge.

    example:

    for N gauges: allgaugedata=[dictionary_1,...,dictionary_N] where
    dictionary_n.keys() = gauge#,x,y,level,t,q1,...qmq+1.
    where level,t,q1...q mq+1 are numpy arrays. 'q mq+1' is the surface elevation.
    """

    fid=open(setgaugefile)
    inp='#'
    while inp == '#':
        inpl=fid.readline()
        inp=inpl[0]

    inp = fid.readline()
    mgauges=int(inp.split()[0])
    gaugelocs=[]
    linesread=0
    while linesread < mgauges :
        row=string.split(fid.readline())
        if row!=[]:
           gaugelocs.append(row)
           linesread=linesread+1

    fid.close()

    data=loadtxt(datafile)

    allgaugedata=[]

    for n in xrange(mgauges) :
        onegaugedata=data[mlab.find(data[:,0]==int(gaugelocs[n][0]))]
        dict={}
        dict['gauge']=int(gaugelocs[n][0])
        dict['x']=float(gaugelocs[n][1])
        dict['y']=float(gaugelocs[n][2])
        onegaugedata = data[mlab.find(data[:,0]==dict['gauge'])]
        dict['level']=onegaugedata[:,1]
        dict['t']=onegaugedata[:,2]
        dict['mq'] = len(onegaugedata[0])-3
        for m in xrange(dict['mq']) :
            dict['q'+str(m+1)]=onegaugedata[:,3 + m]

        allgaugedata.append(dict)

    return allgaugedata

def CheckElevation(vec):
    #print vec
    vec_diff =  np.diff(np.array(vec))
    elChange = vec[-1] - vec[1];
    if (elChange > ELEVATION_THRESHOLD):
        #print 'DOWN'
        if (mlab.find(vec_diff >= 0).size > VECTOR_SIZE*ELEVATION_TREND_PRECENT):
            return elChange,HAND_DOWN
    elif (elChange < -ELEVATION_THRESHOLD):
        #print 'UP'
        if (mlab.find(vec_diff <= 0).size > VECTOR_SIZE*ELEVATION_TREND_PRECENT):
            return elChange,HAND_UP
    return 0,HAND_UNKWON