def modularity_spectrum(G):
    """Return eigenvalues of the modularity matrix of G.

    Parameters
    ----------
    G : Graph
       A NetworkX Graph or DiGraph

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    See Also
    --------
    modularity_matrix

    References
    ----------
    .. [1] M. E. J. Newman, "Modularity and community structure in networks",
       Proc. Natl. Acad. Sci. USA, vol. 103, pp. 8577-8582, 2006.
    """
    from scipy.linalg import eigvals
    if G.is_directed():
        return eigvals(nx.directed_modularity_matrix(G))
    else:
        return eigvals(nx.modularity_matrix(G))

def get_graph_props(g):
    travgl = g.transitivity_avglocal_undirected()
    tru = g.transitivity_undirected()
    d = g.diameter()
    asd = g.assortativity_degree()
    apl = g.average_path_length()
    omega = g.omega()
    density = g.density()
    maxd = g.maxdegree()
    medd = np.median(g.degree())
    plaw = ig.power_law_fit(g.degree())
    spnorm = max(np.abs(spl.eigvals(g.get_adjacency().data)))
    leigs = spl.eigvals(g.laplacian())
    algc = abs(sorted([e for e in leigs if e > 1e-10])[1])
    
    return [travgl,     # avg. local transitivity
            tru,        # global transitivity
            d,          # diameter
            asd,        # degree assortativity
            apl,        # avg. path length
            omega,      # max. clique
            density,
            maxd,       # max. degree
            medd,       # median degree
            plaw.alpha, # power law exponent
            spnorm,     # largest eigenvalue of adj. matrix
            algc,       # 2nd smallest non-zero eigenvalue of laplacian
            ]

    def testMinrealBrute(self):
        for n, m, p in permutations(range(1,6), 3):
            s = matlab.rss(n, p, m)
            sr = s.minreal()
            if s.states > sr.states:
                self.nreductions += 1
            else:
                np.testing.assert_array_almost_equal(
                    np.sort(eigvals(s.A)), np.sort(eigvals(sr.A)))
                for i in range(m):
                    for j in range(p):
                        ht1 = matlab.tf(
                            matlab.ss(s.A, s.B[:,i], s.C[j,:], s.D[j,i]))
                        ht2 = matlab.tf(
                            matlab.ss(sr.A, sr.B[:,i], sr.C[j,:], sr.D[j,i]))
                        try:
                            self.assert_numden_almost_equal(
                                ht1.num[0][0], ht2.num[0][0],
                                ht1.den[0][0], ht2.den[0][0])
                        except Exception as e:
                            # for larger systems, the tf minreal's
                            # the original rss, but not the balanced one
                            if n < 6:
                                raise e

        self.assertEqual(self.nreductions, 2)

def format_eig_svd():
    def format_cplx(z):
        if z.imag < 1e-300:
            return '{0:.4f}'.format(z.real)
        return '{0:.4f}+{1:.4f}i'.format(z.real, z.imag)

    eig12 = sp.eigvals(generate_matrix(12))
    svd12 = sp.svdvals(generate_matrix(12))

    eig25 = sp.eigvals(generate_matrix(25))
    svd25 = sp.svdvals(generate_matrix(25))

    result12 = r'\begin{tabular}{cc}' + '\n'
    result12 += r'    Eigenvalues&Singular values\\' + '\n'
    result12 += '     \\hline\n'
    result25 = copy.copy(result12)
    for k in range(25):
        if k < 12:
            result12 += r'    ${0}$&${1:.4f}$\\'.format(format_cplx(eig12[k]), svd12[k]) + '\n'
        result25 += r'    ${0}$&${1:.4f}$\\'.format(format_cplx(eig25[k]), svd25[k]) + '\n'

    result12 += '\\end{tabular}\n'
    result25 += '\\end{tabular}\n'

    print(result12)

    print(result25)

    def test_dare(self):
        A = matrix([[-0.6, 0],[-0.1, -0.4]])
        Q = matrix([[2, 1],[1, 0]])
        B = matrix([[2, 1],[0, 1]])
        R = matrix([[1, 0],[0, 1]])

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T * X * A - X -
            A.T * X * B * solve(B.T * X * B + R, B.T * X * A) + Q, zeros((2,2)))
        assert_array_almost_equal(solve(B.T * X * B + R, B.T * X * A), G)
        # check for stable closed loop
        lam = eigvals(A - B * G)
        assert_array_less(abs(lam), 1.0)

        A = matrix([[1, 0],[-1, 1]])
        Q = matrix([[0, 1],[1, 1]])
        B = matrix([[1],[0]])
        R = 2

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T * X * A - X -
            A.T * X * B * solve(B.T *  X * B + R, B.T * X * A) + Q, zeros((2,2)))
        assert_array_almost_equal(B.T * X * A / (B.T * X * B + R), G)
        # check for stable closed loop
        lam = eigvals(A - B * G)
        assert_array_less(abs(lam), 1.0)

