Python scipy.unique函数代码示例

def filter_introns(introns, genes, options):
    
    ### build interval trees of all genes starts and ends
    chrms = sp.array([_.strand for _ in genes])
    strands = sp.array([_.chr for _ in genes])
    gene_trees = dict()
    for c in sp.unique(chrms):
        for s in sp.unique(strands):
            gene_trees[(c, s)] = it.IntervalTree()
            c_idx = sp.where((chrms == c) & (strands == s))[0]
            for i in c_idx:
                gene_trees[(c, s)][genes[i].start:genes[i].stop] = i

    ### match all introns agains trees and remove elements overlapping
    ### more than one gene on the same chr/strand
    cnt_tot = 0
    cnt_rem = 0
    strand_list = ['+', '-']
    offset = options.intron_edges['append_new_terminal_exons_len']
    for si, s in enumerate(strand_list):
        for i in range(introns.shape[0]):
            if introns[i, si].shape[0] == 0:
                continue
            k_idx = []
            cnt_tot += introns[i, si].shape[0]
            for j in range(introns[i, si].shape[0]):
                if len(gene_trees[(s, genes[i].chr)].overlap(introns[i, si][j, 0] - offset, introns[i, si][j, 1] + offset)) == 1:
                    k_idx.append(j)
            if len(k_idx) < introns[i, si].shape[0]:
                cnt_rem += (introns[i, si].shape[0] - len(k_idx))
                introns[i, si] = introns[i, si][k_idx, :]
    print('removed %i of %i (%.2f percent) introns overlapping to no or multiple genes' % (cnt_rem, cnt_tot, cnt_rem / float(max(cnt_tot, 1)) * 100))

    return introns

def main():
    args = getArguments(getParser())

    # prepare logger
    logger = Logger.getInstance()
    if args.debug:
        logger.setLevel(logging.DEBUG)
    elif args.verbose:
        logger.setLevel(logging.INFO)

    # check if output image exists
    if not args.force:
        if os.path.exists(args.output):
            logger.warning("The output image {} already exists. Exiting.".format(args.output))
            exit(-1)

    # load input image
    input_data, input_header = load(args.input)

    logger.debug("Old number of regions={}.".format(len(scipy.unique(input_data))))

    # cut and relabel along the required dimension
    logger.info("Cutting and relabeling...")
    dimensions = range(input_data.ndim)
    del dimensions[args.dimension]
    __split_along(input_data, dimensions)

    logger.debug("New number of regions={}.".format(len(scipy.unique(input_data))))

    # save result contour volume
    save(input_data, args.output, input_header, args.force)

    logger.info("Successfully terminated.")

 def _read_sky_logfile(self):
     #TODO : expand to read errors, msgs etc
     # read in the whole sky log file, shouldn't be big
     f = open(self.skylogfile)
     lines = f.readlines()
     f.close()
     dust = [line.split()[1:] for line in lines if line.startswith('dtau_dust')]
     line = [line.split()[1:] for line in lines if line.startswith('dtau_line')]
     dust = _sp.array(dust, dtype='float')
     line = _sp.array(line, dtype='float')
     transitions = _sp.unique(dust[:,0])
     shells = _sp.unique(dust[:,1])
     dtau_dust = dict()
     dtau_line = dict()
     dtau_tot = dict()
     for t in transitions:
         d = []
         l = []
         for s in shells:
             d.append( _sp.mean([i[2] for i in dust if ((i[0]==t) * (i[1]==s))]) )
             l.append( _sp.mean([i[2] for i in line if ((i[0]==t) * (i[1]==s))]) )
         dtau_dust[t] = _sp.copy(d)
         dtau_line[t] = _sp.copy(l)
         dtau_tot[t] = _sp.array(d) + _sp.array(l)
     # create object to store in main class
     class Tau(object):pass
     Tau.dtau_dust = dtau_dust
     Tau.dtau_line = dtau_line
     Tau.dtau_tot = dtau_tot
     Tau.transitions = transitions
     Tau.shells = shells
     self.Tau = Tau

