M 個の整数 1, 2, ..., M の順列を考えます。このとき、i 番目 (先頭要素が 1 番目) の要素が整数 i ではない順列を「完全順列 (derangement)」といいます。今回は 1 から M までの整数値で完全順列を生成する述語を clpfd で作ってみましょう。
プログラムは簡単です。順列を生成する処理で、N 番目の要素は N と等しくない、という制約を加えるだけです。次のリストを見てください。
リスト : 完全順列
:- use_module(library(clpfd)).
check([], _).
check([X | Xs], N) :-
X #\= N,
N1 #= N + 1,
check(Xs, N1).
derangement(N, Xs) :-
length(Xs, N),
Xs ins 1 .. N,
all_different(Xs),
check(Xs, 1),
label(Xs).
述語 check で N 番目の要素が N と等しくない制約を加えています。check の第 1 引数がリストで、第 2 引数の N が要素の位置を表します。制約は簡単で、先頭要素 X が N と等しくないこと、つまり X #\= N とすればいいだけです。
それでは実行してみましょう。
?- derangement(3, Xs). Xs = [2, 3, 1] ; Xs = [3, 1, 2]. ?- derangement(4, Xs). Xs = [2, 1, 4, 3] ; Xs = [2, 3, 4, 1] ; Xs = [2, 4, 1, 3] ; Xs = [3, 1, 4, 2] ; Xs = [3, 4, 1, 2] ; Xs = [3, 4, 2, 1] ; Xs = [4, 1, 2, 3] ; Xs = [4, 3, 1, 2] ; Xs = [4, 3, 2, 1]. ?-
正常に動作していますね。ところで、述語 check は numlist/3 と maplist/3 を使うともっと簡単に定義することができます。
numlist/3 は整数列を生成する述語です。
numlist(Low, High, Xs)
numlist は Low から High までの整数列をリストに格納して返します。簡単な使用例を示しましょう。
?- numlist(1, 1, Xs). Xs = [1]. ?- numlist(1, 0, Xs). false. ?- numlist(1, 8, Xs). Xs = [1, 2, 3, 4, 5, 6, 7, 8].
プログラムの修正は簡単で、check を numlist と maplist に置き換えるだけです。
リスト : check の置き換え
derangement(N, Xs) :-
length(Xs, N),
Xs ins 1 .. N,
all_different(Xs),
numlist(1, N, Ys),
maplist(#\=, Xs, Ys),
label(Xs).
numlist で 1 から N までの整数列 Ys を生成し、それを maplist に渡します。maplist は Xs の要素 (変数) と Ys の要素 (位置) を述語 #\= に渡すだけです。これで完全順列を生成することができます。
「ラテン方陣 (Latin square)」はナンバープレースの枠の条件を無くした方陣です。ラテン方陣の定義を参考文献『数理パズルのはなし』より引用します。 『ラテン方陣を一般的にいうなら、n 行 n 列の正方形の枡に n 種類の記号を n 個ずつ配列して、各行各列に記号の重複のないものを n 次のラテン方陣というのです。』 このラテン方陣をパズルに応用したものがナンバープレースというわけです。
簡単な例を示しましょう。3 次のラテン方陣は次に示す 12 通りになります。
0 1 2 0 1 2 0 2 1 0 2 1 1 0 2 1 0 2
1 2 0 2 0 1 1 0 2 2 1 0 0 2 1 2 1 0
2 0 1 1 2 0 2 1 0 1 0 2 2 1 0 0 2 1
標準形
1 2 0 1 2 0 2 0 1 2 0 1 2 1 0 2 1 0
0 1 2 2 0 1 0 1 2 1 2 0 0 2 1 1 0 2
2 0 1 0 1 2 1 2 0 0 1 2 1 0 2 0 2 1
図 : 3 次のラテン方陣
この中で、最初の行と列の要素を昇順に並べたものを「標準形」といいます。3 次のラテン方陣の場合、標準形は 1 種類しかありません。ラテン方陣は任意の行を交換する、または任意の列を交換してもラテン方陣になります。3 次のラテン方陣の場合、標準形から行または列を交換することで、残りの 11 種類のラテン方陣を生成することができます。
今回は標準形ラテン方陣を求めるプログラムを clpfd で作ってみましょう。プログラムは次のようになります。
リスト : ラテン方陣
:- use_module(library(clpfd)).
% 行列の生成
make_matrix(0, _, []).
make_matrix(N, M, [Xs | Ys]) :-
N > 0,
length(Xs, M),
N1 is N - 1,
make_matrix(N1, M, Ys).
check(N, Xs) :- Xs ins 1..N.
latin(N, [Xs1 | Xs]) :-
make_matrix(N, N, [Xs1 | Xs]),
maplist(check(N), Xs),
maplist(all_different, Xs),
transpose([Xs1 | Xs], [Ys1 | Ys]),
maplist(all_different, Ys),
numlist(1, N, Zs),
maplist(#=, Xs1, Zs),
maplist(#=, Ys1, Zs),
maplist(label, Xs).
述語 make_matrix は再帰定義で N 行 M 列の行列を生成します。行列の要素はすべて自由変数です。なお、再帰定義を使わずに行列を生成することもできます。次のリストを見てください。
リスト : 行列の生成
xlength(N, Xs) :- length(Xs, N).
make_matrix(N, M, Xs) :-
length(Xs, N),
maplist(xlength(M), Xs).
xlength/2 は length/2 の引数の順番を逆にしたものです。make_matrix は、最初に N 個の自由変数を格納したリストを length で Xs に生成します。次に、maplist/2 で Xs に xlenght(M) を適用すると、Xs の要素は自由変数なので、それが xlength(M) に渡されることになります。その結果、自由変数は長さ M のリスト (要素は自由変数) になるわけです。
自由変数の使い方が Prolog らしくてクールですね。M.Hiroi は参考 URL 6 『CLP(FD) Constraint Logic Programming over Finite Domains』でこの方法を知りました。作者様に感謝いたします。これで行列を生成できるのですから、Prolog はやっぱり面白いなあと思いました。
次に、制約を定義します。先頭の行 Xs1 の値は 1 から N に確定するので、それ以降の行 Xs に対して制約を定義します。述語 check は変数の範囲 (1 .. N) を指定します。そして、all_different で変数の値がすべて異なることを指定します。
列の制約を定義する場合、clpfd の述語 transpose/ 2 で転置行列を求めると簡単です。次の図を見てください。
[[1,2,3], [[1,4,7], [4,5,6], = 転置行列 => [2,5,8], [7,8,9]] [3,6,9]]
このように、行列の行と列を入れ替えた行列を「転置行列 (transposed matrix)」といいます。transpose(Xs, Ys) は引数 Xs を行列とみなして、その転置行列を求めます。簡単な使用例を示します。
?- transpose([[1,2],[3,4]], Xs). Xs = [[1, 3], [2, 4]]. ?- transpose([[1,2],[3,4],[5,6]], Xs). Xs = [[1, 3, 5], [2, 4, 6]]. ?- transpose([[1,2,3],[4,5,6],[7,8,9]], Xs). Xs = [[1, 4, 7], [2, 5, 8], [3, 6, 9]]. ?-
先頭の列 Ys1 は値が 1 から N に確定しているので、それ以降の列 Ys に対して、all_different で変数の値がすべて異なることを指定します。あとは、numlist と maplist で、Xs1 と Ys1 に 1 から N までの値を指定します。最後に、Xs に格納されている変数の値を label で決定します。
それでは実行してみましょう。
?- latin(3, Xs), maplist(writeln, Xs). [1,2,3] [2,3,1] [3,1,2] Xs = [[1, 2, 3], [2, 3, 1], [3, 1, 2]] ; false. ?- latin(4, Xs), maplist(writeln, Xs). [1,2,3,4] [2,1,4,3] [3,4,1,2] [4,3,2,1] Xs = [[1, 2, 3, 4], [2, 1, 4, 3], [3, 4, 1, 2], [4, 3, 2, 1]] ; [1,2,3,4] [2,1,4,3] [3,4,2,1] [4,3,1,2] Xs = [[1, 2, 3, 4], [2, 1, 4, 3], [3, 4, 2, 1], [4, 3, 1, 2]] ; [1,2,3,4] [2,3,4,1] [3,4,1,2] [4,1,2,3] Xs = [[1, 2, 3, 4], [2, 3, 4, 1], [3, 4, 1, 2], [4, 1, 2, 3]] ; [1,2,3,4] [2,4,1,3] [3,1,4,2] [4,3,2,1] Xs = [[1, 2, 3, 4], [2, 4, 1, 3], [3, 1, 4, 2], [4, 3, 2, 1]] ; false. ?-
正常に動作していますね。ちなみに、標準形ラテン方陣の総数は次のようになります。
I4 = 4 I5 = 56 I6 = 9408 I7 = 16942080
?- time((findall(Xs, latin(6, Xs), Ys), length(Ys, N), write(N))). 9408 % 21,469,021 inferences, 1.283 CPU in 1.283 seconds (100% CPU, 16729672 Lips) Ys = [[[1, 2, 3, 4, 5, 6], [2, 1, 4, 3, 6, 5], [3, 4, 5, 6, 1|...], [4, 3, 6, 5|...], [5, 6, 1|...], [6, 5|...]], [[1, 2, 3, 4, 5, 6], [2, 1, 4, 3, 6|...], [3, 4, 5, 6|...], [4, 3, 6|...], [5, 6|...], [6|...]], [[1, 2, 3, 4, 5|...], [2, 1, 4, 3|...], [3, 4, 5|...], [4, 3|...], [5|...], [...|...]], [[1, 2, 3, 4|...], [2, 1, 4|...], [3, 4|...], [4|...], [...|...]|...], [[1, 2, 3|...], [2, 1|...], [3|...], [...|...]|...], [[1, 2|...], [2|...], [...|...]|...], [[1|...], [...|...]|...], [[...|...]|...], [...|...]|...], N = 9408. ?-
このプログラムは 6 次の標準形ラテン方陣の総数を求めるのに約 1.3 秒かかりました。プログラムは簡単に作成できましたが、制約条件がゆるいせいか実行時間は遅いですね。高次の標準形ラテン方陣の総数は、簡単に求めることができない非常にハードな問題だといわれています。興味のある方は挑戦してみてください。
今回は N Queens Problem を解いてみましょう。最初に Prolog のプログラムを示します。
リスト : N Queens Problem
not_attack(_, [], _).
not_attack(X, [Y | Ys], N) :-
abs(X - Y) =\= N,
N1 is N + 1,
not_attack(X, Ys, N1).
safe([_]).
safe([X | Xs]) :- not_attack(X, Xs, 1), safe(Xs).
nqueens(N, Ys) :-
numlist(1, N, Xs),
permutation(Xs, Ys),
safe(Ys).
プログラムの説明は拙作のページ「Prolog 入門: 8 クイーン」をお読みください。実行結果は次のようになります。
?- nqueens(8, Xs). Xs = [1, 5, 8, 6, 3, 7, 2, 4] ; Xs = [1, 6, 8, 3, 7, 4, 2, 5] ; Xs = [1, 7, 4, 6, 8, 2, 5, 3] ; Xs = [1, 7, 5, 8, 2, 4, 6, 3] . ?- time((nqueens(10, Q), fail)). % 87,872,830 inferences, 10.289 CPU in 10.289 seconds (100% CPU, 8540201 Lips) false. ?-
このプログラムは単純な生成検定法なので、クイーンの個数が増えると実行時間は極端に遅くなります。
clpfd で N Queens Problem を解く場合、述語 not_attack の演算子 =\= と is を #\= と #= に変更します。プログラムは次のようになります。
リスト ; N Queens Problem (clpfd 版)
:- use_module(library(clpfd)).
not_attack(_, [], _).
not_attack(X, [Y | Ys], N) :-
abs(X - Y) #\= N,
N1 #= N + 1,
not_attack(X, Ys, N1).
safe([_]).
safe([X | Xs]) :- not_attack(X, Xs, 1), safe(Xs).
nqueens(N, Xs) :-
length(Xs, N),
Xs ins 1 .. N,
all_different(Xs),
safe(Xs).
length/2 で無名変数を格納したリスト Xs を生成します。あとは、Xs の順列を生成し、パズルの条件を満たすように述語 safe で制約を記述するだけです。
実行結果は次のようになりました。
?- nqueens(8, Xs), label(Xs). Xs = [1, 5, 8, 6, 3, 7, 2, 4] ; Xs = [1, 6, 8, 3, 7, 4, 2, 5] ; Xs = [1, 7, 4, 6, 8, 2, 5, 3] ; Xs = [1, 7, 5, 8, 2, 4, 6, 3] . ?- time((nqueens(10, Xs), label(Xs), fail)). % 10,570,645 inferences, 0.721 CPU in 0.721 seconds (100% CPU, 14666479 Lips) false. ?-
単純な生成検定法よりも約 14 倍速くなりました。制約を記述するだけでここまで高速化するとは、clpfd のパワーに M.Hiroi も大変驚きました。
ただし、N Queens Problem の場合、無駄な順列を生成しないように枝刈りを入れれば、clpfd 版よりも高速になります。
リスト : N Queens Problem (高速版)
nqueens1(Xs, Q) :- queen_sub(Xs, [], Q).
queen_sub(L, SafeQs, Q) :-
select(X, L, RestQs),
not_attack(X, SafeQs, 1),
queen_sub(RestQs, [X | SafeQs], Q).
queen_sub([], Q, Q).
順列を生成している途中で、クイーンが衝突していないか not_attack でチェックしています。実行結果は次のようになりました。
?- time((nqueens1([1,2,3,4,5,6,7,8,9,10], Xs), fail)). % 1,075,402 inferences, 0.135 CPU in 0.135 seconds (100% CPU, 7974961 Lips) false. ?-
clpfd 版よりも 5 倍以上速くなりました。
なお、解をひとつ求めるだけでよければ、clpfd 版も負けてはいません。次の例を見てください。
?- time((numlist(1, 20, Qs), nqueens1(Qs, Xs))). % 18,869,293 inferences, 2.414 CPU in 2.414 seconds (100% CPU, 7816889 Lips) Qs = [1, 2, 3, 4, 5, 6, 7, 8, 9|...], Xs = [11, 6, 14, 7, 10, 8, 19, 16, 9|...] . ?-
?- time((nqueens(20, Xs), labeling([ff], Xs))). % 166,349 inferences, 0.014 CPU in 0.014 seconds (100% CPU, 12091690 Lips) Xs = [1, 3, 5, 14, 17, 4, 16, 7, 12|...] . ?- time((nqueens(30, Xs), labeling([ff], Xs))). % 249,375 inferences, 0.027 CPU in 0.027 seconds (100% CPU, 9280209 Lips) Xs = [1, 3, 5, 24, 26, 4, 23, 7, 28|...] . ?- time((nqueens(40, Xs), labeling([ff], Xs))). % 366,900 inferences, 0.035 CPU in 0.035 seconds (100% CPU, 10628252 Lips) Xs = [1, 3, 5, 26, 33, 4, 28, 7, 34|...] . ?-
nqueens1 の場合、クイーンの個数が増えると最初の解を求めるのに時間がかかるようになりますが、clpfd 版では labeling/2 のオプションに ff を指定すると、最初の解を高速に求めることができます。興味のある方はいろいろ試してみてください。
「ビンパッキング問題 (Bin Packing Problem)」は、大きさが異なる N 個の品物を大きさ B のビンに詰めるとき、最小のビンの本数と品物の詰め方を求める問題です。ビンは品物を入れる箱と考えてください。品物の大きさは B 以下で、ビンも品物を詰めるだけの本数は十分に用意されているものとします。
ビンパッキング問題は部分和問題と同じく NP 問題になるので、品物の数が多くなると最適解を求めるのが大変難しくなります。そこで今回は、ビンの本数をあらかじめ決めておいて、品物の詰め方を求めるプログラムを作ってみましょう。それでは問題です。
4 人で 8 個の荷物を運びます。各荷物の重さは 3.3 kg, 6.1 kg, 5.8 kg, 4.1 kg, 5.0 kg, 2.1 kg, 6.0 kg, 6.4 kg です。各自の運ぶ荷物の重さの合計が 11 kg 以下になるように荷物を割り当てることはできるでしょうか。できる場合はその割り当て方を求めてください。
出典 : Coprisによる制約プログラミング入門, (田村直之さん)
clpfd は整数しか扱えないので、重さを 10 倍して整数に直して解くことにします。このような問題を解く場合、データ構造に「結合行列」を使うと簡単です。次の図を見てください。
: a b c d e f g h
---+-----------------
A : 1 0 1 0 0 0 0 0
B : 0 0 0 1 0 0 1 0
C : 0 1 0 0 1 0 0 0
D : 0 0 0 0 0 1 0 1
図 : 結合行列
a - h が品物、A = D が人とします。結合行列は人が運ぶ品物を 1 で、運ばない品物を 0 で表します。上図の場合、A が運ぶ品物は a, c で、D が運ぶ品物が f, h となります。
制約の記述も簡単です。一人が運ぶ品物の重さの制限は scalar_product/4 を使えば簡単ですね。もう一つ制約があって、品物を表す列には 1 が一つしかないことです。この場合、clpfd の述語 sum/3 を使うと簡単に制約を表すことができます。
sum(Xs, Pred, N)
sum/3 は Xs の総和を求めて、その値と引数 N が述語 Pred を満たすときに成功します。なお、sum で使用できる Pred は、#=, #\=, #<, #>, #=<, #>= だけです。たとえば、リスト Xs の要素が 0 と 1 だけの場合、sum(Xs, #=, 1) とすれば、Xs には 1 が一つしかないという制約を表すことができます。
プログラムは次のようになります。
リスト : ビンパッキング問題
:- use_module(library(clpfd)).
% N 行 M 列の行列を生成
xlength(N, Xs) :- length(Xs, N).
make_matrix(N, M, Xs) :-
length(Xs, N),
maplist(xlength(M), Xs).
% 行の制約
check_line(Ws, W, Xs) :-
Xs ins 0..1,
scalar_product(Ws, Xs, #=<, W).
% 列の制約
check_column(Xs) :- sum(Xs, #=, 1).
% 解法
solver :-
Ws = [33, 61, 58, 41, 50, 21, 60, 64],
make_matrix(4, 8, Xs),
maplist(check_line(Ws, 110), Xs),
transpose(Xs, Ys),
maplist(check_column, Ys),
maplist(label, Xs),
maplist(writeln, Xs).
make_matrix は N 行 M 列の行列を作ります。要素は自由変数になります。行の制約は述語 check_line で定義します。引数 Ws は品物の重さのリスト、W が重さの上限値、Xs が結合行列の行です。ins で Xs の変数の値が 0 or 1 になることと、重さの合計値が W 以下であることを scalar_product で設定します。
列の制約は述語 check_column で行います。Xs は列を表すリストで、sum(Xs, #=, 1) とするだけです。列を求める場合、clpfd の述語 transpose/ 2 で転置行列を求めると簡単です。あとは、solver で maplist に check_line と check_column を渡して制約を定義し、label で変数の値を求めます。最後に、maplist(writeln, Xs) で行列を表示するだけです。
それでは実行してみましょう。
?- solver. [0,0,0,0,0,1,0,1] [0,0,0,0,1,0,1,0] [0,0,1,1,0,0,0,0] [1,1,0,0,0,0,0,0] true ; [0,0,0,0,0,1,0,1] [0,0,0,0,1,0,1,0] [0,1,0,1,0,0,0,0] [1,0,1,0,0,0,0,0] true ; [0,0,0,0,0,1,0,1] [0,0,0,0,1,0,1,0] [1,0,1,0,0,0,0,0] [0,1,0,1,0,0,0,0] true . ?-
最初の解では、A が 8.5 kg, B が 11.0 kg, C が 9.9 kg, D が 9.4 kg となりました。総重量 (38.8 kg) を 4 人で割ると 9.7 kg になりますが、これを上限値とすると、解はなくなります。4 人の場合、解があるのは 10.8 kg までです。興味のある方はいろいろ試してみてください。
次は下図に示す経路で、A から I までの最短経路を求めてみましょう。なお、今回は各地点間の距離を 1 とすることにします。
Ａ───Ｂ───Ｃ
│ │ │
│ │ │
│ │ │
Ｄ───Ｅ───Ｆ
│ │ │
│ │ │
│ │ │
Ｇ───Ｈ───Ｉ
図 : 経路図
制約プログラミングで経路を求める場合、辺 (edge) を中心に考えると簡単です。グラフには辺に向きがある有向グラフと、向きがない無向グラフの 2 種類ありますが、今回は両方でプログラムを作ってみましょう。
なお、最短経路問題は「ダイクストラのアルゴリズム」を使って高速に解くことができます。興味のある方は拙作のページ「Algorithms with Python: 欲張り法」をお読みください。
最初は無向グラフでプログラムを作ります。基本的な考え方は簡単です。辺を変数に対応させ、ある辺を選んだらその変数の値を 1 に、選ばなかったら 0 にします。辺と変数の対応を下表に示します。
表 : 辺と変数の対応
辺 : 変数
-------+------
A - B : AB
A - D : AD
B - C : BC
B - E : BE
C - F : CF
D - E : DE
D - G : DG
E - F : EF
E - H : EH
F - I : FI
G - H : GH
H - I : HI
あとは、選んだ辺がスタートからゴールまでつながることを表わせばいいわけです。この制約も簡単です。各頂点において、そこに接続されている辺の本数をチェックするだけです。
スタート A とゴール I は、辺が 1 本だけ接続されています。経路上にある頂点、たとえば A - B - C という経路において、A から B へ進むとき辺 AB を通り、B から C へ進むときは辺 BC を通りますね。このように、経路上にある頂点は必ず辺が 2 本接続されることになります。頂点が経路上になければ、辺は接続されていないので、本数は 0 になります。したがって、A, I 以外の頂点につながる辺の本数は 0 \/ 2 で表すことができます。
頂点につながっている辺の本数は、辺を表す変数の値を足し算するだけで求めることができます。プログラムは次のようになります。
リスト : 経路の探索
:- use_module(library(clpfd)).
% 辺の表示
print_edge([], []).
print_edge([1 | Xs], [Y | Ys]) :- writeln(Y), print_edge(Xs, Ys).
print_edge([0 | Xs], [_ | Ys]) :- print_edge(Xs, Ys).
% 制約
check(_, _, []).
check(S, G, [S-Xs | As]) :-
sum(Xs, #=, 1), check(S, G, As).
check(S, G, [G-Xs | As]) :-
sum(Xs, #=, 1), check(S, G, As).
check(S, G, [_-Xs | As]) :-
sum(Xs, #=, K), K in 0 \/ 2, check(S, G, As).
% 探索
solver :-
Vars = [AB, BC, AD, BE, CF, DE,
EF, DG, EH, FI, GH, HI],
Name = ["A-B", "B-C", "A-D", "B-E", "C-F", "D-E",
"E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
Vars ins 0..1,
Adj = [a-[AB, AD], b-[AB, BC, BE], c-[BC, CF],
d-[AD, DE, DG], e-[BE, DE, EF, EH], f-[CF, EF, FI],
g-[DG, GH], h-[EH, GH, HI], i-[FI, HI]],
check(a, i, Adj),
sum(Vars, #=, L),
labeling([min(L)], Vars),
writeln(Vars),
print_edge(Vars, Name).
リスト Vars に辺を表す変数を格納し、辺の名前をリスト Name に格納します。Vars の範囲は ins で 0..1 に指定します。Adj は隣接リストを表します。今回は SWI-Prolog に用意されている連想リスト (pairs) を使います。
連想リストはキー (key) と値 (value) を組にして、それをリストに格納したものです。SWI-Prolog の場合、2 つのデータをハイフン (-) でつなぐと、それを組として扱うことができます。もちろん、パターンマッチングにハイフンを使うこともできます。
簡単な使用例を示しましょう。
?- pairs_keys_values(As, [a,b,c,d], [1,2,3,4]). As = [a-1, b-2, c-3, d-4]. ?- pairs_keys_values(As, [a,b,c,d], [1,2,3,4]), pairs_values(As, X). As = [a-1, b-2, c-3, d-4], X = [1, 2, 3, 4]. ?- pairs_keys_values(As, [a,b,c,d], [1,2,3,4]), pairs_keys(As, X). As = [a-1, b-2, c-3, d-4], X = [a, b, c, d]. ?-
連想リストは述語 pairs_keys_values/3 で生成することができます。第 1 引数が連想リスト、第 2 引数がキーを格納したリスト、第 3 引数が値を格納したリストです。pairs_values/2 は連想リストから値を取り出し、pairs_keys/2 はキーを取り出します。これらの述語は、次のように pairs_keys_values を使って定義することができます。
pairs_keys(As, Xs) :- pairs_keys_values(As, Xs, _). pairs_values(As, Xs) :- pairs_keys_values(As, _, Xs).
連想リストの探索は Prolog の述語、たとえば member を使って簡単に行うことができます。次の例を見てください。
?- member(b-Y, [a-1, b-2, c-3, d-4]). Y = 2 ; false. ?- member(X-3, [a-1, b-2, c-3, d-4]). X = c ; false. ?- member(X-Y, [a-1, b-2, c-3, d-4]). X = a, Y = 1 ; X = b, Y = 2 ; X = c, Y = 3 ; X = d, Y = 4. ?-
member に組 b-Y を渡せば、キー b の値を求めることができます。同様に、X-3 を渡せば値が 3 のキーを求めることができます。キーと値を両方とも変数にすることも可能です。その場合は連想リストに格納されているキーと値を一組ずつ取り出すことができます。
次は制約の定義を説明します。制約は述語 check/3 で定義します。第 1 引数がスタート地点、第 2 引数がゴール地点、第 3 引数が連想リストです。最初の規則はスタート地点の定義です。変数のリスト Xs の総和が 1 であることを sum で指定します。2 番目の規則がゴール地点の定義で、スタートの規則と同じです。最後の規則がそれ以外の頂点の場合です。リスト Xs の総和を変数 K として、その値が 0 \/ 2 であることを指定します。
あとのプログラムは簡単です。最短経路は Vars の総和 L が最小になることなので、labeling のオプションで min(L) を指定します。これでプログラムは完成です。
それでは実行してみましょう。
?- solver. [0,0,1,0,0,0,0,1,0,0,1,1] A-D D-G G-H H-I true ; [0,0,1,0,0,1,0,0,1,0,0,1] A-D D-E E-H H-I true ; [0,0,1,0,0,1,1,0,0,1,0,0] A-D D-E E-F F-I true ; [1,0,0,1,0,0,0,0,1,0,0,1] A-B B-E E-H H-I true ; [1,0,0,1,0,0,1,0,0,1,0,0] A-B B-E E-F F-I true ; [1,1,0,0,1,0,0,0,0,1,0,0] A-B B-C C-F F-I true ; [0,0,1,0,0,0,1,1,1,1,1,0] A-D E-F D-G E-H F-I G-H true . ?-
次に示すように、最短経路は 6 通りあります。
A B C A B C A B C A--B C A--B C A--B--C | | | | | | D E F D--E F D--E--F D E F D E--F D E F | | | | | | G--H--I G H--I G H I G H--I G H I G H I
次は有向グラフを使ってプログラムを作りましょう。基本的な考え方は無向グラフと同じです。辺を変数に対応させ、選んだ辺を 1 に、選ばない辺を 0 にします。有向グラフなので、A から B へ進む辺と、B から A に進む辺を区別します。この場合、隣接リストではなく隣接行列を使うと簡単です。次のリストを見てください。
リスト : 隣接行列
% A B C D E F G H I
adjacent([[0,_,0, _,0,0, 0,0,0], % A
[_,0,_, 0,_,0, 0,0,0], % B
[0,_,0, 0,0,_, 0,0,0], % C
[_,0,0, 0,_,0, _,0,0], % D
[0,_,0, _,0,_, 0,_,0], % E
[0,0,_, 0,_,0, 0,0,_], % F
[0,0,0, _,0,0, 0,_,0], % G
[0,0,0, 0,_,0, _,0,_], % H
[0,0,0, 0,0,_, 0,_,0]]). % I
リストのリストで行列を表します。行が頂点から出て行く辺、列が頂点に入ってくる辺を表します。たとえば、A は B と D につながっているので、A 行の B 番目と D 番目の無名変数が A から出て行く辺、A 列の B 番目と D 番目の無名変数が A に入ってくる辺を表します。
制約の定義ですが、頂点に辺が接続されている場合、入力と出力はともに 1 になります。辺がつながっていない場合、入力と出力はともに 0 になります。つまり、入力 (または出力) の本数が 0 or 1 で、入力と出力の値が等しいことが条件になります。
プログラムは次のようになります。
リスト : 制約の定義
check(_, _, _, [], []).
check(S, S, G, [X | Xs], [Y | Ys]) :-
sum(X, #=, 1), % スタートから出る辺は 1 つ
sum(Y, #=, 0),
N1 #= S + 1,
check(N1, S, G, Xs, Ys).
check(G, S, G, [X | Xs], [Y | Ys]) :-
sum(X, #=, 0),
sum(Y, #=, 1), % ゴールに入る辺は 1 つ
N1 #= G + 1,
check(N1, S, G, Xs, Ys).
check(N, S, G, [X | Xs], [Y | Ys]) :-
K in 0..1,
sum(X, #=, K), % 出力と入力の辺の数が同じ (0 or 1)
sum(Y, #=, K),
N1 #= N + 1,
check(N1, S, G, Xs, Ys).
述語 check/5 の第 1 引数が行 (または列) の番号、第 2, 3 引数の S, G がスタート地点とゴール地点を表す整数、第 3 引数が隣接行列、第 4 引数がその転置行列です。つまり、X が頂点 N から出て行く辺、Y が頂点に入ってくる辺を格納したリストになります。
第 1 引数が S と等しい場合、X の総和が 1 で、Y の総和が 0 になります。G と等しい場合、X の総和が 0 で、Y の総和が 1 になります。それ以外の頂点では、X の総和と Y の総和が等しくて、その値 K は 0 or 1 になります。
最後に経路を探索する述語 solver を定義します。
リスト : 最短経路の探索
solver :-
adjacent(Xs),
flatten(Xs, Ys),
include(var, Ys, Vars), % 変数を集める
Vars ins 0..1,
transpose(Xs, Zs),
check(1, 1, 9, Xs, Zs),
sum(Vars, #=, L), % 経路長 L
labeling([min(L)], Vars),
maplist(writeln, Xs),
print_edge(0, Xs).
隣接行列 adjacent から行列を取り出して、flatten/2 で平坦化し、include/3 で自由変数を取り出します。そして、それらの範囲を ins で 0 .. 1 に設定します。次に、transpose で Xs の転置行列を生成して、述語 check で制約を設定します。あとは labeling で最短経路を探索して、print_edge で選択した辺を表示するだけです。
あとのプログラムは簡単なので説明は割愛します。詳細はプログラムリストをお読みください。
それでは実行してみましょう。
?- solver. [0,0,0,1,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,1,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,1,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,0,0] A -> D D -> G G -> H H -> I true ; [0,0,0,1,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,1,0,0,0,0] [0,0,0,0,0,0,0,1,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,0,0] A -> D D -> E E -> H H -> I true ; [0,0,0,1,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,1,0,0,0,0] [0,0,0,0,0,1,0,0,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] A -> D D -> E E -> F F -> I true ; [0,1,0,0,0,0,0,0,0] [0,0,0,0,1,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,1,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,0,0] A -> B B -> E E -> H H -> I true ; [0,1,0,0,0,0,0,0,0] [0,0,0,0,1,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,1,0,0,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] A -> B B -> E E -> F F -> I true ; [0,1,0,0,0,0,0,0,0] [0,0,1,0,0,0,0,0,0] [0,0,0,0,0,1,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] A -> B B -> C C -> F F -> I true ; [0,0,0,1,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,0,0,0] [0,0,0,0,0,0,1,0,0] [0,0,0,0,0,1,0,0,0] [0,0,0,0,0,0,0,0,1] [0,0,0,0,0,0,0,1,0] [0,0,0,0,1,0,0,0,0] [0,0,0,0,0,0,0,0,0] A -> D D -> G E -> F F -> I G -> H H -> E true . ?-
当然ですが、6 通りの最短経路が出力されます。
リスト : 最短経路の探索 (有向グラフ)
:- use_module(library(clpfd)).
% 隣接行列
% A B C D E F G H I
adjacent([[0,_,0, _,0,0, 0,0,0], % A
[_,0,_, 0,_,0, 0,0,0], % B
[0,_,0, 0,0,_, 0,0,0], % C
[_,0,0, 0,_,0, _,0,0], % D
[0,_,0, _,0,_, 0,_,0], % E
[0,0,_, 0,_,0, 0,0,_], % F
[0,0,0, _,0,0, 0,_,0], % G
[0,0,0, 0,_,0, _,0,_], % H
[0,0,0, 0,0,_, 0,_,0]]). % I
% 制約
check(_, _, _, [], []).
check(S, S, G, [X | Xs], [Y | Ys]) :-
sum(X, #=, 1), % スタートから出る辺は 1 つ
sum(Y, #=, 0),
N1 #= S + 1,
check(N1, S, G, Xs, Ys).
check(G, S, G, [X | Xs], [Y | Ys]) :-
sum(X, #=, 0),
sum(Y, #=, 1), % ゴールに入る辺は 1 つ
N1 #= G + 1,
check(N1, S, G, Xs, Ys).
check(N, S, G, [X | Xs], [Y | Ys]) :-
K in 0..1,
sum(X, #=, K), % 出力と入力の辺の数が同じ (0 or 1)
sum(Y, #=, K),
N1 #= N + 1,
check(N1, S, G, Xs, Ys).
% 辺の表示
print_edge(_, []).
print_edge(N, [X | Xs]) :-
sumlist(X, 0),
N1 is N + 1,
print_edge(N1, Xs).
print_edge(N, [X | Xs]) :-
nth0(M, X, 1),
C1 is N + 0x41,
C2 is M + 0x41,
format('~c -> ~c~n', [C1, C2]),
N1 is N + 1,
print_edge(N1, Xs).
% A から I までの最短経路を求める
solver :-
adjacent(Xs),
flatten(Xs, Ys),
include(var, Ys, Vars), % 変数を集める
Vars ins 0..1,
transpose(Xs, Zs),
check(1, 1, 9, Xs, Zs),
sum(Vars, #=, L), % 経路長 L
labeling([min(L)], Vars),
maplist(writeln, Xs),
print_edge(0, Xs).