def diffusion_matrix(G, nodeList, npyFile, gamma=8):
    """
    compute inverse of Laplacian matrix and symmetrize. (default gamma is arbitrary)
    the higher the influence, the closer we expect function to be, so let distances be reciprocal.
    """
    if not npyFile or not os.path.isfile(npyFile):

        L = np.array(nx.laplacian_matrix(G,nodeList))
        # depending on version of networkx, might get sparse matrix instead. if so, do this:
        if np.shape(L) == ():
            L = np.array(nx.laplacian_matrix(G,nodeList).todense())
        m, n = np.shape(L)
        L = L + (np.eye(m, n)*gamma)
        D = np.linalg.inv(L)
        n = len(nodeList)
        for i in xrange(n):
            for j in xrange(i+1,n):
                D[i][j] = D[j][i] = 1/(min(D[i][j], D[j][i]))

        if npyFile:
            np.save(npyFile, D)
    else:
        D = np.load(npyFile)

    return D

def laplacian_spectrum(G, weight="weight"):
    """Return eigenvalues of the Laplacian 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
    --------
    laplacian_matrix
    """
    try:
        import numpy as np
    except ImportError:
        raise ImportError("laplacian_spectrum() requires NumPy: http://scipy.org/ ")
    return np.linalg.eigvals(nx.laplacian_matrix(G, weight=weight))

def laplacian_spectrum(G, weight='weight'):
    """Return eigenvalues of the Laplacian 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
    --------
    laplacian_matrix
    """
    from scipy.linalg import eigvals
    return eigvals(nx.laplacian_matrix(G,weight=weight).todense())

def eig_vis_opt(G, F, beta):
	"""
		Computes first and second eigenvector of sqrt(C+beta*L)^T CAC sqrt(C+beta*L) matrix for visualization.
		Input:
			* G: graph
			* F: graph signal
			* beta: regularization parameter
		Output:
			* v1: first eigenvector
			* v2: second eigenvector	
	"""
	ind = {}
	i = 0
	
	for v in G.nodes():
		ind[v] = i
		i = i + 1
	
	C = laplacian_complete(networkx.number_of_nodes(G))
	A = weighted_adjacency_complete(G, F, ind)
	CAC = numpy.dot(numpy.dot(C,A), C)
	L = networkx.laplacian_matrix(G).todense()
	
	isqrtCL = sqrtmi( C + beta * L)
	M = numpy.dot(numpy.dot(isqrtCL, CAC), isqrtCL)
	
	(eigvals, eigvecs) = scipy.linalg.eigh(M,eigvals=(0,1))
	x1 = numpy.asarray(numpy.dot(eigvecs[:,0], isqrtCL))[0,:]
	x2 = numpy.asarray(numpy.dot(eigvecs[:,1], isqrtCL))[0,:]

	return x1, x2

    def etape(self, frontier, i_graph):
        """ Calculates the most probable seed within the infected nodes"""

        # Taking the actual submatrix, not the laplacian matrix. The change
        # lies in the total number of connections (The diagonal terms) for the
        # infected nodes connected to uninfected ones in the initial graph
        i_laplacian_matrix = nx.laplacian_matrix(i_graph)
        for i in range(0, len(i_graph.nodes())):
            if frontier.has_node(i_graph.nodes()[i]):
                i_laplacian_matrix[i, i] +=\
                                frontier.node[i_graph.nodes()[i]]['clear']

        # SymPy
        Lm = Matrix(i_laplacian_matrix.todense())
        i = self.Sym2NumArray(Matrix(Lm.eigenvects()[0][2][0])).argmax()

        # NumPy
        # val, vect = linalg.eigh(i_laplacian_matrix.todense())
        #i = vect[0].argmax()

        # SciPY
        # val, vect = eigs(i_laplacian_matrix.rint())
        # i = vect[:, 0].argmax()

        seed = (i_graph.nodes()[i])

        return seed

