def cut_size(G, S, T=None, weight=None):
    """Returns the size of the cut between two sets of nodes.

    A *cut* is a partition of the nodes of a graph into two sets. The
    *cut size* is the sum of the weights of the edges "between" the two
    sets of nodes.

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

    S : sequence
        A sequence of nodes in ``G``.

    T : sequence
        A sequence of nodes in ``G``. If not specified, this is taken to
        be the set complement of ``S``.

    weight : object
        Edge attribute key to use as weight. If not specified, edges
        have weight one.

    Returns
    -------
    number
        Total weight of all edges from nodes in set ``S`` to nodes in
        set ``T`` (and, in the case of directed graphs, all edges from
        nodes in ``T`` to nodes in ``S``).

    Examples
    --------
    In the graph with two cliques joined by a single edges, the natural
    bipartition of the graph into two blocks, one for each clique,
    yields a cut of weight one::

        >>> G = nx.barbell_graph(3, 0)
        >>> S = {0, 1, 2}
        >>> T = {3, 4, 5}
        >>> nx.cut_size(G, S, T)
        1

    Each parallel edge in a multigraph is counted when determining the
    cut size::

        >>> G = nx.MultiGraph(['ab', 'ab'])
        >>> S = {'a'}
        >>> T = {'b'}
        >>> nx.cut_size(G, S, T)
        2

    Notes
    -----
    In a multigraph, the cut size is the total weight of edges including
    multiplicity.

    """
    edges = nx.edge_boundary(G, S, T, data=weight, default=1)
    if G.is_directed():
        edges = chain(edges, nx.edge_boundary(G, T, S, data=weight, default=1))
    return sum(weight for u, v, weight in edges)

 def test_null_graph(self):
     null = nx.null_graph()
     assert_equal(list(nx.edge_boundary(null, [])), [])
     assert_equal(list(nx.edge_boundary(null, [], [])), [])
     assert_equal(list(nx.edge_boundary(null, [1, 2, 3])), [])
     assert_equal(list(nx.edge_boundary(null, [1, 2, 3], [4, 5, 6])), [])
     assert_equal(list(nx.edge_boundary(null, [1, 2, 3], [3, 4, 5])), [])

 def test_null_edge_boundary(self):
     """null nxgraph has empty edge boundaries"""
     null=self.null
     assert_equal(nx.edge_boundary(null,[]),[])
     assert_equal(nx.edge_boundary(null,[],[]),[])
     assert_equal(nx.edge_boundary(null,[1,2,3]),[])
     assert_equal(nx.edge_boundary(null,[1,2,3],[4,5,6]),[])
     assert_equal(nx.edge_boundary(null,[1,2,3],[3,4,5]),[])

 def test_directed(self):
     """Tests the edge boundary of a directed graph."""
     G = nx.DiGraph([(0, 1), (1, 2), (2, 3), (3, 4), (4, 0)])
     S = {0, 1}
     boundary = list(nx.edge_boundary(G, S))
     expected = [(1, 2)]
     assert_equal(boundary, expected)

 def test_multigraph(self):
     """Tests the edge boundary of a multigraph."""
     G = nx.MultiGraph(list(nx.cycle_graph(5).edges()) * 2)
     S = {0, 1}
     boundary = list(nx.edge_boundary(G, S))
     expected = [(0, 4), (0, 4), (1, 2), (1, 2)]
     assert_equal(boundary, expected)

def island_update(topology, new_policy):
    """
    precondition: Assumes that only one island update is performed,
    and no subspace updates have been performed. This assumption is
    forced by our use of VLAN tags instead of MPLS labels
    provides per-packet
    """
    inst.stats.tally_update(new_policy)
    log.info("Island update")

    old_policy = inst.current_abstract_policy

    # Switches which didn't change in new policy
    nops = set( s1 for s1, c1 in old_policy \
                if switch_covered(c1, new_policy[s1]))
    # Everything else
    new =  set(topology.switches()) - nops
    old = set()

    fixpoint = island_fixpoint(topology, new_policy)
    while new:
        additions = fixpoint(new, old)
        old |= new
        new = additions
        
    mods = old

    subpolicy = restrict_policy(mods, new_policy)

    boundary = nx.edge_boundary(topology, mods)
    fake_edge_ports = \
      [topology.node[x]['ports'][y] for (x,y) in boundary \
       if topology.node[y]['isSwitch']]

    # retrieve current data from inst        
    current_internal_policy = inst.current_internal_policy
    current_edge_policy = inst.current_edge_policy
    current_version = inst.current_version
    current_priority = inst.current_priority

    # calculate new version and priority
    new_version = current_version + 1
    new_priority = current_priority - 1

    # Have to manually construct the correct edge policies by
    # distinguishing between "true" edge ports to hosts and "fake"
    # edge ports to other switches running the old version.
    internal_policy, edge_policy = \
      mk_versioned_policies(subpolicy, new_version, new_priority, topology,
                            old_version=current_version, 
                            fake_edge_ports=fake_edge_ports)

    old_internal_policy = restrict_policy(mods, current_internal_policy)
    old_edge_policy = restrict_policy(mods, current_edge_policy)
    
    return UpdateObject(internal_policy, edge_policy,
                        old_internal_policy,
                        old_edge_policy,
                        new_priority, new_version)

 def test_multidigraph(self):
     """Tests the edge boundary of a multdiigraph."""
     edges = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 0)]
     G = nx.MultiDiGraph(edges * 2)
     S = {0, 1}
     boundary = list(nx.edge_boundary(G, S))
     expected = [(1, 2), (1, 2)]
     assert_equal(boundary, expected)

