ECLiPSe CLPで以下の問題を解くプログラムを書く方法を説明します。

問題（Problems and exercises in Operations Research 19ページ 4.1 Compact storage of similar sequences）

DNAマッピングプロセスの際には、非常に長い同じ長さ＆ほとんど同じ並びのDNA配列をコンパクトに保存するという問題に良く遭遇する。
ここでは問題を単純化するためDNA配列を0と1の配列と仮定する。2つのビット列をそれぞれA[i]、B[i]の配列として表現したとき、これらのハミング距離は「|A[i]-B[i]|のビット列全体での総和」で定義される（= 同じ位置のビットが異なる箇所の総数）
以下の6つのビット列のそれぞれのハミング距離は以下となる（問題からコピー）

ハミング距離が大きすぎなければ、1つの完全なビット列のみストレージに記憶し、残りはそのビット列からの派生として差分のみ保存するという方法を採用できる。
上記を実現するため、全ビット列を表現するために最小の差分となるような全域木を求めよ。

解説：
いきなりDNAとか出てきてびっくりしてしまいますが、少し考えると全ノード（ビット列）が接続されたグラフがあり、ノードを結ぶそれぞれのエッジの距離がハミング距離であるグラフを考え、その最小全域木を求める問題ということがわかります。
最小全域木を求めるためにはgraph_algorithmsライブラリ の
minimum_spanning_tree/4　という術語を使用します。

minimum_spanning_tree(+Graph, +DistanceArg, -Tree, -TreeWeight)　の引数の説明
Graph：グラフ
DistanceArg：エッジ e(_,_,Distance) の Distanceの部分に数字で距離を設定する場合は0を指定
Tree：求められた最小全域木(エッジのリスト)。戻り値
TreeWeight：求められた最小全域木の総コスト。戻り値

:-lib(graph_algorithms).

solve(Tree,TreeWeight):-
	make_graph(
		6,
		[ 
			e(1,2,4),
			e(1,3,4),
			e(1,4,5),
			e(1,5,4),
			e(1,6,3),
			e(2,3,4),
			e(2,4,3),
			e(2,5,4),
			e(2,6,5),
			e(3,4,5),
			e(3,5,2),
			e(3,6,5),
			e(4,5,3),
			e(4,6,6),
			e(5,6,5)
		],
		Graph),
	minimum_spanning_tree(Graph, 0, Tree, TreeWeight).

実行結果

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

以下のような全域木となることがわかります（解答からコピー）

ECLiPSe CLPを用いて以下のナップサック問題を解きます。
Example の Knapsack のページ を参考にしてプログラムを作成しました。

問題（Problems and exercises in Operations Research の 18ページ 3.5 Knapsack Branch and Bound ）：

投資銀行の総予算は1400万ユーロで、4種類の投資が可能です。次の表は、投資額(Amount)と各投資による純収入(Net revenue)を示しています（単位は百万ユーロ）。投資を行う場合は、各投資を全額行う必要があります。総純収入を最大化する投資方法を答えよ。

この問題文だけでは読み取れなかったのですが、解答を見ると4種類の投資は「投資しない/１個口のみ投資する」の2択らしく、同じ投資を2個以上するのは禁止のようです。

解き方：
「制約を満たす」だけではなく「利益を最大化する」「コストを最小化する」解を求める際は分枝限定法(branch_and_bound)ライブラリを使用します。

分枝限定法の厳密な定義やアルゴリズムに関しては他のページや書籍を参考にしてください。
ECLiPSe CLPのbranch and bound ライブラリでは以下の手順でコストを最小化する解を得ます。
１．コスト以外の制約を満たす解を得ます。
２．その時のコストを計算します。計算したコスト＝Aとします。
３．新たな制約「コストがA未満であること」を追加して、解を得ます。
４．上記操作を繰り返して解が得られなくなったら、直前に計算したコストが最小コストの解となります。

プログラム：

:- lib(ic).
:- lib(branch_and_bound).

model(Capacity, Values, Amounts, Revenues, SumRevenues) :-
	( foreach(Amount, Amounts), foreach(Value, Values),  param(Capacity) do
		Value #::0..1   % ある投資は最大でも1回しか行えない
	),
	integers(Amounts),
	sum(Amounts*Values) $= SumAmounts,
	SumAmounts $=< Capacity,
	sum(Values*Revenues) $= SumRevenues.

solve(Capacity, Values, Amounts, Revenue, SumRevenues) :-
	model(Capacity, Values, Amounts, Revenue, SumRevenues),
	Cost #= -SumRevenues,
	minimize(search(Values, 0, largest, indomain_reverse_split, complete, []), Cost).

解説：
6行目でValuesに[投資1,投資2,投資3,投資4]のそれぞれの投資数(0 or 1)の制約を設定しています。Valuesの個々の要素が確保されているのもこの部分ですが、動作に関してfor系の説明をしなくてはなりません。いつか改めて記事書きます。

13行目の述語solve内にある16行目の minimize/2 で分枝限定法を行っていて、繰り返し実行するgoal（この例ではsearch/6)と、コスト計算結果が入る変数（この例ではCost）を指定すると、先ほど解説した1～4の一連の動作を自動で行ってくれます。minimizeではCostを最小化する解を得ます（問題文が利益最大化を目指すものなので15行目でコスト変数への代入時に符号を逆転してます）

16行目のsearch/6の引数のうち、いくつかを説明します。
第三引数(largest): 制約変数のリストのどの変数から値を設定していくかを決める際、「一番大きい数値をドメインに持つ変数」から決定していきます（すべての変数のドメインが0 or 1なのでこのケースでは意味ないです。コピペしてきたので…）
第四引数(indomain_reverse_split) :制約変数に値を設定していくときの順番ですが、「ドメイン内の一番大きいものから降順に」値を設定していきます(このケースでは1を設定して、そのあと0を設定する）
第五引数(complete) :もれなくすべての選択肢を試します（オプションで端折るものがあります）

実行結果：

?- solve(14, Values, [5, 7, 4, 3], [16, 22, 12, 8], SumRevenues).
Values = [0, 1, 1, 1]
SumRevenues = 42
Yes (0.02s cpu)

Found a solution with cost -38
Found a solution with cost -42
Found no solution with cost -58.0 .. -43.0

6行目からのFound a solution with cost ～　の部分で繰り返しコストの条件を厳しくしながらサーチしているのがわかります。
Valuesを見ると２，３，４番目の投資を行うと最大利益42を得られるようです。

上記の問題を｢0/1個口の投資」ではなく「0 ～ N個口の投資を行う」という条件に変更した際の解をサーチするプログラムもついでに掲載します。
プログラム：

:- lib(ic).
:- lib(branch_and_bound).

model(Capacity, Values, Amounts, Revenues, SumRevenues) :-
	( foreach(Amount, Amounts), foreach(Value, Values),  param(Capacity) do
	    % 同じ種類の投資を複数回行える条件に変更
	    Bound is Capacity // Amount,
	    Value #:: 0..Bound
	),
	integers(Amounts),
	sum(Amounts*Values) $= SumAmounts,
	SumAmounts $=< Capacity,
	sum(Values*Revenues) $= SumRevenues.

solve(Capacity, Values, Amounts, Revenue, SumRevenues) :-
	model(Capacity, Values, Amounts, Revenue, SumRevenues),
	Cost #= -SumRevenues,
	minimize(search(Values, 0, largest, indomain_reverse_split, complete, []), Cost).

7行目のところが変更部分ですが、// Amountの「//」はコメントではなくて演算子で、　「Capacity / Amount の整数部分」を計算します。

実行結果：

?- solve(14, Values, [5, 7, 4, 3], [16, 22, 12, 8], SumRevenues).
Values = [2, 0, 1, 0]
SumRevenues = 44
Yes (0.00s cpu)

Found a solution with cost -32
Found a solution with cost -40
Found a solution with cost -42
Found a solution with cost -44
Found no solution with cost -144.0 .. -45.0

やはり解は変わり、44が最大利益となり、1番目の投資先に2個、3番目の投資先に1個投資するのが一番利益が高くなるようです。

ついでに、ネットで見かけた複雑そうな30個の荷物の0-1ナップサック問題を解いてみました（プログラムは前半で掲載したもの使いました）
実行結果：

?- time(solve(499887702, Values, [137274936, 989051853, 85168425, 856699603, 611065509, 22345022, 678298936, 616908153, 28801762, 478675378, 706900574, 738510039, 135746508, 599020879, 738084616, 545330137, ...], [128990795, 575374246, 471048785, 640066776, 819841327, 704171581, 536108301, 119980848, 117241527, 325850062, 623319578, 998395208, 475707585, 863910036, 340559411, 122579234, ...], SumRevenues)).
Values = [0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, ...]
SumRevenues = 3673016420
Yes (0.00s cpu)

Found a solution with cost -2593529208
Found a solution with cost -2967709379
Found a solution with cost -3056523353
Found a solution with cost -3188370729
Found a solution with cost -3242744021
Found a solution with cost -3243814047
Found a solution with cost -3325229604
Found a solution with cost -3326299630
Found a solution with cost -3403240143
Found a solution with cost -3535087519
Found a solution with cost -3617573102
Found a solution with cost -3671946394
Found a solution with cost -3673016420
Found no solution with cost -6192818779.0 .. -3673016421.0

Success, times: 0.000000s thread, 0.000+0.031s process, 0.036s real

自分のマシン(Intel(R) Core(TM) i5-3340M CPU @ 2.70GHz 2.70 GHz　メモリ8GB しょぼめのノートPC)では36msecで解けていました。

Example の Knapsack のページ を見ると、0/1のナップサック問題（上記の0個/1個口の投資しか許されない問題）の場合には別の解き方を行っていたり、EPLEXライブラリを使用した書き方もできることがわかります。また、ナップサック問題は分枝限定法ではなく動的計画法を用いて解くのもメジャーな解き方のようです。

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 – [～]の初めの数字はかかるコストが出てます。
続くリストは、どのエッジを通ったかをリストアップしてるのですが、目的地（ゴール）～ノード１（スタート）の順で出てます

以上