def sort_on_centroids(fn_mwk, fn_nev,k_cent_all,pc_sd_all,T_all,sel=SEL_ELEC_A, nmax=N_SNIPPET_MAX, npts=N_SNIPPET_PTS,c_v=VOLT_CONV):
	#this spikes sorts based on how close a spike is to the farthest centroid from center of mass of all the centroids
	sorted_data = {}
	for arg in getspk(fn_mwk, fn_nev, override_elecs=sel):
			# -- preps
			w = np.array(arg['wav']['unsorted'])
			w *= c_v
			ch = arg['ch']
			if ch not in sorted_data.keys():
				sorted_data[ch] = {'good_sp':[],'bad_sp':[],'g_sp_t_abs':[],'g_sp_t_rel':[],'g_sp_t_imgonset':[],'g_sp_im_id':[],'g_sp_im_ind':[],'b_sp_t_abs':[],'b_sp_t_rel':[],'b_sp_t_imgonset':[],'b_sp_im_id':[],'b_sp_im_ind':[]}
			ch_ind = sel.index(ch)
			#pick the centroid and the spikes sd for that channel
			pc_sd = pc_sd_all[ch_ind]
			k_cent = k_cent_all[ch_ind]
			k_dist = np.array([fastnorm(k-k_cent.mean(axis=0)) for k in k_cent])
			T = T_all[ch_ind]
			b_c = np.flatnonzero(k_dist == k_dist.max())
			#compute the distance of the waveform from the the centroids and determine whether it is a good spike or a bad spike.
			w_dist = np.array([fastnorm((np.dot(w,T)/pc_sd)-k_cent[k]) for k in range(len(k_cent))])
			if  np.flatnonzero(w_dist == w_dist.min())==b_c and k_dist.max() > 1.75:
				sorted_data[ch]['good_sp'].append(w)
				sorted_data[ch]['g_sp_t_abs'].append(arg['t_abs'])
				sorted_data[ch]['g_sp_t_rel'].append(arg['t_rel'])
				sorted_data[ch]['g_sp_t_imgonset'].append(arg['t_imgonset'])
				sorted_data[ch]['g_sp_im_id'].append(arg['imgid'])
				sorted_data[ch]['g_sp_im_ind'].append(arg['iimg'])
			else:
				sorted_data[ch]['b_sp_t_abs'].append(arg['t_abs'])
				sorted_data[ch]['b_sp_t_rel'].append(arg['t_rel'])
				sorted_data[ch]['b_sp_t_imgonset'].append(arg['t_imgonset'])
				sorted_data[ch]['b_sp_im_id'].append(arg['imgid'])
				sorted_data[ch]['b_sp_im_ind'].append(arg['iimg'])
				
	return sorted_data

def get_cap_check_indices(i):
    # Assumes that there is a test pulse followed by the stimulus pulses (downward first)
    di = np.diff(i)
    up_idx = np.flatnonzero(di > 0)
    down_idx = np.flatnonzero(di < 0)
    
    return up_idx[2::2], down_idx[1::2]

def align(s_x, s_y, mode="front"):
    """Function to align the input series

    :param s_x: First Series
    :type s_x: :mod:`pandas.Series`

    :param s_y: Second Series
    :type s_y: :mod:`pandas.Series`

    :param mode: Align Front/Back
    :type mode: str
    """

    p_x = np.flatnonzero(s_x)
    p_y = np.flatnonzero(s_y)

    if not len(p_x) or not len(p_y):
        return s_x, s_y, 0

    if mode == "front":
        p_x = p_x[0]
        p_y = p_y[0]

    if mode == "back":
        p_x = p_x[-1]
        p_y = p_y[-1]

    shift = p_x - p_y

    s_x, s_y = shift_series(s_x, s_y, shift)
    return s_x, s_y, shift

    def connection_field_plot_continuous(self, index, afferent=True, density=30):
        weights = numpy.array(self.proj.get('weight', format='list', gather=True))
        if afferent:
            idx = numpy.array(numpy.flatnonzero(weights[:,1].flatten()==index))
            x = self.proj.pre.positions[0][weights[idx,0].astype(int)]
            y = self.proj.pre.positions[1][weights[idx,0].astype(int)]
            w = weights[idx,2]
        else:
            idx = numpy.flatnonzero(weights[:,0]==index)
            x = self.proj.post.positions[0][weights[idx,1].astype(int)]
            y = self.proj.post.positions[1][weights[idx,1].astype(int)]
            w = weights[idx,2]

        xi = numpy.linspace(min(x), max(x), 100)
        yi = numpy.linspace(min(y), max(y), 100)
        zi = griddata(x, y, w, xi, yi)
        pylab.figure()
        #pylab.imshow(zi)
        pylab.scatter(x,y,marker='o',c=w,s=50)
        pylab.xlim(-self.source.parameters.sx/2,self.source.parameters.sx/2)
        pylab.ylim(-self.source.parameters.sy/2,self.source.parameters.sy/2)
        pylab.colorbar()
        pylab.title('Connection field from %s to %s of neuron %d' % (self.source.name,
                                                                     self.target.name,
                                                                     index))

