SWI-Prolog で制約プログラミングを行うとき、述語 maplist/2 を使うとプログラムが簡単になる場合があります。
maplist(Pred, List)
maplist/2 は List の要素に述語 Pred を適用し、すべての要素が成功すれば maplist/2 も成功します。
簡単な例を示しましょう。
?- maplist(=:=(2), [2,2,2,2,2]). true. ?- maplist(=:=(2), [2,2,2,2,1]). false. ?- maplist(=\=(2), [1,3,5,7,9]). true. ?- maplist(=\=(2), [1,2,3,5,7,9]). false. ?-
maplist/2 を使うと制約を簡単に記述することができます。簡単な例題として、拙作のページ「お気楽 Prolog プログラミング入門: 地図の配色問題」を clpfd で解いてみましょう。
┌─────────┐
│ ａ │
├──┬───┬──┤
│ ｂ │ ｃ │ ｄ │
├──┴─┬─┴──┤
│ ｅ │ ｆ │
└────┴────┘
図 : 簡単な地図
出典 : Leon Sterling, Ehud Shapiro,『Prolog の技芸』, 共立出版, 1988
今回は、図に示す簡単な地図を 4 色で塗り分けます。プログラムは次のようになります。
リスト : 地図の配色問題
:- use_module(library(clpfd)).
solver([A, B, C, D, E, F]) :-
[A, B, C, D, E, F] ins 1 .. 4,
maplist(#\=(A), [B, C, D]),
maplist(#\=(B), [A, C, E]),
maplist(#\=(C), [A, B, D, E, F]),
maplist(#\=(D), [A, C, F]),
maplist(#\=(E), [B, C, F]),
maplist(#\=(F), [C, D, E]).
領域を変数 A, B, C, D, E, F に、色を数字 1, 2, 3, 4 に対応させます。あとは、隣の領域の色が自分と異なることを maplist で記述するだけです。つまり、隣接リストを定義するだけで問題を解くことができるわけです。実行結果は次のようになります。
?- solver(Xs), label(Xs). Xs = [1, 2, 3, 2, 1, 4] ; Xs = [1, 2, 3, 2, 4, 1] ; Xs = [1, 2, 3, 4, 1, 2] ; Xs = [1, 2, 3, 4, 4, 1] ; Xs = [1, 2, 3, 4, 4, 2] . ?-
┌─────────┐ │ ■ │ ├──┬───┬──┤ │ ■ │ ■ │ ■ │ ├──┴─┬─┴──┤ │ ■ │ ■ │ └────┴────┘ 図 : 解答の一例
なお、一番最後の制約は記述しなくても解くことができます。ちょっと冗長な記述でも、矛盾していなければ clpfd で解を求めることができます。
次は、もう少し大きな地図で試してみましょう。
┌──────┬──────┐
│ Ａ │ Ｂ │
│ ┌──┬─┴─┬──┐ │
│ │ │ Ｄ │ │ │
│ │ Ｃ ├─┬─┤ Ｅ │ │
│ │ │ │ │ │ │
│ ├──┤Ｇ│Ｈ├──┤ │
│ │ │ │ │ │ │
│ │ Ｆ ├─┴─┤ Ｉ │ │
│ │ │ Ｊ │ │ │
│ ├──┴─┬─┴──┤ │
│ │ Ｋ │ Ｌ │ │
│ └────┴────┴─┤
│ │
└─────────────┘
図 : 簡単な地図 (2)
上の地図を 4 色で塗り分けます。プログラムと実行結果を示します。
リスト : 地図の配色問題 (2)
solver1(Xs) :-
Xs = [A,B,C,D,E,F,G,H,I,J,K,L],
Xs ins 1..4,
maplist(#\=(A), [B,C,D,F,K,L]),
maplist(#\=(B), [A,D,E,I,L]),
maplist(#\=(C), [A,D,F,G]),
maplist(#\=(D), [A,B,C,E,G,H]),
maplist(#\=(E), [B,D,H,I]),
maplist(#\=(F), [A,C,G,J,K]),
maplist(#\=(G), [C,D,F,H,J]),
maplist(#\=(H), [D,E,G,I,J]),
maplist(#\=(I), [B,E,H,J,L]),
maplist(#\=(J), [F,G,H,I,K,L]),
maplist(#\=(K), [A,F,J,L]),
maplist(#\=(L), [A,B,I,J,K]).
?- solver1(Xs), label(Xs), writeln(Xs). [1,2,2,3,1,3,4,2,3,1,2,4] Xs = [1, 2, 2, 3, 1, 3, 4, 2, 3|...] ; [1,2,2,3,1,3,4,2,4,1,2,3] Xs = [1, 2, 2, 3, 1, 3, 4, 2, 4|...] ; [1,2,2,3,1,3,4,2,4,1,4,3] Xs = [1, 2, 2, 3, 1, 3, 4, 2, 4|...] . ?- findall(Xs, (solver1(Xs), label(Xs)), Ys), length(Ys, N). Ys = [[1, 2, 2, 3, 1, 3, 4, 2|...], [1, 2, 2, 3, 1, 3, 4|...], [1, 2, 2, 3, 1, 3|...], [1, 2, 2, 3, 1|...], [1, 2, 2, 3|...], [1, 2, 2|...], [1, 2|...], [1|...], [...|...]|...], N = 216. ?-
解は全部で 216 通りあります。このプログラムでは重複解のチェックを行っていないので、多数の解が出力されます。地域Ａの色を 1 に限定すると 54 通りの解となります。最初に表示される解を下図に示します。
┌──────┬──────┐
│ ● │ ● │
│ ┌──┬─┴─┬──┐ │
│ │ │ ● │ │ │
│ │ ● ├─┬─┤ ● │ │
│ │ │ │ │ │ │
│ ├──┤●│●├──┤ │
│ │ │ │ │ │ │
│ │ ● ├─┴─┤ ● │ │
│ │ │ ● │ │ │
│ ├──┴─┬─┴──┤ │
│ │ ● │ ● │ │
│ └────┴────┴─┤
│ │
└─────────────┘
図 : 解答の一例
部分和問題は、要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題です。たとえば、集合 {2, 3, 5, 8} の場合、総和が 10 となる部分集合は {2, 3, 5} と {2, 8} がありますが、14 となる部分集合はありません。部分集合の総数は、要素数を n とすると \(2^n\) 個になるので、n が大きくなるとナイーブな方法では時間がかかってしまいます。実際には、動的計画法などの手法を使って、現実的な時間で部分和問題を解くことができるといわれています。clpfd でも高速に解くことができるか試してみましょう。
最初にナイーブな方法で部分和問題を解いてみましょう。今回は要素を正整数に限定します。部分和問題は「べき集合」を生成する述語 power_set を作ると簡単です。次のリストを見てください。
リスト : べき集合 power_set([], []). power_set([_ | Xs], Ys) :- power_set(Xs, Ys). power_set([X | Xs], [X | Ys]) :- power_set(Xs, Ys).
べき集合を求める述語 power_set は簡単です。最初の規則は空集合のべき集合は空集合であることを表しています。これが再帰呼び出しの停止条件になります。2 番目の規則は先頭要素を取り除いたリスト Xs のべき集合を求めます。3 番目の規則は、先頭要素を取り除いたリスト Xs のべき集合 Ys を求め、Ys の先頭に X を追加します。これでべき集合を生成することができます。
簡単な実行例を示します。
?- power_set([2,3,5,8], Xs). Xs = [] ; Xs = [8] ; Xs = [5] ; Xs = [5, 8] ; Xs = [3] ; Xs = [3, 8] ; Xs = [3, 5] ; Xs = [3, 5, 8] ; Xs = [2] ; Xs = [2, 8] ; Xs = [2, 5] ; Xs = [2, 5, 8] ; Xs = [2, 3] ; Xs = [2, 3, 8] ; Xs = [2, 3, 5] ; Xs = [2, 3, 5, 8]. ?-
[2, 3, 5, 8] の部分集合は空集合 [ ] を含めて 16 通りあります。この power_set を使うと部分和問題のプログラムは次のようになります。
リスト : 部分和問題
subset_sum(Xs, N, As) :-
power_set(Xs, As),
sum_list(As, N).
部分集合 As の総和を述語 sum_list で求め、それが N と等しいかチェックするだけです。sum_list/2 はリストの総和を求める述語です。簡単な使用例を示します。
?- sum_list([1,2,3,4,5], A). A = 15. ?- sum_list([1,2,3,4,5], 15). true. ?- sum_list([1,2,3,4,5], 16). false. ?-
それでは実行してみましょう。
?- subset_sum([2,3,5,8], 10, A). A = [2, 8] ; A = [2, 3, 5] ; false. ?- subset_sum([2,3,5,8], 14, A). false. ?-
とても簡単ですね。ただし、集合の要素数が多くなると、実行時間がかかるようになります。次のテストプログラムを見てください。
リスト : 簡単なテスト
% フィボナッチ数の生成
make_fibo(1, [1]).
make_fibo(2, [2, 1]).
make_fibo(N, Xs) :-
N > 2,
M is N - 2,
make_fibo_sub(M, [2, 1], Xs).
make_fibo_sub(0, Xs, Ys) :- !, reverse(Xs, Ys).
make_fibo_sub(N, [B, A | Xs], Ys) :-
N1 is N - 1,
C is A + B,
make_fibo_sub(N1, [C, B, A | Xs], Ys).
% 簡単なテスト
test(N, A) :-
make_fibo(N, Xs),
sum_list(Xs, M),
M1 is M - 1,
subset_sum(Xs, M1, A).
述語 make_fibo は N 個のフィボナッチ数列を生成します。簡単な実行例を示します。
?- make_fibo(10, A), writeln(A). [1,2,3,5,8,13,21,34,55,89] A = [1, 2, 3, 5, 8, 13, 21, 34, 55|...]. ?-
要素の総和を M とすると、1 から M までの整数は、要素を組み合わせて必ず作ることができます。これはフィボナッチ数列の面白い特徴です。test は 総和 - 1 となる組み合わせを subset_sum で求め、その実行時間を計測します。結果は次のようになりました。
?- time((test(16, A), writeln(A), fail)). [2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597] % 1,310,807 inferences, 0.139 CPU in 0.139 seconds (100% CPU, 9439672 Lips) false. ?- time((test(17, A), writeln(A), fail)). [2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584] % 2,752,603 inferences, 0.284 CPU in 0.284 seconds (100% CPU, 9680130 Lips) false. ?- time((test(18, A), writeln(A), fail)). [2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181] % 5,767,264 inferences, 0.596 CPU in 0.596 seconds (100% CPU, 9669156 Lips) false. ?- time((test(19, A), writeln(A), fail)). [2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765] % 12,058,725 inferences, 1.232 CPU in 1.232 seconds (100% CPU, 9784185 Lips) false. ?- time((test(20, A), writeln(A), fail)). [2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765,10946] % 25,165,930 inferences, 2.560 CPU in 2.560 seconds (100% CPU, 9829040 Lips) false. ?-
要素がひとつ増えると実行時間は約 2 倍になっていることがわかります。要素数を n とすると、subset_sum の実行時間は \(2^n\) に比例する遅いプログラムなのです。
次は clpfd でプログラムを作ってみましょう。部分和問題は集合 X の要素を \(x_i\) とし、その係数を \(a_i\)すると、次の等式を満たすか判定する問題になります。
係数 \(a_i\) の値が 0 か 1 かを決める問題になります。このような問題を「0-1 整数計画問題」といいます。clpfd にはこのような問題を解くのにぴったりの述語 scalar_product/4 が用意されています。
scalar_product(Xs, As, Pred, N).
scala product は「内積」のことです。たとえば、3 次元のベクトル (x1, y1, z1), (x2, y2, z2) の内積は x1 * x2 + y1 * y2 + z1 * z2 で求めることができます。
scala_product はリスト Xs, As の要素を掛け算して、その総和を求めます。その値と引数 N の関係が述語 Pred を満たすとき、scalar_product は成功します。実をいうと、この述語だけで部分和問題は解けてしまいます。簡単な実行例を示しましょう。
?- length(Xs, 5), scalar_product([1,2,3,4,5], Xs, #=, 15), Xs ins 0..1. Xs = [1, 1, 1, 1, 1]. ?- length(Xs, 5), scalar_product([1,2,3,4,5], Xs, #=, 10), Xs ins 0..1, label(Xs). Xs = [0, 1, 1, 0, 1] ; Xs = [1, 0, 0, 1, 1] ; Xs = [1, 1, 1, 1, 0]. ?-
これをそのままプログラムすると、部分和問題の解法プログラムは次のようになります。
リスト : 部分和問題 (clpfd 版)
subset_sum_clp(Xs, N, As) :-
length(Xs, M),
length(As, M),
As ins 0..1,
scalar_product(Xs, As, #=, N).
% 表示
print_ans([],[]) :- nl.
print_ans([X | Xs], [1 | Ys]) :- format('~d ', [X]), print_ans(Xs, Ys).
print_ans([_ | Xs], [0 | Ys]) :- print_ans(Xs, Ys).
% 簡単なテスト
test_clp(Xs, N) :-
subset_sum_clp(Xs, N, Ans),
label(Ans),
print_ans(Xs, Ans),
fail.
プログラムは簡単なので説明は割愛します。それでは実行してみましょう。
?- time((make_fibo(20, Xs), sum_list(Xs, M), M1 is M - 1, test_clp(Xs, M1))). 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 % 4,094 inferences, 0.001 CPU in 0.001 seconds (100% CPU, 6210558 Lips) false. ?- time((make_fibo(30, Xs), sum_list(Xs, M), M1 is M - 1, test_clp(Xs, M1))). 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 1346269 % 5,876 inferences, 0.002 CPU in 0.002 seconds (100% CPU, 3019527 Lips) false. ?- time((make_fibo(40, Xs), sum_list(Xs, M), M1 is M - 1, test_clp(Xs, M1))). 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 1346269 2178309 3524578 5702887 9227465 14930352 24157817 39088169 63245986 102334155 165580141 % 7,766 inferences, 0.001 CPU in 0.001 seconds (100% CPU, 6273021 Lips) false. ?-
次は最大値 (末尾の要素) - 1 を求めてみましょう。
?- time((make_fibo(20, Xs), last(Xs, M), M1 is M - 1, test_clp(Xs, M1))). 1 3 8 21 55 144 377 987 2584 6765 % 13,100 inferences, 0.002 CPU in 0.002 seconds (100% CPU, 5873913 Lips) false. ?- time((make_fibo(30, Xs), last(Xs, M), M1 is M - 1, test_clp(Xs, M1))). 1 3 8 21 55 144 377 987 2584 6765 17711 46368 121393 317811 832040 % 26,215 inferences, 0.004 CPU in 0.004 seconds (100% CPU, 6070957 Lips) false. ?- time((make_fibo(40, Xs), last(Xs, M), M1 is M - 1, test_clp(Xs, M1))). 1 3 8 21 55 144 377 987 2584 6765 17711 46368 121393 317811 832040 2178309 5702887 14930352 39088169 102334155 % 43,730 inferences, 0.006 CPU in 0.006 seconds (100% CPU, 7906346 Lips) false. ?-
最後に総和 + 1 の値、つまり失敗する場合をためしてみましょう。
?- time((make_fibo(20, Xs), sum_list(Xs, M), M1 is M + 1, test_clp(Xs, M1))). % 2,712 inferences, 0.001 CPU in 0.001 seconds (100% CPU, 5071055 Lips) false. ?- time((make_fibo(30, Xs), sum_list(Xs, M), M1 is M + 1, test_clp(Xs, M1))). % 4,002 inferences, 0.000 CPU in 0.000 seconds (100% CPU, 8200820 Lips) false. ?- time((make_fibo(40, Xs), sum_list(Xs, M), M1 is M + 1, test_clp(Xs, M1))). % 5,292 inferences, 0.001 CPU in 0.001 seconds (100% CPU, 4816601 Lips) false. ?-
このように、clpfd を使うと個数を増やしても高速に解くことができます。clpfd の制約ソルバーはとても優秀ですね。M.Hiroi も大変驚きました。
リスト : 部分和問題
:- use_module(library(clpfd)).
% べき集合
power_set([], []).
power_set([_ | Xs], Ys) :- power_set(Xs, Ys).
power_set([X | Xs], [X | Ys]) :- power_set(Xs, Ys).
% ナイーブな方法
subset_sum(Xs, N, As) :-
power_set(Xs, As),
sum_list(As, N).
% フィボナッチ数の生成
make_fibo(1, [1]).
make_fibo(2, [2, 1]).
make_fibo(N, Xs) :-
N > 2,
M is N - 2,
make_fibo_sub(M, [2, 1], Xs).
make_fibo_sub(0, Xs, Ys) :- !, reverse(Xs, Ys).
make_fibo_sub(N, [B, A | Xs], Ys) :-
N1 is N - 1,
C is A + B,
make_fibo_sub(N1, [C, B, A | Xs], Ys).
% 簡単なテスト
test(N, A) :-
make_fibo(N, Xs),
sum_list(Xs, M),
M1 is M - 1,
subset_sum(Xs, M1, A).
% clpfd 版
subset_sum_clp(Xs, N, As) :-
length(Xs, M),
length(As, M),
As ins 0..1,
scalar_product(Xs, As, #=, N).
% 表示
print_ans([],[]) :- nl.
print_ans([X | Xs], [1 | Ys]) :- format('~d ', [X]), print_ans(Xs, Ys).
print_ans([_ | Xs], [0 | Ys]) :- print_ans(Xs, Ys).
% 簡単なテスト
test_clp(Xs, N) :-
subset_sum_clp(Xs, N, Ans),
label(Ans),
print_ans(Xs, Ans),
fail.
次は「ナップザック問題」を取り上げます。ナップザック (knapsack) とは辞書を引いてみると、ランドセルのような背中にせおう四角形の袋や箱のことを意味します。ここでは物を入れる袋と簡単に考えてください。
ここで、ナップザックの中に品物を詰め込むことを考えてみます。一つのナップザックと複数の品物が与えられたとき、袋に詰めた品物の合計金額が最大になるような選び方を求めることが「ナップザック問題」です。ナップザック問題にはバリエーションがあって、同じ品物をいくつも入れて良い場合と、一つしか入れてはいけない場合があります。後者の場合を「0-1 ナップザック問題」といいます。
ナップザック問題は、部分和問題と同様に NP 問題になります。これは厳密に解を求めようとすると、全ての場合について総当たりで調べるしか方法がなく、データ数が多くなると時間がかかるため、現実的な時間では解答を出すことができないというものです。品物の詰め方が難問の一つ、といわれてもピンとこないと思いますが、ナップザック問題は品物の種類が増えるにしたがって、その組み合わせが爆発的に増えるのです。
ところが、幸いなことにナップザック問題は実用的には解決済みの問題と考えられています。とくに有名なのが「動的計画法」を用いた解法です。ナップザックと品物の大きさを整数値に限定すれば、動的計画法を用いることで厳密解を求めることが可能です。興味のある方は拙作のページ「Algorithms with Python: 動的計画法」をお読みください。
今回は clpfd を使ってナップザック問題を解いてみましょう。
まず最初に、0-1 ナップザック問題を解いてみましょう。0-1 ナップザック問題を数式で表すと次のようになります。
基本的には、部分和問題と同じようなプログラムになりますが、ナップザックに入る範囲で、最大の金額になるような入れ方を求める必要があります。
このような場合、制約ソルバーに最大値を探すためのオプションを指定する必要があります。オプションの指定には述語 labeling/2 を使います。
labeling(Opts, Vars)
第 1 引数のリストにオプションを、第 2 引数に変数を指定します。第 1 引数が空リストの場合、label/1 と同じ動作 (デフォルトの動作) になります。
最大値または最小値を求める場合は、オプション max(Expr) と min(Expr) を使います。引数の Expr には変数だけではなく式を指定してもかまいません。簡単な例を示しましょう。
?- X in 1..3, Y in 1..3, labeling([max(X + Y)], [X, Y]). X = Y, Y = 3 ; X = 2, Y = 3 ; X = 3, Y = 2 ; X = 1, Y = 3 ; X = Y, Y = 2 ; X = 3, Y = 1 ; X = 1, Y = 2 ; X = 2, Y = 1 ; X = Y, Y = 1 ; false. ?- X in 1..3, Y in 1..3, labeling([min(X + Y)], [X, Y]). X = Y, Y = 1 ; X = 1, Y = 2 ; X = 2, Y = 1 ; X = 1, Y = 3 ; X = Y, Y = 2 ; X = 3, Y = 1 ; X = 2, Y = 3 ; X = 3, Y = 2 ; X = Y, Y = 3 ; false. ?-
max(X + Y) を指定すると、X + Y の大きいほうから順番に、逆に min(X + Y) を指定すると小さいほうから順番に値を求めることができます。
それから、labeling/2 と label/1 は自由変数に値を割り当てるとき、リストの左側 (先頭) の変数から行います。たとえば、label([X, Y, Z]) であれば、X, Y, Z の順番で値を決め、バックトラックするときは、逆の順番で行われます。つまり、Z の領域をすべて試したあと、バックトラックして Y の領域を試し、最後に X の領域を試します。
簡単な実行例を示します。
?- X in 1..2, Y in 3..4, Z in 5..6, labeling([], [X, Y, Z]). X = 1, Y = 3, Z = 5 ; X = 1, Y = 3, Z = 6 ; X = 1, Y = 4, Z = 5 ; X = 1, Y = 4, Z = 6 ; X = 2, Y = 3, Z = 5 ; X = 2, Y = 3, Z = 6 ; X = 2, Y = 4, Z = 5 ; X = 2, Y = 4, Z = 6. ?-
オプションに ff を指定すると、領域の小さな変数から値を割り当てます。次の例を見てください。
?- X in 1..3, Y in 4..5, labeling([ff], [X, Y]). X = 1, Y = 4 ; X = 2, Y = 4 ; X = 3, Y = 4 ; X = 1, Y = 5 ; X = 2, Y = 5 ; X = 3, Y = 5. ?-
変数 Y の領域は X の領域よりも小さいので、最初に Y の値を決め、X の領域をすべて試したあと、バックトラックして Y の領域を試します。
それでは実際に簡単な問題を解いてみましょう。
下表に示す品物をサイズ 15 のナップザックに入れるとき、金額が最大となる入れ方を求めてください。
| 品物 | 金額 | サイズ |
|---|---|---|
| A | 4 | 3 |
| B | 5 | 4 |
| C | 6 | 5 |
| D | 8 | 7 |
| E | 10 | 9 |
出典 : Coprisによる制約プログラミング入門, (田村直之さん)
プログラムは次のようになります。
リスト : 0-1 ナップザック問題
:- use_module(library(clpfd)).
solver :-
Ws = [3, 4, 5, 7, 9],
Ps = [4, 5, 6, 8, 10],
Bs = [A, B, C, D, E],
Bs ins 0..1,
scalar_product(Ws, Bs, #=<, 15),
scalar_product(Ps, Bs, #=, V),
labeling([max(V)], [A, B, C, D, E, V]),
writeln(V),
writeln(Bs).
Ws にサイズ、Ps に金額、Bs に品物の選択結果を表す変数 (0 or 1) をセットします。次に、scalar_product でサイズの合計が 15 以下であること、金額の合計を変数 V で表すことを定義します。あとは、labeling でオプションに max(V) を指定するだけです。
それでは実行してみましょう。
?- solver. 18 [1,0,1,1,0] true ; 17 [1,1,0,1,0] true ; 16 [0,0,1,0,1] true . ?- once(solver). 18 [1,0,1,1,0] true. ?-
金額の最大値は 18 で、選択した品物は A, C, D の 3 つです。バックトラックせずに 1 回だけ実行したい場合は述語 once(Goal) を使うと便利です。once は引数の Goal を実行し、成功してもバックトラックはしません。
次は同じ品物をいくつ選んでもよい問題を解いてみましょう。
下表に示す品物をサイズ 10 のナップザックに入れるとき、金額が最大となる入れ方を求めてください。なお、同じ品物を何個入れてもかまいません。
| 品物 | 金額 | サイズ |
|---|---|---|
| A | 6 | 4 |
| B | 4 | 3 |
| C | 1 | 1 |
リスト : ナップザック問題
solver1 :-
Ws = [4, 3, 1],
Ps = [6, 4, 1],
Bs = [A, B, C],
maplist(#=<(0), Bs),
scalar_product(Ws, Bs, #=<, 10),
scalar_product(Ps, Bs, #=, V),
labeling([max(V)], [A, B, C, V]),
writeln(V),
writeln(Bs).
プログラムは簡単で、変数 (係数) A, B, C の範囲が 0 or 1 ではなく、0 からナップザックに入るまでの個数 (上限値) になるだけです。上限値はナップザックの大きさから決定できるので、変数の範囲は 0 以上であることを maplist で指定するだけです。
それでは実行してみましょう。
?- solver1. 14 [1,2,0] true ; 14 [2,0,2] true . 13 [0,3,1] true .
このように、clpfd で簡単に求めることができますが、実をいうと、ナップザックの大きさが増えると、このプログラムでは時間がかかるようになります。0-1 ナップザック問題でも、品物の種類が増えると、さきほどのプログラムでは時間がかかるかもしれません。
それでは、今までのプログラムでは時間がかかる例を示しましょう。
下表に示す品物をサイズ W のナップザックに入れるとき、金額が最大となる入れ方を求めてください。なお、同じ品物を何個入れてもかまいません。
| 品物 | 金額 | サイズ |
|---|---|---|
| A | 91 | 3 |
| B | 120 | 4 |
| C | 610 | 20 |
| D | 930 | 30 |
この問題を今までと同じようにプログラムすると、W が大きくなるにしたがい時間がかかるようになります。
リスト : ナップザック問題の解法 (Bad)
solver_bad(W) :-
Price = [91, 120, 610, 930],
Size = [3, 4, 20, 30],
Vars = [A, B, C, D],
maplist(#=<(0), Vars),
scalar_product(Size, Vars, #=<, W),
scalar_product(Price, Vars, #=, Value),
labeling([max(Value)], [A, B, C, D, Value]),
writeln(Value),
writeln(Vars).
?- time(once(solver_bad(300))). 9300 [0,0,0,10] % 4,294,909 inferences, 0.494 CPU in 0.494 seconds (100% CPU, 8685452 Lips) true. ?- time(once(solver_bad(400))). 12392 [2,1,0,13] % 9,724,644 inferences, 1.135 CPU in 1.135 seconds (100% CPU, 8571058 Lips) true. ?- time(once(solver_bad(500))). 15490 [0,0,1,16] % 16,156,678 inferences, 1.938 CPU in 1.938 seconds (100% CPU, 8336625 Lips) true. ?-
時間がかかる理由は値 V の下限値が 0 から始まることです。最初、変数 A, B, C, D の値はすべて 0 なので、V の値も 0 になります。このあと、変数の値を一つずつ増やしながら試していくので、V の下限値がなかなか大きくならず、探索範囲を狭まることができないのです。
この場合、最初に最適値に近い V の値を求めることができれば、それが制約となって、探索範囲を狭めることができます。今回は「欲張り法」で値を求めることにしましょう。欲張り法の説明は拙作のページ「Algorithms with Python: 欲張り法」をお読みください。
欲張り法でナップザック問題を解く場合、金額 / サイズ を単価と考えて、単価の高い品物からナップザックに入れられるだけ入れていく、という方法がよく使われます。これを実現するには、labeling/2 のオプションに down を指定します。すると、変数の値は領域の大きな値から順番にセットされていきます。簡単な実行例を示します。
?- X in 1 .. 3, Y in 4 .. 5, labeling([down], [X, Y]). X = 3, Y = 5 ; X = 3, Y = 4 ; X = 2, Y = 5 ; X = 2, Y = 4 ; X = 1, Y = 5 ; X = 1, Y = 4. ?-
あとは labeling/2 で変数を指定するとき、単価の高い順に並べれば OK です。これで最初に求められる V の値は欲張り法と同じになります。プログラムは次のようになります。
リスト : ナップザック問題の解法 (Good)
solver(W) :-
% 単価
% A: 30.3, B: 30.0, C: 30.5, D: 31.0
Price = [91, 120, 610, 930],
Size = [3, 4, 20, 30],
Vars = [A, B, C, D],
maplist(#=<(0), Vars),
scalar_product(Size, Vars, #=<, W),
scalar_product(Price, Vars, #=, Value),
labeling([down, max(Value)], [D, C, A, B, Value]),
writeln(Value),
writeln(Vars).
labeling の指定で、オプションに down を指定して、変数を単価の高い順に D, C, A, B と並べるだけです。実行結果は次のようになりました。
?- time(once(solver(300))). 9300 [0,0,0,10] % 295,653 inferences, 0.032 CPU in 0.032 seconds (100% CPU, 9197366 Lips) true. ?- time(once(solver(400))). 12392 [2,1,0,13] % 549,723 inferences, 0.061 CPU in 0.061 seconds (100% CPU, 9037260 Lips) true. ?- time(once(solver(500))). 15490 [0,0,1,16] % 891,345 inferences, 0.102 CPU in 0.102 seconds (100% CPU, 8744429 Lips) true. ?-
最初のプログラムよりもずいぶんと高速に解くことができました。なお、問題のサイズが大きくなると clpfd で解くのは困難になるでしょうが、それほど大きなサイズでなければ、clpfd でも解くことができるのではないかと思いました。興味のある方はいろいろ試してみてください。