#format python

from LuDecomposition import *
from Numeric import *
from ArrayPrinter import *
import pdb

def getMatrixY(aMatrixLower, b):
	matrixY = makeEmptyMatrix(len(b),1)
	for n in range(len(b)):
		matrixY[n][0] = float(b[n][0] - _minusForGetMatrixY(n, aMatrixLower, matrixY)) / float(aMatrixLower[n][n])
	return matrixY

def getMatrixX(aMatrixUpper, y):
	limitedMatrix = len(y)
	matrixX = makeEmptyMatrix(limitedMatrix,1)
	for n in range(limitedMatrix-1, -1,-1):
		matrixX[n][0] = float(y[n][0] - _minusForGetMatrixX(n, aMatrixUpper, matrixX)) / float(aMatrixUpper[n][n])
		#print "x[%d]: y[%d][0] - minus:[%f] / u[%d][%d]:%f : %f"% (n,n,_minusForGetMatrixX(n, aMatrixUpper, matrixX),n,n, aMatrixUpper[n][n], matrixX[n][0])
	return matrixX

def _minusForGetMatrixX(n, aUpperMatrix, aMatrixX):
	totalMinus = 0
	limitedMatrix = len(aMatrixX)
	for t in range(n+1,limitedMatrix):
		totalMinus += aUpperMatrix[n][t] * aMatrixX[t][0]
	return totalMinus

def _minusForGetMatrixY(n, aLowerMatrix, aMatrixY):
	totalMinus = 0
	for t in range(n):
		totalMinus += aLowerMatrix[n][t] * aMatrixY[t][0]
	return totalMinus
	

def makeEmptyMatrix(aRow, aCol):
	emptyMatrix = []
	for i in range(0,aRow):
		emptyRow = []
		for j in range(0,aCol):
			emptyRow.append(0.0)
		emptyMatrix.append(emptyRow)
	return emptyMatrix

def prettyPrintMatrix(aMatrix):
	print array2string(array(aMatrix))
def printSeparator():
	print "*-" * 35

if __name__=="__main__":
	a =	[[4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
			[1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
			[0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
			[0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0],
			[0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0],
			[0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0],
			[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0],
			[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0],
			[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0]]
	b = [[5.0], [6.0], [6.0], [6.0], [6.0], [6.0], [6.0], [6.0], [5.0]]
	l, u = LuDecomposition(a).perform()
	printSeparator()
	print "LuDecomposition - Matrix L:"
	prettyPrintMatrix(l)
	#pdb.set_trace()
	printSeparator()
	matrixY = getMatrixY(l, b)
	print "Calculated - Matrix Y"
	prettyPrintMatrix(matrixY)
	printSeparator()
	print "LuDecomposition - Matrix U:"
	prettyPrintMatrix(u)
	printSeparator()
	print "Final - Matrix X"
	matrixX = getMatrixX(u, matrixY)
	prettyPrintMatrix(matrixX)