from numpy import *
import matplotlib.pyplot as plt

# Chapter 1

def trifactor(a, d, c):
    """
    Compute the entries in the LU-factorization of a tridiagonal matrix A using algorithm 1.8. 
    Returns the sequence l on the sub-diagonal of L, and the sequence u on the diagonal of U. 
    
    a: The n-1 entries on the sub-diagonal of A
    d: The n entries on the diagonal of A
    c: The n-1 entries on the super-diagonal of A
    """
    u, l = zeros_like(d), zeros_like(a)
    u[0] = d[0]
    for k in range(len(a)):
        l[k] = a[k]/u[k]
        u[k+1] = d[k+1] - l[k]*c[k]
    return l, u

def trisolve(l, u, c, b):
    """
    Compute the solution to a system on the form A*x=b with A a tridiagonal matrix, based on the LU-factorization of A 
    as given by the function trifactor (using algorithm 1.9). 
    The input vector b will be updated with the solution x (i.e. in-place).
    
    l: The entries on the sub-diagonal of L in the LU-factorization of A
    u: The entries on the diagonal of U in the LU-factorization of A
    c: The entries on the super-diagonal of U in the LU-factorization of A (this equals the entries on the super-diagonal of A)
    b: The right hand side
    """
    n = shape(b)[0]
    for k in range(1, n):
        b[k] -= l[k-1]*b[k - 1]
    b[n-1] /= u[n-1]
    for k in range(n-2,-1,-1):
        b[k] -= c[k]*b[k + 1]
        b[k] /= u[k]

def splineint(a, b, y, mu1, munp1):
    """
    Use algorithm 1.19 to solve the D_2 spline problem.
    Returns the left end points of all sub-inetrvals, and a matrix holding all coefficients in all the piecewise cubic polynomials 
    (one row for each piecewise polynomial)
    
    a: The left end point of the interval
    b: The rigth end point of the interval
    y: The values at the interpolation points
    mu1 and munp1: Initial conditions/moments.
    """
    n = len(y) - 1
    h = (b - a)/float(n)
    c, d = ones(n - 2), 4*ones(n - 1)
    l ,u = trifactor(c, d, c)
    b1 = (6/h**2)*(y[2:(n + 1)] - 2*y[1:n] + y[:(n - 1)])
    b1[0] -= mu1
    b1[n - 2] -= munp1
    trisolve(l, u, c, b1)
    mu2ton = b1
    mu = zeros(n + 1)
    mu[0] = mu1
    mu[1:n] = mu2ton
    mu[n] = munp1
    
    C = array(zeros((n, 4))) # Apply lemma 1.3 to obtain the coefficient matrix
    C[:, 0] = y[:n]
    C[:, 1] = (y[1:(n+1)] - y[:n])/h - h*mu[:n]/3-h*mu[1:(n+1)]/6
    C[:, 2] = mu[:n]/2
    C[:, 3] = (mu[1:(n + 1)] - mu[0:n])/(6*h)
    C = matrix(C)
    return linspace(a, b - h, n), C


def findsubintervals (t ,x):
    """
    Use algorithm 1.21 to return an index set for the points in x, saying which intervals they belong to.
    
    t: the start points of the intervals
    x: The points to find the interval indices for.
    """
    k, m = len(t), len(x)
    if k<2:
        return zeros(m,1)
    else:
        j = concatenate([t, x]).argsort()
        i = nonzero(j >= k)
        arr = arange(0,m)
        arr =  i - arr - 1
        arr = arr[0]
        return arr

def splineval(x,C,X):
    """
    Use algorithm 1.22 to compute and return the values of a cubic spline (x, C) at points X. (x, C) can be obtained with the function splineint. 
     
    x: The start points of the intervals.
    C: A matrix holding all coefficients in all the piecewise cubic polynomials (one row for each piecewise polynomial)
    X: The points we want to evaluate for the spline
    """
    m = len(X)
    i = findsubintervals(x,X) 
    G = zeros(m)
    for j in range(m):
        k = i[j]
        t = X[j] - x[k]
        G[j]=C[k,:]* t**array([[0],[1],[2],[3]])
    return G

# Chapter 2

