Python scipy.nonzero函数代码示例

    def getGenoIndex(self,pos_start=None,pos_end=None,chrom=None,windowsize=0):
        """computes 0-based genotype index from position of cumulative position.
        Positions can be given in one out of two ways:
        - position (pos_start-pos_end on chrom)
        - cumulative position (pos_cum_start-pos_cum_end)
        If all these are None (default), then all genotypes are returned

        Args:
            pos_start:       position based selection (start position)
            pos_end:       position based selection (end position)
            chrom:      position based selection (chromosome)

        Returns:
            idx_start:         genotype index based selection (start index)
            idx_end:         genotype index based selection (end index)
        """
        if (pos_start is not None) & (pos_end is not None):
            assert pos_start[0] == pos_end[0], "chromosomes don't match between start and end position"
            I = self.position["chrom"]==pos_start[0]
            I = I & (self.position["pos"]>=(pos_start[1]-windowsize)) & (self.position["pos"]<(pos_end[1]+windowsize))
            I =  sp.nonzero(I)[0]
            idx_start = I.min()
            idx_end = I.max()
        elif chrom is not None:
            I = self.position["chrom"]==chrom
            I =  sp.nonzero(I)[0]
            if I.size==0:
                return None
            idx_start = I.min()
            idx_end = I.max()
        else:
            idx_start=None
            idx_end=None
        return idx_start,idx_end

def graph_connected_components(graph_mat):
    """takes graph as matrix and return list of connected components
          
       arguments:
       graph_mat - matrix of graph
              
       output:
       list_of_comp - list of components, list_of_comp[i] - list of node numbers in (i-1)-th component"""
    component_of_graph = scipy.zeros(graph_mat.shape[0],dtype = scipy.int8) # component_of_graph[i] is the number of component of i-th node.
    cur_comp = 1 #the number of current component (see below)
    try:
        tmp_nodes_to_process = [scipy.nonzero(component_of_graph==0)[0][0]] #node indexes to process
    except IndexError:
            #exceptional situation, when graph_mat is empty
            return []

    #kind of breadth first search
    while(len(tmp_nodes_to_process)>0): #while there is nodes to process
        cur_node = tmp_nodes_to_process.pop() #take one from array
        lnodes_numbers = scipy.nonzero(graph_mat[cur_node,:])[0] #take indexes of all of it linked nodes
        #and choose that corresponds to the non processed nodes of them, the node is non processed if its component is zero
        lnodes_numbers = scipy.extract(component_of_graph[lnodes_numbers] == 0,lnodes_numbers)
        tmp_nodes_to_process +=lnodes_numbers.tolist()
        component_of_graph[lnodes_numbers] = cur_comp
        # if there is no linked nodes, start processing of new connected component, and next unprocessed node
        if (len(tmp_nodes_to_process) == 0):
            cur_comp+=1
            tmp_arr = scipy.nonzero(component_of_graph==0)[0]
            if (len(tmp_arr)>0):tmp_nodes_to_process = [tmp_arr[0]] 

    list_of_comp = [] #collect list
    for i in range(cur_comp+1):
        tmp_arr=scipy.nonzero(component_of_graph==(i+1))[0].tolist()
        if (len(tmp_arr)>0):list_of_comp+=[tmp_arr]
    return list_of_comp

 def cut(self):
   average = sp.sum(sp.absolute(self.data))/sp.size(self.data)
   head = sp.nonzero(sp.absolute(self.data)>average)[0][5]
   bottom = sp.nonzero(sp.absolute(self.data)>average)[0][-1]
   self.data = self.data[head:bottom]
   self.duration_list = self.duration_list[head:bottom]
   self.duration = self.duration_list[-1] - self.duration_list[0]

 def setup_Q(self):
     self.Q = scipy.zeros((self.nperiods, self.ndists, self.ndists))
     #D = self.D * self.Dmask
     for p in range(self.nperiods):
         for i, d1 in enumerate(self.dists):
             s1 = sum(d1)
             if s1 > 0:
                 for j, d2 in enumerate(self.dists):
                     s2 = sum(d2)
                     xor = scipy.logical_xor(d1, d2)
                     # only consider transitions between dists that are
                     # 1 step (dispersal or extinction) apart
                     if sum(xor) == 1:
                         dest = scipy.nonzero(xor)[0]
                         #prior = self.dist_priors[i]
                         if s1 < s2: # dispersal
                             rate = 0.0
                             for src in scipy.nonzero(d1)[0]:
                                 rate += self.D[p,src,dest] * \
                                         self.Dmask[p,src,dest]
                             # for each area in d1, add rate of
                             # dispersal to dest
                         else: # extinction
                             rate = self.E[p,dest]
                         #self.Q[i][j] = (prior * rate)
                         self.Q[p,i,j] = rate
         self.set_Qdiag(p)

    def getGenoIndex(self,pos0=None,pos1=None,chrom=None,pos_cum0=None,pos_cum1=None):
        """computes 0-based genotype index from position of cumulative position. 
        Positions can be given in one out of two ways: 
        - position (pos0-pos1 on chrom)
        - cumulative position (pos_cum0-pos_cum1)
        If all these are None (default), then all genotypes are returned

        Args:
            pos0:       position based selection (start position)
            pos1:       position based selection (stop position)
            chrom:      position based selection (chromosome)
            pos_cum0:   cumulative position based selection (start position)
            pos_cum1:   cumulative position based selection (stop position)
        
        Returns:
            i0:         genotype index based selection (start index)
            i1:         genotype index based selection (stop index)
        """
        if (pos0 is not None) & (pos1 is not None) & (chrom is not None):
            I = self.gneoChrom==chrom
            I = I & (self.genoPos>=p0) & (self.genoPos<p1)
            I = SP.nonzero(I)[0]
            i0 = I.min()
            i1 = I.max()
        elif (pos_cum0 is not None) & (pos_cum1 is not None):
            I = (self.genoPos_cum>=pos_cum0) & (self.genoPos_cum<pos_cum1)
            I = SP.nonzero(I)[0]
            if I.size==0:
                return None
            i0 = I.min()
            i1 = I.max()
        else:
            i0=None
            i1=None
        return i0,i1

