def polyhedra(self, wm):
        '''Iterates through the polyhedra that make up the closest volume to a certain vertex'''
        for p, facerow in enumerate(self.connected):
            faces = facerow.indices
            pts, polys = _ptset(), _quadset()
            if len(faces) > 0:
                poly = np.roll(self.polys[faces[0]], -np.nonzero(self.polys[faces[0]] == p)[0][0])
                assert pts[wm[p]] == 0
                assert pts[self.pts[p]] == 1
                pts[wm[poly[[0, 1]]].mean(0)]
                pts[self.pts[poly[[0, 1]]].mean(0)]

                for face in faces:
                    poly = np.roll(self.polys[face], -np.nonzero(self.polys[face] == p)[0][0])
                    a = pts[wm[poly].mean(0)]
                    b = pts[self.pts[poly].mean(0)]
                    c = pts[wm[poly[[0, 2]]].mean(0)]
                    d = pts[self.pts[poly[[0, 2]]].mean(0)]
                    e = pts[wm[poly[[0, 1]]].mean(0)]
                    f = pts[self.pts[poly[[0, 1]]].mean(0)]

                    polys((0, c, a, e))
                    polys((1, f, b, d))
                    polys((1, d, c, 0))
                    polys((1, 0, e, f))
                    polys((f, e, a, b))
                    polys((d, b, a, c))

            yield pts.points, np.array(list(polys.triangles))

 def compute_sketch(self):
     start_time = time.time()
     if self.sketch is not None:
         return self.sketch
     # basically, want to init an empty csr matrix 
     if self.randomized and (self.m > 100 * self.l):
         print "using sparse sketch"
         # lets use the sparse version of randomized sketch here 
         return self.compute_sparse_sketch()
     else:
         self._sketch_func = self._fast_rand_sketch
     # what do we do differently here? we need to iterate over the nzrow_inds,
     mat_b = np.zeros([self.l + self.b_size, self.m])
     # other way: np.where(~mat_b.any(axis=1))[0]
     # zero_rows = np.nonzero([round(s, 7) == 0.0 for s in np.sum(mat_b, axis = 1)])[0].tolist()
     zero_rows = np.nonzero([round(s, 7) == 0.0 for s in np.sum(mat_b[:self.l, :], axis = 1)])[0]
     zero_rows = np.hstack((zero_rows, np.arange(self.l, self.l + self.b_size))).tolist()
     # iterate through the nzrow_inds
     for i in self.nzrow_inds:
         mat_b[zero_rows[0], :] = self.mat.getrow(i).todense()
         zero_rows.remove(zero_rows[0])
         if len(zero_rows) == 0:
             print "sketching ", i
             self._sketch_func(mat_b)
             zero_rows = np.nonzero([round(s, 7) == 0.0 for s in np.sum(mat_b[:self.l, :], axis = 1)])[0]
             zero_rows = np.hstack((zero_rows, np.arange(self.l, self.l + self.b_size))).tolist()
     self._sketch_func(mat_b)
     self.sketch = mat_b[:self.l, :]
     self.sketching_time = time.time() - start_time 
     return self.sketch

	def init_from_ms(self, ms_run):
		'''
		Initialize the population from the results of an ms run to avoid having to do
		lengthy 'burn-in' phase at beginning of each simulation.
		'''

		if not isinstance(ms_run, ms.MsReader):
			raise TypeError("Argument 'ms_sample' must be an MsRun object.")

		sys.stderr.write("Initializing population from ms sample with header:\n{}\n".format(ms_run.header))

		# read in a simulation
		ms_sample = ms_run.next()

		# initialize position of mutations
		chrlen = self.chrlen[0]
		pos = np.array(ms_sample.positions, dtype = np.float32)*chrlen # scale by chromosome length

		assert(len(ms_sample.samples) >= 2*self.size)

		for i in range(0, self.size):
			alleles = np.array([ int(x) for x in ms_sample.samples[i] ])
			derived = np.nonzero(alleles)[0]
			if len(derived):
				self.chroms[0][i] = pos[ np.nonzero(alleles)[0] ]
			else:
				self.chroms[0][i] = np.ndarray((0,), dtype = np.float32)

