def precision_and_recall(actual,predicted,cls):
    c = (actual == cls)
    si = sp.argsort(-c)
    tp = sp.cumsum(sp.single(predicted[si] == cls))
    fp = sp.cumsum(sp.single(predicted[si] != cls))
    rec = tp /sp.sum(predicted == cls)
    prec = tp / (fp + tp)
    return prec,rec

def pca(dat, npca=None, verbose = False):
    if isinstance(dat, sp.ndarray):
        dat = pd.DataFrame(dat)
        names = []
        for i in range(dat.shape[1]):
            names.append("x"+str(i+1))
        dat.columns = names
    names = list(dat.columns)
    nr = dat.shape[0]
    nc = dat.shape[1]
    r = sp.corrcoef(dat, rowvar=False)
    heikin = dat.mean(axis=0)
    bunsan = dat.var(axis=0, ddof=1)
    sd = sp.sqrt(bunsan)
    eval, evec = linalg.eig(r)
    eval = sp.real(eval)
    rank = rankdata(eval, method="ordinal")
    rank = nc+1-rank
    eval2 = eval.copy()
    evec2 = evec.copy()
    for i in range(nc):
        j = sp.where(rank == i+1)[0][0]
        eval[i] = eval2[j]
        evec[:, i] = evec2[:, j]
    contr = eval/nc*100
    cum_contr = sp.cumsum(contr)
    fl = (sp.sqrt(eval)*evec)
    for i in range(nc):
        dat.ix[:, i] = (dat.ix[:, i]-heikin[i]) / sd[i]
    fs = sp.dot(dat, evec*sp.sqrt(nr/(nr-1)))
    if npca is None:
        npca = sp.sum(eval >= 1)
    eval = eval[0:npca]
    cont = eval/nc
    cumc = sp.cumsum(cont)
    fl = fl[:, 0:npca]
    rcum = sp.sum((fl ** 2), axis=1)
    if verbose:
        print("            ", end="")
        for j in range(npca):
            print("{0:>8s}".format("PC"+str(j+1)), end="")
        print("  Contribution")
        for i in range(nc):
            print("{0:>12s}".format(names[i]), end="")
            for j in range(npca):
                print(" {0:7.3f}".format(fl[i, j]), end="")
            print(" {0:7.3f}".format(rcum[i]))
        print("  Eigenvalue", end="")
        for j in range(npca):
            print(" {0:7.3f}".format(eval[j]), end="")
        print("\nContribution", end="")
        for j in range(npca):
            print(" {0:7.3f}".format(cont[j]), end="")
        print("\nCum.contrib.", end="")
        for j in range(npca):
            print(" {0:7.3f}".format(cumc[j]), end="")
        print()
    return {"r":r, "fl":fl, "eval":eval, "fs":fs[:, 0:npca]}

def fastunwrap(thetaArray, discont = scipy.pi):
    # takes an array of theta values
    # returns the data in unwrapped form (unwrapping over the axis == 1)
    diff = scipy.zeros_like(thetaArray)
    diff[1:,:] = scipy.diff(thetaArray, axis = 0)
    upSteps = diff > discont
    downSteps = diff < -discont
    shift = scipy.cumsum(upSteps, axis = 0) - scipy.cumsum(downSteps, axis = 0)
    return thetaArray - 2.0*discont*shift

def linspace_weighed(a, b, n, points):
    """Positions 'n' points in range ['a', 'b'] such that space around
       points[:][0] has an additional weight points[:][1] and
       half-width points[:][2].

       The shape of the weight is ``(|x - x0|/w + 1)**(-2)``, so that if the
       range is infinite, w is indeed the half-width of the distribution.
       points[:][1] describes the relative weights of such peaks."""

    x = linspace(a, b, 5*n)
    density = _n.zeros([len(x)], _n.float_)

    for point in points:
        point = list(point)
        
        if point[0] < a:
            point[0] = a
        if point[0] > b:
            point[0] = b

        density_shape = 1 / (abs(x - point[0])/point[2] + 1)**2
        base_weight = scipy.integrate.trapz(density_shape, x)

        density += (point[1] / base_weight) * density_shape

    if len(points) == 0:
        density[:] = 1

    cumdensity = scipy.cumsum(density) - density[0]
    cumdensity /= cumdensity[-1]
    interpolant = scipy.interpolate.interp1d(cumdensity, x)

    y = linspace(0, 1, n)
    return interpolant(y)