def rforwardsolve(A, b, d):
    """
    Solve the system A*x = b for a lower triangular d-banded matrix A, using algorithm 2.6 (row-oriented forward substitution)
    The input vector b will be updated with the solution x (i.e. in-place).
    
    A: A lower triangular matrix
    b: The right hand side
    d: The band-width
    """
    n = len(b)
    b[0] /= A[0, 0]
    for k in range(1,n):
        lk = array([0,k-d]).max()
        b[k] = b[k] - dot(A[k, lk:k],b[lk:k])
        b[k] /= A[k, k]
    
def rbacksolve(A, b, d):
    """
    Solve the system A*x = b for an upper triangular d-banded matrix A, using algorithm 2.7 (row-oriented backward substitution)
    The input vector b will be updated with the solution x (i.e. in-place).
    
    A: An upper triangular matrix
    b: The right hand side
    d: The band-width
    """
    n = len(b)
    b[n - 1] /= A[n - 1,n - 1]
    for k in range(n-2,-1,-1):
        uk = array([n, k + d + 1]).min()
        b[k] = b[k] - dot(A[k,(k+1):uk], b[(k+1):uk])
        b[k] /= A[k,k]
    
def cforwardsolve(A, b, d):
    """
    Solve the system A*x = b for a lower triangular d-banded matrix A, using algorithm 2.9 (column-oriented forward substitution)
    The input vector b will be updated with the solution x (i.e. in-place).
    
    A: A lower triangular matrix
    b: The right hand side
    d: The band-width
    """
    A = matrix(A)
    n = len(b)
    for k in range(n-1):
        b[k] /= A[k, k]
        uk = array([n, k + d + 1]).min()
        b[(k+1):uk] -= A[(k+1):uk, k]*b[k]
    b[n - 1] /= A[n - 1,n - 1] 
    
def cbacksolve(A, b, d):
    """
    Solve the system A*x = b for a upper triangular d-banded matrix A, using algorithm 2.10 (column-oriented backward substitution)
    The input vector b will be updated with the solution x (i.e. in-place).
    
    A: An upper triangular matrix
    b: The right hand side
    d: The band-width
    """
    A = matrix(A)
    n = len(b)
    for k in range(n - 1,0,-1):
        b[k] /= A[k, k]
        lk = array([0, k - d]).max()
        b[lk:k] -= (A[lk:k, k]*b[k])
    b[0] /= A[0,0]
    
# Algorithm 2.18
def L1U(A, d):
    """
    Compute the L1U-factorization of the d-banded matrix A.
    Returns a lower triangular matrix L with with ones on the diagonal, and an upper triangular matrix U. 
    Note that, for d=1, the algorithm is equivalent to trisolve.
    
    A: A matrix
    d: The band-width
    """
    n = shape(A)[0]
    L = eye(n)
    U = matrix(zeros((n,n))); U[0,0] = A[0,0]
    for k in range(1,n):
        km = array([0, k - d]).max()
        if km < k:
            L[k, km:k] = A[k, km:k]
            rforwardsolve(U[km:k, km:k].T, L[k, km:k].T, d) # L
        U[km:(k + 1), k] = A[km:(k + 1), k]
        rforwardsolve(L[km:(k + 1), km:(k + 1)], U[km:(k + 1), k], d) # U
    return L, U

def L1Uv2(A, d):
    """
    Compute the L1U-factorization of the d-banded matrix A.
    This is more advanced than L1U in the sense that it avoids allocation of new arrays: The result is stored directly in A (i.e. in-place), i.e.
    overwriting the previous A. The entries in A and U are at the same places in A, and the diagonal with ones in L is not stored in A.
    In order to achieve this, we can't use the function rforwardsolve, as L1U uses it (due to the special form with ones on the diagonal).
    
    A: A matrix
    d: The band-width
    """
    n = shape(A)[0]
    for k in range(1,n):
        km = array([0, k - d]).max() # First index of r we need to update
        for r in range(km, k - 1):
            A[k, r] /= A[r, r]
            uk = array([k, r + d + 1]).min() # last index not included
            A[k, (r + 1):uk] -= A[r, (r + 1):uk]*A[k, r]
        A[k, k - 1] /= A[k - 1,k - 1] 
        for r in range(km, k):
            uk = array([k + 1, r + d + 1]).min() # last index not included
            A[(r + 1):uk, k] -= A[(r + 1):uk, r]*A[r, k]
                   
