Prime Python Philanthropy

With each passing day, I realize that I love programming more an more. Part of it comes from the fact that every line of code is a little puzzle: how it interacts with other lines, itself, variables, functions elsewhere in the code, etc. Getting everything to work together is such an accomplishment, something akin to (at least for me) finishing a race, winning a game, or doing well on an exam. I take almost any opportunity to program, even if it’s just a few short lines to decide what to eat for dinner.

Recently, I’ve been helping a few of my friends with their Python homework for CSE 231. I never took the course, my roommate took it last semester, but we both just like flexing our programming muscles every now and again. We are both in a programming course right now (him in CSE 232, me in MTH 451), but just going back and writing little Python programs is always fun. Today, I wrote a short bit of code that added two binary numbers via string manipulation: completely unnecessary in the real world, but fun to program in its own right.

Of course, my own programming homework is always more difficult than what the CSE 231 homework is. This week, for instance, we needed to write three functions to solve a system of linear equations using three different Gaussian elimination schemes (no pivoting, partial pivoting, and scaled partial pivoting). My first stab at it did not go well, since I started with the actual numerical solving schemes than tackling the machinery underneath. After fusing around for a while, I scrapped it and started fresh, first writing basic matrix manipulation functions (switch rows, multiply row by a constant, etc.), then gradually moved up to actually solving the equations.

# Matrix manipulation functions
def exchange_rows(A, i, j):
  A[i], A[j] = A[j], A[i]
  return A

def multiply_row(A, i, mult):
  for q in xrange(len(A[i])):
    A[i][q] = A[i][q]*mult
  return A

def add_rows(A, i, j):
  for q in xrange(len(A[i])):
    A[j][q] = A[j][q] + A[i][q]
  return A

# Pivot Checking
def check_pivot(A, i):
  if A[i][i] == 0:
    A = exchange_rows(A, i, i+1)
  else:
    pass
  return A

# Generic entry elimination functions
def forward_substitution(A, i):
  A = multiply_row(A, i, 1/float(A[i][i]))
  for q in xrange(i+1, len(A)):
    if A[q][i] != 0:
      mult = -1.0*float(A[q][i])
      A = multiply_row(A, i, mult)
      A = add_rows(A, i, q)
      A = multiply_row(A, i, 1.0/mult)
    else:
      pass
  return A

def backward_substitution(A, i):
  for q in xrange(i, 0, -1):
    if A[q-1][i] != 0:
      mult = -1.0*float(A[q-1][i])
      A = multiply_row(A, i, mult)
      A = add_rows(A, i, q-1)
      A = multiply_row(A, i, 1.0/mult)
    else:
      pass
  return A

# Three solver methods
def gauss_np(A):
  for index in xrange(len(A)):
    A = check_pivot(A, index)
    A = forward_substitution(A, index)
  for index in xrange(len(A)-1, -1, -1):
    A = backward_substitution(A, index)
  return A

def gauss_pp(A):
  startrow = 0
  for index in xrange(startrow, len(A)):
    maxindex = 0
    maxvalue = 1
    for i in xrange(index, len(A)):
      if A[i][index] > maxvalue:
        maxindex = i
        maxvalue = A[i][index]
    A = exchange_rows(A, index, maxindex)
    A = forward_substitution(A, index)
    startrow += 1
  for index in xrange(len(A)-1, -1, -1):
    A = backward_substitution(A, index)
  return A

def gauss_spp(A):
  list = []
  for j in xrange(len(A)):
    maxvalue = -1000
    for i in xrange(len(A[0])-1):
      if A[j][i] > maxvalue:
        maxvalue = A[j][i]
    list.append(maxvalue)
  for index in xrange(len(A)):
    temp = list
    for i in xrange(len(temp)):
      temp[i] = abs(A[i][index]) / float(temp[i])
    max = 0
    ind = i
    for i in xrange(index, len(temp)):
      if temp[i] > max:
        max = temp[i]
        ind = i
    temp[index], temp[i] = temp[i], temp[index]
    A = exchange_rows(A, index, i)
    A = forward_substitution(A, index)
  for index in xrange(len(A)-1, -1, -1):
    A = backward_substitution(A, index)
  return A

# Main function commands
def grab_input():
  A = input("Input Matrix A: ")
  b = input("Input Vector b: ")
  option = raw_input("(1-Gaussian, no pivot : 2-Gaussian, partial pivot : Gaussian, scaled partial pivot)\nSelect Solving Method: ")
  if len(b) != len(A[0]):
    print "Matrix/Vector dimensions mismatched."
    A = "1"
  return A, b, option

def prep_matrix(A, b):
  for q in xrange(len(A)):
    A[q].append(b[q])
  return A

def solve_function(A, option):
  if option == "1":
    A = gauss_np(A)
  elif option == "2":
    A = gauss_pp(A)
  elif option == "3":
    A = gauss_spp(A)
  return A

def grab_solution(A):
  x = []
  for q in xrange(len(A)):
    x.append(A[q][-1])
  return x

def main():
  A, b, option = grab_input()
  A = prep_matrix(A, b)
  A = solve_function(A, option)
  x = grab_solution(A)
  print x

if __name__ == "__main__":
  main()

For the solver functions, I wrote gauss_np(), then gauss_pp(), then gauss_spp(). You can kind of tell that order since the latter ones are way less elegant than the rest of the code. I guess that happens after fusing with the program for a few hours straight, but at the same time I enjoyed trying to track down those little bugs and index errors and what-not that crept in while typing the commands.

It’s a scavenger hunt, academic contest, tabletop puzzle, and basketball game all rolled into one!

Advertisements

Comments are closed.