def _hausdoff_ab(a, b, one_ring_neighbour):
    """
    Compute hausdoff distance of h(a,b)
    part unit of function hausdoff_distance
    
    Parameters:
    -----------
    a: array with 1 label
    b: array with 1 label
    one_ring_neighbour: one ring neighbour matrix

    Return:
    -------
    h: hausdoff(a,b)

    """
    a = np.array(a)
    b = np.array(b)
    h = 0
    for i in np.flatnonzero(a):
        hd = np.inf
        for j in np.flatnonzero(b):
            d = surf_dist(i,j, one_ring_neighbour)    
            if d<hd:
                hd = copy.deepcopy(d)
        if hd>h:
            h = hd
    return h

    def retrieve_coverage_information(self, cov_src=""):
        u"""
        Retrieve coverage information, i.e. whether locations to have spectra
        extracted have valid data and are not covered by clouds, margins etc.
        Information is retrieved from an external file.
        """
        self.coverage = dict()
        if cov_src and os.path.isfile(cov_src):
            for img_id in self.image_data:
                self.coverage[img_id] = list()
            try:
                data = np.loadtxt(cov_src)
            except:
                data = np.genfromtxt(cov_src, converters={0: lambda s: s.strip()})

            for d in data:
                # extracting row_id, usually a location id, i.e. for a plot
                try:
                    row_id = int(d[0])
                except:
                    row_id = d[0]
                # retrieving non-zero elements for each row, i.e. indicators
                # for visibility and coverage for each image id
                try:
                    nz = np.flatnonzero(d[1:])
                except:
                    # this line is necessary to deal with zero-rank arrays
                    # that are generated by np.genfromtxt above if the
                    # original location ids cannot be converted to float,
                    # i.e. if they contain letters
                    nz = np.flatnonzero(np.array(d.tolist()[1:]))
                for n in nz:
                    self.coverage[n].append(row_id)

    def __ComputeHistogram__(self, x, x_range):
        """Function to compute the histogram of the data
        
        Parameters
        ----------
        x: ndarray
            Array containing the data
        x_range: tuple
            Tuple containing the minimum and maxumum to consider
            to build the histogram of the data x

        Returns
        -------
        pdf: 1d-array 
            Return the PDF of x using the range provided by the user
        bin_edge: 1d-array
            Return the bins corresponding to the PDF
    
        """        
        # Compute the histogram for the data x with unit bins
        pdf_rel, bin_edges_rel = np.histogram(x, bins=(np.max(x) - np.min(x)), density=True)

        # We need to translate the pdf depending of the range given by x_range
        ### Create an array with unit bins depending of x_range
        ### We need max - min + 1 bins
        pdf_abs = np.zeros((x_range[1] - x_range[0],))
        bin_edges_abs = np.array(range(x_range[0], x_range[1] + 1))
        ### Copy the relative pdf at the right position
        pdf_abs[np.flatnonzero(bin_edges_abs==bin_edges_rel[0])[0] : np.flatnonzero(bin_edges_abs==bin_edges_rel[-1])[0]] = pdf_rel[:]

        return (pdf_abs, bin_edges_abs)

def do_pca_analysis(profiles, lens, name='', plot=False):
	L = np.array(0.446*(lens-np.mean(lens)), dtype='float64')
	pr = []
	for i,p in enumerate(profiles):
		mask = np.isnan(p)
		p[mask] = np.interp(np.flatnonzero(mask), np.flatnonzero(~mask), p[~mask])
		av, va = moving_average(np.log(p+0.001), 46, 100)
		pr.append(av)
	y = np.array(pr)
	pca = PCA(n_components=2)
	pca.fit(y)
	print pca.explained_variance_ratio_
	yp = pca.transform(y)
	m,b,r,p,_ = stats.linregress(L, yp[:,0])
	p1 = [p]
	r1 = [r]
	for _ in xrange(300):
		sample = np.random.choice(L.shape[0], L.shape[0], replace=True)
		m,b,r,p,_ = stats.linregress(L[~sample], yp[~sample,0])
		p1.append(p)
		r1.append(r)
	m,b,r,p,_ = stats.linregress(L, yp[:,1])
	p2 = [p]
	r2 = [r]
	for _ in xrange(300):
		sample = np.random.choice(L.shape[0], L.shape[0], replace=True)
		m,b,r,p,_ = stats.linregress(L[~sample], yp[~sample,1])
		p2.append(p)
		r2.append(r)
	if plot:
		plot_pca(y, pca, yp, L, name)
	return r1, p1, r2, p2, L.shape[0], name, np.std(L)

def fixgaps(x):
    """FIXGAPS: Linearly interpolates gaps in a time series
     YOUT=FIXGAPS(YIN) linearly interpolates over NaN in the input time
     series (may be complex), but ignores trailing and leading NaNs.
     R. Pawlowicz 11/6/99
     Version 1.0
    """
    

    #find nans
    bd = np.isnan(x)

    #early exit if there are no nans  
    if not bd.any():
        return x
    
    #find nonnans index numbers
    gd = np.flatnonzero(~bd)

    #ignore leading and trailing nans
    bd[:gd.min()]=False
    bd[(gd.max()+1):]=False
    
    #interpolate nans
    x[bd] = np.interp(np.flatnonzero(bd),gd,x[gd])

    return x

def preprocessRawData(raw_data,y_raw):
    # if age >= 38, labeled with 1
    # else, labeled with -1.
    labelvec = np.array(y_raw)
    y = np.ones(len(labelvec))
    neg = labelvec < 38
    y[neg] = -1 

    # positive examples: larger than or equal to 38
    # negative examples: smaller than 38
    num_pos = len(np.flatnonzero(y > 0))
    num_neg = len(np.flatnonzero(y < 0))
    #print('Number of positive/negative examples = %d/%d' % (num_pos, num_neg))

    headers = list(raw_data.columns.values) # get the features' name
    raw_feat = np.array(raw_data[headers]) # feature matrix without age feature  

    avg = np.mean(raw_feat,axis = 0)
    std_dev = np.std(raw_feat, axis = 0)


    X = (raw_feat-avg)/std_dev # scaled features matrix [-1,1]
    #print X.shape # X is N x n, where N is 342 subjects, n is 7 features
    #print y.shape # y is N x 1, where N is 342 subjects
    return X,y

def get_square_stim_characteristics(i, t, no_test_pulse=False):
    '''
    Identify the start time, duration, amplitude, start index, and
    end index of a square stimulus.
    This assumes that there is a test pulse followed by the stimulus square.
    '''

    di = np.diff(i)
    up_idx = np.flatnonzero(di > 0)
    down_idx = np.flatnonzero(di < 0)

    idx = 0 if no_test_pulse else 1

    # second square is the stimulus
    if up_idx[idx] < down_idx[idx]: # positive square
        start_idx = up_idx[idx] + 1 # shift by one to compensate for diff()
        end_idx = down_idx[idx] + 1
    else: # negative square
        start_idx = down_idx[idx] + 1
        end_idx = up_idx[idx] + 1

    stim_start = float(t[start_idx])
    stim_dur = float(t[end_idx] - t[start_idx])
    stim_amp = float(i[start_idx])

    return (stim_start, stim_dur, stim_amp, start_idx, end_idx)

def plotTrialPath(jsonData, trialNum, preTime=0, postTime=0, bPreRelStart=True, bPostRelEnd=True, bFolded=False, bLabels=True):
	if type(trialNum) == type(int()):
		trialNum = [trialNum]
	for tn in trialNum:
		trial = jsonData['trials'][tn]
		trialPath = getTrialPath(jsonData, tn, preTime, postTime, bPreRelStart, bPostRelEnd, bFolded)	
		if trial['bAvoidedShock']:
			pyplot.plot(trialPath[:,1],trialPath[:,2], 'g')
		else:
			pyplot.plot(trialPath[:,1],trialPath[:,2], 'r')
		#plot the LED
		if bFolded or trial['side']==0:
			pyplot.plot(0,35/2.0,marker='o',mec='r',mfc='r',ms=10)
			pyplot.plot(80,35/2.0,marker='o',mec='r',mfc='None',ms=10)
		else:
			pyplot.plot(0,35/2.0,marker='o',mec='r',mfc='None',ms=10)
			pyplot.plot(80,35/2.0,marker='o',mec='r',mfc='r',ms=10)	
		#plot dot where LED turned on.
		ndxLED = np.flatnonzero(trialPath[:,0]>0)
		if len(ndxLED)>0:
			pyplot.plot(trialPath[ndxLED[0],1],trialPath[ndxLED[0],2],marker='o',mec='r',mfc='r',ms=5)  
		#plot dot where Shock starts.
		ndxShock = np.flatnonzero(trialPath[:,0]>jsonData['parameters']['LEDTimeMS']/1000.0 - .1)
		if not trial['bAvoidedShock'] and len(ndxShock)>0:
			pyplot.plot(trialPath[ndxShock[0],1],trialPath[ndxShock[0],2],marker='o',mec='y',mfc='y',ms=5)  
		pyplot.text(trialPath[-1,1],trialPath[-1,2], str(tn))
	pyplot.xlim([0,80])
	pyplot.ylim([0,35])
	pyplot.axvline(x=getEscapePosition(jsonData),color='k',ls='--') 
	if bLabels:
		pyplot.xlabel('mm')
		pyplot.ylabel('mm')

def splint(xa, ya, y2a, x):
    '''spline interpolation'''
    try:
        len(x)
    except TypeError:
        x = np.array([x])

    try:
        klo, khi = np.array([
            (np.flatnonzero(xa < i)[-1], np.flatnonzero(xa > i)[0])
            for i in x
            ]).transpose()
    except IndexError:
        raise ValueError(
            'Input values must be between %s and %s'
            % (np.exp(xa[0]), np.exp(xa[-1]))
            )

    h = xa[khi] - xa[klo]
    if any(h <= 0):
        raise ValueError, 'xa input must be strictly increasing'
    a = (xa[khi] - x) / h
    b = (x - xa[klo]) / h

    res = (
        a * ya[klo]
        + b * ya[khi]
        + ((a**3 - a) * y2a[klo] + (b**3 - b) * y2a[khi]) * (h**2) / 6
        )
    return res

    def readPulsar(self, psr, psrname):
        print("WARNING: readPulsar has been deprecated!")
        psr.name = psrname

        # Read the content of the par/tim files in a string
        psr.parfile_content = str(self.getData(psrname, 'parfile', required=False))
        psr.timfile_content = str(self.getData(psrname, 'timfile', required=False))

        # Read the timing model parameter descriptions
        psr.ptmdescription = map(str, self.getData(psrname, 'tmp_name'))
        psr.ptmpars = np.array(self.getData(psrname, 'tmp_valpre'))
        psr.ptmparerrs = np.array(self.getData(psrname, 'tmp_errpre'))
        psr.flags = map(str, self.getData(psrname, 'efacequad', 'Flags'))

        # Read the position of the pulsar
        if self.hasField(psrname, 'raj'):
            psr.raj = np.float(self.getData(psrname, 'raj'))
        else:
            rajind = np.flatnonzero(np.array(psr.ptmdescription) == 'RAJ')
            psr.raj = np.array(self.getData(psrname, 'tmp_valpre'))[rajind]

        if self.hasField(psrname, 'decj'):
            psr.decj = np.float(self.getData(psrname, 'decj'))
        else:
            decjind = np.flatnonzero(np.array(psr.ptmdescription) == 'DECJ')
            psr.decj = np.array(self.getData(psrname, 'tmp_valpre'))[decjind]

        # Obtain residuals, TOAs, etc.
        psr.toas = np.array(self.getData(psrname, 'TOAs'))
        psr.toaerrs = np.array(self.getData(psrname, 'toaErr'))
        psr.prefitresiduals = np.array(self.getData(psrname, 'prefitRes'))
        psr.residuals = np.array(self.getData(psrname, 'postfitRes'))
        psr.detresiduals = np.array(self.getData(psrname, 'prefitRes'))
        psr.freqs = np.array(self.getData(psrname, 'freq'))
        psr.Mmat = np.array(self.getData(psrname, 'designmatrix'))

