beyn.py

import solver as Solver
import numpy  as np
import numpy.matlib

def simple(operator, origin, radius,
           nodes=100, maxIt=10, lStart=1, lStep=1, rankTol=1e-4, verbose=True):
    """Solves an eigenvalue problem using Beyn's algorithm (simple version)

    Keyword arguments:
    operator -- the solver defining the operator to use
    origin -- the origin (in the complex plane) of the above circular contour
    radius -- the radius of the circular contour used to search the eigenvalues
    nodes -- the number of nodes for the trapezoidal integration rule (optional)
    lStart -- the number of columns used for A0 when algorithm starts (optional)
    lStep -- the step used for increasing the number of columns of A0 (optional)
    rankTol -- the tolerance for the rank test
    verbose -- should I be verbose? (optional)

    Returns the computed eigenvalues
    """

    # Display the parameter used
    if(verbose): display(nodes, maxIt, lStart, lStep, rankTol, origin, radius)

    # Initialise A0 search
    myPath = path(nodes, origin, radius)
    hasK   = False
    it     = 0
    m      = operator.size()
    l      = lStart
    k      = -1

    # Search A0
    if(verbose): print "Searching A0..."
    while(not hasK and it != maxIt):
        if(verbose): print " # Iteration: " + str(it+1)

        vHat = randomMatrix(m ,l)                       # Take a random VHat
        A0   = integrate(operator, myPath, 0, vHat)     # Compute A0
        k    = np.linalg.matrix_rank(A0, tol=rankTol)   # Rank test

        if(k == 0):                                     # Rank is zero?
            raise RuntimeError(zeroRankErrorStr())      #  -> Stop
        elif(k == m):                                   # Maximum rank reached?
            raise RuntimeError(maxRankErrorStr())       #  -> Stop
        elif(k < l):                                    # Found a null SV?
            hasK = True                                 #  -> We have A0
        else:                                           # Matrix is full rank?
            l = l + lStep                               #  -> Increase A0 size

        it += 1                                         # Keep on searching A0

    # Check if maxIt was reached
    if(it == maxIt):
        raise RuntimeError(maxItErrorStr())
    else:
        if(verbose): print "Constructing linear EVP..."

    # Compute V, S and Wh
    #  NB: For SVD(A) = V*S*Wh, numpy computes {v, s, w}, such that:
    #      v = V; diag(s) = S and w = Wh
    V, S, Wh = np.linalg.svd(A0, full_matrices=False, compute_uv=1)

    # Extract V0, W0 and S0Inv
    V0    = np.delete(V,  l-1, 1)
    W0    = np.delete(Wh, l-1, 0).H
    S0Inv = np.matrix(np.diag(1/np.delete(S, l-1, 0)))

    # Compute A1 and B
    A1 = integrate(operator, myPath, 1, vHat)
    B  = V0.H * A1 * W0 * S0Inv

    # Eigenvalues of B
    if(verbose): print "Solving linear EVP..."
    myLambda, QHat = numpy.linalg.eig(B)

    # Done
    if(verbose): print "Done!"
    return myLambda


## Import only simple (other functions are just helpers)
__all__ = ['simple']


## Helper functions
def path(nodes, origin, radius):
    """Returns a list with the coordinates of a circular contour

    Keyword arguments:
    nodes -- the number of nodes used to discretise the contour
    radius -- the radius of the circular contour
    origin -- the origin (in the complex plane) of the circular contour

    The returned list contains nodes+1 elements,
    such that the first and the last are identical (up the machine precision)
    """

    step         = 1.0 / nodes
    nodesPlusOne = nodes + 1

    path = list()
    for i in range(nodesPlusOne):
        path.append(origin + radius * np.exp(1j * 2 * np.pi * i * step))

    return path


def integrate(operator, path, order, vHat):
    """Computes the countour integral of Beyn's method, that is matrix A_p

    Keyword arguments:
    operator -- the solver defining the operator to use
    path -- the path to integrate on
    order -- the order of Beyn's integral (that is the 'p' defining matrix A_p)
    vHat -- the RHS matrix defing Beyn's integral
    """

    # Initialise I
    I = np.matlib.zeros(vHat.shape, dtype=complex)

    # Initialise integration loop
    F1  = multiSolve(operator, vHat, path[0])
    F1 *= np.power(path[0], order)

    # Integration loop
    pathSizeMinus = len(path) - 1;
    for i in range(pathSizeMinus):
        F2  = multiSolve(operator, vHat, path[i + 1])
        F2 *= np.power(path[i + 1], order)

        I += (F1 + F2) * (path[i + 1] - path[i])
        F1 = F2

    # Done
    return I


def multiSolve(solver, B, w):
    """Solves for multiple RHS

    Keyword arguments:
    solver -- the solver to use
    B -- the matrix of RHS
    w -- the complex frequency to use
    """

    # Number of solve
    nSolve = B.shape[1]

    # Initialise X
    size = solver.size()
    X    = np.matlib.empty((size, nSolve), dtype=complex)

    # Loop and solve
    for i in range(nSolve):
        b = B[:, i]
        solver.solve(b, w)
        X[:, i] = solver.solution()

    # Done
    return X


def randomMatrix(n, m):
    """Returns a random complex matrix of the given size (n x m)"""
    return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j


def display(nodes, maxIt, lStart, lStep, rankTol, origin, radius):
    print "Beyn's contour integral method (simple)"
    print "---------------------------------------"
    print " # Nodes used for the trapezoidal rule:" + " " + str(nodes)
    print " # Maximum number of iterations:       " + " " + str(maxIt)
    print " # Initial size of col(A0):            " + " " + str(lStart)
    print " # Step size for col(A0):              " + " " + str(lStep)
    print " # Rank test tolerance:                ",
    print format(rankTol, '.2e')
    print "---------------------------------------"
    print " # Cirular path origin:                ",
    print "(" + format(np.real(origin).tolist(), '+.2e') + ")",
    print "+",
    print "(" + format(np.imag(origin).tolist(), '+.2e') + ")j"
    print " # Cirular path radius:                ",
    print format(radius, '+.2e')
    print "---------------------------------------"


def zeroRankErrorStr():
    """Returns a string explaining the probable reason of a zero rank"""
    return ("Found a rank of zero: " +
            "the contour is probably enclosing no eigenvalues")


def maxRankErrorStr():
    """Returns a string explaining the probable reason of a maximal rank"""
    return ("Maximal rank found: " +
            "the contour is probably enclosing to many eigenvalues")


def maxItErrorStr():
    """Returns a string claiming: the maximum number of iterations is reached"""
    return "Maximum number iterations is reached!"