Optimizing pathfinding in Constraint Logic Programming with Prolog
https://stackoverflow.com/questions/8458945
-
13-03-2021 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
I am working on a small prolog application to solve the Skyscrapers and Fences puzzle.
An unsolved puzzle:
A solved puzzle:
When I pass the program already solved puzzles it is quick, almost instantaneous, to validate it for me. When I pass the program really small puzzles (2x2, for example, with modified rules, of course), it is also quite fast to find a solution.
The problem is on computing puzzles with the "native" size of 6x6. I've left it running for 5 or so hours before aborting it. Way too much time.
I've found that the part that takes the longest is the "fences" one, not the "skyscrapers". Running "skyscrapers" separately results in a fast solution.
Here's my algorithm for fences:
- Vertices are represented by numbers, 0 means the path doesn't pass through that particular vertex, > 1 represents that vertex's order in the path.
- Constrain each cell to have the appropriate amount of lines surrounding it.
- That means that two vertexes are connected if they have sequential numbers, e.g., 1 -> 2, 2 -> 1, 1 ->
Max,Max-> 1 (Maxis the number for the last vertex in the path. computed viamaximum/2)
- That means that two vertexes are connected if they have sequential numbers, e.g., 1 -> 2, 2 -> 1, 1 ->
- Make sure each non-zero vertex has at least two neighboring vertices with sequential numbers
- Constrain
Maxto be equal to(BoardWidth + 1)^2 - NumberOfZeros(BoardWidth+1is the number of vertices along the edge andNumberOfZerosis computed viacount/4). - Use
nvalue(Vertices, Max + 1)to make sure the number of distinct values inVerticesisMax(i.e. the number of vertices in the path) plus1(zero values) - Find the first cell containing a
3and force the path to begin and end there, for efficiency purposes
What can I do to increase efficiency? Code is included below for reference.
skyscrapersinfences.pro
:-use_module(library(clpfd)).
:-use_module(library(lists)).
:-ensure_loaded('utils.pro').
:-ensure_loaded('s1.pro').
print_row([]).
print_row([Head|Tail]) :-
write(Head), write(' '),
print_row(Tail).
print_board(Board, BoardWidth) :-
print_board(Board, BoardWidth, 0).
print_board(_, BoardWidth, BoardWidth).
print_board(Board, BoardWidth, Index) :-
make_segment(Board, BoardWidth, Index, row, Row),
print_row(Row), nl,
NewIndex is Index + 1,
print_board(Board, BoardWidth, NewIndex).
print_boards([], _).
print_boards([Head|Tail], BoardWidth) :-
print_board(Head, BoardWidth), nl,
print_boards(Tail, BoardWidth).
get_board_element(Board, BoardWidth, X, Y, Element) :-
Index is BoardWidth*Y + X,
get_element_at(Board, Index, Element).
make_column([], _, _, []).
make_column(Board, BoardWidth, Index, Segment) :-
get_element_at(Board, Index, Element),
munch(Board, BoardWidth, MunchedBoard),
make_column(MunchedBoard, BoardWidth, Index, ColumnTail),
append([Element], ColumnTail, Segment).
make_segment(Board, BoardWidth, Index, row, Segment) :-
NIrrelevantElements is BoardWidth*Index,
munch(Board, NIrrelevantElements, MunchedBoard),
select_n_elements(MunchedBoard, BoardWidth, Segment).
make_segment(Board, BoardWidth, Index, column, Segment) :-
make_column(Board, BoardWidth, Index, Segment).
verify_segment(_, 0).
verify_segment(Segment, Value) :-
verify_segment(Segment, Value, 0).
verify_segment([], 0, _).
verify_segment([Head|Tail], Value, Max) :-
Head #> Max #<=> B,
Value #= M+B,
maximum(NewMax, [Head, Max]),
verify_segment(Tail, M, NewMax).
exactly(_, [], 0).
exactly(X, [Y|L], N) :-
X #= Y #<=> B,
N #= M +B,
exactly(X, L, M).
constrain_numbers(Vars) :-
exactly(3, Vars, 1),
exactly(2, Vars, 1),
exactly(1, Vars, 1).
iteration_values(BoardWidth, Index, row, 0, column) :-
Index is BoardWidth - 1.
iteration_values(BoardWidth, Index, Type, NewIndex, Type) :-
\+((Type = row, Index is BoardWidth - 1)),
NewIndex is Index + 1.
solve_skyscrapers(Board, BoardWidth) :-
solve_skyscrapers(Board, BoardWidth, 0, row).
solve_skyscrapers(_, BoardWidth, BoardWidth, column).
solve_skyscrapers(Board, BoardWidth, Index, Type) :-
make_segment(Board, BoardWidth, Index, Type, Segment),
domain(Segment, 0, 3),
constrain_numbers(Segment),
observer(Type, Index, forward, ForwardObserver),
verify_segment(Segment, ForwardObserver),
observer(Type, Index, reverse, ReverseObserver),
reverse(Segment, ReversedSegment),
verify_segment(ReversedSegment, ReverseObserver),
iteration_values(BoardWidth, Index, Type, NewIndex, NewType),
solve_skyscrapers(Board, BoardWidth, NewIndex, NewType).
build_vertex_list(_, Vertices, BoardWidth, X, Y, List) :-
V1X is X, V1Y is Y, V1Index is V1X + V1Y*(BoardWidth+1),
V2X is X+1, V2Y is Y, V2Index is V2X + V2Y*(BoardWidth+1),
V3X is X+1, V3Y is Y+1, V3Index is V3X + V3Y*(BoardWidth+1),
V4X is X, V4Y is Y+1, V4Index is V4X + V4Y*(BoardWidth+1),
get_element_at(Vertices, V1Index, V1),
get_element_at(Vertices, V2Index, V2),
get_element_at(Vertices, V3Index, V3),
get_element_at(Vertices, V4Index, V4),
List = [V1, V2, V3, V4].
build_neighbors_list(Vertices, VertexWidth, X, Y, [NorthMask, EastMask, SouthMask, WestMask], [NorthNeighbor, EastNeighbor, SouthNeighbor, WestNeighbor]) :-
NorthY is Y - 1,
EastX is X + 1,
SouthY is Y + 1,
WestX is X - 1,
NorthNeighborIndex is (NorthY)*VertexWidth + X,
EastNeighborIndex is Y*VertexWidth + EastX,
SouthNeighborIndex is (SouthY)*VertexWidth + X,
WestNeighborIndex is Y*VertexWidth + WestX,
(NorthY >= 0, get_element_at(Vertices, NorthNeighborIndex, NorthNeighbor) -> NorthMask = 1 ; NorthMask = 0),
(EastX < VertexWidth, get_element_at(Vertices, EastNeighborIndex, EastNeighbor) -> EastMask = 1 ; EastMask = 0),
(SouthY < VertexWidth, get_element_at(Vertices, SouthNeighborIndex, SouthNeighbor) -> SouthMask = 1 ; SouthMask = 0),
(WestX >= 0, get_element_at(Vertices, WestNeighborIndex, WestNeighbor) -> WestMask = 1 ; WestMask = 0).
solve_path(_, VertexWidth, 0, VertexWidth) :-
write('end'),nl.
solve_path(Vertices, VertexWidth, VertexWidth, Y) :-
write('switch row'),nl,
Y \= VertexWidth,
NewY is Y + 1,
solve_path(Vertices, VertexWidth, 0, NewY).
solve_path(Vertices, VertexWidth, X, Y) :-
X >= 0, X < VertexWidth, Y >= 0, Y < VertexWidth,
write('Path: '), nl,
write('Vertex width: '), write(VertexWidth), nl,
write('X: '), write(X), write(' Y: '), write(Y), nl,
VertexIndex is X + Y*VertexWidth,
write('1'),nl,
get_element_at(Vertices, VertexIndex, Vertex),
write('2'),nl,
build_neighbors_list(Vertices, VertexWidth, X, Y, [NorthMask, EastMask, SouthMask, WestMask], [NorthNeighbor, EastNeighbor, SouthNeighbor, WestNeighbor]),
L1 = [NorthMask, EastMask, SouthMask, WestMask],
L2 = [NorthNeighbor, EastNeighbor, SouthNeighbor, WestNeighbor],
write(L1),nl,
write(L2),nl,
write('3'),nl,
maximum(Max, Vertices),
write('4'),nl,
write('Max: '), write(Max),nl,
write('Vertex: '), write(Vertex),nl,
(Vertex #> 1 #/\ Vertex #\= Max) #=> (
((NorthMask #> 0 #/\ NorthNeighbor #> 0) #/\ (NorthNeighbor #= Vertex - 1)) #\
((EastMask #> 0 #/\ EastNeighbor #> 0) #/\ (EastNeighbor #= Vertex - 1)) #\
((SouthMask #> 0 #/\ SouthNeighbor #> 0) #/\ (SouthNeighbor #= Vertex - 1)) #\
((WestMask #> 0 #/\ WestNeighbor #> 0) #/\ (WestNeighbor #= Vertex - 1))
) #/\ (
((NorthMask #> 0 #/\ NorthNeighbor #> 0) #/\ (NorthNeighbor #= Vertex + 1)) #\
((EastMask #> 0 #/\ EastNeighbor #> 0) #/\ (EastNeighbor #= Vertex + 1)) #\
((SouthMask #> 0 #/\ SouthNeighbor #> 0) #/\ (SouthNeighbor #= Vertex + 1)) #\
((WestMask #> 0 #/\ WestNeighbor #> 0) #/\ (WestNeighbor #= Vertex + 1))
),
write('5'),nl,
Vertex #= 1 #=> (
((NorthMask #> 0 #/\ NorthNeighbor #> 0) #/\ (NorthNeighbor #= Max)) #\
((EastMask #> 0 #/\ EastNeighbor #> 0) #/\ (EastNeighbor #= Max)) #\
((SouthMask #> 0 #/\ SouthNeighbor #> 0) #/\ (SouthNeighbor #= Max)) #\
((WestMask #> 0 #/\ WestNeighbor #> 0) #/\ (WestNeighbor #= Max))
) #/\ (
((NorthMask #> 0 #/\ NorthNeighbor #> 0) #/\ (NorthNeighbor #= 2)) #\
((EastMask #> 0 #/\ EastNeighbor #> 0) #/\ (EastNeighbor #= 2)) #\
((SouthMask #> 0 #/\ SouthNeighbor #> 0) #/\ (SouthNeighbor #= 2)) #\
((WestMask #> 0 #/\ WestNeighbor #> 0) #/\ (WestNeighbor #= 2))
),
write('6'),nl,
Vertex #= Max #=> (
((NorthMask #> 0 #/\ NorthNeighbor #> 0) #/\ (NorthNeighbor #= 1)) #\
((EastMask #> 0 #/\ EastNeighbor #> 0) #/\ (EastNeighbor #= 1)) #\
((SouthMask #> 0 #/\ SouthNeighbor #> 0) #/\ (SouthNeighbor #= 1)) #\
((WestMask #> 0 #/\ WestNeighbor #> 0) #/\ (WestNeighbor #= 1))
) #/\ (
((NorthMask #> 0 #/\ NorthNeighbor #> 0) #/\ (NorthNeighbor #= Max - 1)) #\
((EastMask #> 0 #/\ EastNeighbor #> 0) #/\ (EastNeighbor #= Max - 1)) #\
((SouthMask #> 0 #/\ SouthNeighbor #> 0) #/\ (SouthNeighbor #= Max - 1)) #\
((WestMask #> 0 #/\ WestNeighbor #> 0) #/\ (WestNeighbor #= Max - 1))
),
write('7'),nl,
NewX is X + 1,
solve_path(Vertices, VertexWidth, NewX, Y).
solve_fences(Board, Vertices, BoardWidth) :-
VertexWidth is BoardWidth + 1,
write('- Solving vertices'),nl,
solve_vertices(Board, Vertices, BoardWidth, 0, 0),
write('- Solving path'),nl,
solve_path(Vertices, VertexWidth, 0, 0).
solve_vertices(_, _, BoardWidth, 0, BoardWidth).
solve_vertices(Board, Vertices, BoardWidth, BoardWidth, Y) :-
Y \= BoardWidth,
NewY is Y + 1,
solve_vertices(Board, Vertices, BoardWidth, 0, NewY).
solve_vertices(Board, Vertices, BoardWidth, X, Y) :-
X >= 0, X < BoardWidth, Y >= 0, Y < BoardWidth,
write('process'),nl,
write('X: '), write(X), write(' Y: '), write(Y), nl,
build_vertex_list(Board, Vertices, BoardWidth, X, Y, [V1, V2, V3, V4]),
write('1'),nl,
get_board_element(Board, BoardWidth, X, Y, Element),
write('2'),nl,
maximum(Max, Vertices),
(V1 #> 0 #/\ V2 #> 0 #/\
(
(V1 + 1 #= V2) #\
(V1 - 1 #= V2) #\
(V1 #= Max #/\ V2 #= 1) #\
(V1 #= 1 #/\ V2 #= Max)
)
) #<=> B1,
(V2 #> 0 #/\ V3 #> 0 #/\
(
(V2 + 1 #= V3) #\
(V2 - 1 #= V3) #\
(V2 #= Max #/\ V3 #= 1) #\
(V2 #= 1 #/\ V3 #= Max)
)
) #<=> B2,
(V3 #> 0 #/\ V4 #> 0 #/\
(
(V3 + 1 #= V4) #\
(V3 - 1 #= V4) #\
(V3 #= Max #/\ V4 #= 1) #\
(V3 #= 1 #/\ V4 #= Max)
)
) #<=> B3,
(V4 #> 0 #/\ V1 #> 0 #/\
(
(V4 + 1 #= V1) #\
(V4 - 1 #= V1) #\
(V4 #= Max #/\ V1 #= 1) #\
(V4 #= 1 #/\ V1 #= Max)
)
) #<=> B4,
write('3'),nl,
sum([B1, B2, B3, B4], #= , C),
write('4'),nl,
Element #> 0 #=> C #= Element,
write('5'),nl,
NewX is X + 1,
solve_vertices(Board, Vertices, BoardWidth, NewX, Y).
sel_next_variable_for_path(Vars,Sel,Rest) :-
% write(Vars), nl,
findall(Idx-Cost, (nth1(Idx, Vars,V), fd_set(V,S), fdset_size(S,Size), fdset_min(S,Min), var_cost(Min,Size, Cost)), L),
min_member(comp, BestIdx-_MinCost, L),
nth1(BestIdx, Vars, Sel, Rest),!.
var_cost(0, _, 1000000) :- !.
var_cost(_, 1, 1000000) :- !.
var_cost(X, _, X).
%build_vertex_list(_, Vertices, BoardWidth, X, Y, List)
constrain_starting_and_ending_vertices(Vertices, [V1,V2,V3,V4]) :-
maximum(Max, Vertices),
(V1 #= 1 #/\ V2 #= Max #/\ V3 #= Max - 1 #/\ V4 #= 2 ) #\
(V1 #= Max #/\ V2 #= 1 #/\ V3 #= 2 #/\ V4 #= Max - 1 ) #\
(V1 #= Max - 1 #/\ V2 #= Max #/\ V3 #= 1 #/\ V4 #= 2 ) #\
(V1 #= 2 #/\ V2 #= 1 #/\ V3 #= Max #/\ V4 #= Max - 1 ) #\
(V1 #= 1 #/\ V2 #= 2 #/\ V3 #= Max - 1 #/\ V4 #= Max ) #\
(V1 #= Max #/\ V2 #= Max - 1 #/\ V3 #= 2 #/\ V4 #= 1 ) #\
(V1 #= Max - 1 #/\ V2 #= 2 #/\ V3 #= 1 #/\ V4 #= Max ) #\
(V1 #= 2 #/\ V2 #= Max - 1 #/\ V3 #= Max #/\ V4 #= 1 ).
set_starting_and_ending_vertices(Board, Vertices, BoardWidth) :-
set_starting_and_ending_vertices(Board, Vertices, BoardWidth, 0, 0).
set_starting_and_ending_vertices(Board, Vertices, BoardWidth, BoardWidth, Y) :-
Y \= BoardWidth,
NewY is Y + 1,
solve_path(Board, Vertices, BoardWidth, 0, NewY).
set_starting_and_ending_vertices(Board, Vertices, BoardWidth, X, Y) :-
X >= 0, X < BoardWidth, Y >= 0, Y < BoardWidth,
build_vertex_list(_, Vertices, BoardWidth, X, Y, List),
get_board_element(Board, BoardWidth, X, Y, Element),
(Element = 3 ->
constrain_starting_and_ending_vertices(Vertices, List)
;
NewX is X + 1,
set_starting_and_ending_vertices(Board, Vertices, BoardWidth, NewX, Y)).
solve(Board, Vertices, BoardWidth) :-
write('Skyscrapers'), nl,
solve_skyscrapers(Board, BoardWidth),
write('Labeling'), nl,
labeling([ff], Board), !,
write('Setting domain'), nl,
NVertices is (BoardWidth+1)*(BoardWidth+1),
domain(Vertices, 0, NVertices),
write('Starting and ending vertices'), nl,
set_starting_and_ending_vertices(Board, Vertices, BoardWidth),
write('Setting maximum'), nl,
maximum(Max, Vertices),
write('1'),nl,
Max #> BoardWidth + 1,
write('2'),nl,
Max #< NVertices,
count(0, Vertices, #=, NZeros),
Max #= NVertices - NZeros,
write('3'),nl,
write('Calling nvalue'), nl,
ValueCount #= Max + 1,
nvalue(ValueCount, Vertices),
write('Solving fences'), nl,
solve_fences(Board, Vertices, BoardWidth),
write('Labeling'), nl,
labeling([ff], Vertices).
main :-
board(Board),
board_width(BoardWidth),
vertices(Vertices),
solve(Board, Vertices, BoardWidth),
%findall(Board,
% labeling([ff], Board),
% Boards
%),
%append(Board, Vertices, Final),
write('done.'),nl,
print_board(Board, 6), nl,
print_board(Vertices, 7).
utils.pro
get_element_at([Head|_], 0, Head).
get_element_at([_|Tail], Index, Element) :-
Index \= 0,
NewIndex is Index - 1,
get_element_at(Tail, NewIndex, Element).
reverse([], []).
reverse([Head|Tail], Inv) :-
reverse(Tail, Aux),
append(Aux, [Head], Inv).
munch(List, 0, List).
munch([_|Tail], Count, FinalList) :-
Count > 0,
NewCount is Count - 1,
munch(Tail, NewCount, FinalList).
select_n_elements(_, 0, []).
select_n_elements([Head|Tail], Count, FinalList) :-
Count > 0,
NewCount is Count - 1,
select_n_elements(Tail, NewCount, Result),
append([Head], Result, FinalList).
generate_list(Element, NElements, [Element|Result]) :-
NElements > 0,
NewNElements is NElements - 1,
generate_list(Element, NewNElements, Result).
generate_list(_, 0, []).
s1.pro
% Skyscrapers and Fences puzzle S1
board_width(6).
%observer(Type, Index, Orientation, Observer),
observer(row, 0, forward, 2).
observer(row, 1, forward, 2).
observer(row, 2, forward, 2).
observer(row, 3, forward, 1).
observer(row, 4, forward, 2).
observer(row, 5, forward, 1).
observer(row, 0, reverse, 1).
observer(row, 1, reverse, 1).
observer(row, 2, reverse, 2).
observer(row, 3, reverse, 3).
observer(row, 4, reverse, 2).
observer(row, 5, reverse, 2).
observer(column, 0, forward, 2).
observer(column, 1, forward, 3).
observer(column, 2, forward, 0).
observer(column, 3, forward, 2).
observer(column, 4, forward, 2).
observer(column, 5, forward, 1).
observer(column, 0, reverse, 1).
observer(column, 1, reverse, 1).
observer(column, 2, reverse, 2).
observer(column, 3, reverse, 2).
observer(column, 4, reverse, 2).
observer(column, 5, reverse, 2).
board(
[
_, _, 2, _, _, _,
_, _, _, _, _, _,
_, 2, _, _, _, _,
_, _, _, 2, _, _,
_, _, _, _, _, _,
_, _, _, _, _, _
]
).
vertices(
[
_, _, _, _, _, _, _,
_, _, _, _, _, _, _,
_, _, _, _, _, _, _,
_, _, _, _, _, _, _,
_, _, _, _, _, _, _,
_, _, _, _, _, _, _,
_, _, _, _, _, _, _
]
).
La solution
I also, like twinterer, enjoyed this puzzle. But being a principiant, I had first to discover an appropriate strategy, for both skyscrapes and fences part, and then deeply debugging the latter, cause a copy variables problem that locked me many hours.
Once solved the bug, I faced the inefficiency of my first attempt. I reworked in plain Prolog a similar schema, just to verify how inefficient it was.
At least, I understood how use CLP(FD) more effectively to model the problem (with help from the twinterer' answer), and now the program is fast (0,2 sec). So now I can hint you about your code: the required constraints are far simpler than those you coded: for the fences part, i.e. with a buildings placement fixed, we have 2 constraints: number of edges where height > 0, and linking the edges together: when an edge is used, the sum of adjacents must be 1 (on both sides).
Here is the last version of my code, developed with SWI-Prolog.
/* File: skys.pl
Author: Carlo,,,
Created: Dec 11 2011
Purpose: questions/8458945 on http://stackoverflow.com
http://stackoverflow.com/questions/8458945/optimizing-pathfinding-in-constraint-logic-programming-with-prolog
*/
:- module(skys, [skys/0, fences/2, draw_path/2]).
:- [index_square,
lambda,
library(clpfd),
library(aggregate)].
puzzle(1,
[[-,2,3,-,2,2,1,-],
[2,-,-,2,-,-,-,1],
[2,-,-,-,-,-,-,1],
[2,-,2,-,-,-,-,2],
[1,-,-,-,2,-,-,3],
[2,-,-,-,-,-,-,2],
[1,-,-,-,-,-,-,2],
[-,1,1,2,2,2,2,-]]).
skys :-
puzzle(1, P),
skyscrapes(P, Rows),
flatten(Rows, Flat),
label(Flat),
maplist(writeln, Rows),
fences(Rows, Loop),
writeln(Loop),
draw_path(7, Loop).
%% %%%%%%%%%%
% skyscrapes part
% %%%%%%%%%%
skyscrapes(Puzzle, Rows) :-
% massaging definition: separe external 'visibility' counters
first_and_last(Puzzle, Fpt, Lpt, Wpt),
first_and_last(Fpt, -, -, Fp),
first_and_last(Lpt, -, -, Lp),
maplist(first_and_last, Wpt, Lc, Rc, InnerData),
% InnerData it's the actual 'playground', Fp, Lp, Lc, Rc are list of counters
maplist(make_vars, InnerData, Rows),
% exploit symmetry wrt rows/cols
transpose(Rows, Cols),
% each row or col contains once 1,2,3
Occurs = [0-_, 1-1, 2-1, 3-1], % allows any grid size leaving unspecified 0s
maplist(\Vs^global_cardinality(Vs, Occurs), Rows),
maplist(\Vs^global_cardinality(Vs, Occurs), Cols),
% apply 'external visibility' constraint
constraint_views(Lc, Rows),
constraint_views(Fp, Cols),
maplist(reverse, Rows, RRows),
constraint_views(Rc, RRows),
maplist(reverse, Cols, RCols),
constraint_views(Lp, RCols).
first_and_last(List, First, Last, Without) :-
append([[First], Without, [Last]], List).
make_vars(Data, Vars) :-
maplist(\C^V^(C \= (-) -> V #= C ; V in 0..3), Data, Vars).
constraint_views(Ns, Ls) :-
maplist(\N^L^
( N \= (-)
-> constraint_view(0, L, Rs),
sum(Rs, #=, N)
; true
), Ns, Ls).
constraint_view(_, [], []).
constraint_view(Top, [V|Vs], [R|Rs]) :-
R #<==> V #> 0 #/\ V #> Top,
Max #= max(Top, V),
constraint_view(Max, Vs, Rs).
%% %%%%%%%%%%%%%%%
% fences part
% %%%%%%%%%%%%%%%
fences(SkyS, Ps) :-
length(SkyS, D),
% allocate edges
max_dimensions(D, _,_,_,_, N),
N1 is N + 1,
length(Edges, N1),
Edges ins 0..1,
findall((R, C, V),
(nth0(R, SkyS, Row), nth0(C, Row, V), V > 0),
Buildings),
maplist(count_edges(D, Edges), Buildings),
findall((I, Adj1, Adj2),
(between(0, N, I), edge_adjacents(D, I, Adj1, Adj2)),
Path),
maplist(make_path(Edges), Path, Vs),
flatten([Edges, Vs], Gs),
label(Gs),
used_edges_to_path_coords(D, Edges, Ps).
count_edges(D, Edges, (R, C, V)) :-
cell_edges(D, (R, C), Is),
idxs0_to_elems(Is, Edges, Es),
sum(Es, #=, V).
make_path(Edges, (Index, G1, G2), [S1, S2]) :-
idxs0_to_elems(G1, Edges, Adj1),
idxs0_to_elems(G2, Edges, Adj2),
nth0(Index, Edges, Edge),
[S1, S2] ins 0..3,
sum(Adj1, #=, S1),
sum(Adj2, #=, S2),
Edge #= 1 #<==> S1 #= 1 #/\ S2 #= 1.
%% %%%%%%%%%%%%%%
% utility: draw a path with arrows
% %%%%%%%%%%%%%%
draw_path(D, P) :-
forall(between(1, D, R),
( forall(between(1, D, C),
( V is (R - 1) * D + C - 1,
U is (R - 2) * D + C - 1,
( append(_, [V, U|_], P)
-> write(' ^ ')
; append(_, [U, V|_], P)
-> write(' v ')
; write(' ')
)
)),
nl,
forall(between(1, D, C),
( V is (R - 1) * D + C - 1,
( V < 10
-> write(' ') ; true
),
write(V),
U is V + 1,
( append(_, [V, U|_], P)
-> write(' > ')
; append(_, [U, V|_], P)
-> write(' < ')
; write(' ')
)
)),
nl
)
).
% convert from 'edge used flags' to vertex indexes
%
used_edges_to_path_coords(D, EdgeUsedFlags, PathCoords) :-
findall((X, Y),
(nth0(Used, EdgeUsedFlags, 1), edge_verts(D, Used, X, Y)),
Path),
Path = [(First, _)|_],
edge_follower(First, Path, PathCoords).
edge_follower(C, Path, [C|Rest]) :-
( select(E, Path, Path1),
( E = (C, D) ; E = (D, C) )
-> edge_follower(D, Path1, Rest)
; Rest = []
).
The output:
[0,0,2,1,0,3]
[2,1,3,0,0,0]
[0,2,0,3,1,0]
[0,3,0,2,0,1]
[1,0,0,0,3,2]
[3,0,1,0,2,0]
[1,2,3,4,5,6,13,12,19,20,27,34,41,48,47,40,33,32,39,46,45,38,31,24,25,18,17,10,9,16,23,
22,29,30,37,36,43,42,35,28,21,14,7,8,1]
0 1 > 2 > 3 > 4 > 5 > 6
^ v
7 > 8 9 < 10 11 12 < 13
^ v ^ v
14 15 16 17 < 18 19 > 20
^ v ^ v
21 22 < 23 24 > 25 26 27
^ v ^ v
28 29 > 30 31 32 < 33 34
^ v ^ v ^ v
35 36 < 37 38 39 40 41
^ v ^ v ^ v
42 < 43 44 45 < 46 47 < 48
As I mentioned, my first attempt was more 'procedural': it draws a loop, but the problem I was unable to solve is basically that the cardinality of vertices subset must be known before, being based on the global constraint all_different. It painfully works on a reduced 4*4 puzzle, but I stopped it after some hours on the 6*6 original. Anyway, learning from scratch how to draw a path with CLP(FD) has been rewarding.
t :-
time(fences([[0,0,2,1,0,3],
[2,1,3,0,0,0],
[0,2,0,3,1,0],
[0,3,0,2,0,1],
[1,0,0,0,3,2],
[3,0,1,0,2,0]
],L)),
writeln(L).
fences(SkyS, Ps) :-
length(SkyS, Dt),
D is Dt + 1,
Sq is D * D - 1,
% min/max num. of vertices
aggregate_all(sum(V), (member(R, SkyS), member(V, R)), MinVertsT),
MinVerts is max(4, MinVertsT),
MaxVerts is D * D,
% find first cell with heigth 3, for sure start vertex
nth0(R, SkyS, Row), nth0(C, Row, 3),
% search a path with at least MinVerts
between(MinVerts, MaxVerts, NVerts),
length(Vs, NVerts),
Vs ins 0 .. Sq,
all_distinct(Vs),
% make a loop
Vs = [O|_],
O is R * D + C,
append(Vs, [O], Ps),
% apply #edges check
findall(rc(Ri, Ci, V),
(nth0(Ri, SkyS, Rowi),
nth0(Ci, Rowi, V),
V > 0), VRCs),
maplist(count_edges(Ps, D), VRCs),
connect_path(D, Ps),
label(Vs).
count_edges(Ps, D, rc(R, C, V)) :-
V0 is R * D + C,
V1 is R * D + C + 1,
V2 is (R + 1) * D + C,
V3 is (R + 1) * D + C + 1,
place_edges(Ps, [V0-V1, V0-V2, V1-V3, V2-V3], Ts),
flatten(Ts, Tsf),
sum(Tsf, #=, V).
place_edges([A,B|Ps], L, [R|Rs]) :-
place_edge(L, A-B, R),
place_edges([B|Ps], L, Rs).
place_edges([_], _L, []).
place_edge([M-N | L], A-B, [Y|R]) :-
Y #<==> (A #= M #/\ B #= N) #\/ (A #= N #/\ B #= M),
place_edge(L, A-B, R).
place_edge([], _, []).
connect(X, D, Y) :-
D1 is D - 1,
[R, C] ins 0 .. D1,
X #= R * D + C,
( C #< D - 1, Y #= R * D + C + 1
; R #< D - 1, Y #= (R + 1) * D + C
; C #> 0, Y #= R * D + C - 1
; R #> 0, Y #= (R - 1) * D + C
).
connect_path(D, [X, Y | R]) :-
connect(X, D, Y),
connect_path(D, [Y | R]).
connect_path(_, [_]).
Thanks you for such interesting question.
MORE EDIT:here the main miss piece of code for the complete solution (index_square.pl)
/* File: index_square.pl
Author: Carlo,,,
Created: Dec 15 2011
Purpose: indexing square grid for FD mapping
*/
:- module(index_square,
[max_dimensions/6,
idxs0_to_elems/3,
edge_verts/4,
edge_is_horiz/3,
cell_verts/3,
cell_edges/3,
edge_adjacents/4,
edge_verts_all/2
]).
%
% index row : {D}, left to right
% index col : {D}, top to bottom
% index cell : same as top edge or row,col
% index vert : {(D + 1) * 2}
% index edge : {(D * (D + 1)) * 2}, first all horiz, then vert
%
% {N} denote range 0 .. N-1
%
% on a 2*2 grid, the numbering schema is
%
% 0 1
% 0-- 0 --1-- 1 --2
% | | |
% 0 6 0,0 7 0,1 8
% | | |
% 3-- 2 --4-- 3 --5
% | | |
% 1 9 1,0 10 1,1 11
% | | |
% 6-- 4 --7-- 5 --8
%
% while on a 4*4 grid:
%
% 0 1 2 3
% 0-- 0 --1-- 1 --2-- 2 --3-- 3 --4
% | | | | |
% 0 20 21 22 23 24
% | | | | |
% 5-- 4 --6-- 5 --7-- 6 --8-- 7 --9
% | | | | |
% 1 25 26 27 28 29
% | | | | |
% 10--8 --11- 9 --12--10--13--11--14
% | | | | |
% 2 30 31 32 33 34
% | | | | |
% 15--12--16--13--17--14--18--15--19
% | | | | |
% 3 35 36 37 38 39
% | | | | |
% 20--16--21--17--22--18--23--19--24
%
% | |
% --+-- N --+--
% | |
% W R,C E
% | |
% --+-- S --+--
% | |
%
% get range upper value for interesting quantities
%
max_dimensions(D, MaxRow, MaxCol, MaxCell, MaxVert, MaxEdge) :-
MaxRow is D - 1,
MaxCol is D - 1,
MaxCell is D * D - 1,
MaxVert is ((D + 1) * 2) - 1,
MaxEdge is (D * (D + 1) * 2) - 1.
% map indexes to elements
%
idxs0_to_elems(Is, Edges, Es) :-
maplist(nth0_(Edges), Is, Es).
nth0_(Edges, I, E) :-
nth0(I, Edges, E).
% get vertices of edge
%
edge_verts(D, E, X, Y) :-
S is D + 1,
edge_is_horiz(D, E, H),
( H
-> X is (E // D) * S + E mod D,
Y is X + 1
; X is E - (D * S),
Y is X + S
).
% qualify edge as horizontal (never fail!)
%
edge_is_horiz(D, E, H) :-
E >= (D * (D + 1)) -> H = false ; H = true.
% get 4 vertices of cell
%
cell_verts(D, (R, C), [TL, TR, BL, BR]) :-
TL is R * (D + 1) + C,
TR is TL + 1,
BL is TR + D,
BR is BL + 1.
% get 4 edges of cell
%
cell_edges(D, (R, C), [N, S, W, E]) :-
N is R * D + C,
S is N + D,
W is (D * (D + 1)) + R * (D + 1) + C,
E is W + 1.
% get adjacents at two extremities of edge I
%
edge_adjacents(D, I, G1, G2) :-
edge_verts(D, I, X, Y),
edge_verts_all(D, EVs),
setof(E, U^V^(member(E - (U, V), EVs), E \= I, (U == X ; V == X)), G1),
setof(E, U^V^(member(E - (U, V), EVs), E \= I, (U == Y ; V == Y)), G2).
% get all edge_verts/4 for grid D
%
edge_verts_all(D, L) :-
( edge_verts_all_(D, L)
-> true
; max_dimensions(D, _,_,_,_, S), %S is (D + 1) * (D + 2) - 1,
findall(E - (X, Y),
( between(0, S, E),
edge_verts(D, E, X, Y)
), L),
assert(edge_verts_all_(D, L))
).
:- dynamic edge_verts_all_/2.
%% %%%%%%%%%%%%%%%%%%%%
:- begin_tests(index_square).
test(1) :-
cell_edges(2, (0,1), [1, 3, 7, 8]),
cell_edges(2, (1,1), [3, 5, 10, 11]).
test(2) :-
cell_verts(2, (0,1), [1, 2, 4, 5]),
cell_verts(2, (1,1), [4, 5, 7, 8]).
test(3) :-
edge_is_horiz(2, 0, true),
edge_is_horiz(2, 5, true),
edge_is_horiz(2, 6, false),
edge_is_horiz(2, 9, false),
edge_is_horiz(2, 11, false).
test(4) :-
edge_verts(2, 0, 0, 1),
edge_verts(2, 3, 4, 5),
edge_verts(2, 5, 7, 8),
edge_verts(2, 6, 0, 3),
edge_verts(2, 11, 5, 8).
test(5) :-
edge_adjacents(2, 0, A, B), A = [6], B = [1, 7],
edge_adjacents(2, 9, [2, 6], [4]),
edge_adjacents(2, 10, [2, 3, 7], [4, 5]).
test(6) :-
cell_edges(4, (2,1), [9, 13, 31, 32]).
:- end_tests(index_square).
Autres conseils
A quick glance over your program suggests that you use reification quite heavily. Unfortunately, such formulations imply weak consistency in current systems like SICStus.
Often, however, things can be formulated more compactly leading to better consistency. Here is one example which you might adapt to your needs.
Say, you want to express that (X1,Y1) and (X2,Y2) are horizontal or vertical neighbors. You could say ( X1+1 #= X2 #/\ Y1 #= Y2 ) #\ ... for each possiblity (and check if your health insurance covers RSI).
Or you can say abs(X1-X2)+abs(Y1-Y2) #= 1. In the olden tymes SICStus Prolog used to have a symmetric difference (--)/2 for that, but I assume you are using version 4.
Above formulation maintains interval consistency (at least I conclude this from the examples I tried):
| ?- X1 in 1..9, abs(X1-X2)+abs(Y1-Y2) #= 1.
X1 in 1..9,
X2 in 0..10,...
So the X2 is readily constrained!
There might be situations (as you indicate in your response) where you need the reified form to maintain other constraints. In this case you might consider to post both.
Leaf through the manual, there are several combinatorial constraints that might be interesting too. And as a quick fix, smt/1 might help (new in 4.2.0). Would be interested to hear about this...
Another possibility might be to use another implementation: For example library(clpfd) of YAP or SWI.
What a nice little puzzle! In order to understand the properties, I implemented a solution in ECLiPSe. It can be found here: http://pastebin.com/eZbgjgFA (Don't worry if you see loops in the code: these can be easily translated into standard Prolog predicates. There's other stuff, though, that's not so easily translated from ECLiPSe into Sicstus)
Execution time is faster than what you report, but could probably be better:
?- snf(L).
L = [[]([]([](0,0,1,1),[](1,1,0,0),[](0,1,0,1),[](0,1,0,0),[](0,1,0,0),[](0,1,1,1)),
[]([](1,1,0,0),[](0,0,1,0),[](1,1,1,0),[](1,0,0,1),[](0,0,1,0),[](1,1,0,1)),
[]([](1,0,0,0),[](0,0,1,1),[](1,0,0,0),[](0,1,1,1),[](1,0,0,0),[](0,1,1,0)),
[]([](1,0,1,0),[](1,1,0,1),[](0,0,1,0),[](1,1,0,0),[](0,0,0,1),[](0,0,1,0)),
[]([](1,0,0,0),[](0,1,1,1),[](1,0,1,0),[](1,0,1,0),[](1,1,1,0),[](1,0,1,0)),
[]([](1,0,1,1),[](1,1,0,0),[](0,0,1,0),[](1,0,1,1),[](1,0,1,0),[](1,0,1,1))),
...]
Yes (40.42s cpu, solution 1, maybe more)
No (52.88s cpu)
What you see in the answer is the matrix of edges. Each inner term denotes for a field in the puzzle, which edge is active (left,up,right,down). I edited out the rest.
I used eight arrays in total: the HxWx4 array of edges (0/1), a (H+1)x(W+1) array of active edges per field vertex (0/2), a HxW array of sums of active edges (0..3), a HxW array of buildings (0/1), two [H,W]x3 arrays of building heights, and two [H,W]x3 arrays of building positions.
The requirement that there must be only one path is not put as a constraint, but merely executed as a check after a potential solution is found during labeling.
The constraints are:
the sum array must contain for each field the sum of the active edges for that field
touching edges of neighbouring fields must contain the same value
the vertex points must have two active edges connected to them, or none
in each column/row, exactly three buildings must be placed. Some of the buildings are placed by the definition of the puzzle
each building height in a row/column must be different
the building height corresponds to the sum of active edges at this position
the number of visible buildings is specified by the definition of the puzzle. This restricts the order in which the buildings can appear in the row/column.
the positions of the buildings in a row/column must be given in ascending order
once the position of the first/second/third building is known, we can infer some positions where no building can be placed.
With this set of constraints, we are now ready to label. Labeling is done in two steps, which speeds up the solving process.
In a first step, only the building positions are labeled. This is the most restricted part, and if we find a solution here, the rest is much easier.
In the second step, all other variables are labeled. For both steps, I chose "first fail" as labeling strategy, i.e., label variables with smallest domain first.
Without solving the building positions first, the program takes much longer (I always stopped it after some minutes). Since I didn't have a second puzzle instance available, I'm not sure that the search strategy will feasible in all instances, though
Looking through your program again, it seems that you follow a similar strategy of placing the buildings first. However, you iterate between setting constraints and labeling. This is not efficient. In CLP, you should always place the constraints upfront (unless the constraints really depend on the current state of the partial solution), and only when the constraints are posted you search for a solution. This way, you can detect failure regarding all constraints during your search. Otherwise, you may find a partial solution that fulfils the set of constraints that you posted so far, only to find out that you cannot complete the solution once the other constraints are added.
Also, if you have different sets of variables, experiment with the order in which the variables are labeled. There's no universal recipe for that, though.
Hope this helps!
