""" Bifurcation point classes.  Each class locates and processes bifurcation points.

    * _BranchPointFold is a version based on BranchPoint location algorithms
    * BranchPoint: Branch process is broken (can't find alternate branch -- see MATCONT notes)

    Drew LaMar, March 2006
"""

from __future__ import absolute_import, print_function

from .misc import *
from PyDSTool.common import args
from .TestFunc import DiscreteMap, FixedPointMap

from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, \
    subtract, divide, transpose, eye, real, imag, \
    conjugate, average
from scipy import optimize, linalg
from numpy import dot as matrixmultiply
from numpy import array, float, complex, int, float64, complex64, int32, \
    zeros, divide, subtract, reshape, argsort, nonzero

#####
_classes = ['BifPoint', 'BPoint', 'BranchPoint', 'FoldPoint', 'HopfPoint',
            'BTPoint', 'ZHPoint', 'CPPoint',
            'BranchPointFold', '_BranchPointFold', 'DHPoint',
            'GHPoint', 'LPCPoint', 'PDPoint', 'NSPoint', 'SPoint']

__all__ = _classes
#####

class BifPoint(object):
    def __init__(self, testfuncs, flagfuncs, label='Bifurcation', stop=False):
        self.testfuncs = []
        self.flagfuncs = []
        self.found = []
        self.label = label
        self.stop = stop
        self.data = args()

        if not isinstance(testfuncs, list):
            testfuncs = [testfuncs]
        if not isinstance(flagfuncs, list):
            flagfuncs = [flagfuncs]

        self.testfuncs.extend(testfuncs)
        self.flagfuncs.extend(flagfuncs)
        self.tflen = len(self.testfuncs)

    def locate(self, P1, P2, C):
        pointlist = []

        for i, testfunc in enumerate(self.testfuncs):
            if self.flagfuncs[i] == iszero:
                for ind in range(testfunc.m):
                    X, V = testfunc.findzero(P1, P2, ind)
                    pointlist.append((X,V))

        X = average([point[0] for point in pointlist], axis=0)
        V = average([point[1] for point in pointlist], axis=0)
        C.Corrector(X,V)

        return X, V

    def process(self, X, V, C):
        data = args()
        data.X = todict(C, X)
        data.V = todict(C, V)

        self.found.append(data)

    def info(self, C, ind=None, strlist=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        if C.verbosity >= 1:
            print(self.label + ' Point found ')
        if C.verbosity >= 2:
            print('========================== ')
            for n, i in enumerate(ind):
                print(n, ': ')
                Xd = self.found[i].X
                for k, j in Xd.items():
                    print(k, ' = ', j)
                print('')
                if hasattr(self.found[i], 'eigs'):
                    print('Eigenvalues = \n')
                    for x in self.found[i].eigs:
                        print('     (%f,%f)' % (x.real, x.imag))
                    print('\n')
                if strlist is not None:
                    for string in strlist:
                        print(string)
                    print('')


class SPoint(BifPoint):
    """Special point that represents user-selected free parameter values."""
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'S', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)
        self.info(C, -1)
        return True


class BPoint(BifPoint):
    """Special point that represents boundary of computational domain."""
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'B', stop=stop)

    def locate(self, P1, P2, C):
        # Find location that triggered testfunc and initialize testfunc to that index
        val1 = (P1[0]-self.testfuncs[0].lower)*(self.testfuncs[0].upper-P1[0])
        val2 = (P2[0]-self.testfuncs[0].lower)*(self.testfuncs[0].upper-P2[0])
        ind = nonzero(val1*val2 < 0)
        self.testfuncs[0].ind = ind
        self.testfuncs[0].func = self.testfuncs[0].one

        X, V = BifPoint.locate(self, P1, P2, C)

        # Set testfunc back to monitoring all
        self.testfuncs[0].ind = None
        self.testfuncs[0].func = self.testfuncs[0].all

        return X, V

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        BifPoint.info(self, C, ind)



