Numerical Analysis in Python

numerical.py

import numpy as np

# numpy.dot
# For 2-D arrays it is equivalent to matrix multiplication, and for
# 1-D arrays to inner product of vectors (without complex conjugation).
# Ex: np.dot([[1, 0], [0, 1]], [[4, 1], [2, 2]])

# numpy.vdot
# The vdot(a, b) function handles complex numbers differently than dot(a, b).
# If the first argument is complex the complex conjugate of the first argument is
# used for the calculation of the dot product. Does not perform a matrix product!
# Ex: np.vdot([1+2j,3+4j], [5+6j,7+8j])
# Ex: np.mat([1+2j,3+4j]) * np.map([[5+6j],[7+8j]])
# numpy.tensordot
# Compute tensor dot product along specified axes for arrays >= 1-D.
# Ex:
# >>> a = np.arange(60.).reshape(3,4,5)
# >>> b = np.arange(24.).reshape(4,3,2)
# >>> c = np.tensordot(a,b, axes=([1,0],[0,1]))

# np.cross
# Return the cross product of two (arrays of) vectors.
# Ex:
# >>> np.cross([1, 2, 3], [4, 5, 6])

# numpy.meshgrid
# Return coordinate matrices from coordinate vectors.
# >>> nx, ny = (3, 2)
# >>> x = np.linspace(0, 1, nx)
# >>> y = np.linspace(0, 1, ny)
# >>> xv, yv = meshgrid(x, y)
# >>> xv, yv = meshgrid(x, y, sparse=True)

# np.array (numpy.ndarray)
# An array object represents a multidimensional, homogeneous array of fixed-size items.
# Ex:
# >>> x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
# >>> x[1:7:2]

# numpy.allclose
# Returns True if two arrays are element-wise equal within a tolerance.

# np.matrix
# The matrix objects are a subclass of the numpy arrays (ndarray).

# x.dtype
# The ype of elements

# numpy.mat(data, dtype=None)
# Unlike matrix, asmatrix does not make a copy if the input is already a matrix or an ndarray. Equivalent to matrix(data, copy=False).
# >>> x = np.array([[1, 2], [3, 4]])
# >>> m.dtype
# >>> m = np.asmatrix(x)
# >>> m.dtype

# Insertint data into arrays:
# fill_diagonal(a, val[, wrap])		Fill the main diagonal
# place(arr, mask, vals)	Change elements of an array based on conditional and input values.

# Iterating over arrays:
# nditer			Efficient multi-dimensional iterator object to iterate over arrays.
# ndenumerate(arr)	Multidimensional index iterator.
# ndindex(*shape)	An N-dimensional iterator object to index arrays.
# flatiter			Flat iterator object to iterate over arrays.

# numpy.linalg.svd
# Singular Value Decomposition.
# The decomposition is performed using LAPACK routine _gesdd
# >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
# Reconstruction based on full SVD:
# >>> U, s, V = np.linalg.svd(a, full_matrices=True)
# Reconstruction based on reduced SVD:
# >>> U, s, V = np.linalg.svd(a, full_matrices=False)

# numpy.empty
# Return a new array of given shape and type, without initializing entries.
# empty, unlike zeros, does not set the array values to zero, and may therefore
# be marginally faster.
# Ex:
# >>> np.empty([2, 2])
# >>> np.empty([2, 2], dtype=int)

# numpy.linalg.norm
# Matrix or vector norm.
# Ex:
# >>> from numpy import linalg as la
# >>> la.norm()

# numpy.seterr
# Set how floating-point errors are handled.
# >>> old_settings = np.seterr(all='print')
# >>> np.geterr()
# Ex:
# >>> np.int16(32000) * np.int16(3)
# >>> m = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
# >>> n = m.reshape((3, 3))
# >>> la.norm(b, 'fro')
# >>> la.norm(a, np.inf)
# >>> la.norm(a, 1)
# >>> la.norm(a, 2)

from numpy import linalg as la
from pprint import pprint
import json

class Householder:
    """QR factorization of a m x n matrix
    with Householder transformations"""

    def __init__(self, M):
        self.M = M

    def QR():
        return M

