カテゴリー: PROLOG

ECLiPSe CLPで記述するソルバーの構成はおおまかに以下の部分に分かれます。
１．制約の記述
２．解のサーチ
３．解の表示

２．の解のサーチというのは、具体的には以下の2つのループの入れ子の構造となります。
　ループ１：複数の制約変数の中から1つを選択する
　ループ２：変数にドメイン内の値を一つ割り当てる

このとき、ループ１で「どのような順番で変数を選択するか」ループ２で「どのような順番で値を決定するか」という２つの制御をおこなうことができます。この記事はループ2の値決定の順番をいろいろ試す記事となります。

ECLiPSe CLPではループ１，ループ２は例えば以下のような方法で記述可能です(他にもあると思いますが）
・labeling/1を使用
　ループ１：リストの先頭の変数から選択
　ループ２：ドメイン内の数値の小さい値から大きい値へと順に値を決定
・search/6を使用（ループ１、ループ２の制御を一括で指定）
・以下の2つの述語を使用
　ループ１：delete/5を使用して値を決める変数を決定
　ループ２：indomain/2を使用してどのように変数に値を設定するかを決定

本記事ではsearch/6に関して、「どのような順で変数に値を割り当てるか」というのをオプションを変えながら実験したいと思います。

以下のテスト用プログラムで実験します。

:-lib(ic).

search_test(Choice):-
	A #:: 1..10,
	search([A], 0, input_order, Choice, complete, []),
	write(A),write(','),
	fail.

Choice の部分を色々変更することにより実験します。解説の部分は直訳なので自分でも意味わかってない箇所あります。特に「テストされた値を削除」「ドメイン分割」のところが良くわかってません。誰か教えて下さい。

①indomain
組み込みのindomain/1を使用します。値は昇順で試行されます。失敗しても以前にテストされた値は削除されません。

?- search_test(indomain).
No (0.00s cpu)
1,2,3,4,5,6,7,8,9,10,

②indomain_min
昇順で試行されます。失敗すると以前にテストされた値が削除されます。値はindomainの場合と同じ順序でテストされますが、バックトラックがより早く発生する可能性があります。

?- search_test(indomain_min).
No (0.00s cpu)
1,2,3,4,5,6,7,8,9,10,

③indomain_max
値は降順で試行されます。 失敗すると以前にテストされた値が削除されます。

?- search_test(indomain_max).
No (0.00s cpu)
10,9,8,7,6,5,4,3,2,1,

④indomain_reverse_min
indomain_minと同様ですが、選択は逆の順序で試行されます。つまり最小の値が最初にドメインから削除され、それが失敗した場合にのみ、値が割り当てられます。

?- search_test(indomain_reverse_min).
No (0.00s cpu)
10,9,8,7,6,5,4,3,2,1,

⑤indomain_reverse_max
indomain_maxと同様ですが、選択は逆の順序で試行されます。 つまり、最大値が最初にドメインから削除され、それが失敗した場合にのみ、値が割り当てられます。

?- search_test(indomain_reverse_max).
No (0.00s cpu)
1,2,3,4,5,6,7,8,9,10,

⑥indomain_middle
値はドメインの中央から開始して試行されます。 失敗すると以前にテストされた値が削除されます。

?- search_test(indomain_middle).
No (0.00s cpu)
5,6,4,7,3,8,2,9,1,10,

⑦indomain_median
値はドメインの中央値から開始して試行されます。 失敗すると以前にテストされた値が削除されます。

?- search_test(indomain_median).
No (0.00s cpu)
6,7,5,8,4,9,3,10,2,1,

⑧indomain_split
ドメインの下半分を最初に試行して、連続するドメイン分割によって試行されます。 失敗すると、試行された間隔が削除されます。 これは、indomainまたはindomain_minと同じ順序で値を列挙しますが、それより早く失敗する可能性があります（注：解サーチ時はなるべく早く失敗したほうが処理時間が速くなる。早く失敗≒変数の決定の入れ子のループのより外側のループ回数が減ることと同等）

?- search_test(indomain_split).
No (0.00s cpu)
1,2,3,4,5,6,7,8,9,10,

⑨indomain_reverse_split
ドメインの上半分を最初に試行して、連続するドメイン分割によって試行されます。 失敗すると、試行された間隔が削除されます。 これは、indomain_maxと同じ順序で値を列挙しますが、それより早く失敗する可能性があります。

?- search_test(indomain_reverse_split).
No (0.00s cpu)
10,9,8,7,6,5,4,3,2,1,

⑩indomain_random
値はランダムな順序で試行されます。 バックトラック時に、以前に試行された値が削除されます。 このルーチンを使用すると、別の呼び出しで異なる順序で乱数が作成されるため、再現性のない結果が生じる可能性があります。 この方法では、組み込みのrandom / 1を使用して乱数を作成し、seed / 1を使用して別の実行で同じ番号生成シーケンスを強制することができます。

?- search_test(indomain_random).
No (0.00s cpu)

何度も実行（その都度結果が異なる）
8,2,1,7,3,9,10,4,5,6,
4,6,3,10,5,2,1,7,9,8,
7,10,9,4,1,2,3,8,6,5,
10,4,7,6,3,8,5,9,2,1,
10,3,5,7,8,4,1,2,9,6,

⑪indomain_interval
ドメインが複数の間隔で構成されている場合は、最初に間隔の選択で分岐します。 1つの間隔では、ドメイン分割を使用します。

?- search_test(indomain_interval).
No (0.00s cpu)
1,2,3,4,5,6,7,8,9,10,

ECLiPSe CLPのEPLEXライブラリを使用して混合整数線形計画問題を解きます。

混合整数計画問題とは、線形計画問題の変数のすべてもしくは一部が整数という条件がついた問題です。

Problems and exercises in Operations Research の 3.2 Gomory cutsの問題を解きます。

ほぼ線形計画問題と同じですが、x(x1とx2)に整数であるという条件がついています。

線形計画法の時と比べたプログラムの書き方の違いなのですがなんと整数の制約を持つ変数にintegers制約を追加するだけでできてしまいます。（下記のプログラムの13行目）
線形計画問題と混合整数計画問題の解き方の違いとして、混合整数計画問題のほうは分枝限定法を使用したり複雑な方法になるのですが、それらはライブラリの中で行ってくれるため使用する側はほとんど手間が変わりません。

プログラム：

:- lib(eplex).
:- eplex_instance(prob). % declare an eplex instance

solve(Cost, Vars) :-
	% create the problem variables and set their range
	Vars = [X1,X2],
	prob: (Vars $:: 0.0..1.0Inf),
	% post the constraints for the problem to the eplex instance
	prob: (-4*X1 + 6*X2 $=< 9),
	prob: (X1 + X2 $=< 4),
	prob: (X1 $>= 0),
	prob: (X2 $>= 0),
	prob: integers(Vars),
	
	% set up the external solver with the objective function
	prob: eplex_solver_setup(min(X1 - 2*X2)),
	prob: eplex_solve(Cost). % Solve problem using external solver

実行結果：

?- solve(Cost, [X1, X2]).
Cost = -3.0
X1 = X1{0.0 .. 1.7976931348623157e+308 @ 1.0}
X2 = X2{0.0 .. 1.7976931348623157e+308 @ 2.0}
Yes (0.00s cpu)

X1 = 1
X2 = 2
が解答です。かなり便利です！

ECLiPSe CLPのEPLEXライブラリを用いて線形計画法の問題を解く方法を説明します。

EPLEXは外部の数理計画法(Mathematical Programming)ソルバーで、線形計画法(Linear Programming = LP)や混合整数計画法(Mixed Integer Programming = MIP)の問題を解くことができます。

線形計画法とは複数の変数に1次方程式・1次不等式で記述されたいくつかの制約があるときに、目的関数(これも一次式)を最大化もしくは最小化する解を得るためのアルゴリズムです。

混合整数計画法は線形計画法の変数の一部もしくは全部が「整数である」という制約を持つ問題です（こっちのほうが難しいようです）

EPLEXライブラリに関してはTutorial の Chapter 16 The Eplex Libraryに使用方法が書いてあります。EPLEXはECLiPSe CLPのネイティブなライブラリではなく、外部（サードパーティ）のソルバーです。Eplexライブラリを用いることにより他のネイティブなライブラリと同じように使用可能です。また、ECLiPSe CLPをインストールした際に自動的にインストールされるので、特に設定の必要はなく初めから使えます。

世に線形計画法のソルバーはたくさんあるので、別にECLiPSe CLPを使用しなくてもよいのですが、ECLiPSeを使えば他の便利な制約論理プログラミングの制約が同時に使用できるので、のちのち良いことがあると思います。

問題は Problems and exercises in Operations Research の Chapter 2 Linear programming 2.1 Graphical solutionを解いてみます。

問題①

Tは転置行列です。
cx = 16*X1 + 25*X2 を最小にするような X1, X2の値を求めよという問題です。

プログラム：

:- lib(eplex).
:- eplex_instance(prob). % declare an eplex instance

solve(Cost, Vars) :-
	% create the problem variables and set their range
	Vars = [X1,X2],
	prob: (Vars $:: 0.0..1.0Inf),
	
	% post the constraints for the problem to the eplex instance
	prob: (1*X1+7*X2 $>= 4),
	prob: (1*X1+5*X2 $>= 5),
	prob: (2*X1+3*X2 $>= 9),
	prob: (X1 $>= 0),
	prob: (X2 $>= 0),
	
	% set up the external solver with the objective function
	prob: eplex_solver_setup(min(16*X1+25*X2)),
	prob: eplex_solve(Cost). % Solve problem using external solver

実行結果：

?- solve(Cost, [X1, X2]).
Cost = 72.142857142857139
X1 = X1{0.0 .. 1.7976931348623157e+308 @ 4.2857142857142856}
X2 = X2{0.0 .. 1.7976931348623157e+308 @ 0.14285714285714293}
Yes (0.00s cpu)

解説：
ほとんどそのまま制約を書いてあるだけですが、ポイントは
①2行目、eplex_instance()でprobという名前のインスタンスを設定し、EPLEXの機能を使うときはそのインスタンスを指定する。
②17行目、eplex_solver_setupで目的関数の設定を行う
③実行結果のX1{0.0 .. 1.7976931348623157e+308 @ 4.2857142857142856}の部分は　変数名{最小値..最大値 @ 解}　の意味
くらいです。あとは見ればわかると思います。

問題②：

プログラム：

:- lib(eplex).
:- eplex_instance(prob). % declare an eplex instance

solve(Cost, Vars) :-
	% create the problem variables and set their range
	Vars = [X1,X2],
	prob: (Vars $:: 0.0..1.0Inf),
	
	% post the constraints for the problem to the eplex instance
	prob: (	2*X1 + X2 $=< 4),
	prob: (-2*X1 + X2 $=< 2),
	prob: (	  X1 - X2 $=< 1),
	prob: (X1 $>= 0),
	prob: (X2 $>= 0),
	
	% set up the external solver with the objective function
	prob: eplex_solver_setup(max(3*X1+2*X2)),
	prob: eplex_solve(Cost). % Solve problem using external solver

実行結果：

?- solve(Cost, [X1, X2]).
Cost = 7.5
X1 = X1{0.0 .. 1.7976931348623157e+308 @ 0.5}
X2 = X2{0.0 .. 1.7976931348623157e+308 @ 3.0}
Yes (0.00s cpu)

問題③

プログラム：

:- lib(eplex).
:- eplex_instance(prob). % declare an eplex instance

solve(Cost, Vars) :-
	% create the problem variables and set their range
	Vars = [X1,X2,X3],
	prob: (Vars $:: 0.0..1.0Inf),
	
	% post the constraints for the problem to the eplex instance
	prob: (2*X1 + 3*X3 $=1),
	prob: (3*X1+2*X2 - X3 $= 5),
	prob: (X1 $>= 0),
	prob: (X2 $>= 0),
	prob: (X3 $>= 0),
	
	% set up the external solver with the objective function
	prob: eplex_solver_setup(min(X1 - 2*X2)),
	prob: eplex_solve(Cost). % Solve problem using external solver

実行結果：

?- solve(Cost, [X1, X2, X3]).
Cost = -5.333333333333333
X1 = X1{0.0 .. 1.7976931348623157e+308 @ 0.0}
X2 = X2{0.0 .. 1.7976931348623157e+308 @ 2.6666666666666665}
X3 = X3{0.0 .. 1.7976931348623157e+308 @ 0.33333333333333331}
Yes (0.00s cpu)

カンタンですね！

ECLiPSe CLPで機械の買い換え問題（Problems and exercises in Operations Researchの1.5 Renewal plan）を解くプログラムの説明をします。

問題：
ある会社が12000ユーロの機械を購入し使用しています。機械には毎年維持費(costs)がかかります（維持費は年々上がっていきます）。
使用中の機械は中古で販売(gain)し新しく買いなおすこともできます（中古機械の販売価格は年々下がっていきます）
維持費と中古機械の販売価格は以下の表のとおりになっています（1keuroは1000ユーロ）。
5年間の総運用コストを最小限に抑える機械の買い換え計画を決定しなさい。

行間を読むとこの問題は機械の生産する利益に関しては触れられていないのでこれは毎年一定で考慮する必要なしということのようです。

難しくて解けなかったので答えと解説を見ました。説明します。
「i年目に購入した機械をj年目まで使用し、買い替えたときのコスト」をCijという変数で表し、1～6年目時点を表す6つのノードがCijで結ばれた有向グラフを考える（6は5年間の期限終了時を表す）
Cijの計算方法は以下となります。

Cij = 12000 + [(j-i)年間にかかった維持費] – [(j-i)年使用した機械の中古販売価格]

たとえば　2年目に購入した機械を4年目まで使用し買い替えたコストだと
C24 = 12000 + (2000+4000) – 6000 = 12000
となります。これはノード2とノード4を結ぶエッジとなります。
そのように計算した有向グラフは以下のようになります（回答からコピー）

このようなグラフの、ノード１からノード６へ至る最小コストの経路が答えとなります。

single_pair_short_path/6を使用すると、最小コストの経路をすべて取得できます。

single_pair_short_path(+Graph, +DistanceArg, +SourceNode, +SinkNode, +Tolerance, -Path)
　+Tolerance を 0にするとすべての最小コストの経路が取得され、当然それらのコストはみな同じとなります。
　+Tolerance に0より大きい数値を設定すると、最小＋その値　までのコストの経路も抽出されるようです。

プログラム

:-lib(graph_algorithms).

solve(Path):-
	make_graph(6,
		[ 
			e(1,2,7),
			e(1,3,12),
			e(1,4,21),
			e(1,5,31),
			e(1,6,44),
			e(2,3,7),
			e(2,4,12),
			e(2,5,21),
			e(2,6,31),
			e(3,4,7),
			e(3,5,12),
			e(3,6,21),
			e(4,5,7),
			e(4,6,12),
			e(5,6,7)
		],Graph),
	single_pair_short_path(Graph, 0, 1, 6, 0, Path).

実行結果：

?- solve(Path).
Path = 31 - [e(5, 6, 7), e(3, 5, 12), e(1, 3, 12)]
Yes (0.00s cpu, solution 1, maybe more)
Path = 31 - [e(4, 6, 12), e(3, 4, 7), e(1, 3, 12)]
Yes (0.00s cpu, solution 2, maybe more)
Path = 31 - [e(4, 6, 12), e(2, 4, 12), e(1, 2, 7)]
Yes (0.00s cpu, solution 3)

コストは31(=31000ユーロ)が最小で、以下の買い換え計画があるようです。
3年目と5年目で買い替える
3年目と4年目で買い替える
2年目と4年目で買い替える

ECLiPSe CLPでグラフの最小カットを求める方法を説明します。

グラフのカットというのは、最大流問題で用いた、エッジに許容流量の設定のあるスタートとゴールのノードを持つ有向グラフのすべてのノードを、「スタート側」「ゴール側」の2種類のノードに分けることをいいます。

ノードをどのように分けてもよいので、カットの仕方には複数通りの方法があります。
グラフに|V|個のノードがあったとして、スタートとゴールのノードのみどちら側か決まっているので　2^(|V|-2)通り　のカットの方法があります。

最小カットというのは、複数あるカットのうち、「カット容量が最小となるもの(そしてその容量)」のことを言います。
カット容量とは、ノードとノードを結ぶすべてのエッジのうち「スタート側に属するノード」から「ゴール側に属するノード」の方向のエッジの許容流量の総計のことを言います（スタート側→スタート側、ゴール側→ゴール側、ゴール側→スタート側 のエッジの容量は数えない）

最小カット容量は、スタート→ゴールの流れのボトルネックのようなもので、この値がそのまま最大流と等しくなるとのことです。
（最小カット容量＝最大フロー）

前置き長くなりましたが上記最小カットを求めます。
library(all_min_cuts) の述語 all_min_cuts/8 を使います。

all_min_cuts(+Graph, +CapacityArg, +SourceNode, +SinkNode, -MaxFlowValue, -MaxFlowEdges, -MinCuts, -MinCutEdges)

Grapha – グラフ
CapacityArg – エッジ e(_,_,A) のAのどこが最大流量か設定する
SourceNode – スタートのノード番号
SinkNode – ゴールのノード番号
MaxFlowValue – 最大フロー(=最小カット容量)(戻り値)
MaxFlowEdges – 「フロー　- エッジ」のリスト(戻り値)
MinCuts – 最小カット容量をもつカットのリスト（複数あれば複数）(戻り値)
MinCutEdges – すべての最小カット組（それぞれのカット組はスタート側-ゴール側を分けるエッジのリスト）(戻り値)

①以下の図でノード１がスタート、ノード7がゴールとして最小カットを求めます。

プログラム：

:-lib(graph_algorithms).
:-lib(all_min_cuts).

solve(MaxFlowValue, MaxFlowEdges, MinCuts, MinCutEdges):-
	make_graph(
		7,
		[ 
			e(1,2,6),
			e(1,3,6),
			e(2,3,1),
			e(2,4,3),
			e(2,5,3),
			e(3,5,7),
			e(4,5,1),
			e(4,7,1),
			e(5,7,4),
			e(5,6,5),
			e(6,7,4)
		],
		Graph),
	all_min_cuts(Graph, 0, 1, 7, MaxFlowValue, MaxFlowEdges, MinCuts, MinCutEdges).

実行結果（見やすいように手で改行入れてます）：

?- solve(MaxFlowValue, MaxFlowEdges, MinCuts, MinCutEdges).
MaxFlowValue = 9
MaxFlowEdges = 
[
4 - e(6, 7, 4), 
4 - e(5, 7, 4), 
4 - e(5, 6, 5), 
1 - e(4, 7, 1), 
1 - e(4, 5, 1), 
4 - e(3, 5, 7), 
3 - e(2, 5, 3), 
2 - e(2, 4, 3), 
4 - e(1, 3, 6), 
5 - e(1, 2, 6)
]
MinCuts = [[2, 5, 3, 6, 1, 4]]
MinCutEdges = 
[
[
e(4, 7, 1), 
e(6, 7, 4), 
e(5, 7, 4)
]
]
Yes (0.00s cpu)

最小カットは９で
ノード　4-7、6-7、5-7　のエッジがスタート側とゴール側を分けるエッジのようです。

②以下のグラフで最小カットを求めます。

プログラム（便宜的にスタートノードとゴールノードをそれぞれノード5、ノード6に置き換え）：

:-lib(graph_algorithms).
:-lib(all_min_cuts).

solve(MaxFlowValue, MaxFlowEdges, MinCuts, MinCutEdges):-
	make_graph(
		6,
		[ 
			e(5,1,3),
			e(5,2,10),
			e(1,4,4),
			e(1,3,4),
			e(2,1,7),
			e(2,3,2),
			e(4,2,3),
			e(4,3,8),
			e(4,6,3),
			e(3,6,8)
		],
		Graph),
	all_min_cuts(Graph, 0, 5, 6, MaxFlowValue, MaxFlowEdges, MinCuts, MinCutEdges).

実行結果：

?- solve(MaxFlowValue, MaxFlowEdges, MinCuts, MinCutEdges).
MaxFlowValue = 10
MaxFlowEdges = 
[
	7 - e(5, 2, 10), 
	3 - e(5, 1, 3), 
	3 - e(4, 6, 3), 
	1 - e(4, 3, 8), 
	7 - e(3, 6, 8), 
	2 - e(2, 3, 2), 
	5 - e(2, 1, 7), 
	4 - e(1, 4, 4), 
	4 - e(1, 3, 4)
]
MinCuts = [[2, 5, 1]]
MinCutEdges = 
[
	[e(1, 4, 4), e(1, 3, 4), e(2, 3, 2)]
]
Yes (0.02s cpu)

最大フロー(最小カット)10、
ノード1-4、1-3、2-3　を結ぶエッジがスタート側とゴール側を結ぶ

ECLiPSe CLPで最大流問題を解く方法の説明をします。
最大流問題って自分は知らなかったのですが、グラフが与えられ、ノードからノードを結ぶそれぞれのエッジに転送量の上限（転送するのはデータでも物でも何でも良い）が決められているとき、あるノードからもう一つのノードへの最大の転送量を求めるという問題です。

電子回路のキルヒホッフの法則のように、各ノードに入る量とそのノードから出ていく量は一致している必要があります（始点と終点は除く）つまりあるノードに注目したとき、入ってくる量の許容量が6で出ていく量の許容量が4の場合は少ない方の4に合わされてしまうことになります（ボトルネックとなる）

自分の勝手な解釈ですが、いくつかのノードが接続された形の水道管ネットワークみたいのがあったとして、ノード同士をつなぐ水道管の最大流量が決まっているとき、入口Aから出口Bに流せる最大の水の量を求めるみたいな感じかなとイメージしました。

以下のようなグラフ（Problems and exercises in Operations Research 10ページ目 1.3 Maximum flow の図）が与えられたときにノード1からノード7への最大流を求めるプログラムを書きます。
エッジの数値はノード間の距離ではなくエッジの許容流量を表していることに注意です。

library(max_flow)を使用します。
またlibrary(graph_algorithms)も使用します。

グラフの作成はmake_graph/3を使用します（使い方は以前の記事参照）

最大流問題を解く述語はmax_flow/5やmax_flow/7ですが、今回は得られる情報が多いmax_flow/7を使用します。

max_flow(
　　+Graph,
　　+CapacityArg,
　　+SourceNode,
　　+SinkNode,
　　-MaxFlowValue,
　　-MaxFlowEdges,
　　-MaxFlowEdgesGraph)

引数の意味は以下となります。
　1.Graph　生成したグラフ
　2.CapacityArg エッジの許容量をどこに書いたか指定
　　（e(_,_,Capacity)のCapacityの部分に数値で記載した場合は0)
　3.SourceNode 開始ノード番号
　4.SinkNode 終了ノード番号
　5.MaxFlowValue 最大流量（戻り値）
　6.MaxFlowEdges 使用されたエッジとそのエッジに流される流量
　　（流量ゼロのエッジは出力されない）　（戻り値）
　7.MaxFlowEdgesGraph 使用されたエッジをグラフの形で出力
　　（流量ゼロのエッジは出力されない）　（戻り値）