# Chapter 3

def LDL(A, d):
    """
    Find the LDL factorization of the hermitian matrix A, using algorithm 3.4
    
    A: A hermitian matrix
    d: The band-width
    """
    n = shape(A)[0]
    L = array(eye(n))
    dg = zeros(n)
    dg[0] = A[0, 0]
    for k in range(1, n):
        m = reshape(array(A[:k, k].copy()), k)
        rforwardsolve(L[:k, :k], m, d)
        L[k, :k] = m/dg[:k]
        dg[k] = A[k, k] - dot(L[k, :k], m)
    return L, dg 

def bandcholesky(A, d):
    """
    Find the Cholesky factorization of the positive definite matrix A, using algorithm 3.13
    
    A: A positive definite matrix
    d: The band-width
    """
    L, dg = LDL(A, d)
    return matrix(L)*diag(sqrt(dg))
    
def bandsemicholeskyL(A, d):
    """
    Find the Cholesky factorization of the positive semidefinite matrix A, using algorithm 3.24 
    The input matrix A will be changed by this method
    
    A: A positive semidefinite matrix
    d: The band-width
    """
    n = shape(A)[0]
    for k in range(n):
        if A[k,k] > 0:
            kp=array([n, k + 1 + d]).min();
            A[k,k] = sqrt(A[k,k])
            A[(k+1):kp, k] =A [(k+1):kp, k]/A[k, k]
            for j in range(k+1, kp):
                A[j:kp, j] = A[j:kp, j] - A[j, k]*A[j:kp, k]
        else:
            A[k:kp, k] = 0
    return tril(A)
    
# Chapter 4

def housegen(x): 
    """
    Return u and a in Theorem 4.16, using algorithm 4.17
    
    x: The vector to reflect on ae_1
    """
    a = linalg.norm(x)
    if a == 0:
        u=x; u[0]=sqrt(2); return u, a
    if x[0] == 0:
        r = 1
    else:
        r =x[0]/abs(x[0])
    u = conj(r)*x/a
    u[0]=u[0]+1
    u=u/sqrt(u[0])
    a=-r*a
    return u, a

def housetriang(A, B):
    """
    Return R, C so that R is upper trapezoidal and resulting from a series of Householder transformations to A, and C so that it results from applying the same series of Householder transformations to B.
    This is algorithm 4.21

    A: The matrix to apply Householder transformations to, in order to transform it to upper trapezoidal form.
    B: The matrix to apply the same Householder transformations to. 
    """
    m, n = shape(A); r = shape(B)[1] ; A=hstack([A,B]); 
    minval = array([n, m - 1]).min()
    for k in range(minval):
        v, A[k, k] = housegen(A[k:m, k])
        v = matrix(reshape(v, (m - k, 1)))
        C = A[k:m, (k+1):(n+r)] ; A[k:m, (k + 1):(n + r)] = C - v*(v.T*C)
    R = triu(A[:, :n]); C = A[:, n:(n + r)]
    return R, C

def rothesstri(A, b):
    """
    Solve the system Ax=b where A is in upper Hessenberg form using Givens rotations, using algorithm 4.36
    
    A: A matrix in upper Hessenberg form.
    b: The right hand side.
    """
    n = shape(A)[0]
    A = hstack([A, b])
    for k in range(n-1):
        r = linalg.norm([ A[k , k] , A[k + 1, k] ])
        if r>0:
            c=A[k, k]/r; s=A[k + 1, k]/r
            A[[k, k + 1],(k + 1):(n + 1)]=[[c, s],[-s, c]]*A[[k, k + 1],(k + 1):(n + 1)]
        A[k, k] = r; A[k+1,k] = 0
    z = A[:, n].copy()
    rbacksolve(A[:, :n], z, n)
    return z
      
# testcode

def _test_chap1(): 
    print 'Testing all 5 algorithms from chapter 1 with a spline plot'
    xplot = linspace(0,1,101) 
    xknot = linspace(0,1,11)
    yknot = xknot**4
    xknot, C = splineint(0, 1, yknot, 0, 12) # xknot is reduced with one in size
    G = splineval(xknot, C, xplot)
    plt.plot(xplot,G)
    plt.show()

