小屋畑ソフトウェアサービス ブログページ

カテゴリー: ECLiPse

  • ECLiPSe CLPで2次元の板取り問題(2d cutting stock problem)を解く

    以下のGitHubページにソースコードを掲載しましたがこちらにも書きます。
    GitHubリンク

    2次元板取り問題というのは、例えば10×10のサイズの紙があったとして、幅と高さが以下のサイズの紙片を切り抜きたいときどのように切り抜くことが可能かを解く問題です。
    3×3
    2×1
    6×7
    5×2

    使い方

    eclipse -f cuttingstock2d.ecl -e "solve_cutting_stock_2d(入力CSVファイルパス,出力CSVファイルパス)"

    （ECLiPSeの実行ファイルのパスを前もって環境変数PATHに設定しておく必要あり)

    入力CSVファイルフォーマット(すべての数字は整数の必要あり):
    1行目: 紙の幅, 紙の高さ
    2行目: 紙片1の幅, 紙片1の高さ
    3行目: 紙片2の幅, 紙片2の高さ
    …

    出力CSVファイルフォーマット:(紙片が90度回転している可能性あり、y座標は下方向に増加)
    1行目: 紙片1の左のx座標,紙片1の上のy座標,紙片1の幅, 紙片1の高さ
    2行目: 紙片2の左のx座標,紙片2の上のy座標,紙片2の幅, 紙片2の高さ
    …

    例
    in.csvの中身:

    10,10
3,3
2,1
6,7
5,2

    実行コマンド:

    eclipse -f cuttingstock2d.ecl -e "solve_cutting_stock_2d('in.csv','out.csv')"

    標準出力:

      1  1  1  4  4  4  4  4------
  1  1  1  4  4  4  4  4------
  1  1  1---------------------
  2  3  3  3  3  3  3---------
  2  3  3  3  3  3  3---------
---  3  3  3  3  3  3---------
---  3  3  3  3  3  3---------
---  3  3  3  3  3  3---------
---  3  3  3  3  3  3---------
---  3  3  3  3  3  3---------

    out.csv の出力:

    0,0,3,3
0,3,1,2
1,3,6,7
3,0,5,2

    注意：このプログラムは内部的にバックトラックにより別解を求めることが可能ですが、上記の使用方法では最初の解しか得られません。ECLiPSe外部からの使用で別解が欲しい場合はPHPやPythonから呼ぶライブラリの使用を検討してください

    ソースコード（cuttingstock2d.ecl）

    :- lib(ic).
:- lib(csv).

