#! /usr/bin/env python
"""Perform some linear algebra and vector operations

Date: November 14, 2008

Author: Bob Wimmer

"""

#from math import *   #This imports all math functions and constants
#from scipy.integrate import odeint
#from numpy import *
#from numpy.linalg import inv,solve
#import Gnuplot, Gnuplot.funcutils
#import sys, time

from numpy import matrix
from scipy.linalg import inv, det

from numpy import allclose, arange, eye, linalg, ones
from scipy import linsolve, sparse


def main():
    A=matrix([[1,1,1],[4,4,3],[7,8,5]]) # 3 lines 3 rows
    b = matrix([1,2,1]).transpose()     # 3 lines 1 rows.
    print det(A)
    print inv(A)*b


    #Now let's have a look at large sparse matrices
    #It's very ineffint to use ordinary inversion techniques for such matrices
    A = eye(10)
    print A
    Asp = sparse.lil_matrix(A)
    print Asp
    #The lil format is LIst of Lists
    Asp2 = sparse.lil_matrix((10,10))
    Asp2.setdiag(ones(10))
    print Asp2

    #Now let's time both operations for large matrices
    from time import time
    t = time(); A = eye(1000); t1=time()-t; print t1
    t=time(); Asp  = sparse.lil_matrix(A);  t2=time()-t; print t2
    #so conversion to sparse takes some time

    t = time(); Asp2 = sparse.lil_matrix((1000,1000)); t3=time()-t; print t3
    t = time(); Asp2.setdiag(ones(1000)); t4=time()-t; print t4
    #Wow! that's a lot faster!
    #ratio is
    print "The second method is " + str((t1+t2)/(t3+t4)) + " times faster!"

    #Now use this to solve a trivial system of equations, A*x = b
    b = arange(1,1001)
    t = time(); x = linalg.solve(A,b); t1 = time()-t; print t1
    t = time(); xsp = linsolve.spsolve(Asp2,b); t2 = time()-t; print t2
    #and check whether x and b are close (they should be)
    print allclose(x,b)
    print allclose(xsp,b)
  






    
    #To perform inversions, etc., it is more efficient to convert lil to another format
    
    



main()