def Capillary_Pressure_Curve(net,
                             fluid,
                             capillary_pressure='capillary_pressure',
                             pore_volume='volume',
                             throat_volume='volume',
                             fig=None):
  r"""
  Plot drainage capillary pressure curve

  Parameters
  ----------
  net : OpenPNM Network Object
      The network for which the graphs are desired
  fig : Matplotlib figure object
      Canvas on which to draw plots

  """
  if type(fluid)==str: fluid = net.find_object_by_name(fluid)
  try:
    PcPoints = sp.unique(fluid.get_throat_data(prop=capillary_pressure))
  except KeyError:
    raise Exception('Capillary pressure simulation has not been run')

  PcPoints = sp.unique(fluid.get_throat_data(prop=capillary_pressure))
  Snwp = sp.zeros_like(PcPoints)
  Ps = sp.r_[0:net.num_pores('internal')]
  for i in range(1,sp.size(PcPoints)):
      Pc = PcPoints[i]
      Snwp[i] = sum((fluid.get_throat_data(prop=capillary_pressure)[Ps]<Pc)*(net.get_throat_data(prop=throat_volume)[Ps]))/sum(net.get_throat_data(prop=throat_volume)[Ps])
  
  if fig==None: fig = plt.figure()
  ax = fig.add_subplot(111)
  ax.plot(PcPoints,Snwp,'r.-')
  ax.set_xlabel('Capillary Pressure')
  ax.set_ylabel('Fluid Saturation')

def batch_metrics(unit_list, threshold, t_ref, t_cen):
    ''' This here function runs metrics on a batch of data.  Pass in units from the catalog.
    '''
    
    from scipy import unique
    
    samp_rate = 30000.
    n_samples = 30
    n_chans = 4
    
    # Find common Sessions
    sessions = unique([unit.session for unit in unit_list])
    
    for session in sessions:
        
        units = session.units
        tetrodes = unique([unit.tetrode for unit in units])
        
        for tetrode in tetrodes:
            data = load_spikes(session.path, tetrode,  samp_rate, n_samples, n_chans)
            f_p, f_n = metrics(data, threshold, t_ref, t_cen, session.duration)
            # Doing this because sometimes there is no cluster 0 sometimes
            f_p.setdefault(1)
            f_n.setdefault(1)
            units = [ unit for unit in session.units if unit.tetrode == tetrode] 
            for unit in units:
                unit.falsePositive = f_p[unit.cluster]
                unit.falseNegative = f_n[unit.cluster]

    def test_set_boundary_conditions_bctypes(self):
        self.alg.setup(invading_phase=self.water,
                       defending_phase=self.air,
                       trapping=True)
        Ps = sp.random.randint(0, self.net.Np, 10)

        self.alg.set_boundary_conditions(pores=Ps, bc_type='inlets')
        assert sp.sum(self.alg['pore.inlets']) == sp.size(sp.unique(Ps))
        self.alg['pore.inlets'] = False

        self.alg.set_boundary_conditions(pores=Ps, bc_type='outlets')
        assert sp.sum(self.alg['pore.outlets']) == sp.size(sp.unique(Ps))
        self.alg['pore.outlets'] = False

        self.alg.set_boundary_conditions(pores=Ps, bc_type='residual')
        assert sp.sum(self.alg['pore.residual']) == sp.size(sp.unique(Ps))
        self.alg['pore.residual'] = False

        flag = False
        try:
            self.alg.set_boundary_conditions(pores=Ps, bc_type='bad_type')
        except:
            flag = True
        assert flag

        flag = False
        try:
            self.alg.set_boundary_conditions(bc_type=None, mode='bad_type')
        except:
            flag = True
        assert flag

    def _check_trapping(self, inv_val):
        r"""
        Determine which pores and throats are trapped by invading phase.  This
        method is called by ``run`` if 'trapping' is set to True.
        """
        # Generate a list containing boolean values for throat state
        Tinvaded = self['throat.inv_Pc'] < sp.inf
        # Add residual throats, if any, to list of invaded throats
        Tinvaded = Tinvaded + self['throat.residual']
        # Invert logic to find defending throats
        Tdefended = ~Tinvaded
        [pclusters, tclusters] = self._net.find_clusters2(mask=Tdefended,
                                                          t_labels=True)
        # See which outlet pores remain uninvaded
        outlets = self['pore.outlets']*(self['pore.inv_Pc'] == sp.inf)
        # Identify clusters connected to remaining outlet sites
        def_clusters = sp.unique(pclusters[outlets])
        temp = sp.in1d(sp.unique(pclusters), def_clusters, invert=True)
        trapped_clusters = sp.unique(pclusters)[temp]
        trapped_clusters = trapped_clusters[trapped_clusters >= 0]

        # Find defending clusters NOT connected to the outlet pores
        pmask = np.in1d(pclusters, trapped_clusters)
        # Store current applied pressure in newly trapped pores
        pinds = (self['pore.trapped'] == sp.inf) * (pmask)
        self['pore.trapped'][pinds] = inv_val

        # Find throats on the trapped defending clusters
        tinds = self._net.find_neighbor_throats(pores=pinds,
                                                mode='intersection')
        self['throat.trapped'][tinds] = inv_val
        self['throat.entry_pressure'][tinds] = 1000000

