from collections import deque
from heapq import heappush, heappop

demand, capacity, weight = 'demand', 'capacity', 'weight'
def _dijkstra(G, source, get_weight, pred=None, paths=None, cutoff=None, target=None):
    G_succ = G.succ if G.is_directed() else G.adj

    push = heappush
    pop = heappop
    dist = {}  # dictionary of final distances
    seen = {source: 0}
    c = count()
    fringe = []  # use heapq with (distance,label) tuples
    push(fringe, (0, next(c), source))
    while fringe:
        (d, _, v) = pop(fringe)
        if v in dist:
            continue  # already searched this node.
        dist[v] = d
        if v == target:
            break

        for u, e in G_succ[v].items():
            cost = get_weight(v, u, e)
            if cost is None:
                continue
            vu_dist = dist[v] + get_weight(v, u, e)
            if cutoff is not None:
                if vu_dist > cutoff:
                    continue
            if u in dist:
                if vu_dist < dist[u]:
                    raise ValueError('Contradictory paths found:',
                                     'negative weights?')
            elif u not in seen or vu_dist < seen[u]:
                seen[u] = vu_dist
                push(fringe, (vu_dist, next(c), u))
                if paths is not None:
                    paths[u] = paths[v] + [u]
                if pred is not None:
                    pred[u] = [v]
            elif vu_dist == seen[u]:
                if pred is not None:
                    pred[u].append(v)

    if paths is not None:
        return (dist, paths)
    if pred is not None:
        return (pred, dist)
    return dist

    
class Graph:
    def connected_components(G):
        seen={}
        for v in G:
            if v not in seen:
                c = sp_length(G, v)
                yield list(c)
                seen.update(c)
            
    def bidirectional_dijkstra(G, source, target, weight='weight'):
        if source == target:
            return (0, [source])
        push = heappush
        pop = heappop
        # Init:   Forward             Backward
        dists  = [{},                {}]  # dictionary of final distances
        paths  = [{source: [source]}, {target: [target]}]  # dictionary of paths
        fringe = [[],                []]  # heap of (distance, node) tuples for
                                          # extracting next node to expand
        seen   = [{source: 0},        {target: 0}]  # dictionary of distances to
                                                    # nodes seen
        c = count()
        # initialize fringe heap
        push(fringe[0], (0, next(c), source))
        push(fringe[1], (0, next(c), target))
        # neighs for extracting correct neighbor information
        if G.is_directed():
            neighs = [G.successors_iter, G.predecessors_iter]
        else:
            neighs = [G.neighbors_iter, G.neighbors_iter]
        # variables to hold shortest discovered path
        #finaldist = 1e30000
        finalpath = []
        dir = 1
        while fringe[0] and fringe[1]:
            # choose direction
            # dir == 0 is forward direction and dir == 1 is back
            dir = 1 - dir
            # extract closest to expand
            (dist, _, v) = pop(fringe[dir])
            if v in dists[dir]:
                # Shortest path to v has already been found
                continue
            # update distance
            dists[dir][v] = dist  # equal to seen[dir][v]
            if v in dists[1 - dir]:
                # if we have scanned v in both directions we are done
                # we have now discovered the shortest path
                return (finaldist, finalpath)
    
            for w in neighs[dir](v):
                if(dir == 0):  # forward
                    if G.is_multigraph():
                        minweight = min((dd.get(weight, 1)
                                         for k, dd in G[v][w].items()))
                    else:
                        minweight = G[v][w].get(weight, 1)
                    vwLength = dists[dir][v] + minweight  # G[v][w].get(weight,1)
                else:  # back, must remember to change v,w->w,v
                    if G.is_multigraph():
                        minweight = min((dd.get(weight, 1)
                                         for k, dd in G[w][v].items()))
                    else:
                        minweight = G[w][v].get(weight, 1)
                    vwLength = dists[dir][v] + minweight  # G[w][v].get(weight,1)
    
                if w in dists[dir]:
                    if vwLength < dists[dir][w]:
                        raise ValueError(
                            "Contradictory paths found: negative weights?")
                elif w not in seen[dir] or vwLength < seen[dir][w]:
                    # relaxing
                    seen[dir][w] = vwLength
                    push(fringe[dir], (vwLength, next(c), w))
                    paths[dir][w] = paths[dir][v] + [w]
                    if w in seen[0] and w in seen[1]:
                        # see if this path is better than than the already
                        # discovered shortest path
                        totaldist = seen[0][w] + seen[1][w]
                        if finalpath == [] or finaldist > totaldist:
                            finaldist = totaldist
                            revpath = paths[1][w][:]
                            revpath.reverse()
                            finalpath = paths[0][w] + revpath[1:]
        raise NetworkXNoPath("No path between %s and %s." % (source, target))

