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

Question

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

给定一个 10x10 方格的网格，您的任务是只访问每个方格一次。在每个步骤中，您可以

• 水平或垂直跳过 2 个方块或
• 对角跳过 1 个正方形

换句话说（接近我执行以下）中，从1至100个这样的标记10×10网格数，在坐标每平方(X, Y)为1或等于一个大于“前”方在(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

我使用 CLP(FD) 在 SWI-Prolog 中的实现不能很好地扩展。事实上，除非预先指定了几乎两行，否则它甚至无法解决问题的 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).

Answer 1

首先：

这里发生了什么？

要查看发生了什么，以下是PostScript定义，可让我们将搜索可视化：

/n 5 定义

340 n div dup 刻度
-0.9 0.1 翻译 % 为线条笔划留出空间

/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 显示
    grestore } 绑定定义

/clr { gsave coords 1 setgray 0 0 1 1 4 copy rectfill
     0 setgray 0.02 setlinewidth rectstroke grestore} 绑定 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 数量
1 2 厘
2 3 4 数量

你得到：

PostScript 是一种用于可视化搜索过程的出色编程语言，我还建议您查看postscript以获取更多信息。

我们可以轻松修改您的程序以发出合适的 PostScript 指令，让我们直接观察搜索。我强调了相关的补充：

constrain_entries([], _Square, _X, _Y)。
constrain_entries([E|Es], Square, X, Y) :-
    冻结（E，后记（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使程序更加普遍。

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

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

例如，我们在突出显示的位置看到很多回溯：

然而，根本问题已经完全在别处：

突出显示的方块都不是有效的移动！

由此，我们看到您当前的公式没有——至少不是足够早——让求解器识别何时无法完成解决方案的部分分配！这是个坏消息，因为无法识别不一致的分配通常会导致不可接受的性能。例如，为了纠正 1 → 3 转换（这种转换永远不会发生，但在这种情况下已经是首选之一），求解器必须在枚举后回溯大约 8 个方格——如一个非常粗略的估计——25 ⁸  =  152587890625 个部分解，然后从棋盘中的第二个位置重新开始。

在约束文献中，这种回溯被称为thrashing。意思是同一个原因反复失败。

这怎么可能？您的模型似乎是正确的，可用于检测解决方案。那挺好的！然而，一个好的约束公式不仅可以识别解决方案，还可以快速检测到解决方案无法完成的部分分配。这使求解器能够有效地修剪搜索，而正是在这一重要方面，您当前的公式不足。其原因之一与您正在使用的具体化约束中的约束传播有关。特别地，请考虑以下查询：

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

直觉上，我们期望B = 0. 但事实并非如此！相反，我们得到：

X inf..1\/3..sup,
X+1#=_3840,
_3840#=3#B，
B 在 0..1。

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

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

修复核心问题

我们应该消除回溯事业的核心。为了修剪搜索，我们必须更早地识别不一致的（部分）分配。

直观地说，我们正在寻找连接的游览，并且希望在明确游览无法以预期方式继续时立即回溯。

为了完成我们想要的，我们至少有两个选择：

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

选项 1：分配策略

CLP(FD) 约束的一个主要吸引力在于它们让我们将任务描述与搜索分离。使用 CLP(FD) 约束时，我们经常通过label/1或执行搜索labeling/2。然而，我们可以以任何我们想要的方式自由地为变量赋值。如果我们遵循（正如您所做的那样）将“约束发布”部分放入其自己的谓词（称为核心关系）的良好实践，这将非常容易。

例如，这是一个自定义分配策略，可确保游览始终保持连接：

分配（Vs）：-
        长度（Vs，N），
        numlist(1, N, Ns),
        地图列表（成员_（Vs），Ns）。

member_(Es, E) :- 成员(E, Es)。

通过这个策略，我们从头开始得到一个 5×5 实例的解决方案：

?- number_puzzle_(Square, Vars), time( allocation(Vars) )。
% 5,030,522 推理，0.913 秒内 0.907 CPU（99% CPU，5549133 唇）
Square = square(row(1, 8, 5, 2, 11), ...),
变量 = [1, 8, 5, 2, 11, 16, 21, 24, 15|...] 

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

从标准的标签策略来看， min标签选项在这种情况下实际上与此策略非常相似，并且确实也为 5×5 情况找到了解决方案：

?- number_puzzle_(Square, Vars), time(labeling([ min ], Vars))。
% 22,461,798 推理，4.142 CPU 在 4.174 秒内（99% CPU，5422765 Lips）
Square = square(row(1, 8, 5, 2, 11), ...),
变量 = [1, 8, 5, 2, 11, 16, 21, 24, 15|...] 。

然而，即使是拟合分配策略也不能完全补偿弱约束传播。对于 10×10 实例，使用min选项进行一些搜索后，板子看起来像这样：

请注意，我们还必须调整nPostScript 代码中的 值以按预期进行可视化。

理想情况下，我们应制定任务，以这样的方式，我们从强大的传播中受益，然后也用好配置策略。

选项 2：改造

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

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

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

?- Vs = [_,_,_], 电路（Vs），标签（Vs）。
Vs = [2, 3, 1] ;
Vs = [3, 1, 2]。

我circuit/1用来描述也是封闭式旅行的解决方案。这意味着，如果我们找到这样的解决方案，那么我们可以通过从找到的巡回赛中最后一个方格的有效移动重新从头开始：

n_tour(N, Vs) :-
        L#= N*N,
        长度（Vs，L），
        后继者（Vs，N，1），
        电路（Vs）。

后继者（[]，_，_）。
后继者([V|Vs], N, K0) :-
        findall(Num, n_k_next(N, K0, Num), [Next|Nexts]),
        foldl(num_to_dom, Nexts, Next, Dom),
        V 在 Dom，
        K1 #= K0 + 1,
        后继者（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] 插入 1..N，
        X #= X0 + DX，
        Y#= Y0 + DY，
        n_x_y_k（N，X，Y，下一个），
        标签（[DX，DY]）。

请注意现在如何将可接受的后继表示为域元素，从而减少了约束的数量并完全消除了具体化的需要。最重要的是，现在在搜索过程中的每个点都会自动考虑并强制执行预期的连通性。谓词n_x_y_k/4将(X,Y)坐标与列表索引相关联。您可以通过更改n_k_next/3. 我离开推广到开放旅游是一个挑战。

以下是额外的定义，让我们以更易读的形式打印解决方案：

:- set_prolog_flag(double_quotes, chars)。

print_tour(Vs) :-
        长度（Vs，L），
        L#=N*N，N#>0，
        长度（Ts，N），
        Tour_enumeration(Vs, N, Es),
        短语（格式字符串（Ts，0, 4），Fs），
        地图列表（格式（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])),
        格式字符串（休息，N，我）。

Tour_enumeration(Vs, N, Es) :-
        长度（Es，N），
        地图列表（相同长度（Es），Es），
        追加（Es，Ls），
        foldl(vs_enumeration(Vs, Ls), Vs, 1-1, _)。

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

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

?- n_tour(10, Vs),
   时间（标签（[ff]，Vs）），
   print_tour(Vs)。
% 5,642,756 推理，0.996 秒内 0.988 CPU（99% CPU，5710827 唇）
   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) 公式，我将其他公式视为挑战。