def find_bridge_edges(raw, squelched, groupsize):
    groups = build_affinity_groups(squelched, groupsize)
    for (left, right) in itertools.combinations(groups, 2):
        left_name = name_subgraph(left)
        right_name = name_subgraph(right)
        boundary = nx.edge_boundary(raw, left.nodes(), right.nodes())
        if boundary:
            yield left_name, right_name, boundary

	def check_state(self, min_size=None):

		# check nodes of G are all in partition assignment
		V = set(self.g.nodes())
		nodes = set(self.nodes.keys())
		assert V == nodes

		# check each nodes's partition contains the node
		assigned_partitions = set([])
		for u in V:
			i = self.nodes[u]
			assert u in set(self.partition[i])
			assigned_partitions.add(i)

		assert assigned_partitions == set(self.partition.keys())

		# check that the set of nodes in partitioning is exactly V
		tmp = reduce(lambda x,y: x+y, self.partition.values())
		nodes = set(tmp)
		assert len(tmp) == len(nodes)
		assert V == nodes

		# check min size of partitions
		partition_sizes = map(len, self.partition.values())
		if min_size != None:
			assert min(partition_sizes) >= min_size
		assert sum(partition_sizes) == len(V)

		assert set(self.partition_graph.nodes()) == set(self.partition.keys())

		for i in self.partition_graph:
			for j in self.partition_graph:

				if i < j:
					continue

				if i == j:
					nodes_i = set(self.partition[i])
					neighbors = []
					for u in nodes_i:
						neighbors += self.g.neighbors(u)
					count = 0
					for nbr in neighbors:
						count += nbr in nodes_i
					count /= 2  # each edge is double counted
				else: 
					nodes_i = self.partition[i]
					nodes_j = self.partition[j]
					count = len(networkx.edge_boundary(self.g, nodes_i, nodes_j))

				if self.partition_graph.has_edge(i,j):
					stored_count = self.partition_graph.get_edge_data(i,j,key=0)["count"]
					assert count == stored_count, \
					"Mismatch: edges(%d,%d)=%d stored count=%d" % (i,j,count, stored_count)
				else:
					assert count == 0, \
					"Mismatch: edges(%d,%d)=%d but no count stored" % (i,j, count)