プログラム：

:-lib(graph_algorithms).
:-lib(max_flow).

solve(MaxFlowValue,MaxFlowEdges,MaxFlowEdgesGraph):-
	make_graph(
		7,
		[ 
			e(1,2,6),
			e(1,3,6),
			e(2,3,1),
			e(2,4,3),
			e(2,5,3),
			e(3,5,7),
			e(4,5,1),
			e(4,7,1),
			e(5,7,4),
			e(5,6,5),
			e(6,7,4)
		],
		Graph),
	max_flow(Graph, 0, 1, 7, MaxFlowValue,MaxFlowEdges,MaxFlowEdgesGraph).

実行結果（見やすくなるよう手で改行入れてます）：

?- solve(MaxFlow, Edges, Graph).
MaxFlow = 9
Edges = 
[
4 - e(6, 7, 4), 
4 - e(5, 7, 4), 
4 - e(5, 6, 5), 
1 - e(4, 7, 1), 
1 - e(4, 5, 1), 
4 - e(3, 5, 7), 
3 - e(2, 5, 3), 
2 - e(2, 4, 3), 
4 - e(1, 3, 6), 
5 - e(1, 2, 6)
]
Graph = 
graph(7, 
[](
[e(1, 2, 6), e(1, 3, 6)], 
[e(2, 4, 3), e(2, 5, 3)], 
[e(3, 5, 7)], 
[e(4, 5, 1), e(4, 7, 1)], 
[e(5, 6, 5), e(5, 7, 4)], 
[e(6, 7, 4)], []),
 _2193, _2194, _2195, _2196, _2197, _2198)
Yes (0.00s cpu)

最大9の量流せるようですね。

一応注意点ですがこのlibrary(max_flow)はステータスが「prototype」となっていることに気を付けてください。もしバグ見つかったらECLiPSeの開発者に報告すると喜ばれます。手間でしたら私（koyahataアットマークkoyahatataku.com)に報告くだされば私が報告します。

以上

ECLiPSe CLPを使ってグラフの最短経路木(shortest path tree)を求める方法を説明します。

前回の例ではshortest_paths/4を使用して最短経路木を求めましたが、今回はshortest_paths_bellman_ford/4を使用しBellman-Ford’sアルゴリズムで求めます。

特徴としては、前回の例ではコストは負の数が許されなかったのですが、Bellman-Ford’sアルゴリズムでは負のコストも許されます。マニュアルには記載されていませんが前回のshortest_paths/4はおそらくダイクストラ法を使っているのではと予想しています（ソースにもコメントなし）。

以下のようなグラフ（Problems and exercises in Operations Research 9ページ目 1.2 Bellman-Ford’s algorithm の図）で最短経路木を求めます。

shortest_paths_bellman_ford/4 の引数はコストに負数が許されること以外はshortest_paths/4と同じなので説明を割愛します（前回記事を参照）

プログラム：

:-lib(graph_algorithms).

