最小全域木問題を解いてみる.
まず,描画などのためにNetworkXとMatplotlibをそれぞれimportする.
import networkx as nx
import matplotlib.pyplot as plt
最小全域木問題の問題例として頂点集合$\{1,2,3,4,5\}$,枝集合$\{\{1,2\}, \ \{1,4\}, \ \{2,3\}, \ \{2,4\}, \ \{2,5\}, \ \{3,5\}, \ \{4,5\}\}$の無向グラフを用意する.
枝にはweightも付加する. また,描画のため,頂点には座標を付加する.
G = nx.Graph()
G.add_node(1, position=(0, 1))
G.add_node(2, position=(1, 2))
G.add_node(3, position=(3, 2))
G.add_node(4, position=(1, 0))
G.add_node(5, position=(3, 0))
G.add_edge(1, 2, weight=2)
G.add_edge(1, 4, weight=1)
G.add_edge(2, 3, weight=3)
G.add_edge(2, 4, weight=1)
G.add_edge(2, 5, weight=5)
G.add_edge(3, 5, weight=4)
G.add_edge(4, 5, weight=3)
なお,この問題例は
久保幹雄,組合せ最適化とアルゴリズム,共立出版,2000年.
の図1.8からの引用である.
枝重み付き無向グラフを描画する関数を以下に定義する.
def draw_graph_with_weight(G, node_position, title='A network', with_labels=True):
plt.title(title)
nx.draw_networkx(G, pos=node_position, with_labels=with_labels, node_color='w')
weight = nx.get_edge_attributes(G, 'weight')
edge_labels = nx.draw_networkx_edge_labels(G, pos=node_position, edge_labels=weight)
return
先程用意した問題例のグラフを描画して確認する.
%matplotlib inline
node_position = nx.get_node_attributes(G, 'position')
draw_graph_with_weight(G, node_position, 'An instance of the minimum spanning tree problem')
NetworkXには最小全域木問題を解いてくれる関数minimum_spanning_treeが用意されているので,それを使ってみる.
なお,この関数minimum_spanning_treeではKruskalのアルゴリズムが使われているようだ.
mst = nx.minimum_spanning_tree(G)
mst # 最適解を確認する.
<networkx.classes.graph.Graph at 0x10a0becf8>
NetworkXの関数minimum_spanning_treeは最適解をGraph形式で返す.
描画して確認する.
%matplotlib inline
plt.figure(figsize=(16,4))
plt.subplot(1, 2, 1)
draw_graph_with_weight(G, node_position, 'An instance of the minimum spanning tree problem')
plt.subplot(1, 2, 2)
draw_graph_with_weight(mst, node_position, 'A minimum spanning tree')
最小重みであることを確認するのは簡単ではないかもしれないが,全域木であることは確認できる.
次に,NetworkXをデータ構造として使いつつ,Kruskalのアルゴリズムを実装する.
まず,union findデータ構造は使わず,閉路ができてしまうか否かの判定はNetworkXの関数node_connected_componentを使う. この関数node_connected_componentは指定した頂点を含む連結成分を返す関数なので,事実上深さ優先探索(あるいは幅優先探索)で閉路ができてしまうか否かを判定するのと同じである(と思う).
def kruskal_simple(G):
edge = []
for u, v in G.edges():
edge.append((G.edge[u][v]['weight'], (u, v)))
edge.sort()
T = nx.Graph()
for v in G.nodes():
T.add_node(v, position=G.node[v]['position'])
for w, e in edge:
if len(T.edges()) - 1 == len(G.nodes()):
break
u, v = e
if v in nx.node_connected_component(T, u):
continue
T.add_edge(u, v, weight=w)
return T
実装したKruskalのアルゴリズムを実行してみる.
T = kruskal_simple(G)
得られた最小全域木と問題例を描画してみる.
%matplotlib inline
plt.figure(figsize=(16,4))
plt.subplot(1, 2, 1)
draw_graph_with_weight(G, node_position, 'An instance of the minimum spanning tree problem')
plt.subplot(1, 2, 2)
draw_graph_with_weight(T, node_position, 'A minimum spanning tree')
確かに最小全域木が求められているようだ.
次に,union findデータ構造の1つであるunion find treeを用いて,Kruskalのアルゴリズムを実装する.
def make_initial_union_find_tree(G):
union_find_tree = nx.DiGraph()
for v in G.nodes():
union_find_tree.add_node(v, size=1)
return union_find_tree
def name_and_size(v, union_find_tree):
pred = union_find_tree.predecessors(v)
if len(pred) == 0:
return v, union_find_tree.node[v]['size']
else:
return name_and_size(pred[0], union_find_tree)
def merge_tree(u, v, union_find_tree):
u_root_name, u_tree_size = name_and_size(u, union_find_tree)
v_root_name, v_tree_size = name_and_size(v, union_find_tree)
if u_root_name == v_root_name:
return
if u_tree_size > v_tree_size:
union_find_tree.add_edge(u_root_name, v_root_name)
union_find_tree.node[u_root_name]['size'] += v_tree_size
else:
union_find_tree.add_edge(v_root_name, u_root_name)
union_find_tree.node[v_root_name]['size'] += u_tree_size
return
def kruskal_union_find(G):
union_find_tree = make_initial_union_find_tree(G)
edge = []
for u, v in G.edges():
edge.append((G.edge[u][v]['weight'], (u, v)))
edge.sort()
T = nx.Graph()
for v in G.nodes():
T.add_node(v, position=G.node[v]['position'])
for w, (u, v) in edge:
if len(T.edges()) - 1 == len(G.nodes()):
break
if name_and_size(u, union_find_tree) == name_and_size(v, union_find_tree):
continue
T.add_edge(u, v, weight=w)
merge_tree(u, v, union_find_tree)
return T
T = kruskal_union_find(G)
得られた最小全域木と問題例を描画してみる.
%matplotlib inline
plt.figure(figsize=(16,4))
plt.subplot(1, 2, 1)
draw_graph_with_weight(G, node_position, 'An instance of the minimum spanning tree problem')
plt.subplot(1, 2, 2)
draw_graph_with_weight(T, node_position, 'A minimum spanning tree')
確かに最小全域木が得られているようだ.
次に,union find treeをより理解するため,Kruskalのアルゴリズムにおける解とunion find treeの変化を描画してみる.
先程定義した関数kruskal_union_findを,途中経過観察可能なように再定義する.
def kruskal_union_find2(G):
union_find_tree = [make_initial_union_find_tree(G)]
edge = []
for u, v in G.edges():
edge.append((G.edge[u][v]['weight'], (u, v)))
edge.sort()
T = [nx.Graph()]
for v in G.nodes():
T[-1].add_node(v, position=G.node[v]['position'])
for w, (u, v) in edge:
if len(T[-1].edges()) - 1 == len(G.nodes()):
break
if name_and_size(u, union_find_tree[-1]) == name_and_size(v, union_find_tree[-1]):
continue
T.append(T[-1].copy())
T[-1].add_edge(u, v, weight=w)
union_find_tree.append(union_find_tree[-1].copy())
merge_tree(u, v, union_find_tree[-1])
return T, union_find_tree
再定義したkruskal_union_findを実行してみる.
T, UFT = kruskal_union_find2(G)
実行結果(最小全域木の生成過程とunion find treeの成長過程)を描画する関数を以下に定義する.
def draw_forest(instance, forest, node_position, title='A forest', with_labels=True):
nx.draw_networkx(instance, node_position, with_labels=with_labels, node_color='w', alpha=0.5)
plt.title(title)
nx.draw_networkx(forest, pos=node_position, with_labels=with_labels, edge_color='r', node_color='w')
weight = nx.get_edge_attributes(instance, 'weight')
edge_labels = nx.draw_networkx_edge_labels(G, pos=node_position, edge_labels=weight)
return
def draw_forest_and_union_find_tree(instance, forest, union_find_tree, with_labels=True):
node_position = nx.get_node_attributes(instance, 'position')
total_phase = len(forest)
plt.figure(figsize=(16, 4 * total_phase))
for i in range(total_phase):
plt.subplot(total_phase, 2, 2 * i + 1)
draw_forest(instance, forest[i], node_position, title='A forest {0}'.format(i), with_labels=False)
plt.subplot(total_phase, 2, 2 * i + 2)
draw_graph_with_weight(union_find_tree[i], node_position, title='A union find tree {0}'.format(i), with_labels=with_labels)
return
Kruskalのアルゴリズムの途中経過を描画してみる.
draw_forest_and_union_find_tree(G, T, UFT, node_position)
小さな問題例ではunion find treeの変化がわかりにくいかもしれないので,もう少し大きな例を作って見てみる.
枝重みを一様ランダムに決めた格子状のグラフを用意する.
import random
def make_random_grid_instance(m, n, random_seed=1):
random.seed(random_seed)
grid_graph = nx.grid_2d_graph(m, n)
for i in range(m):
for j in range(n):
grid_graph.node[(i, j)]['position'] = (i, j)
for u, v in grid_graph.edges():
grid_graph.edge[u][v]['weight'] = random.randint(1, 99)
return grid_graph
g33 = make_random_grid_instance(3, 3)
%matplotlib inline
plt.figure(figsize=(8,4))
g33_node_position = nx.get_node_attributes(g33, 'position')
draw_graph_with_weight(g33, g33_node_position, 'An instance of the minimum spanning tree problem', with_labels=False)
t33, u33 = kruskal_union_find2(g33)
draw_forest_and_union_find_tree(g33, t33, u33, with_labels=False)
import time
import numpy as np
def mst_time(max_size=10):
log = []
for m in range(2, max_size + 1):
gmn = make_random_grid_instance(m, m)
size = len(gmn.nodes())
start_time = time.time()
mst_mn = kruskal_simple(gmn)
simple_kruskal_time = time.time() - start_time
start_time = time.time()
mst_mn = kruskal_union_find(gmn)
kruskal_union_find_time = time.time() - start_time
log.append((size, simple_kruskal_time, kruskal_union_find_time))
return log
log = np.array(mst_time(40))
%matplotlib inline
plt.figure(figsize=(16, 8))
plt.title("CPU time of Kruskal's algorithm")
plt.xlabel('instance size (# of nodes)')
plt.ylabel('CPU time')
plt.yscale('linear')
plt.plot(log[:,0], log[:,1], 'o-', label='simple Kruskal')
plt.plot(log[:,0], log[:,2], 'o-', label='Kruskal with UFT')
z = plt.legend(loc='upper left')
ランダムに問題例を作ると,単純な実装でもそんなに悪くない.
単純な実装に対して意地悪な問題例を作って試してみる.
def make_intractable_instance(n):
g = nx.Graph()
g.add_node(0, position=((1+n)/2,1))
for i in range(1, n + 1):
g.add_node(i, position=(0,i))
g.add_edge(0, 1, weight=1)
for i in range(1, n):
g.add_edge(i, i + 1, weight=i*2)
g.add_edge(0, i + 1, weight=i*2+1)
return g
def mst_time(max_size=10):
log = []
for n in range(1, max_size + 1, 100):
g = make_intractable_instance(n)
size = len(g.nodes())
start_time = time.time()
mst_n = kruskal_simple(g)
simple_kruskal_time = time.time() - start_time
start_time = time.time()
mst_n = kruskal_union_find(g)
kruskal_union_find_time = time.time() - start_time
log.append((size, simple_kruskal_time, kruskal_union_find_time))
return log
log = np.array(mst_time(2000))
%matplotlib inline
plt.figure(figsize=(16, 8))
plt.title("CPU time of Kruskal's algorithm")
plt.xlabel('instance size (# of nodes)')
plt.ylabel('CPU time')
plt.yscale('linear')
plt.plot(log[:,0], log[:,1], 'o-', label='simple Kruskal')
plt.plot(log[:,0], log[:,2], 'o-', label='Kruskal with UFT')
z = plt.legend(loc='upper left')