def __entrofy(X, k, w=None, q=None, pre_selects=None):
    '''See entrofy() for documentation'''

    n_participants, n_attributes = X.shape

    if w is None:
        w = np.ones(n_attributes)

    if q is None:
        q = 0.5 * np.ones(n_attributes)

    assert 0 < k <= n_participants
    assert not np.any(w < 0)
    assert np.all(q >= 0.0) and np.all(q <= 1.0)
    assert len(w) == n_attributes
    assert len(q) == n_attributes

    if k == n_participants:
        return np.arange(n_participants)

    # Initialization
    y = np.zeros(n_participants, dtype=bool)

    if pre_selects is None:
        # Select one at random
        pre_selects = np.random.choice(n_participants, size=1)

    y[pre_selects] = True

    # Where do we have missing data?
    Xn = np.isnan(X)

    while True:
        i = y.sum()
        if i >= k:
            break

        # Initialize the distribution vector
        p = np.nanmean(X[y], axis=0)
        p[np.isnan(p)] = 0.0

        # Compute the candidate distributions
        p_new = (p * i + X) / (i + 1.0)

        # Wherever X is nan, propagate the old p since we have no new information
        p_new[Xn] = (Xn * p)[Xn]

        # Compute marginal gain for each candidate
        delta = obj(p_new, w, q) - obj(p, w, q)

        # Knock out the points we've already taken
        delta[y] = -np.inf

        # Select the top score.  Break near-ties randomly.
        target_score = delta.max()
        target_score = target_score - 1e-3 * np.abs(target_score)
        new_idx = np.random.choice(np.flatnonzero(delta >= target_score))
        y[new_idx] = True

    return obj(np.nanmean(X[y], axis=0), w, q), np.flatnonzero(y)