class BranchPoint(BifPoint):
    """May only work for EquilibriumCurve ... (needs fixing)"""
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'BP', stop=stop)

    def __locate_newton(self, X, C):
        """x[0:self.dim] = (x,alpha)
           x[self.dim] = beta
           x[self.dim+1:2*self.dim] = p
        """
        J_coords = C.CorrFunc.jac(X[0:C.dim], C.coords)
        J_params = C.CorrFunc.jac(X[0:C.dim], C.params)

        return r_[C.CorrFunc(X[0:C.dim]) + X[C.dim]*X[C.dim+1:], \
                  matrixmultiply(transpose(J_coords),X[C.dim+1:]), \
                  matrixmultiply(transpose(X[C.dim+1:]),J_params), \
                  matrixmultiply(transpose(X[C.dim+1:]),X[C.dim+1:]) - 1]

    def locate(self, P1, P2, C):
        # Initiliaze p vector to eigenvector with smallest eigenvalue
        X, V = P1
        X2, V2 = P2

        J_coords = C.CorrFunc.jac(X, C.coords)

        W, VL = linalg.eig(J_coords, left=1, right=0)
        ind = argsort([abs(eig) for eig in W])[0]
        p = real(VL[:,ind])

        initpoint = zeros(2*C.dim, float)
        initpoint[0:C.dim] = X
        initpoint[C.dim+1:] = p

        X = optimize.fsolve(self.__locate_newton, initpoint, C)
        self.data.psi = X[C.dim+1:]
        X = X[0:C.dim]

        V = 0.5*(V+V2)

        return X, V

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        # Finds the new branch
        J_coords = C.CorrFunc.jac(X, C.coords)
        J_params = C.CorrFunc.jac(X, C.params)


        singular = True
        perpvec = r_[1,zeros(C.dim-1)]
        d = 1
        while singular and d <= C.dim:
            try:
                v0 = linalg.solve(r_[c_[J_coords, J_params],
                                  [perpvec]], \
                                  r_[zeros(C.dim-1),1])
            except:
                perpvec = r_[0., perpvec[0:(C.dim-1)]]
                d += 1
            else:
                singular = False

        if singular:
            raise PyDSTool_ExistError("Problem in _compute: Failed to compute tangent vector.")
        v0 /= linalg.norm(v0)
        V = sign([x for x in v0 if abs(x) > 1e-8][0])*v0

        A = r_[c_[J_coords, J_params], [V]]
        W, VR = linalg.eig(A)
        W0 = [ind for ind, eig in enumerate(W) if abs(eig) < 5e-5]
        V1 = real(VR[:,W0[0]])

        H = C.CorrFunc.hess(X, C.coords+C.params, C.coords+C.params)
        c11 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V, V) for i in range(H.shape[0])])
        c12 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V, V1) for i in range(H.shape[0])])
        c22 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V1, V1) for i in range(H.shape[0])])

        beta = 1
        alpha = -1*c22/(2*c12)
        V1 = alpha*V + beta*V1
        V1 /= linalg.norm(V1)

        self.found[-1].eigs = W
        self.found[-1].branch = todict(C, V1)

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        for n, i in enumerate(ind):
            strlist.append('branch angle = ' + repr(matrixmultiply(tocoords(C, self.found[i].V), \
                tocoords(C, self.found[i].branch))))

        X = tocoords(C, self.found[-1].X)
        V = tocoords(C, self.found[-1].V)
        C._preTestFunc(X, V)
        strlist.append('Test function #1: ' + repr(self.testfuncs[0](X,V)[0]))

        BifPoint.info(self, C, ind, strlist)



class FoldPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'LP', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        # Compute normal form coefficient
        # NOTE: These are for free when using bordering technique!)
        # NOTE: Does not agree with MATCONT output! (if |p| = |q| = 1, then it does)
        J_coords = C.CorrFunc.jac(X, C.coords)
        W, VL, VR = linalg.eig(J_coords, left=1, right=1)
        minW = min(abs(W))
        ind = [(abs(eig) < minW+1e-8) and (abs(eig) > minW-1e-8) for eig in W].index(True)
        p, q = real(VL[:,ind]), real(VR[:,ind])
        p /= matrixmultiply(p,q)

        B = C.CorrFunc.hess(X, C.coords, C.coords)
        self.found[-1].a = abs(0.5*matrixmultiply(p,[bilinearform(B[i,:,:], q, q) for i in range(B.shape[0])]))
        self.found[-1].eigs = W

        numzero = len([eig for eig in W if abs(eig) < 1e-4])
        if numzero > 1:
            if C.verbosity >= 2:
                print('Fold-Fold!\n')
            del self.found[-1]
            return False
        elif numzero == 0:
            if C.verbosity >= 2:
                print('False positive!\n')
            del self.found[-1]
            return False

        if C.verbosity >= 2:
            print('\nChecking...')
            print('  |q| = %f' % linalg.norm(q))
            print('  <p,q> = %f' % matrixmultiply(p,q))
            print('  |Aq| = %f' % linalg.norm(matrixmultiply(J_coords,q)))
            print('  |transpose(A)p| = %f\n' % linalg.norm(matrixmultiply(transpose(J_coords),p)))

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        for n, i in enumerate(ind):
            strlist.append('a = ' + repr(self.found[i].a))

        BifPoint.info(self, C, ind, strlist)



class HopfPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'H', stop=stop)

    def process(self, X, V, C):
        """Tolerance for eigenvalues a possible problem when checking for neutral saddles."""
        BifPoint.process(self, X, V, C)

        J_coords = C.CorrFunc.jac(X, C.coords)
        eigs, LV, RV = linalg.eig(J_coords,left=1,right=1)

        # Check for neutral saddles
        found = False
        for i in range(len(eigs)):
            if abs(imag(eigs[i])) < 1e-5:
                for j in range(i+1,len(eigs)):
                    if C.verbosity >= 2:
                        if abs(eigs[i]) < 1e-5 and abs(eigs[j]) < 1e-5:
                            print('Fold-Fold point found in Hopf!\n')
                        elif abs(imag(eigs[j])) < 1e-5 and abs(real(eigs[i]) + real(eigs[j])) < 1e-5:
                            print('Neutral saddle found!\n')
            elif abs(real(eigs[i])) < 1e-5:
                for j in range(i+1, len(eigs)):
                    if abs(real(eigs[j])) < 1e-5 and abs(real(eigs[i]) - real(eigs[j])) < 1e-5:
                        found = True
                        w = abs(imag(eigs[i]))
                        if imag(eigs[i]) > 0:
                            p = conjugate(LV[:,j])/linalg.norm(LV[:,j])
                            q = RV[:,i]/linalg.norm(RV[:,i])
                        else:
                            p = conjugate(LV[:,i])/linalg.norm(LV[:,i])
                            q = RV[:,j]/linalg.norm(RV[:,j])

        if not found:
            del self.found[-1]
            return False

        direc = conjugate(1/matrixmultiply(conjugate(p),q))
        p = direc*p

        # Alternate way to compute 1st lyapunov coefficient (from Kuznetsov [4])

        #print (1./(w*w))*real(1j*matrixmultiply(conjugate(p),b1)*matrixmultiply(conjugate(p),b3) + \
        #   w*matrixmultiply(conjugate(p),trilinearform(D,q,q,conjugate(q))))

        self.found[-1].w = w
        self.found[-1].l1 = firstlyapunov(X, C.CorrFunc, w, J_coords=J_coords, p=p, q=q, check=(C.verbosity==2))
        self.found[-1].eigs = eigs

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        for n, i in enumerate(ind):
            strlist.append('w = ' + repr(self.found[i].w))
            strlist.append('l1 = ' + repr(self.found[i].l1))

        BifPoint.info(self, C, ind, strlist)


# Codimension-2 bifurcations

class BTPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'BT', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.CorrFunc.sysfunc.jac(X, C.coords)

        W, VL, VR = linalg.eig(J_coords, left=1, right=1)

        self.found[-1].eigs = W

        if C.verbosity >= 2:
            if C.CorrFunc.testfunc.data.B.shape[1] == 2:
                b = matrixmultiply(transpose(J_coords), C.CorrFunc.testfunc.data.w[:,0])
                c = matrixmultiply(J_coords, C.CorrFunc.testfunc.data.v[:,0])
            else:
                b = C.CorrFunc.testfunc.data.w[:,0]
                c = C.CorrFunc.testfunc.data.v[:,0]
            print('\nChecking...')
            print('  <b,c> = %f' % matrixmultiply(transpose(b), c))
            print('\n')

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        BifPoint.info(self, C, ind)



class ZHPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'ZH', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.CorrFunc.sysfunc.jac(X, C.coords)

        W, VL, VR = linalg.eig(J_coords, left=1, right=1)

        self.found[-1].eigs = W

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        BifPoint.info(self, C, ind)



class CPPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'CP', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.CorrFunc.sysfunc.jac(X, C.coords)
        B = C.CorrFunc.sysfunc.hess(X, C.coords, C.coords)

        W, VL, VR = linalg.eig(J_coords, left=1, right=1)

        q = C.CorrFunc.testfunc.data.C/linalg.norm(C.CorrFunc.testfunc.data.C)
        p = C.CorrFunc.testfunc.data.B/matrixmultiply(transpose(C.CorrFunc.testfunc.data.B),q)

        self.found[-1].eigs = W

        a = 0.5*matrixmultiply(transpose(p), reshape([bilinearform(B[i,:,:], q, q) \
                for i in range(B.shape[0])],(B.shape[0],1)))[0][0]
        if C.verbosity >= 2:
            print('\nChecking...')
            print('  |a| = %f' % a)
            print('\n')

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        BifPoint.info(self, C, ind)


class BranchPointFold(BifPoint):
    """Check Equilibrium.m in MATCONT"""
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'BP', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        pind = self.testfuncs[0].pind

        # Finds the new branch
        J_coords = C.CorrFunc.jac(X, C.coords)
        J_params = C.CorrFunc.jac(X, C.params)

        A = r_[c_[J_coords, J_params[:,pind]]]
        #A = r_[c_[J_coords, J_params], [V]]

        W, VR = linalg.eig(A)
        W0 = [ind for ind, eig in enumerate(W) if abs(eig) < 5e-5]
        tmp = real(VR[:,W0[0]])
        V1 = r_[tmp[:-1], 0, 0]
        V1[len(tmp)-1+pind] = tmp[-1]

        """NEED TO FIX THIS!"""
        H = C.CorrFunc.hess(X, C.coords+C.params, C.coords+C.params)
        # c11 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V, V) for i in range(H.shape[0])])
        # c12 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V, V1) for i in range(H.shape[0])])
        # c22 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V1, V1) for i in range(H.shape[0])])

        # beta = 1
        # alpha = -1*c22/(2*c12)
        # V1 = alpha*V + beta*V1
        # V1 /= linalg.norm(V1)

        self.found[-1].eigs = W
        self.found[-1].branch = None
        self.found[-1].par = C.freepars[self.testfuncs[0].pind]
        # self.found[-1].branch = todict(C, V1)

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        #for n, i in enumerate(ind):
        #    strlist.append('branch angle = ' + repr(matrixmultiply(tocoords(C, self.found[i].V), \
        #        tocoords(C, self.found[i].branch))))

        X = tocoords(C, self.found[-1].X)
        V = tocoords(C, self.found[-1].V)
        C._preTestFunc(X, V)
        strlist.append('Test function #1: ' + repr(self.testfuncs[0](X,V)[0]))

        BifPoint.info(self, C, ind, strlist)

