カテゴリー: ECLiPse

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バージョンがあり、こちらのほうが記述が簡潔なのでそのうち解析してみたいです。