scp.py

import networkx as nx
import random


# Note NetworkX, flow must have positive capacity
# Hence we must build edges in both directions
# This is undesirable but we can't do much about it
# We assume in the graph, if (a,b) is an edge, then (b,a) is also an edge
# We also assume the default orientation is (a,b) where a<b

from itertools import chain, combinations, product

def powerset(iterable):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return list(map(frozenset,chain.from_iterable(combinations(s, r) for r in range(len(s)+1))))

def normalize(i,j):
    return min(i,j), max(i,j)

# Input: Undirected graph graph with weights, set of requests
def SCP(G, A, outside_edge):
    must_use = 0
    apsp = dict(nx.all_pairs_dijkstra_path_length(G))
    for (i,j) in A:
        must_use+=apsp[i][j]

    # Find the spanning tree
    T = nx.Graph(G)
    for x,y in outside_edge:
        T.remove_edge(x,y)
    
    # print(nx.is_tree(T))
    #print("T", T.edges())
    
    F = []
    # find the edges index the fundamental cycle basis
    for (x,y) in G.edges():
        i,j = normalize(x,y)
        if not T.has_edge(i,j):
            F.append((i,j))

    #print("F", F)

    # find the fundamental cycle basis
    C = dict()
    for x in range(len(F)):
        (i,j) = F[x]
        path = nx.shortest_path(T, j, i)
        path_edges = list(zip(path, path[1:]))
        cycle = [(i,j)]+path_edges
        C[x] = cycle
    #print("C", C)

    # find the minimum cost flow satisfies all requests
    # find the flow in the edge space of G
    # note we would need an standard orientation, which would be (i,j) where i<j
    f = min_circulation(G, A)
    #print("flow",f)

    p = len(F)

    if p==0:
        return

    lamb = list([0]*p)
    for i in range(len(F)):
        lamb[i] = f[F[i]]
    # print(lamb)
    
    # try all possible result, the Delta of the cycle basis coefficient from -p to p
    # for each coordinate
    # basically we want to know the optimal of (f + \sum_{i=1}^p Delta_1C_1)   

    results = []
    for Delta in map(list,product(range(-p,p+1), repeat=p)):
        circulation = circulation_from_fundamental_basis(Delta, C)
        new_flow = flow_sum(circulation, f)
        results.append((opt_given_flow(new_flow, A, G) + must_use,Delta))

    close_results = [(x,y) for x,y in results if abs(max(y))<=1 and abs(min(y))<=1]

    if min(close_results)[0] != min(results)[0]:
        print("FAIL")
        print(min(results),print(G.edges(data=True)),print(A))
    #print(results)
    return min(results)

# compute sum of two flow
def flow_sum(h, g):
    f = dict()
    for x in h.keys():
        f[x] = 0
    for x in g.keys():
        f[x] = 0
    for x in h.keys():
        f[x] += h[x]
    for x in g.keys():
        f[x] += g[x]
    return f

# given fundamental basis and coefficients
# compute \sum_{i} Delta_i*C_i
def circulation_from_fundamental_basis(Delta, C):
    f = dict()
    for i in range(len(Delta)):
        for (x,y) in C[i]:
            if x<y:
                if (x,y) not in f:
                    f[(x,y)]=0
                f[(x,y)]+=Delta[i]
            else:
                if (y,x) not in f:
                    f[(y,x)]=0
                f[(y,x)]-=Delta[i]
    return f

# Input: given flow, demand pairs and graph, compute a connected flow
def opt_given_flow(f, A, G):
    # compute the cost of the flow
    flow_cost = sum([abs(f[(i,j)])*G[i][j]['weight'] for (i,j) in f.keys()])

    # we build a new graph by adding edges in A
    # and also flow of non-zero value
    H = nx.Graph()
    H.add_nodes_from(G.nodes())
    H.add_edges_from(A)
    for (i,j) in f.keys():
        if f[(i,j)] != 0:
            H.add_edge(i,j)
    connected_components = list(nx.connected_components(H))

    if len(connected_components)==1:
        return flow_cost
    #print("connected components", connected_components)

    # contract all the components, resulting contracted graph is H
    Gprime = nx.MultiGraph(G)
    # print("Gprime",Gprime.edges(data=True))
    Hprime = nx.quotient_graph(Gprime, connected_components, weight='weight')
    # print("Hprime",Hprime.edges(data=True))

    # build a simple graph out of H
    H = nx.Graph()
    for u,v,data in Hprime.edges(data=True):
        w = data[0]['weight']
        if H.has_edge(u,v):
            H[u][v]['weight'] = min(H[u][v]['weight'],w)
        else:
            H.add_edge(u, v, weight=w)

    # print(H.edges(data=True))

    # terminals are contracted vertices
    terminals = [x for x in H.nodes() if len(x)>1]
    steiner_tree_cost = 2 * steiner_tree(H, terminals)
    return steiner_tree_cost + flow_cost

