Python LinearAlgebra.inverse方法代碼示例

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
    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)

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
 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"

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
 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))

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
def inversedot_woodbury(m, v):
    a = zeros((n, 11), Float)
    for i in range(n):
        for j in range(max(-7, -i), min(4, n - i)):
            a[i, j + 7] = m[i, i + j]
    print a
    al, indx, d = band.bandec(a, 7, 3)
    VtZ = identity(4, Float)
    Z = zeros((n, 4), Float)
    for i in range(4):
        u = zeros(n, Float)
        for j in range(4):
            u[j] = m[j, n - 4 + i]
        band.banbks(a, 7, 3, al, indx, u)
        for k in range(n):
            Z[k, i] = u[k]
        #Z[:,i] = u
        for j in range(4):
            VtZ[j, i] += u[n - 4 + j]
    print Z
    print VtZ
    H = la.inverse(VtZ)
    print H
    band.banbks(a, 7, 3, al, indx, v)
    return(v - dot(Z, dot(H, v[n - 4:])))

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
 def inverse(self):
     """
     @returns: the inverse
     @rtype: L{quaternion}
     """
     import LinearAlgebra
     inverse = LinearAlgebra.inverse(self.asMatrix())
     return quaternion(inverse[:, 0])

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
def inversedot(m, v):
    return dot(la.inverse(m), v)
    n, nn = m.shape
    if 1:
        for i in range(n):
            sys.stdout.write('% ')
            for j in range(n):
                if m[i, j] > 0: sys.stdout.write('+ ')
                elif m[i, j] < 0: sys.stdout.write('- ')
                else: sys.stdout.write('  ')
            sys.stdout.write('\n')

    cyclic = False
    for i in range(4):
        for j in range(n - 4, n):
            if m[i, j] != 0:
                cyclic = True
    print '% cyclic:', cyclic
    if not cyclic:
        a = zeros((n, 11), Float)
        for i in range(n):
            for j in range(max(-5, -i), min(6, n - i)):
                a[i, j + 5] = m[i, i + j]
        for i in range(n):
            sys.stdout.write('% ')
            for j in range(11):
                if a[i, j] > 0: sys.stdout.write('+ ')
                elif a[i, j] < 0: sys.stdout.write('- ')
                else: sys.stdout.write('  ')
            sys.stdout.write('\n')
        al, indx, d = band.bandec(a, 5, 5)
        print a
        band.banbks(a, 5, 5, al, indx, v)
        return v
    else:
        #return inversedot_woodbury(m, v)
        bign = 3 * n
        a = zeros((bign, 11), Float)
        u = zeros(bign, Float)
        for i in range(bign):
            u[i] = v[i % n]
            for j in range(-7, 4):
                a[i, j + 7] = m[i % n, (i + j + 7 * n) % n]
        #print a
        if 1:
            for i in range(bign):
                sys.stdout.write('% ')
                for j in range(11):
                    if a[i, j] > 0: sys.stdout.write('+ ')
                    elif a[i, j] < 0: sys.stdout.write('- ')
                    else: sys.stdout.write('  ')
                sys.stdout.write('\n')
        #print u
        al, indx, d = band.bandec(a, 5, 5)
        band.banbks(a, 5, 5, al, indx, u)
        #print u
        return u[n + 2: 2 * n + 2]

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
 def __getattr__(self, attr):
     if attr == 'A':
         return squeeze(self.array)
     elif attr == 'T':
         return Matrix(Numeric.transpose(self.array))
     elif attr == 'H':
         if len(self.array.shape) == 1:
             self.array.shape = (1,self.array.shape[0])
         return Matrix(Numeric.conjugate(Numeric.transpose(self.array)))
     elif attr == 'I':
         return Matrix(LinearAlgebra.inverse(self.array))
     elif attr == 'real':
         return Matrix(self.array.real)
     elif attr == 'imag':
         return Matrix(self.array.imag)
     elif attr == 'flat':
         return Matrix(self.array.flat)
     else:
         raise AttributeError, attr + " not found."

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
    def ComputeComplianceMatrix(self):
        print "In BoxMinimizer.ComputeComplianceMatrix"
        bm = self.bulkModulus
        mu = self.shearModulus
        Lambda = bm - (2./3) * mu

        ecMatrix = num.zeros((6,6),num.Float)

        for idx in xrange(0,3):
            ecMatrix[idx, idx] = Lambda + 2*mu

        ecMatrix[0,1] = Lambda
        ecMatrix[1,0] = Lambda
        ecMatrix[0,2] = Lambda
        ecMatrix[2,0] = Lambda
        ecMatrix[1,2] = Lambda
        ecMatrix[2,1] = Lambda
        for idx in xrange(3,6):
            ecMatrix[idx, idx] = mu


        self.compliance = LA.inverse(ecMatrix)

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
    m[4 * i + 5][4 * i + 0] = 0
    m[4 * i + 5][4 * i + 1] = 0
    m[4 * i + 5][4 * i + 2] = 1
    m[4 * i + 5][4 * i + 3] = .5
    m[4 * i + 5][4 * i + 4] = 0
    m[4 * i + 5][4 * i + 5] = 0
    m[4 * i + 5][4 * i + 6] = -1
    m[4 * i + 5][4 * i + 7] = .5