def modularity(g, comms):
    """Comput modularity: Community-centric version."""
    Q = 0.0
    M = float(g.number_of_edges())
    for c, nodes in comms.iteritems():
        E_in = len(networkx.edge_boundary(g, nodes, nodes))/2
        assert E_in/2 == E_in//2
        K_in = sum(g.degree(n) for n in nodes)
        Q += E_in/(M*1) - (K_in/(2*M))**2
    return Q

    def test_path_edge_boundary(self):
        """Check edge boundaries in path nxgraph."""
        P10=self.P10

        assert_equal(nx.edge_boundary(P10,[]),[])
        assert_equal(nx.edge_boundary(P10,[],[]),[])
        assert_equal(nx.edge_boundary(P10,[1,2,3]),[(3, 4)])
        assert_equal(sorted(nx.edge_boundary(P10,[4,5,6])),[(4, 3), (6, 7)])
        assert_equal(sorted(nx.edge_boundary(P10,[3,4,5,6,7])),[(3, 2), (7, 8)])
        assert_equal(nx.edge_boundary(P10,[8,9,10]),[(8, 7)])
        assert_equal(sorted(nx.edge_boundary(P10,[4,5,6],[9,10])),[])
        assert_equal(nx.edge_boundary(P10,[1,2,3],[3,4,5]) ,[(2, 3), (3, 4)])

 def test_path_graph(self):
     P10 = cnlti(nx.path_graph(10), first_label=1)
     assert_equal(list(nx.edge_boundary(P10, [])), [])
     assert_equal(list(nx.edge_boundary(P10, [], [])), [])
     assert_equal(list(nx.edge_boundary(P10, [1, 2, 3])), [(3, 4)])
     assert_equal(sorted(nx.edge_boundary(P10, [4, 5, 6])),
                  [(4, 3), (6, 7)])
     assert_equal(sorted(nx.edge_boundary(P10, [3, 4, 5, 6, 7])),
                  [(3, 2), (7, 8)])
     assert_equal(list(nx.edge_boundary(P10, [8, 9, 10])), [(8, 7)])
     assert_equal(sorted(nx.edge_boundary(P10, [4, 5, 6], [9, 10])), [])
     assert_equal(list(nx.edge_boundary(P10, [1, 2, 3], [3, 4, 5])),
                  [(2, 3), (3, 4)])

    def _pre_init(self, pa_name, group, dgraph, fd, boundary_params):
        """Return a tuple of the form (pa_inputs, pa_outputs, renames)
        for the PseudoAssembly that would be created given the nodes in
        group and the given graph.
        """

        # First, find our group boundary
        self._orig_group_nodes = list(group) + list(boundary_params)
        allnodes = dgraph.find_prefixed_nodes(self._orig_group_nodes)
        out_edges = nx.edge_boundary(dgraph, allnodes)
        in_edges = nx.edge_boundary(dgraph,
                                    set(dgraph.nodes()).difference(allnodes))
        solver_states = []
        if fd is False:
            for comp in group:

                # Keep any node marked 'solver_state'. Note, only inputs can
                # be solver_states.
                solver_states.extend([node for node in dgraph.find_prefixed_nodes([comp])
                                      if 'solver_state' in dgraph.node[node]])

        pa_inputs = edges_to_dict(in_edges).values()
        pa_inputs.extend(solver_states)
        pa_outputs = set([a[0] for a in out_edges])

        renames = {}

        # Add pseudoassy inputs
        for varpath in list(flatten_list_of_iters(pa_inputs)) + \
                       list(pa_outputs):
            varname = to_PA_var(varpath, pa_name)
            if varpath in dgraph:
                renames[varpath] = varname
                old = dgraph.base_var(varpath)
                if old != varpath:
                    renames[old] = to_PA_var(old, pa_name)

        # make boundary params outputs of the PA
        pa_outputs.update(boundary_params)

        return pa_inputs, pa_outputs, renames

 def test_complete_graph(self):
     K10 = cnlti(nx.complete_graph(10), first_label=1)
     ilen = lambda iterable: sum(1 for i in iterable)
     assert_equal(list(nx.edge_boundary(K10, [])), [])
     assert_equal(list(nx.edge_boundary(K10, [], [])), [])
     assert_equal(ilen(nx.edge_boundary(K10, [1, 2, 3])), 21)
     assert_equal(ilen(nx.edge_boundary(K10, [4, 5, 6, 7])), 24)
     assert_equal(ilen(nx.edge_boundary(K10, [3, 4, 5, 6, 7])), 25)
     assert_equal(ilen(nx.edge_boundary(K10, [8, 9, 10])), 21)
     assert_equal(sorted(nx.edge_boundary(K10, [4, 5, 6], [9, 10])),
                  [(4, 9), (4, 10), (5, 9), (5, 10), (6, 9), (6, 10)])
     assert_equal(sorted(nx.edge_boundary(K10, [1, 2, 3], [3, 4, 5])),
                  [(1, 3), (1, 4), (1, 5), (2, 3), (2, 4),
                   (2, 5), (3, 4), (3, 5)])

def boundary_nodes(graph, nodes):
    # TODO: move to utils
