Cubic Spline/1002/Tri Diagonal.py

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)