def continuous_phase(phase, axis=0, center=False):
    """Add and subtract 2 pi such that the phase in the array is
       as continuous as possible, along first or given axis. Optionally,
       it also centers the phase data so that the average is smallest."""

    phase = _n.array(phase, copy=0)

    rowshape = list(phase.shape)
    
    if len(rowshape) > 0:
        rowshape[axis] = 1

        slip = _n.concatenate([ _n.zeros(rowshape),
                                scipy.diff(phase, axis=axis) ],
                              axis=axis)
        slip = _n.around(slip/(2*_n.pi))
        cumslip = scipy.cumsum(slip, axis=axis)

        phase = phase - 2*_n.pi*cumslip
    else:
        pass

    if center:
        offset = _n.around(scipy.average(phase, axis=axis)/(2*_n.pi))
        offset = _n.reshape(offset, rowshape)
        offset = _n.repeat(offset, cumslip.shape[axis], axis=axis)
        phase = phase - 2*_n.pi*offset
    
    return phase

 def eliminatePercentileTails(self, mskDds, loPercentile=10.0, hiPercentile=90.0):
     """
     Trims lower and/or upper image histogram tails by replacing :samp:`mskDds`
     voxel values with :samp:`mskDds.mtype.maskValue()`. 
     """
     rootLogger.info("Eliminating percentile tails...")
     rootLogger.info("Calculating element frequencies...")
     elems, counts = elemfreq(mskDds)
     rootLogger.info("elems:\n%s" % (elems,))
     rootLogger.info("counts:\n%s" % (counts,))
     cumSumCounts = sp.cumsum(counts, dtype="float64")
     percentiles = 100.0*(cumSumCounts/float(cumSumCounts[-1]))
     percentileElems = elems[sp.where(sp.logical_and(percentiles > loPercentile, percentiles < hiPercentile))]
     loThresh = percentileElems[0]
     hiThresh = percentileElems[-1]
     rootLogger.info("Masking percentiles range (%s,%s) = (%s,%s)" % (loPercentile, hiPercentile, loThresh, hiThresh))
     mskDds.asarray()[...] = \
         sp.where(
             sp.logical_and(
                 sp.logical_and(mskDds.asarray() >= loThresh, mskDds.asarray() <= hiThresh),
                 mskDds.asarray() != mskDds.mtype.maskValue()
             ),
             mskDds.asarray(),
             mskDds.mtype.maskValue()
         )
     rootLogger.info("Done eliminating percentile tails.")

def recombine(G,rates):
    """
    Performs recombination of genotype matrix G corresponding to given rates

    input:      G       Nx2 matrix of integers where each columns give
                        haplotypes for N markers
                rates   (N-1)x1 vector of floats with recombination rates

    Element rates[i] is the rate of recombination between markers i and i+1

    Example:
    #>>> from numpy import random, ones, zeros, column_stack
    >>> G = column_stack((zeros((50,1)),ones((50,1))))
    >>> rates = 0.2*ones((G.shape[0],1))
    >>> random.seed(seed=0)
    >>> recombine(G,rates).T
    array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
             0.,  0.,  1.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,
             0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.,
             0.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,  0.,  0.],
           [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
             1.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,
             1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,  1.,  1.,
             1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.,  0.,  1.,  1.]])
   
    """
   
    #draw uniform numbers for all marker interval
    crossover = mod(cumsum(random.uniform(size=(G.shape[0],1))<rates),2)
   
    #generate recombined G
    recG = array([ [g[c],g[c-1]] for g,c in zip(G[1:G.shape[0],:],crossover) ])
    recG = vstack((G[0,:],recG))
    
    return recG