def householder(a):
    from numpy import dot,diagonal,outer,identity
    from math import sqrt
    n = len(a)
    for k in range(n-2):
        u = a[k+1:n,k]
        uMag = sqrt(dot(u,u))
        if u[0] < 0.0: uMag = -uMag
        u[0] = u[0] + uMag
        h = dot(u,u)/2.0
        v = dot(a[k+1:n,k+1:n],u)/h
        g = dot(u,v)/(2.0*h)
        v = v - g*u
        a[k+1:n,k+1:n] = a[k+1:n,k+1:n] - outer(v,u) \
                         - outer(u,v)
        a[k,k+1] = -uMag
    return diagonal(a),diagonal(a,1)

def computeP(a):
    from numpy import dot,diagonal,outer,identity
    from math import sqrt
    n = len(a)
    p = identity(n)*1.0
    for k in range(n-2):
        u = a[k+1:n,k]
        h = dot(u,u)/2.0
        v = dot(p[1:n,k+1:n],u)/h
        p[1:n,k+1:n] = p[1:n,k+1:n] - outer(v,u)
    return p

class Orthogonalization:
    """Functionality related to orthogonalization of matrices"""
    def __init__(self, M):
    	self.M = M

    def gr1(self):
    	return True

    def gramschmidt(self):
    	return self.M

    def householder(self,x):
        x[0] = x[0] + np.sign(x[0]) * la.norm(x)
        print "x"
        print x
        x /= la.norm(x)
        return x

    def qr(self):
        R = np.matrix(self.M).copy()
        Qt = np.eye(R.shape[0])

        return Qt.T, R

    def qr2(self):
        m, n = self.M.shape
        Q = np.eye(m)
        for i in range(n - (m == n)):
            H = np.eye(m)
            H[i:, i:] = self.make_householder(self.M[i:, i])
            Q = np.dot(Q, H)
            self.M = np.dot(H, self.M)
        return Q, self.M

    def make_householder(self,a):
        v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
        v[0] = 1
        H = np.eye(a.shape[0])
        H -= (2 / np.dot(v, v)) * np.dot(v[:, None], v[None, :])
        return H

class Matrix:
    def __init__(self, rows, cols):
        #self.M = np.zeros(shape=(rows,cols))
        #self.M = np.arange(9).reshape((rows,cols))
        self.M = np.random.rand(rows,cols)
        #self.M = np.eye(rows)
        self.propsdict = {'size': {'rows': rows, 'cols': cols}}
#    def __init__(self, M):
 #       self.M = np.array(M)
#        self.propsdict = {'Name': 'Zara', 'Age': 7, 'Class': 'First'}

    def printme(self):
        print self.M

    def printme2(self):
        # Return an iterator yielding pairs of
        # array coordinates and values.
        for index, x in np.ndenumerate(self.M):
            print index, x

    def getdiag(self):
        # Let the row count be the size of index vector

        # Note: there is also a np.diag(), ex:
        # np.diag(self.M)

        # Better (using numpy.flatiter):
        # diag_indices = range(0,self.M.size,self.M.shape[1]+1)
        # a.flat(diag_indices)

        diag_indices = np.diag_indices(self.M.shape[0])
        return self.M[diag_indices]

    def triliu(self):
        # Return the indices for the lower-triangle of an (n, m) array.
        il1 = np.tril_indices(self.M.shape[0])
        # Return the indices for the upper-triangle of an (n, m) array.
        iu1 = np.triu_indices(self.M.shape[0])
        # Take elements from an array along an axis.
        indices = [0, 1, 4]
        np.take(a, indices)

        # Column wise iteration using self.M.flat
        # Ex: flat_itr = self.M.flat
        # This is an row wise iterator

        # From documentation:
        # Iteration is done in row-major, C-style order
        # (the last index varying the fastest).
        # The iterator can also be indexed using basic
        # slicing or advanced indexing.

        flat_itr = self.M.flat
        for item in flat_itr:
                print item



    def info(self):
        rows, cols = self.M.shape
        self.propsdict["det"] = la.det(self.M)
        self.propsdict["norms"] = {'max(1)': la.norm(self.M,1),
#                                   'min(-1)': la.norm(self.M,-1),
                                   'euclid(2)': la.norm(self.M,2),
                                   'frob': la.norm(self.M,'fro'),
                                   'inf': la.norm(self.M,np.inf)}
        self.propsdict["cond"] = la.cond(self.M)
        self.propsdict["orthogonality"] = la.norm(self.M.T*self.M - np.eye(max(rows, cols)),2)
        U, s, V = la.svd(self.M)
        self.propsdict["singular_vals"] = str(s)
        w, v = la.eig(self.M)
        self.propsdict["eigenvalues"] = str(w)
        self.propsdict["rank"] = la.matrix_rank(self.M)
        print json.dumps(self.propsdict, sort_keys=True, indent=4, separators=(',', ': '))

def testing( matrix ) :
	a = np.arange(6).reshape(2,3)
	#for x in np.nditer(a, op_flags=['readwrite']):

from math import sqrt
import logging

class HouseholderTransform:
    def __init__(self):
        logging.debug( "New instance of HouseholderTransform" )

class GivensRotation:
    def __init__(self):
        logging.debug( "New instance of GivensRotation" )

# Gram-Schmidt
def gramschmidt(M):
    V = M
    rows, cols = M.shape

    for j in range(cols):
       for i in range(rows):
          D = sum([V[p][i]*V[p][j] for p in range(rows)])

          for p in range(rows):
            V[p][j] -= (D * V[p][i])

       invnorm = 1.0 / sqrt(sum([(V[p][j])**2 for p in range(rows)]))

       for p in range(rows):
           V[p][j] *= invnorm

    return V

def gs(X, row_vecs = False, norm = True):
    if not row_vecs:
        X = X.T
    Y = X[0:1,:].copy()
    for i in range(1, X.shape[0]):
        proj = np.diag((X[i,:].dot(Y.T)/np.linalg.norm(Y,axis=1)**2).flat).dot(Y)
        Y = np.vstack((Y, X[i,:] - proj.sum(0)))
    if norm:
        Y = np.diag(1/np.linalg.norm(Y,axis=1)).dot(Y)
    if row_vecs:
        return Y
    else:
        return Y.T

class TestSuit:
    """Test suit with test cases"""
    def __init__(self):
    	return True
    def Test01(self):
        return 'hello world'

if __name__ == "__main__":
    #o = Orthogonalization([[12.0,-51,4],[6,167,-68],[-4,24,-41]])
    #print np.matrix(o.M)

    #Q,R = o.qr()
    #print Q

    X, Y = householder(np.matrix([[12.01,1,4],[6,17,-8]]))
    print X
    print Y

    logging.basicConfig(level=logging.DEBUG)

    # We don't need more right now
    np.set_printoptions(precision=4)
    json.encoder.FLOAT_REPR = lambda f: ("%.4f" % f)

    #m = Matrix([[12.0001,-51,4],[6,167,-68],[-4,24,-41]])
    m = Matrix(3,3)
    m.info()
    m.printme2()
    m.printme()
    print m.getdiag()
    ht = HouseholderTransform()
    gr = GivensRotation()

    if 1==2:
        h = Householder([[1,2,3],[0,1,2],[0,2,3]])
        print h.M

        no = np.matrix([[0.997,0,0],[0,1.007,3],[0.021,0,0.989]])

#    cm = np.array([1, 2, 3], dtype=complex)
#    print cm
        print la.norm(h.M,2)
#    print la.norm(cm,2)
    #### ||A^TA - I|| -- Orthogonality of a matrix
        diff = (no.T *  no - np.eye(3))

        print "----"
        print no
#    print no.T
 #   print diff
        print la.norm(diff,2) ## euclidian/2-norm
        print la.norm(diff,1) ## max
        print la.norm(diff,np.inf)
        print la.norm(diff,'fro')
        U, s, V = np.linalg.svd(no, full_matrices=True)
#    print U
        print s
        print np.absolute(s[:1]*s[:1] - 1)
        print "----"
        print la.cond(no,2) ### == quotient btw max and min eigenvalue of A^TA
        print np.absolute(s[0]/s[2])
    #print np.absolute(s*s - np.eye(3))
#    print V