Optimizing pathfinding in Constraint Logic Programming with Prolog
Asked Answered
B

3

11

I am working on a small prolog application to solve the Skyscrapers and Fences puzzle.

An unsolved puzzle:

Skyscrapers in fences puzzle (unsolved)

A solved puzzle:

Skyscrapers in fences puzzle (solved)

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 (Max is the number for the last vertex in the path. computed via maximum/2)
  • Make sure each non-zero vertex has at least two neighboring vertices with sequential numbers
  • Constrain Max to be equal to (BoardWidth + 1)^2 - NumberOfZeros (BoardWidth+1 is the number of vertices along the edge and NumberOfZeros is computed via count/4).
  • Use nvalue(Vertices, Max + 1) to make sure the number of distinct values in Vertices is Max (i.e. the number of vertices in the path) plus 1 (zero values)
  • Find the first cell containing a 3 and 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(
    [
        _, _, _, _, _, _, _,
        _, _, _, _, _, _, _,
        _, _, _, _, _, _, _,
        _, _, _, _, _, _, _,
        _, _, _, _, _, _, _,
        _, _, _, _, _, _, _,
        _, _, _, _, _, _, _
    ]
).
Beale answered 10/12, 2011 at 18:44 Comment(0)
H
13

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
        https://mcmap.net/q/969314/-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).
Hockey answered 16/12, 2011 at 6:6 Comment(20)
This is very, very nice! Thank you for working on this and posting your solution. I wonder if CLP(FD)'s automaton/3 or automaton/8 constraint might also be useful to describe a path in this problem.Tonie
I wonder also! I tried to understand it, but the lacks of examples or explanations has put me off. The underling mechanics it's so complex, but I'll try again!Hockey
And: The definition of idxs0_to_elems/3 is missing (presumably it is in the file "index_square"). And: In your initial attempt, it seems that you can use the circuit/1 constraint from library(clpfd) to describe the loop.Tonie
When you publish the full code, will you include the modules index_square and lambda? I'd like to trace a run in your solution in order to understand the search strategy, and to see whether I can replicate it in ECLiPSe.Mclin
index_square.pl indeed contains some reusable code developed for the latter approach: i think could be reused, so moved it. I'll post it. Lambda it's from Ulrich Neumerkle, but because caused the problem that locked me initially (unexpected copying), i'll remove it. circuit/1 indeed should be the way to go, but it needs the entire graph. The main problem was that just a subset is involved.Hockey
@chac: I meant try circuit/1 on the subset. If necessary each subset. Yes it's inefficient, I only meant that if you try each subset anyway (it sounded like you did this originally) you might as well try circuit/1 on each subset instead of formulating essentially the same thing with other constraints.Tonie
+1. Nice solution! Chac: How was your problem of "unexpected copying" related to library(lambda)?Fraction
@mat: you are right. While developing, I reinvented the wheel. I was more concentrated on atomaton/8, but that eluded me...Hockey
@false: simply put: I losed the constraints on output variables, but I had to learn to inspect the CLP 'store' to understand what's was going on. Partly because I'm a rather enthusiast user of lambda (I tend to use it everywhere, cause I enjoy functional and compact code), but lambda it's somewhat unfriendly when debugging: tracing it's a pain, when acquainted to SWI-Prolog debugger, OTOH no breakpoints inside allowed, and stepping into it's very tedious and distracting (I find almost useless, best use writeln).Hockey
I never use a de-bugger, because before finding my errors I start to de-bug the de-bugger. Please look at this definition of $ and @. $ means trace call/exit. @ means there is at least a single answer. These are my only friends. Even when debugging foreign code, the first I do is to "install" this code. That is, paste it into the file.Fraction
I've tried Diadem.pl some time ago, as suggested by Ulrich when i prompted on SWI-Prolog mailing list, but i prefer an interactive debugger. Using Diadem I miss a 'report' phase, that allows filter the relevant info. You know, Prolog terms tends to be rather clumsy, and the main problem is to know where the bug is...The SWI-Prolog debugger it perfectible (for instance, a very simple change could make debugging DCG much more friendly, now it's a pain) but it's usable. CLP(FD) add complexity,and despite the support from debugger, I had some hard time to understand what was going on.Hockey
Did you really use @? I very much doubt it. You probably used $ only. In my code, I frequently leave @ in and it really helps.Fraction
let us continue this discussion in chatHockey
@CapelliC: could this solution be transfered to the original skyscraper puzzle?Hileman
@LasseAKarlsen: sorry, I don't understand. Could you elaborate your question ?Hockey
@CapelliC: I'm trying to implement a solution to skyscraper puzzles in Prolog, and thought perhaps your solution to this puzzle could be useful.Hileman
@LasseAKarlsen: yes, the code above is parametric vs grid dimension, should work. Try it with the obvious changes, and post a question if you have problems. Skyscrapes solver it's just the first part. Fences were more difficult to solve.Hockey
@CapelliC: I really appreciate it. From what I can tell you are using constraint propagation in order to solve the skyscraper part of the puzzle, and also some kind of depth-first search? With the original skyscraper puzzle, you fill the squares with digits 1 to N, N being the size of the board. How can I expand your solution so that this will be encounted for?Hileman
@LasseAKarlsen: no, I don't remember using depth search. Seems that the scyscrapes puzzle could be a bit simpler than this one. all_different with the additional constraint (see constraint_views/2 in my code) could be sufficient.Hockey
let us continue this discussion in chatHileman
F
7

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.

Fraction answered 10/12, 2011 at 21:3 Comment(2)
Thank you for your answer! Thing is, in solve_vertices/5, for example, I need to use reification because I must force the number of successes to be equal to C: sum([B1, B2, B3, B4], #= , C).Beale
You can post both. For other comments, see above.Fraction
M
5

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!

Mclin answered 13/12, 2011 at 11:46 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.