"""
This code was written to help make sure i'm doing my
Linear Algebra homework correctly.

To the best of my knowledge, it multiplies, adds,
and finds inverses for 2-dimensional matrixes.

Scott Hurring - scott at hurring dot com
"""

class MatrixException(Exception):
	pass

class Matrix2D:
	"""
	Simple class for representing and working with 2D matrixes.
	"""
	matrix = []

	def __init__(self, rows, cols):
		"""
		Create a new matrix
		"""
		self.rows = rows
		self.cols = cols
		self.create()

	def create(self, val=0):
		"""
		Create the matrix.
		"""
		self.A = []
		for row in range(self.rows):
			self.A.insert(row, [])
			for col in range(self.cols):
				self.A[row].insert(col, val)

	def ones(self):
		"""
		Set matrix to all '1's.
		"""
		self.create(1)

	def zeros(self):
		"""
		Set matrix to all '1's.
		"""
		self.create(0)

	def identity(self):
		"""
		Make this matrix an identity matrix
		"""
		if self.rows != self.cols:
			raise MatrixException("Cannot make identity matrix")

		self.zeros()
		for row in range(self.rows):
			self.set(row,row, 1)

	def output(self):
		"""
		Print the matrix.
		"""
		for row in self.A:
			for col in row:
				print "%+1.2f " % (col),
			print

	def getmatrix(self):
		"""
		Get the entire matrix.
		"""
		return self.A

	def setmatrix(self, B):
		"""
		Set the entire matrix in one shot
		"""
		self.A = B

	def get(self, row, col):
		"""
		Return an element from the matrix
		"""
		return self.A[row][col]

	def set(self, row,col, value):
		"""
		Set an element of the matrix.
		"""
		self.A[row][col] = value

	def add(self, B):
		"""
		Add this matrix with another.
		"""
		(m,n) = (B.rows, B.cols)

		if m != self.rows or n != self.cols:
			raise MatrixException("add: Invalid dimensions (%i,%i)" % (m,n))

		matrix = Matrix2D(m, n)
		for row in range(m):
			for col in range(n):
				matrix.set(row,col, self.get(row,col)+B.get(row,col))
		return matrix

	def summate(self, B, i, j, n):
		"""
		Summate the value of A.summate(B, i, j)
		#(AB)ij
		"""
		value = 0
		for k in range(n):
			#print "value +=A[%i][%i] * B[%i][%i]" % (i,k, k,j)
			value += self.get(i,k) * B.get(k,j)
		#print
		return value

	def mult(self, B):
		"""
		Multiply this matrix against another.
		"""
		m = self.rows
		(n, l) = (B.rows, B.cols)

		if n != self.cols:
			raise MatrixException("mult: Invalid dimensions (%i,%i)" % (n,l))

		matrix = Matrix2D(m, l)
		for row in range(m):
			for col in range(l):
				value = self.summate(B, row,col, n)
				matrix.set(row,col, value)
		return matrix


	def append_identity(self):
		"""
		Append the identity matrix to the end of this matrix.
		Before finding inverse.
		"""
		B = Matrix2D(self.rows, self.cols*2)
		for row in range(self.rows):
			for col in range(self.cols):
				B.set(row,col, self.get(row,col))

				if col == row: val = 1
				else: val = 0
				B.set(row,col+self.cols, val)
		return B

	def unshift_identity(self):
		"""
		Remove the identity matrix from the beginning of this matrix.
		After finding inverse.
		"""
		cols = self.cols / 2
		B = Matrix2D(self.rows, cols)
		for row in range(self.rows):
			for col in range(cols):
				B.set(row,col, self.get(row,col+cols))
		return B

	def is_identity(self):
		"""
		Is this matrix the identity matrix?
		"""
		for row in range(self.rows):
			for col in range(self.cols):
				if row == col:
					if self.get(row,col) != 1:
						return False
				elif self.get(row,col) != 0:
					return False
		return True

	def inverse(self):
		"""
		Compute the inverse of this matrix.
		Return result as a new matrix.
		"""
		# Gauss-Jordan inverse computation
		B = self.append_identity()

		print "Upper Triangular"
		B.upper_triangular()
		B.output()

		print "Pivots to 1"
		B.pivots_to_1()
		B.output()

		print "Lower Triangular"
		B.lower_triangular()
		B.output()

		C = B.unshift_identity()

		return C

	def apply_to_elements(self, func):
		"""
		Apply a function to every element in the list.
		"""
		for row in range(self.rows):
			for col in range(self.cols):
				v = func(self.get(row,col))
				#print "[",row,"][",col,"] =", v
				self.set(row,col, v)

	def pivots_to_1(self):
		"""
		Divide entire row by pivot to get pivot = 1.
		"""
		(m, n) = (self.rows, self.cols)
		for i in range(m):
			pivot = self.get(i,i)
			if pivot == 0: break
			for c in range(n):
				v = (1.*self.get(i,c)) / (1.*pivot)
				#print "pivot=",pivot, "v=", 1.0*self.get(i,c), "/", 1.*pivot, "=", v
				self.set(i,c, v)

	def upper_triangular(self):
		"""
		Put this matrix into upper-triangular form.
		Go TOP->BOTTOM, eliminating all numbers UNDER pivots.
		"""
		(m, n) = (self.rows, self.cols)
		# For each pivot
		for i in range(m-1):
			pivot = self.get(i,i)
			# Skip pivots that are "0"
			if pivot == 0: continue

			# For each row underneath this pivot
			for j in range(i+1,m):
				if self.get(j,i) == 0: continue
				mult = -1.0 * (1.*self.get(j,i)) / (1.*pivot)
				# If this is non-zero below a pivot
				self.subtract_pivot(mult, i, j)

	def lower_triangular(self):
		"""
		Put this matrix into lower-triangular form.
		Go BOTTOM->TOP eliminating all numbers ABOVE pivots.
		"""
		(m, n) = (self.rows, self.cols)
		for i in range(m-1, 0, -1):
			pivot = self.get(i,i)
			# Skip pivots that are "0"
			if pivot == 0: continue

			# For each row above this pivot
			for j in range(i-1, -1, -1):
				if self.get(j,i) == 0: continue
				mult = -1.0 * (1.*self.get(j,i)) / (1.*pivot)
				# If this is non-zero above a pivot
				self.subtract_pivot(mult, i, j)
				self.output()

	def subtract_pivot(self, mult, i,j):
		"""
		Subtract (row i)*mult from (row j).
		"""
		print "mult = ", mult

		# Muliplier_row = Pivot_row * mult
		mrow = []
		for c in range(self.cols):
			mrow.insert(c, self.get(i,c)*mult)

		#print "[",i,"][",i,"] [",j,"] row=", self.get(j,i) ,"; mult=",mult, " :: ", mrow

		# Now subtract mrow from A[row+1]
		for c in range(self.cols):
			v = self.get(j,c)+mrow[c]
			#print self.get(j,c), "+", mrow[c], "=", v
			self.set(j,c, v)

#
#