def _generate_clique_alt(variables, obj, inequalities, equalities):
    n_dim = len(variables)
    rmat = spmatrix(1.0, range(n_dim), range(n_dim))
    for support in get_support(variables, obj):
        nonzeros = np.nonzero(support)[0]
        value = random.random()
        for i in nonzeros:
            for j in nonzeros:
                rmat[i, j] = value
    for polynomial in flatten([inequalities, equalities]):
        support = np.any(get_support(variables, polynomial), axis=0)
        nonzeros = np.nonzero(support)[0]
        value = random.random()
        for i in nonzeros:
            for j in nonzeros:
                rmat[i, j] = value
    rmat = rmat + 5*n_dim*spmatrix(1.0, range(n_dim), range(n_dim))
    # compute symbolic factorization using AMD ordering
    symb = cp.symbolic(rmat, p=amd.order)
    ip = symb.ip
    # symb = cp.symbolic(rmat)
    # ip = range(n_dim)
    cliques = symb.cliques()
    R = np.zeros((len(cliques), n_dim))
    for i, clique in enumerate(cliques):
        for j in range(len(clique)):
            R[i, ip[cliques[i][j]]] = 1
    return R

 def _get_ind_under_point(self, event):
     'get the index of the vertex under point if within epsilon tolerance'
     try:
         x, y = zip(*self._poly.xy)
         
         # display coords
         xt, yt = self._poly.get_transform().numerix_x_y(x, y)
         d = np.sqrt((xt-event.x)**2 + (yt-event.y)**2)
         indseq = np.nonzero(np.equal(d, np.amin(d)))
         ind = indseq[0]
     
         if d[ind]>=self._epsilon:
             ind = None
     
         return ind
     except:
         # display coords
         xy = np.asarray(self._poly.xy)
         xyt = self._poly.get_transform().transform(xy)
         xt, yt = xyt[:, 0], xyt[:, 1]
         d = np.sqrt((xt-event.x)**2 + (yt-event.y)**2)
         indseq = np.nonzero(np.equal(d, np.amin(d)))[0]
         ind = indseq[0]
         
         if d[ind]>=self._epsilon:
             ind = None
         
         return ind

def intersubjectconsensus():
    """Compute inter-subjects clustering consensus.

    """
    base_dir = r'/nfs/h1/workingshop/huanglijie/uni_mul_analysis'
    db_dir = os.path.join(base_dir, 'multivariate', 'detection', 'mvpcluster')

    n_clusters = 60

    mask_file = os.path.join(base_dir, 'multivariate', 'detection',
                             'mask.nii.gz')
    mask = nib.load(mask_file).get_data()

    for n in range(1, n_clusters):
        n += 1
        merged_file = os.path.join(db_dir, 'merged_cluster_'+str(n)+'.nii.gz')
        merged_data = nib.load(merged_file).get_data()
        n_subjs = merged_data.shape[3]
        mtx = np.zeros((n_subjs, n_subjs))
        for i in range(n_subjs):
            for j in range(n_subjs):
                data_i = merged_data[..., i]
                data_j = merged_data[..., j]
                vtr_i = data_i[np.nonzero(mask)]
                vtr_j = data_j[np.nonzero(mask)]
                tmp = metrics.adjusted_mutual_info_score(vtr_i, vtr_j)
                mtx[i, j] = tmp
        outfile = os.path.join(db_dir, 'consensus_'+str(n)+'.csv')
        np.savetxt(outfile, mtx, delimiter=',')

def rect(time, t0, t1, height=1, noise = 0.0):
	"""Rectangular signal of given height and width t1-t0

    Parameters
    ----------
        time : np.ndarray of shape (N,)
        	time instants (equidistant)
        t0 : float
        	time instant of rect lhs
        t1 : float
        	time instant of rect rhs
        height : float
         	signal maximum
        noise :float, optional
        	std of simulated signal noise

    Returns
    -------
        x : np.ndarray of shape (N,)
         	signal amplitudes at time instants
    """

	x = np.zeros((len(time),))
	x[np.nonzero(time > t0)] = height
	x[np.nonzero(time > t1)] = 0.0

	if noise > 0:
		x = x + np.random.randn(len(time)) * noise
	return x

