Bipartite Matching¶
Bipartite Matching問題を解くアルゴリズムを実装してみる。Figure 8のような二部グラフをまずは作成する。
In [46]:
# Figure 8のグラフを作成
# もろもろをインポート
import networkx as nx
import matplotlib.pyplot as plt
from collections import deque
G = nx.DiGraph()
edgelist = [('v1','w1'),('v1','w2'),('v2','w1'),('v2','w2'),('v3','w2'),('v4','w3'),('v4','w4'),('v5','w5'),('v5','w6')]
G.add_edges_from(edgelist)
for e in G.edges():
G[e[0]][e[1]]['capacity'] = 1
G[e[0]][e[1]]['flow'] = 0
# 描画
pos={'v1':(0,10),'v2':(0,8),'v3':(0,6),'v4':(0,4),'v5':(0,2),'w1':(2,10),'w2':(2,8),'w3':(2,6),'w4':(2,4),'w5':(2,2),'w6':(2,0)}
edge_labels = {(i, j): (w['capacity'], str((w['flow'])) ) for i, j, w in G.edges(data=True)}
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_color='r')
nx.draw_networkx_labels(G, pos)
nx.draw(G, pos)
plt.axis('off')
plt.show()
In [47]:
# V,Wを設定。Gに加えてV,Wを与えれば二部グラフは定まる
V = ['v1','v2','v3','v4','v5']
W = ['w1','w2','w3','w4','w5','w6']
In [48]:
def add_st(G,V,W):
'''Just add a source s and a sink t to the original bipartite graph G.
Parameters
----------
G : bipartite graph
V,W : the lists of nodes, two groups of G
Return
------
H : modified graph
'''
H = G.copy()
H.add_nodes_from(['s','t'])
for p in V:
H.add_edge('s', p)
H['s'][p]['flow'] = 0
H['s'][p]['capacity'] = 1
for p in W:
H.add_edge(p, 't')
H[p]['t']['flow'] = 0
H[p]['t']['capacity'] = 1
return H
In [19]:
G0 = add_st(G,V,W)
In [20]:
# s,tを加えたグラフを描画
pos={'s':(-1,5),'t':(3,5),'v1':(0,10),'v2':(0,8),'v3':(0,6),'v4':(0,4),'v5':(0,2),'w1':(2,10),'w2':(2,8),'w3':(2,6),'w4':(2,4),'w5':(2,2),'w6':(2,0)}
edge_labels = {(i, j): (w['capacity'], str((w['flow'])) ) for i, j, w in G0.edges(data=True)}
nx.draw_networkx_edge_labels(G0, pos, edge_labels=edge_labels, font_color='r')
nx.draw_networkx_labels(G0, pos)
nx.draw(G0, pos)
plt.axis('off')
plt.show()
まず、push-relabel algorithmを用いてbipartite matchingを解いてみる。
In [49]:
# ch03で作成したものたち
import numpy as np
def makeResidualGraph(G):
'''
Input: a graph G
Output: its residual graph Gf
'''
Gf = G.copy()
edgeList = G.edges()
for edge in edgeList:
# Initialize flow
Gf[edge[0]][edge[1]]['flow'] = 0
# 逆向きのedgeがないものは追加
if not (edge[1], edge[0]) in edgeList:
Gf.add_edge(edge[1],edge[0])
Gf[edge[1]][edge[0]]['capacity'] = Gf[edge[0]][edge[1]]['flow']
Gf[edge[1]][edge[0]]['flow'] = 0
return Gf
def Initialize(G, s, t):
'''
Inputs:
G: a residual graph
V: verteces
height
excess
E: edges
flow(preflow)
capacity
NB: flow <= capacity
s: a start point
t: an end point
Output:
None
Initialize the height and excess of each vertex.
'''
for v in G.nodes():
G.node[v]['excess'] = 0
for v in G.nodes():
# 始点以外の点について
if v != s:
G.node[v]['height'] = 0
# G.node[v]['excess'] = 0
for p in G.neighbors(v):
G[v][p]['flow'] = 0
# 始点について
else:
# 始点の高さをnに
G.node[v]['height'] = G.number_of_nodes() # n
# 始点から出る枝のpreflowをcapacityギリギリに
# このとき、始点に隣接する点のexcessを更新
for p in G.neighbors(v):
G[v][p]['flow'] = G[v][p]['capacity']
# backward edgeのcapacityの更新
G[p][v]['capacity'] = G[v][p]['flow']
G.node[p]['excess'] = G[v][p]['flow']
def Push(G, v, w, forwardEdges):
'''
G: a graph(a residual graph)
v, w: vertices of G such that h(v) = h(w) + 1
forwardEdges: a list of the edges of the original graph
'''
# 更新量diffを計算
residualCapacity = G[v][w]['capacity'] - G[v][w]['flow']
diff = np.min([G.node[v]['excess'], residualCapacity])
''' for debugging
print(v)
print(w)
print('diff = ' + str(diff))
print('res.cap. = ' + str(residualCapacity))
print('excess of ' + v + ': ' + str(G.node[v]['excess']))
print('excess of ' + w + ': ' + str(G.node[w]['excess']))
'''
# (v,w), (w,v)を更新, p.3を参考に
if (v,w) in forwardEdges:
G[v][w]['flow'] += diff
G[w][v]['capacity'] += diff
G.node[v]['excess'] -= diff
G.node[w]['excess'] += diff
else:
# このとき、(w,v)がforward edge、(v,w)がbackward edge
G[w][v]['flow'] -= diff
G[v][w]['capacity'] -= diff
G.node[w]['excess'] += diff
G.node[v]['excess'] -= diff
def loopCondition(G, s, t):
nodelist = G.nodes()
nodelist.remove(s)
nodelist.remove(t)
for v in nodelist:
if(G.node[v]['excess'] > 0):
return True
return False
def PushRelabel(G, s, t):
'''
Inputs:
G: a graph
s: a starting point
t: an end point
Output:
the graph G with its maximum flow
'''
# Forward Edges を記録
forwardEdges = G.edges()
# s,tを除いたnodeのlistを作る
nodeList = G.nodes()
nodeList.remove(s)
nodeList.remove(t)
# residual graph の作成
Gf = makeResidualGraph(G)
# Initialization
Initialize(Gf, s, t)
# Main Loop
while(loopCondition(Gf, s, t)):
# 高さが最も高く、かつexcess > 0の点vを探す
height = -100 # 適当に負の値を設定
for p in nodeList:
if(Gf.node[p]['excess'] > 0 and Gf.node[p]['height'] > height):
v = p
height = Gf.node[p]['height']
# h(v) = h(w) + 1 を満たす点wを探す
# vのneighborsの中から探す(そうじゃないとkey eroorに)
w = None
for p in Gf.neighbors(v):
if(Gf.node[v]['height'] == Gf.node[p]['height'] + 1 and (Gf[v][p]['capacity'] - Gf[v][p]['flow']) > 0):
w = p
break
# そのような点wがない場合、increment h(v)
if w == None:
Gf.node[v]['height'] += 1
# ある場合
else:
Push(Gf, v, w, forwardEdges)
# もともと無かったedgeを消去
for edge in Gf.edges():
if not edge in forwardEdges:
Gf.remove_edge(edge[0],edge[1])
return Gf
In [50]:
def biMatch_pr(G,V,W):
'''Calculate the perfect matching of (G,V,W) by using push-relabel algorithm.
Parameters
----------
G : bipartite graph
V,W : the lists of nodes, two groups of G
Return
------
G : modified graph with perfect matching
'''
H = add_st(G,V,W)
I = PushRelabel(H, 's', 't')
I.remove_node('s')
I.remove_node('t')
return I
In [51]:
H = biMatch_pr(G,V,W)
v4 w3 diff = 1 res.cap. = 1 excess of v4: 1 excess of w3: 0 w3 t diff = 1 res.cap. = 1 excess of w3: 1 excess of t: 0 v5 w5 diff = 1 res.cap. = 1 excess of v5: 1 excess of w5: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w5 t diff = 1 res.cap. = 1 excess of w5: 1 excess of t: 1 w2 t diff = 1 res.cap. = 1 excess of w2: 1 excess of t: 2 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 t diff = 1 res.cap. = 1 excess of w1: 1 excess of t: 3 v1 w1 diff = 1 res.cap. = 1 excess of v1: 1 excess of w1: 0 w1 v2 diff = 1 res.cap. = 1 excess of w1: 1 excess of v2: 0 v2 w2 diff = 1 res.cap. = 1 excess of v2: 1 excess of w2: 0 w2 v3 diff = 1 res.cap. = 1 excess of w2: 1 excess of v3: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w2 v2 diff = 1 res.cap. = 1 excess of w2: 1 excess of v2: 0 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 v1 diff = 1 res.cap. = 1 excess of w1: 1 excess of v1: 0 v1 w1 diff = 1 res.cap. = 1 excess of v1: 1 excess of w1: 0 w1 v2 diff = 1 res.cap. = 1 excess of w1: 1 excess of v2: 0 v2 w2 diff = 1 res.cap. = 1 excess of v2: 1 excess of w2: 0 w2 v3 diff = 1 res.cap. = 1 excess of w2: 1 excess of v3: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w2 v2 diff = 1 res.cap. = 1 excess of w2: 1 excess of v2: 0 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 v1 diff = 1 res.cap. = 1 excess of w1: 1 excess of v1: 0 v1 w1 diff = 1 res.cap. = 1 excess of v1: 1 excess of w1: 0 w1 v2 diff = 1 res.cap. = 1 excess of w1: 1 excess of v2: 0 v2 w2 diff = 1 res.cap. = 1 excess of v2: 1 excess of w2: 0 w2 v3 diff = 1 res.cap. = 1 excess of w2: 1 excess of v3: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w2 v2 diff = 1 res.cap. = 1 excess of w2: 1 excess of v2: 0 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 v1 diff = 1 res.cap. = 1 excess of w1: 1 excess of v1: 0 v1 w1 diff = 1 res.cap. = 1 excess of v1: 1 excess of w1: 0 w1 v2 diff = 1 res.cap. = 1 excess of w1: 1 excess of v2: 0 v2 w2 diff = 1 res.cap. = 1 excess of v2: 1 excess of w2: 0 w2 v3 diff = 1 res.cap. = 1 excess of w2: 1 excess of v3: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w2 v2 diff = 1 res.cap. = 1 excess of w2: 1 excess of v2: 0 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 v1 diff = 1 res.cap. = 1 excess of w1: 1 excess of v1: 0 v1 w1 diff = 1 res.cap. = 1 excess of v1: 1 excess of w1: 0 w1 v2 diff = 1 res.cap. = 1 excess of w1: 1 excess of v2: 0 v2 w2 diff = 1 res.cap. = 1 excess of v2: 1 excess of w2: 0 w2 v3 diff = 1 res.cap. = 1 excess of w2: 1 excess of v3: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w2 v2 diff = 1 res.cap. = 1 excess of w2: 1 excess of v2: 0 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 v1 diff = 1 res.cap. = 1 excess of w1: 1 excess of v1: 0 v1 w1 diff = 1 res.cap. = 1 excess of v1: 1 excess of w1: 0 w1 v2 diff = 1 res.cap. = 1 excess of w1: 1 excess of v2: 0 v2 w2 diff = 1 res.cap. = 1 excess of v2: 1 excess of w2: 0 w2 v3 diff = 1 res.cap. = 1 excess of w2: 1 excess of v3: 0 v3 w2 diff = 1 res.cap. = 1 excess of v3: 1 excess of w2: 0 w2 v2 diff = 1 res.cap. = 1 excess of w2: 1 excess of v2: 0 v2 w1 diff = 1 res.cap. = 1 excess of v2: 1 excess of w1: 0 w1 v1 diff = 1 res.cap. = 1 excess of w1: 1 excess of v1: 0 v1 s diff = 1 res.cap. = 1 excess of v1: 1 excess of s: 0
In [52]:
# 描画
pos={'v1':(0,10),'v2':(0,8),'v3':(0,6),'v4':(0,4),'v5':(0,2),'w1':(2,10),'w2':(2,8),'w3':(2,6),'w4':(2,4),'w5':(2,2),'w6':(2,0)}
edge_labels = {(i, j): (w['capacity'], str((w['flow'])) ) for i, j, w in H.edges(data=True)}
nx.draw_networkx_edge_labels(H, pos, edge_labels=edge_labels, font_color='r')
nx.draw_networkx_labels(H, pos)
nx.draw(H, pos)
plt.axis('off')
plt.show()
Ford-Fulkersonでも求めてみる
In [39]:
#ch02で作ったやつ
import networkx as nx
import matplotlib.pyplot as plt
from collections import deque
import numpy as np
def bfsFlowPath(G, s, e):
'''
Search a path from s to t all of whose points have strictly positive capacity.
Inputs:
G: a graph
s: a start point
e: an end point
each edge has two attributes:
capacity: capacity
flow: its current flow which should be no more than its capacity
Output:
path: a list of edges which represents a path from s to e.
At each edge of the path, its current flow is strictly less than its capacity.
In case there is no path from s to t, return None.
'''
# 過去に訪れた点を記録
# sは最初から入れておく
past = [s]
# pathを記録するためのリスト
path = []
# 全ての点のsからの距離の初期値を無限大に
for p in G.nodes():
G.node[p]['dist'] = float('inf')
# node s の距離を0に
G.node[s]['dist'] = 0
# sに隣接する点をqueueに
# queueには、今後訪れるべき点が格納される
queue = deque()
for p in G.successors(s):
# current flow < capacity となるedgeだけをpathの候補に
# flow < capacity となるedge以外は存在しないものとして扱うのと同じ
if G[s][p]['flow'] < G[s][p]['capacity']:
queue.append(p)
# あとで条件分岐に用いる
numberOfSuccessorsOfSource = len(queue)
# sに隣接する点の距離を1に
for p in queue:
G.node[p]['dist'] = 1
# BFSを用いてflow < capacityを満たすsからeへのpathがあるのか調べる
# pastに過去に訪れた点を格納
while len(queue)>0:
v = queue.popleft()
if v == e: break
else:
past.append(v)
for p in G.successors(v):
# (過去に訪れていない and flow < capacity)を満たすedge
if ( (not p in past) and ( G[v][p]['flow'] < G[v][p]['capacity']) ):
if ( not p in queue):
queue.append(p)
if G.node[p]['dist'] > G.node[v]['dist'] + 1:
G.node[p]['dist'] = G.node[v]['dist'] + 1
# sからeへのpathが存在しない場合はNoneを返す
if numberOfSuccessorsOfSource == 0:
v = s
if v != e or v == s:
# print ('There is no path.')
return None
# 以下、sからeへのpathが存在する場合
# 終点から遡ってpathを形成する
pp = e
while (1):
if pp == s: break
pred = G.predecessors(pp)
count = 0
for p in pred:
# ここに、flow < capacity の条件を追加
if ( G.node[p]['dist'] == G.node[pp]['dist']-1 and G[p][pp]['flow'] < G[p][pp]['capacity']):
path.insert(0, (p,pp))
pp = p
break
else:
count += 1
# 条件を満たすedgeがない
if count == len(pred):
break
return path
def makeResidualGraph(G):
'''
Input: a graph G
Output: its residual graph Gf
'''
Gf = G.copy()
edgeList = G.edges()
for edge in edgeList:
# Initialize flow
Gf[edge[0]][edge[1]]['flow'] = 0
# 逆向きのedgeがないものは追加
if not (edge[1], edge[0]) in edgeList:
Gf.add_edge(edge[1],edge[0])
Gf[edge[1]][edge[0]]['capacity'] = Gf[edge[0]][edge[1]]['flow']
Gf[edge[1]][edge[0]]['flow'] = 0
return Gf
def fordFulkerson(G, s, t):
'''
Inputs:
G: a graph
s: a source vertex
t: a sink vertex
Outputs:
the graph G whose flow was modified by Ford-Fulkerson algorithm
In case there is no path from s to t, return None.
'''
# initialize flows
for e in G.edges():
G[e[0]][e[1]]['flow'] = 0
# Forward edgesを記録
forwardEdges = G.edges()
# Residual Graphの作成
Gf = makeResidualGraph(G)
# そもそもGにおいてsからtへのパスがあるのか確認
path = bfsFlowPath(G, s, t)
if path == None:
print ("There is no path from " + str(s) + " to "+ str(t) )
return None
# Gにおいてsからtへのパスがある場合
while(1):
# これ以上パスがみつからない場合は最適なのでそこでループを打ち切る
path = bfsFlowPath(Gf, s, t)
if path == None:
break
# まだパスがある（最適でない）場合
else:
# path上のedgeについて、capacity - flowの最小値を調べる
diff = float('inf')
for edge in path:
diff = np.min([diff, Gf[edge[0]][edge[1]]['capacity'] - Gf[edge[0]][edge[1]]['flow'] ])
# path上のedgeのflowを更新
for edge in path:
if edge in forwardEdges:
Gf[edge[0]][edge[1]]['flow'] += diff
# このとき、backward edgeのcapacityを更新する必要あり？
Gf[edge[1]][edge[0]]['capacity'] += diff
else:
Gf[edge[0]][edge[1]]['flow'] -= diff
Gf[edge[1]][edge[0]]['capacity'] -= diff
# もともと無かったedgeを消去
for edge in Gf.edges():
if not edge in forwardEdges:
Gf.remove_edge(edge[0],edge[1])
return Gf
In [35]:
def biMatch_ff(G,V,W):
'''Calculate the perfect matching of (G,V,W) by using Ford-Fulkerson algorithm.
Parameters
----------
G : bipartite graph
V,W : the lists of nodes, two groups of G
Return
------
G : modified graph with perfect matching
'''
H = add_st(G,V,W)
I = fordFulkerson(H, 's', 't')
I.remove_node('s')
I.remove_node('t')
return I
In [36]:
H = biMatch_ff(G,V,W)
In [37]:
print(H['v3']['w2']['flow'])
1.0
In [38]:
# 描画
pos={'v1':(0,10),'v2':(0,8),'v3':(0,6),'v4':(0,4),'v5':(0,2),'w1':(2,10),'w2':(2,8),'w3':(2,6),'w4':(2,4),'w5':(2,2),'w6':(2,0)}
edge_labels = {(i, j): (w['capacity'], str((w['flow'])) ) for i, j, w in H.edges(data=True)}
nx.draw_networkx_edge_labels(H, pos, edge_labels=edge_labels, font_color='r')
nx.draw_networkx_labels(H, pos)
nx.draw(H, pos)
plt.axis('off')
plt.show()
うまく動いている
In [40]:
# NetworkXのライブラリを用いて求めるとこうなる
flow_value, flows = nx.maximum_flow(G0,'s','t')
print('maximum flow: {}'.format(flow_value))
caps = nx.get_edge_attributes(G0, 'capacity')
for u in nx.topological_sort(G0):
for v, flow in sorted(flows[u].items()):
print('({}, {}): {}/{}'.format(u, v, flow, caps[(u, v)]))
maximum flow: 4 (s, v1): 0/1 (s, v2): 1/1 (s, v3): 1/1 (s, v4): 1/1 (s, v5): 1/1 (v3, w2): 1/1 (v1, w1): 0/1 (v1, w2): 0/1 (v4, w3): 1/1 (v4, w4): 0/1 (w3, t): 1/1 (w4, t): 0/1 (v2, w1): 1/1 (v2, w2): 0/1 (w2, t): 1/1 (w1, t): 1/1 (v5, w5): 0/1 (v5, w6): 1/1 (w5, t): 0/1 (w6, t): 1/1
In [ ]: