Thema anzeigen - Lauffähiger Python Code für Kurs 1142, KE 4, Fernuni Hagen :: Astronomy and Space Flight, Astronomie und Raumfahrt und dies und das... ::

# simple gauss elimination
# Written by Rainer Gerhards 2009-05-26 * <rgerhards@adiscon.com>
# WARNING: no syntax/semantic checks done so far on input stream!
# Needs the numpy package for matrix operations
#

import sys, operator
import numpy
#import numpy.linalg

def getNewWords():
"""get list of words from a single line of stdin"""
while True:
   line = sys.stdin.readline()
   if line == "":
      return
   if line[0] == '#' or line == "\n":
      continue
   # a bit ugly, but works...
   mylist = []
   for item in line.strip().split():
      if item != '':
         mylist.append(item)
   return mylist

def initMat(m, rows, cols):
m['rows'] = int(rows)
m['cols'] = int(cols)
m['mat'] = numpy.zeros((rows, cols), numpy.float_, 'C')

def fillDataMat(m):
words = getNewWords()
for i in range(m['rows']):
   if words[0] == "enddata":
      return
   for j in range(m['cols']):
      m['mat'][i, j] = float(words[j])
   words = getNewWords()

def matCreate(name):
"""create a new matrix object"""
m = {}
m['name'] = name
words = getNewWords()
while words:
   if words[0] == "dim":
      initMat(m, int(words[1]), int(words[2]))
   elif words[0] == "data":
      fillDataMat(m)
   elif words[0] == "endmatrix":
      return m
   else:
      print "unexpected token in input: '" + words[0] + "'"
   words = getNewWords()

def matPrint(M, name):
"""print matrix """
print "Matrix " + name + "" + str(M.shape) + ":"
print numpy.array_str(M)
print

def matPrintAll():
for name in allMat:
   matPrint(allMat[name], name)

def matDup(name, newname):
""" duplicate matrix """
new = {}
new['name'] = newname
allMat[newname] = new
mat = allMat[name]
M = mat['mat'].copy()
new['mat'] = M
new['rows'] = mat['rows']
new['cols'] = mat['cols']

# helpers mentioned, but not provided, in 1142 script
def swaprows(M, i, k):
if verbose:
   print "swap rows " + str(i) + "," + str(k)
for j in range(M.shape[1]):
   x = M[i, j]
   M[i, j] = M[k, j]
   M[k, j] = x

def swapcolumns(M, index, i, k):
""" swap columns in a matrix """
""" Note that we must also swap variable indexes, as we rename them """
if verbose:
   print "swap columns " + str(i) + "," + str(k)
x = index[k]
index[k] = index[i]
index[i] = x
for j in range(M.shape[0]):
   x = M[j, i]
   M[j, i] = M[j, k]
   M[j, k] = x

def findpivot(M, index, k):
""" find pivot element """
# k is a_kk which would be a pivot element, but is zero
for i in range(k, M.shape[0]): # all remaining rows
   for j in range(k, M.shape[1]): # all remaining cols
      if(M[i, j] != 0):
         return (1, i, j)
# if we reach this point, all elements were zero and there
# is no pivot element
return (0, k, k)

def isSolvable(M):
""" test if linear system of equations is solvable """
m = M.shape[0]
n = M.shape[1] - 1 # result vector b (Ax=b) inside matrix
print m, n
if m > n:
   return False; # irregular!
# check if rank A = rank A|b
i = m - 1;
while i >= 0 and M[i, i] == 0:
   if M[i, n] != 0:
      if verbose:
         print "rank different in " + str(i+1) + ": " + \
            str(M[i,i]) + "/" + str(M[i,n])
      return False; # rank different!
   i = i - 1
return True;
# end helpers mentioned in 1142 script

### functions imported from FeU Hagen Script 1142, Chapter 5.4
# lightly adopted for this framework

def gausselim(A, i, j):
""" gausselimination on the given matrix for a given rowset """
# i: start line, j: column in question
if verbose:
   print "doing gausselim at " + str(i) + "," + str(j)
for k in range(i + 1, A.shape[0]):
   A_kj=A[k,j]
   A_ij=A[i,j]
   Q_kj=-A_kj/A_ij
   A[k,:]=A[k,:]+Q_kj*A[i,:]
if verbose:
   matPrint(A, 'gausselim result')

# not yet done!
def compute_x(C,ind,d):
# C=[A|b], wobei A eine obere Dreiecksmatrix mit m <= n ist.
# Wir loesen das System Ax=b unter Beruecksichtigung der
# Permutation index fuer x.

## added rgerhards ##
solvable = isSolvable(C)
n = C.shape[0]
m = C.shape[1]
x = numpy.zeros((n, 1), numpy.float_, 'C')
## end rgerhards ##

if solvable:
   x[ind[d]]=C[d,n]/C[d,d]
   for i in range(1,d+1):
      x[ind[d-i]]=C[d-i,n]
      for j in range(0,i):
         x[ind[d-i]]=x[ind[d-i]]-C[d-i,d-j]*x[ind[d-j]]
         x[ind[d-i]]=x[ind[d-i]]/C[d-i,d-i]
return solvable, x

def gaussalg(matA, matB):
A = matA['mat']
m = matA['rows']
n = matA['cols']
d=min(m, n)
index=range(0,n) # Permutation der Variablen
C=numpy.concatenate((A,matB['mat']),1) # Verschmelze A,b nebeneinander
for k in range(d):
   if verbose:
      print "pivot element at " + str(k) + " is: " + str(C[k,k])
   if C[k,k]==0:
      pivotfound,i,j = findpivot(C,index,k)
      if pivotfound == 0:
         break
      else:
         if i != k: # vertausche Zeilen i and k
            swaprows(C,i,k)
         if j != k: # vertausche Spalten
            swapcolumns(C,index,j,k) # und Variablen j and k
   gausselim(C,k,k)
matA['mat'] = C
matA['index'] = index
solvable, x = compute_x(C,index,k)
return solvable, x

### end functions imported from FeU Hagen Script 1142, Chapter 5.4

def solve(matA, matB):
""" (try to) solve a system of linear equations """
solvable, x = gaussalg(matA, matB)
print "system is solvable: " + str(solvable)
if solvable:
   matPrint(x, 'result matrix')

# ---------------------------------------------------------------------- #
# main program
print "Gauss 0.1, written by Rainer Gerhards (20090526)"
verbose = False
allMat = {}
words = getNewWords()
while words:
cmd = words[0].lower()
if cmd == "verbose":
      if words[1] == "on":
         verbose = True
      else:
         verbose = False
elif cmd == "matrix":
   allMat[words[1]] = matCreate(words[1])
elif cmd == "printmatrix":
   matPrint(allMat[words[1]]['mat'], words[1])
elif cmd == "printallmatrix":
   matPrintAll()
elif cmd == "duplicatematrix":
   matDup(words[1], words[2])
elif cmd == "issolvable":
   print "matrix " + words[1] + " is solvable: " + \
      str(isSolvable(allMat[words[1]]['mat']))
elif cmd == "gausselim":
   gausselim(allMat[words[1]]['mat'], int(words[2])-1, int(words[3])-1)
elif cmd == "gaussalg":
   gaussalg(allMat[words[1]], allMat[words[2]])
   print "associated index: " + str(allMat[words[1]]['index'])
elif cmd == "solve":
   solve(allMat[words[1]], allMat[words[2]])
elif cmd == "swaprows":
   swaprows(allMat[words[1]]['mat'], int(words[2])-1, int(words[3])-1)
elif cmd == "swapcols":
   index = range(0, allMat[words[1]]['cols'])
   swapcolumns(allMat[words[1]]['mat'], index,
         int(words[2])-1, int(words[3])-1)
   print "associated index: " + str(index)
else:
   print "unexpected token in input: '" + cmd + "'"
words = getNewWords()