# Input: undirected graph, terminal vertices, output the weight
def steiner_tree(G, T):
    weight = 999999999
    # remove degree 1 Steiner vertices
    S = frozenset(G.nodes()) - frozenset(T)

    # remove all degree 1 steiner vertex
    while len([x for x in S&frozenset(G.nodes()) if G.degree[x]<=1])!=0:
        for x in S&frozenset(G.nodes()):
            if G.degree[x]<=1:
                G.remove_node(x)

    # remove all degree 2 steiner vertex
    while len([x for x in S&frozenset(G.nodes()) if G.degree[x]==2])!=0:
        for x in S&frozenset(G.nodes()):
            if G.degree[x]==2:
                new_weight = sum([t['weight'] for (_,_,t) in G.edges(x,data=True)])
                [u,v] = list(G[x])
                G.add_edge(u,v,weight=new_weight)
                G.remove_node(x)

    S = S&frozenset(G.nodes())

    # now S only contain Stiner vertex of degree at least 3
    # then do 2^k times MST, k = len(S)
    for X in powerset(S):
        H = nx.induced_subgraph(G, frozenset(G.nodes())-X)
        # first check connectivity
        if nx.is_connected(H):
            T = nx.minimum_spanning_tree(H)
            # obtain tree weight
            weight = min(weight,sum([data['weight'] for a,b,data in T.edges(data=True)]))
    return weight

# given a graph and some requests, compute the min circulation satisfies the requests
def min_circulation(G, A):
    G = G.to_directed()
    demand = dict()
    for v in list(G.nodes()):
        G.nodes[v]["demand"] = 0 
    for (i,j) in A:
        G.nodes[i]["demand"] -= 1
        G.nodes[j]["demand"] += 1

    flowDict = nx.min_cost_flow(G, demand='demand', capacity='capacity', weight='weight')

    f = dict()
    for (x,y) in G.edges():
        i = min(x,y)
        j = max(x,y)
        f[(i,j)] = flowDict[i][j] - flowDict[j][i]

    return f



### 3 item graph

G = nx.Graph()

basepath = [(0,1),(1,2),(2,3),(3,4),(4,5)]
path1 = [(1,6),(6,7),(7,8),(8,9),(9,10),(10,11),(11,3)]
path2 = [(2,12),(12,13),(13,14),(14,15),(15,16),(16,17),(17,4)]
path3 = [(0,18),(18,19),(19,20),(20,21),(21,22),(22,23),(23,5)]

G.add_edges_from(basepath)
G.add_edges_from(path1)
G.add_edges_from(path2)
G.add_edges_from(path3)

outside_edge = [(1,6),(2,12),(0,18)]

#G = nx.Graph()
#G.add_edge(0,1,weight=11)
#G.add_edge(1,2,weight=4)
#G.add_edge(2,3,weight=11)
#G.add_edge(3,5,weight=11)
#G.add_edge(4,5,weight=3)
#G.add_edge(0,4,weight=15)
#G.add_edge(0,3,weight=16)

#A = [(4,5),(1,2)]

#SCP(G,A)

# G = nx.Graph()
# G.add_edge(0,1)
# G.add_edge(1,2)
# G.add_edge(2,3)
# G.add_edge(3,0)
# G.add_edge(3,4)
# G.add_edge(4,5)
# G.add_edge(5,6)
# G.add_edge(6,1)
# G.add_edge(6,7)
# G.add_edge(7,8)
# G.add_edge(8,9)
# G.add_edge(9,4)


for iter in range(100000):
    n = 24
    #G = nx.random_unlabeled_tree(n)
    # print(G.edges())


    # randomly add 3 edges
    #for i in range(3):
    #    a = random.randrange(0,n)
    #    b = random.randrange(0,n)
    #    G.add_edge(a,b)

    # random weights
    for x,y in G.edges():
        G[x][y]['weight'] = int(random.random()*100000)

    for iter in range(100):
        A = []
        # random requests
        for i in range(10):
            a = random.randrange(0,n)
            b = random.randrange(0,n)
            A.append((a,b))
        SCP(G,A,outside_edge)