solve(Path):-
	make_graph(
		6,
		[ 
			e(1,2,2),
			e(2,4,-1),
			e(1,3,4),
			e(3,5,2),
			e(5,2,-1),
			e(5,6,5),
			e(4,5,0),
			e(6,4,3),
			e(3,6,1)
		],
		Graph),
	shortest_paths(Graph,0,1,Path).

実行結果（見やすくするために改行を手で入れてます）：

?- solve(Path).
Path = 
[](0 - [], 
2 - [e(1, 2, 2)], 
4 - [e(1, 3, 4)], 
1 - [e(2, 4, -1), e(1, 2, 2)], 
1 - [e(4, 5, 0), e(2, 4, -1), e(1, 2, 2)], 
5 - [e(3, 6, 1), e(1, 3, 4)])
Yes (0.00s cpu)

出力結果の見方は前回記事と同じ

総コストが負となるループがある場合（この例ではノード2->4->5) 何度もそのループを回ることでいくらでもコストを下げられるので「一番コストが低い経路」というのは存在しないのですが、ここらへんの説明はドキュメントにはありませんでした。常識的に考えるとノードを重複して通らない前提で一番コストが低い経路を返しているとは思いますが…

本記事で説明した「あるスタート地点からすべてのノードへの経路（1対N）」ではなく「あるスタート地点からあるゴールまでの経路（1対1）」を求めるsingle_pair_という接頭語が付いた述語もありますので興味のある方は調べてみてください。

ECLiPSe CLPでグラフの最短経路木を求めるやり方の紹介をします。

graph_algorithmsライブラリが使用できます。

まずmake_graph/3でグラフを作成します。引数の意味は、
第一引数：ノードの個数
第二引数：エッジ群の情報
第三引数：生成されたグラフ（戻り値）

例えば以下のグラフ（Problems and exercises in Operations Research 9ページ目 1.1 Dijkstra’s algorithm の図）の場合：

	make_graph(
		6,
		[ 
			e(1,2,2),
			e(2,4,3),
			e(1,3,4),
			e(3,5,2),
			e(5,2,8),
			e(5,6,5),
			e(4,5,0),
			e(6,4,3),
			e(3,6,1)
		],
		Graph)

みたいな感じで定義します。e のとこの意味は
e(ノード１番号,ノード２番号,コスト）
みたいな感じです（ノード１→ノード２の向きが存在します）
このコストのところはもっと複雑な情報を設定できるようです。
無向グラフの場合は逆向きエッジの定義も必要となります。

最短経路木なのですが、shortest_paths/4を使用して得ることが出来ます。
引数の意味は
第一引数：先ほど作ったグラフ
第二引数：e(_,_,EdgeData)で定義したエッジの3番目の引数EdgeDataのどの引数を距離（コスト）として使用するか（今回は0)
第三引数：スタート地点のノード番号
第四引数：生成された経路情報