def effective_resistance_project(G, beacons):
    from numpy.linalg import pinv
    projection = np.zeros((G.number_of_nodes() - len(beacons), len(beacons)))
    L = nx.laplacian_matrix(G)
    B = nx.incidence_matrix(G).T
    B_e = B.copy()
    
    L_pseudo = pinv(L)
    for i in xrange(B.shape[0]):
        min_ace = np.min(np.where(B[i,:] ==1)[1])
        B_e[i, min_ace] = -1
    
    for i,beacon in enumerate(beacons):
        node_index = 0
        for j,node in enumerate(G.nodes()):
            if node in beacons:
                continue
                
            battery = np.zeros((B_e.shape[1],1))
            battery[i] = 1
            battery[node_index] = -1

            p = L_pseudo * battery
            projection[node_index][i] = abs(p[i] - p[j])
            node_index += 1 
    return projection

def MinCut(G):

	#Calcula a matriz laplaciana do grafo G
	#Opcionalmente, pode-se usar a laplaciana normalizada
	#lap = nx.normalized_laplacian(G)
	lap = nx.laplacian_matrix(G)
	eigenValues, eigenVectors = la.eigh(lap)

	orthoVector = []
	
	#pega-se entao os componentes orthonormais dos
	#autovetores e cria-se um novo vetor
	for vectors in eigenVectors:
		orthoVector.append(vectors[1])

	#para o Ratio-cut, usa-se a mediana para dividir
	#o grafo.
	#med = np.median(eigenVectors[1])

	nodesleft = []
	nodesright = []

	#divide-se entao o grafo em 2 componentes, baseado no sinal
	#do vetor orthonormal. Compara-se a lista de nodos com o vetor.
	#Se o valor for maior que zero, vai pra uma componente, caso contrario,
	#vai pra outra.
	for node, vec in zip(G.nodes(), orthoVector):
		if(vec > 0):
			nodesleft.append(node)
		else:
			nodesright.append(node)
	
	return (nodesleft, nodesright)

def cheb_spectral_cut(CAC, start, F, G, beta, k, n, ind):
	"""
		Fast spectral cut implementation using chebyshev polynomials.
		Input:
			* CAC: C*A*C where C is the Laplacian of a complete graph and A is a pairwise squared difference matrix
			* start: initialization
			* F: graph signal
			* G: graph
			* L: graph laplacian matrix
			* beta: regularization parameter
			* k: max edges cut
			* n: number of polynomials
			* ind: vertex index vertex: unique integer
		Output:
			* res: dictionary with following fields:
				- x: indicator vector
				- size: number of edges cut
				- score: cut score
				- energy: cut energy
	"""
	L = networkx.laplacian_matrix(G)
	M = chebyshev_approx_2d(n, beta, CAC, L)
	
	eigvec = power_method(-M, start, 10)
	x = chebyshev_approx_1d(n, beta, eigvec, L)
	
	(x, score, size, energy) = sweep_opt(x, beta, F, G, k, ind)
	
	res = {}
	res["x"] = numpy.array(x)
	res["size"] = size
	res["score"] = score
	res["energy"] = energy
	
	return res

def create_laplacian_matrix(G):
	row = []
	column = []
	value = []

	for t in range(G.num_snaps()):
		Lg = networkx.laplacian_matrix(G.snap(t))
		for (i,j) in zip(*scipy.nonzero(Lg)):
			row.append(G.size()*t + i)
			column.append(G.size()*t + j)
			
			if i != j:
				value.append(Lg[i,j])
			else:
				if t > 0 and t < G.num_snaps() - 1:
					value.append(Lg[i,j] + 2 * G.swap_cost())
				else:	
					value.append(Lg[i,j] + 1 * G.swap_cost())
	
	for t in range(G.num_snaps()-1):
		for v in range(G.size()):
			row.append(t*G.size() + v)
			column.append((t+1)*G.size() + v)
			value.append(-1 * G.swap_cost())
			
			column.append(t*G.size() + v)
			row.append((t+1)*G.size() + v)
			value.append(-1 * G.swap_cost())


	sz = G.num_snaps() * G.size()
	return scipy.sparse.csr_matrix((value, (row, column)), shape=(sz, sz), dtype=float)

 def test_abbreviation_of_method(self):
     G = nx.path_graph(8)
     A = nx.laplacian_matrix(G)
     sigma = 2 - sqrt(2 + sqrt(2))
     ac = nx.algebraic_connectivity(G, tol=1e-12, method='tracemin')
     assert_almost_equal(ac, sigma)
     x = nx.fiedler_vector(G, tol=1e-12, method='tracemin')
     check_eigenvector(A, sigma, x)