def _test_chap2():
    print 'Testing all algorithms from chapter 2'
    
    n = 10
    d = 9
    b = random.random((n,1))
    
    print 'Testing rforwardsolve'
    A = matrix(zeros((n, n)))
    for k in range(n):
        for s in range(d+1):
            if k+s < n:
                A[k + s, k] = random.random()
    bcopy = b.copy()
    rforwardsolve(A, bcopy, d)
    xver = linalg.solve(A,b)
    assert abs(xver-bcopy).max() < 1E-4, 'rforwardsolve failed'
    
    print 'Testing cforwardsolve'
    bcopy = b.copy()
    cforwardsolve(A, bcopy, d)
    assert abs(xver-bcopy).max() < 1E-4, 'cforwardsolve failed'
    
    b = random.random((n,1))
    print 'Testing rbacksolve'
    A = zeros((n, n))
    for k in range(n):
        for s in range(d + 1):
            if k+s < n:
                A[k, k + s] = random.random()
    bcopy = b.copy()
    rbacksolve(A, bcopy, d)
    xver = linalg.solve(A,b)
    assert abs(xver-bcopy).max() < 1E-4, 'rbacksolve failed'
    
    print 'Testing cbacksolve'
    bcopy = b.copy()
    cbacksolve(A, bcopy, d)
    assert abs(xver-bcopy).max() < 1E-4, 'cbacksolve failed'
    
    print 'Testing L1U'
    A = matrix(zeros((n, n)))
    for k in range(n):
        for s in range(n):
            if abs(k-s) <= d:
                A[k, s] = random.random()
    x = matrix(random.random((n, 1)))
    L, U = L1U(A, d)
    assert abs(L*U*x -A*x).max() < 1E-4, 'L1U failed'
        
    print 'Testing L1Uv2'
    Acopy = A.copy()
    L1Uv2(Acopy, d)
    for k in range(n):
        assert abs(Acopy[k,k] - U[k, k]) < 1E-4, 'L1Uv2 failed'
        for l in range(k):
            assert abs(Acopy[k,l] - L[k, l]) < 1E-4, 'L1Uv2 failed'
            assert abs(Acopy[l,k] - U[l, k]) < 1E-4, 'L1Uv2 failed'
 
def _test_chap3():
    print 'Testing LDL'
    n = 10
    d = 9
    A = matrix(zeros((n, n)))
    for k in range(n):
        for s in range(k + 1):
            if abs(k-s) <= d:
                A[k, s] = random.random()
                A[s, k] = A[k, s]
    L, dg = LDL(A, d)
    L = matrix(L)
    dg = diag(dg)
    newmatr = L*matrix(dg)*L.T
    assert abs(A - newmatr).max() < 1E-4, 'LDL failed'
    
    print 'Testing bandcholesky'
    A = matrix(zeros((n, n)))
    for k in range(n):
        for s in range(k + 1):
            if abs(k-s) <= d:
                A[s, k] = random.random()
    A = matrix(A)*matrix(A).T
    L = bandcholesky(A, d)
    L = matrix(L)
    newmatr = L*L.T
    assert abs(A - newmatr).max() < 1E-4, 'bandcholesky failed'
    
    print 'Testing bandsemicholeskyL'
    Acopy = A.copy()
    L = bandsemicholeskyL(Acopy, d)
    L = matrix(L)
    newmatr = L*L.T
    assert abs(A - newmatr).max() < 1E-4, 'bandsemicholeskyL failed'
    
def _test_chap4():
    print 'Testing housegen and housetriang'
    A = array([[1,3,1],[1,3,7],[1,-1,-4],[1,-1,2]])
    R, C = housetriang(A,eye(4))
    assert abs(matrix(linalg.inv(C))*R - A).max() < 1E-4, 'housetriang failed'
    
    n = 10
    print 'rothesstri'
    A = matrix(zeros((n, n)))
    for k in range(n):
        for s in range(k + 2):
            if s < n:
                A[s, k] = random.random()
    b = random.random((n,1))
    xexact = linalg.solve(A, b)
    x = rothesstri(A, b)
    assert abs(x - xexact).max() < 1E-4, 'rothesstri failed'
    
if __name__=='__main__':
    _test_chap1()
    _test_chap2()
    _test_chap3()
    _test_chap4()