def env_init(self):
        """
            Based on the levin model, the dispersion probability is initialized.
        """
        self.dispersionModel = InvasiveUtility.Levin
        notDirectedG = networkx.Graph(self.simulationParameterObj.graph)
        adjMatrix = adjacency_matrix(notDirectedG)

        edges = self.simulationParameterObj.graph.edges()
        simulationParameterObj = self.simulationParameterObj
        if self.dispersionModel == InvasiveUtility.Levin:
            parameters = InvasiveUtility.calculatePath(notDirectedG,adjMatrix, edges, simulationParameterObj.downStreamRate,
                simulationParameterObj.upStreamRate)
            C = (1 - simulationParameterObj.upStreamRate * simulationParameterObj.downStreamRate) / (
                (1 - 2 * simulationParameterObj.upStreamRate) * (1 - simulationParameterObj.downStreamRate))
            self.dispertionTable = np.dot(1 / C, parameters)
            self.germinationObj = GerminationDispersionParameterClass(1, 1)
        #calculating the worst case fully invaded rivers cost
        worst_case = repmat(1, 1, self.simulationParameterObj.nbrReaches * self.simulationParameterObj.habitatSize)[0]
        cost_state_unit = InvasiveUtility.get_unit_invaded_reaches(worst_case,
            self.simulationParameterObj.habitatSize) * self.actionParameterObj.costPerReach
        stateCost = cost_state_unit + InvasiveUtility.get_invaded_reaches(
            worst_case) * self.actionParameterObj.costPerTree
        stateCost = stateCost + InvasiveUtility.get_empty_slots(worst_case) * self.actionParameterObj.emptyCost
        costAction = InvasiveUtility.get_budget_cost_actions(repmat(3, 1, self.simulationParameterObj.nbrReaches)[0],
            worst_case, self.actionParameterObj)
        networkx.adjacency_matrix(self.simulationParameterObj.graph)
        return "VERSION RL-Glue-3.0 PROBLEMTYPE non-episodic DISCOUNTFACTOR " + str(
            self.discountFactor) + " OBSERVATIONS INTS (" + str(
            self.simulationParameterObj.nbrReaches * self.simulationParameterObj.habitatSize) + " 1 3) ACTIONS INTS (" + str(
            self.simulationParameterObj.nbrReaches) + " 1 4) REWARDS (" + str(self.Bad_Action_Penalty)+" "+str(
            -1 * (costAction + stateCost)) + ") EXTRA "+str(self.simulationParameterObj.graph.edges()) + " BUDGET "+str(self.actionParameterObj.budget) +" by Majid Taleghan."

def dsd_matrix(G, nodeList, npyFile, LMsetSize=50, centralityFunc='degree', **kwargs):
    """
    any kwargs, if specified, will be passed into centrality function call.
    """
    # if npy path not entered, or file does not exist, compute D
    if not npyFile or not os.path.isfile(npyFile):
        #
        # construct hemat
        #
        adjMatrix = np.array(nx.adjacency_matrix(G,nodeList))
        if np.shape(adjMatrix) == ():
            adjMatrix = np.array(nx.adjacency_matrix(G,nodeList).todense())
        HEmatrix = dsd.hematrix(adjMatrix)

        # construct DSD
        LMset = get_LMset(G, nodeList, LMsetSize, centralityFunc, **kwargs)
        D = dsd.DSD(HEmatrix,LMset)

        if npyFile:
            try:
                np.save(npyFile, D)
            except IOError:
                os.makedirs(npyFile[:npyFile.rfind('/')])
                np.save(npyFile, D)
    # otherwise just load and return it
    else:
        D = np.load(npyFile)
    return D

def nsim_bvd04(G1, G2, max_iter=100, eps=1e-4):
    """
    Algorithm to calculate node-node similarity matrix of
    two directed graphs.

    Return
    ------
    A 2d similarity matrix of |V1| x |V2|.

    Reference
    ---------
    Blondel, Vincent D. et al. "A Measure of Similarity between Graph Vertices:
    Applications to Synonym Extraction and Web Searching." SIAM Review (2004)
    """
    N = len(G1.nodes())
    M = len(G2.nodes())
    A = nx.adjacency_matrix(G1).todense()
    B = nx.adjacency_matrix(G2).todense()
    nsim_prev = np.zeros((M, N))
    nsim = np.ones((M, N))

    for i in range(max_iter):
        if np.allclose(nsim, nsim_prev, atol=eps):
            break

        nsim_prev = np.copy(nsim)
        nsim = np.dot(np.dot(B, nsim_prev), A.T) + \
            np.dot(np.dot(B.T, nsim_prev), A)

        fnorm = np.linalg.norm(nsim, ord='fro')
        nsim = nsim / fnorm

    print("Converge after %d iterations (eps=%f)." % (i, eps))

    return nsim.T