def eig_vis_nc(G):
	L = networkx.laplacian_matrix(G).todense()
	(eigvals, eigvecs) = scipy.linalg.eigh(L,eigvals=(1,2))

	x1 = numpy.asarray(eigvecs[:,0])
	x2 = numpy.asarray(eigvecs[:,1])
	
	return x1, x2

def ematrix(G,nodes):
    L = np.array(nx.laplacian_matrix(G,nodes))
    # not using eig.  gives very strange results for non-invertible matrices
    # which don't agree with matlab's eig.  svd however gives sensible results
    u,s,vt = np.linalg.svd(L)
    s = 1/s
    s[np.where(np.abs(s)>1e9)]=0
    s = np.sqrt(s)
    return np.dot(u,np.diag(s))

def my_algebraic_connectivity(graph, normalise=False):
    if normalise:
        eigvals, eigvecs = sp.sparse.linalg.eigsh(nx.normalized_laplacian_matrix(graph).asfptype(), 2, which='SA')
        a = eigvals[1]
    else:
        eigvals, eigvecs = sp.sparse.linalg.eigsh(nx.laplacian_matrix(graph).asfptype(), 2, which='SA')
        a = eigvals[1]
    if a < MACHINE_EPSILON: a = 0.0
    return a

 def test_cycle(self):
     G = nx.cycle_graph(8)
     A = nx.laplacian_matrix(G)
     sigma = 2 - sqrt(2)
     for method in self._methods:
         ac = nx.algebraic_connectivity(G, tol=1e-12, method=method)
         assert_almost_equal(ac, sigma)
         x = nx.fiedler_vector(G, tol=1e-12, method=method)
         check_eigenvector(A, sigma, x)

 def test_two_nodes(self):
     G = nx.Graph()
     G.add_edge(0, 1, weight=1)
     A = nx.laplacian_matrix(G)
     for method in self._methods:
         assert_almost_equal(nx.algebraic_connectivity(
             G, tol=1e-12, method=method), 2)
         x = nx.fiedler_vector(G, tol=1e-12, method=method)
         check_eigenvector(A, 2, x)
     G = nx.MultiGraph()
     G.add_edge(0, 0, spam=1e8)
     G.add_edge(0, 1, spam=1)
     G.add_edge(0, 1, spam=-2)
     A = -3 * nx.laplacian_matrix(G, weight='spam')
     for method in self._methods:
         assert_almost_equal(nx.algebraic_connectivity(
             G, weight='spam', tol=1e-12, method=method), 6)
         x = nx.fiedler_vector(G, weight='spam', tol=1e-12, method=method)
         check_eigenvector(A, 6, x)

def write_adj_mat(G, fileobj=sys.stdout):
    """Write G to a sparse matrix format that Julia and Matlab can read."""
    lapmatrix = nx.laplacian_matrix(G)
    norm_lapl = nx.normalized_laplacian_matrix(G)
    adjmatrix = nx.adjacency_matrix(G)
    mdict = {'laplacian': lapmatrix,
             'norm_lapl': norm_lapl,
             'adjacency': adjmatrix}
    sio.savemat(fileobj, mdict)
    return mdict

 def test_problematic_graph_issue_2381(self):
     G = nx.path_graph(4)
     G.add_edges_from([(4, 2), (5, 1)])
     A = nx.laplacian_matrix(G)
     sigma = 0.438447187191
     for method in self._methods:
         ac = nx.algebraic_connectivity(G, tol=1e-12, method=method)
         assert_almost_equal(ac, sigma)
         x = nx.fiedler_vector(G, tol=1e-12, method=method)
         check_eigenvector(A, sigma, x)

