Python numpy.vdot函数代码示例

def berrycurvature(H,kx,ky,mu,d=10**-6):
    '''
    Calculate the Berry curvature of the occupied bands for a Hamiltonian with the given chemical potential using the Kubo formula.

    Parameters
    ----------
    H : callable
        Input function which returns the Hamiltonian as a 2D array.
    kx,ky : float
        The two parameters which specify the 2D point at which the Berry curvature is to be calculated.
        They are also the input parameters to be conveyed to the function H.
    mu : float
        The chemical potential.
    d : float, optional
        The spacing to be used to calculate the derivatives.

    Returns
    -------
    float
        The calculated Berry curvature for function H at point kx,ky with chemical potential mu.
    '''
    result=0
    Vx=(H(kx+d,ky)-H(kx-d,ky))/(2*d)
    Vy=(H(kx,ky+d)-H(kx,ky-d))/(2*d)
    Es,Evs=eigh(H(kx,ky))
    for n in range(Es.shape[0]):
        for m in range(Es.shape[0]):
            if Es[n]<=mu and Es[m]>mu:
                result-=2*(np.vdot(np.dot(Vx,Evs[:,n]),Evs[:,m])*np.vdot(Evs[:,m],np.dot(Vy,Evs[:,n]))/(Es[n]-Es[m])**2).imag
    return result

    def testMCMCFitting(self):
        print "starting mcmc"
        true, data = self.generate_data()
        model = HDPMixtureModel(3,100,100,1)
        model.seed = 1
        start = time()
        r = model.fit(data, verbose=10)
#        end = time() - start
#        print r.pis.shape
#        print r.pis, r.pis[0]
#        print r[0],r[0].pis
        print r.mus
        self.assertEqual(len(r), 2, 'results are the wrong length: %d' %len(r))
        diffs = {}
        #print 'r.mus:', r.mus()
        for i in gen_mean:
            diffs[i] = np.min(np.abs(r[0].mus-gen_mean[i]),0)
            #print r[0].mus[0], np.min(np.abs(r[0].mus[0]-gen_mean[i]),0)
            #diffs[i] = np.abs(r[0].mus()[i]-gen_mean[i])
            #print i, gen_mean[i],r[0].mus()[i], diffs[i], np.vdot(diffs[i],diffs[i])
            self.assertLessEqual( np.vdot(diffs[i],diffs[i]),1, 
                                  'diff to large: %f' % np.vdot(diffs[i], diffs[i]))
        
        for i in gen_mean:
            diffs[i] = np.min(np.abs(r[1].mus-gen_mean[i]),0)
            #diffs[i] = np.abs(r[1].mus()[i]-gen_mean[i])
            #print i, gen_mean[i],r[1].mus()[i], diffs[i], np.vdot(diffs[i],diffs[i])
            self.assertLessEqual( np.vdot(diffs[i],diffs[i]),1, 
                                  'diff to large: %f' % np.vdot(diffs[i], diffs[i]))

    def test_integral_moment_first(self):
        lmax = 5 #XXX use meshgrid above l=5

        # Test first-order moments of the spherical harmonics
        for l1,m1 in lmiter(lmax, comm=world):
            s1_L = Y(l1, m1, theta_L, phi_L)
            for l2,m2 in lmiter(lmax):
                s2_L = Y(l2, m2, theta_L, phi_L)

                # Note that weights times surface area make up for sin(theta)
                v_ex = 4 * np.pi * np.vdot(s1_L, np.cos(phi_L) \
                    * np.sin(theta_L) * s2_L * weight_L)
                v_ey = 4 * np.pi * np.vdot(s1_L, np.sin(phi_L) \
                    * np.sin(theta_L) * s2_L * weight_L)
                v_ez = 4 * np.pi * np.vdot(s1_L, np.cos(theta_L) \
                    * s2_L * weight_L)

                v0_ex = intYY_ex(l1, m1, l2, m2)
                v0_ey = intYY_ey(l1, m1, l2, m2)
                v0_ez = intYY_ez(l1, m1, l2, m2)

                self.assertAlmostEqual(v_ex, v0_ex, 12, '%s != %s (l1=%2d, ' \
                    'm1=%2d, l2=%2d, m2=%2d)' % (v_ex,v0_ex,l1,m1,l2,m2))
                self.assertAlmostEqual(v_ey, v0_ey, 12, '%s != %s (l1=%2d, ' \
                    'm1=%2d, l2=%2d, m2=%2d)' % (v_ey,v0_ey,l1,m1,l2,m2))
                self.assertAlmostEqual(v_ez, v0_ez, 12, '%s != %s (l1=%2d, ' \
                    'm1=%2d, l2=%2d, m2=%2d)' % (v_ez,v0_ez,l1,m1,l2,m2))