def elam_spline(xin, yin, yspl_in, xout):
    """
    interpolate values from Elam photoabsorption and
    scattering tables, according to Elam, and following
    standard interpolation methods.  Calc borrowed from D. Dale.

    Parameters:
        xin (ndarray): x values for interpolation data
        yin (ndarray): y values for interpolation data
        yspl_in (ndarray): spline coefficients (second derivatives of y) for
                       interpolation data
        xout(float or ndarray): x values to be evaluated at

    Returns:
        ndarray: interpolated values
    """
    x = as_ndarray(xout)
    x[np.where(x < min(xin))] = min(xin)
    x[np.where(x > max(xin))] = max(xin)

    lo, hi = np.array([(np.flatnonzero(xin < e)[-1],
                        np.flatnonzero(xin > e)[0])
                       for e in x]).transpose()

    diff = xin[hi] - xin[lo]
    if any(diff <= 0):
        raise ValueError('x must be strictly increasing')
    a = (xin[hi] - x) / diff
    b = (x - xin[lo]) / diff
    return (a * yin[lo] + b * yin[hi] +
            (diff*diff/6) * ((a*a - 1) * a * yspl_in[lo] +
                             (b*b - 1) * b * yspl_in[hi]))

 def time_of_state_transition(self,prevStateID,nextStateID):
     '''Returns the time of the transition from prevStateID to nextStateID on each trial.
     If prevStateID is an empty string it calculates transition from any state.
     It assumes that the states have the same ID on all trials.
     If the transition did not occur, it returns NaN for that trial.
     If it occurs more than once, it returns the first time it happens.
     '''
     # -- If input is a state name, find its ID --
     if isinstance(nextStateID,str):
         nextStateIDval = self.find_stateID(nextStateID)
     else:
         nextStateIDval = nextStateID
     if len(prevStateID):
         USE_PREVSTATE = True
         if isinstance(prevStateID,str):
             prevStateIDval = self.find_stateID(prevStateID)
         else:
             prevStateIDval = prevStateID
     else:
         USE_PREVSTATE = False
     transitionEventTimes = np.empty(self['nTrials'])
     transitionEventTimes.fill(np.nan)
     eventsIndList = events_each_trial(self['RawEventsTrialID'])
     for trialInd,eventsInd in enumerate(eventsIndList):
         if USE_PREVSTATE:
             transitions = (self['RawEvents'][eventsInd,0]==prevStateIDval) &\
                           (self['RawEvents'][eventsInd,3]==nextStateIDval)
             transitionEventInds = np.flatnonzero(transitions)
         else:
             transitionEventInds = np.flatnonzero(self['RawEvents'][eventsInd,3]==nextStateIDval)
         if len(transitionEventInds):
             # Note that we store only the first transition
             eventOfInterest = eventsInd[transitionEventInds[0]]
             transitionEventTimes[trialInd] = self['RawEvents'][eventOfInterest,2]
     return transitionEventTimes

    def get_feature_mask(self, mask):
        """
        Parameters
        ----------
        mask: (N, F) or (F,) ndarray of bool
            True for unobserved features.

        Returns
        -------
        feature_mask: (N, D) or (D,) ndarray of bool
            True for unobserved features.
        """
        assert(mask.dtype == bool)
        if mask.ndim == 1:
            assert(mask.shape[0] == self.F)
            feature_mask = np.ones(self.D).astype(bool)
            for ind in np.flatnonzero(~mask):
                feature_mask[slice(*self.feature_bounds[ind])] = False
        else:
            assert(mask.shape[1] == self.F)
            N = mask.shape[0]
            feature_mask = np.ones((N, self.D)).astype(bool)
            for i in xrange(N):
                for ind in np.flatnonzero(~mask[i]):
                    feature_mask[i, slice(*self.feature_bounds[ind])] = False
        return feature_mask

def gen_crossed_logit_pandas(nc, cs, s1, s2):

    np.random.seed(3799)

    a = np.kron(np.arange(nc), np.ones(cs))
    b = np.kron(np.ones(cs), np.arange(nc))
    fe = np.ones(nc * cs)

    vc = np.zeros(nc * cs)
    for i in np.unique(a):
        ii = np.flatnonzero(a == i)
        vc[ii] += s1*np.random.normal()
    for i in np.unique(b):
        ii = np.flatnonzero(b == i)
        vc[ii] += s2*np.random.normal()

    lp = -0.5 * fe + vc
    pr = 1 / (1 + np.exp(-lp))
    y = 1*(np.random.uniform(size=nc*cs) < pr)

    ident = np.zeros(2*nc, dtype=np.int)
    ident[nc:] = 1

    df = pd.DataFrame({"fe": fe, "a": a, "b": b, "y": y})

    return df

  def _calculateSegmentActivity(connections, activeInput, connectedPermanence,
                                activationThreshold, minThreshold):
    """
    Calculate the active and matching segments for this timestep.

    @param connections (SparseMatrixConnections)
    @param activeInput (numpy array)

    @return (tuple)
    - activeSegments (numpy array)
      Dendrite segments with enough active connected synapses to cause a
      dendritic spike

    - matchingSegments (numpy array)
      Dendrite segments with enough active potential synapses to be selected for
      learning in a bursting column

    - potentialOverlaps (numpy array)
      The number of active potential synapses for each segment.
      Includes counts for active, matching, and nonmatching segments.
    """

    # Active
    overlaps = connections.computeActivity(activeInput, connectedPermanence)
    activeSegments = np.flatnonzero(overlaps >= activationThreshold)

    # Matching
    potentialOverlaps = connections.computeActivity(activeInput)
    matchingSegments = np.flatnonzero(potentialOverlaps >= minThreshold)

    return (activeSegments,
            matchingSegments,
            potentialOverlaps)