def normalizeHistogram(data):
    histogram = scipy.ndimage.histogram(data.astype("f"), 0, 255, 256)
    cumulatedHistogram = scipy.cumsum(histogram)
    nch = cumulatedHistogram.astype("f")/len(data.flat)
    inch = (nch*255).astype("i")
    normalize = scipy.vectorize(lambda i: inch[i])
    return normalize(data)

def pct_sigma(array):
    """
    Get normal quantiles

    Parameters
    ----------
    x : array_like
        distribtion of values

    Returns
    -------
    sigma : ndarray
      normal quantile
    pct : ndarray
      percentile
    y : ndarray
      value
    """
    qrank = lambda x: ((x - 0.3175)/(x.max() + 0.365))

    y = array.copy()
    y = y[~np.isnan(y)]
    y = np.sort(y)

    if y.size == 0:
        blank = np.zeros(y.shape)
        return blank, blank, blank

    n = sp.ones(len(y))
    cs = sp.cumsum(n)
    pct = qrank(cs)
    sigma = sps.norm.ppf(pct)

    return sigma, pct, y

    def initialize(self, state, chain):
        params = {}
        for key in self.scan_range.keys():
            # Check for single range
            if len(self.scan_range[key]) == 2:
                params[key] = sp.rand() * (self.scan_range[key][1] - self.scan_range[key][0]) + self.scan_range[key][0]
            else:
                # calculate weights of sub_regions
                sub_size = sp.array([])
                # Determine weights of region
                for i in range(0, len(self.scan_range[key]), 2):
                    sub_size = sp.append(sub_size, self.scan_range[key][i + 1] - self.scan_range[key][i])
                    self.range_weight[key] = sub_size / float(sp.sum(sub_size))

                # sample region based on size
                i_sel = 2 * sp.searchsorted(sp.cumsum(self.range_weight[key]), sp.rand())
                # sample point
                params[key] = (
                    sp.rand() * (self.scan_range[key][i_sel + 1] - self.scan_range[key][i_sel])
                    + self.scan_range[key][i_sel]
                )

        # params=dict([(key,sp.rand()*(self.scan_range[key][1]-self.scan_range[key][0])+self.scan_range[key][0]) for key in self.scan_range.keys() if type(self.scan_range[key])==list])

        # Add constant parameters
        for key in self.constants.keys():
            params[key] = self.constants[key]

        for key in self.functions.keys():
            params[key] = self.functions[key](params)

        modelid = "%i%01i" % (self.rank, 0) + "%i" % chain.accepted

        return params, modelid