def decTY1(raw_8, raw_16=None, raw_32=None):
    """
    Modified byte offset decompressor used in Oxford Diffraction images
    @param raw_8,raw_16,raw_32: strings containing raw data with integer of the given size
    @return numpy.ndarray 
    """
    data = numpy.fromstring(raw_8, dtype="uint8").astype(int)
    data -= 127
    if raw_32 is not None:
        int32 = numpy.fromstring(raw_32, dtype="int32").astype(int)
        exception32 = numpy.nonzero(data == 128)
    if raw_16 is not None:
        int16 = numpy.fromstring(raw_16, dtype="int16").astype(int)
        exception16 = numpy.nonzero(data == 127)
        data[exception16] = int16
    if raw_32:
        data[exception32] = int32
    summed = data.cumsum()
    smax = summed.max()
    if (smax > (2 ** 31 - 1)):
        bytecode = "int64"
    elif (smax > (2 ** 15 - 1)):
        bytecode = "int32"
    elif (smax > (2 ** 7 - 1)):
        bytecode = "int16"
    else:
        bytecode = "int8"
    return summed.astype(bytecode)

def d2init(p):

    # useful indexing for 2nd derivs

    # for difference terms
    n = (p - 2) * (p - 1) // 2
    Id = zeros((n, 2), dtype=int)
    n = -1
    for i in range(p - 1):
        for j in range(i + 1, p - 1):
            n += 1
            Id[n, :] = [i, j]

    # for permutation terms
    Ip = loop(zeros(p - 3, dtype=int), p, 0)
    loop(reset=True)

    # indexing for sums and products
    Jd = zeros((p - 1, p - 2), dtype=int)  # Jd[i,:] lists the rows of Id with i
    Jp = zeros((p - 1, p - 2), dtype=int)  # Jp[i,:] lists the rows of Ip without i
    mask = zeros(Ip.shape, dtype=int)
    for i in range(p - 1):
        (Jd[i, :], x) = nonzero(Id == i)
        (Jp[i, :],) = nonzero((Ip == i).choose(mask, 1).sum(axis=1) == 0)

    return Id, Jd, Ip, Jp

    def _hessian_main(self, params):
        params_infl = params[:self.k_inflate]
        params_main = params[self.k_inflate:]

        y = self.endog
        w = self.model_infl.predict(params_infl)
        w = np.clip(w, np.finfo(float).eps, 1 - np.finfo(float).eps)
        score = self.score(params)
        zero_idx = np.nonzero(y == 0)[0]
        nonzero_idx = np.nonzero(y)[0]

        mu = self.model_main.predict(params_main)

        hess_arr = np.zeros((self.k_exog, self.k_exog))

        coeff = (1 + w[zero_idx] * (np.exp(mu[zero_idx]) - 1))

        #d2l/dp2
        for i in range(self.k_exog):
            for j in range(i, -1, -1):
                hess_arr[i, j] = ((
                    self.exog[zero_idx, i] * self.exog[zero_idx, j] *
                    mu[zero_idx] * (w[zero_idx] - 1) * (1 / coeff -
                    w[zero_idx] * mu[zero_idx] * np.exp(mu[zero_idx]) /
                    coeff**2)).sum() - (mu[nonzero_idx] * self.exog[nonzero_idx, i] *
                    self.exog[nonzero_idx, j]).sum())

        return hess_arr

def cleaningSineTracks(tfreq, minTrackLength=3):
	"""
	Delete short fragments of a collection of sinusoidal tracks 
	tfreq: frequency of tracks
	minTrackLength: minimum duration of tracks in number of frames
	returns tfreqn: output frequency of tracks
	"""

	if tfreq.shape[1] == 0:                                 # if no tracks return input
		return tfreq
	nFrames = tfreq[:,0].size                               # number of frames
	nTracks = tfreq[0,:].size                               # number of tracks in a frame
	for t in range(nTracks):                                # iterate over all tracks
		trackFreqs = tfreq[:,t]                               # frequencies of one track
		trackBegs = np.nonzero((trackFreqs[:nFrames-1] <= 0)  # begining of track contours
								& (trackFreqs[1:]>0))[0] + 1
		if trackFreqs[0]>0:
			trackBegs = np.insert(trackBegs, 0, 0)
		trackEnds = np.nonzero((trackFreqs[:nFrames-1] > 0)   # end of track contours
								& (trackFreqs[1:] <=0))[0] + 1
		if trackFreqs[nFrames-1]>0:
			trackEnds = np.append(trackEnds, nFrames-1)
		trackLengths = 1 + trackEnds - trackBegs              # lengths of trach contours
		for i,j in zip(trackBegs, trackLengths):              # delete short track contours
			if j <= minTrackLength:
				trackFreqs[i:i+j] = 0
	return tfreq

