针对数字板拼图优化的 CLP(FD) 求解器

Optimized CLP(FD) solver for number board puzzle

考虑来自https://puzzling.stackexchange.com/questions/20238/explore-the-square-with-100-hops的问题:

Given a grid of 10x10 squares, your task is to visit every square exactly once. In each step, you may

  • skip 2 squares horizontally or vertically or
  • skip 1 square diagonally

换句话说(更接近我下面的实现),用从 1 到 100 的数字标记一个 10x10 网格,使得坐标 (X, Y) 处的每个方格都是 1 或等于 "previous" 广场 (X, Y-3)(X, Y+3)(X-3, Y)(X+3, Y)(X-2, Y-2)(X-2, Y+2)(X+2, Y-2)(X+2, Y+2).

这看起来像是一个简单的约束编程问题,Z3 可以通过简单的声明性规范在 30 秒内解决它:https://twitter.com/johnregehr/status/1070674916603822081

我在 SWI-Prolog 中使用 CLP(FD) 的实现不能很好地缩放。事实上,它甚至无法解决问题的 5x5 实例,除非预先指定几乎两行:

?- number_puzzle_(_Square, Vars), Vars = [1,24,14,2,25, 16,21,5,8,20 |_], time(once(labeling([], Vars))).
% 10,063,059 inferences, 1.420 CPU in 1.420 seconds (100% CPU, 7087044 Lips)
_Square = square(row(1, 24, 14, 2, 25), row(16, 21, 5, 8, 20), row(13, 10, 18, 23, 11), row(4, 7, 15, 3, 6), row(17, 22, 12, 9, 19)),
Vars = [1, 24, 14, 2, 25, 16, 21, 5, 8|...].

?- number_puzzle_(_Square, Vars), Vars = [1,24,14,2,25, 16,21,5,8,_ |_], time(once(labeling([], Vars))).
% 170,179,147 inferences, 24.152 CPU in 24.153 seconds (100% CPU, 7046177 Lips)
_Square = square(row(1, 24, 14, 2, 25), row(16, 21, 5, 8, 20), row(13, 10, 18, 23, 11), row(4, 7, 15, 3, 6), row(17, 22, 12, 9, 19)),
Vars = [1, 24, 14, 2, 25, 16, 21, 5, 8|...].

?- number_puzzle_(_Square, Vars), Vars = [1,24,14,2,25, 16,21,5,_,_ |_], time(once(labeling([], Vars))).
% 385,799,962 inferences, 54.939 CPU in 54.940 seconds (100% CPU, 7022377 Lips)
_Square = square(row(1, 24, 14, 2, 25), row(16, 21, 5, 8, 20), row(13, 10, 18, 23, 11), row(4, 7, 15, 3, 6), row(17, 22, 12, 9, 19)),
Vars = [1, 24, 14, 2, 25, 16, 21, 5, 8|...].

(这是在装有 SWI-Prolog 6.0.0 的旧机器上。在装有 SWI-Prolog 7.2.3 的较新机器上,它的运行速度大约是它的两倍,但这还不足以击败明显的指数复杂性。)

这里使用的部分方案来自https://www.nurkiewicz.com/2018/09/brute-forcing-seemingly-simple-number.html

所以,我的问题是:我怎样才能加快以下 CLP(FD) 程序的速度?

额外感谢的其他问题:是否有一个特定的标记参数可以显着加快此搜索速度,如果有,我如何做出有根据的猜测它可能是哪一个?

:- use_module(library(clpfd)).

% width of the square board
n(5).


% set up a term square(row(...), ..., row(...))
square(Square, N) :-
    length(Rows, N),
    maplist(row(N), Rows),
    Square =.. [square | Rows].

row(N, Row) :-
    functor(Row, row, N).


% Entry is the entry at 1-based coordinates (X, Y) on the board. Fails if X
% or Y is an invalid coordinate.
square_coords_entry(Square, (X, Y), Entry) :-
    n(N),
    0 < Y, Y =< N,
    arg(Y, Square, Row),
    0 < X, X =< N,
    arg(X, Row, Entry).


