M.Hiroi's Home Page

Prolog Programming

制約論理プログラミング超入門


Copyright (C) 2016-2023 Makoto Hiroi
All rights reserved.

●経路の探索 (ハミルトン路)

あるグラフにおいて、すべての頂点をちょうど一回ずつ通る経路を「ハミルトン路 (Hamiltonian path)」といいます。また、すべての頂点を一回ずつ通って出発点に戻る経路を「ハミルトン閉路」といいます。任意のグラフがハミルトン路を持っているか否かを判定する方法は、現在では知られていません (NP 問題といわれています)。

今回は前回の経路図で A から I までのハミルトン路を求めてみましょう。経路図を再掲します。

A───B───C
│      │      │
│      │      │
│      │      │
D───E───F
│      │      │
│      │      │
│      │      │
G───H───I

    図 : 経路図

●プログラムの作成

前回のプログラムは、スタートとゴール以外の頂点の制約条件で、つながる辺の本数を 0 or 2 に指定しました。今回はすべての頂点を通るので、辺の本数を 2 に指定すれば解けるような気がします。実際にプログラムを作って確かめてみましょう。

リスト : ハミルトン路を求める (No Good)

:- use_module(library(clpfd)).

% 辺の表示
print_edge([], []).
print_edge([1 | Xs], [Y | Ys]) :- writeln(Y), print_edge(Xs, Ys).
print_edge([0 | Xs], [_ | Ys]) :- print_edge(Xs, Ys).