def truncate_hist1( self, xmin, xmax ):
    buf   = get_buffer_hist1( self )
    sbuf  = get_err_buffer_hist1( self )
    edges, fixed = get_bin_edges_axis( self.GetXaxis(), type=True )

    e1 = numpy.fabs(edges[:-1]-xmin)<1.e-9
    e2 = numpy.fabs(edges[1:]-xmax)<1.e-9
    assert numpy.any( e1 ) and numpy.any( e2 ), 'Invalid new histogram limits'
    i1 = numpy.nonzero( e1 )[0][0]
    i2 = numpy.nonzero( e2 )[0][-1]+1

    if fixed:
        newhist = self.__class__( self.GetName(), self.GetTitle(), i2-i1, xmin, xmax )
    else:
        newhist = self.__class__( self.GetName(), self.GetTitle(), i2-i1, edges[i1:i2] )

    newbuf = get_buffer_hist1( newhist )
    if sbuf is None:
        newsbuf = None
    else:
        newhist.Sumw2()
        newsbuf = get_err_buffer_hist1( newhist )

    newbuf[:] = buf[i1:i2]
    if not sbuf is None:
        newsbuf[:] = sbuf[i1:i2]

    newhist.SetEntries( newhist.Integral() )

    return newhist

def petscKron(A,B):
    dim = A.shape[0]*B.shape[0] # length of resulting matrix
    
    # Used to get indexes where values are non-zero
    Br,Bc = np.nonzero(B)
    Ar,Ac = np.nonzero(A)

    # Need to have values on first axis
    Ar = np.asarray(Ar).ravel(); Ac = np.asarray(Ac).ravel()
    Br = np.asarray(Br).ravel(); Bc = np.asarray(Bc).ravel()

    # Distance between each 'block'
    n = B.shape[1]
    
    # create petsc resulting matrix
    K = PETSc.Mat().createAIJ([dim,dim])
    K.setFromOptions(); K.setUp()
    start,end = K.getOwnershipRange()

    for i in xrange(len(Ar)): # Go through each non-zero value in A
        # br,bc are used to track which 'block' we're in (in result matrix)
        br,bc = n*Ar[i], n*Ac[i]

        for j in xrange(len(Br)): # Go through non-zero values in B
            # kr,kc used to see where to put the number in K (the indexs)
            kr = (Br[j]+br).astype(np.int32)
            kc = (Bc[j]+bc).astype(np.int32)

            if start <= kr < end: # Make sure we're in the correct processor
                K[kr, kc] = A[Ar[i],Ac[i]] * B[Br[j],Bc[j]]

    K.assemble()
    return K

def zpeaki(source,order=1,fpeak=fhigh):
    '''
        寻找n阶高/低点
        返回值为高点数据序列，以及该高点最大跨度的坐标(即计算该高/低点所需用到的最远的未来数据的坐标)
        order默认为1,小于1当作1
        返回值中第一个是高/低点非0,其余为0的序列 sh
                第二个是该高低点的最远未来数据的坐标序列 si
                其中 sh[np.nonzero(sh)]为高点序列, si[np.nonzero(sh)]为坐标序列,sif.time[si[np.nonzero(sh)]]为坐标的影响时间序列
    '''
    tsx1 = fpeak(source)
    sx1 = np.select([tsx1!=0],[source],0)
    icovered = rollx(np.arange(len(source)),-1)
    if order <= 1:
        return sx1,np.select([tsx1],[icovered],0)
    icursx = np.nonzero(tsx1)[0]
    for i in xrange(1,order):   #必然进入循环
        sxx = source[icursx]
        tsxx = fpeak(sxx)
        icovered[icursx] = rollx(icovered[icursx],-1)   #当前高/低点的计算范围,即之前顶点的范围左转一位(排除掉不是顶点的)
        icursx = icursx[np.nonzero(tsxx)[0]]
    osx = np.zeros_like(source)
    osx[icursx] = source[icursx]
    iz = np.zeros_like(source)
    iz[icursx] = icovered[icursx]   #去掉icovered之中不必要的那些数字
    return osx,iz