第二引数の説明は以下を読んでください
DistanceArg refers to the graph's EdgeData information that was specified when the graph was constructed. If EdgeData is a simple number, then DistanceArg should be 0 and EdgeData will be taken as the length of the edge. If EdgeData is a compound data structure, DistanceArg should be a number between 1 and the arity of that structure and determines which argument of the EdgeData structure will be interpreted as the edge's length. Important: the distance information in EdgeData must be a non-negative number, and the numeric type (integer, float, etc) must be the same in all edges.

If DistanceArg is given as -1, then any EdgeData is ignored and the length of every edge is assumed to be equal to 1.

意訳：
DistanceArgには
・e(_,_,EdgeData)のEdgeDataが1つの数値で表されている場合は0を指定
・すべてのエッジの距離が1固定の場合は-1を指定
・「1～EdgeDataのアリティ」の範囲内で指定した場合は復号項(compound term)の何番目に距離が設定されてるかを指定したことになる
距離（コスト）は非負の数値のみ許可される。

今回の定義でノード１からすべてのノードへの最短経路を取得する場合は、以下のように呼びます

shortest_paths(Graph,0,1,Path)

プログラム全体と結果は以下のようになります。

プログラム：

:-lib(graph_algorithms).

solve(Path):-
	make_graph(
		6,
		[ 
			e(1,2,2),
			e(2,4,3),
			e(1,3,4),
			e(3,5,2),
			e(5,2,8),
			e(5,6,5),
			e(4,5,0),
			e(6,4,3),
			e(3,6,1)
		],
		Graph),
	shortest_paths(Graph,0,1,Path).