m[n * 4 + 2][2] = 1
m[n * 4 + 3][3] = 1

m[0][n * 4 + 2] = 1
m[1][n * 4 + 3] = 1

def printarr(m):
    for j in range(n * 4 + 4):
        for i in range(n * 4 + 4):
            print '%6.1f' % m[j][i],
        print ''

sys.output_line_width = 160
#print array2string(m, precision = 3)
mi = la.inverse(m)
#printarr(mi)
print ''
for j in range(n + 1):
    for k in range(4):
        print '%7.2f' % mi[j * 4 + k][(n / 2) * 4 + 2],
    print ''

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
def polyfitw(x, y, w, ndegree, return_fit=0):
   """
   Performs a weighted least-squares polynomial fit with optional error estimates.

   Inputs:
      x: 
         The independent variable vector.

      y: 
         The dependent variable vector.  This vector should be the same 
         length as X.

      w: 
         The vector of weights.  This vector should be same length as 
         X and Y.

      ndegree: 
         The degree of polynomial to fit.

   Outputs:
      If return_fit==0 (the default) then polyfitw returns only C, a vector of 
      coefficients of length ndegree+1.
      If return_fit!=0 then polyfitw returns a tuple (c, yfit, yband, sigma, a)
         yfit:  
            The vector of calculated Y's.  Has an error of + or - Yband.

         yband: 
            Error estimate for each point = 1 sigma.

         sigma: 
            The standard deviation in Y units.

         a: 
            Correlation matrix of the coefficients.

   Written by:   George Lawrence, LASP, University of Colorado,
                 December, 1981 in IDL.
                 Weights added, April, 1987,  G. Lawrence
                 Fixed bug with checking number of params, November, 1998, 
                 Mark Rivers.  
                 Python version, May 2002, Mark Rivers
   """
   n = min(len(x), len(y)) # size = smaller of x,y
   m = ndegree + 1         # number of elements in coeff vector
   a = Numeric.zeros((m,m),Numeric.Float)  # least square matrix, weighted matrix
   b = Numeric.zeros(m,Numeric.Float)    # will contain sum w*y*x^j
   z = Numeric.ones(n,Numeric.Float)     # basis vector for constant term

   a[0,0] = Numeric.sum(w)
   b[0] = Numeric.sum(w*y)

   for p in range(1, 2*ndegree+1):     # power loop
      z = z*x   # z is now x^p
      if (p < m):  b[p] = Numeric.sum(w*y*z)   # b is sum w*y*x^j
      sum = Numeric.sum(w*z)
      for j in range(max(0,(p-ndegree)), min(ndegree,p)+1):
         a[j,p-j] = sum

   a = LinearAlgebra.inverse(a)
   c = Numeric.matrixmultiply(b, a)
   if (return_fit == 0):
      return c     # exit if only fit coefficients are wanted

   # compute optional output parameters.
   yfit = Numeric.zeros(n,Numeric.Float)+c[0]   # one-sigma error estimates, init
   for k in range(1, ndegree +1):
      yfit = yfit + c[k]*(x**k)  # sum basis vectors
   var = Numeric.sum((yfit-y)**2 )/(n-m)  # variance estimate, unbiased
   sigma = Numeric.sqrt(var)
   yband = Numeric.zeros(n,Numeric.Float) + a[0,0]
   z = Numeric.ones(n,Numeric.Float)
   for p in range(1,2*ndegree+1):     # compute correlated error estimates on y
      z = z*x		# z is now x^p
      sum = 0.
      for j in range(max(0, (p - ndegree)), min(ndegree, p)+1):
         sum = sum + a[j,p-j]
      yband = yband + sum * z      # add in all the error sources
   yband = yband*var
   yband = Numeric.sqrt(yband)
   return c, yfit, yband, sigma, a

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
cell, atoms = get_info(filename)