def aff_cercle_visi(lon, lat, dlon, dlat, col, fig):
    """
    Affichage des cercles de visibilité

    lon, lat : position utilisateur (float)
    dlon, dlat : vecteurs
    col : param affichage des cercles
    """
    from params import CRD

    s = len(dlon)
    lon_vis = np.zeros((2, s))
    lon_vis[0, :] = lon + dlon
    lon_vis[1, :] = lon - dlon
    #if min(lon_vis < 0) or max(lon_vis > 2 * np.pi):

    dlat2 = dlat[(s-1)::-1]  # np.array (250,)

    #lon_vis_array = np.array([lon_vis[0], lon_vis[1], lon_vis[0, 0]]) * CRD
    #dlat_array = (np.array([dlat, dlat2, dlat[0]]) + lat) * CRD
    #plt.plot(lon_vis_array, dlat_array, col)

    latv = dlat + lat

    ## bordures
    ## haut du demi-cercle
    indlat = np.nonzero(latv > np.pi / 2)[0]  # 1 dimension
    latv[indlat] = np.pi - latv[indlat]
    lon_vis[:, indlat] = lon_vis[:, indlat] + np.pi
    #plt.plot(lon_vis[0, indlat]*CRD, latv[indlat]*CRD, 'g-')
    #toto = np.nonzero(latv >= np.pi)[0].size
    #if toto:
    #    print("*"*5, toto)

    ## bas du demi-cercle
    indlat = np.nonzero(latv < -np.pi / 2)[0]  # 1 dimension
    latv[indlat] = -np.pi - latv[indlat]
    lon_vis[:, indlat] = lon_vis[:, indlat] + np.pi
    lon_vis[1, :] = lon_vis[1, (s-1)::-1]
    #plt.plot(lon_vis[1, indlat]*CRD, latv[indlat]*CRD, 'b-')
    #toto = np.nonzero(latv <= -np.pi)[0].size
    #if toto:
    #    print(toto, "*"*5)

    ## côtés du demi-cercle
    lon_vis = lon_vis + 2 * np.pi
    lon_vis = lon_vis % (2 * np.pi)

    latv2 = latv[(s-1)::-1]

    #plt.plot(lon_vis[0, :]*CRD, latv*CRD, col)
    #plt.plot(lon_vis[1, :]*CRD, latv2*CRD, col)
    cercled = lon_vis[0, :]*CRD, latv*CRD
    cercleg = lon_vis[1, :]*CRD, latv2*CRD

    #plt.plot(np.array([lon_vis[0, 0], lon_vis[1, 0]])*CRD,
    #         np.array([latv2[0], latv2[0]])*CRD, col)

    #plt.show()
    return cercleg, cercled

def output_wiggle(bins, binsize, norm_factor, by_strand, name, extra_trackline = ""):
    """write all non-empty bins to bedgraph format strings; always includes
    minimal track line; Output is in 1-based wiggle format."""
    if not by_strand:
        print "track type=wiggle_0 alwaysZero=on visibility=full maxHeightPixels=100:80:50 " \
                + ("name='%s'" % name) + extra_trackline
        for chrom in sorted(bins.keys()):
            print "variableStep chrom=%s span=%d" % (chrom, binsize)
            non_zero_bins = numpy.nonzero(bins[chrom] > 0)
            result = numpy.column_stack((non_zero_bins[0] * binsize + 1,
                bins[chrom][non_zero_bins] * norm_factor))
            numpy.savetxt(sys.stdout, result, "%d\t%.8f")
    else:
        for strand in (0, 1):
            if strand == 0:
                nf = norm_factor
            else:
                nf = -norm_factor
            print "track type=wiggle_0 alwaysZero=on visibility=full maxHeightPixels=100:80:50 " \
                    + ("name='%s[%s]'" % (name, strand and '-' or '+')) + extra_trackline
            for chrom in sorted(bins.keys()):
                print "variableStep chrom=%s span=%d" % (chrom, binsize)
                non_zero_bins = numpy.nonzero(bins[chrom][strand] > 0)
                result = numpy.column_stack((non_zero_bins[0] * binsize + 1,
                    bins[chrom][strand][non_zero_bins] * nf))
                numpy.savetxt(sys.stdout, result, "%d\t%.8f")