#TODO: use networkx boundary nodes directly: does the same thing
    """ returns nodes at boundary of G based on edge_boundary from networkx """
    graph = unwrap_graph(graph)
    nodes = list(nodes)
    nbunch = list(unwrap_nodes(nodes))
    # find boundary
    b_edges = nx.edge_boundary(graph, nbunch)  # boundary edges
    internal_nodes = [s for (s, t) in b_edges]
    assert(all(n in nbunch for n in internal_nodes))  # check internal

    return wrap_nodes(graph, internal_nodes)

    def test_k10_edge_boundary(self):
        """Check edge boundaries in K10"""
        K10=self.K10

        assert_equal(nx.edge_boundary(K10,[]),[])
        assert_equal(nx.edge_boundary(K10,[],[]),[])
        assert_equal(len(nx.edge_boundary(K10,[1,2,3])),21)
        assert_equal(len(nx.edge_boundary(K10,[4,5,6,7])),24)
        assert_equal(len(nx.edge_boundary(K10,[3,4,5,6,7])),25)
        assert_equal(len(nx.edge_boundary(K10,[8,9,10])),21)
        assert_equal(sorted(nx.edge_boundary(K10,[4,5,6],[9,10])),
                     [(4, 9), (4, 10), (5, 9), (5, 10), (6, 9), (6, 10)])
        assert_equal(nx.edge_boundary(K10,[1,2,3],[3,4,5]),
                     [(1, 3), (1, 4), (1, 5), (2, 3), (2, 4), 
                      (2, 5), (3, 4), (3, 5)])

	def sample_graph(self, enforce_min_degree=False, enforce_connected_comps=False):
		g2 = networkx.Graph()

		for i in self.state.partition.keys():
			edges_within = self.state.edge_count(i)
			nodes_i = self.state.get_nodes(i)
			g2.add_nodes_from(nodes_i)
			assert edges_within <= len(nodes_i) * (len(nodes_i)-1)/2 and edges_within >= 0
			while edges_within > 0:
				u = random.choice(nodes_i)
				v = random.choice(nodes_i)
				if u == v or g2.has_edge(u,v):
					continue
				else:
					g2.add_edge(u,v)
					edges_within -= 1

			for j in self.state.neighbors(i):
				if i <= j:
					continue
				edges_between = self.state.edge_count(i,j)
				nodes_j = self.state.get_nodes(j)
				assert edges_between <= len(nodes_i) * len(nodes_j) and edges_between >= 0
				while edges_between > 0:
					u = random.choice(nodes_i)
					v = random.choice(nodes_j)
					if g2.has_edge(u,v):
						continue
					else:
						g2.add_edge(u,v)
						edges_between -= 1

		if enforce_min_degree:
			self.min_degree(g2)
		if enforce_connected_comps:
			self.conn_comps(g2)

		if debug_sampling:
			#import pdb
			print "Checking sampled graph"
			for i in self.state.partition:
				nodes_i = self.state.get_nodes(i)
				assert self.state.edge_count(i) == len(networkx.subgraph(g2, nodes_i).edges())
				for j in self.state.partition:
					if i <= j:
						continue
					nodes_j = self.state.get_nodes(j)
					assert self.state.edge_count(i,j) == len(networkx.edge_boundary(g2, nodes_i, nodes_j))
		return g2

def boundary_nodes(G, nodes):
    #TODO: move to utils
    """ returns nodes at boundary of G
    TODO: check works for both directed and undirected graphs
    based on edge_boundary from networkx """
    import autonetkit.ank as ank_utils
    graph = ank_utils.unwrap_graph(G)
    nodes = list(nodes)
    nbunch = list(ank_utils.unwrap_nodes(nodes))
    # find boundary
    b_edges = nx.edge_boundary(graph, nbunch) # boundary edges
    internal_nodes = [s for (s, t) in b_edges]
    assert(all(n in nbunch for n in internal_nodes)) # check internal 

    return ank_utils.wrap_nodes(G, internal_nodes)

    def best_edge(component):
        """Returns the optimum (minimum or maximum) edge on the edge
        boundary of the given set of nodes.

        A return value of ``None`` indicates an empty boundary.

        """
        # TODO In Python 3.4 and later, we can just do
        #
        #     boundary = nx.edge_boundary(G, component, data=weight)
        #     return opt(boundary, key=lambda e: e[-1][weight], default=None)
        #
        # which is better because it doesn't require creating a list.
        boundary = list(nx.edge_boundary(G, component, data=True))
        if not boundary:
            return None
        return opt(boundary, key=lambda e: e[-1][weight])

def deltaQ(g, c1, c2):
    """Change of modularity if c1 and c2 were merged"""
    # The formula for modulary change in the paper is deceiving.
    # Derive it yourself:
    #  deltaQ = E12/M - 2*K1*K2/(2M)^2
    M = float(g.number_of_edges())
    E12 = len(networkx.edge_boundary(g, comms[c1], comms[c2]))
    #E1 = g.subgraph(comms[c1]).number_of_edges()
    #E2 = g.subgraph(comms[c2]).number_of_edges()
    K1 = sum(g.degree(n) for n in comms[c1])
    K2 = sum(g.degree(n) for n in comms[c2])
    #dQ = 2*(eij - ai*aj)
    dQ = E12/M - 2*K1*K2/(2*M)**2
    return dQ

    # Newman method:
    e12 = E12/M
    a1 = (K1 - E1)/M
    a2 = (K2 - E2)/M
    return e12 - a1*a2