% Constraint is a CLP(FD) constraint term relating variable Var and the
% previous variable at coordinates (X, Y). X and Y may be arithmetic
% expressions. If X or Y is an invalid coordinate, this predicate succeeds
% with a trivially false Constraint.
square_var_coords_constraint(Square, Var, (X, Y), Constraint) :-
    XValue is X,
    YValue is Y,
    (   square_coords_entry(Square, (XValue, YValue), PrevVar)
    ->  Constraint = (Var #= PrevVar + 1)
    ;   Constraint = (0 #= 1) ).


% Compute and post constraints for variable Var at coordinates (X, Y) on the
% board. The computed constraint expresses that Var is 1, or it is one more
% than a variable located three steps in one of the cardinal directions or
% two steps along a diagonal.
constrain_entry(Var, Square, X, Y) :-
    square_var_coords_constraint(Square, Var, (X - 3, Y), C1),
    square_var_coords_constraint(Square, Var, (X + 3, Y), C2),
    square_var_coords_constraint(Square, Var, (X, Y - 3), C3),
    square_var_coords_constraint(Square, Var, (X, Y + 3), C4),
    square_var_coords_constraint(Square, Var, (X - 2, Y - 2), C5),
    square_var_coords_constraint(Square, Var, (X + 2, Y - 2), C6),
    square_var_coords_constraint(Square, Var, (X - 2, Y + 2), C7),
    square_var_coords_constraint(Square, Var, (X + 2, Y + 2), C8),
    Var #= 1 #\/ C1 #\/ C2 #\/ C3 #\/ C4 #\/ C5 #\/ C6 #\/ C7 #\/ C8.


% Compute and post constraints for the entire board.
constrain_square(Square) :-
    n(N),
    findall(I, between(1, N, I), RowIndices),
    maplist(constrain_row(Square), RowIndices).

constrain_row(Square, Y) :-
    arg(Y, Square, Row),
    Row =.. [row | Entries],
    constrain_entries(Entries, Square, 1, Y).

constrain_entries([], _Square, _X, _Y).
constrain_entries([E|Es], Square, X, Y) :-
    constrain_entry(E, Square, X, Y),
    X1 is X + 1,
    constrain_entries(Es, Square, X1, Y).


% The core relation: Square is a puzzle board, Vars a list of all the
% entries on the board in row-major order.
number_puzzle_(Square, Vars) :-
    n(N),
    square(Square, N),
    constrain_square(Square),
    term_variables(Square, Vars),
    Limit is N * N,
    Vars ins 1..Limit,
    all_different(Vars).

首先:

这是怎么回事?

要查看发生了什么,这里有 PostScript 定义,可以让我们可视化 搜索:

/n 5 def

340 n div dup scale
-0.9 0.1 translate % leave room for line strokes

/Palatino-Roman 0.8 selectfont

/coords { n exch sub translate } bind def

/num { 3 1 roll gsave coords 0.5 0.2 translate
    5 string cvs dup stringwidth pop -2 div 0 moveto show
    grestore } bind def

/clr { gsave coords 1 setgray 0 0 1 1 4 copy rectfill
     0 setgray 0.02 setlinewidth rectstroke grestore} bind def

1 1 n { 1 1 n { 1 index clr } for pop } for

这些定义为您提供了两个程序:

  • clr清除一个方块
  • num 在正方形上显示数字。

例如,如果您将这些定义保存到 tour.ps,然后使用以下命令调用 PostScript 解释器 Ghostscript

gs -r72 -g350x350 tour.ps

然后输入以下指令:

1 2 3 num
1 2 clr
2 3 4 num

你得到:

PostScript 是一种用于可视化搜索过程的出色编程语言,我还建议查看 以了解更多信息。

我们可以很容易地修改您的程序以发出合适的 PostScript 指令,让我们直接观察搜索。我强调相关的补充:

constrain_entries([], _Square, _X, _Y).
constrain_entries([E|Es], Square, X, Y) :-
    freeze(E, postscript(X, Y, E)),
    constrain_entry(E, Square, X, Y),
    X1 #= X + 1,
    constrain_entries(Es, Square, X1, Y).

postscript(X, Y, N) :- format("~w ~w ~w num\n", [X,Y,N]).
postscript(X, Y, _) :- format("~w ~w clr\n", [X,Y]), false.

我也冒昧把(is)/2改成了(#=)/2,让程序更通用。

假设您在 tour.ps 中保存了 PostScript 定义,在 tour.pl 中保存了 Prolog 程序,下面对 SWI-Prolog 和 Ghostscript 的调用说明了这种情况:

swipl -g "number_puzzle_(_, Vs), label(Vs)" tour.pl | gs -g350x350 -r72 tour.ps -dNOPROMPT

比如我们在高亮位置看到很多回溯:

然而,本质问题已经完全存在于其他地方:

None 个突出显示的方块是有效的移动!

据此,我们看到您当前的公式没有——至少不够早——让求解器识别何时无法完成对解决方案的部分分配!这是坏消息,因为未能识别不一致的分配通常会导致不可接受的性能。例如,为了纠正 1 → 3 过渡(永远不会以这种方式发生,但已经是在这种情况下做出的第一个选择之一),求解器将不得不回溯大约 8 个正方形,在枚举之后——作为一个非常粗略的估计——258 = 152587890625 部分解,然后只从棋盘的第二个位置重新开始。

在约束文献中,这种回溯称为颠簸。这意味着重复失败由于同样的原因

这怎么可能?您的模型似乎是正确的,可用于检测解决方案。那挺好的!但是,一个好的约束公式不仅可以识别解决方案,还可以快速检测 部分分配无法 完成的解决方案。这就是允许求解器有效修剪搜索的原因,正是在这个重要方面,您当前的公式不足。原因之一与您正在使用的 具体约束 中的约束传播有关。特别是考虑以下查询:

?- (X + 1 #= 3) #<==> B, X #\= 2.

直觉上,我们期望 B = 0。但事实并非如此!相反,我们得到:

X in inf..1\/3..sup,
X+1#=_3840,
_3840#=3#B,
B in 0..1.

因此,求解器不会非常强烈地传播具体化的相等性。 也许应该!只有来自 Prolog 从业者的足够反馈 才能判断约束求解器的这个区域是否应该改变,可能会换一点速度为了更强的传播。此反馈的高度相关性是我建议只要有机会就使用 CLP(FD) 约束的原因之一,即每次推理 整数.[=49 时=]

对于这种特殊情况,我可以告诉您,从这个意义上讲,使求解器更强大并没有太大的区别。您基本上最终得到了一个核心问题仍然存在的电路板版本,其中有许多转换(其中一些在下面突出显示)在任何解决方案中都不会发生:

解决核心问题

我们应该从根本上消除回溯的原因。要修剪搜索,我们必须更早地识别不一致的(部分)分配。

直觉上,我们正在搜索连接的游览,并且希望在明确无法以预期方式继续游览时立即返回。

要完成我们想要的,我们至少有两个选择:

  1. 更改分配策略以考虑连通性
  2. 以更强烈地考虑连通性的方式对问题建模。

选项 1:分配策略

CLP(FD) 约束的一个主要吸引力在于它们让我们解耦 任务描述与搜索。在使用 CLP(FD) 约束时,我们经常通过 label/1labeling/2 执行搜索。然而,我们可以自由地以我们想要的任何方式给变量赋值。如果我们像您所做的那样遵循将 "constraint posting" 部分放入其自己的谓词(称为 核心关系.

的良好做法,这将非常容易

例如,这里有一个自定义分配策略,可确保游览始终保持连接状态:

allocation(Vs) :-
        length(Vs, N),
        numlist(1, N, Ns),
        maplist(member_(Vs), Ns).

member_(Es, E) :- member(E, Es).

通过这个策略,我们从头开始得到了 5×5 实例的解决方案:

?- number_puzzle_(Square, Vars), time(allocation(Vars)).
% 5,030,522 inferences, 0.907 CPU in 0.913 seconds (99% CPU, 5549133 Lips)
Square = square(row(1, 8, 5, 2, 11), ...),
Vars = [1, 8, 5, 2, 11, 16, 21, 24, 15|...] 

此策略有多种修改方式值得一试。例如,当多个正方形是可接受的时,我们可以尝试通过考虑正方形的剩余域元素的数量来做出更智能的选择。我将尝试此类改进作为一种挑战。

从标准的标签策略来看,min标签选项在这种情况下实际上与该策略非常相似,实际上它也找到了 5×5 情况的解决方案:

?- number_puzzle_(Square, Vars), time(labeling([min], Vars)).
% 22,461,798 inferences, 4.142 CPU in 4.174 seconds (99% CPU, 5422765 Lips)
Square = square(row(1, 8, 5, 2, 11), ...),
Vars = [1, 8, 5, 2, 11, 16, 21, 24, 15|...] .

然而,即使是合适的分配策略也不能完全补偿弱约束传播。对于 10×10 的实例,在使用 min 选项进行一些搜索后,棋盘看起来像这样:

请注意,我们还必须调整 PostScript 代码中 n 的值,以按预期将其可视化。

理想情况下,我们应该以从强传播中获益的方式来制定任务,然后使用好的分配策略。

选项 2:改造

一个好的 CLP 公式传播尽可能强烈(在可接受的时间内)。因此,我们应该努力使用允许求解器更容易推理任务最重要要求的约束。在这个具体案例中,这意味着我们应该尝试为当前表示为具体化约束的析取的内容找到更合适的公式,如上所示,不允许太多传播。理论上,约束求解器可以自动识别这种模式。然而,这对于许多用例来说是不切实际的,因此我们有时必须通过手动尝试几种有希望的公式来进行试验。尽管如此,在这种情况下:如果应用程序程序员有足够的反馈,这种情况更有可能得到改进和处理!

我现在使用 CLP(FD) 约束 circuit/1 来表明我们正在寻找 哈密顿电路在特定的图表中。该图表示为整数变量列表,其中每个元素表示其 后继者 在列表中的位置。

例如,一个包含 3 个元素的列表恰好包含 2 个哈密顿回路:

?- Vs = [_,_,_], circuit(Vs), label(Vs).
Vs = [2, 3, 1] ;
Vs = [3, 1, 2].

我使用 circuit/1 来描述也是 封闭式游览 的解决方案。这意味着,如果我们找到这样的解决方案,那么我们可以通过从找到的路线中的最后一个方格开始的有效移动从头开始:

n_tour(N, Vs) :-
        L #= N*N,
        length(Vs, L),
        successors(Vs, N, 1),
        circuit(Vs).

successors([], _, _).
successors([V|Vs], N, K0) :-
        findall(Num, n_k_next(N, K0, Num), [Next|Nexts]),
        foldl(num_to_dom, Nexts, Next, Dom),
        V in Dom,
        K1 #= K0 + 1,
        successors(Vs, N, K1).

num_to_dom(N, D0, D0\/N).

n_x_y_k(N, X, Y, K) :- [X,Y] ins 1..N, K #= N*(Y-1) + X.

n_k_next(N, K, Next) :-
        n_x_y_k(N, X0, Y0, K),
        (   [DX,DY] ins -2 \/ 2
        ;   [DX,DY] ins -3 \/ 0 \/ 3,
            abs(DX) + abs(DY) #= 3
        ),
        [X,Y] ins 1..N,
        X #= X0 + DX,
        Y #= Y0 + DY,
        n_x_y_k(N, X, Y, Next),
        label([DX,DY]).

请注意现在如何将可接受的后继者表示为域 元素 ,减少了约束的数量并完全消除了对具体化的需要。最重要的是,现在在搜索过程中的每个点都会自动考虑并强制执行预期的连通性。谓词 n_x_y_k/4(X,Y) 坐标与列表索引相关联。您可以通过更改 n_k_next/3 轻松地使该程序适应其他任务(例如,骑士之旅)。我把对 open 旅行的概括作为一个挑战。

这里有一些额外的定义,可以让我们以更易读的形式 打印 解决方案:

:- set_prolog_flag(double_quotes, chars).

print_tour(Vs) :-
        length(Vs, L),
        L #= N*N, N #> 0,
        length(Ts, N),
        tour_enumeration(Vs, N, Es),
        phrase(format_string(Ts, 0, 4), Fs),
        maplist(format(Fs), Es).

format_(Fs, Args, Xs0, Xs) :- format(chars(Xs0,Xs), Fs, Args).

format_string([], _, _) --> "\n".
format_string([_|Rest], N0, I) -->
        { N #= N0 + I },
        "~t~w~", call(format_("~w|", [N])),
        format_string(Rest, N, I).

tour_enumeration(Vs, N, Es) :-
        length(Es, N),
        maplist(same_length(Es), Es),
        append(Es, Ls),
        foldl(vs_enumeration(Vs, Ls), Vs, 1-1, _).

vs_enumeration(Vs, Ls, _, V0-E0, V-E) :-
        E #= E0 + 1,
        nth1(V0, Ls, E0),
        nth1(V0, Vs, V).

在具有强传播的公式中,预定义的 ff 搜索策略通常是一个很好的策略。事实上,它可以让我们在商用机器上几秒钟内解决整个任务,即原始的 10×10 实例:

?- n_tour(10, Vs),
   time(labeling([ff], Vs)),
   print_tour(Vs).
% 5,642,756 inferences, 0.988 CPU in 0.996 seconds (99% CPU, 5710827 Lips)
   1  96  15   2  97  14  80  98  13  79
  93  29  68  94  30  27  34  31  26  35
  16   3 100  17   4  99  12   9  81  45
  69  95  92  28  67  32  25  36  33  78
  84  18   5  83  19  10  82  46  11   8
  91  23  70  63  24  37  66  49  38  44
  72  88  20  73   6  47  76   7  41  77
  85  62  53  86  61  64  39  58  65  50
  90  22  71  89  21  74  42  48  75  43
  54  87  60  55  52  59  56  51  40  57
Vs = [4, 5, 21, 22, 8, 3, 29, 26, 6|...] 

为了获得最佳性能,我建议您也尝试使用其他 Prolog 系统。商业级 CLP(FD) 系统的效率通常是购买 Prolog 系统的重要原因。

另请注意,这绝不是任务的唯一有前途的 Prolog 甚至 CLP(FD) 公式,我将考虑其他公式作为挑战。