def are_isomorphic(G, H):
    """Check whether two graphs G and H are isomorphic.
    
    Note: This function is brute force and very slow.
    
    args:
        G: a networkx Graph
        H: a networkx Graph
    
    returns:
    	True if G and H are isomorphic.
    	False if G and H are not isomorphic.
    """
    n = len(G.nodes())
    m = len(H.nodes())
    if n != m:
        return False
    if sorted(G.degree().values()) != sorted(H.degree().values()):
        return False
    else:
        a_g = nx.adjacency_matrix(G)
        vertex_perms = list(permutations(H.nodes(), m))
        for i in vertex_perms:
            a_h = nx.adjacency_matrix(H, i)
            if (a_h == a_g).all():
                #print(list(zip(G.nodes(), i)), "is an isomorphism") 
                return True
        return False

def lesion_graph_degree_thresh(graph, threshold):
    """
    Remove vertices from a graph with degree greater than or equal to
    threshold.

    Parameters:
    -----------
        graph: NetworkX graph to be lesioned.
        threshold: Degree above which to remove nodes.

    Returns:
    --------
        G: NetworkX graph
        A: Adjacency matrix for graph

    """
    # Error checking
    G = deepcopy(graph)

    assert threshold >= 0, " In percolation, `threshold` must be >= 0."
    # Check if lesioning is necessary for threshold
    if threshold > max(G.degree().values()):
        return G, nx.adjacency_matrix(G)

    # Identify all node indices >= threshold
    node_inds = np.where(np.asarray(G.degree().values()) >= threshold)[0]
    # Eliminate these nodes
    G.remove_nodes_from(node_inds)

    if G.order() > 0:
        return G, nx.adjacency_matrix(G)
    else:
        #print 'Graph completely lesioned.'
        return None, None

    def similarity(self, G=None, H=None, iters=1000):
        """ Returns the graph similarity based on

        :param G: networkx graph of original graph (default: self.G)
        :param H: networkx graph of inferred graph (default: self.H)
        :param iter: number of iterations (default: 20)
        :return: float

        """
        if G is None:
            G = self.G
        if H is None:
            H = self.H

        n = len(G)

        gA = nx.adjacency_matrix(G)
        hA = nx.adjacency_matrix(H)
        s = np.identity(n)   # initial condition
        for i in range(iters):
            temp = (np.kron(gA, hA) + np.kron(gA.T, hA.T)) * s + \
                np.identity(n) * 0.0000001
            s = temp / np.linalg.norm(temp)

        a = np.trace(s)

        temp = (np.kron(gA, hA) + np.kron(gA.T, hA.T)) * s
        s = temp / np.linalg.norm(temp)

        a += np.trace(s)
        return a / 2

def lesion_graph_degree(graph, num_lesions):
    """
    Remove vertices from a graph according to degree.

    Args:
        graph: NetworkX graph to be lesioned.
        num_lesions: Number of top degree nodes to remove.

    Returns:
        G: NetworkX graph
        A: Adjacency matrix for graph

    """
    # Error checking
    G = deepcopy(graph)
    if num_lesions == 0:
        return G, nx.adjacency_matrix(G)

    assert num_lesions >= 0 and num_lesions < graph.order, 'Attempting to\
        remove too many/few nodes'

    for l in range(num_lesions):
        # Identify node to cut
        node_i, node_d = max(G.degree().items(),
                             key=lambda degree: degree[1])
        G.remove_node(node_i)
        #print (node_i, node_d)

    if G.order() > 0:
        return G, nx.adjacency_matrix(G)
    else:
        print 'Graph completely lesioned.'
        return None, None

 def update_nl(self, dist):
     self.G.remove_edges_from(self.G.edges())
     print "dr:"
     print self.dr()
     for row in list(enumerate(self.dr())):
         for col in list(enumerate(row[1])):
             if col[1] > dist:
                 self.G.add_edge(col[0], row[0], x=self.x[col[0]]-self.x[row[0]], y=self.y[col[0]]-self.y[row[0]], z=0., r=col[1])
     print nx.adjacency_matrix(self.G)

 def createMultiplex(self):
     print self.politicsGraph.nodes().__len__()
     sortedPoliticalComentators = sorted(self.politicsGraph.nodes())
     politicalComentatorsAdjMat = nx.adjacency_matrix(self.politicsGraph, sortedPoliticalComentators)
     sortedChurchComentators = sorted(self.churchGraph.nodes())
     churchComentatorsAdjMat = nx.adjacency_matrix(self.churchGraph, sortedChurchComentators)
     otherComentators = sorted(self.othersGraph.nodes())
     otherAdjMat = nx.adjacency_matrix(self.othersGraph, otherComentators)
     self.addLayerToGraph(politicalComentatorsAdjMat, sortedPoliticalComentators, 'L1', 1)
     self.addLayerToGraph(churchComentatorsAdjMat, sortedChurchComentators, 'L2', 2)
     self.addLayerToGraph(otherAdjMat, otherComentators, 'L3', 3)