class _BranchPointFold(BifPoint):
    """Check Equilibrium.m in MATCONT"""
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'BP', stop=stop)

    def __locate_newton(self, X, C):
        """Note:  This is redundant!! B is a column of A!!!  Works for now, though..."""
        pind = self.testfuncs[0].pind

        J_coords = C.CorrFunc.jac(X[0:C.dim], C.coords)
        J_params = C.CorrFunc.jac(X[0:C.dim], C.params)

        A = c_[J_coords, J_params[:,pind]]
        B = J_params[:,pind]

        return r_[C.CorrFunc(X[0:C.dim]) + X[C.dim]*X[C.dim+1:], \
                  matrixmultiply(transpose(A),X[C.dim+1:]), \
                  matrixmultiply(transpose(X[C.dim+1:]),B), \
                  matrixmultiply(transpose(X[C.dim+1:]),X[C.dim+1:]) - 1]

    def locate(self, P1, P2, C):
        # Initiliaze p vector to eigenvector with smallest eigenvalue
        X, V = P1

        pind = self.testfuncs[0].pind

        J_coords = C.CorrFunc.jac(X, C.coords)
        J_params = C.CorrFunc.jac(X, C.params)

        A = r_[c_[J_coords, J_params[:,pind]]]

        W, VL = linalg.eig(A, left=1, right=0)
        ind = argsort([abs(eig) for eig in W])[0]
        p = real(VL[:,ind])

        initpoint = zeros(2*C.dim, float)
        initpoint[0:C.dim] = X
        initpoint[C.dim+1:] = p

        X = optimize.fsolve(self.__locate_newton, initpoint, C)
        self.data.psi = X[C.dim+1:]
        X = X[0:C.dim]

        return X, V

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        pind = self.testfuncs[0].pind

        # Finds the new branch
        J_coords = C.CorrFunc.jac(X, C.coords)
        J_params = C.CorrFunc.jac(X, C.params)

        A = r_[c_[J_coords, J_params[:,pind]]]
        #A = r_[c_[J_coords, J_params], [V]]

        W, VR = linalg.eig(A)
        W0 = [ind for ind, eig in enumerate(W) if abs(eig) < 5e-5]
        tmp = real(VR[:,W0[0]])
        V1 = r_[tmp[:-1], 0, 0]
        V1[len(tmp)-1+pind] = tmp[-1]

        """NEED TO FIX THIS!"""
        H = C.CorrFunc.hess(X, C.coords+C.params, C.coords+C.params)
        c11 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V, V) for i in range(H.shape[0])])
        c12 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V, V1) for i in range(H.shape[0])])
        c22 = matrixmultiply(self.data.psi,[bilinearform(H[i,:,:], V1, V1) for i in range(H.shape[0])])

        beta = 1
        alpha = -1*c22/(2*c12)
        V1 = alpha*V + beta*V1
        V1 /= linalg.norm(V1)

        self.found[-1].eigs = W
        self.found[-1].branch = None
        self.found[-1].par = C.freepars[self.testfuncs[0].pind]
        self.found[-1].branch = todict(C, V1)

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        #for n, i in enumerate(ind):
        #    strlist.append('branch angle = ' + repr(matrixmultiply(tocoords(C, self.found[i].V), \
        #        tocoords(C, self.found[i].branch))))

        X = tocoords(C, self.found[-1].X)
        V = tocoords(C, self.found[-1].V)
        C._preTestFunc(X, V)
        strlist.append('Test function #1: ' + repr(self.testfuncs[0](X,V)[0]))

        BifPoint.info(self, C, ind, strlist)


class DHPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'DH', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.CorrFunc.sysfunc.jac(X, C.coords)
        eigs, LV, RV = linalg.eig(J_coords,left=1,right=1)

        self.found[-1].eigs = eigs

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        BifPoint.info(self, C, ind)



class GHPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'GH', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.CorrFunc.sysfunc.jac(X, C.coords)
        eigs, LV, RV = linalg.eig(J_coords,left=1,right=1)

        # Check for neutral saddles
        found = False
        for i in range(len(eigs)):
            if abs(imag(eigs[i])) < 1e-5:
                for j in range(i+1,len(eigs)):
                    if C.verbosity >= 2:
                        if abs(eigs[i]) < 1e-5 and abs(eigs[j]) < 1e-5:
                                print('Fold-Fold point found in Hopf!\n')
                        elif abs(imag(eigs[j])) < 1e-5 and abs(real(eigs[i]) + real(eigs[j])) < 1e-5:
                                print('Neutral saddle found!\n')
            elif abs(real(eigs[i])) < 1e-5:
                for j in range(i+1, len(eigs)):
                    if abs(real(eigs[j])) < 1e-5 and abs(real(eigs[i]) - real(eigs[j])) < 1e-5:
                        found = True
                        w = abs(imag(eigs[i]))
                        if imag(eigs[i]) > 0:
                            p = conjugate(LV[:,j]/linalg.norm(LV[:,j]))
                            q = RV[:,i]/linalg.norm(RV[:,i])
                        else:
                            p = conjugate(LV[:,i]/linalg.norm(LV[:,i]))
                            q = RV[:,j]/linalg.norm(RV[:,j])

        if not found:
            del self.found[-1]
            return False

        direc = conjugate(1/matrixmultiply(conjugate(p),q))
        p = direc*p

        # Alternate way to compute 1st lyapunov coefficient (from Kuznetsov [4])

        #print (1./(w*w))*real(1j*matrixmultiply(conjugate(p),b1)*matrixmultiply(conjugate(p),b3) + \
        #   w*matrixmultiply(conjugate(p),trilinearform(D,q,q,conjugate(q))))

        self.found[-1].w = w
        self.found[-1].l1 = firstlyapunov(X, C.CorrFunc.sysfunc, w, J_coords=J_coords, p=p, q=q, check=(C.verbosity==2))
        self.found[-1].eigs = eigs

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        for n, i in enumerate(ind):
            strlist.append('w = ' + repr(self.found[i].w))
            strlist.append('l1 = ' + repr(self.found[i].l1))

        BifPoint.info(self, C, ind, strlist)



# Discrete maps

class LPCPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'LPC', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.sysfunc.jac(X, C.coords)

        W, VL, VR = linalg.eig(J_coords, left=1, right=1)

        self.found[-1].eigs = W

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        X = tocoords(C, self.found[-1].X)
        V = tocoords(C, self.found[-1].V)
        C._preTestFunc(X, V)
        strlist.append('Test function #1: ' + repr(self.testfuncs[0](X,V)[0]))
        strlist.append('Test function #2: ' + repr(self.testfuncs[1](X,V)[0]))

        BifPoint.info(self, C, ind, strlist)

class PDPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'PD', stop=stop)

    def process(self, X, V, C):
        """Do I need to compute the branch, or will it always be in the direction of freepar = constant?"""
        BifPoint.process(self, X, V, C)

        F = DiscreteMap(C.sysfunc, period=2*C.sysfunc.period)
        FP = FixedPointMap(F)

        J_coords = FP.jac(X, C.coords)
        J_params = FP.jac(X, C.params)

        # Locate branch of double period map
        W, VL = linalg.eig(J_coords, left=1, right=0)
        ind = argsort([abs(eig) for eig in W])[0]
        psi = real(VL[:,ind])

        A = r_[c_[J_coords, J_params], [V]]
        W, VR = linalg.eig(A)
        W0 = argsort([abs(eig) for eig in W])[0]
        V1 = real(VR[:,W0])

        H = FP.hess(X, C.coords+C.params, C.coords+C.params)
        c11 = matrixmultiply(psi,[bilinearform(H[i,:,:], V, V) for i in range(H.shape[0])])
        c12 = matrixmultiply(psi,[bilinearform(H[i,:,:], V, V1) for i in range(H.shape[0])])
        c22 = matrixmultiply(psi,[bilinearform(H[i,:,:], V1, V1) for i in range(H.shape[0])])

        beta = 1
        alpha = -1*c22/(2*c12)
        V1 = alpha*V + beta*V1
        V1 /= linalg.norm(V1)

        J_coords = C.sysfunc.jac(X, C.coords)
        W = linalg.eig(J_coords, right=0)

        self.found[-1].eigs = W
        self.found[-1].branch_period = 2*C.sysfunc.period
        self.found[-1].branch = todict(C, V1)

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        strlist = []
        for n, i in enumerate(ind):
            strlist.append('Period doubling branch angle = ' + repr(matrixmultiply(tocoords(C, self.found[i].V), \
                tocoords(C, self.found[i].branch))))

        BifPoint.info(self, C, ind, strlist)



class NSPoint(BifPoint):
    def __init__(self, testfuncs, flagfuncs, stop=False):
        BifPoint.__init__(self, testfuncs, flagfuncs, 'NS', stop=stop)

    def process(self, X, V, C):
        BifPoint.process(self, X, V, C)

        J_coords = C.sysfunc.jac(X, C.coords)

        eigs, VL, VR = linalg.eig(J_coords, left=1, right=1)

        # Check for nonreal multipliers
        found = False
        for i in range(len(eigs)):
            for j in range(i+1,len(eigs)):
                if abs(imag(eigs[i])) > 1e-10 and \
                   abs(imag(eigs[j])) > 1e-10 and \
                   abs(eigs[i]*eigs[j] - 1) < 1e-5:
                    found = True

        if not found:
            del self.found[-1]
            return False

        self.found[-1].eigs = eigs

        self.info(C, -1)

        return True

    def info(self, C, ind=None):
        if ind is None:
            ind = list(range(len(self.found)))
        elif isinstance(ind, int):
            ind = [ind]

        BifPoint.info(self, C, ind)

