Python numpy.divide函数代码示例

def multistate_distribution(data, parameters, limit, 
                            normalize_likelihood_level_cell_counts = True):

    data_grandpa, data_parent, data_children = data
    sigma, b, a_grandpa, a_parent, a_children = parameters
    
    normalization_factor = normalize(sigma, a_grandpa, b, limit)
    grandpa_dist = [steady_state_distribution(x, sigma, a_grandpa, b, normalization_factor) for x in data_grandpa]
    
    normalization_factor = normalize(sigma, a_parent, b, limit)
    parent_dist = [steady_state_distribution(x, sigma, a_parent, b, normalization_factor) for x in data_parent]
    
    normalization_factor = normalize(sigma, a_children, b, limit)
    children_dist = [steady_state_distribution(x, sigma, a_children, b, normalization_factor) for x in data_children]
    
    grandpa_dist = np.array(grandpa_dist, dtype = float)
    parent_dist = np.array(parent_dist, dtype = float)
    children_dist = np.array(children_dist, dtype = float)
    
    if normalize_likelihood_level_cell_counts:
        grandpa_dist = np.divide(grandpa_dist, float(data_grandpa.size))
        parent_dist = np.divide(parent_dist, float(data_parent.size))
        children_dist = np.divide(children_dist, float(data_children.size))
        
    return grandpa_dist, parent_dist, children_dist

def ratio_err(top,bottom,top_low,top_high,bottom_low,bottom_high):
    #uses simple propagation of errors (partial derivatives)
    #note it returns errorbars, not interval

    #-make sure input is numpy arrays-
    top = np.array(top)
    top_low = np.array(top_low)
    top_high = np.array(top_high)
    bottom = np.array(bottom)
    bottom_low = np.array(bottom_low)
    bottom_high = np.array(bottom_high)

    #-calculate errorbars-
    top_errlow = np.subtract(top,top_low)
    top_errhigh = np.subtract(top_high,top)
    bottom_errlow = np.subtract(bottom,bottom_low)
    bottom_errhigh = np.subtract(bottom_high,bottom)

    #-calculate ratio_low-
    ratio_low  = np.sqrt( np.square(np.divide(top_errlow,bottom)) + np.square( np.multiply(np.divide(top,np.square(bottom)),bottom_errlow)) )
    #-calculate ratio_high-
    ratio_high = np.sqrt( np.square(np.divide(top_errhigh,bottom)) + np.square( np.multiply(np.divide(top,np.square(bottom)),bottom_errhigh)) )
#    ratio_high = ((top_errhigh/bottom)**2.0 + (top/(bottom**2.0))*bottom_errhigh)**2.0)**0.5

    # return two vectors, err_low and err_high
    return ratio_low,ratio_high

def HS(im1, im2, alpha, ite,):

	#set up initial velocities
	uInitial = np.zeros([im1.shape[0],im1.shape[1]])
	vInitial = np.zeros([im1.shape[0],im1.shape[1]])

	# Set initial value for the flow vectors
	u = uInitial
	v = vInitial

	# Estimate derivatives
	[fx, fy, ft] = computeDerivatives(im1, im2)

	# Averaging kernel
	kernel=np.matrix([[1/12, 1/6, 1/12],[1/6, 0, 1/6],[1/12, 1/6, 1/12]])

	print fx[100,100],fy[100,100],ft[100,100]

	# Iteration to reduce error
	for i in range(ite):
		# Compute local averages of the flow vectors
		uAvg = cv2.filter2D(u,-1,kernel)
		vAvg = cv2.filter2D(v,-1,kernel)

		uNumer = (fx.dot(uAvg) + fy.dot(vAvg) + ft).dot(ft)
		uDenom = alpha + fx**2 + fy**2
		u = uAvg - np.divide(uNumer,uDenom)

		# print np.linalg.norm(u)

		vNumer = (fx.dot(uAvg) + fy.dot(vAvg) + ft).dot(ft)
		vDenom = alpha + fx**2 + fy**2
		v = vAvg - np.divide(vNumer,vDenom)
	return (u,v)