def x_input_to_states(xinput, CORR_VAL_OUT=0, PARALLEL = False):
    sinput = np.zeros(xinput.shape)
    num_samples = xinput.shape[0]
    num_sensors = xinput.shape[1]
    if num_samples < num_sensors:
        print '[WARN] number of samplesa are smaller than number of sensors'

    print 'Mapping', xinput.shape, ' marix to discrete states '

    for k, samples in enumerate(xinput.T):
        obs = samples[:,np.newaxis]
        label, opt_num_cluster, model, score, score_err_sum = state_retrieval(obs, max_num_cluster=6, est_method='kmean', PARALLEL=PARALLEL)
        high_peak_label_idx = np.argmax(model.cluster_centers_)
        low_peak_label_idx = np.argmin(model.cluster_centers_)
        high_peak_idx = np.nonzero(label == high_peak_label_idx)[0]
        sinput[high_peak_idx,k] = 1
        low_peak_idx = np.nonzero(label == low_peak_label_idx)[0]
        sinput[low_peak_idx, k] = -1

    corr_state_val = []
    if CORR_VAL_OUT == 1:
        print 'Compute Correlation Score....'
        for k,(row1, row2) in enumerate(zip(sinput.T, xinput.T)):
            corr_state_val.append(round(stats.pearsonr(row1, row2)[0],3))
    corr_state_val = np.array(corr_state_val)
    return sinput, corr_state_val

def main():
	rfile = sys.argv[1]

	csvfile = open(rfile, 'rb')
	dat = csv.reader(csvfile, delimiter=',')

	X  = []
	Y = []

	for i, row in enumerate(dat):
		if i > 0:
			X.append(float(row[0]))
			Y.append(int(row[1]))

	X = np.array(X)
	Y = np.array(Y)

	class1 = np.array(X[np.nonzero(Y == 1)[0]])
	class2 = np.array(X[np.nonzero(Y == 2)[0]])

	print("computing...")
	#build GMM for two classes
	model1 = build_models1(class1)
	#model2 = build_models2(class2)

	print("Here are the models!")

	print(model1[0])
	#print(model2[0])


	plt.plot(range(1,20+1), model1[1], 'ro')
	#plt.plot(range(1,20+1), model2[1], 'ro')
	plt.show() 

def CCML(ccpart,ccmat,mu,r,nu,s,row_alpha,col_alpha):
	lp = []
	
	ccmat = np.array(ccmat)

	state = ccpart.states[1]

	# loop through the states
	for state in ccpart.states:
		all_cols = state['col_parts']
		all_rows = state['row_parts']

		K = max(all_cols)+1
				
		lp_temp = lcrp(all_cols,col_alpha)
		for view in range(K):
			row_part = all_rows[view,:]
			lp_temp += lcrp(row_part,row_alpha)
			cols_view = np.nonzero(all_cols==view)[0]
			for col in cols_view:
				for cat in range(row_part.max()+1):
					X = ccmat[np.nonzero(row_part==cat)[0],col]
					lp_temp += NGML(X,mu,r,nu,s)

		lp.append(lp_temp);

	# return the normalized probabilities
	return lp-sp.misc.logsumexp(lp)

def dup_idx(arr):
    """
    Return the indices of all duplicated array elements.

    Parameters
    ----------

    arr : array-like object

    Returns
    -------

    idx : NumPy array
        An array containing the indices of the duplicated elements

    Examples
    --------

    >>> from root_numpy import dup_idx
    >>> dup_idx([1, 2, 3, 4, 5])
    array([], dtype=int64)
    >>> dup_idx([1, 2, 3, 4, 5, 5])
    array([4, 5])
    >>> dup_idx([1, 2, 3, 4, 5, 5, 1])
    array([0, 4, 5, 6])

    """
    _, b = np.unique(arr, return_inverse=True)
    return np.nonzero(np.logical_or.reduce(
        b[:, np.newaxis] == np.nonzero(np.bincount(b) > 1),
        axis=1))[0]