def set_aw(sys,poles):
    """Divide in controller in input and feedback part
       for anti-windup

    Usage
    =====
    [sys_in,sys_fbk]=set_aw(sys,poles)

    Inputs
    ------

    sys: controller
    poles : poles for the anti-windup filter

    Outputs
    -------
    sys_in, sys_fbk: controller in input and feedback part
    """
    sys=ss(sys);
    den_old=poly(eigvals(sys.A))
    den = poly(poles)
    tmp= tf(den_old,den,sys.Tsamp)
    tmpss=tf2ss(tmp)
    sys_in=ss(tmp*sys)
    sys_in.Tsamp=sys.Tsamp
    sys_fbk=ss(1-tmp)
    sys_fbk.Tsamp=sys.Tsamp
    return sys_in, sys_fbk

def bb_step(sys,X0=None,Tf=None,Ts=0.001):
    """Plot the step response of the continous system sys

    Call:
    y=bb_step(sys [,Tf=final time] [,Ts=time step])

    Parameters
    ----------
    sys : Continous System in State Space form
    X0: Initial state vector (not used yet)
    Ts  : sympling time
    Tf  : Final simulation time
 
    Returns
    -------
    Nothing

    """
    if Tf==None:
        vals = eigvals(sys.A)
        r = min(abs(real(vals)))
        if r < 1e-10:
            r = 0.1
        Tf = 7.0 / r
    sysd=c2d(sys,Ts)
    dstep(sysd,Tf=Tf)

    def testMinreal(self, verbose=False):
        """Test a minreal model reduction"""
        #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
        A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
        #B = [0.3, -1.3; 0.1, 0; 1, 0]
        B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
        #C = [0, 0.1, 0; -0.3, -0.2, 0]
        C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
        #D = [0 -0.8; -0.3 0]
        D = [[0., -0.8], [-0.3, 0.]]
        # sys = ss(A, B, C, D)

        sys = ss(A, B, C, D)
        sysr = minreal(sys)
        self.assertEqual(sysr.states, 2)
        self.assertEqual(sysr.inputs, sys.inputs)
        self.assertEqual(sysr.outputs, sys.outputs)
        np.testing.assert_array_almost_equal(
            eigvals(sysr.A), [-2.136154, -0.1638459])

        s = tf([1, 0], [1])
        h = (s+1)*(s+2.00000000001)/(s+2)/(s**2+s+1)
        hm = minreal(h)
        hr = (s+1)/(s**2+s+1)
        np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0])
        np.testing.assert_array_almost_equal(hm.den[0][0], hr.den[0][0])

def adjacency_spectrum(G, weight='weight'):
    """Return eigenvalues of the adjacency matrix of G.

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    adjacency_matrix
    """
    from scipy.linalg import eigvals
    return eigvals(nx.adjacency_matrix(G,weight=weight).todense())

def ar_model_check_stable(A):
    """check if this AR model is stable

    :Parameters:
        A : ndarray
            The coefficient matrix of the model
    """

    # inits and checks
    m, p = A.shape
    p /= m
    if p != round(p):
        raise ValueError('bad inputs!')

    # check for stable model
    A1 = N.concatenate((
        A,
        N.concatenate((
            N.eye((p - 1) * m),
            N.zeros(((p - 1) * m, m))
        ), axis=1)
    ))
    lambdas = NL.eigvals(A1)
    rval = True
    if (N.absolute(lambdas) > 1).any():
        rval = False
    del A1, lambdas
    return rval

def compare_eigen_methods():
    """ timing of different eigensolver methods"""
    import scipy.linalg as linalg 
    print '\n *** diagonalize a real symmetric band matrix by different methods ***\n'
    mt=mytimer.create(10)    
    try:
        N=float(sys.argv[1])
        S=float(sys.argv[2])
    except:
        sys.exit('supply N S (command line arguments): matrix dimension N and and number of super-diagonals S ')

    np.random.seed(1)
    a=np.random.random((N,N))

    ab=GMatrix(np.random.random((S+1,N)),storage='upper_banded')
    mt[0].start('lapack symmetric upper banded')
    print ab.store+'\n',ab.eigvals()[:5]
    mt[0].stop()

    a=ab.restore('full')
    mt[1].start('lapack symmetric full')
    print a.store+'\n',a.eigvals()[:5]
    mt[1].stop()

    print 'lingalg'
    mt[2].start('linalg general')
    print np.sort(linalg.eigvals(a.m).real)[:5]
    mt[2].stop()

    mytimer.table()

    def test_timescales(self):
        P_dense=self.bdc.transition_matrix()
        P=self.bdc.transition_matrix_sparse()
        ev=eigvals(P_dense)
        """Sort with decreasing magnitude"""
        ev=ev[np.argsort(np.abs(ev))[::-1]]
        ts=-1.0/np.log(np.abs(ev))

        """k=None"""
        with self.assertRaises(ValueError):
            tsn=timescales(P)

        """k is not None"""
        tsn=timescales(P, k=self.k)
        self.assertTrue(np.allclose(ts[1:self.k], tsn[1:]))

        """k is not None, ncv is not None"""
        tsn=timescales(P, k=self.k, ncv=self.ncv)
        self.assertTrue(np.allclose(ts[1:self.k], tsn[1:]))
        

        """tau=7"""      

        """k is not None"""
        tsn=timescales(P, k=self.k, tau=7)
        self.assertTrue(np.allclose(7*ts[1:self.k], tsn[1:]))