def show_graph_features(G,circular=False):
    for node,adj_list in G.adjacency_iter():
        print node, ': ',adj_list   
        
    print  nx.adjacency_matrix(G)
    
    #N\u00f3s:
    print 'N\u00famero de n\u00f3s: ',G.number_of_nodes()
    print 'N\u00f3s: \\n','\\t',G.nodes()  
    #Arestas:
    print 'N\u00famero de arestas: ',G.number_of_edges()
    print 'Arestas: \\n','\\t',G.edges(data=True)

def graph_diff(g0, g1):
    m0 = nx.adjacency_matrix(g0)
    m1 = nx.adjacency_matrix(g1)
    delta = 0
    for i in range(0, m0.shape[0]):
        for j in range(0, m0.shape[1]):
            if i > j:
                # print m0[(i,j)],int(m1[(i,j)]>0)
                # if(m0[(i,j)]==1 and int(m1[(i,j)]>0)!=1):
                if m0[(i, j)] != int(m1[(i, j)] > 0):
                    delta += 1
    return delta

def cal_exact_Nd_simple(H, random_weight=False):
    """return (Nd, N-Rank) """
    G = H.copy()
    N = H.number_of_nodes()
    try:
        G2 = nx.adjacency_matrix(G, weight = 'weight')
        #~ print 'weight'
    except:
        G2 = nx.adjacency_matrix(G)
    if random_weight:
        G2 = np.array(G2)*np.random.random((N, N))
    rank_G = mranksvd(G2)
    return max(1, N-rank_G), N-rank_G

 def GTrieInsert(self, graph, label=None, states=False):
     if not self.root.isLeaf() and self.null_graph:
         self.insertRecursive(networkx.Graph(), [], networkx.adjacency_matrix(networkx.Graph()).todense(),
                              self.root, 0, label, states)
     components = networkx.connected_components(graph.to_undirected()) \
         if networkx.is_directed(graph) else networkx.connected_components(graph)
     component_len = [1 for x in components if len(x) > 1]
     if len(list(components)) > 1 and sum(component_len) > 1:
         print "Illegal Graph Insert: Graph has more than one connnected component."
         return
     cannonGraph = self.GTCannon(graph.copy())
     matrix = networkx.adjacency_matrix(cannonGraph).todense()
     conditions = self.utility.symmetryConditions(cannonGraph)
     self.insertRecursive(cannonGraph, conditions, matrix, self.root, 0, label, states)

def init():
    global op1, op2, grouping, network, adj_mat
    op1 = [(-1 + 2 * random.random()) for agents in range(no_of_agents)]
    op2 = [(-1 + 2 * random.random()) for agents in range(no_of_agents)]
    grouping = [random.randint(1,m) for agents in range(no_of_agents)]
    network = nx.erdos_renyi_graph(no_of_agents,p)
    adj_mat = nx.adjacency_matrix(network)

 def test_eigenvector_v_katz_random(self):
     G = nx.gnp_random_graph(10,0.5, seed=1234)
     l = float(max(eigvals(nx.adjacency_matrix(G).todense())))
     e = nx.eigenvector_centrality_numpy(G)
     k = nx.katz_centrality_numpy(G, 1.0/l)
     for n in G:
         assert_almost_equal(e[n], k[n])

 def busmap_by_spectral_clustering(network, n_clusters, **kwds):
     lines = network.lines.loc[:,['bus0', 'bus1']].assign(weight=1./network.lines.x).set_index(['bus0','bus1'])
     G = OrderedGraph()
     G.add_nodes_from(network.buses.index)
     G.add_edges_from((u,v,dict(weight=w)) for (u,v),w in lines.itertuples())
     return pd.Series(sk_spectral_clustering(nx.adjacency_matrix(G), n_clusters, **kwds) + 1,
                      index=network.buses.index)