check1(_, _, []).
check1(S, G, [S-Xs | As]) :-
    sum(Xs, #=, 1), check1(S, G, As).
check1(S, G, [G-Xs | As]) :-
    sum(Xs, #=, 1), check1(S, G, As).
check1(S, G, [_-Xs | As]) :-
    sum(Xs, #=, 2), check1(S, G, As).

solver1 :-
    Vars = [AB, BC, AD, BE, CF, DE,
            EF, DG, EH, FI, GH, HI],
    Name = ["A-B", "B-C", "A-D", "B-E", "C-F", "D-E",
            "E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
    Vars ins 0..1,
    Adj = [a-[AB, AD], b-[AB, BC, BE], c-[BC, CF],
           d-[AD, DE, DG], e-[BE, DE, EF, EH], f-[CF, EF, FI],
           g-[DG, GH], h-[EH, GH, HI], i-[FI, HI]],
    check1(a, i, Adj),
    label(Vars),
    writeln(Vars),
    print_edge(Vars, Name).

述語 check1 において、最後の規則でつながる辺の本数を 2 に指定します。あとは、label で Vars の値を求めるだけです。

実行結果は次のようになりました。

?- solver1.
[0,1,1,1,1,0,0,1,1,1,1,0]
B-C
A-D
B-E
C-F
D-G
E-H
F-I
G-H
true ;
[0,1,1,1,1,0,1,1,0,0,1,1]
B-C
A-D
B-E
C-F
E-F
D-G
G-H
H-I
true ;
false.

?-
  A  B--C    A  B--C
  |  |  |    |  |  |
  D  E  F    D  E--F
  |  |  |    |
  G--H  I    G--H--I

最初の解はハミルトン路になっていますが、次の解は閉路 (B - C - E - F) ができるため、ハミルトン路にはなりません。このように、頂点につながる辺の本数だけでは条件不足になるのです。

●閉路のチェック

閉路ができるのを防ぐため、今回は経路に沿って頂点に番号を順番に付けていくことにします。このとき、隣り合った頂点の差分は 1 になることに注意してください。閉路の場合、たとえば B - C - F - E で番号を 2, 3, 4, 5 と付けると、B と E は隣り合っているにもかかわらず差分が 3 になるため、この経路は除外されることになります。

プログラムは次のようになります。

リスト : ハミルトン路を求める (Good)

check2(E, X-Y) :- E #==> abs(X - Y) #= 1.

solver2 :-
    Vars = [AB, BC, AD, BE, CF, DE,
            EF, DG, EH, FI, GH, HI],
    Name = ["A-B", "B-C", "A-D", "B-E", "C-F", "D-E",
            "E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
    Vars ins 0..1,
    Adj = [a-[AB, AD], b-[AB, BC, BE], c-[BC, CF],
           d-[AD, DE, DG], e-[BE, DE, EF, EH], f-[CF, EF, FI],
           g-[DG, GH], h-[EH, GH, HI], i-[FI, HI]],
    check1(a, i, Adj),
    % 経路に沿って頂点に番号をつける
    Node = [A, B, C, D, E, F, G, H, I],
    Node ins 1..9,
    Edge = [A-B, B-C, A-D, B-E, C-F, D-E,
            E-F, D-G, E-H, F-I, G-H, H-I],
    A #= 1,
    I #= 9,
    maplist(check2, Vars, Edge),
    label(Vars),
    writeln(Vars),
    writeln(Node),
    print_edge(Vars, Name).

リスト Node に頂点の表す変数を格納します。範囲は 1..9 で、A を 1 に、I を 9 に設定します。そして、辺の両端にある頂点の変数を組で表し、それを連想リスト Edge に格納します。辺の値が 1 ならば、両端の変数の差分が 1 であることを、述語 check2 と maplist で設定します。

●実行結果

それでは実行してみましょう。

?- solver2.
[0,1,1,1,1,0,0,1,1,1,1,0]
[1,6,7,2,5,8,3,4,9]
B-C
A-D
B-E
C-F
D-G
E-H
F-I
G-H
true ;
[1,1,0,0,1,1,1,1,0,0,1,1]
[1,2,3,6,5,4,7,8,9]
A-B
B-C
C-F
D-E
E-F
D-G
G-H
H-I
true ;
false.
  A  B--C    A--B--C
  |  |  |          |
  D  E  F    D--E--F
  |  |  |    |
  G--H  I    G--H--I

2 通りの経路を求めることができました。

●有向グラフの場合

有向グラフの場合でも同様な問題が生じます。次のプログラムのように、スタートとゴール以外の頂点で、入力の辺と出力の辺が 1 であることを条件にすると、閉路が生じてしまいます。

リスト : ハミルトン路を求める (No Good)

% 制約条件
check1(_, _, _, [], []).
check1(S, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 1),   % スタートは出る辺しかない
    sum(Y, #=, 0),
    N1 #= S + 1,
    check1(N1, S, G, Xs, Ys).
check1(G, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 0),   % ゴールは入る辺しかない
    sum(Y, #=, 1),
    N1 #= G + 1,
    check1(N1, S, G, Xs, Ys).
check1(N, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 1),   % 出る辺と入る辺の数が同じ (1)
    sum(Y, #=, 1),
    N1 #= N + 1,
    check1(N1, S, G, Xs, Ys).

% A から I までのハミルトン路を求める
solver1 :-
    adjacent(Xs),
    flatten(Xs, Ys),
    include(var, Ys, Vars),     % 変数を集める  
    Vars ins 0..1,
    transpose(Xs, Zs),
    check1(1, 1, 9, Xs, Zs),
    label(Vars),
    print_edge(0, Xs).

プログラムの実行結果を示します。

?- solver1.
A -> D
B -> E
C -> F
D -> G
E -> B
F -> C
G -> H
H -> I
true ;
A -> D
B -> E
C -> B
D -> G
E -> F
F -> C
G -> H
H -> I
true ;
A -> D
B -> C
C -> F
D -> G
E -> B
F -> I
G -> H
H -> E
true .
  A    B    C    A    B ← C   A    B → C
  ↓   ↓↑  ↓↑   ↓    ↓    ↑   ↓    ↑    ↓
  D    E    F    D    E → F   D    E    F
  ↓                ↓               ↓    ↑    ↓
  G → H → I    G → H → I   G → H    I

無向グラフと違って、2 つの頂点を行き来する閉路もあります。

この場合も経路に沿って頂点に番号をつけることで、閉路の生成を防ぐことができます。プログラムは次のようになります。

リスト : ハミルトン路を求める (Good)

check3(_, [], []).
check3(X, [Y | Ys], [N | Ns]) :-
    var(Y), !, Y #==> X + 1 #= N, check3(X, Ys, Ns).
check3(X, [_ | Ys], [_ | Ns]) :-
    check3(X, Ys, Ns).

check2([], [], _).
check2([N | Ns], [X | Xs], Node) :-
    check3(N, X, Node),
    check2(Ns, Xs, Node).

solver2 :-
    adjacent(Xs),
    flatten(Xs, Ys),
    include(var, Ys, Vars),     % 変数を集める  
    Vars ins 0..1,
    transpose(Xs, Zs),
    check1(1, 1, 9, Xs, Zs),
    length(Node, 9),
    Node ins 1..9,
    check2(Node, Xs, Node),
    element(1, Node, 1),
    element(9, Node, 9),
    label(Vars),
    writeln(Node),
    print_edge(0, Xs).

頂点を表す変数は length で生成し、述語 check2 で制約を設定します。第 1, 3 引数が頂点を表す変数のリスト、第 2 引数が隣接行列です。check2 では、変数と行を順番に取り出して check3 に渡します。check3 の第 1 引数が頂点で、第 2 引数が頂点につながる辺を格納したリスト、第 3 引数がもう一方の頂点になります。辺 Y が自由変数の場合、Y を選択したときは X から N に進むので X + 1 #= N を満たせばいいわけです。

あとは Node の先頭要素の値が 1 で、末尾要素の値が 9 であることを clpfd の述語 element/3 で指定します。

element(N, Xs, V)

element は nth1 と同様にリストの要素を 1 から数えます。そして、Xs の N 番目の要素が V であれば成功します。

簡単な例を示しましょう。

?- element(3, [10, 20, 30, 40, 50], 30).
true.

?- element(4, [10, 20, 30, 40, 50], N).
N = 40.

?- element(N, [10, 20, 30, 40, 50], 50).
N = 5.

?- element(N, [10, 20, 30, 40, 50], M), label([N, M]).
N = 1,
M = 10 ;
N = 2,
M = 20 ;
N = 3,
M = 30 ;
N = 4,
M = 40 ;
N = 5,
M = 50.

?- length(Xs, 3), Xs ins 1..3, element(1, Xs, 3), label(Xs).
Xs = [3, 1, 1] ;
Xs = [3, 1, 2] ;
Xs = [3, 1, 3] ;
Xs = [3, 2, 1] ;
Xs = [3, 2, 2] ;
Xs = [3, 2, 3] ;
Xs = [3, 3, 1] ;
Xs = [3, 3, 2] ;
Xs = [3, 3, 3].

●実行結果

実行結果は次のようになりました。

?- solver2.
[1,6,7,2,5,8,3,4,9]
A -> D
B -> C
C -> F
D -> G
E -> B
F -> I
G -> H
H -> E
true ;
[1,2,3,6,5,4,7,8,9]
A -> B
B -> C
C -> F
D -> G
E -> D
F -> E
G -> H
H -> I
true ;
false.

?-

当然ですが 2 通りの経路が出力されます。


●プログラムリスト1

リスト : 経路の探索 (無向グラフ)

:- use_module(library(clpfd)).

% 辺の表示
print_edge([], []).
print_edge([1 | Xs], [Y | Ys]) :- writeln(Y), print_edge(Xs, Ys).
print_edge([0 | Xs], [_ | Ys]) :- print_edge(Xs, Ys).

check(_, _, []).
check(S, G, [S-Xs | As]) :-
    sum(Xs, #=, 1), check(S, G, As).
check(S, G, [G-Xs | As]) :-
    sum(Xs, #=, 1), check(S, G, As).
check(S, G, [_-Xs | As]) :-
    sum(Xs, #=, K), K in 0 \/ 2, check(S, G, As).

solver :-
    Vars = [AB, BC, AD, BE, CF, DE,
            EF, DG, EH, FI, GH, HI],
    Name = ["A-B", "B-C", "A-D", "B-E", "C-F", "D-E",
            "E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
    Vars ins 0..1,
    Adj = [a-[AB, AD], b-[AB, BC, BE], c-[BC, CF],
           d-[AD, DE, DG], e-[BE, DE, EF, EH], f-[CF, EF, FI],
           g-[DG, GH], h-[EH, GH, HI], i-[FI, HI]],
    check(a, i, Adj),
    sum(Vars, #=, L),
    labeling([min(L)], Vars),
    writeln(Vars),
    print_edge(Vars, Name).

% ハミルトン路を求める
check1(_, _, []).
check1(S, G, [S-Xs | As]) :-
    sum(Xs, #=, 1), check1(S, G, As).
check1(S, G, [G-Xs | As]) :-
    sum(Xs, #=, 1), check1(S, G, As).
check1(S, G, [_-Xs | As]) :-
    sum(Xs, #=, 2), check1(S, G, As).

solver1 :-
    Vars = [AB, BC, AD, BE, CF, DE,
            EF, DG, EH, FI, GH, HI],
    Name = ["A-B", "B-C", "A-D", "B-E", "C-F", "D-E",
            "E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
    Vars ins 0..1,
    Adj = [a-[AB, AD], b-[AB, BC, BE], c-[BC, CF],
           d-[AD, DE, DG], e-[BE, DE, EF, EH], f-[CF, EF, FI],
           g-[DG, GH], h-[EH, GH, HI], i-[FI, HI]],
    check1(a, i, Adj),
    label(Vars),
    writeln(Vars),
    print_edge(Vars, Name).

check2(E, X-Y) :- E #==> abs(X - Y) #= 1.

solver2 :-
    Vars = [AB, BC, AD, BE, CF, DE,
            EF, DG, EH, FI, GH, HI],
    Name = ["A-B", "B-C", "A-D", "B-E", "C-F", "D-E",
            "E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
    Vars ins 0..1,
    Adj = [a-[AB, AD], b-[AB, BC, BE], c-[BC, CF],
           d-[AD, DE, DG], e-[BE, DE, EF, EH], f-[CF, EF, FI],
           g-[DG, GH], h-[EH, GH, HI], i-[FI, HI]],
    check1(a, i, Adj),
    %
    Node = [A, B, C, D, E, F, G, H, I],
    Node ins 1..9,
    Edge = [A-B, B-C, A-D, B-E, C-F, D-E,
            E-F, D-G, E-H, F-I, G-H, H-I],
    A #= 1,
    I #= 9,
    maplist(check2, Vars, Edge),
    label(Vars),
    writeln(Vars),
    writeln(Node),
    print_edge(Vars, Name).

●プログラムリスト2

リスト : 経路の探索 (有向グラフ)

:- use_module(library(clpfd)).

%          A B C  D E F  G H I
adjacent([[0,_,0, _,0,0, 0,0,0],   % A
          [_,0,_, 0,_,0, 0,0,0],   % B
          [0,_,0, 0,0,_, 0,0,0],   % C
          [_,0,0, 0,_,0, _,0,0],   % D
          [0,_,0, _,0,_, 0,_,0],   % E
          [0,0,_, 0,_,0, 0,0,_],   % F
          [0,0,0, _,0,0, 0,_,0],   % G
          [0,0,0, 0,_,0, _,0,_],   % H
          [0,0,0, 0,0,_, 0,_,0]]). % I

% 制約条件
check(_, _, _, [], []).
check(S, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 1),   % スタートは出る辺しかない
    sum(Y, #=, 0),
    N1 #= S + 1,
    check(N1, S, G, Xs, Ys).
check(G, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 0),   % ゴールは入る辺しかない
    sum(Y, #=, 1),
    N1 #= G + 1,
    check(N1, S, G, Xs, Ys).
check(N, S, G, [X | Xs], [Y | Ys]) :-
    K in 0..1,
    sum(X, #=, K),   % 出る辺と入る辺の数が同じ (0 or 1)
    sum(Y, #=, K),
    N1 #= N + 1,
    check(N1, S, G, Xs, Ys).

% 辺の表示
print_edge(_, []).
print_edge(N, [X | Xs]) :-
    sumlist(X, 0),
    N1 is N + 1,
    print_edge(N1, Xs).
print_edge(N, [X | Xs]) :-
    nth0(M, X, 1),
    C1 is N + 0x41,
    C2 is M + 0x41,
    format('~c -> ~c~n', [C1, C2]),
    N1 is N + 1,
    print_edge(N1, Xs).

% A から I までの最短経路を求める
solver :-
    adjacent(Xs),
    flatten(Xs, Ys),
    include(var, Ys, Vars),     % 変数を集める  
    Vars ins 0..1,
    transpose(Xs, Zs),
    check(1, 1, 9, Xs, Zs),
    sum(Vars, #=, L),           % 経路長 L
    labeling([min(L)], Vars),
    maplist(writeln, Xs),
    print_edge(0, Xs).

% 制約条件
check1(_, _, _, [], []).
check1(S, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 1),   % スタートは出る辺しかない
    sum(Y, #=, 0),
    N1 #= S + 1,
    check1(N1, S, G, Xs, Ys).
check1(G, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 0),   % ゴールは入る辺しかない
    sum(Y, #=, 1),
    N1 #= G + 1,
    check1(N1, S, G, Xs, Ys).
check1(N, S, G, [X | Xs], [Y | Ys]) :-
    sum(X, #=, 1),   % 出る辺と入る辺の数が同じ (1)
    sum(Y, #=, 1),
    N1 #= N + 1,
    check1(N1, S, G, Xs, Ys).

% ハミルトン路を求める (No Good)
solver1 :-
    adjacent(Xs),
    flatten(Xs, Ys),
    include(var, Ys, Vars),     % 変数を集める  
    Vars ins 0..1,
    transpose(Xs, Zs),
    check1(1, 1, 9, Xs, Zs),
    label(Vars),
    maplist(writeln, Xs),
    print_edge(0, Xs).

check3(_, [], []).
check3(X, [Y | Ys], [N | Ns]) :-
    var(Y), !, Y #==> X + 1 #= N, check3(X, Ys, Ns).
check3(X, [_ | Ys], [_ | Ns]) :-
    check3(X, Ys, Ns).

check2([], [], _).
check2([N | Ns], [X | Xs], Node) :-
    check3(N, X, Node),
    check2(Ns, Xs, Node).

% ハミルトン路を求める (Good)
solver2 :-
    adjacent(Xs),
    flatten(Xs, Ys),
    include(var, Ys, Vars),     % 変数を集める  
    Vars ins 0..1,
    transpose(Xs, Zs),
    check1(1, 1, 9, Xs, Zs),
    length(Node, 9),
    Node ins 1..9,
    check2(Node, Xs, Node),
    element(1, Node, 1),
    element(9, Node, 9),
    label(Vars),
    maplist(writeln, Xs),
    writeln(Node),
    print_edge(0, Xs).

●経路の探索 (ハミルトン閉路)

次は下図に示す経路図でハミルトン閉路を求めてみましょう。

A───B───C
│      │    /│
│      │  /  │
│      │/    │
D───E───F
│    /│      │
│  /  │      │
│/    │      │
G───H───I

   図 : 経路図

今回求める経路は閉路なので、すべての頂点に辺が 2 本つながることになります。無向グラフでプログラムすると、次のようになります。

リスト : ハミルトン閉路を求める

:- use_module(library(clpfd)).

% 辺の表示
print_edge([], []).
print_edge([1 | Xs], [Y | Ys]) :- writeln(Y), print_edge(Xs, Ys).
print_edge([0 | Xs], [_ | Ys]) :- print_edge(Xs, Ys).

% 制約
check1(Xs) :- sum(Xs, #=, 2).
check2(E, X-Y) :-
    E #==> (abs(X - Y) #= 1 #\/ X #= 1).  % [補足]

solver :-
    Vars = [AB, BC, AD, BE, CE, CF, DE,
            EG, EF, DG, EH, FI, GH, HI],
    Name = ["A-B", "B-C", "A-D", "B-E", "C-E", "C-F", "D-E",
            "E-G", "E-F", "D-G", "E-H", "F-I", "G-H", "H-I"],
    Vars ins 0..1,
    Adj = [[AB, AD], [AB, BC, BE], [BC, CE, CF],
           [AD, DE, DG], [BE, CE, DE, EG, EF, EH], [CF, EF, FI],
           [DG, EG, GH], [EH, GH, HI], [FI, HI]],
    maplist(check1, Adj),
    Node = [A, B, C, D, E, F, G, H, I],
    Node ins 1..9,
    A #= 1,
    B #= 2,
    D #= 9,
    AB #= 1,
    AD #= 1,
    all_different(Node),  % 必要
    Edge = [A-B, B-C, A-D, B-E, C-E, C-F, D-E,
            E-G, E-F, D-G, E-H, F-I, G-H, H-I],
    maplist(check2, Vars, Edge),
    label(Vars),
    label(Node),          % 必要
    writeln(Vars),
    writeln(Node),
    print_edge(Vars, Name).

制約を表す述語 check1 は、頂点に 2 本の辺がつながることを表します。述語 check2 は、隣り合った頂点の差分が 1 で、そうでなければ出発点に戻る辺であることを表します。[補足] あとは、頂点 A, B, D を 1, 2, 9 に、辺 AB と AD を 1 に限定して解を求めます。

ただし、これだけでは不十分で、頂点の値はすべて異なることを all_differnet で指定して、label(Node) で Node の値を求める必要があります。無向グラフの場合、辺には方向がないので、隣り合った頂点の値は差分でチェックすることになります。これでは制約が不足するのかもしれません。これは有向グラフのプログラムを作るときに確かめてみましょう。

実行結果は次のようになりました。

?- solver.
[1,0,1,1,1,1,0,0,0,1,0,1,1,1]
[1,2,4,9,3,5,8,7,6]
A-B
A-D
B-E
C-E
C-F
D-G
F-I
G-H
H-I
true ;
[1,1,1,0,0,1,0,1,0,1,1,1,0,1]
[1,2,3,9,7,4,8,6,5]
A-B
B-C
A-D
C-F
E-G
D-G
E-H
F-I
H-I
true ;
[1,1,1,0,0,1,1,1,0,0,0,1,1,1]
[1,2,3,9,8,4,7,6,5]
A-B
B-C
A-D
C-F
D-E
E-G
F-I
G-H
H-I
true ;
[1,1,1,0,1,0,0,0,1,1,0,1,1,1]
[1,2,3,9,4,5,8,7,6]
A-B
B-C
A-D
C-E
E-F
D-G
F-I
G-H
H-I
true.

?-
  A--B  C    A--B--C    A--B--C    A--B--C
  |  | /|    |     |    |     |    |    /
  |  |/ |    |     |    |     |    |   /
  D  E  F    D  E  F    D--E  F    D  E--F
  |     |    | /|  |      /   |    |     |
  |     |    |/ |  |     /    |    |     |
  G--H--I    G  H--I    G--H--I    G--H--I

4 通りの経路が出力されました。

●有向グラフのプログラム

次は有向グラフを使ってプログラムを作ってみましょう。

リスト : ハミルトン閉路を求める (有向グラフ)

:- use_module(library(clpfd)).

%
% 隣接行列
%          A B C  D E F  G H I
adjacent([[0,_,0, _,0,0, 0,0,0],   % A
          [_,0,_, 0,_,0, 0,0,0],   % B
          [0,_,0, 0,_,_, 0,0,0],   % C
          [_,0,0, 0,_,0, _,0,0],   % D
          [0,_,_, _,0,_, _,_,0],   % E
          [0,0,_, 0,_,0, 0,0,_],   % F
          [0,0,0, _,_,0, 0,_,0],   % G
          [0,0,0, 0,_,0, _,0,_],   % H
          [0,0,0, 0,0,_, 0,_,0]]). % I

% 辺の表示
print_edge(_, []).
print_edge(N, [X | Xs]) :-
    sumlist(X, 0),
    N1 is N + 1,
    print_edge(N1, Xs).
print_edge(N, [X | Xs]) :-
    nth0(M, X, 1),
    C1 is N + 0x41,
    C2 is M + 0x41,
    format('~c -> ~c~n', [C1, C2]),
    N1 is N + 1,
    print_edge(N1, Xs).

% 制約
check3(_, [], []).
check3(X, [Y | Ys], [N | Ns]) :-
    var(Y), !, Y #==> (X + 1 #= N #\/ N #= 1),  % [補足]
    check3(X, Ys, Ns).
check3(X, [_ | Ys], [_ | Ns]) :-
    check3(X, Ys, Ns).

check2([], [], _).
check2([N | Ns], [X | Xs], Node) :-
    check3(N, X, Node),
    check2(Ns, Xs, Node).

check1([], []).
check1([X | Xs], [Y | Ys]) :-
    sum(X, #=, 1),   % 出る辺と入る辺の数が同じ (1)
    sum(Y, #=, 1),
    check1(Xs, Ys).

solver :-
    adjacent(Xs),
    flatten(Xs, Ys),
    include(var, Ys, Vars),     % 変数を集める  
    Vars ins 0..1,
    transpose(Xs, Zs),
    check1(Xs, Zs),
    length(Node, 9),
    Node ins 1..9,
    check2(Node, Xs, Node),
    element(1, Node, 1),    % A
    element(2, Node, 2),    % B
    element(4, Node, 9),    % D
    label(Vars),
    writeln(Node),
    print_edge(0, Xs).

制約を表す述語 check1 は頂点に辺が 2 本つながることを表します。述語 check2 と check3 は、X から N へ進む辺 Y を選択したら、N の値は X + 1 になるか、出発点 (N = 1) に戻ることを表しています。[補足] あとは、solver で頂点 A, B, D の値を 1, 2, 9 に設定して、label で辺の値を求めます。

有向グラフの場合、頂点に対して all_differnet を指定する必要はありません。また、頂点は label で探索しなくても、値は一意に定まるようです。どうやら、有向グラフのほうが厳しい制約になるみたいです。

実行結果は次のようになりました。

?- solver.
[1,2,4,9,3,5,8,7,6]
A -> B
B -> E
C -> F
D -> A
E -> C
F -> I
G -> D
H -> G
I -> H
true ;
[1,2,3,9,7,4,8,6,5]
A -> B
B -> C
C -> F
D -> A
E -> G
F -> I
G -> D
H -> E
I -> H
true ;
[1,2,3,9,8,4,7,6,5]
A -> B
B -> C
C -> F
D -> A
E -> D
F -> I
G -> E
H -> G
I -> H
true ;
[1,2,3,9,4,5,8,7,6]
A -> B
B -> C
C -> E
D -> A
E -> F
F -> I
G -> D
H -> G
I -> H
true ;
false.

?-


当然ですが 4 通りの経路が出力されます。


●補足

今回のように、辺を選んだときに進行方向の頂点の番号が +1 になるか、またはスタート地点 (1) に戻る、という条件を設定した場合、頂点の個数が多くなると時間がとてもかかるようになります。スタート地点と次の地点、そして最後の地点を決定できる場合は、最後の地点からスタート地点に戻る辺を制約から除外したほうが速くなります。詳しい説明は拙作のページ「騎士の周遊」をお読みください。

●完全グラフ

任意の 2 頂点間に辺があるグラフを「完全グラフ (complete graph)」といいます。n 個の頂点の完全グラフを Kn と表します。Kn の辺の数は n(n - 1) / 2 になります。たとえば、3 個の頂点 A, B, C の完全グラフは三角形になり、4 個の頂点 A, B, C, D の完全グラフは四角形に対角線を描いた図になります。

  A────B
  │\    /│
  │  \/  │
  │  /\  │
  │/    \│
  C────D

 図 : 完全グラフ

clpfd には完全グラフのハミルトン閉路を求める述語 circuit/1 が用意されています。

circuit(Xs)

circuit の引数はリストで、リストの長さが頂点の個数に対応します。circuit は頂点に 1 から n までの番号を割り当てます。リストの添字と要素の値が頂点を表していて、添字の頂点から要素の頂点へ進む辺があると考えてください。

簡単な例を示しましょう。

?- length(Xs, 4), circuit(Xs), label(Xs).
Xs = [2, 3, 4, 1] ;
Xs = [2, 4, 1, 3] ;
Xs = [3, 1, 4, 2] ;
Xs = [3, 4, 2, 1] ;
Xs = [4, 1, 2, 3] ;
Xs = [4, 3, 1, 2] ;
false.

?-

これを図に示すと次のようになります。

 1  2  3  4 |         辺         :       閉路
------------+--------------------+------------------
 2  3  4  1 : 1-2, 2-3, 3-4, 4-1 : 1 - 2 - 3 - 4 - 1
 2  4  1  3 : 1-2, 2-4, 3-1, 4-3 : 1 - 2 - 4 - 3 - 1
 3  1  4  2 : 1-3, 2-1, 3-4, 4-2 : 1 - 3 - 4 - 2 - 1
 3  4  2  1 : 1-3, 2-4, 3-2, 4-1 : 1 - 3 - 2 - 4 - 1
 4  1  2  3 : 1-4, 2-1, 3-2, 4-3 : 1 - 4 - 3 - 2 - 1
 4  3  1  2 : 1-4, 2-3, 3-1, 4-2 : 1 - 4 - 2 - 3 - 1

  1────2      1────2      1        2
  │        │        \    /        │\    /│
  │        │          \/          │  \/  │
  │        │          /\          │  /\  │
  │        │        /    \        │/    \│
  4────3      4────3      4        3

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

逆回りを別経路と考えると 6 通りの経路になります。ところで、circuit を使うと「巡回セールスマン問題 (TSP)」を解くこともできます。

●巡回セールスマン問題

セールスマンは下図に示す各都市を 1 回ずつもれなく訪問して帰ってこなくてはなりません。このとき、一番短い巡路 (全ての都市を含む単純な閉路) を見つけるのが「巡回セールスマン問題 (Traveling Salesperson Problem, TSP)」です。巡路の長さは自由に定義してかまいません。たとえば、時間であったり料金であってもいいのですが、今回はオーソドックスに「距離」と定義することにしましょう。

    A              B
    ●              ●

        ●      ●
        E      F
    ●              ●
    C              D
 
        都市の配置

図 : 巡回セールスマン問題

上図の場合、都市の個数は 6 つあるので、巡路の総数は (6 - 1)! = 120 通りしかありません。これならば総当りで簡単に答えを求めることができます。ところが、都市の個数が多くなると簡単ではありません。n 個の都市の場合は (n - 1)! の巡路を調べなくてはいけません。巡路の総数が n の階乗に比例して増えるというのは、まさに爆発的に増えるので、n が増えると厳密解を求めるのは大変困難になります。

今回は clpfd の述語 circuit を使って、巡回セールスマン問題を解いてみましょう。

●プログラムの作成

最初は上図の配置でプログラムを作りましょう。都市を頂点とみなして、お互いの距離を行列で表すと、次のようになります。

  :    座標
--+------------
A : (  0,   0)
B : (400,   0)
C : (  0, 200)
D : (400, 200)
E : (100, 100)
F : (300, 100)


          頂点間の距離

   :  A    B    C    D    E    F
---+------------------------------
 A :   0  400  200  447  141  316 
 B : 400    0  447  200  316  141 
 C : 200  447    0  400  141  316 
 D : 447  200  400    0  316  141 
 E : 141  316  141  316    0  200 
 F : 316  141  316  141  200    0 

述語 circuit を使う場合、頂点 (A - F) を整数 (1 - 6) に割り当てます。circuit が生成するリストの要素は、次へ進む頂点を表わしています。たとえば、リストを [A, B, C, D, E, F] とすると、経路は 1 -> A, 2 -> B, 3 -> C, 4 -> D, 5 -> E, 6 -> 7 になるので、頂点間の距離は element を使って簡単に求めることができます。

プログラムは次のようになります。

リスト : 巡回セールスマン問題

:- use_module(library(clpfd)).

tsp :-
    Node = [A, B, C, D, E, F],
    element(A, [  0, 400, 200, 447, 141, 316], L1),
    element(B, [400,   0, 447, 200, 316, 141], L2),
    element(C, [200, 447,   0, 400, 141, 316], L3),
    element(D, [447, 200, 400,   0, 316, 141], L4),
    element(E, [141, 316, 141, 316,   0, 200], L5),
    element(F, [316, 141, 316, 141, 200,   0], L6), 
    L #= L1 + L2 + L3 + L4 + L5 + L6,
    circuit(Node),
    labeling([min(L)], Node),
    writeln(L),
    writeln(Node).

変数 (A - F) の値は、その頂点から次へ進む頂点を表しています。したがって、element(A, [ ... ], L1) とすれば、1 -> A の距離を変数 L1 に求めることができます。あとは、L1 から L6 を総和を変数 L に求め、circuit で Node のハミルトンを閉路を生成し、その最小値を labeling で探索すればいいわけです。

●実行結果

実行結果は次のようになりました。

?- tsp.
1282
[2,4,1,6,3,5]
true ;
1282
[3,1,5,2,6,4]
true ;
1282
[3,6,4,2,1,5]
true ;
1282
[5,4,1,3,6,2]
true ;
1314
[3,4,5,6,2,1]
true .

?-

最小値は 1282 で、逆周りを含めて 4 通りの経路になりました。

なお、circuit を使う方法は、都市の個数が多くなると時間がとてもかかるようになります。また、頂点の並べ方によっても実行時間は左右されます。ご参考までに、頂点の座標を与えて最短閉路を求めるプログラムを示します。興味のある方はいろいろためしてみてください。

リスト : 巡回セールスマン問題 (2)

:- use_module(library(clpfd)).

% Algorithms with Python: 巡回セールスマン問題[4] より
% 規則的なデータ (data16)
point_data(1, [20-20, 120-20, 220-20, 320-20,
               70-120, 170-120, 270-120, 370-120, 
               20-220, 120-220, 220-220, 320-220,
               70-320, 170-320, 270-320, 370-320]).

% 乱数データ (r15)
% 並べ方を変更 (元のままでは時間がかかりすぎる)
point_data(2, [21-17, 38-85, 41-57, 61-239, 97-23,
               336-15, 384-43, 242-178, 183-419, 425-145,
               321-276, 385-338, 365-355, 365-399, 396-356]).

make_line(_,[],[]).
make_line(X1-Y1, [X2-Y2 | Ys], [Z | Zs]) :-
    Dx is abs(X1 - X2),
    Dy is abs(Y1 - Y2),
    Z is round(sqrt(Dx * Dx + Dy * Dy)),
    make_line(X1-Y1, Ys, Zs).

make_matrix([], _, []).
make_matrix([X | Xs], Ys, [Zs1 | Zs]) :-
    make_line(X, Ys, Zs1),
    make_matrix(Xs, Ys, Zs).

solver(Q) :-
    point_data(Q, Ps),
    make_matrix(Ps, Ps, Xs),
    length(Ps, N),
    length(Vs, N),    % 頂点
    length(Cs, N),    % 距離 (コスト)
    maplist(element, Vs, Xs, Cs),
    sum(Cs, #=, L),
    circuit(Vs),
    labeling([min(L)], Vs),
    writeln(L),
    writeln(Vs).
?- time(solver(1)).
1672
[2,3,4,7,1,5,8,12,10,11,6,16,9,13,14,15]
% 19,359,902 inferences, 1.812 CPU in 1.812 seconds (100% CPU, 10683220 Lips)
true .

?- time(solver(2)).
1607
[3,4,2,9,1,5,6,10,14,7,8,11,15,13,12]
% 236,161,365 inferences, 20.870 CPU in 20.870 seconds (100% CPU, 11315851 Lips)
true .

?-

初版 2016 年 6 月 19 日
改訂 2023 年 5 月 14 日