Python scipy.where函数代码示例

def est_condprob2(data, val, given):
    """Calculate the probability of P(X|Y,Z)

    est_condprob2(data, 'A', ['M', 'LC'])"""

    if not isinstance(given, list):
        raise IndexError("Given must be a list or tuple of givens")
    elif len(given) != 2:
        raise IndexError("I need multiple givens!  Give me more...give me more!")

    gcols = []
    for g in given:
        if g in ['M', 'F']:
            gcols.append(1)
        elif g in ['LC', 'SC', 'T']:
            gcols.append(2)
        elif g in ['A', 'B', 'C']:
            gcols.append(0)

    if val in ['M', 'F']:
        vcol = 1
    elif val in ['LC', 'SC', 'T']:
        vcol = 2
    elif val in ['A', 'B', 'C']:
        vcol = 0

    datsize = data.shape[0]
    needed = [val, given[0], given[1]]
    t = sp.where([sp.all(data[i]==needed) for i in range(datsize)])[0]

    t2 = sp.where([sp.all(data[i,1:]==given) for i in range(datsize)])[0]
    
    return float(t.size)/t2.size

    def cross_validation(self, x, y, v=5, sig_r=2.0 ** sp.arange(-8, 0), mu_r=10.0 ** sp.arange(-15, 0)):
        # Get parameters
        n = x.shape[0]
        ns = sig_r.size
        nm = mu_r.size
        err = sp.zeros((ns, nm))

        # Initialization of the indices for the cross validation
        cv = CV()
        cv.split_data_class(y, v=v)

        for i in range(ns):
            for j in range(nm):
                for k in range(v):
                    model_temp = KDA()
                    model_temp.train(x[cv.it[k], :], y[cv.it[k]], sig=sig_r[i], mu=mu_r[j])
                    yp = model_temp.predict(x[cv.iT[k], :], x[cv.it[k], :], y[cv.it[k]])
                    yp.shape = y[cv.iT[k]].shape
                    t = sp.where(yp != y[cv.iT[k]])[0]
                    err[i, j] += float(t.size) / yp.size
                    del model_temp
        err /= v
        t = sp.where(err == err.min())
        self.sig = sig_r[t[0][0]]
        self.mu = mu_r[t[1][0]]
        return sig_r[t[0][0]], mu_r[t[1][0]], err