結果（見づらいのでPathのとこ改行しました。本当は1行でダダっと出ます）：

?- solve(Path).
Path = 
[](
0 - [], 
2 - [e(1, 2, 2)], 
4 - [e(1, 3, 4)], 
5 - [e(2, 4, 3), e(1, 2, 2)], 
5 - [e(4, 5, 0), e(2, 4, 3), e(1, 2, 2)], 
5 - [e(3, 6, 1), e(1, 3, 4)]
)
Yes (0.00s cpu)

結果のPathの見方なのですが、
最初の 0 – [] （おそらく自分自身へ行く経路）のあと、
ノード2へ行く最短パス、ノード3へ行く最短パス、、、ノード6へ行く最短パス
がそれぞれ表示されてます。
2 – [～]、4 – [～]の初めの数字はかかるコストが出てます。
続くリストは、どのエッジを通ったかをリストアップしてるのですが、目的地（ゴール）～ノード１（スタート）の順で出てます

以上

ECLiPSe CLPのexampleページにあるnonogram solverを解析していて、大体理解したので説明します。（内部で regular expression constraintを使用しています。）
nonogramは日本でいうところのお絵かきロジック、ピクロスです。

ルールの説明は割愛します。

プログラムの説明も細かいところ端折ってエッセンスの説明だけします。
このプログラムの一番重要な部分は、それぞれの行列のヒントのリスト(例えば[3,3,2]など)から盤面を表すリストの制約を生成するline述語で、この内容がわかれば8割方理解したことになります。