def neb_modify_force(band):
	peak_image = numpy.argmax(band.Energy)
	for i in xrange(1, band.nimages-1):
		# "improved" tangent calculation
		Forwards = putil.separation_vector_raw(band.Pos[(i+1)*3*band.NAtoms:(i+2)*3*band.NAtoms], band.Pos[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], band.PBC, band.Dim)
		Backwards = putil.separation_vector_raw(band.Pos[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], band.Pos[(i-1)*3*band.NAtoms:(i)*3*band.NAtoms], band.PBC, band.Dim)
		mForwards = putil.magnitude(Forwards)
		mBackwards = putil.magnitude(Backwards)
	
		if(band.Energy[i]>band.Energy[i-1] and band.Energy[i+1] > band.Energy[i]):
			norm = putil.normalise(Forwards)
		elif(band.Energy[i] < band.Energy[i-1] and band.Energy[i+1] < band.Energy[i]):
			norm = putil.normalise(Backwards)
		else:
			norm = putil.normalise(Backwards * math.fabs(band.Energy[i] - band.Energy[i-1]) + Forwards * math.fabs(band.Energy[i+1] - band.Energy[i]))

			
		# if climbing then move uphill independent of the rest of the band
		if(i==peak_image):
			band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms] -= 2 * norm * numpy.vdot(band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], norm)
			
		else:
			# remove parallel component of real force
			band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms] -= norm * numpy.vdot(band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms], norm)
			# add in force due to the springs
			band.Force[(i)*3*band.NAtoms:(i+1)*3*band.NAtoms] += norm*(mForwards - mBackwards)*band.springk
			
		#band.mForce = numpy.absolute(band.Force).max()
		band.mForce = putil.magnitude(band.Force)
		band.peak_image = peak_image
		

		
	return band	