class DiGraph:
    def __init__(self, edges):
        self.edges = edges
        self.nodes = None
        pass

    def nodes_iter(self):
        return iter(self.nodes)
    def edges(self):
        return self.edges

    def add_nodes_from(self, nodes):
        self.nodes
        
def _initial_tree_solution(G):
    
    H = DiGraph((edge for edge in G.edges() if edge[2].get(capacity, 1) > 0))
    demand_nodes = (node for node in G.nodes_iter() if node[1].get(demand, 0) != 0)
    H.add_nodes_from(demand_nodes)
    r = H.nodes()[0]

    T = DiGraph()
    y = {r: 0}
    artificial_edges = []
    flow_cost = 0

    n = H.number_of_nodes()
    try:
        max_weight = max(abs(d[weight]) for u, v, d in H.edges() if weight in d)
    except ValueError:
        max_weight = 0
    huge_weight = 1 + n * max_weight

    for v, d in H.nodes()[1:]:
        v_demand = d.get(demand, 0)
        if v_demand >= 0:
            if not (r, v) in H.edges():
                H.add_edge(r, v, {weight: huge_weight, 'flow': v_demand})
                artificial_edges.append((r, v))
                y[v] = H[r][v].get(weight, 0)
                T.add_edge(r, v)
                flow_cost += v_demand * H[r][v].get(weight, 0)

            else: # (r, v) in H.edges()
                if (not capacity in H[r][v]
                    or v_demand <= H[r][v][capacity]):
                    H[r][v]['flow'] = v_demand
                    y[v] = H[r][v].get(weight, 0)
                    T.add_edge(r, v)
                    flow_cost += v_demand * H[r][v].get(weight, 0)

                else: # existing edge does not have enough capacity
                    new_label = generate_unique_node()
                    H.add_edge(r, new_label, {weight: huge_weight, 'flow': v_demand})
                    H.add_edge(new_label, v, {weight: huge_weight, 'flow': v_demand})
                    artificial_edges.append((r, new_label))
                    artificial_edges.append((new_label, v))
                    y[v] = 2 * huge_weight
                    y[new_label] = huge_weight
                    T.add_edge(r, new_label)
                    T.add_edge(new_label, v)
                    flow_cost += 2 * v_demand * huge_weight

        else:
            if not (v, r) in H.edges():
                H.add_edge(v, r, {weight: huge_weight, 'flow': -v_demand})
                artificial_edges.append((v, r))
                y[v] = -H[v][r].get(weight, 0)
                T.add_edge(v, r)
                flow_cost += -v_demand * H[v][r].get(weight, 0)

            else:
                if (not capacity in H[v][r] or -v_demand <= H[v][r][capacity]):
                    H[v][r]['flow'] = -v_demand
                    y[v] = -H[v][r].get(weight, 0)
                    T.add_edge(v, r)
                    flow_cost += -v_demand * H[v][r].get(weight, 0)
                else: 
                    new_label = generate_unique_node()
                    H.add_edge(v, new_label, {weight: huge_weight, 'flow': -v_demand})
                    H.add_edge(new_label, r, {weight: huge_weight, 'flow': -v_demand})
                    artificial_edges.append((v, new_label))
                    artificial_edges.append((new_label, r))
                    y[v] = -2 * huge_weight
                    y[new_label] = -huge_weight
                    T.add_edge(v, new_label)
                    T.add_edge(new_label, r)
                    flow_cost += 2 * -v_demand * huge_weight

    return H, T, y, artificial_edges, flow_cost, r


def _find_entering_edge(H, c):
    capacity = 'capacity'
    new_edge = ()
    for u, v, d in H.edges_iter():
        if d.get('flow', 0) == 0:
            if c[(u, v)] < 0:
                new_edge = (u, v)
                break
        else:
            if capacity in d:
                if (d.get('flow', 0) == d[capacity]
                    and c[(u, v)] > 0):
                    new_edge = (u, v)
                    break
    return new_edge


