はじめに
栗林公園の概要や最速で巡るモチベーションについては、前回の記事で詳しく説明しているので、先にそちらをご覧ください。
その1では、栗林公園の北西部に絞って、全ての道を1回以上通るルートの中で最短になる手順を求めました。庭園全体で求めることができなかった理由として、マッチング問題の計算が非常に重く、園内全てを対象にすると計算時間が爆発的に増えてしまうためです。
この記事では、マッチング問題を解くアルゴリズムを変えることで計算量を改善して、園内全域を最短ルートを求めることを目指します。
手法
問題を解く考え方やコードの大部分は前回の記事のものを再利用して、マッチングアルゴリズムと地図データの取得に関する部分だけを変えます。マッチング問題を解くアルゴリズムは以下の3種類について実装します。１は準備運動、２は本命、３はおまけの感覚です。
- 動的計画法（bitDP）による解法
- Blossom algorithm
- PuLPを利用した解法
1.動的計画法
アルゴリズムの概要
この動的計画法では、以下の値を管理することで全てのマッチングの中で、距離の総和が最小の解を求めることができます。ただし、$x$は次数が奇数の頂点全体の集合の部分集合です。
$dp[x] := (xに含まれる要素をマッチングしたときの距離の最小値)$
マッチングは必ず2つの頂点をペアにするため、$x$の要素数が偶数の時のみ$dp[x]$は意味のある値を持ちます。遷移式は、2頂点間の距離を$dist(i, j)$とすると、以下のように書けます。
$dp[S \cup \{ i, j\}] = dp[S] + dist(i, j)$
$dp[\emptyset] = 0$からスタートして、徐々に大きな集合を考えていくことで、最終的に全ての頂点をマッチングした時の最小値をもとめることができます。このアルゴリズムの優れている点は、$[(1, 2), (3, 4)]$と$[(1, 3), (2, 4)]$のような同じ頂点を使った複数のマッチングについて、それらの中で距離が最小の結果だけを覚えておくことで無駄を省いて計算できる点です。全探索のアルゴリズムと比較して、このような重複を排除して計算を進めることができるため、計算量が改善できるというカラクリです。
実装では集合をそのまま管理するのは少し面倒なので、集合に含まれている要素のビットを1、含まれていない要素のビットを0で表現した2進数に変換する形で、集合を数値におきかえて管理します。このテクニックを用いた動的計画法のことを「bitDP」と呼びます。
実際にこの遷移式を用いて計算をすると、計算量は$O(N^2 2^N)$となり、前回の$O(N!)$のアルゴリズムと比べて$N$が大きな値でも高速に計算できることが分かります。
データ
計算量を踏まえると、$N=20$以下なら現実的な時間で解を求めることができそうです。そこで、このアルゴリズムでは公園の南側の以下の範囲で最適解を求めます。
実装
マッチングを求める部分を以下のコードに差し替えることで、最適なマッチングを求めることができます。距離の大きさをfloat型のままで計算すると、少し時間がかかってしまうためint型に変換して計算しています。
INF = 1 << 60
N = len(even_nodes)
dp = [INF] * (1 << N)
dp[0] = 0
m = [0] * (1<<N) # 動的計画法による結果の復元に使用するリスト
for i in tqdm.tqdm(range(1 << N)):
if dp[i] == INF:
continue
v = []
for j in range(N):
if (i >> j) % 2 == 0:
v.append(j)
for j in range(len(v)-1):
for k in range(j+1, len(v)):
vj = v[j]
vk = v[k]
n = i + (1<<vj) + (1<<vk)
if dp[n] > dp[i] + int(even_dists[vj][vk]*10000):
dp[n] = dp[i] + int(even_dists[vj][vk]*10000)
m[n] = [vj, vk] # 状態nでマッチングしたペアを保存
pairs = []
n = (1<<N) - 1
while n>0:
pairs.append((even_nodes[m[n][0]], even_nodes[m[n][1]]))
j, k = m[n]
n -= 1<<j
n -= 1<< k
結果
総移動距離は1776mでした。マッチングの計算にかかる時間は、81秒と予想に反して少し長かったです。公園の1/3くらいの大きさに広げることができましたが、グラフの大きさに対して指数的に計算時間が増えることを考えると、これ以上の拡張は厳しそうです。
2. Blossom algorithm
アルゴリズムの概要
Blossomアルゴリズムは、最大重み一般マッチング問題を$O(N^3)$で解くことができるアルゴリズムです。Edmondsのアルゴリズムとも呼ばれます。これまでのアルゴリズムとは異なり多項式オーダーで計算できるため、非常に高速で計算できますが、アルゴリズムの内容はかなり複雑で自力で実装するのも非常に困難です。
AtCoderでは、問題としてマッチング問題が扱われること自体少なく、出題されても二部グラフのマッチングです。過去5年の間でAtCoder Begginger Contestにおいて、このアルゴリズムが要求された例はないレベルであり、筆者自身も完全には理解できていないので、Blossomアルゴリズムの解説はここでは割愛します。
参考
GabowのEdmonds' Algorithm(一般マッチング O(VElogV))の解説と実装 - Qiita
Blossom algorithm - Wikipedia
データ
公園全体の交差点（頂点）の数はおおよそ300点なので、Blossomアルゴリズムを用いれば、十分高速に最適なマッチングを求めることができそうです。
実装
自力では実装できないため、Library Checkerからお借りしました。このアルゴリズムは最大重みマッチング求めるものなので、辺の重みを十分大きな値から距離を引き算した値にして、最大マッチングが総距離最小マッチングとなるように変換しています。また、float型だと計算が重くなってしまうため、値を切り捨ててint型で計算するようにしています。float型で計算した時と異なる結果が出る可能性がありますが、それでも合計距離の誤差は1m以内に収まるのでここでは許容することにします。
実装
INF = 1 << 60
from collections import deque
class GeneralWeightedMatching():
def __init__(self, n):
self.n = n
self.nx = n
self.m = 2 * n + 1
self.u = [0] * self.m * self.m
self.v = [0] * self.m * self.m
self.w = [0] * self.m * self.m
self.match = [0] * self.m
self.slack = [0] * self.m
self.flower = [[] for _ in range(self.m)]
self.flower_from = [0] * self.m * self.m
self.label = [0] * self.m
self.root = [0] * self.m
self.par = [0] * self.m
self.col = [0] * self.m
self.vis = [0] * self.m
self.que = deque()
self.t = 0
for u in range(1, self.m):
for v in range(1, self.m):
self.u[u * self.m + v] = u
self.v[u * self.m + v] = v
def dist(self, u, v):
u, v = self.u[u * self.m + v], self.v[u * self.m + v]
return self.label[u] + self.label[v] - self.w[u * self.m + v] * 2
def add_edge(self, u, v, w):
u += 1; v += 1
self.w[u * self.m + v] = max(self.w[u * self.m + v], w)
self.w[v * self.m + u] = max(self.w[v * self.m + u], w)
def update_slack(self, u, x):
if not self.slack[x] or self.dist(u, x) < self.dist(self.slack[x], x):
self.slack[x] = u
def set_slack(self, x):
self.slack[x] = 0
for u in range(1, self.n + 1):
if self.w[u * self.m + x] > 0 and self.root[u] != x and self.col[self.root[u]] == 0:
self.update_slack(u, x)
def que_push(self, x):
stack = [x]
while stack:
x = stack.pop()
if x <= self.n:
self.que.append(x)
continue
for i, fi in enumerate(self.flower[x]):
stack.append(fi)
def set_root(self, x, b):
stack = [x]
while stack:
x = stack.pop()
self.root[x] = b
if x <= self.n:
continue
for i, fi in enumerate(self.flower[x]):
stack.append(fi)
def get_pr(self, b, xr):
f = self.flower[b]
pr = f.index(xr)
if pr % 2 == 1:
f = self.flower[b] = f[0:1] + f[1:][::-1]
return len(f) - pr
else:
return pr
def set_match(self, u, v):
self.match[u] = self.v[u * self.m + v]
if u <= self.n:
return
xr = self.flower_from[u * self.m + self.u[u * self.m + v]]
pr = self.get_pr(u, xr)
f = self.flower[u]
for i in range(pr):
self.set_match(f[i], f[i ^ 1])
self.set_match(xr, v)
self.flower[u] = f[pr:] + f[:pr]
def augment(self, u, v):
xnv = self.root[self.match[u]]
self.set_match(u, v)
while xnv:
self.set_match(xnv, self.root[self.par[xnv]])
u, v = self.root[self.par[xnv]], xnv
xnv = self.root[self.match[u]]
self.set_match(u, v)
def get_lca(self, u, v):
self.t += 1
while u or v:
if not u:
u, v = v, u
continue
if self.vis[u] == self.t:
return u
self.vis[u] = self.t
u = self.root[self.match[u]]
if u: u = self.root[self.par[u]]
u, v = v, u
return 0
def add_blossom(self, u, lca, v):
b = self.n + 1
while b <= self.nx and self.root[b]:
b += 1
if b > self.nx:
self.nx += 1
self.label[b] = 0
self.col[b] = 0
self.match[b] = self.match[lca]
f = self.flower[b] = []
f.append(lca)
x = u
while x != lca:
f.append(x)
y = self.root[self.match[x]]
f.append(y)
self.que_push(y)
x = self.root[self.par[y]]
f = self.flower[b] = f[0:1] + f[1:][::-1]
x = v
while x != lca:
f.append(x)
y = self.root[self.match[x]]
f.append(y)
self.que_push(y)
x = self.root[self.par[y]]
self.set_root(b, b)
for x in range(1, self.nx + 1):
self.w[b * self.m + x] = self.w[x * self.m + b] = 0
for x in range(1, self.n + 1):
self.flower_from[b * self.m + x] = 0
for i, xs in enumerate(f):
for x in range(1, self.nx + 1):
if self.w[b * self.m + x] == 0 or self.dist(xs, x) < self.dist(b, x):
self.u[b * self.m + x] = self.u[xs * self.m + x]
self.u[x * self.m + b] = self.u[x * self.m + xs]
self.v[b * self.m + x] = self.v[xs * self.m + x]
self.v[x * self.m + b] = self.v[x * self.m + xs]
self.w[b * self.m + x] = self.w[xs * self.m + x]
self.w[x * self.m + b] = self.w[x * self.m + xs]
for x in range(1, self.n + 1):
if self.flower_from[xs * self.m + x]:
self.flower_from[b * self.m + x] = xs
self.set_slack(b)
def expand_blossom(self, b):
f = self.flower[b]
for i, fi in enumerate(f):
self.set_root(fi, fi)
xr = self.flower_from[b * self.m + self.u[b * self.m + self.par[b]]]
pr = self.get_pr(b, xr)
f = self.flower[b]
for i in range(0, pr, 2):
xs = f[i]
xns = f[i + 1]
self.par[xs] = self.u[xns * self.m + xs]
self.col[xs] = 1
self.col[xns] = 0
self.slack[xs] = 0
self.set_slack(xns)
self.que_push(xns)
self.col[xr] = 1
self.par[xr] = self.par[b]
for i in range(pr + 1, len(f)):
xs = f[i]
self.col[xs] = -1
self.set_slack(xs)
self.root[b] = 0
def on_found_edge(self, u, v):
eu = self.u[u * self.m + v]
ev = self.v[u * self.m + v]
u = self.root[eu]
v = self.root[ev]
if self.col[v] == -1:
self.par[v] = eu
self.col[v] = 1
nu = self.root[self.match[v]]
self.slack[v] = self.slack[nu] = 0
self.col[nu] = 0
self.que_push(nu)
elif self.col[v] == 0:
lca = self.get_lca(u, v)
if not lca:
self.augment(u, v)
self.augment(v, u)
return 1
else:
self.add_blossom(u, lca, v)
return 0
def matching(self):
for i in range(self.nx + 1):
self.col[i] = -1
self.slack[i] = 0
self.que.clear()
for x in range(1, self.nx + 1):
if self.root[x] == x and not self.match[x]:
self.par[x] = 0
self.col[x] = 0
self.que_push(x)
if not self.que:
return 0
while True:
while self.que:
u = self.que.popleft()
if self.col[self.root[u]] == 1:
continue
for v in range(1, self.n + 1):
if self.w[u * self.m + v] and self.root[u] != self.root[v]:
if self.dist(u, v) == 0:
if self.on_found_edge(u, v):
return 1
else:
self.update_slack(u, self.root[v])
d = INF
for b in range(self.n + 1, self.nx + 1):
if self.root[b] == b and self.col[b] == 1:
d = min(d, self.label[b] // 2)
for x in range(1, self.nx + 1):
if self.root[x] == x and self.slack[x]:
if self.col[x] == -1:
d = min(d, self.dist(self.slack[x], x))
elif self.col[x] == 0:
d = min(d, self.dist(self.slack[x], x) // 2)
for u in range(1, self.n + 1):
if self.col[self.root[u]] == 0:
if self.label[u] <= d:
return 0
self.label[u] -= d
elif self.col[self.root[u]] == 1:
self.label[u] += d
for b in range(self.n + 1, self.nx + 1):
if self.root[b] == b:
if self.col[b] == 0:
self.label[b] += d * 2
elif self.col[b] == 1:
self.label[b] -= d * 2
self.que.clear()
for x in range(1, self.nx + 1):
if self.root[x] == x and self.slack[x] and self.root[self.slack[x]] != x and self.dist(self.slack[x], x) == 0:
if self.on_found_edge(self.slack[x], x):
return 1
for b in range(self.n + 1, self.nx + 1):
if self.root[b] == b and self.col[b] == 1 and self.label[b] == 0:
self.expand_blossom(b)
return 0
def solve(self):
cnt = 0
ans = 0
for u in range(self.n + 1):
self.root[u] = u
self.flower[u].clear()
w_max = 0
for u in range(1, self.n + 1):
for v in range(1, self.n + 1):
self.flower_from[u * self.m + v] = u if u == v else 0
w_max = max(w_max, self.w[u * self.m + v])
for u in range(1, self.n + 1):
self.label[u] = w_max
while self.matching():
cnt += 1
for u in range(1, self.n + 1):
if self.match[u] and self.match[u] < u:
ans += self.w[u * self.m + self.match[u]]
for i in range(self.n):
self.match[i] = self.match[i + 1] - 1
return ans, cnt
N = len(dist)
G = GeneralWeightedMatching(N)
for i in range(N):
for j in range(N):
if i == j:
continue
G.add_edge(i, j, 10**9 - int(dist[i][j]*10000))
_ = G.solve()
pair = []
for i in range(N):
if G.match[i] > i:
pair.append((even[i], even[G.match[i]]))
結果
頂点のマッチング結果は以下のようになりました。
[(0, 1), (2, 65), (3, 173), (4, 138), (5, 139), (6, 135), (7, 26), (8, 92), (9, 10), (11, 81), (12, 54), (13, 14), (15, 79), (16, 61), (17, 18), (19, 20), (21, 117), (22, 23), (25, 113), (28, 112), (29, 69), (31, 200), (32, 180), (33, 174), (34, 62), (36, 74), (37, 38), (39, 60), (40, 41), (44, 57), (45, 46), (47, 133), (48, 50), (51, 109), (53, 134), (58, 59), (63, 64), (66, 67), (70, 73), (71, 141), (72, 75), (76, 111), (78, 176), (82, 97), (83, 84), (85, 175), (88, 184), (89, 110), (93, 100), (94, 123), (95, 96), (104, 116), (105, 161), (114, 115), (136, 137), (140, 142), (143, 144), (152, 153)]
マッチングの組の数から、期待通り全ての頂点がマッチングされています。前回のルート構築のコードを用いて最適ルートを可視化すると以下のようになります。結果の描画に際して、これまでのコードではある2頂点間を結ぶ道が複数ある場合に、距離が短い方しか通らないという問題が起こるため、1回目は距離が長い方を通り、2回目以降は距離が短い方を通るというルールを設けて描画しています1。
総移動距離は9724mです。徒歩（時速5km）で観覧した場合、約２時間かかるという結果でした。2回通る部分の総距離は1430mであるので、栗林公園の歩道の総和の約8.3kmに対して、余分に歩くことになる距離が1/6程度であることが分かります。
3. PuLPを利用した解法
PuLPは、pythonで混合整数線形計画問題を扱うことができるライブラリで、複数の連続(整数)変数の一次式で表現される目的関数や制約条件の最適化を行うことができます。このライブラリを用いると、目的関数の最小値とその時の変数の値をある程度高速に求解できるので、これまでの2つの解法のように、マッチング問題を解くための専用というわけではなく、巡回セールスマン問題やナップザック問題なども定式化の工夫で解くことができます。この節では一般マッチング問題を一次式の形に変形してPuLPで解いてもらおうという狙いになっています。
定式化
一般マッチング問題を、複数の変数を使った一次式の形で表現します。変数$x$を以下のように定義します。
\begin{eqnarray}
x_{ij}
=
\begin{cases}
1 & ( i,j\text{ がペア} ) \\
0 & ( i,j\text{ がペアでない} )
\end{cases}
\end{eqnarray}
ただし、奇数頂点数を$N$として$1 \leqq i,j \leqq N$を満たします。行列$x$は、マッチングのペアの成分にのみ$1$が入っているという意味になります。
最小化すべき目的関数は以下の通りになります。ペア間の距離の総和が最小になるようなマッチングを求めることが目的であるため、$x$を用いてその値を表現したものとなります。
\text{Objective} = \sum_{1 \leqq i,j \leqq N} dist(i, j) \times x_{ij}
行列$x$の制約条件は以下の4つです。$x$が一般マッチングの解として満たす必要がある性質を条件としています。条件1,2は「マッチングの相手が1つであること」、条件3は「自分自身とマッチングしないこと」、条件4は「マッチングの相方から見た相方が自分であること」を課しています。
- $x$の各行の和が$1$であること。$\sum_{j=1}^N x_{ij} = 1$ $(1 \leqq i \leqq N)$
- $x$の各列の和が$1$であること。$\sum_{i=1}^N x_{ij} = 1$ $(1 \leqq j \leqq N)$
- $x$の対角成分は$0$であること。$x_{ij} = 0$ $(i = j, 1 \leqq i, j \leqq N)$
- $x$の対象行列であること。$x_{ij} = x_{ji}$ $(1 \leqq i, j \leqq N)$
実装と結果
定式化で定義した目的関数と制約条件をPuLPの形式で設定して、solveを実行すると解を求めてくれます。
PuLPによる実装例
N = len(even_dists)
model = pp.LpProblem("matching", pp.LpMinimize)
x = [[pp.LpVariable("x(%s, %s)"%(i, j), cat="Binary") for i in range(N)] for j in range(N)]
objective = pp.lpSum([pp.lpDot(w, xs) for w, xs in zip(even_dists, x)])
model += objective
# 制約1: 行和が1
for i in range(N):
model += pp.lpSum(x[i][j] for j in range(N)) == 1
# 制約2: 列和が1
for j in range(N):
model += pp.lpSum(x[j][i] for j in range(N)) == 1
# 制約3: 対象行列
# 制約4: 対角成分が0
for i in range(N):
for j in range(i, N):
if i == j:
model += (x[i][j] == 0)
else:
model += (x[i][j] == x[j][i])
status = model.solve()
# 結果の出力
print("Status: {}".format(pp.LpStatus[status]))
print("Optimal Value: {}".format(objective.value() / 2))
# マッチングの結果の出力
result = []
for i in range(N):
for j in range(i, N):
if x[i][j].value() == 1:
result.append((even_nodes[i], even_nodes[j]))
print("Result:")
print(result)
Status: Optimal
Optimal Value: 1429.0725243928862
Result:
[(0, 1), (2, 65), (3, 173), (4, 138), (5, 139), (6, 135), (7, 26), (8, 92), (9, 10), (11, 81), (12, 54), (13, 14), (15, 79), (16, 61), (17, 18), (19, 20), (21, 117), (22, 23), (25, 113), (28, 112), (29, 69), (31, 200), (32, 180), (33, 174), (34, 62), (36, 74), (37, 38), (39, 60), (40, 41), (44, 57), (45, 46), (47, 133), (48, 50), (51, 109), (53, 134), (58, 59), (63, 64), (66, 67), (70, 73), (71, 141), (72, 75), (76, 111), (78, 176), (82, 97), (83, 84), (85, 175), (88, 184), (89, 110), (93, 100), (94, 123), (95, 96), (104, 116), (105, 161), (114, 115), (136, 137), (140, 142), (143, 144), (152, 153)]
2のBlossom algorithmで求めた結果と同じ結果が得られたので、結果の正しさが確認できたと思います。実行時間は7.0秒とBlossom algorithmと比べて少し長めでした。
まとめ
今回は、香川県高松市にある日本最大の文化財庭園「栗林公園」を最速で巡るために、公園内の道をすべて通りつつ、最短で回覧する道順を構成するプログラムを作成しました。特に、後編ではマッチングの最適化アルゴリズムに工夫を入れて、計算量を下げることで計算時間の短縮を行い、プログラムの適応範囲を公園全体で広げることができました。
園内全ての道を一度以上と通りたい場合は、総移動距離が約9.7kmで、休憩なしで2時間かかることが分かりました。ぜひ高松観光のスケジュールを組む時の参考にしていただけたらと思います。
最後に、この問題は「すべての道を一度以上通る最適なルート」を計算することができ、一般的な配送の最適化問題のような「すべての点を通る最適なルート」とは本質的に異なる問題であるため、当然有用となるシチュエーションも異なります。例としては道路の点検業務のルートを作るときとかに利用可能だと思います。
さらに、都会の住宅地などでは「すべての道を通る」という要請を「すべての道沿いにある家を通る」とおおよそ読み換えることができるので、ポスティングや国勢調査の訪問でも有効だと思われます。
-
前回の記事で説明したように、一筆書きの順序には任意性がありますので、どの順番で通っても問題はないです。 ↩