[((1327, 1333), 'PyDSTool.common.args', 'args', ([], {}), '()\n', (1331, 1333), False, 'from PyDSTool.common import args\n'), ((1936, 1986), 'numpy.average', 'average', (['[point[0] for point in pointlist]'], {'axis': '(0)'}), '([point[0] for point in pointlist], axis=0)\n', (1943, 1986), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((1999, 2049), 'numpy.average', 'average', (['[point[1] for point in pointlist]'], {'axis': '(0)'}), '([point[1] for point in pointlist], axis=0)\n', (2006, 2049), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((2144, 2150), 'PyDSTool.common.args', 'args', ([], {}), '()\n', (2148, 2150), False, 'from PyDSTool.common import args\n'), ((4054, 4078), 'numpy.nonzero', 'nonzero', (['(val1 * val2 < 0)'], {}), '(val1 * val2 < 0)\n', (4061, 4078), False, 'from numpy import array, float, complex, int, float64, complex64, int32, zeros, divide, subtract, reshape, argsort, nonzero\n'), ((5679, 5716), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(0)'}), '(J_coords, left=1, right=0)\n', (5689, 5716), False, 'from scipy import optimize, linalg\n'), ((5779, 5797), 'numpy.real', 'real', (['VL[:, (ind)]'], {}), '(VL[:, (ind)])\n', (5783, 5797), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((5816, 5839), 'numpy.zeros', 'zeros', (['(2 * C.dim)', 'float'], {}), '(2 * C.dim, float)\n', (5821, 5839), False, 'from numpy import array, float, complex, int, float64, complex64, int32, zeros, divide, subtract, reshape, argsort, nonzero\n'), ((5914, 5965), 'scipy.optimize.fsolve', 'optimize.fsolve', (['self.__locate_newton', 'initpoint', 'C'], {}), '(self.__locate_newton, initpoint, C)\n', (5929, 5965), False, 'from scipy import optimize, linalg\n'), ((6850, 6865), 'scipy.linalg.norm', 'linalg.norm', (['v0'], {}), '(v0)\n', (6861, 6865), False, 'from scipy import optimize, linalg\n'), ((6984, 6997), 'scipy.linalg.eig', 'linalg.eig', (['A'], {}), '(A)\n', (6994, 6997), False, 'from scipy import optimize, linalg\n'), ((7078, 7098), 'numpy.real', 'real', (['VR[:, (W0[0])]'], {}), '(VR[:, (W0[0])])\n', (7082, 7098), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((7569, 7584), 'scipy.linalg.norm', 'linalg.norm', (['V1'], {}), '(V1)\n', (7580, 7584), False, 'from scipy import optimize, linalg\n'), ((8819, 8856), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (8829, 8856), False, 'from scipy import optimize, linalg\n'), ((9036, 9056), 'numpy.dot', 'matrixmultiply', (['p', 'q'], {}), '(p, q)\n', (9050, 9056), True, 'from numpy import dot as matrixmultiply\n'), ((10738, 10775), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (10748, 10775), False, 'from scipy import optimize, linalg\n'), ((13289, 13326), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (13299, 13326), False, 'from scipy import optimize, linalg\n'), ((14418, 14455), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (14428, 14455), False, 'from scipy import optimize, linalg\n'), ((15103, 15140), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (15113, 15140), False, 'from scipy import optimize, linalg\n'), ((16427, 16440), 'scipy.linalg.eig', 'linalg.eig', (['A'], {}), '(A)\n', (16437, 16440), False, 'from scipy import optimize, linalg\n'), ((16522, 16542), 'numpy.real', 'real', (['VR[:, (W0[0])]'], {}), '(VR[:, (W0[0])])\n', (16526, 16542), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((19133, 19163), 'scipy.linalg.eig', 'linalg.eig', (['A'], {'left': '(1)', 'right': '(0)'}), '(A, left=1, right=0)\n', (19143, 19163), False, 'from scipy import optimize, linalg\n'), ((19226, 19244), 'numpy.real', 'real', (['VL[:, (ind)]'], {}), '(VL[:, (ind)])\n', (19230, 19244), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((19263, 19286), 'numpy.zeros', 'zeros', (['(2 * C.dim)', 'float'], {}), '(2 * C.dim, float)\n', (19268, 19286), False, 'from numpy import array, float, complex, int, float64, complex64, int32, zeros, divide, subtract, reshape, argsort, nonzero\n'), ((19361, 19412), 'scipy.optimize.fsolve', 'optimize.fsolve', (['self.__locate_newton', 'initpoint', 'C'], {}), '(self.__locate_newton, initpoint, C)\n', (19376, 19412), False, 'from scipy import optimize, linalg\n'), ((19841, 19854), 'scipy.linalg.eig', 'linalg.eig', (['A'], {}), '(A)\n', (19851, 19854), False, 'from scipy import optimize, linalg\n'), ((19936, 19956), 'numpy.real', 'real', (['VR[:, (W0[0])]'], {}), '(VR[:, (W0[0])])\n', (19940, 19956), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((20529, 20544), 'scipy.linalg.norm', 'linalg.norm', (['V1'], {}), '(V1)\n', (20540, 20544), False, 'from scipy import optimize, linalg\n'), ((21696, 21733), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (21706, 21733), False, 'from scipy import optimize, linalg\n'), ((22325, 22362), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (22335, 22362), False, 'from scipy import optimize, linalg\n'), ((24873, 24910), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (24883, 24910), False, 'from scipy import optimize, linalg\n'), ((26064, 26101), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(0)'}), '(J_coords, left=1, right=0)\n', (26074, 26101), False, 'from scipy import optimize, linalg\n'), ((26166, 26184), 'numpy.real', 'real', (['VL[:, (ind)]'], {}), '(VL[:, (ind)])\n', (26170, 26184), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((26243, 26256), 'scipy.linalg.eig', 'linalg.eig', (['A'], {}), '(A)\n', (26253, 26256), False, 'from scipy import optimize, linalg\n'), ((26319, 26336), 'numpy.real', 'real', (['VR[:, (W0)]'], {}), '(VR[:, (W0)])\n', (26323, 26336), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((26769, 26784), 'scipy.linalg.norm', 'linalg.norm', (['V1'], {}), '(V1)\n', (26780, 26784), False, 'from scipy import optimize, linalg\n'), ((26844, 26873), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'right': '(0)'}), '(J_coords, right=0)\n', (26854, 26873), False, 'from scipy import optimize, linalg\n'), ((27795, 27832), 'scipy.linalg.eig', 'linalg.eig', (['J_coords'], {'left': '(1)', 'right': '(1)'}), '(J_coords, left=1, right=1)\n', (27805, 27832), False, 'from scipy import optimize, linalg\n'), ((8990, 9008), 'numpy.real', 'real', (['VL[:, (ind)]'], {}), '(VL[:, (ind)])\n', (8994, 9008), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((9007, 9025), 'numpy.real', 'real', (['VR[:, (ind)]'], {}), '(VR[:, (ind)])\n', (9011, 9025), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((15181, 15220), 'scipy.linalg.norm', 'linalg.norm', (['C.CorrFunc.testfunc.data.C'], {}), '(C.CorrFunc.testfunc.data.C)\n', (15192, 15220), False, 'from scipy import optimize, linalg\n'), ((13556, 13616), 'numpy.dot', 'matrixmultiply', (['J_coords', 'C.CorrFunc.testfunc.data.v[:, (0)]'], {}), '(J_coords, C.CorrFunc.testfunc.data.v[:, (0)])\n', (13570, 13616), True, 'from numpy import dot as matrixmultiply\n'), ((15275, 15312), 'numpy.transpose', 'transpose', (['C.CorrFunc.testfunc.data.B'], {}), '(C.CorrFunc.testfunc.data.B)\n', (15284, 15312), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((6318, 6334), 'numpy.zeros', 'zeros', (['(C.dim - 1)'], {}), '(C.dim - 1)\n', (6323, 6334), False, 'from numpy import array, float, complex, int, float64, complex64, int32, zeros, divide, subtract, reshape, argsort, nonzero\n'), ((9722, 9736), 'scipy.linalg.norm', 'linalg.norm', (['q'], {}), '(q)\n', (9733, 9736), False, 'from scipy import optimize, linalg\n'), ((9773, 9793), 'numpy.dot', 'matrixmultiply', (['p', 'q'], {}), '(p, q)\n', (9787, 9793), True, 'from numpy import dot as matrixmultiply\n'), ((10887, 10900), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (10891, 10900), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((12040, 12052), 'numpy.conjugate', 'conjugate', (['p'], {}), '(p)\n', (12049, 12052), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((13482, 13501), 'numpy.transpose', 'transpose', (['J_coords'], {}), '(J_coords)\n', (13491, 13501), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((22474, 22487), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (22478, 22487), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23635, 23647), 'numpy.conjugate', 'conjugate', (['p'], {}), '(p)\n', (23644, 23647), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((5293, 5312), 'numpy.transpose', 'transpose', (['J_coords'], {}), '(J_coords)\n', (5302, 5312), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((5362, 5386), 'numpy.transpose', 'transpose', (['X[C.dim + 1:]'], {}), '(X[C.dim + 1:])\n', (5371, 5386), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((9840, 9867), 'numpy.dot', 'matrixmultiply', (['J_coords', 'q'], {}), '(J_coords, q)\n', (9854, 9867), True, 'from numpy import dot as matrixmultiply\n'), ((11324, 11337), 'numpy.real', 'real', (['eigs[i]'], {}), '(eigs[i])\n', (11328, 11337), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((13821, 13833), 'numpy.transpose', 'transpose', (['b'], {}), '(b)\n', (13830, 13833), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((15381, 15393), 'numpy.transpose', 'transpose', (['p'], {}), '(p)\n', (15390, 15393), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((18647, 18659), 'numpy.transpose', 'transpose', (['A'], {}), '(A)\n', (18656, 18659), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((18709, 18733), 'numpy.transpose', 'transpose', (['X[C.dim + 1:]'], {}), '(X[C.dim + 1:])\n', (18718, 18733), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((22919, 22932), 'numpy.real', 'real', (['eigs[i]'], {}), '(eigs[i])\n', (22923, 22932), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((5431, 5455), 'numpy.transpose', 'transpose', (['X[C.dim + 1:]'], {}), '(X[C.dim + 1:])\n', (5440, 5455), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((9943, 9962), 'numpy.transpose', 'transpose', (['J_coords'], {}), '(J_coords)\n', (9952, 9962), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((18771, 18795), 'numpy.transpose', 'transpose', (['X[C.dim + 1:]'], {}), '(X[C.dim + 1:])\n', (18780, 18795), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((27997, 28010), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (28001, 28010), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((28049, 28062), 'numpy.imag', 'imag', (['eigs[j]'], {}), '(eigs[j])\n', (28053, 28062), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((6550, 6566), 'numpy.zeros', 'zeros', (['(C.dim - 1)'], {}), '(C.dim - 1)\n', (6555, 6566), False, 'from numpy import array, float, complex, int, float64, complex64, int32, zeros, divide, subtract, reshape, argsort, nonzero\n'), ((11560, 11573), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (11564, 11573), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11602, 11615), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (11606, 11615), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23155, 23168), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (23159, 23168), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23197, 23210), 'numpy.imag', 'imag', (['eigs[i]'], {}), '(eigs[i])\n', (23201, 23210), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11422, 11435), 'numpy.real', 'real', (['eigs[j]'], {}), '(eigs[j])\n', (11426, 11435), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11653, 11674), 'numpy.conjugate', 'conjugate', (['LV[:, (j)]'], {}), '(LV[:, (j)])\n', (11662, 11674), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11672, 11695), 'scipy.linalg.norm', 'linalg.norm', (['LV[:, (j)]'], {}), '(LV[:, (j)])\n', (11683, 11695), False, 'from scipy import optimize, linalg\n'), ((11733, 11756), 'scipy.linalg.norm', 'linalg.norm', (['RV[:, (i)]'], {}), '(RV[:, (i)])\n', (11744, 11756), False, 'from scipy import optimize, linalg\n'), ((11816, 11837), 'numpy.conjugate', 'conjugate', (['LV[:, (i)]'], {}), '(LV[:, (i)])\n', (11825, 11837), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11835, 11858), 'scipy.linalg.norm', 'linalg.norm', (['LV[:, (i)]'], {}), '(LV[:, (i)])\n', (11846, 11858), False, 'from scipy import optimize, linalg\n'), ((11896, 11919), 'scipy.linalg.norm', 'linalg.norm', (['RV[:, (j)]'], {}), '(RV[:, (j)])\n', (11907, 11919), False, 'from scipy import optimize, linalg\n'), ((23017, 23030), 'numpy.real', 'real', (['eigs[j]'], {}), '(eigs[j])\n', (23021, 23030), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23328, 23351), 'scipy.linalg.norm', 'linalg.norm', (['RV[:, (i)]'], {}), '(RV[:, (i)])\n', (23339, 23351), False, 'from scipy import optimize, linalg\n'), ((23491, 23514), 'scipy.linalg.norm', 'linalg.norm', (['RV[:, (j)]'], {}), '(RV[:, (j)])\n', (23502, 23514), False, 'from scipy import optimize, linalg\n'), ((11173, 11186), 'numpy.imag', 'imag', (['eigs[j]'], {}), '(eigs[j])\n', (11177, 11186), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11452, 11465), 'numpy.real', 'real', (['eigs[i]'], {}), '(eigs[i])\n', (11456, 11465), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11468, 11481), 'numpy.real', 'real', (['eigs[j]'], {}), '(eigs[j])\n', (11472, 11481), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((22764, 22777), 'numpy.imag', 'imag', (['eigs[j]'], {}), '(eigs[j])\n', (22768, 22777), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23047, 23060), 'numpy.real', 'real', (['eigs[i]'], {}), '(eigs[i])\n', (23051, 23060), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23063, 23076), 'numpy.real', 'real', (['eigs[j]'], {}), '(eigs[j])\n', (23067, 23076), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((23266, 23289), 'scipy.linalg.norm', 'linalg.norm', (['LV[:, (j)]'], {}), '(LV[:, (j)])\n', (23277, 23289), False, 'from scipy import optimize, linalg\n'), ((23429, 23452), 'scipy.linalg.norm', 'linalg.norm', (['LV[:, (i)]'], {}), '(LV[:, (i)])\n', (23440, 23452), False, 'from scipy import optimize, linalg\n'), ((11203, 11216), 'numpy.real', 'real', (['eigs[i]'], {}), '(eigs[i])\n', (11207, 11216), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((11219, 11232), 'numpy.real', 'real', (['eigs[j]'], {}), '(eigs[j])\n', (11223, 11232), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((22794, 22807), 'numpy.real', 'real', (['eigs[i]'], {}), '(eigs[i])\n', (22798, 22807), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n'), ((22810, 22823), 'numpy.real', 'real', (['eigs[j]'], {}), '(eigs[j])\n', (22814, 22823), False, 'from numpy import Inf, NaN, isfinite, r_, c_, sign, mod, subtract, divide, transpose, eye, real, imag, conjugate, average\n')]