def adjointness_error(op, its=100):
    """Check adjointness of op.A and op.As for 'its' instances of random data.

    For random unit-normed x and y, this finds the error in the adjoint
    identity <Ax, y> == <x, A*y>:
        err = abs( vdot(A(x), y) - vdot(x, Astar(y)) ).

    The type and shape of the input to A are specified by inshape and indtype.

    Returns a vector of the error magnitudes.

    """
    inshape = op.inshape
    indtype = op.indtype
    outshape = op.outshape
    outdtype = op.outdtype

    x = get_random_normal(inshape, indtype)

    errs = np.zeros(its, dtype=indtype)
    for k in xrange(its):
        x = get_random_normal(inshape, indtype)
        x = x/np.linalg.norm(x)
        y = get_random_normal(outshape, outdtype)
        y = y/np.linalg.norm(y)
        ip_A = np.vdot(op.A(x), y)
        ip_Astar = np.vdot(x, op.As(y))
        errs[k] = np.abs(ip_A - ip_Astar)

    return errs

    def braket(self, other, space="g"):
        """
        Returns the scalar product <u1|u2> of the periodic part of two wavefunctions 
        computed in G-space or r-space, depending on the value of space.

        Args:
            other: 
                Other wave (right-hand side)
            space: 
                Integration space. Possible values ["g", "gsphere", "r"]
                if "g" or "r" the scalar product is computed in G- or R-space on the FFT box.
                if space="gsphere" the integration is done on the G-sphere. Note that
                this option assumes that self and other have the same list of G-vectors. 
        """
        space = space.lower()

        if space == "g":  
            ug1_mesh = self.gsphere.tofftmesh(self.mesh, self.ug)
            ug2_mesh = other.gsphere.tofftmesh(self.mesh, other.ug)
            return np.vdot(ug1_mesh, ug2_mesh)

        elif space == "gsphere":
            return np.vdot(self.ug, other.ug)

        elif space == "r":
            return np.vdot(self.ur, other.ur) * self.mesh.dv

        else:
            raise ValueError("Wrong space: %s" % space)

    def  __compute_velocities_at_collision(obj1, obj2):

        v_obj1 = np.array([obj1.velocity.x, obj1.velocity.y], dtype=np.float64)
        v_obj2 = np.array([obj2.velocity.x, obj2.velocity.y], dtype=np.float64)

        v_normal = np.array([v_obj1[0] - v_obj2[0], v_obj1[1] - v_obj2[1]])
        #print v_obj1, v_obj2
        #print v_normal

        if not np.linalg.norm(v_normal):
            #print "No colision..."
            return obj1.velocity, obj2.velocity

        v_unit_normal = np.divide(v_normal, np.linalg.norm(v_normal))
        v_unit_tangent = np.array([-1*v_unit_normal[1], v_unit_normal[0]])

        v_obj1_normal = np.vdot(v_unit_normal, v_obj1)
        v_obj1_tangent = np.vdot(v_unit_tangent, v_obj1)
        v_obj2_normal = np.vdot(v_unit_normal, v_obj2)
        v_obj2_tangent = np.vdot(v_unit_tangent, v_obj2)

        v_obj1_new_normal = (v_obj1_normal*(obj1.mass - obj2.mass) + 2*obj2.mass*v_obj2_normal)/(obj1.mass + obj2.mass)
        v_obj2_new_normal = (v_obj2_normal*(obj2.mass - obj1.mass) + 2*obj1.mass*v_obj1_normal)/(obj1.mass + obj2.mass)

        vec_obj1_normal = np.dot(v_obj1_new_normal, v_unit_normal)
        vec_obj1_tangent = np.dot(v_obj1_tangent, v_unit_tangent)
        vec_obj2_normal = np.dot(v_obj2_new_normal, v_unit_normal)
        vec_obj2_tangent = np.dot(v_obj2_tangent, v_unit_tangent)

        v_obj1_final = np.add(vec_obj1_normal, vec_obj1_tangent)
        v_obj2_final = np.add(vec_obj2_normal, vec_obj2_tangent)

        return Velocity(v_obj1_final[0], v_obj1_final[1]), Velocity(v_obj2_final[0], v_obj2_final[1])

    def step(self,r,f):
        r = np.array(r)
        f = np.reshape(f,(-1,3))
        if self.v is None:
            self.v = np.zeros((len(f), 3))
        else:
            try:
                vf = np.vdot(self.v,f)
            except ValueError:
                self.v = np.zeros((len(f), 3))
                vf = np.vdot(self.v,f)

            if vf > 0.0:
                self.v = ((1.0 - self.a) * self.v + 
                          (self.a * f / 
                           np.sqrt(np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))))
                if self.Nsteps > self.Nmin:
                    self.dt = min(self.dt * self.finc, self.dtmax)
                    self.a *= self.fa
                self.Nsteps += 1
            else:
                self.v *= 0.0
                self.a = self.astart
                self.dt *= self.fdec
                self.Nsteps = 0

        self.v += self.dt * f
        dr = self.dt * self.v
        dr = self.determine_step(dr)
        return r + dr.reshape(r.shape)

def plot_ritz(A, n, iters):
    ''' Plot the relative error of the Ritz values of `A'.
    '''
    Amul = A.dot
    b = np.random.rand(A.shape[0])
    Q = np.empty((len(b), iters+1), dtype = np.complex128)
    H = np.zeros((iters+1, iters), dtype = np.complex128)
    Q[:, 0] = b / la.norm(b)
    eigvals = np.sort(abs(la.eig(A)[0]))[::-1]
    eigvals = eigvals[:n]
    abs_err = np.zeros((iters,n))

    for j in xrange(iters):
        Q[:, j+1] = Amul(Q[:, j])
        for i in xrange(j+1):
            H[i,j] = np.vdot(Q[:,i].conjugate(), (Q[:, j+1]))
            Q[:,j+1] = Q[:,j+1] - H[i,j] * (Q[:,i])

        H[j+1, j] = np.sqrt(np.vdot(Q[:, j+1], Q[:, j+1].conjugate()))
        Q[:,j+1] = Q[:,j+1] / H[j+1, j]

        if j < n:
            rit = np.zeros(n, dtype = np.complex128)
            rit[:j+1] = np.sort(la.eig(H[:j+1, :j+1])[0])[::-1]
            abs_err[j,:] = abs(eigvals - rit) / abs(eigvals)
        else:
            rit = np.sort(la.eig(H[:j+1,:j+1])[0])[::-1]
            rit = rit[:n]
            abs_err[j,:] = abs(eigvals - rit) / abs(eigvals)

    for i in xrange(n):
        plt.semilogy(abs_err[:,i])
    plt.show()