def __interpolateBetweenBinaryObjects(obj1, obj2, slices):
    """
    Takes two binary objects and puts slices slices in-between them, each of which
    contains a smooth binary transition between the objects.
    @note private inner function
    """
    if not obj1.shape == obj2.shape:
        raise AttributeError(
            "The two supplied objects have to be of the same shape, not {} and {}.".format(obj1.shape, obj2.shape)
        )

    # constant
    offset = 0.5  # must be a value smaller than the minimal distance possible
    temporal_dimension = 3

    # get all voxel position
    obj1_voxel = scipy.nonzero(obj1)
    obj2_voxel = scipy.nonzero(obj2)

    # get smallest pairwise distances between all object voxels
    distances = cdist(scipy.transpose(obj1_voxel), scipy.transpose(obj2_voxel))

    # keep for each True voxel of obj1 only the smallest distance to a True voxel in obj2
    min_distances = distances.min(1)

    # test if all seems to work
    if len(min_distances) != len(obj1_voxel[0]):
        raise Exception("Invalid number of minimal distances received.")

    # replace True voxels in obj1 with their respective distances to the True voxels in obj2
    thr_obj = obj1.copy()
    thr_obj = thr_obj.astype(scipy.float_)
    thr_obj[obj1_voxel] = min_distances
    thr_obj[obj1_voxel] += offset  # previous steps distances include zeros, therefore this is required

    # compute the step size for each slice that is added
    maximum = min_distances.max()
    step = maximum / float(slices + 1)
    threshold = maximum

    # control step: see if thr_obj really corresponds to obj1
    if not scipy.all(thr_obj.astype(scipy.bool_) == obj1.astype(scipy.bool_)):
        raise Exception("First created object does not correspond to obj1.")

    # assemble return volume
    return_volume = [thr_obj.astype(scipy.bool_)]  # corresponds to obj1
    for _ in range(slices):
        threshold -= step
        # remove all value higher than the threshold
        thr_obj[thr_obj > threshold] = 0
        # add binary volume to list /makes a copy)
        return_volume.append(thr_obj.astype(scipy.bool_))

    # add last slice (corresponds to es obj2 slice)
    thr_obj[thr_obj > offset] = 0
    return_volume.append(thr_obj.astype(scipy.bool_))

    # return binary scipy array
    return scipy.rollaxis(scipy.asarray(return_volume, dtype=scipy.bool_), 0, temporal_dimension + 1)