def align_chain_to_seq(sequence,chain,verbose=False):
	#Build Polypeptides from the chains
	polypeptides = build_polypeptides(chain)
	
	#Can't be broken out into another function, because we need seq_lens
	contiguous_seqs = [single_pp.get_sequence().tostring() for single_pp in polypeptides]
	ATOM_joined_seq = ''.join(contiguous_seqs)
	
	seq_lens = [0] + [len(single_pp) for single_pp in polypeptides]


	#Figuring all of this out took days...
	#I am so tired of dealing with mapping various numberings around
	#I wish Biopython, especially Bio.pairwise2 had better documentation
	breaks = set(S.cumsum(seq_lens) )#TODO : Tear hair out GYAAAAA
	
	nogaps = lambda x,y: -2000 -200*y #There really should not be inserts with respect to the database sequence.

	def specificgaps(x,y):
		if x in breaks:#very minor penalty for gaps at breaks in the PDB structure, consider using 0
			return (0 -y)
		else:
			return (-2000 -200*y)#strongly discourage gaps anywhere else.
	
	alignments = __PW.align.globalxc(sequence.seq.tostring(),ATOM_joined_seq,nogaps,specificgaps)
	
	if verbose:
		#some output?
		for a in alignments:
			__stderr.write( __PW.format_alignment(*a) )
			__stderr.write('\n')

	return alignments

 def __init__(self,layers,gridOpts):
   ''' Initialize the grid using the given layers and grid options.
   '''
   segments = []
   qStart   =  scipy.inf
   qEnd     = -scipy.inf
   for layer in layers:
     if layer.isQuantum:
       d1 = dn = gridOpts.dzQuantum
       segments += [self.get_dz_segment(d1,dn,layer.thickness)]
       qStart = min(qStart,sum([len(seg) for seg in segments[:-1]]))
       qEnd   = max(qEnd,  sum([len(seg) for seg in segments]))
     elif gridOpts.useFixedGrid:
       d1 = dn = gridOpts.dz
       segments += [self.get_dz_segment(d1,dn,layer.thickness)]
     elif layer.thickness*gridOpts.dzCenterFraction > gridOpts.dzEdge:
       d1 = dn = gridOpts.dzEdge
       dc = gridOpts.dzCenterFraction*layer.thickness
       segments += [self.get_dz_segment(d1,dc,layer.thickness/2),
                    self.get_dz_segment(dc,dn,layer.thickness/2)]
     else:
       d1 = dn = gridOpts.dzEdge
       segments += [self.get_dz_segment(d1,dn,layer.thickness)]
   self.dz       = scipy.concatenate(segments)
   self.z        = scipy.concatenate(([0],scipy.cumsum(self.dz)))
   self.zr       = (self.z[:-1]+self.z[1:])/2
   self.znum     = len(self.z)
   self.rnum     = len(self.zr)
   self.gridOpts = gridOpts
   self.qIndex   = scipy.arange(qStart,qEnd+1)   # Wavefunction index
   self.qrIndex  = scipy.arange(qStart,qEnd)     # Quantum region index   

def impz(b, a=1):
    """Plot step and impulse response of an FIR filter.

    b : float
        Forward terms of the FIR filter.
    a : float
        Feedback terms of the FIR filter. (Default value = 1)

    From http://mpastell.com/2010/01/18/fir-with-scipy/

    Returns
    -------
    None

    """
    l = len(b)
    impulse = np.repeat(0., l)
    impulse[0] = 1.
    x = np.arange(0, l)
    response = sp.lfilter(b, a, impulse)
    plt.subplot(211)
    plt.stem(x, response)
    plt.ylabel('Amplitude')
    plt.xlabel(r'n (samples)')
    plt.title(r'Impulse response')
    plt.subplot(212)
    step = sp.cumsum(response)
    plt.stem(x, step)
    plt.ylabel('Amplitude')
    plt.xlabel(r'n (samples)')
    plt.title(r'Step response')
    plt.subplots_adjust(hspace=0.5)

 def apply_flow(self,flowrate):
     r'''
     Convert the invaded sequence into an invaded time for a given flow rate
     considering the volume of invaded pores and throats.
     
     Parameters
     ----------
     flowrate : float
         The flow rate of the injected fluid
         
     Returns
     -------
     Creates a throat array called 'invasion_time' in the Algorithm 
     dictionary
     
     '''
     P12 = self._net['throat.conns']  # List of throats conns
     a = self['throat.invasion_sequence']  # Invasion sequence
     b = sp.argsort(self['throat.invasion_sequence'])
     P12_inv = self['pore.invasion_sequence'][P12]  # Pore invasion sequence
     # Find if the connected pores were invaded with or before each throat
     P1_inv = P12_inv[:,0] == a
     P2_inv = P12_inv[:,1] == a
     c = sp.column_stack((P1_inv,P2_inv))  
     d = sp.sum(c,axis=1,dtype=bool)  # List of Pores invaded with each throat
     # Find volume of these pores
     P12_vol = sp.zeros((self.Nt,))
     P12_vol[d] = self._net['pore.volume'][P12[c]]
     # Add invaded throat volume to pore volume (if invaded)
     T_vol = P12_vol + self._net['throat.volume']
     # Cumulative sum on the sorted throats gives cumulated inject volume
     e = sp.cumsum(T_vol[b]/flowrate)
     t = sp.zeros((self.Nt,))
     t[b] = e  # Convert back to original order
     self._phase['throat.invasion_time'] = t

    def blobs(shape, porosity, blobiness=8):
        """
        Generates an image containing amorphous blobs

        Parameters
        ----------
        shape : list
            The size of the image to generate in [Nx, Ny, Nz] where N is the
            number of voxels

        blobiness : scalar
            Controls the morphology of the image.  A higher number results in
            a larger number of smaller blobs.

        porosity : scalar
            The porosity of the final image.  This number is approximated by
            the method so the returned result may not have exactly the
            specified value.

        """
        if sp.size(shape) == 1:
            shape = sp.full((3, ), int(shape))
        [Nx, Ny, Nz] = shape
        sigma = sp.mean(shape)/(4*blobiness)
        mask = sp.rand(Nx, Ny, Nz)
        mask = spim.gaussian_filter(mask, sigma=sigma)
        hist = sp.histogram(mask, bins=1000)
        cdf = sp.cumsum(hist[0])/sp.size(mask)
        xN = sp.where(cdf >= porosity)[0][0]
        im = mask <= hist[1][xN]
        return im