def bilinear_concentric_potential(r_g, dr_g, f_g, ft_g, l1, l2, alpha, rfilt=None):
    """Calculate corrections for concentric functions and potentials::

                 /     _   _a    _   _a    ~   _   _a  ~ _   _a   _
        v      = | f  (r - R ) V(r - R ) - f  (r - R ) V(r - R ) dr
         m1,m2   /  L1,L2                    L1,L2

    where f(r) and ft(r) are bilinear product of two localized functions which
    are radial splines times real spherical harmonics (l1,m1) or (l2,m2) and::

          _       1       _ -1              ~ _    erf(alpha*r)  _ -1
        V(r) = --------- |r|        ^       V(r) = ------------ |r|
               4 pi eps0                            4 pi eps0

    Note that alpha (and rfilt) should conform with the cutoff radius.
    """
    work_g = erf(alpha*r_g)

    if rfilt is None:
        M = np.vdot(f_g - ft_g * work_g, r_g * dr_g)
    else:
        M = np.vdot((f_g - ft_g * work_g)[r_g>=rfilt], \
            (r_g * dr_g)[r_g>=rfilt])

        # Replace 1/r -> (3-r^2/rfilt^2)/(2*rfilt) for r < rfilt
        M += np.vdot((f_g - ft_g * work_g)[r_g<rfilt], \
            (r_g**2/(2*rfilt) * (3-(r_g/rfilt)**2) * dr_g)[r_g<rfilt])

    v_mm = np.empty((2*l1+1, 2*l2+1), dtype=float)
    for m1 in range(2*l1+1):
        for m2 in range(2*l2+1):
            v_mm[m1,m2] = M * intYY(l1, m1-l1, l2, m2-l2)
    return v_mm

    def calc_orbital_params(self):
        #Calculates Keplerian orbital parameters based on the state vectors
        #This method should be called after each change to velocity vector

        #a = - mu / (v^2 - 2 * mu/r)
        v_2 = numpy.vdot(self.v, self.v)
        r_ = numpy.linalg.norm(self.r)
        # TODO: this won't work for parabolic trajectories (as the denominator is 0)!
        self.a = -EARTH_MU / (v_2 - 2 * EARTH_MU / r_)

        #T = 2*Pi*sqrt(a^3/ni)
        # TODO: again, orbital period is not defined for non-bound orbits
        self.T = 2.0 * numpy.pi * numpy.sqrt(self.a**3 / EARTH_MU)

        #Calculate specific relative angular momentum h = r X v
        h = numpy.cross(self.r, self.v)
        h_2 = numpy.vdot(h, h)
        h_ = numpy.sqrt(h_2)

        #Calculate eccentricity vector e = (v X h) / EARTH_MU - r/|r|
        e = numpy.cross(self.v, h) / EARTH_MU - self.r/r_
        self.e = numpy.linalg.norm(e)

        i_rad = numpy.arccos(h[2] / h_) #h[2] = hz
        self.i = numpy.degrees(i_rad)
        #However, some soruces state that if hz < 0 then inclination is 180 deg - i; should check this
        #n is the vector pointing to the ascending node
        n = numpy.array((-h[1], h[0], 0))
        n_ = numpy.linalg.norm(n)
        if i_rad == 0.0:
            o_rad = 0.0
        else:
            if n[1] >= 0.0: #ie. if h[0] >= 0
                o_rad = numpy.arccos(n[0] / n_)
            else:
                o_rad = 2 * numpy.pi - numpy.arccos(n[0] / n_)
        self.o = numpy.degrees(o_rad)

        #Calculate ni (true anomaly)
        q = numpy.vdot(self.r, self.v) #What the hell is q?
        ni_x = h_2 / (r_ * EARTH_MU) - 1.0
        ni_y = h_ * q / (r_ * EARTH_MU)
        self.ni = numpy.degrees(numpy.arctan2(ni_y, ni_x))

        if self.e == 0.0:
            #For circular orbit w is 0 by convention
            self.w = 0.0
        else:
            if n_ == 0.0:
                #For equatorial orbit
                self.w = numpy.degrees(numpy.arctan2(e[1], e[0]))
            else:
                self.w = numpy.degrees(numpy.arccos(numpy.vdot(n, e) / (n_ * self.e)))
        if e[2] < 0.0:
            self.w = 360.0 - self.w
        if self.w < 0.0:
            self.w = 360.0 + self.w

        self.rp = self.a * (1.0 - self.e) #Periapsis distance
        self.ra = self.a * (1.0 + self.e) #Apoapsis distance