def move_nodes(d, lfun, nodes):#, dely):
    " Move nodes one timestep (projecting lost points to the boundary) "

    # get delauney 
    dely = get_dely(d, nodes)

    # force constant 
    dt = .1520015
    deps = h0 * sqrt(finfo(float64).eps)
    restlength_factor = 1.400025
    
    # bars and midpoints
    bars, barmids = dely[-2:]
    barvecs = nodes[:,bars[:,1]] - nodes[:,bars[:,0]]
    barls = sqrt((barvecs**2).sum(0))
    u = barvecs / barls       # unit vectors
    
    # force from each bar
    restlen = restlength_factor * lfun(*barmids)

    # print('restlen: {0} \nbarls : {1}'.format(restlen.shape, barls.shape))
    logic = restlen > barls
    f = where(logic, restlen-barls, 0.0)
    # f = where(f<h0/2.0, f, h0/2.0)
    # 
    ns = nodes.shape
    spmat = sparse.csc_matrix
    # print(ns)
    # print(u[0].shape)
    # print(f.shape)
    # print(bars[:,0].shape)

    dp =  (-spmat((u[0]*f, (0*bars[:,0], bars[:,0])), shape=ns).todense() +
            spmat((u[0]*f, (0*bars[:,1], bars[:,1])), shape=ns).todense())
    dp += (-spmat((u[1]*f, (0*bars[:,0]+1, bars[:,0])), shape=ns).todense() +
            spmat((u[1]*f, (0*bars[:,1]+1, bars[:,1])), shape=ns).todense())

    nodes = array(nodes + dt*dp)

    # project boundary points back into the domain
    d_ = d(*nodes)
    ix = nonzero(d_>0)
    some_out = True
    count = 0
    while some_out:
        gradx = 1.0/deps * (d(nodes[0,ix]+deps, nodes[1,ix]) - d_[ix])
        grady = 1.0/deps * (d(nodes[0,ix], nodes[1,ix] + deps) - d_[ix])
        norm = sqrt(gradx**2 + grady**2)
        nodes[0,ix] -= d_[ix]*gradx / norm
        nodes[1,ix] -= d_[ix]*grady / norm
        d_ = d(*nodes)
        ix = nonzero(d_>geps)
        some_out = ix[0].size
        count+=1
        if count>5: #raise ValueError("counted "+str(ix[0].size)+" nodes oob")
            print("counted ",str(ix[0].size)," nodes oob")
            break
    return nodes

 def __call__(self, dx):
     """
     hat function evaluation
     """
     x = convarg(dx) - self.c
     y = sc.zeros_like(x)
     i = sc.nonzero((x>self.a)&(x<=.0))
     y[i] = 1. - x[i]/self.a
     i = sc.nonzero((x>.0)&(x<self.b))
     y[i] = 1. - x[i]/self.b
     return(y)

def gen_single_trial(interval_lengths, rates):
    """ Generates a single spike train with intervals of length
    `interval_lengths` and the firing rates given in `rates`.
    """
    boundaries = sp.ones(len(interval_lengths) + 1) * pq.s
    boundaries[1:] = [l.rescale(boundaries.units) for l in interval_lengths]
    rates = rates[sp.nonzero(boundaries[1:])]
    boundaries = boundaries[sp.nonzero(boundaries)]
    boundaries[0] = 0.0 * pq.s
    boundaries = sp.cumsum(boundaries)
    return stools.st_concatenate([stg.gen_homogeneous_poisson(
        rate, t_start=boundaries[i], t_stop=boundaries[i + 1])
        for i, rate in enumerate(rates)])

def subgraph_with_center(graph_mat, node_names, center_name):
    """takes number of center node, and returns new subgraph, that consists of all nodes connected with center.
        
        arguments:
        graph_mat,node_names - graph description as matrix (see matrix_from_tuple_list doc for details)
        center_num - name of the center node

        output:
        subgraph,sub_node_names - subgraph description as matrix (see matrix_from_tuple_list doc for details)"""
    center_num = scipy.nonzero(node_names==center_name)[0][0]
    center_friends_num = scipy.nonzero(graph_mat[center_num,:])[0] #indexes of nodes that linked with  central node including itself
    subgraph = graph_mat[center_friends_num,:][:,center_friends_num] # FIXME we consider part of graph which consists of nodes linked with center only 
    sub_node_names = node_names[center_friends_num]
    return (subgraph,sub_node_names)

def dup_fig(x, y):
    colors = [None, 'black','blue','brown','purple','orange','cyan','gray','yellow','black','red','green']
    k = len(sp.unique(y))
    for i in range(1, k+1):
        inds = sp.nonzero(y == i)
        plt.plot(x[inds, 0], x[inds, 1], 'o', color=colors[i])
    plt.show()

    def getPhenotypes(self,phenotype_IDs=None,phenotype_query=None,sample_idx=None,center=True,intersection=False):
        """load Phenotypes

        Args:
            phenotype_IDs:      names of phenotypes to load
            phenotype_query:    string hoding a pandas query (e.g. "(environment==1) & (phenotype_ID=='growth')"
                                selects all phenotypes that have a phenotype_ID equal to growth under environment 1.
            sample_idx:         Boolean sample index for subsetting individuals
            center:             Boolean: mean center (and mean-fill in missing values if intersection==False)? (default True)
            impute:             imputation of missing values (default: True)
            intersection:       restrict observation to those obseved in all phenotypes? (default: False)

        Returns:
            phenotypes:     [N x P] scipy.array of phenotype values for P phenotypes
            sample_idx_intersect:        index of individuals in phenotypes after filtering missing valuesS
        """
        if phenotype_IDs is not None:
            I = SP.array([SP.nonzero(self.phenotype_ID==n)[0][0] for n in phenotype_IDs])
        elif phenotype_query is not None:
            try:
                I = self.index_frame.query(phenotype_query).values[:,0]

                #if there are no results we won't actually get an exception, we just get an
                #empty response
                if len(I) == 0:
                    print "query '%s' yielded no results!" % (phenotype_query)
                    I = SP.zeros([0],dtype="int")

            except Exception, arg:

                print "query '%s' yielded no results: %s" % (phenotype_query, str(arg))

                I = SP.zeros([0],dtype="int")

	def __init__(self, grid, fArray, zDrawsSorted):
		assert(len(grid) == len(fArray))
		(self.grid, self.fArray) = (grid, fArray)
		self.zDraws = zDraws		
		self.slopes = scipy.zeros(len(grid) - 1)
		self.dx = grid[1] - grid[0]
		for i in range(len(grid) - 1):
			self.slopes[i] = (fArray[i+1] - fArray[i]) / self.dx
		# set up sums
		self.cellSums = scipy.zeros(len(grid) + 1)
		self.boundaryIndices = [len(zDraws)] * len(grid)
		for (i, x) in enumerate(grid):
			indices = scipy.nonzero(self.zDraws >= x)[0]
			if (len(indices) > 0):
				self.boundaryIndices[i] = indices[0]
		self.cellSums[0] = scipy.sum(self.zDraws[0:self.boundaryIndices[0]])
		for i in range(1, len(self.cellSums)-1):
			self.cellSums[i] = scipy.sum(self.zDraws[self.boundaryIndices[i-1] : self.boundaryIndices[i]])
		self.cellSums[-1] = scipy.sum(self.zDraws[self.boundaryIndices[-1] : ])
		
		diff = scipy.sum(self.zDraws) - scipy.sum(self.cellSums)
		print("diff: %f" % diff)
		for i in range(len(grid)):
			if (self.boundaryIndices[i] < len(self.zDraws)):
				print("grid point %f, boundary %f" % (self.grid[i], self.zDraws[self.boundaryIndices[i]]))
			else:
				print("grid point %f, no draws to right" % self.grid[i])

    def errorApproximation(self, ratio, dim=20):

        self.buildMatrix()

        sumNonzeros = (self.vxm !=0).sum()
        numTest = int(ratio*sumNonzeros)

        elementList = []

        nonZeroTuple = sp.nonzero(self.vxm)

        for x in range(int(numTest)):
            rInt = sp.random.randint(0,nonZeroTuple[0].size)
            randrow = nonZeroTuple[0][rInt]
            randcolumn = nonZeroTuple[1][rInt]

            valElementIndex = [randrow,randcolumn]
            elementList.append(valElementIndex)

        self.modvxm = sp.copy(self.vxm)

        for x in elementList:
            self.modvxm[x[0],x[1]] = 0

        self.modvmx = self.fillAverages(vxm = self.modvxm)
        self.newmodvxm = self.predict(dim,vxm=self.modvxm)

        sqDiff = 0
        for x in elementList:
            sqDiff += sp.square(self.newmodvxm[x[0],x[1]] - self.vxm[x[0],x[1]])
        self.rmse = sp.sqrt(sqDiff/len(elementList))

    def _generate_pores(self):
        r"""
        Generate the pores (coordinates, numbering and types)
        """
        self._logger.info("generate_pores: Create specified number of pores")

        #Find non-zero elements in image
        template = self._template
        Np = np.sum(template > 0)
        #Add pores to data and ifo
        pind = np.arange(0, Np)
        self.set_pore_info(label='all', locations=pind)
        self.set_pore_data(prop='numbering', data=pind)  # Remove eventually

        
        img_ind = np.ravel_multi_index(sp.nonzero(template), dims=sp.shape(template), order='F')
        self.set_pore_data(prop='voxel_index', data=img_ind)

        #This voxel_to_pore map is messy but works
        temp = sp.prod(sp.shape(template))*sp.ones(np.prod(sp.shape(template),),dtype=sp.int32)
        temp[img_ind] = pind
        self._voxel_to_pore_map = temp

        coords = self._Lc*(0.5 + np.transpose(np.nonzero(template)))
        self.set_pore_data(prop='coords', data=coords)
        self._logger.debug("generate_pores: End of method")