line述語にはヒントのリストと盤面1行分(or1列分)のリストを渡します。
line述語内でfromtoを駆使してヒントからTableに怪しげなリストを生成しています。

以下の例では、tkeclipseのコンソール？で
length(Xs,15),line([3,2,3],Xs).
と入力してステップ実行した結果です。これは、ピクロスの問題でいうところの、ヒントとして「3,2,3」が与えられた15マスの1つの行の制約をつける処理に相当します。

このfromtoの説明も今回は割愛します。このプログラムではあまり良い使い方してないので他の例で学んだほうが良いです。fromtoなどのLogical Loopはかなり重要なのでいつか改めて解説記事書くと思います。

このTableなのですが、fromtoの操作が終了した時点では図のような内容になっています（右の四角内）

並べ替えて書くと以下の内容です（ここら辺ぐちゃぐちゃに登録してるのがこのプログラムの良くない≒テキトーなところなので見習わなくて良いです）
[0,1,1]
[1,1,2]
[1,2,3]
[1,3,4]
[0,4,5]
[0,5,5]
[1,5,6]
[1,6,7]
[0,7,8]
[0,8,8]
[1,8,9]
[1,9,10]
[1,10,11]
[0,11,11]

このリストなのですが、オートマトンの状態遷移表を意味していて、盤面のリストを左から読み込んでいった際の状態の変遷を記述してます。
3つある数字の意味は
[読み込んだ数字(0 or 1), 現在の状態, 次の状態]
を表していて、
「現在の状態　で　0または1　を読み込んだら　次はどの状態になるか」
をすべて書き出しています。

状態を11個持つ以下のようなオートマトンが出来たことになります。便宜的に状態にqをつけました。

状態はq1から開始し、矢印に書いてある数字を読み込んだ際に次の状態に遷移します。
矢印がない数字が入力された際はヒントの条件を破ってしまっているのでオートマトン不受理（エラー、制約違反）となります。
q11の状態が最終状態（受理≒読み込み終わったときにこの状態であればOK）となります。

状態遷移表が完成したのち、それをregularという術語に渡して実際の制約を作成しています。

regularの説明です。
Qsは状態のリストです。この例では15個数字を読み込むので、それぞれの読み込み時に対応したものとなります。作成したオートマトンから、上記の例だと最初の状態がq1、最後Qnがq11となります。
与えられた盤面の1行(列)のリストに「現在の状態」「次の状態」の2つの変数を追加して
[マス目の数字, 現在の状態, 次の状態]
が1要素となるリスト（この例だと要素は15個）を作成しています。
状態の部分にはQsの要素を順番に指定して以下のように隣の要素同士を関連付けています。

そして、それぞれの要素に関し
(Goal infers ac)@Module

で実際の制約をつけるのですが、これが実際どのようにコールされているかというと、
member([マス目の数字, 現在の状態, 次の状態], List) infers ac
という呼ばれ方です。Listの中身はオートマトンの状態遷移表です。
この操作で、[マス目の数字, 現在の状態, 次の状態]の3つの組が、オートマトンの状態遷移表の要素のいずれかの組み合わせに一致する（それ以外の組はありえない）という制約を設定しています。

member(A,L)はご存知の通り「AはL内の要素の一つ」を表すprologの組み込み述語で、ECLiPSeの制約ではありません。
ここで使用されている 「infers」というのがこのブログでも何度か触れたpropiaという仕組みで、任意のprologの述語（組み込み述語でも自作述語でも良い）をECLiPSeの制約に変換してしまうというすごい強力な仕組みです。memberは本来バックトラックを発生させてしまうため制約としては使用できないのですが、propiaを使用することにより「AはL内の要素」という意味の制約として機能するようになります。
infers ac の acの部分は制約伝搬の際にどの程度絞り込むか（候補の枝狩り）の動作に関連します。
acはarc consistency の略でアーク整合のことで、これは制約プログラミングの用語です。自分はぼんやりとしか把握してないのですがそのうちちゃんと理解したら記事書こうと思います。
infersのマニュアルページ

ちなみにECLiPSeは主にアーク整合アルゴリズムAC-3とAC-5を使用しているらしいです。

制約を作るときにオートマトンを使用するというのは面白いアイデアで今後何らかの参考になるかもしれません。

Nonogram Solverは今回解説したオートマトン使用したものの他にECLiPSeメイン開発者のJoachim Schimpfさんが作成したgecodeバージョンがあり、こちらのほうが記述が簡潔なのでそのうち解析してみたいです。