Python LinearAlgebra类代码示例

def linePlaneIntersectionNumeric(p1, p2, p3, p4, p5):
  if not useNumeric:
    return linePlaneIntersection(p1, p2, p3, p4, p5)
  if useNumpy:
    top = [
        [1., 1., 1., 1.],
        [p1[0], p2[0], p3[0], p4[0]], [p1[1], p2[1], p3[1], p4[1]],
        [p1[2], p2[2], p3[2], p4[2]]]
    topDet = numpy.linalg.det(top)
    bottom = [
        [1., 1., 1., 0.], [p1[0], p2[0], p3[0], p5[0]-p4[0]],
        [p1[1], p2[1], p3[1], p5[1]-p4[1]], [p1[2], p2[2], p3[2], p5[2]-p4[2]]]
    botDet = numpy.linalg.det(bottom)
  else:  # actually use numeric
    top = Matrix.Matrix(
        [[1., 1., 1., 1.], [p1[0], p2[0], p3[0], p4[0]], [p1[1], p2[1],
         p3[1], p4[1]], [p1[2], p2[2], p3[2], p4[2]]])
    topDet = LinearAlgebra.determinant(top)
    bottom = Matrix.Matrix(
        [[1., 1., 1., 0.], [p1[0], p2[0], p3[0], p5[0]-p4[0]], [p1[1],
         p2[1], p3[1], p5[1]-p4[1]], [p1[2], p2[2], p3[2], p5[2]-p4[2]]])
    botDet = LinearAlgebra.determinant(bottom)
  if topDet == 0.0 or botDet == 0.0:
    return False
  t = -topDet/botDet
  x = p4[0] + (p5[0]-p4[0]) * t
  y = p4[1] + (p5[1]-p4[1]) * t
  z = p4[2] + (p5[2]-p4[2]) * t
  return [x, y, z]

 def testgeneralizedInverse (self):
     "Test LinearAlgebra.generalized_inverse"
     import LinearAlgebra
     ca = Numeric.array([[1,1-1j],[0,1]])
     cai = LinearAlgebra.inverse(ca)
     cai2 = LinearAlgebra.generalized_inverse(ca)
     self.failUnless(eq(cai, cai2))

    def __init__(self, atoms, constraints):
	self.atoms = atoms
	natoms = len(self.atoms)
	nconst = reduce(operator.add, map(len, constraints))
	b = Numeric.zeros((nconst, natoms), Numeric.Float)
	c = Numeric.zeros((nconst,), Numeric.Float)
	i = 0
	for cons in constraints:
	    cons.setCoefficients(self.atoms, b, c, i)
	    i = i + len(cons)
	u, s, vt = LinearAlgebra.singular_value_decomposition(b)
	self.rank = 0
	for i in range(min(natoms, nconst)):
	    if s[i] > 0.:
		self.rank = self.rank + 1
	self.b = b
	self.bi = LinearAlgebra.generalized_inverse(b)
	self.p = Numeric.identity(natoms)-Numeric.dot(self.bi, self.b)
	self.c = c
	self.bi_c = Numeric.dot(self.bi, c)
	c_test = Numeric.dot(self.b, self.bi_c)
	if Numeric.add.reduce((c_test-c)**2)/nconst > 1.e-12:
	    Utility.warning("The charge constraints are inconsistent."
			    " They will be applied as a least-squares"
			    " condition.")

    def render(self, rstate):
        # Set up the texture matrix as the modelview inverse
        m = glGetFloatv(GL_MODELVIEW_MATRIX)
        glMatrixMode(GL_TEXTURE)
        m = LinearAlgebra.inverse(m)
        m[3][0] = 0
        m[3][1] = 0
        m[3][2] = 0
        glLoadMatrixf(m)
        glMatrixMode(GL_MODELVIEW)

        # Set up texture coordinate generation
        glTexEnvfv(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
        glTexGenfv(GL_S, GL_TEXTURE_GEN_MODE, GL_REFLECTION_MAP_EXT)
        glTexGenfv(GL_T, GL_TEXTURE_GEN_MODE, GL_REFLECTION_MAP_EXT)
        glTexGenfv(GL_R, GL_TEXTURE_GEN_MODE, GL_REFLECTION_MAP_EXT)
        glEnable(GL_TEXTURE_GEN_S)
        glEnable(GL_TEXTURE_GEN_T)
        glEnable(GL_TEXTURE_GEN_R)

        # We're blending on top of existing polygons, so use the same tricks as the decal pass
        DecalRenderPass.render(self, rstate)

        # Clean up
        glMatrixMode(GL_TEXTURE)
        glLoadIdentity()
        glMatrixMode(GL_MODELVIEW)
        glDisable(GL_TEXTURE_GEN_S)
        glDisable(GL_TEXTURE_GEN_T)
        glDisable(GL_TEXTURE_GEN_R)

 def __pow__(self, other):
     shape = self.array.shape
     if len(shape)!=2 or shape[0]!=shape[1]:
         raise TypeError, "matrix is not square"
     if type(other) in (type(1), type(1L)):
         if other==0:
             return Matrix(identity(shape[0]))
         if other<0:
             result=Matrix(LinearAlgebra.inverse(self.array))
             x=Matrix(result)
             other=-other
         else:
             result=self
             x=result
         if other <= 3:
             while(other>1):
                 result=result*x
                 other=other-1
             return result
         # binary decomposition to reduce the number of Matrix
         #  Multiplies for other > 3.
         beta = _binary(other)
         t = len(beta)
         Z,q = x.copy(),0
         while beta[t-q-1] == '0':
             Z *= Z
             q += 1
         result = Z.copy()
         for k in range(q+1,t):
             Z *= Z
             if beta[t-k-1] == '1':
                 result *= Z
         return result
     else:
         raise TypeError, "exponent must be an integer"

def plotit(xs,ys,title,legends):
    #
    # Do it 
    #
    num_points=0
    for x in xs:
        num_points=num_points+len(x)
    print num_points,'data points'
    #x=[1,2,3,4,5,6,7,8,9]
    #y=[2,4,6,8,10,12,14,16,18]
    mat_fix=[]
    vec_fix=[]
    for x in xs:
        for point in x:
            mat_fix.append([point,1.0])
    for y in ys:
        for point in y:
            vec_fix.append(point)
    import LinearAlgebra
    sols,rsq,rank,junk=LinearAlgebra.linear_least_squares(Numeric.array(mat_fix),
                                                          Numeric.array(vec_fix))
    slope=sols[0]
    intercept=sols[1]
    print rsq
    rsq=float(rsq[0])
    print 'Slope: %.2f, Intercept: %.2f, R^2: %.2f' %(slope,intercept,rsq)
    file=dislin_driver.graf_mult3(xs,ys,title,'Simple E','PBE_ene',legends)
    return

    def testLinearLeastSquares2(self):
        """
        From bug #503733. Failing with dlapack_lite
        """
        import LinearAlgebra
        d = [0.49910197] + [0.998203938] * 49
        d1 = [0.000898030454] * 50
        def tridiagonal(sz):
            G = Numeric.zeros((sz,sz), Numeric.Float64)
            for i in range(sz):
                G[i,i] = d[i]
            for i in range(0,sz-1):
                G[i+1,i] = G[i,i+1] = d1[i]
            return G

        yfull = Numeric.array(
            [4.37016668e+18, 4.09591905e+18, 3.82167167e+18, 4.12952660e+18,
             2.60084719e+18, 2.05944452e+18, 1.69850960e+18, 1.51450383e+18,
             1.39419275e+18, 1.25264986e+18, 1.18187857e+18, 1.16772440e+18,
             1.17126300e+18, 1.13941580e+18, 1.17834000e+18, 1.20664860e+18,
             1.18895580e+18, 1.18895580e+18, 1.21726440e+18, 1.24557296e+18,
             1.22434149e+18, 1.23495719e+18, 1.24203436e+18, 1.22434160e+18,
             1.23495720e+18, 1.21372580e+18, 1.22434160e+18, 1.21018740e+18,
             1.22080300e+18, 1.15357020e+18, 1.19957160e+18, 1.18187880e+18,
             1.19249440e+18, 1.18895579e+18, 1.28449704e+18, 1.27742021e+18,
             1.30218984e+18, 1.30926700e+18, 1.25972716e+18, 1.15003156e+18,
             1.17126296e+18, 1.15710876e+18, 1.10756882e+18, 1.20311006e+18,
             1.29511286e+18, 1.28449726e+18, 1.29157446e+18, 1.44373273e+18,])
        for i in range(20, 40):
            G = tridiagonal(i)
            y = yfull[:i]
            A = LinearAlgebra.linear_least_squares(G, y)[0]
            total = Numeric.add.reduce(y)
            residual = Numeric.absolute(y - Numeric.dot(G, A))
            assert_eq(0.0, Numeric.add.reduce(residual)/total)

    def screwMotion(self):
	import LinearAlgebra
	axis, angle = self.rotation().axisAndAngle()
	d = self.vector*axis
	x = d*axis-self.vector
	r0 = Numeric.dot(LinearAlgebra.generalized_inverse(
	                    self.tensor.array-Numeric.identity(3)), x.array)
	return VectorModule.Vector(r0), axis, angle, d

 def inverse(self):
     """
     @returns: the inverse
     @rtype: L{quaternion}
     """
     import LinearAlgebra
     inverse = LinearAlgebra.inverse(self.asMatrix())
     return quaternion(inverse[:, 0])

 def testEigenvalues(self):
     """
     From bug #545259 -- illegal value to DGEEV routine
     """
     import LinearAlgebra
     a = Numeric.array([[2,3], [4,5]])
     v = LinearAlgebra.eigenvalues(a)
     assert_eq(v, Numeric.array([3.5-math.sqrt(57)/2.,
                                 3.5+math.sqrt(57)/2.]))

 def testHeigenvalues(self):
     """
     From bug #464980.
     """
     import LinearAlgebra
     a = Numeric.array([[1.0, 0.01j], [-0.01j, 1.0]])
     v = LinearAlgebra.eigenvalues(a)
     assert_eq(v, [1.01+0.j, 0.99+0.j])
     Hv = LinearAlgebra.Heigenvalues(a)
     assert_eq(v, [1.01+0.j, 0.99+0.j])

	def calcEigenPairs( self, RtR ):
		"""
		Calculate the corresponding eigenpairs for RtR, and sort them accordingly:
		M{m1 >= m2 >= m3}; set M{v3 = v1 x v2} to ensure a RHS
		(CORRECT)
		
		@postcondition: The eigen pairs are calculated, sorted such that M{m1 >= m2 >= m3} and
		M{v3 = v1 x v2}.
		
		@param RtR: 3x3 Matrix of M{R^t * R}.
		@type  RtR: 3x3 Matrix
		@return : Eigenpairs for the RtR matrix.
		@rtype  : List of stuff
		
		"""
		
		eVal, eVec = LinearAlgebra.eigenvectors(RtR)

		# This is cool.  We sort it using Numeric.sort(eVal)
		# then we reverse it using nifty-crazy ass notation [::-1].
		eVal2 = Numeric.sort(eVal)[::-1]
		eVec2 = [[],[],[]] #Numeric.zeros((3,3), Numeric.Float64)
				
		# Map the vectors to their appropriate owners		
		if ( eVal2[0] == eVal[0]):
			eVec2[0] = eVec[0]
			if ( eVal2[1] == eVal[1] ):
				eVec2[1] = eVec[1]
				eVec2[2] = eVec[2]
			else:
				eVec2[1] = eVec[2]
				eVec2[2] = eVec[1]
		elif( eVal2[0] == eVal[1]):
			eVec2[0] = eVec[1]
			if ( eVal2[1] == eVal[0] ):
				eVec2[1] = eVec[0]
				eVec2[2] = eVec[2]
			else:
				eVec2[1] = eVec[2]
				eVec2[2] = eVec[0]
		elif( eVal2[0] == eVal[2]):
			eVec2[0] = eVec[2]
			if ( eVal2[1] == eVal[1] ):
				eVec2[1] = eVec[1]
				eVec2[2] = eVec[0]
			else:
				eVec2[1] = eVec[0]
				eVec2[2] = eVec[1]

		eVec2[2][0] = eVec2[0][1]*eVec2[1][2] - eVec2[0][2]*eVec2[1][1]
		eVec2[2][1] = eVec2[0][2]*eVec2[1][0] - eVec2[0][0]*eVec2[1][2]
		eVec2[2][2] = eVec2[0][0]*eVec2[1][1] - eVec2[0][1]*eVec2[1][0]
		
		return [eVal2, eVec2]

def eig(object,**kw):
    """We try to adhere to the numpy way of doing things, so we need to do a transpose on the results
    that we get back from Numeric"""
    global _LinearAlgebra, _numpy_linalg
    if _LinearAlgebra:
        eigval,eigvec =  _LinearAlgebra.eigenvectors(object,**kw)
        eigvec = objects.numeric.transpose(eigvec)
        return eigval,eigvec
    elif _numpy_linalg:
        return _numpy_linalg.eig(object,**kw)
    else:
        raise AttributeError("No linear algebra functionality to deal with an eigenvectors.")

 def testLinearLeastSquares(self):
     """
     From bug #503733.
     """
     # XXX not positive on this yet
     import LinearAlgebra
     from RandomArray import seed, random
     seed(7,19)
     (n, m) = (180, 35)
     yp = random((n,m))
     y  = random(n)
     x, residuals, rank, sv = LinearAlgebra.linear_least_squares(yp, y)
     # shouldn't segfault.
     assert rank == m

def fit_polynomial(x_vector, data_vector):
    """Fit a polynomial of degree `degree' at the points 
    in `x_vector' to the `data_vector' by interpolation.
    """

    import Numeric as num
    import LinearAlgebra as la

    degree = len(x_vector)-1

    vdm = vandermonde(x_vector, degree)
    result = la.solve_linear_equations(vdm, num.array(data_vector))
    result = list(result)
    result.reverse()
    return Polynomial(result)