def _find_leaving_edge(H, T, cycle, new_edge, capacity = 'capacity', reverse=False):
    eps = False
    leaving_edge = ()

    # If cycle is a digon.
    if len(cycle) == 3:
        u, v = new_edge
        if reverse:
            if H[u][v].get('flow', 0) > H[v][u].get('flow', 0):
                return (v, u), H[v][u].get('flow', 0)
            else:
                return (u, v), H[u][v].get('flow', 0)
        else:
            uv_residual = H[u][v].get(capacity, 0) - H[u][v].get('flow', 0)
            vu_residual = H[v][u].get(capacity, 0) - H[v][u].get('flow', 0)
            if (uv_residual > vu_residual):
                return (v, u), vu_residual
            else:
                return (u, v), uv_residual

    # Find the forward edge with the minimum value for capacity - 'flow'
    # and the reverse edge with the minimum value for 'flow'.
    for index, u in enumerate(cycle[:-1]):
        edge_capacity = False
        edge = ()
        v = cycle[index + 1]
        if (u, v) in T.edges() + [new_edge]: #forward edge
            if capacity in H[u][v]: # edge (u, v) has finite capacity
                edge_capacity = H[u][v][capacity] - H[u][v].get('flow', 0)
                edge = (u, v)
        else: #reverse edge
            edge_capacity = H[v][u].get('flow', 0)
            edge = (v, u)

        # Determine if edge might be the leaving edge.
        if edge:
            if leaving_edge:
                if edge_capacity < eps:
                    eps = edge_capacity
                    leaving_edge = edge
            else:
                eps = edge_capacity
                leaving_edge = edge

    if not leaving_edge:
        raise ValueError("No leaving edge")

    return leaving_edge, eps


def _create_flow_dict(G, H):
    flow_dict = dict([(u, {}) for u in G])

    for u in G.nodes_iter():
        for v in G.neighbors(u):
            if H.has_edge(u, v):
                flow_dict[u][v] = H[u][v].get('flow', 0)
            else:
                flow_dict[u][v] = 0
    return flow_dict


def network_simplex(G):
    demand, capacity,  weight = 'demand', 'capacity', 'weight'
    # Fix an arbitrarily chosen root node and find an initial tree solution.
    H, T, y, artificial_edges, flow_cost, r = \
            _initial_tree_solution(G, demand = demand, capacity = capacity, weight = weight)

    # Initialize the reduced costs.
    c = {}
    for u, v, d in H.edges_iter():
        c[(u, v)] = d.get(weight, 0) + y[u] - y[v]

    # Main loop.
    while True:
        new_edge = _find_entering_edge(H, c, capacity = capacity)
        if not new_edge:
            break # Optimal basis found. Main loop is over.
        cycle_cost = abs(c[new_edge])

        # Find the cycle created by adding new_edge to T.
        path1 = shortest_path(T.to_undirected(), r, new_edge[0])
        path2 = shortest_path(T.to_undirected(), r, new_edge[1])
        join = r
        for index, node in enumerate(path1[1:]):
            if index + 1 < len(path2) and node == path2[index + 1]:
                join = node
            else:
                break
        path1 = path1[path1.index(join):]
        path2 = path2[path2.index(join):]
        cycle = []
        if H[new_edge[0]][new_edge[1]].get('flow', 0) == 0:
            reverse = False
            path2.reverse()
            cycle = path1 + path2
        else: # new_edge is at capacity
            reverse = True
            path1.reverse()
            cycle = path2 + path1

        # Find the leaving edge. Will stop here if cycle is an infinite
        # capacity negative cost cycle.
        leaving_edge, eps = _find_leaving_edge(H, T, cycle, new_edge, capacity=capacity, reverse=reverse)

        # Actual augmentation happens here. If eps = 0, don't bother.
        if eps:
            flow_cost -= cycle_cost * eps
            if len(cycle) == 3:
                if reverse:
                    eps = -eps
                u, v = new_edge
                H[u][v]['flow'] = H[u][v].get('flow', 0) + eps
                H[v][u]['flow'] = H[v][u].get('flow', 0) + eps
            else:
                for index, u in enumerate(cycle[:-1]):
                    v = cycle[index + 1]
                    if (u, v) in T.edges() + [new_edge]:
                        H[u][v]['flow'] = H[u][v].get('flow', 0) + eps
                    else: # (v, u) in T.edges():
                        H[v][u]['flow'] -= eps

        # Update tree solution.
        T.add_edge(*new_edge)
        T.remove_edge(*leaving_edge)

        # Update distances and reduced costs.
        if new_edge != leaving_edge:
            forest = DiGraph(T)
            forest.remove_edge(*new_edge)
            R, notR = connected_component_subgraphs(forest.to_undirected())
            if r in notR.nodes(): # make sure r is in R
                R, notR = notR, R
            if new_edge[0] in R.nodes():
                for v in notR.nodes():
                    y[v] += c[new_edge]
            else:
                for v in notR.nodes():
                    y[v] -= c[new_edge]
            for u, v in H.edges():
                if u in notR.nodes() or v in notR.nodes():
                    c[(u, v)] = H[u][v].get(weight, 0) + y[u] - y[v]


    for u, v in artificial_edges:
        if H[u][v]['flow'] != 0:
            raise ValueError("No flow satisfying all demands.")
        H.remove_edge(u, v)

    for u in H.nodes():
        if not u in G:
            H.remove_node(u)

    flow_dict = _create_flow_dict(G, H)

    return flow_cost, flow_dict