"""
Very simple 2-dimensional Matrix class for 2d matrix operations.
I wrote this code to help with my Linear Algebra homework one night
and then reworked it a few months later to take better advantage
of closures, lambda funcs.

The purpose of this code was to help me understand linear
algebra better and brush up on my python skills... i can't 
vouch for it's completeness or correctness.

It also helps a lot if you're familiar with Python, becuase
a lot of the power in this package lies in the ability to 
paa:ss in lambda expressions or functions to the matrix.

@examples
	# Construct a 3x5 matrix
	A = Matrix2D.create(3,5)
	# Construct a 5x5 identity matrix
	I = Matrix2D.identity(5,5)

	# Multiply A times B
	A = Matrix2D.init_list(2, [1,2, 3,4])
	B = Matrix2D.init_list(2, [5,6, 7,8])
	C = A * B	# equivalent to: C = A.multiply(B)
	C.output()	# [[8, 5], [20, 13]]

	#
	I = B / B	# eqivalent to: I = B * B.inverse()
	
	# Multiply by scalar 2
	f = lambda x: x*2
	A.apply_lambda(f)
	A.output()	# [[1,4], [6,8]]
	
	# Round everything to 2 decimal places
	f = lambda x: round(x, 2)
	A.apply_lambda(f)

	# Blank the diagonal with 0's
	# leave everything else untouched
	def f(row,col, value):
		if (row==col): return 0
	A.apply_func(f)
	A.output()	# [[0,4], [6,0]]

@author Scott Hurring [scott at hurring dot com]
	http://hurring.com/code/python/matrix2d/
@version v0.3 - Feb 08, 2006
"""

class MatrixException(Exception):
	pass
	