def test_linear_solvers():
    pn = OpenPNM.Network.Cubic([1, 40, 30], spacing=0.0001)
    geom = OpenPNM.Geometry.Toray090(network=pn,
                                     pores=pn.pores(),
                                     throats=pn.throats())
    air = OpenPNM.Phases.Air(network=pn)
    phys_air = OpenPNM.Physics.Standard(network=pn,
                                        phase=air,
                                        pores=pn.pores(),
                                        throats=pn.throats())

    BC1_pores = pn.pores(labels=['left'])
    BC2_pores = pn.pores(labels=['right'])

    alg_1 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_1.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=1,
                                  pores=BC1_pores)
    alg_1.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_1.run(iterative_solver='gmres')

    alg_2 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_2.set_boundary_conditions(bctype='Neumann',
                                  bcvalue=-1e-11,
                                  pores=BC1_pores)
    alg_2.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_2.run(iterative_solver='cg')

    alg_3 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_3.set_boundary_conditions(bctype='Neumann_group',
                                  bcvalue=-3e-10,
                                  pores=BC1_pores)
    alg_3.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_3.run()

    alg_4 = OpenPNM.Algorithms.FickianDiffusion(network=pn, phase=air)
    alg_4.set_boundary_conditions(bctype='Neumann_group',
                                  bcvalue=-3e-10,
                                  pores=BC1_pores)
    alg_4.set_boundary_conditions(bctype='Dirichlet',
                                  bcvalue=0,
                                  pores=BC2_pores)
    alg_4.setup()
    alg_4.solve()

    assert round(sp.absolute(alg_1.rate(BC1_pores))[0], 16) == round(sp.absolute(alg_1.rate(BC2_pores))[0], 16)
    assert round(sp.absolute(alg_2.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_2['pore.'+air.name+'_bcval_Neumann']))[0]*len(BC1_pores), 16)
    assert round(sp.absolute(alg_3.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_3['pore.'+air.name+'_bcval_Neumann_group']))[0], 16)
    assert round(sp.absolute(alg_4.rate(BC2_pores))[0], 16) == round(sp.absolute(sp.unique(alg_4['pore.'+air.name+'_bcval_Neumann_group']))[0], 16)

    assert round(sp.absolute(sp.sum(alg_1.rate(BC1_pores,mode='single'))),16) == round(sp.absolute(alg_1.rate(BC1_pores))[0],16)
    assert round(sp.absolute(sp.sum(alg_2.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_2.rate(BC2_pores))[0],16)
    assert round(sp.absolute(sp.sum(alg_3.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_3.rate(BC2_pores))[0],16)
    assert round(sp.absolute(sp.sum(alg_4.rate(BC2_pores,mode='single'))),16) == round(sp.absolute(alg_4.rate(BC2_pores))[0],16)

 def _do_outer_iteration_stage(self):
     #Generate curve from points
     for inv_val in self._inv_points:
         #Apply one applied pressure and determine invaded pores
         logger.info('Applying capillary pressure: '+str(inv_val))
         self._do_one_inner_iteration(inv_val)
     #Store results using networks' get/set method
     self['pore.inv_Pc'] = self._p_inv
     self['throat.inv_Pc'] = self._t_inv
     #Find invasion sequence values (to correspond with IP algorithm)
     self._p_seq = sp.searchsorted(sp.unique(self._p_inv),self._p_inv)
     self._t_seq = sp.searchsorted(sp.unique(self._t_inv),self._t_inv)
     self['pore.inv_seq'] = self._p_seq
     self['throat.inv_seq'] = self._t_seq
     #Calculate Saturations
     v_total = sp.sum(self._net['pore.volume'])+sp.sum(self._net['throat.volume'])
     sat = 0.
     self['pore.inv_sat'] = 1.
     self['throat.inv_sat'] = 1.
     for i in range(self._npts):
         inv_pores = sp.where(self._p_seq==i)[0]
         inv_throats = sp.where(self._t_seq==i)[0]
         new_sat = (sum(self._net['pore.volume'][inv_pores])+sum(self._net['throat.volume'][inv_throats]))/v_total
         sat += new_sat
         self['pore.inv_sat'][inv_pores] = sat
         self['throat.inv_sat'][inv_throats] = sat

def __compute_affiliation(label_image, mask_image, bounding_boxes):
    """
    Computes which regions of the supplied label_image belong to the mask_image's foreground
    respectively background. When a region belongs to both, it is assigned to the foreground
    if more voxels belong to the foreground than in the background and vice-versa.
    In the case of equal affiliation, the region is assigned to the background.
    @return fg_ids, bg_ids
    """
    # simple extraction
    fg_ids = list(scipy.unique(label_image[mask_image]))
    bg_ids = list(scipy.unique(label_image[~mask_image]))
    # decide for overlapping regions whether they are 50 or more in fg or in bg
    for rid in set(fg_ids) & set(bg_ids):
        relevant_region_label_image = label_image[bounding_boxes[rid - 1]]
        relevant_region_mask_image = mask_image[bounding_boxes[rid - 1]]
        fg_part = 0
        bg_part = 0
        for affiliation, rid2 in zip(relevant_region_mask_image.ravel(), relevant_region_label_image.ravel()):
            if rid2 == rid:
                if affiliation: fg_part += 1
                else: bg_part += 1
        #fg_part = relevant_region_label_image[relevant_region_mask_image]
        #bg_part = relevant_region_label_image[~relevant_region_mask_image]
        if fg_part > bg_part: # if more voxels of region rid in fg than in bg
            bg_ids.remove(rid)
        else:
            fg_ids.remove(rid)
    # debug line, can be removed if the above code is final
    if 0 != len(set(fg_ids) & set(bg_ids)): raise Exception('Error making fg and bg ground truth distinct.') 
    return fg_ids, bg_ids

def fit_dispersion(counts, disp_raw, disp_conv, sf, CFG, dmatrix1):

    mean_count = sp.mean(counts / sf, axis=1)[:, sp.newaxis]
    index = sp.where(disp_conv)[0]

    lowerBound = sp.percentile(sp.unique(disp_raw[index]), 1)
    upperBound = sp.percentile(sp.unique(disp_raw[index]), 99)

    idx = sp.where((disp_raw > lowerBound) & (disp_raw < upperBound))[0]

    matrix = sp.ones((idx.shape[0], 2), dtype='float')
    matrix[:, 0] /= mean_count[idx].ravel()

    modGamma = sm.GLM(disp_raw[idx], matrix, family=sm.families.Gamma(sm.families.links.identity))
    res = modGamma.fit()
    Lambda = res.params

    disp_fitted = disp_raw.copy()
    ok_idx = sp.where(~sp.isnan(disp_fitted))[0]
    disp_fitted[ok_idx] = Lambda[0] / mean_count[ok_idx] + Lambda[1]

    if sp.sum(disp_fitted > 0) > 0:
        print "Found dispersion fit"

    if CFG['diagnose_plots']:
        plot.mean_variance_plot(counts=counts,
                                disp=disp_fitted,
                                matrix=dmatrix1,
                                figtitle='Fitted Dispersion Estimate',
                                filename=os.path.join(CFG['plot_dir'], 'dispersion_fitted.pdf'),
                                CFG=CFG)

    return (disp_fitted, Lambda, idx)

 def plot_setup(p):
     pylab.ylabel("Throughput gain [\%]")
     pylab.xscale('log', basex=2)
     pylab.xticks(
         list(scipy.unique(group['symbols'])),
         list(scipy.unique(group['symbols'])))
     plotter.set_markers(p)
     plotter.set_slave_info(slavename)

    def find_neighbor_throats(self,pnums,flatten=True,mode='union'):
        r"""
        Returns a list of throats neighboring the given pore(s)

        Parameters
        ----------
        pnums : array_like
            Indices of pores whose neighbors are sought
        flatten : boolean, optional
            If flatten is True (default) a 1D array of unique throat ID numbers
            is returned. If flatten is False the returned array contains arrays
            of neighboring throat ID numbers for each input pore, in the order
            they were sent.
        mode : string, optional
            Specifies which neighbors should be returned.  The options are: 
            
            * 'union' : All neighbors of the input pores

            * 'intersection' : Only neighbors shared by all input pores 
            
            * 'not_intersection' : Only neighbors not shared by any input pores

        Returns
        -------
        neighborTs : 1D array (if flatten is True) or ndarray of arrays (if
            flatten if False)

        Examples
        --------
        >>> pn = OpenPNM.Network.Cubic(name='doc_test').generate(divisions=[5,5,5],lattice_spacing=[1])
        >>> pn.find_neighbor_throats(pnums=[0,1])
        array([0, 1, 2, 3, 4, 5])
        >>> pn.find_neighbor_throats(pnums=[0,1],flatten=False)
        array([array([0, 1, 2]), array([0, 3, 4, 5])], dtype=object)
        """
        #Test for existance of incidence matrix
        try:
            neighborTs = self.incidence_matrix['lil']['connections'].rows[[pnums]]
        except:
            self._logger.info('Creating incidence matrix, please wait')
            self.create_incidence_matrix()
            neighborTs = self.incidence_matrix['lil']['connections'].rows[[pnums]]
        if flatten:
            #All the empty lists must be removed to maintain data type after hstack (numpy bug?)
            neighborTs = [sp.asarray(x) for x in neighborTs if x]
            neighborTs = sp.hstack(neighborTs)
            #Remove references to input pores and duplicates
            if mode == 'not_intersection':
                neighborTs = sp.unique(sp.where(sp.bincount(neighborTs)==1)[0])
            elif mode == 'union':
                neighborTs = sp.unique(neighborTs)
            elif mode == 'intersection':
                neighborTs = sp.unique(sp.where(sp.bincount(neighborTs)>1)[0])
        else:
            for i in range(0,sp.size(pnums)):
                neighborTs[i] = sp.array(neighborTs[i])
        return sp.array(neighborTs,ndmin=1)

 def plot_setup(p):
     pylab.ylabel("Throughput" + " [" + list(group['unit'])[0] + "]")
     pylab.yscale('log')
     pylab.xscale('log', basex=2)
     pylab.xticks(
         list(scipy.unique(group['symbols'])),
         list(scipy.unique(group['symbols'])))
     plotter.set_markers(p)
     plotter.set_legend_columns(3)

 def test_no_late_filling(self):
     mip = op.algorithms.Porosimetry(network=self.net)
     mip.setup(phase=self.hg)
     mip.set_inlets(pores=self.net.pores('left'))
     mip.run()
     assert len(sp.unique(mip['pore.invasion_pressure'])) > 1
     assert len(sp.unique(mip['pore.invasion_sequence'])) > 1
     assert len(sp.unique(mip['throat.invasion_pressure'])) > 1
     assert len(sp.unique(mip['throat.invasion_sequence'])) > 1