def projection_shifts_init(A, E, B, shift_options):
    """Find starting shift parameters for low-rank ADI iteration using
    Galerkin projection on spaces spanned by LR-ADI iterates.

    See [PK16]_, pp. 92-95.

    Parameters
    ----------
    A
        The |Operator| A from the corresponding Lyapunov equation.
    E
        The |Operator| E from the corresponding Lyapunov equation.
    B
        The |VectorArray| B from the corresponding Lyapunov equation.
    shift_options
        The shift options to use (see :func:`lyap_solver_options`).

    Returns
    -------
    shifts
        A |NumPy array| containing a set of stable shift parameters.
    """
    for i in range(shift_options['init_maxiter']):
        Q = gram_schmidt(B, atol=0, rtol=0)
        shifts = spla.eigvals(A.apply2(Q, Q), E.apply2(Q, Q))
        shifts = shifts[np.real(shifts) < 0]
        if shifts.size == 0:
            # use random subspace instead of span{B} (with same dimensions)
            if shift_options['init_seed'] is not None:
                np.random.seed(shift_options['init_seed'])
                np.random.seed(np.random.random() + i)
            B = B.space.make_array(np.random.rand(len(B), B.space.dim))
        else:
            return shifts
    raise RuntimeError('Could not generate initial shifts for low-rank ADI iteration.')

	def __init__(self,para):	
		self.parms = Parameters()
		self.parms.N = para.N
		self.parms.Nin = para.Nin
		self.parms.Nout= para.Nout
		self.parms.inp_sc= para.inp_sc
		self.parms.spectr_rad= para.spectr_rad
		self.parms.leak_rate= para.leak_rate
		self.parms.inpBias = para.inpBias

		self.parms.Pscaling = para.Pscaling
		self.parms.OutScaling = para.OutScaling

		self.parms.alpha = para.alpha

		self.parms.t1 = para.t1
		self.parms.dt = para.dt

		self.parms.y_d = para.y_d

		self.P = np.eye(self.parms.N)*self.parms.Pscaling#100
		self.inpBias = self.parms.inpBias * np.random.randn(self.parms.N,1)

		wtemp = np.random.randn(self.parms.N,self.parms.N)
		self.W = self.parms.spectr_rad*wtemp/max(abs(linalg.eigvals(wtemp)))
		self.V = np.random.randn(self.parms.N,self.parms.Nin)*np.dot(np.ones((self.parms.N,1)),self.parms.inp_sc) 

		self.Woutp = np.random.randn(self.parms.Nout,self.parms.N)

		self.state_A = np.random.rand(self.parms.N,1)*2.-1.
		self.state_B = self.state_A
		self.orgstate = self.state_A

def _run_hamiltonian(verbose=True):
    c = classicalHamiltonian()
    if verbose:
        print(c.potential(array([-0.5, 0.5])))
        print(c.potential(array([-0.5, 0.0])))
        print(c.potential(array([0.0, 0.0])))

    xopt = optimize.fmin(c.potential, c.initialposition(), xtol=1e-10)

    hessian = nd.Hessian(c.potential)

    H = hessian(xopt)
    true_H = np.array([[5.23748385e-12, -2.61873829e-12],
                       [-2.61873829e-12, 5.23748385e-12]])
    error_estimate = np.NAN
    if verbose:
        print(xopt)
        print('H', H)
        print('H-true_H', np.abs(H - true_H))
        # print('error_estimate', info.error_estimate)

        eigenvalues = linalg.eigvals(H)
        normal_modes = c.normal_modes(eigenvalues)

        print('eigenvalues', eigenvalues)
        print('normal_modes', normal_modes)
    return H, error_estimate, true_H