class Matrix2D(object):
	"""
	Simple class for representing and working with 2D matrixes.
	"""
	dbg = 0
	A = []
	
	def __init__(self, rows,cols, matrix):
		self.rows = rows
		self.cols = cols
		self.setmatrix(matrix)

	#############################################################
	# Static constructor methods

	def init_value(rows, cols, value=0):
		return Matrix2D.create(rows,cols, value)
	
	def init_list(cols, matrix):
		rows = len(matrix) / cols
		A = []
		for row in range(rows):
			A.insert(row, [])
			for col in range(cols):
				value = matrix[(cols*row)+col]
				A[row].insert(col, value)
		return Matrix2D(rows,cols, A)

	def create(rows,cols, value=0):
		"""
		Create a matrix.
		"""
		A = []
		for row in range(rows):
			A.insert(row, [])
			for col in range(cols):
				A[row].insert(col, value)
		return Matrix2D(rows,cols, A)

	def ones(rows,cols):
		return Matrix2D.create(rows,cols, 1)

   	def zeros(rows,cols):
		return Matrix2D.create(rows,cols, 0)

	def identity(rows,cols):
		M = Matrix2D.create(rows,cols, 0)
		for row in range(M.rows):
			M.set(row,row, 1)
		return M	
		
	init_value = staticmethod(init_value)
	init_zero = init_value
	init_list = staticmethod(init_list)
	create = staticmethod(create)
	ones = staticmethod(ones)
	zeros = staticmethod(zeros)
	identity = staticmethod(identity)

	##################################################################
	# Return new matrices based on the current one

	def clone(self):
		"""
		Return a clone of this matrix.
		"""
		return Matrix2D(self.rows,self.cols, self.getmatrix())

	##################################################################
	# Operate on this matrix's elements
	
	def set_all(self, value):
		"""
		Set all elements of a created matrix to a 'value'.
		"""
		for row in range(self.rows):
			for col in range(self.cols):
				self.A[row][col] = value

	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
		
	##################################################################
	# Operations with another matrix
	# A op B -> C
	
	def __add__(self, B): return self.add(B)
	def add(self, B):
		"""
		Add this matrix to another. C=A+B
		"""
		(m,n) = (B.rows, B.cols)

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

		C = Matrix2D.create(m,n)
		def f(row,col, value):
			return self.get(row,col)+B.get(row,col)	
		C.apply_func(f)
		
		#for (row,col) in C.walk():
		#	C.set(row,col, self.get(row,col)+B.get(row,col))
		return C

	def __mul__(self, B): return self.multiply(B)
	def multiply(self, B):
		"""
		Multiply this matrix against another. C=AB
		"""
		m = self.rows
		(n, l) = (B.rows, B.cols)

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

		C = Matrix2D.create(m,l)
		def f(row,col, value):
			return self.dot(B, row,col)
		C.apply_func(f)
		
		#for (row,col) in C.walk():
		#	value = self.dot(B, row,col)
		#	C.set(row,col, value)
		return C
		
	def __div__(self, B): return self * B.inverse()

	####################################################################
	# Options on this matrix that return a new matrix
	# A.op() -> B
	
	def inverse(self):
		"""
		Compute and Return the Gauss-Jordan inverse of this matrix. A->A^-1
		"""
		self.debug("Inverse", 1)
		B = self.append_identity()
		B.reduce()
		B = B.unshift_identity()
		return B

	def reduce(self):
		"""
		Reduce this matrix to 'reduced echelon' form.  A->R
		"""
		self.debug("Reduce", 1)
		#B = self.clone()
		self.upper_triangular()
		self.pivots_to_1()
		self.lower_triangular()

	def transpose(self):
		"""
		Transpose this matrix.  A->At
		"""
		self.debug("Transpose", 1)
		B = Matrix2D.create(self.cols,self.rows)
		
		def f(row,col, value):
			return self.get(col,row)
		self.apply_func(f)
		
		# replaced with apply_func code
		#for row in range(self.rows):
		#	for col in range(self.cols):
		#		B.set(col,row, self.get(row,col))
		
		self.debug(B.getmatrix(), 2)
		return B
		
	def append_identity(self):
		"""
		Return new matrix with identity appended on end. 
		Used to find inverse
		A -> [A I]
		"""
		self.debug("Append identity", 1)
		B = Matrix2D.create(self.rows,self.cols*2)
		
		def f(row,col, value):
			if col >= self.cols: return (int)(row==col-self.cols)
			else: return self.get(row,col)
		B.apply_func(f)
		
		# replaced with apply_func code
		#for row in range(self.rows):
		#	for col in range(self.cols):
		#		B.set(row,col, self.get(row,col))
		#		val = (int)(col == row)
		#		B.set(row,col+self.cols, val)
		
		self.debug(B, 2)
		return B
		
	def unshift_identity(self):
		"""
		Return this matrix without the identity matrix in front.
		Used after finding inverse.
		[I A] -> A
		"""
		self.debug("Unshift identity", 1)
		cols = self.cols / 2
		B = Matrix2D.create(self.rows,cols)
		I = Matrix2D.create(self.rows,cols)
		for row in range(self.rows):
			for col in range(cols):
				I.set(row,col, self.get(row,col))
				B.set(row,col, self.get(row,col+cols))
		if not I.is_identity():
			raise Exception("Expected identity in front, but found ", I.getmatrix())
		self.debug(B.getmatrix(), 2)
		return B

	##################################################################
	# Functions that mutate this matrix in-place

	def upper_triangular(self):
		"""
		Put this matrix into upper-triangular form.
		Go TOP->BOTTOM, eliminating all numbers UNDER pivots.
		"""
		self.debug("Upper triangular", 1)
		self.debug(self.A, 2)
		(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_row(mult, i, j)
		self.debug(self.getmatrix(), 2)
		
	def lower_triangular(self):
		"""
		Put this matrix into lower-triangular form.
		Go BOTTOM->TOP eliminating all numbers ABOVE pivots.
		"""
		self.debug("Lower triangular")
		(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_row(mult, i, j)
				#self.output()
		self.debug(self.getmatrix(), 2)

	def pivots_to_1(self):
		"""
		Divide entire row by pivot to set pivot=1.
		"""
		self.debug("Pivots to 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)
		self.debug(self.getmatrix(), 2)

	def subtract_row(self, mult, i,j):
		"""
		Subtract (row i)*mult from (row j)
		or: (row j) + (row i)*(-mult).
		"""
		self.debug("Subtract pivot ([%i][%i] mult=%0.4f)" % (i,j,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)

		self.debug(self.getmatrix(), 2)
		
	def apply_func(self, func):
		"""
		Apply a function to every element in this matrix.
		if func returns None, the value is left unchanged
		* func(row,col, value)
		"""
		for (row,col) in self.walk():
			fvalue = func(row,col, self.get(row,col))
			if fvalue != None:
				self.set(row,col, fvalue)
				
	def apply_lambda(self, func):
		"""
		Apply a simple lambda to every element in this matrix.
		"""
		for (row,col) in self.walk():
			self.set(row,col, func(self.get(row,col)))
	
	def walk(self):
		"""
		Walk through all elements.
		"""
		for row in range(self.rows):
			for col in range(self.cols):
				yield(row,col)
		
	####################################################################
	# Questions about the matrix
	
	def is_identity(self):
		"""
		Is this matrix the identity matrix?
		"""
		for (row,col) in self.walk():
			value = round(self.get(row,col), 5)
			if row == col:
				if value != 1:
					return False
			elif value != 0:
				return False
		return True
	
	def round(self, n=2):
		f = lambda x: round(x, n)
		self.apply_lambda(f)
		
	def cleanup(self):
		self.round(2)
		def f(row,col, value):
			# cleanup "-0.00" stuff
			if not value: return abs(value)
			return round(value, 2)
		self.apply_func(f)
		
	##################################################################
	# Debug, output, etc... 
	
	def output(self, text=""):
		"""
		Print the matrix.
		"""
		print "Matrix(%i,%i): %s" % (self.rows,self.cols,text)
		for (row,col) in self.walk():
			if col == 0: print "[",
			print "%10.2f |" % self.get(row,col),
			if col == self.cols-1: print 
			
	def debug(self, text, level=3):
		if self.dbg >= level: print text
#
#






## 	def dot(self, B, i=1, j=1):
## 		"""
## 		Dot product of 1xn and nx1 vectors
## 		"""
## 		At = self.transpose()
## 		C = At.summate(B)
	
	def dot(self, B, i=0,j=0):
		"""
		Dot product of row 'i' of A and column 'j' of B
		Sum (for k=range(B.rows) of: A[i][k]*B[k][j]
		scalar = A(row i) dot B(col j)
		"""
		value = 0
		for k in range(B.rows):
			#print "value +=A[%i][%i] * B[%i][%i]" % (i,k, k,j)
			value += self.get(i,k) * B.get(k,j)
		return value