print "Old cell:", cell
print toparams(cell)

#
# The transformation matrix is hardwired in this version.
#
matrix = Numeric.array([[ 0., -1., 1.], [ 0., 0., 1.], [ -2., 2. , -1.]])
emat = Numeric.array([[ 0., -1., 1., 0.28488], 
                      [ 0., 0., 1.,0.12343], 
                      [ -2., 2. , -1., 0.0],
                      [ 0.0, 0.0, 0.0, 1.0] ])

invemat = LinearAlgebra.inverse(emat)

# Extend matrix to deal with translations, in the crystallographic way
#
print atoms
natoms = atoms.shape[0]
ones= 1 + Numeric.zeros((natoms,1))
print ones.shape
extatom = Numeric.concatenate((atoms,ones),axis=1)
print extatom

print "---------------"

newcell = Numeric.dot(matrix,cell)
print newcell
print toparams(newcell)

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
 def inverse(self):
     "Returns the inverse."
     import LinearAlgebra
     inverse = LinearAlgebra.inverse(self.asMatrix())
     return Quaternion(inverse[:, 0])

# 需要導入模塊: import LinearAlgebra [as 別名]
# 或者: from LinearAlgebra import inverse [as 別名]
def solve(path):
    if path[0][2] == '{': closed = 0
    else: closed = 1
    dxs = []
    dys = []
    chords = []
    for i in range(len(path) - 1):
        dxs.append(path[i + 1][0] - path[i][0])
        dys.append(path[i + 1][1] - path[i][1])
        chords.append(hypot(dxs[-1], dys[-1]))
    nominal_th = []
    nominal_k = []
    if not closed:
        nominal_th.append(atan2(dys[0], dxs[0]))
        nominal_k.append(0)
    for i in range(1 - closed, len(path) - 1 + closed):
        x0, y0, t0 = path[(i + len(path) - 1) % len(path)]
        x1, y1, t1 = path[i]
        x2, y2, t2 = path[(i + 1) % len(path)]
        dx = float(x2 - x0)
        dy = float(y2 - y0)
        ir2 = dx * dx + dy * dy
        x = ((x1 - x0) * dx + (y1 - y0) * dy) / ir2
        y = ((y1 - y0) * dx - (x1 - x0) * dy) / ir2
        th = fit_arc(x, y) + atan2(dy, dx)
        bend_angle = mod_2pi(atan2(y2 - y1, x2 - x1) - atan2(y1 - y0, x1 - x0))
        k = 2 * bend_angle/(hypot(y2 - y1, x2 - x1) + hypot(y1 - y0, x1 - x0))
        print '% bend angle', bend_angle, 'k', k
        if t1 == ']':
            th = atan2(y1 - y0, x1 - x0)
            k = 0
        elif t1 == '[':
            th = atan2(y2 - y1, x2 - x1)
            k = 0
        nominal_th.append(th)
        nominal_k.append(k)
    if not closed:
        nominal_th.append(atan2(dys[-1], dxs[-1]))
        nominal_k.append(0)
    print '%', nominal_th
    print '0 0 1 setrgbcolor .5 setlinewidth'
    plot_path(path, nominal_th, nominal_k)
    plot_ks(path, nominal_th, nominal_k)
    th = nominal_th[:]
    k = nominal_k[:]
    n = 8
    for i in range(n):
        ev = make_error_vec(path, th, k)
        m = make_matrix(path, th, k)
        #print m
        #print 'inverse:'
        #print la.inverse(m)
        v = dot(la.inverse(m), ev)
        #print v
        for j in range(len(path)):
            th[j] += 1. * v[2 * j]
            k[j] -= 1. * .5 * v[2 * j + 1]
        if i == n - 1:
            print '0 0 0 setrgbcolor 1 setlinewidth'
        elif i == 0:
            print '1 0 0 setrgbcolor'
        elif i == 1:
            print '0 0.5 0 setrgbcolor'
        elif i == 2:
            print '0.3 0.3 0.3 setrgbcolor'
        plot_path(path, th, k)
        plot_ks(path, th, k)
    print '% th:', th
    print '% k:', k