def one_d_search(G, F, k, ind):
	"""
		Cut computation. Perform 1-D search for beta using golden search.
		Input:
			* G: graph
			* F: graph signal
			* k: max edges to be cut
			* n: number of chebyshev polynomials
			* ind: vertex index vertex: unique integer
		Output:
			* cut
	"""
	C = laplacian_complete(networkx.number_of_nodes(G))
	A = weighted_adjacency_complete(G,F, ind)
	CAC = numpy.dot(numpy.dot(C,A), C)
	start = numpy.ones(networkx.number_of_nodes(G))
	L = networkx.laplacian_matrix(G).todense()

	#Upper and lower bounds for search
	a = 0.
	b = 1000.
	c=b-gr*(b-a)
	d=a+gr*(b-a)
	
	#Tolerance
	tol = 1.

	resab = {}
	resab["size"] = k + 1
	
	#golden search
	while abs(c-d)>tol or resab["size"] > k:      
		resc = spectral_cut(CAC, L, C, A, start, F, G, c, k, ind)
		resd = spectral_cut(CAC, L, C, A, start, F, G, d, k, ind)
		
		if resc["size"] <= k: 
			if resc["score"] > resd["score"]: 
				start = numpy.array(resc["x"])
				b = d
				d = c
				c=b-gr*(b-a)
			else:
				start = numpy.array(resd["x"])
				a=c
				c=d  
				d=a+gr*(b-a)
		else:
				start = numpy.array(resc["x"])
				a=c
				c=d  
				d=a+gr*(b-a)
		
		resab = spectral_cut(CAC, L, C, A, start, F, G, (b+a) / 2, k, ind)
	
	return resab

def spectral_clustering(A, args):
    '''
    Function uses dimension reduction by spectral clustering 
    to cluster with some cluster algorithm 
    '''
    if type(A) == np.ndarray:
        print 'Building graph...'
        G = nx.from_numpy_matrix(A)
        laplacian = nx.laplacian_matrix(G)
    elif type(A) == scsp.csc.csc_matrix:
        print 'Building graph...'
        G = nx.from_scipy_sparse_matrix(A)
        laplacian = nx.laplacian_matrix(G)
    print 'Start finding reduction...'
    max_eigvals = scspl.eigs(laplacian * 1.0, return_eigenvectors=False, 
                             k=args['dim'] + 10)
    alpha = np.amax(max_eigvals)
    modif_laplacian = alpha * scsp.eye(laplacian.shape[0], 
                                       laplacian.shape[1]) - laplacian
    lap_eigval, lap_eigvec = scspl.eigs(modif_laplacian * 1.0, 
                                        k=20)
    if args['show_eigenvalues']:
        print 'Show eigenvalues...'
        plt.scatter(np.array(range(len(lap_eigval))) + 1, 
                    np.sort(-lap_eigval.real + alpha, reverse = True), 
                    marker='.', s=90)
        plt.show()
    U = lap_eigvec[:, 0:args['dim']].real
    print 'Start clustering...'
    if U.shape[0] < 1000:
        clustering = skcl.KMeans(n_clusters=args['number_of_clusters'])
        clustering.fit(U)
    else:
        clustering = skcl.MiniBatchKMeans(n_clusters=args['number_of_clusters'])
        clustering.fit(U)
    center_person_cluster_id = clustering.labels_[args['index']]
    center_person_cluster = [i for i in xrange(len(clustering.labels_)) 
                             if clustering.labels_[i] == center_person_cluster_id]
    
    cluster = [[i for i in xrange(len(clustering.labels_)) if clustering.labels_[i] == j] 
                for j in xrange(args['number_of_clusters'])]
    return center_person_cluster, cluster

def eigenratio (G):
    if nx.is_connected(G):
        # Calculate the eigenrato (lambda_N/lambda_2)
        L = nx.laplacian_matrix(G)
        eigenvalues, eigenvectors = np.linalg.eig(L)
        idx = eigenvalues.argsort()   
        eigenvalues = eigenvalues[idx]
        return eigenvalues[-1] / eigenvalues[1]
    else:
        # If the network is not connected it is not valid
        return float('inf')