def  visResults (m, result_dir, varying_para_values = r_values, xlabel= 'Radius Scale Factor', basenum = 0):

        # plot the raw results
        plotResults(m,result_dir, varying_para_values,xlabel, 'Fitted')

        # plot the normalized results
        norm_m = np.zeros(m.shape)
        for i in range(len(varying_para_values)):
            norm_m[:,:,i] = np.divide(m[:,:,i], m[:,:,basenum]) # r=1.0 is the original model fitting results
        plotResults(norm_m,result_dir,varying_para_values,xlabel,'Normalized')

        for i in range(len(varying_para_values)):
            norm_m[:,:,i] = 100 * np.divide(m[:,:,i]- m[:,:,basenum], m[:,:,basenum]) # r=1.0 is the original model fitting results
        plotResults(norm_m,result_dir,varying_para_values,xlabel,'Difference')



        CmTotal = np.zeros((sample_size, len(varying_para_values)))
        for i in range(len(varying_para_values)):
            CmTotal[:,i] = m[:,4,i] * (m[:,1,i]+ m[:,2,i])  # Cm * (A1+A2) p=0.0 is the original model fitting results
        df_CmTotal = pd.DataFrame(CmTotal, columns = varying_para_values)
        my_box_plot(df_CmTotal, result_dir+ "/total_capacitance.png", xlabel,'Total Capacitance')



        RmUnit = np.zeros((sample_size, len(varying_para_values)))
        for i in range(len(varying_para_values)):
            RmUnit[:,i] = m[:,5,i] / (m[:,1,i]+ m[:,2,i])  # Rm / (A1+A2)
        df_RmUnit = pd.DataFrame(RmUnit, columns = varying_para_values)
        my_box_plot(df_RmUnit, result_dir+ "/total_Rm.png", xlabel,'Unit Membrane Resistance')

def CalcGamma(MeanVec,VarVec,wVec,X):

    #gamma = np.zeros(shape=(8,X.shape[0]),dtype='float128')
    for a in range(X.shape[0]):
        summ = 0
        for i in range(NoM):
            power = np.square(np.subtract(X[a],MeanVec[i]))
##            print 'poweris \n',power
##            print 'v o\n',VarVec[i]
            denp = np.multiply(-2,VarVec[i])
##            print denp
            power = np.divide(power,denp)
##            print power
            power = np.sum(power)
##            print 'power is \n',power

            power = exp(power)
##            print power
            prodVarVec = np.prod(VarVec[i])
            den = 1/(2*math.pi)**(NoM/2)*np.sqrt(prodVarVec)
##            print den
            gamma[i][a] = wVec[i]*np.divide(power,den)
##            print gamma
            summ = summ + gamma[i][a]
##            print gamma[i][a]
        for i in range(NoM):
            gamma[i][a] = gamma[i][a]/summ
        
    return gamma

    def __EM(self):
        old_log_like = -np.inf
        threshold = 1e-15
        probability = 0
        while True:
            # E step
            probability = self.__probability()
            expectation = np.multiply(probability, self.prior)
            expectation = np.divide(expectation, expectation.sum(axis=1))

            # M step: updata parameters
            sumk = expectation.sum(axis=0)
            self.prior = sumk / self.x.shape[0]
            self.mean = np.diag(np.array(np.divide(1, sumk)).flatten()) * \
                        expectation.T * self.x
            for i in range(self.k):
                x_shift = self.x - self.mean[i, :]
                self.sigma[:, :, i] = x_shift.T * \
                    np.diag(np.array(expectation[:, i]).flatten()) * x_shift /\
                    sumk[0, i]

            new_log_like = np.log(probability * self.prior.T).sum()
            if np.abs(new_log_like - old_log_like) < threshold:
                break
            old_log_like = new_log_like
        return probability

 def sw_sums(a, b):
     abw = apply_scale(w, a, b)
     np.divide(abw, 1 + abw, out = abw)
     abw[np.isnan(abw)] = 1
     swr = abw.sum(1, keepdims = True)
     swc = abw.sum(0, keepdims = True)
     return swr, swc

    def normalize(self, mode="integral"):
        """
        Force normalization of filter kernel.

        Parameters
        ----------
        mode : {'integral', 'peak'}
            One of the following modes:
                * 'integral' (default)
                    Kernel normalized such that its integral = 1.
                * 'peak'
                    Kernel normalized such that its peak = 1.
        """
        # There are kernel that sum to zero and
        # the user should be warned in this case
        if np.isinf(self._normalization):
            warnings.warn(
                "Kernel cannot be normalized because the " "normalization factor is infinite.", AstropyUserWarning
            )
            return
        if np.abs(self._normalization) > MAX_NORMALIZATION:
            warnings.warn(
                "Normalization factor of kernel is " "exceptionally large > {0}.".format(MAX_NORMALIZATION),
                AstropyUserWarning,
            )
        if mode == "integral":
            self._array *= self._normalization
        if mode == "peak":
            np.divide(self._array, self._array.max(), self.array)
            self._normalization = 1.0 / self._array.sum()