def calc_random_pattern_exponential(tau,tStart=0,tDuration=1000,tTotal=1000,channel=0,name=''):
    """Calculate a sequence of random state changes with intervals drawn from an
    exponential distribution.

    Parameters
    ----------
    tau: float
        the scale parameter of the exponential distribution (in ms)

    tStart: int
        the time point at which the random changes start (in ms)
        
    tDuration: int
        the total length of the time section in which random changes can occur
        (in ms)
        
    tTotal: int
        the total length of the pattern (in ms)
        
    channel: int
        the channel to switch
        
    name: str
        the name of the pattern
        
    Returns
    -------
    Pattern: RIOpattern
        The generated RIOpattern instance
        
    """
    
    # a little ugly but guarantees to run
    change_times = [0]
    while sp.cumsum(change_times)[-1] < tDuration:
        change_times.append(stats.distributions.expon.rvs(scale=tau))
    
    change_times = sp.cumsum(change_times[1:-1]).astype('int32') + tStart
    
    # to make sure it ends latest at tStart+tDuration
    if len(change_times) % 2 == 1:
        change_times = sp.concatenate((change_times,[tStart+tDuration]))
    
    state_vec = change_times2state_vec(change_times,tTotal)
    Pattern = RIOpattern(name=name, Pulses=States2RIOpulses(state_vec,channel), total_duration=tTotal)
    return Pattern, state_vec

def shapesToSections(shapes):
    sections = []
    sections.append(shapes['syn_wo'][0] * shapes['syn_wo'][1])
    sections.append(shapes['inter_ws'][0] * shapes['inter_ws'][1])
    sections.append(shapes['syn_wi'][0] * shapes['syn_wi'][1])
    sections.append(shapes['wi'][0] * shapes['wi'][1])
    sections.append(shapes['wo'][0] * shapes['wo'][1])
    return cumsum(sections)

def calc_rmat(rho):
    L = len(rho)
    mat = SP.zeros([L,L])
    cr = SP.cumsum(rho)
    for l1 in range(L):
        for l2 in range(l1, L):
            mat[l1,l2] = 0.5*(1-SP.exp(-2*(cr[l2] - cr[l1])))
    return mat

def calculateThreshold(image, coveragePercent):
    import scipy
    data = image.data
    histogram = scipy.histogram(data, len(scipy.unique(data)))
    cumsum = scipy.cumsum(histogram[0])
    targetValue = cumsum[-1] * coveragePercent
    index = scipy.argmin(scipy.absolute(cumsum - targetValue))
    threshold = histogram[1][index]
    return threshold * image.unit

 def prandom_walk_from_here(self,length=10,edge_filter=None):
     import scipy
     if (length<=0):
         return [self]
     if (not edge_filter):
         edge_filter=lambda x,y:1
     sm=map(lambda e:edge_filter(e,length)*e.strength(),self.outedges())
     rv=(random.random()*sum(sm))
     idx=sum((scipy.cumsum(sm)<rv).astype(int))
     return [self]+(self.outedges()[idx]).target().prandom_walk_from_here(length-1,edge_filter)