def z2_vanderbilt(h,nk=30,nt=100,nocc=None,full=False):
  """ Calculate Z2 invariant according to Vanderbilt algorithm"""
  out = [] # output list
  path = np.linspace(0.,1.,nk) # set of kpoints
  fo = open("WANNIER_CENTERS.OUT","w")
  if full:  ts = np.linspace(0.,1.0,nt,endpoint=False)
  else:  ts = np.linspace(0.,0.5,nt,endpoint=False)
  wfall = [[occ_states2d(h,np.array([k,t,0.,])) for k in path] for t in ts] 
  # select a continuos gauge for the first wave
  for it in range(len(ts)-1): # loop over ts
    wfall[it+1][0] = smooth_gauge(wfall[it][0],wfall[it+1][0]) 
  for it in range(len(ts)): # loop over t points
    row = [] # empty list for this row
    t = ts[it] # select the t point
    wfs = wfall[it] # get set of waves 
    for i in range(len(wfs)-1):
      wfs[i+1] = smooth_gauge(wfs[i],wfs[i+1]) # transform into a smooth gauge
#      m = uij(wfs[i],wfs[i+1]) # matrix of wavefunctions
    m = uij(wfs[0],wfs[len(wfs)-1]) # matrix of wavefunctions
    evals = lg.eigvals(m) # eigenvalues of the rotation 
    x = np.angle(evals) # phase of the eigenvalues
    fo.write(str(t)+"    ") # write pumping variable
    row.append(t) # store
    for ix in x: # loop over phases
      fo.write(str(ix)+"  ")
      row.append(ix) # store
    fo.write("\n")
    out.append(row) # store
  fo.close()
  return np.array(out).transpose() # transpose the map

def _default_response_frequencies(A, n):
    """Compute a reasonable set of frequency points for bode plot.

    This function is used by `bode` to compute the frequency points (in rad/s)
    when the `w` argument to the function is None.

    Parameters
    ----------
    A : ndarray
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    w : ndarray
        The 1-D array of length `n` of frequency samples (in rad/s) at which
        the response is to be computed.
    """
    vals = linalg.eigvals(A)
    minpole = min(abs(real(vals)))
    maxpole = max(abs(real(vals)))
    # A reasonable frequency range is two orders of magnitude before the
    # minimum pole (slowest) and two orders of magnitude after the maximum pole
    # (fastest).
    w = numpy.logspace(numpy.log10(minpole) - 2, numpy.log10(maxpole) + 2, n)
    return w

def run_hamiltonian(hessian, verbose=True):
    c = ClassicalHamiltonian()

    xopt = optimize.fmin(c.potential, c.initialposition(), xtol=1e-10)

    hessian.fun = c.potential
    hessian.full_output = True

    h, info = hessian(xopt)
    true_h = np.array([[5.23748385e-12, -2.61873829e-12],
                       [-2.61873829e-12, 5.23748385e-12]])
    eigenvalues = linalg.eigvals(h)
    normal_modes = c.normal_modes(eigenvalues)

    if verbose:
        print(c.potential([-0.5, 0.5]))
        print(c.potential([-0.5, 0.0]))
        print(c.potential([0.0, 0.0]))
        print(xopt)
        print('h', h)
        print('h-true_h', np.abs(h - true_h))
        print('error_estimate', info.error_estimate)

        print('eigenvalues', eigenvalues)
        print('normal_modes', normal_modes)
    return h, info.error_estimate, true_h

def _run_hamiltonian(verbose=True):
    c = classicalHamiltonian()
    if verbose:
        print(c.potential(array([-0.5, 0.5])))
        print(c.potential(array([-0.5, 0.0])))
        print(c.potential(array([0.0, 0.0])))

    xopt = optimize.fmin(c.potential, c.initialposition(), xtol=1e-10)
    # Important to restrict the step in order to avoid the discontinutiy at
    # x=[0,0]
    # hessian = nd.Hessian(c.potential, step_max=1.0, step_nom=np.abs(xopt))
    step = nd.MaxStepGenerator(step_max=2, step_ratio=4, num_steps=16)
    hessian = nd.Hessian(c.potential, step=step, method='central',
                         full_output=True)
    # hessian = algopy.Hessian(c.potential) # Does not work
    # hessian = scientific.Hessian(c.potential) # does not work
    H, info = hessian(xopt)
    true_H = np.array([[5.23748385e-12, -2.61873829e-12],
                       [-2.61873829e-12, 5.23748385e-12]])
    if verbose:
        print(xopt)
        print('H', H)
        print('H-true_H', np.abs(H-true_H))
        print('error_estimate', info.error_estimate)

        eigenvalues = linalg.eigvals(H)
        normal_modes = c.normal_modes(eigenvalues)

        print('eigenvalues', eigenvalues)
        print('normal_modes', normal_modes)
    return H, info.error_estimate, true_H

 def test_simple_complex(self):
     a = [[1,2,3],[1,2,3],[2,5,6+1j]]
     w = eigvals(a)
     exact_w = [(9+1j+sqrt(92+6j))/2,
                0,
                (9+1j-sqrt(92+6j))/2]
     assert_array_almost_equal(w,exact_w)