Python numpy.nonzero函数代码示例

    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