path finding - Optimizing pathfinding in Constraint Logic Programming with Prolog

Question

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 (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) 
        ) 
  

Answer 1

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 +