def _get_node_plot_props(G, node_class=None, max_energy=None, active_node_color=None, active_edge_color=None,
                         dead_node_color=None):
    """
    `node_class` - Generic | Internal | Sensory | Motor
    `node_size` - proportional to the sum of the presynaptic connections it makes with other nodes.
    `node_colors` - function of excitatory/inhibitory, energy_value, firing/inactive

    """
    cm = CMAP_DIFF  # Shade from red (inhibitory) to green (excitatory)
    nodes = G.nodes(node_class)
    adj_matrix = nx.adjacency_matrix(G)
    node_pos = nx.get_node_attributes(G.subgraph(nodes), 'pos')
    edge_width = np.array([d['weight'] for (u, v, d) in G.edges(data=True) if u in nodes])
    firing_nc = colors.hex2color(active_node_color) if active_node_color is not None \
        else list(colors.hex2color(RENDER_NODE_PROPS['Firing']['node_face_color']))
    dead_nc = colors.hex2color(dead_node_color) if dead_node_color is not None \
        else list(colors.hex2color(RENDER_NODE_PROPS['Dead']['node_face_color']))
    _ = firing_nc.append(1.)
    _ = dead_nc.append(1.)

    node_colors = _get_node_colors(G, cm, node_class=node_class, max_energy=max_energy, firing_node_color=firing_nc,
                                   dead_node_color=dead_nc)

    if node_class is not None:
        min_ns, max_ns = RENDER_NODE_PROPS[node_class]['min_node_size'], RENDER_NODE_PROPS[node_class]['max_node_size']
        node_shape = RENDER_NODE_PROPS[node_class]['shape']
        node_size = np.array([np.maximum(adj_matrix[i].sum(), .01) for i, n_id in enumerate(G.nodes())
                              if G.node[n_id]['node_class'] == node_class])  # proportional to the number of connections
    else:
        node_shape, node_size = RENDER_NODE_PROPS['Default']['shape'], adj_matrix.sum(axis=1)
        min_ns, max_ns = RENDER_NODE_PROPS['Default']['min_node_size'], RENDER_NODE_PROPS['Default']['max_node_size']
    node_size = min_ns + (max_ns - min_ns) * (node_size - node_size.min()) / (node_size.max() - node_size.min()) \
        if node_size.max() > node_size.min() else max_ns * np.ones_like(node_size)
    return node_pos, node_colors, node_shape, node_size, edge_width

def hits(G,max_iter=100,tol=1.0e-6):
    M=nx.adjacency_matrix(G,nodelist=G.nodes())
    (n,m)=M.shape # should be square
    A=M.T*M # authority matrix
    x=scipy.ones((n,1))/n  # initial guess
    # power iteration on authority matrix
    i=0
    while True:
        xlast=x
        x=A*x
        x=x/x.sum()
        # check convergence, l1 norm
        err=scipy.absolute(x-xlast).sum()
        if err < tol:
            break
        if i>max_iter:
            raise NetworkXError(\
            "HITS: power iteration failed to converge in %d iterations."%(i+1))
        i+=1

    a=np.asarray(x).flatten()
    h=np.asarray(M*a).flatten()
    hubs=dict(zip(G.nodes(),h/h.sum()))
    authorities=dict(zip(G.nodes(),a/a.sum()))
    return hubs,authorities

    def loglikelihood(self, net):
        """Returns the log likelihood of the given network.
        
        Similar to the loglikelihood method of a Conditional Probability
        Distribution.  
        
        """

        adjmat = nx.adjacency_matrix(net)

        # if any of the mustexist or mustnotexist constraints are violated,
        # return negative infinity
        if (not (adjmat | self.mustexist).all()) or \
           (adjmat & self.mustnotexist).any():
            return NEGINF

        # if any custom constraints are violated, return negative infinity
        if self.constraints and not all(c(adjmat) for c in self.constraints):
            return NEGINF

        loglike = 0.0
        if self.energy_matrix != None:
            energy = N.sum(adjmat * self.energy_matrix) 
            loglike = -self.weight * energy

        return loglike

def coarse_grain_W(num_intervals, num_eigenvectors, g, sparse = True):
    """ Produces W_tilde := R*W*K, where W is the stochastic
        matrix of the original graph, and R,K, are intermediary
        matrices defined in the paper.

        Has an optional arguments to use non-sparse matrices, which are
        (minorly) faster for small graphs.
    """

    
    A = nx.adjacency_matrix(g)
    num_nodes = A.shape[0]
    A = A / np.sum(A, 0)    #stochastic matrix -- don't need A anymore
    A = np.nan_to_num(A)
    eigenvalues,left_eigenvectors = eig(A, left = True, right = False)
    
    if sparse == True:
        A = s.csr_matrix(A)
        groups = make_groups(eigenvalues, left_eigenvectors, num_intervals , num_eigenvectors)
        R = make_sparse_R(groups, num_nodes)
        K = make_sparse_K(groups, num_nodes, g)
        return np.dot(R, np.dot(A, K))
    
    else:
        groups = make_groups(eigenvalues, left_eigenvectors, num_intervals , num_eigenvectors)
        R = make_R(groups, num_nodes)
        K = make_K(groups, num_nodes, g)
        return np.dot(R, np.dot(A, K))