def CGSolve(u0x,u0y,lamb,mu,b,epsilon,dfx,dfy) :
    # Solves JTJ[ux,uy]=b
    #lambd,mu,epsilon,dfx,dfy are needed in the computation of JTJ
    # [u0x,u0y] is the starting point of the iteration algo
    nitmax=100;
    ux=u0x;
    uy=u0y;
    # Computes JTJu
    Ax,Ay=JTJ(ux,uy,dfx,dfy,lamb,mu,epsilon);
    rx=b[0]-Ax
    ry=b[1]-Ay
    px=rx
    py=ry
    rsold=np.linalg.norm(rx)**2+np.linalg.norm(ry)**2
    for i in range(nitmax) :
        Apx,Apy=JTJ(px,py,dfx,dfy,lamb,mu,epsilon);
        alpha=rsold/(np.vdot(rx[:],Apx[:])+np.vdot(ry[:],Apy[:]))
        ux=ux+alpha*px
        uy=uy+alpha*py
        rx=rx-alpha*Apx
        ry=ry-alpha*Apy
        rsnew=np.linalg.norm(rx)**2+np.linalg.norm(ry)**2
        if np.sqrt(rsnew)<1e-10 :
            return [ux,uy]
        px=rx+rsnew/rsold*px
        py=ry+rsnew/rsold*py
        rsold=rsnew
    return [ux,uy]

def berryphase(H,path,ns):
    '''
    Calculate the Berry phase of some bands of a Hamiltonian along a certain path.

    Parameters
    ----------
    H : callable
        Input function which returns the Hamiltonian as a 2D array.
    path : iterable of dict
        The path along which to calculate the Berry phase.
    ns : iterable of int
        The sequences of bands whose Berry phases are wanted.

    Returns
    -------
    1d ndarray
        The wanted Berry phase of the selected bands.
    '''
    ns=np.array(ns)
    for i,parameters in enumerate(path):
        new=eigh(H(**parameters))[1][:,ns]
        if i==0:
            result=np.ones(len(ns),new.dtype)
            evs,old=new,new
        else:
            for j in range(len(ns)):
                result[j]*=np.vdot(old[:,j],new[:,j])
            old=new
    else:
        for j in range(len(ns)):
            result[j]*=np.vdot(old[:,j],evs[:,j])
    return np.angle(result)/np.pi

    def step(self, f):
        coords = self.coords
        if self.v is None:
            self.v = np.zeros((len(coords)))
        else:
            vf = np.vdot(f, self.v)
            if vf > 0.0:
                self.v = (1.0 - self.a) * self.v + self.a * f / np.sqrt(
                    np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))
                if self.Nsteps > self.Nmin:
                    self.dt = min(self.dt * self.finc, self.dtmax)
                    self.a *= self.fa
                self.Nsteps += 1
            else:
                self.v[:] *= 0.0
                self.a = self.astart
                self.dt *= self.fdec
                self.Nsteps = 0

        self.v += self.dt * f
        dr = self.dt * self.v
        if False:  # how do we determine maxstep?
            normdr = np.sqrt(np.vdot(dr, dr))
        else:
            normdr = max(np.abs(dr))
        if normdr > self.maxstep:
            dr = self.maxstep * dr / normdr
        self.coords = coords + dr

    def step(self,f):
        atoms = self.atoms
        if self.v is None:
            self.v = np.zeros((len(atoms), 3))
        else:
            vf = np.vdot(f, self.v)
            if vf > 0.0:
                self.v = (1.0 - self.a) * self.v + self.a * f / np.sqrt(
                    np.vdot(f, f)) * np.sqrt(np.vdot(self.v, self.v))
                if self.Nsteps > self.Nmin:
                    self.dt = min(self.dt * self.finc, self.dtmax)
                    self.a *= self.fa
                self.Nsteps += 1
            else:
                self.v[:] *= 0.0
                self.a = self.astart
                self.dt *= self.fdec
                self.Nsteps = 0

        self.v += self.dt * f
        dr = self.dt * self.v
        normdr = np.sqrt(np.vdot(dr, dr))
        if normdr > self.maxmove:
            dr = self.maxmove * dr / normdr
        r = atoms.get_positions()
        atoms.set_positions(r + dr)
        self.dump((self.v, self.dt))