solve_cutting_stock_2d(InCSV,OutCSV):-

	csv_read(InCSV, CsvRows, []),
	CsvRows = [[CWidth,CHeight]|Panels],

	(foreach(CsvPanel,Panels),fromto(PanelVars,[Panel|RestPanels],RestPanels,[]),param(CWidth,CHeight) do
		(
			[CsvWidth,CsvHeight] = CsvPanel,
			Width#>=0,
			Height#>=0,
			(Width #= CsvWidth and Height #= CsvHeight) or 
			(Width #= CsvHeight and Height #= CsvWidth),
			LeftX #:: 0..CWidth,
			TopY #:: 0..CHeight,
			0 #=<LeftX, LeftX #=< CWidth-Width,
			0 #=<TopY, TopY #=< CHeight-Height,
			Panel = [](LeftX,TopY,Width,Height)
		)
	),
	(fromto(PanelVars,[PickedupPanel|RestPanelVars],RestPanelVars,[]) do 
		(foreach(CheckPanel,RestPanelVars),param(PickedupPanel) do
			PickedupPanel = [](LeftX1,TopY1,Width1,Height1),
			CheckPanel    = [](LeftX2,TopY2,Width2,Height2),
			( LeftX2+Width2 #=< LeftX1 or LeftX1 + Width1 #=< LeftX2) or
			( TopY2+Height2 #=< TopY1 or TopY1 + Height1 #=< TopY2)

		)
	),
	term_variables(PanelVars,Flatten),
	labeling(Flatten),

	%visualization, csv out
	open(OutCSV, write, OutStrm),
	dim(VisDim,[CWidth,CHeight]),
	( foreach([](Left,Top,Width,Height),PanelVars),count(Idx,1,_),param(VisDim,OutStrm) do

	
		write(OutStrm,Left),write(OutStrm,","),
		write(OutStrm,Top),write(OutStrm,","),
		write(OutStrm,Width),write(OutStrm,","),
		writeln(OutStrm,Height),
		( for(Y,Top+1,Top+Height),param(VisDim,Idx,Width,Left)  do
			( for(X, Left+1, Left+Width),param(VisDim,Y,Idx) do
				%VisDim[X,Y] = Idx
				arg([X,Y], VisDim, Idx)
			)
		)
	),
	( for(Y,1,CHeight),param(VisDim,CWidth)  do
		(for(X, 1, CWidth),param(VisDim,Y) do
			arg([X,Y], VisDim, Elem),
			nonvar(Elem)->printf("%3d",[Elem]);printf("---",[])
		),
		nl
	),
	close(OutStrm).

    このプログラムではラベリング（制約変数への値の割り当て）をlabeling/1の使用のみとしてサボっています（値を割り当てる変数の選択はリストの先頭から順に、割り当てる値はドメインの最小値から１つづつ増やしていく）。
    ・ラベリングする変数の選択
    ・変数への値の割り当て方
    を細やかに制御することにより探索速度アップの余地があります。
    例えば「よりサイズの大きい紙片の位置から先に決定」のようなラベリング方法にすると速度の改善がみられるかもしれません。

    制約論理プログラミングでは「制約の記述」と「ラベリング順の記述」が完全に分離しているため、思いついたラベリング順のアイデアを簡単に試すことができます。labeling/1をコールするのが一番手を抜いた方法となります。

    ラベリング順の制御に関しては以下を参考にしてみてください。
    labelingに関して
    ECLiPSe CLPのsearch/6のChoiceの実験
    ECLiPSe tutorial.pdfのTree Search MethodsのChapter
    ECLiPSe e-learning course の Chapter 6 Search Strategies N-Queen Puzzle

  • ECLiPSe CLPで美術館の警備員配置問題を解く

    ECLiPSe CLPで以下の美術館の警備員配置問題を解きます。

    問題（Problems and exercises in Operations Researchの4.7 Museum Guards）:

    以下のような部屋構成の美術館に警備員を配置せよ。
    OR_4_7easy
    条件：
    ・コスト削減のため警備員は各部屋をつなぐドアの位置に配置し2つの部屋を同時に警備する。
    ・どの部屋も少なくとも1人の警備員が見張っているようにする。
    ・上記条件を満たしつつ、なるべく少ない警備員で警備する。

    プログラム：

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

solve(Vars):-
	Graph = 
	[](
		[](a,b,N1),
		[](a,g,N2),
		[](a,f,N3),
		[](a,c,N4),
		[](b,h,N5),
		[](c,d,N6),
		[](d,e,N7),
		[](d,f,N8),
		[](g,i,N9),
		[](h,i,N10),
		[](i,j,N11)
	),
	
	Nodes = [a,b,c,d,e,f,g,h,i,j],
	Vars=[N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11],
	Vars#::0..1,
		
	(
		foreach(Node,Nodes),param(Graph)
		do
		(
			(
				foreacharg([](P1,P2,P3),Graph), fromto(0,In,Out,OutFormula), param(Node)
				do
				(
					P1==Node->
						Out = In+P3;
						(
							P2==Node->
								Out=In+P3;Out=In
						)
				)
			),
			ic:(eval(OutFormula)#>=1)
		)
	),
	ic:(Sum #= sum(Vars)),
	bb_min(labeling(Vars),Sum,bb_options{timeout:60,solutions:all}).

    実行結果：

    ?- solve(Ans).
Ans = [0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1]
Yes (0.00s cpu, solution 1, maybe more)
Ans = [0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1]
Yes (0.00s cpu, solution 2, maybe more)
Ans = [0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1]
Yes (0.00s cpu, solution 3, maybe more)
Ans = [0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1]
Yes (0.00s cpu, solution 4, maybe more)
Ans = [0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1]
Yes (0.00s cpu, solution 5, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1]
Yes (0.00s cpu, solution 6, maybe more)
Ans = [0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1]
Yes (0.00s cpu, solution 7, maybe more)
No (0.00s cpu)

Found a solution with cost 6
Found no solution with cost 2.0 .. 5.0

    解答配置（の一つ）
    OR_4_7easy_ans

    解説：
    まず部屋をノード、ドアをエッジと考えたグラフを作成し、エッジごとに警備員がいるかいないかを表す変数を宣言します。警備員がいれば１、いなければ０。

    Nodesに全部屋のリストを用意しておいて(20行目）、
    それをひとつづつ指定したループを処理します(25行目）
    ループ中ではエッジ配列の中で指定した部屋を持つものを探し、警備員有無の変数をすべて足し合わせます(29行目からのループ）。ある部屋に注目したとき、そのすべてのドアに対応する警備員変数を足し合わせた式を作ります。少なくとも1つのドアには警備員がいるはずなのでこの式に「1以上」の制約をつけます（40行目）

    警備員の数を最小にするのが目標なので警備員有無変数の配列Varsの合計が最小になるように指定し分枝限定法(branch_and_bound)を用いてVarsの値を求めます。
    実行結果を見るとわかる通り、最小コストは６で、7つの解があるようです。

    もう一つ、以下のもっと複雑な部屋割りの場合も解いてみます。
    OR4_7difficult

    プログラム：

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

solve(Vars):-
	Graph = 
	[](
		[](s,r,N1),
		[](s,t,N2),
		[](r,q,N3),
		[](q,p,N4),
		[](n,l,N5),
		[](l,m,N6),
		[](l,k,N7),
		[](t,p,N8),
		[](p,o,N9),
		[](o,k,N10),
		[](k,i,N11),
		[](j,i,N12),
		[](h,g,N13),
		[](g,e,N14),
		[](x,u,N15),
		[](u,w,N16),
		[](d,c,N17),
		[](i,e,N18),
		[](c,e,N19),
		[](x,y,N20),
		[](u,z,N21),
		[](z,a,N22),
		[](a,b,N23),
		[](b,c,N24),
		[](e,f,N25),
		[](p,u,N26)
	),
	
	Nodes = [a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,w,x,y,z],
	Vars=[N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15,N16,N17,N18,N19,N20,N21,N22,N23,N24,N25,N26],
	Vars#::0..1,
		
	(
		foreach(Node,Nodes),param(Graph)
		do
		(
			(
				foreacharg([](P1,P2,P3),Graph), fromto(0,In,Out,OutFormula), param(Node)
				do
				(
					P1==Node->
						Out = In+P3;
						(
							P2==Node->
								Out=In+P3;Out=In
						)
				)
			),
			ic:(eval(OutFormula)#>=1)
		)
	),
	ic:(Sum #= sum(Vars)),
	bb_min(labeling(Vars),Sum,bb_options{timeout:60,solutions:all}),
	(
		fromto(Vars,[First|Rest],Rest,[])
		do
		write(First),write(",")
	),nl.

    実行結果：

    ?- solve(Ans).
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.00s cpu, solution 1, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 2, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 3, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 4, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 5, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 6, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 7, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 8, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 9, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 10, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 11, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.02s cpu, solution 12, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 13, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 14, maybe more)
Ans = [0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 15, maybe more)
Ans = [0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 16, maybe more)
Ans = [0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 17, maybe more)
Ans = [0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 18, maybe more)
Ans = [1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 19, maybe more)
Ans = [1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 20, maybe more)
Ans = [1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 21, maybe more)
Ans = [1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 22, maybe more)
Ans = [1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 23, maybe more)
Ans = [1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 24, maybe more)
Ans = [1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 25, maybe more)
Ans = [1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 26, maybe more)
Ans = [1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, ...]
Yes (0.03s cpu, solution 27, maybe more)
No (0.03s cpu)

Found a solution with cost 14
Found no solution with cost 8.0 .. 13.0
0,1,1,0,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,1,
0,1,1,0,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,1,
0,1,1,0,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,1,
0,1,1,0,1,1,0,0,1,0,1,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
0,1,1,0,1,1,0,0,1,0,1,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
0,1,1,0,1,1,0,0,1,0,1,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
0,1,1,0,1,1,0,0,1,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
0,1,1,0,1,1,0,0,1,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
0,1,1,0,1,1,0,0,1,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
0,1,1,0,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
0,1,1,0,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
0,1,1,0,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
0,1,1,0,1,1,1,0,1,0,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
0,1,1,0,1,1,1,0,1,0,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
0,1,1,0,1,1,1,0,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
0,1,1,1,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
0,1,1,1,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
0,1,1,1,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
1,0,0,1,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
1,0,0,1,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
1,0,0,1,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
1,0,1,0,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
1,0,1,0,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
1,0,1,0,1,1,0,1,0,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,
1,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,1,0,1,1,0,
1,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,0,1,1,0,1,0,
1,1,0,1,1,1,0,0,0,1,0,1,1,0,0,1,1,0,0,1,1,0,1,0,1,0,

    解答配置（の一つ）
    OR_4_7difficult_ans

    最小コストは１４、全部で27通りの解があるようです。

  • PrologのLogical Loopの解説

    ECLiPSe CLPでソルバーを作るときに欠かせないLogical Loopの説明をしたいと思います。

    Prologでプログラムを書いたことがない人は信じられないと思いますが、なんとPure Prologには繰り返し処理を書くための構文がなくて、ループごとに再帰処理を使った述語を定義して繰り返し処理を行います。

    例えばLstというリストに入っている要素をすべて表示するという処理を書きたい場合、以下のようにdsp_lst/1という術語を定義して繰り返しを書いて行います：

    main:-
	Lst=[a,b,c,d,e,f],
	dsp_lst(Lst).

dsp_lst([]).
dsp_lst([First|Rest]):-
	write(First),nl,
	dsp_lst(Rest).

    実行結果：

    ?- main.
Yes (0.00s cpu)
a
b
c
d
e
f

    このように繰り返しの為だけに定義した述語をauxiliary predicate（補助述語）と呼ぶようです。よくPrologのコードで出てくるのですが、述語名のうしろにアンダーバーとかを1つつけて、ループ用の補助述語を定義したりします。

    上記かなり愚かしい仕様に思えますが、恐らく以下のような理由からだと推測します。
    ①Prologの述語は本来「動作」「処理」ではなく、引数同士の「関係」を記述するためのもの（関数ではなく述語と呼ぶゆえん。一階述語論理の述語）
    ②例えば以下のような述語があった場合、

    pred(A,B,C):-
	term1(A,B,C),
	term2(A,C),
	term3(A,B).

    term1～term3のそれぞれの引数同士の関係が「すべて」成り立つとき、pred(A,B,C)が真となる。本来カンマは「処理の区切り」ではなく「複数の述語どうしを結ぶ条件の&」を表すという原則がある（orは;もしくは同名の別の述語定義)。極論を言えばterm1からterm3の順番をバラバラにしても成り立つはず（実際は大いに記述の順番に依存しますが…理想とする姿がそこという意味）、単純に複数の述語を&で結んだものに過ぎないので繰り返しとかいう概念がそもそもない。

    （余談ですが制約論理プログラミングの制約「のみ」で書かれた述語は述語内の制約の順番をバラバラに並び替えても基本的に問題なく動作します。Prologの理想に近い）

    このような仕様に一石を投じたのがECLiPSe CLPのメイン開発者のJoachim Schimpfさんの書いた以下の論文です。
    Logical Loops
    この論文によってPrologに自然な形でのループ処理の記述方法が提示され、ECLiPSeやSicstus Prologなどで採用されています。

    Logical Loopの書き方は、上記の補助述語でループを書いていた人に自然に受け入れられるような仕様となっています。逆にPrologで補助述語を使ったループを書いたことがない人は使いこなすまで時間がかかるかもしれません。再帰を使用してループを記述することに慣れている必要があります。

    例えば上記のLstというリスト内の文字をすべて表示するプログラムをLogical Loopで書くと、以下のようになります。

    main:-
	Lst=[a,b,c,d,e,f],
	(
		fromto(Lst,[First|Rest],Rest,[]) do
		(
			write(First),nl
		)
	).

    上記のように
    1.fromtoなどの繰り返し制御用の記述
    2.do
    3.(繰り返す処理)

    の3つの部分を書きます。
    補助述語との比較で理解すると良いのですが、
    「1.fromtoなどの繰り返し制御用の記述」→述語の引数の部分
    「3.繰り返す処理」→述語の内部の処理
    となります。
    補助述語の書き方をそのまま置き換えたような扱いなので、「繰り返す処理」の部分は別述語のように独立となっていて、ループの外で使われている変数などは明示的に指定しない限り一切使用することはできません。使いたい場合はparamという構文を使用して引数として渡す（後述）必要があります。

    fromtoの4つの引数の意味は(First,In,Out,Last)となっていて、やはり補助述語と関連付けて理解すると良いのですが
    In→補助述語の引数
    Out→補助述語内で再帰呼び出しする際に指定する引数
    Last→補助述語の引数がこの値でコールされたらループ終了
    という感じです。
    一番初めに補助述語をコールする際に指定する引数がFirstとなります。

    駆け足になりますがfromto以外で使用できる構文を列挙していきます。これらを複数同時に指定することもできます。

    foreach(X,List)
    リストList内の要素が順にXに設定され処理が繰り返されます。リストを生成することもできます。Xはループ内局所変数です。

    foreacharg(X,StructOrArray)
    StructOrArrayの部分はECLiPseのArray（[](1,2,3,4)みたいな記述)や複合項 fukugoukou(1,2,3,4) を指定します。
    最初の引数から順にXに設定され処理が繰り返されます。termを生成するためには使用できません。Xはループ内局所変数です。

    foreacharg(X,StructOrArray,Idx)
    引数が2個のものと同じ動きですがIdxに何番目の要素がXに入っているかの数字が入ります（arg(Idx,StructOrArray,X)がtrueとなる)。XとIdxはループ内局所変数です。

    foreachelem(X,Array)
    foreachargと同じ動きですが2次元配列などでも対応しますArray=[]([](1,2,3),[](4,5,6))などの場合Xに1,2,3,4,5,6が順に設定される。Xはループ内局所変数です。Arrayを生成するためには使用できません。

    foreachelem(X,Array,Idx)
    引数2個バージョンと同じ動きですがIdxにインデックス位置が設定されます。XとIdxはループ内局所変数です。

    foreachindex(Idx,Array)
    foreachelemと同じ動きですがXがなくIdxのみ設定されます。Idxはループ内局所変数です。

    for(I,MinExpr,MaxExpr)
    C言語などのfor(I=MinExpr;I<=MaxExpr;I++)と同じ動き。MaxExprは必ず値を設定しておく必要あり。Iはループ内局所変数です。 for(I,MinExpr,MaxExpr,Increment) 引数3個バージョンと同じ動きだがIncrementで何個づつIが増えるかを設定可能 multifor(List,MinList,MaxList) for/3と同じだがIに相当する変数をリストとして複数設定する。MinListとMaxListも対応する下限上限をリストで設定する。 Listはループ内局所変数 multifor(List,MinList,MaxList,IncrementList) 引数3つバージョンと同じだがIncrementListに増加分を設定できる count(I,Min,Max) forと同じのようだがMaxは値設定しなくてOK。ループ終了後にカウントした値など把握できる。 param(Var1,Var2,...) ループ内で使用したい変数を指定する。ここで指定しない限りループの外で使用している変数など使用できない。 補助述語に引数として渡すイメージ。 loop_name(Name) デバッグ用にループの名前を付けることができる。補助述語の名前のようなイメージ。他の述語と重複する名前は許されない。 最後に、先ほど「fromtoなどの記述は複数記述できる」と書いたのですが、これが面白くて、例えばリストを生成する用途で使用したりなどもできます。 リストList=[1,2,3,4,5]があったとして、ListSqにそれぞれの要素を2乗したリストを設定するようなプログラムを書いてみます。 

    main(LstSq):-
	Lst=[1,2,3,4,5],
	(
		fromto(Lst,[First|Rest],Rest,[]),fromto(LstSq,[FirstSq|RestSq],RestSq,[]) do
		(
			FirstSq is First*First
		)
	).

    実行結果：

    ?- main(LstSq).
LstSq = [1, 4, 9, 16, 25]
Yes (0.00s cpu)

    上記動作もfromtoの各引数に関して「補助述語に引数として渡している」イメージを持つと何が起こっているのかわかるのではないかと思います。

    以下のような補助述語が生成されてると考えると良いと思います！

    main(LstSq):-
	Lst=[1,2,3,4,5],
	auxiliary_predicate(Lst,LstSq).

auxiliary_predicate([],[]).
auxiliary_predicate([First|Rest],[FirstSq|RestSq]):-
	FirstSq is First*First,
	auxiliary_predicate(Rest,RestSq).

  • ECLiPSe CLPで輸送問題を解く

    ECLiPSe CLPで以下の輸送問題を解きます。

    問題（Problems and exercises in Operations Research 21ページ 4.5 Transportation ）：
    イタリアの運送会社は、6つの店舗(Verona,Perugia,Rome,Pescara,Taranto,Lamezia)から5つの港(Genoa,Venice,Ancona,Naples,Bari)に空のコンテナを運ばなくてはならない。

    6つの店舗が持っている空のコンテナの数は以下：
    4_5empty_containers

    5つの港がそれぞれ求めているコンテナの数は以下：
    4_5container_demand

    コンテナはトラックで送られるが、各トラックは最大で2個のコンテナしか運べない。
    また、輸送距離1kmあたり30ユーロのコストがかかる。
    各店舗から各港までの距離は以下の通り：
    4_5distance

    コストが最小となるような輸送計画をもとめよ。

    プログラム：

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

solve(Cost, AryL) :-
    % create the problem variables
    dim(AryC,[6,5]),	% containers
    dim(AryL,[6,5]),	% lorries
    CostConstant = [](
		[]( 290, 115, 355, 715, 810),
	    []( 380, 340, 165, 380, 610),
	    []( 505, 530, 285, 220, 450),
	    []( 655, 450, 155, 240, 315),
	    [](1010, 840, 550, 305,  95),
	    [](1072,1097, 747, 372, 333)
	  ),
	term_variables(temp(AryC,AryL), VarList),
    prob: (VarList $:: 0.0..1.0Inf),
    prob: integers(VarList),
    
    % Cost objective function(minimize this)
    (for(Row,1,6),fromto(0,RowSumIn,RowSumOut,WholeCost),param(CostConstant,AryL)
    	do
    	(
	    	for(Col,1,5),fromto(RowSumIn,ColIn,ColOut,RowSumOut),param(Row,CostConstant,AryL) do
	    	(
	    		Wk is CostConstant[Row,Col],
	    		ColOut = ColIn + Wk*AryL[Row,Col]
	    	)
    	)
   	),
   	% Empty containers constraint
    (
    	for(Row,1,6),foreach(EmptyContainer,[10,12,20,24,18,40]),param(AryC) do
    	(
    		(
	    		for(Col,1,5),fromto(0,ColIn,ColIn+AryC[Row,Col],ColCost),param(AryC,Row) do true
	    	),
	    	prob: (EmptyContainer $>= ColCost)
    	)
   	),
   	% Container demand constraint
   	(
   		for(Col,1,5),foreach(ContainerDemand,[20,15,25,33,21]),param(AryC) do
   		(
   			(
	   			for(Row,1,6),fromto(0,RowIn,RowIn+AryC[Row,Col],RowCost),param(AryC,Col) do true
	   		),
   			prob: (ContainerDemand $=< RowCost)
   		)
   	),
   	
   	% Lorries and Containers constraint 
    (
	    foreachelem(C,AryC),foreachelem(L,AryL)
	    do
	    prob:( C $=< 2*L )
	),
	
    % Set up the external solver with the objective function
    prob: eplex_solver_setup(min(30*WholeCost)),
    prob: eplex_solve(Cost), % Solve problem using external solver 
    ( foreachelem(L,AryL) 
    	do 
    	prob:eplex_var_get( L,  typed_solution, L )
    )
	.

    実行結果：

    ?- solve(Cost, Vars).
Cost = 468510.0
Vars = []([](0, 5, 0, 0, 0), [](3, 3, 1, 0, 0), [](7, 0, 0, 3, 0), [](0, 0, 12, 0, 0), [](0, 0, 0, 1, 9), [](0, 0, 0, 13, 2))
Yes (0.02s cpu)

    解答を整理すると以下となります。
    4_5solution

    解説：
    EPLexライブラリを用いて解きました。

    実は、これとほぼ同じ種類の問題がECLiPSeのTutorialにもあります（EPLexライブラリの解説の章の16.1.3 Example formulation of an MP Problem）
    この問題との違いはトラック1つにつきコンテナ2つを運ぶという点だけです。

    解説します。運ぶコンテナの数と実際にそれを運ぶトラックの数とを分けて考えることにします。
    コンテナの数の配列をAryC、トラックの数をAryLとします。
    CostConstantに輸送コストを設定しておきます(8行目)
    コンテナの数、トラックの数には 0以上であること、整数であることという制約をつけます（17,18行目）
    21行目からはコストの目的関数を生成してます。これを最小にするのが目標となります。この処理を抜けた段階で「WholeCost=コスト関数」の式が得られていて、これを後ほど最小化するようにEPLexに設定します(60行目)
    31行目からは各店舗から送られるコンテナ数の上限が持っているコンテナ数であるという制約をつけています。
    41行目からは各港に届くコンテナが要望以上であることを示す制約をつけています。
    53行目からはコンテナとトラックの個数の関係の制約をつけています（各店舗-港間で移動するコンテナ数<=トラック数*2) 60行目で最適化するコスト関数を設定しています。 61行目でEPLexに問題を解くよう指示しています。 62行目からは、AryCとAryLの変数がまだEPLexの範囲用の値を持っているのを、解いた答えの値へ変更しています。

  • ECLiPSe でユーザー定義の算術演算子を定義する

    ECLiPSe CLPでユーザー定義の算術演算子を定義して is で計算してもらえるようにする方法です。

    ちなみにisで既定で計算対象となる述語は ここ を参照してください。

    ここでは演算子の右と左の数字を直角三角形の高さと幅として、斜辺の長さを計算してくれる my_opという演算子を定義して実際に使ってます。
    プログラム：

    :- export op(200 ,yfx ,my_op).

my_op(Number1,Number2,Result):-
	Result is sqrt(Number1 * Number1 + Number2 * Number2).

    実行結果：

    ?- A is 3 my_op 4.
A = 5.0
Yes (0.00s cpu)

    演算子の木構造や優先度に関しては以下の記事など参考になります。
    [PROLOG] compound term
    [Prolog]演算子のカンマとそうでないカンマ、ドット演算子と演算子でないドット（ピリオド）

    opのyfxは3つ以上数字がある場合に木構造が左に伸びてくか右に伸びてくかを決定する部分で、
    例えば引き算で 6 – 2 -1 と記述した場合、
    (6-2)-1 なのか
    6-(2-1) なのかで答え変わってきますが、想定しているのは前者です。
    このような演算子は yfxで定義します。

    200の数字は演算子の優先順位です。200はECLiPSeでは+や-と同じ優先度です。小さいほど優先度が高い。ここは深く考えないでこの数字にしてしまいましたがちゃんと使う場合はまじめな検討が必要です。current_op/3という術語で演算子の優先度や結合の情報(yfxなど)が確認できます。

    SWI-PROLOGでも多分同様のことできるとは思うのですが、やり方わかりませんでした。情報求ム。

  • ECLiPSe CLPで最小全域木を求める

    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つのビット列のそれぞれのハミング距離は以下となる（問題からコピー）
    OR_Excercise_4_1

    ハミング距離が大きすぎなければ、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)

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

  • ECLiPSeで分枝限定法(branch and bound)ライブラリを用いてナップサック問題を解く

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

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

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

    この問題文だけでは読み取れなかったのですが、解答を見ると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のsearch/6のChoiceの実験

    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 ライブラリを用いて混合整数線形計画問題を解く

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

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

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

    OR_excersize_3_2
    ほぼ線形計画問題と同じですが、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 ライブラリを用いて線形計画法の問題を解く

    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を解いてみます。

    問題①
    OR_excersize_2_1
    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}の部分は　変数名{最小値..最大値 @ 解}　の意味
    くらいです。あとは見ればわかると思います。

    問題②：
    OR_excersize_2_2

    プログラム：

    :- 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)

    問題③
    OR_excersize_2_3

    プログラム：

    :- 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)

    カンタンですね！