def compute_ndvi(im, r=0, ir=1, NODATA=-10000):
    """The function computes the NDVI of a multivalued image. It checks if
    there is NODATA value or division per zeros.
    
    Args:
    im: the image to process
    r: the number of the band that corresponds to the red band.
    ir: the number of the band that corresponds to the infra-red band.
    NODATA: the value of the NODATA
    
    Returns:
    ndvi =  the ndvi of the image
    """

    ## Get the size fo the image
    [nl, nc, nb] = im.shape

    ## Be sure that we can do 'floating operation'
    imf = im.astype(sp.float64)
    ndvi = sp.empty((nl, nc))

    if nb < 2:
        print "Two bands are needed to compute the NDVI"
        return None
    else:
        den = imf[:, :, ir - 1] + imf[:, :, r - 1]  # Pre compute the denominator
        t = sp.where((den > 0) & (imf[:, :, 1] != NODATA))
        ndvi[t] = (imf[t[0], t[1], ir - 1] - imf[t[0], t[1], r - 1]) / den[t]  # compute the ndvi on the safe samples

        if len(t[0]) < nl * nc:
            t = sp.where((den == 0) | (imf[:, :, 1] == NODATA))  # check for problematic pixels
            ndvi[t] = NODATA

        imf = []
        return ndvi

 def _generate_masked_mesh(self, cell_mask=None):
     r"""
     Generates the mesh based on the cell mask provided
     """
     #
     if cell_mask is None:
         cell_mask = sp.ones(self.data_map.shape, dtype=bool)
     #
     # initializing arrays
     self._edges = sp.ones(0, dtype=str)
     self._merge_patch_pairs = sp.ones(0, dtype=str)
     self._create_blocks(cell_mask)
     #
     # building face arrays
     mapper = sp.ravel(sp.array(cell_mask, dtype=int))
     mapper[mapper == 1] = sp.arange(sp.count_nonzero(mapper))
     mapper = sp.reshape(mapper, (self.nz, self.nx))
     mapper[~cell_mask] = -sp.iinfo(int).max
     #
     boundary_dict = {
         'bottom':
             {'bottom': mapper[0, :][cell_mask[0, :]]},
         'top':
             {'top': mapper[-1, :][cell_mask[-1, :]]},
         'left':
             {'left': mapper[:, 0][cell_mask[:, 0]]},
         'right':
             {'right': mapper[:, -1][cell_mask[:, -1]]},
         'front':
             {'front': mapper[cell_mask]},
         'back':
             {'back': mapper[cell_mask]},
         'internal':
             {'bottom': [], 'top': [], 'left': [], 'right': []}
     }
     #
     # determining cells linked to a masked cell
     cell_mask = sp.where(~sp.ravel(cell_mask))[0]
     inds = sp.in1d(self._field._cell_interfaces, cell_mask)
     inds = sp.reshape(inds, (len(self._field._cell_interfaces), 2))
     inds = inds[:, 0].astype(int) + inds[:, 1].astype(int)
     inds = (inds == 1)
     links = self._field._cell_interfaces[inds]
     #
     # adjusting order so masked cells are all on links[:, 1]
     swap = sp.in1d(links[:, 0], cell_mask)
     links[swap] = links[swap, ::-1]
     #
     # setting side based on index difference
     sides = sp.ndarray(len(links), dtype='<U6')
     sides[sp.where(links[:, 1] == links[:, 0]-self.nx)[0]] = 'bottom'
     sides[sp.where(links[:, 1] == links[:, 0]+self.nx)[0]] = 'top'
     sides[sp.where(links[:, 1] == links[:, 0]-1)[0]] = 'left'
     sides[sp.where(links[:, 1] == links[:, 0]+1)[0]] = 'right'
     #
     # adding each block to the internal face dictionary
     inds = sp.ravel(mapper)[links[:, 0]]
     for side, block_id in zip(sides, inds):
         boundary_dict['internal'][side].append(block_id)
     self.set_boundary_patches(boundary_dict, reset=True)

 def simulate(self,x0,lambd):
     Dt = self.param['Dt']
     # Dt needs to be a multiple of param['Dt']
     dt = self.param['dt']
     D = lambd[1]
     a = lambd[0]
     N = self.param['N']
     drift = self.param['drift']
     x = scipy.array(x0)
     
     tstart = 0
     tcur = tstart
     while (tcur < tstart + Dt + dt/2 ):
         tcur += dt
         # the random number
         dW=self.rand.normal(loc=0.,scale=scipy.sqrt(2*D*dt),size=N)
       #  if tcur == dt:    #only print random number for first time step
       #      print 'dW =',  dW            
         
         # the process
         drift_term = a * drift(x)
         x=x+drift_term*dt+dW
         # and reflecting boundary conditions
         scipy.where(x>self.domain[1],2*self.domain[1]-x,x)
         scipy.where(x<self.domain[0],2*self.domain[0]-x,x)
     return x

    def populate_out_of_dip_theta(self, n, dip):

        out_of_dip = asarray(self.populate_distribution(
            self.out_of_dip_theta_dist, n))

        (errorIndexes,) = where((out_of_dip > (175 - dip)) &
                                (out_of_dip < (185 - dip)))

        if len(errorIndexes) > 0:
            for i in errorIndexes:
                blnBadNum = True
                count = 0
                while blnBadNum:
                    newNum = self.populate_distribution(
                        self.out_of_dip_theta_dist, 1)
                    if ((newNum[0] <= (175 - dip)) | (newNum[0] >= (185 - dip))):
                        blnBadNum = False
                    count = count + 1
                    if count > 1000:
                        msg = "Bad out of dip theta range in fault \
                                     source file"
                        raise IOError(msg)
                out_of_dip[i] = newNum[0]
        (errorIndexes,) = where((out_of_dip > (175 - dip)) &
                                (out_of_dip < (185 - dip)))

        if len(errorIndexes) > 0:
            msg = "Bad out of dip theta range in fault \
                                     source file"
            raise IOError(msg)

        return out_of_dip

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 getEncodedData(filename,encoding="additive",phenotype_id=None,maf=0.0):
    f = h5py.File(filename,'r')
    phenotype_id = str(phenotype_id)
    if not phenotype_id==None:
        sample_ids = f['Genotype/sample_ids'][:]
        p_sample_ids = f['Phenotypes'][phenotype_id]['sample_ids'][:]
        y = f['Phenotypes'][phenotype_id]['y'][:]
        ind = sp.where(~sp.isnan(y))[0]
        y = y[ind]
        p_sample_ids = p_sample_ids[ind]
        ind = (sp.reshape(sample_ids,(sample_ids.shape[0],1))==p_sample_ids).nonzero()
        raw = f['Genotype/raw'][:]
        raw = raw[ind[0],:]
        [encoded,maf_v] = encodeHeterozygousData(raw)
        ind = sp.where(maf_v>=maf)[0]
        encoded = encoded[:,ind]
        identifiers = f['Genotype/identifiers'][:]
        identifiers = identifiers[ind]
        maf_v = maf_v[ind]
        f.close()
        return [encoded,maf_v,identifiers]
    if encoding=="additive":
        if 'encoded_additive' in f['Genotype'].keys():
            encoded = f['Genotype/encoded_additive'][:]
            maf_v = f['Genotype/global_maf'][:]
        else:
            [encoded,maf_v] = encodeHeterozygousData(f['Genotype/raw'][:])
    identifiers = f['Genotype/identifiers'][:]
    f.close()
    return [encoded,maf_v,identifiers]

    def newEpisode(self):
        if self.learning:
            params = ravel(self.explorationlayer.module.params)
            target = ravel(sum(self.history.getSequence(self.history.getNumSequences()-1)[2]) / 500)
        
            if target != 0.0:
                self.gp.addSample(params, target)
                if len(self.gp.trainx) > 20:
                    self.gp.trainx = self.gp.trainx[-20:, :]
                    self.gp.trainy = self.gp.trainy[-20:]
                    self.gp.noise = self.gp.noise[-20:]
                    
                self.gp._calculate()
                        
                # get new parameters where mean was highest
                max_cov = diag(self.gp.pred_cov).max()
                indices = where(diag(self.gp.pred_cov) == max_cov)[0]
                pick = indices[random.randint(len(indices))]
                new_param = self.gp.testx[pick]
            
                # check if that one exists already in gp training set
                if len(where(self.gp.trainx == new_param)[0]) > 0:
                    # add some normal noise to it
                    new_param += random.normal(0, 1, len(new_param))

                self.explorationlayer.module._setParameters(new_param)

            else:
                self.explorationlayer.drawRandomWeights()
        
        # don't call StateDependentAgent.newEpisode() because it randomizes the params
        LearningAgent.newEpisode(self)

 def cryptoInternal(self, data, base):
     addresses = scipy.array(range(base, base + (len(data) * 2), 2), scipy.uint32)
     for mask, xorVal in self.XOR_TABLE1:
         data = scipy.where((addresses & mask) == mask, data ^ xorVal, data)
     for mask, xorVal in self.XOR_TABLE2:
         data = scipy.where((addresses & mask) != 0,    data ^ xorVal, data)
     return data

 def smart_threshold(self):
     self.median = numpy.median(self.data)
     self.std = numpy.std(self.data)
     blank = scipy.where(self.data < self.median+0.25*self.std)
     signal = scipy.where(self.data > self.median+0.25*self.std)
     self.data[blank] = 0.0
     self.data[signal] = 1.0

 def maskLowStddVoxels(self, dds, nMeanDds, nStddDds):
     unique = np.unique(sp.where(nStddDds.subd.asarray() <= 1.0/3, dds.subd.asarray(), dds.mtype.maskValue()))
     unique = unique[sp.where(unique != dds.mtype.maskValue())]
     if (dds.mpi.comm != None):
         unique = dds.mpi.comm.allreduce(unique.tolist(), op=mpi.SUM)
         unique = np.unique(unique)
     rootLogger.info("Unique constant stdd values = %s" % (unique,))
     rootLogger.info("Creating mask from unique constant values...")
     mskDds = mango.zeros_like(dds, mtype="segmented")
     for uVal in unique:
         mskDds.asarray()[...] = sp.where(dds.asarray() == uVal, 1, mskDds.asarray())
     rootLogger.info("Done creating mask from unique constant values.")
 
     rootLogger.info("Labeling connected constant zero-stdd regions...")
     mskDds.updateHaloRegions()
     mskDds.mirrorOuterLayersToBorder(False)
     self.writeIntermediateDds("_000ZeroStddForLabeling", mskDds)
     lblDds = mango.image.label(mskDds, 1)
     rootLogger.info("Done labeling connected constant stdd regions.")
     self.writeIntermediateDds("_000ZeroStdd", lblDds)
     
     countThresh = 0.01 * sp.product(lblDds.shape)
     rootLogger.info("Eliminating large clusters...")
     lblDds = mango.image.eliminate_labels_by_size(lblDds, minsz=int(countThresh), val=lblDds.mtype.maskValue())
     self.writeIntermediateDds("_000ZeroStddLargeEliminated", lblDds)
 
     rootLogger.info("Assigning mask values...")
     mskDds.subd.asarray()[...] = \
         sp.where(lblDds.subd.asarray() == lblDds.mtype.maskValue(), True, False)
     self.writeIntermediateDds("_000ZeroStddMskd", mskDds)
     del lblDds
     for tmpDds in [dds, nMeanDds, nStddDds]:
         tmpDds.subd.asarray()[...] = \
             sp.where(mskDds.subd.asarray(), tmpDds.mtype.maskValue(), tmpDds.subd.asarray())

 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 step(self, *args):
     """First update the step size, then actually take a step along the gradient."""
     g = self.model.grad(*args);
     
     # Update the weighted Exponential sq avg.
     self.sqExpAvgGrad *= self.exponentAvgM;
     self.sqExpAvgGrad += (1-self.exponentAvgM) * g**2;
     self.sqExpAvgGrad[:] = where(self.sqExpAvgGrad < EPSILON, EPSILON, self.sqExpAvgGrad);
     
     # Uodate the muVect
     possUpdate = 1 + self.qLearningRate * g * self.expAvgGrad / self.sqExpAvgGrad
     #log.debug('max(possUpdate): %.4f,  min(possUpdate): %.4f' % (max(possUpdate), min(possUpdate)))
     ## Keep this from going negative.
     possUpdate = where(possUpdate < 0.001, 0.001, possUpdate);
     self.muVect *= possUpdate
     
     # Do something to cap the update rate.  This is allowing the step rate to overpower the decay completely
     self.muVect = where(self.muVect > self.maxMuVect, self.maxMuVect, self.muVect);
     
     # Then update the exponential average
     self.expAvgGrad *= self.exponentAvgM;
     self.expAvgGrad += (1-self.exponentAvgM) * g;
     
     self.model.params -= self.muVect * g
     Trainer.step(self,*args)

 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