def get_gaussian_weight_patch(gauss_shape=(19, 19), gauss_sigma_frac=.3,
                              gauss_norm_01=True):
    r"""
    2d gaussian image useful for plotting

    Returns:
        ndarray: patch

    CommandLine:
        python -m vtool.coverage_kpts --test-get_gaussian_weight_patch

    Example:
        >>> # ENABLE_DOCTEST
        >>> from vtool.coverage_kpts import *  # NOQA
        >>> # build test data
        >>> # execute function
        >>> patch = get_gaussian_weight_patch()
        >>> # verify results
        >>> result = str(patch)
        >>> print(result)
    """
    # Perdoch uses roughly .95 of the radius
    radius = gauss_shape[0] / 2.0
    sigma = gauss_sigma_frac * radius
    # Similar to SIFT's computeCircularGaussMask in helpers.cpp
    # uses smmWindowSize=19 in hesaff for patch size. and 1.6 for sigma
    # Create gaussian image to warp
    patch = ptool.gaussian_patch(shape=gauss_shape, sigma=sigma)
    if gauss_norm_01:
        np.divide(patch, patch.max(), out=patch)
    return patch

def ndcg_multi(X, Y, Ks):
	assert(X.size == Y.size and all(X.indices == Y.indices) and all(X.indptr == Y.indptr))
	n = Y.shape[1]
	res = zeros(len(Ks))
	nvalid = 0
	Xdata = X.data
	Ydata = Y.data
	indices = Y.indices
	indptr = Y.indptr
	for i in xrange(n):
		[j0, j1] = [indptr[i], indptr[i + 1]]
		if j0 == j1: # skip empty column
			continue
		nvalid += 1
		Xi = Xdata[j0:j1]
		Yi = Ydata[j0:j1]
		I = argsort(-Xi)
		Yi_pred = numpy.exp(Yi[I])-1.0
		Yi_best = numpy.exp(-(sort(-Yi)))-1.0
		Wi = numpy.log(numpy.exp(1) + arange(j1 - j0))
		Yi_pred = numpy.divide(Yi_pred, Wi)
		Yi_best = numpy.divide(Yi_best, Wi)
		for k in xrange(len(Ks)):
			K = Ks[k]
			Ki = min([K, j1 - j0])
			res[k] += sum(Yi_pred[0:Ki]) / sum(Yi_best[0:Ki])
	assert(nvalid > 0)
	res /= nvalid
	return res

    def get_summed_cohp_by_label_and_orbital_list(self, label_list, orbital_list, divisor=1):
        """
        Returns a COHP object that includes a summed COHP divided by divisor

        Args:
            label_list: list of labels for the COHP that should be included in the summed cohp
            orbital_list: list of orbitals for the COHPs that should be included in the summed cohp (same order as label_list)
            divisor: float/int, the summed cohp will be divided by this divisor
        Returns:
            Returns a COHP object including a summed COHP
        """
        # check if cohps are spinpolarized or not
        first_cohpobject = self.get_orbital_resolved_cohp(label_list[0], orbital_list[0])
        summed_cohp = first_cohpobject.cohp.copy()
        summed_icohp = first_cohpobject.icohp.copy()
        for ilabel, label in enumerate(label_list[1:], 1):
            cohp_here = self.get_orbital_resolved_cohp(label, orbital_list[ilabel])
            summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp_here.cohp.copy()[Spin.up]], axis=0)
            if Spin.down in summed_cohp:
                summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp_here.cohp.copy()[Spin.down]], axis=0)
            summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], cohp_here.icohp.copy()[Spin.up]], axis=0)
            if Spin.down in summed_icohp:
                summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp_here.icohp.copy()[Spin.down]], axis=0)

        divided_cohp = {}
        divided_icohp = {}
        divided_cohp[Spin.up] = np.divide(summed_cohp[Spin.up], divisor)
        divided_icohp[Spin.up] = np.divide(summed_icohp[Spin.up], divisor)
        if Spin.down in summed_cohp:
            divided_cohp[Spin.down] = np.divide(summed_cohp[Spin.down], divisor)
            divided_icohp[Spin.down] = np.divide(summed_icohp[Spin.down], divisor)

        return Cohp(efermi=first_cohpobject.efermi, energies=first_cohpobject.energies, cohp=divided_cohp,
                    are_coops=first_cohpobject.are_coops,
                    icohp=divided_icohp)

    def _process_sample (self, ap1, ap2, ap3, triple, tflags):
        """We have computed one independent phase closure triple in one timeslot.

        """
        # Frequency-resolved:
        np.divide (triple, np.abs (triple), triple)
        phase = np.angle (triple)

        self.ap_spec_stats_by_ddid[self.cur_ddid].accum (ap1, phase, tflags + 0.)
        self.ap_spec_stats_by_ddid[self.cur_ddid].accum (ap2, phase, tflags + 0.)
        self.ap_spec_stats_by_ddid[self.cur_ddid].accum (ap3, phase, tflags + 0.)

        # Frequency-averaged:
        triple = np.dot (triple, tflags) / tflags.sum ()
        phase = np.angle (triple)

        self.global_stats_by_time.accum (self.cur_time, phase)

        self.ap_stats_by_ddid[self.cur_ddid].accum (ap1, phase)
        self.ap_stats_by_ddid[self.cur_ddid].accum (ap2, phase)
        self.ap_stats_by_ddid[self.cur_ddid].accum (ap3, phase)
        self.bp_stats_by_ddid[self.cur_ddid].accum ((ap1, ap2), phase)
        self.bp_stats_by_ddid[self.cur_ddid].accum ((ap1, ap3), phase)
        self.bp_stats_by_ddid[self.cur_ddid].accum ((ap2, ap3), phase)

        self.ap_time_stats_by_ddid[self.cur_ddid].accum (self.cur_time, ap1, phase)
        self.ap_time_stats_by_ddid[self.cur_ddid].accum (self.cur_time, ap2, phase)
        self.ap_time_stats_by_ddid[self.cur_ddid].accum (self.cur_time, ap3, phase)

def getGammaAngle(appf,cAtom,oAtom,hAtom):
    # first determine the nAtom
    aminoGroup = appf.select('resnum ' + str(cAtom.getResnum()))
    for at in aminoGroup:
        if(at.getName() == 'N'):
            nAtom = at
        # get coordinates
    cCoords = cAtom.getCoords()
    oCoords = oAtom.getCoords()
    hCoords = hAtom.getCoords()
    nCoords = nAtom.getCoords()
    # get necessary vectors
    oc = np.subtract(oCoords,cCoords)
    nc = np.subtract(nCoords,cCoords)
    ho = np.subtract(hCoords,oCoords)
    n1 = np.cross(oc,nc)
    n1_unit = np.divide(n1,np.linalg.norm(n1))
    # get projection of H-O in O-C direction
    oc_unit = np.divide(oc,np.linalg.norm(oc))
    #print oc_unit
    hproj = np.dot(ho,oc_unit)
    # get projection of H-O onto N-C-O plane
    out = np.dot(ho,n1_unit)
    n2 = np.cross(np.multiply(n1_unit,out),oc)
    #print n2
    ho_ip = np.subtract(ho,np.multiply(n1_unit,out))
    test = np.dot(n2,ho_ip)
    #print test
    ang = hproj/np.linalg.norm(ho_ip)
    ang = math.acos(ang)
    ang = ang*180/math.pi
    #if(test < 0):
    #    ang = ang * -1
    return ang

def compile_stats():
    logging.info('Loading data...')

    for i, filename in enumerate(os.listdir(DATA_DIR)):
        if i > MAX_FILES:
            break
        full_name = DATA_DIR + "/" + filename
        print full_name
        with open(full_name) as f:
            data = np.load(f)
            if len(data["input"].shape) == 3 and len(data["output"].shape) == 3:
                X = data["input"]
                y = data["output"]
                if COLLECT_ACTION_PCT:
                    probs = np.divide(np.sum(y, (0,1)).astype(np.float64),
                                      np.sum(y, (0,1,2)))
                    player_to_stats[filename] = probs
                elif COLLECT_VPIP_PFR:
                    # Assumes actions are (fold, check, call, raise)
                    actions = X[:,:,11:15]
                    if (np.sum(actions, (0,1,2)) < 100):
                        continue
                    probs = np.divide(np.sum(actions, (0,1)).astype(np.float64),
                                      np.sum(actions, (0,1,2)))
                    pfr = probs[3]
                    vpip = probs[3] + probs[2]
                    player_to_stats[filename] = np.array([vpip, pfr])

def _generate(l, k, g, beta, M, e, A, mu, intercept):

    p = beta.shape[0]

    if intercept:
        gradL1 = grad_l1(beta[1:, :])
        gradL2 = grad_l2_squared(beta[1:, :])
        gradGLmu = grad_glmu(beta[1:, :], A, mu)
    else:
        gradL1 = grad_l1(beta)
        gradL2 = grad_l2_squared(beta)
        gradGLmu = grad_glmu(beta, A, mu)

    alpha = -(l * gradL1 + k * gradL2 + g * gradGLmu)
    Mte = np.dot(M.T, e)
    if intercept:
        alpha = np.divide(alpha, Mte[1:, :])
    else:
        alpha = np.divide(alpha, Mte)

    X = np.ones(M.shape)
    if intercept:
        for i in xrange(p - 1):
            X[:, i + 1] = M[:, i + 1] * alpha[i, 0]
    else:
        for i in xrange(p):
            X[:, i] = M[:, i] * alpha[i, 0]

    y = np.dot(X, beta) - e

    return X, y