37
votes

I'm trying to solve a constraint processing problem in prolog.

I need to pack 4 squares of 5x5,4x4,3x3 and 2x2 in a grid of 10x10. They may not overlap.

My variables look like this:

Name: SqX(i), i=1..10, domain: 1..10

Where X is either 5,4,3 or 2. The index i represents the row, the domain the column in the grid.

My first constraints try to define the width and height of the squares. I formulate it as such:

Constraint: SqX(i) > SqX(j)-X /\ i>j-X, range: i>0 /\ j>0

So that the possible points are constrained to be within X rows and columns from each other. Prolog however, stops on these constraints and gives the following result:

Adding constraint "(Sq5_I > Sq5_J-5) /\ (I>J-5)" for values:
        I=1, J=1, 
        I=1, J=2, 
        I=1, J=3, 
        I=1, J=4, 
        I=1, J=5, 
        I=1, J=6, 
=======================[ End Solutions ]=======================

So it stops there, not even checking the other squares. My constraints are most likely too tight, but I can't see why or how. Any suggestions?

5

5 Answers

19
votes

For each square, define X and Y variables that denote the upper left corner. These variable will have domains 1..10-L, where L is the length of the square. If you set the domain to 1..10, the squares may be placed partly outside your 10x10 rectangle.

Then you can post constraints for each pair of rectangles (X,Y) and (X1,Y1) that state that if they overlap on the x axis, they must not overlap on the y axis, and vice versa:

(((X  #=< X1) and (X+L   #> X1)) => ((Y+L #=< Y1) or (Y1+L1 #=< Y))),
(((X1 #=< X)  and (X1+L1 #> X))  => ((Y+L #=< Y1) or (Y1+L1 #=< Y))),
(((Y  #=< Y1) and (Y+L   #> Y1)) => ((X+L #=< X1) or (X1+L1 #=< X))),
(((Y1 #=< Y)  and (Y1+L1 #> Y))  => ((X+L #=< X1) or (X1+L1 #=< X)))

(your particular constraint syntax may vary)

21
votes

Since version 3.8.3 SICStus Prolog offers a number of dedicated placement constraints that match your packing problem nicely. In particular, as your packing problem is two-dimensional, you should consider using the disjoint2/1 constraint.

The following code snippet uses disjoint2/1 to express that rectangles are non-overlapping. The main relation is area_boxes_positions_/4.

:- use_module(library(clpfd)).
:- use_module(library(lists)).

area_box_pos_combined(W_total*H_total,W*H,X+Y,f(X,W,Y,H)) :-
    X #>= 1,
    X #=< W_total-W+1,
    Y #>= 1,
    Y #=< H_total-H+1.

positions_vars([],[]).
positions_vars([X+Y|XYs],[X,Y|Zs]) :-
    positions_vars(XYs,Zs).

area_boxes_positions_(Area,Bs,Ps,Zs) :-
    maplist(area_box_pos_combined(Area),Bs,Ps,Cs),
    disjoint2(Cs),
    positions_vars(Ps,Zs).

On to some queries! First, your initial packing problem:

?- area_boxes_positions_(10*10,[5*5,4*4,3*3,2*2],Positions,Zs),
   labeling([],Zs).
Positions = [1+1,1+6,5+6,5+9],
Zs        = [1,1,1,6,5,6,5,9] ? ...

Next, let's minimize the total area that is required for placing all squares:

?- domain([W,H],1,10),
   area_boxes_positions_(W*H,[5*5,4*4,3*3,2*2],Positions,Zs),
   WH #= W*H,
   minimize(labeling([ff],[H,W|Zs]),WH).
W         = 9,
H         = 7,
Positions = [1+1,6+1,6+5,1+6],
Zs        = [1,1,6,1,6,5,1,6],
WH        = 63 ? ...

Visualizing solutions

What do individual solutions actually look like? ImageMagick can produce nice little bitmaps...

Here's some quick-and-dirty code for dumping the proper ImageMagick command:

:- use_module(library(between)).
:- use_module(library(codesio)).

drawWithIM_at_area_name_label(Sizes,Positions,W*H,Name,Label) :-
    Pix = 20,

    % let the ImageMagick command string begin
    format('convert -size ~dx~d xc:skyblue', [(W+2)*Pix, (H+2)*Pix]),

    % fill canvas 
    format(' -stroke none -draw "fill darkgrey rectangle ~d,~d ~d,~d"', 
           [Pix,Pix, (W+1)*Pix-1,(H+1)*Pix-1]),

    % draw grid
    drawGridWithIM_area_pix("stroke-dasharray 1 1",W*H,Pix),

    % draw boxes
    drawBoxesWithIM_at_pix(Sizes,Positions,Pix),

    % print label
    write( ' -stroke none -fill black'),
    write( ' -gravity southwest -pointsize 16 -annotate +4+0'),
    format(' "~s"',[Label]),

    % specify filename
    format(' ~s~n',[Name]).

Above code for drawWithIM_at_area_name_label/5 relies on two little helpers:

drawGridWithIM_area_pix(Stroke,W*H,P) :-   % vertical lines
    write(' -strokewidth 1 -fill none -stroke gray'),
    between(2,W,X),
    format(' -draw "~s path \'M ~d,~d L ~d,~d\'"', [Stroke,X*P,P, X*P,(H+1)*P-1]),
    false.
drawGridWithIM_area_pix(Stroke,W*H,P) :-   % horizontal lines
    between(2,H,Y),
    format(' -draw "~s path \'M ~d,~d L ~d,~d\'"', [Stroke,P,Y*P, (W+1)*P-1,Y*P]),
    false.
drawGridWithIM_area_pix(_,_,_).

drawBoxesWithIM_at_pix(Sizes,Positions,P) :-
    Colors = ["#ff0000","#00ff00","#0000ff","#ffff00","#ff00ff","#00ffff"],
    write(' -strokewidth 2 -stroke white'),
    nth1(N,Positions,Xb+Yb),
    nth1(N,Sizes,    Wb*Hb),
    nth1(N,Colors,   Color),
    format(' -draw "fill ~sb0 roundrectangle ~d,~d ~d,~d ~d,~d"',
           [Color, Xb*P+3,Yb*P+3, (Xb+Wb)*P-3,(Yb+Hb)*P-3, P/2,P/2]),
    false.
drawBoxesWithIM_at_pix(_,_,_).

Using the visualizers

Let's use the following two queries to produce some still images.

?- drawWithIM_at_area_name_label([5*5,4*4,3*3,2*2],[1+1,6+1,6+5,1+6],9*7,
                                 'dj2_9x7.gif','9x7').

?- drawWithIM_at_area_name_label([5*5,4*4,3*3,2*2],[1+1,1+6,5+6,5+9],10*10,
                                 'dj2_10x10.gif','10x10').

Let's use the following hack-query to produce an image for each solution of the placement of above rectangles on a board of size 9*7:

?- retractall(nSols(_)), 
   assert(nSols(1)), 
   W=9,H=7,
   Boxes = [5*5,4*4,3*3,2*2],
   area_boxes_positions_(W*H,Boxes,Positions,Zs),
   labeling([],Zs), 
   nSols(N), 
   retract(nSols(_)), 
   format_to_codes('dj2_~5d.gif',[N],Name),
   format_to_codes('~dx~d: solution #~d',[W,H,N],Label),
   drawWithIM_at_area_name_label(Boxes,Positions,W*H,Name,Label),
   N1 is N+1,
   assert(nSols(N1)),
   false.

Next, execute all ImageMagick commands output by above queries.

At last, build an animation of the solution set of the third query using ImageMagick:

$ convert -delay 15  dj2_0.*.gif   dj2_9x7_allSolutions_1way.gif 
$ convert dj2_9x7_allSolutions_1way.gif -coalesce -duplicate 1,-2-1 \
          -quiet -layers OptimizePlus -loop 0 dj2_9x7_allSolutions.gif

Results

First, one solution for board size 10*10: 10x10: one solution

Second, one solution for a board of minimum size (9*7): 9x7: one solution

Last, all solutions for a board of minimum size (9*7): 9x7: all solutions


Edit 2015-04-14

Since version 7.1.36 the SWI-Prolog clpfd library supports the constraint disjoint2/1.

Edit 2015-04-22

Here's a sketch of an alternative implementation based on the tuples_in/2 constraint:

  1. For each pair of boxes determine all positions at which these two would be non-overlapping.
  2. Encode the valid combinations as lists of tuples.
  3. For each pair of boxes post one tuples_in/2 constraint.

As a private proof-of-concept, I implemented some code following that idea; like @CapelliC in his answer, I get 169480 distinct solutions for the boxes and board-size the OP stated.

The runtime is comparable to the other clp(FD) based answers; in fact it is very competitive for small boards (10*10 and smaller), but gets a lot worse with larger board sizes.

Please acknowledge that, for the sake of decency, I refrain from posting the code:)

6
votes

There are already several great solutions posted here (+1 for all!), using CLP(FD) constraints.

In addition, I would like to show one conceptually different way to solve such placement and covering tasks, using CLP(B) constraints.

The idea is to consider each possible placement of a tile as a set of TRUE values at specific elements on the grid, where each grid element corresponds to one column of a matrix, and each possible placement of a tile corresponds to one row. The task is then to select a set of rows of said matrix in such a way that each grid element is covered at most once, or in other words, there is at most one TRUE value in each column of the submatrix consisting of the selected rows.

In this formulation, the selection of rows — and hence the placement of tiles at specific positions — is indicated by Boolean variables, one for each row of the matrix.

Here is the code I would like to share, it works in SICStus Prolog and SWI with at most small changes:

:- use_module(library(clpb)).
:- use_module(library(clpfd)).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   The tiles we have available for placement.

   For example, a 2x2 tile is represented in matrix form as:

       [[1,1],
        [1,1]]

   1 indicates which grid elements are covered when placing the tile.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

tile(5*5).
tile(4*4).
tile(3*3).
tile(2*2).

tile_matrix(Rows) :-
        tile(M*N),
        length(Rows, M),
        maplist(length_list(N), Rows),
        append(Rows, Ls),
        maplist(=(1), Ls).

length_list(L, Ls) :- length(Ls, L).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Describe placement of tiles as SAT constraints.

   Notice the use of Cards1 to make sure that each tile is used
   exactly once. Remove or change this constraint if a shape can be
   used multiple times, or can even be omitted.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

placement(M, N, Vs, *(Cs) * *(Cards1)) :-
        matrix(M, N, TilesRows),
        pairs_keys_values(TilesRows, Tiles, Rows),
        same_length(Rows, Vs),
        pairs_keys_values(TilesVs0, Tiles, Vs),
        keysort(TilesVs0, TilesVs),
        group_pairs_by_key(TilesVs, Groups),
        pairs_values(Groups, SameTiles),
        maplist(card1, SameTiles, Cards1),
        Rows = [First|_],
        phrase(all_cardinalities(First, Vs, Rows), Cs).

card1(Vs, card([1], Vs)).

all_cardinalities([], _, _) --> [].
all_cardinalities([_|Rest], Vs, Rows0) -->
        { maplist(list_first_rest, Rows0, Fs, Rows),
          pairs_keys_values(Pairs0, Fs, Vs),
          include(key_one, Pairs0, Pairs),
          pairs_values(Pairs, Cs) },
        [card([0,1], Cs)],
        all_cardinalities(Rest, Vs, Rows).

key_one(1-_).

list_first_rest([L|Ls], L, Ls).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   We build a matrix M_ij, where each row i describes what placing a
   tile at a specific position looks like: Each cell of the grid
   corresponds to a unique column of the matrix, and the matrix
   entries that are 1 indicate the grid positions that are covered by
   placing one of the tiles at the described position. Therefore,
   placing all tiles corresponds to selecting specific rows of the
   matrix such that, for the selected rows, at most one "1" occurs in
   each column.

   We represent each row of the matrix as Ts-Ls, where Ts is the tile
   that is used in each case.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

matrix(M, N, Ms) :-
        Squares #= M*N,
        length(Ls, Squares),
        findall(Ts-Ls, line(N, Ts, Ls), Ms).

line(N, Ts, Ls) :-
        tile_matrix(Ts),
        length(Ls, Max),
        phrase((zeros(0,P0),tile_(Ts,N,Max,P0,P1),zeros(P1,_)), Ls).

tile_([], _, _, P, P) --> [].
tile_([T|Ts], N, Max, P0, P) -->
        tile_part(T, N, P0, P1),
        { (P1 - 1) mod N >= P0 mod N,
          P2 #= min(P0 + N, Max) },
        zeros(P1, P2),
        tile_(Ts, N, Max, P2, P).

tile_part([], _, P, P) --> [].
tile_part([L|Ls], N, P0, P) --> [L],
        { P1 #= P0 + 1 },
        tile_part(Ls, N, P1, P).

zeros(P, P)  --> [].
zeros(P0, P) --> [0], { P1 #= P0 + 1 }, zeros(P1, P).

The following query illustrates which grid elements are covered (1), where each row corresponds to the placement of one of the rectangles:

?- M = 7, N = 9, placement(M, N, Vs, Sat), sat(Sat),
  labeling(Vs), matrix(M, N, Ms), pairs_values(Ms, Rows),
  pairs_keys_values(Pairs0, Vs, Rows),
  include(key_one, Pairs0, Pairs1), pairs_values(Pairs1, Covers),
  maplist(writeln, Covers).
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,0,0,0]
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,0,0,1,1,1,1,0,0,0,0,0,1,1,1,1]
[0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
[0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
M = 7,
N = 9,
etc.

corresponding to the solution:

solution for original problem

Such a CLP(B) formulation is typically less scalable than a CLP(FD) version, also because there are more variables involved. However, it also has a few advantages:

One significant advantage is that it is readily generalized to a version of the task where some or all of the shapes can be used multiple times. For example, in the version above, we can simply change card1/2 to:

custom_cardinality(Vs, card([0,1,2,3,4,5,6,7], Vs)).

and obtain a version where each tile can be used up to 7 times, and can even be omitted entirely (due to the inclusion of 0).

Second, we can easily turn this into a solution for an exact cover problem, which means that each grid element is covered by one of the shapes, by simple changing card([0,1], Cs) to card([1], Cs) in all_cardinalities//3.

Together with the other modification, here is a covering for a 4x4 grid using four 2x2 rectangles:

[1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0]
[0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0]
[0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0]
[0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1]

A third advantage of the CLP(B) formulation is that the number of solutions can be computed without enumerating the solutions explicitly. For example, for the original task:

?- placement(7, 9, Vs, Sat), sat_count(Sat, Count).
Count = 68.

These 68 solutions are already beautifully illustrated by @repeat.

For comparison, here is the number of solutions where each shape can be used between 0 and 7 times:

?- placement(7, 9, Vs, Sat), time(sat_count(Sat, Count)).
% 157,970,727 inferences, 19.165 CPU in 19.571 seconds
...
Count = 17548478.

The same on a 10x10 grid, computed in about 6 minutes (~ 2 billion inferences):

?- placement(10, 10, Vs, Sat), sat_count(Sat, Count).
Count = 140547294509.

And on an 11x11 grid, computed in about half an hour (~ 9 billion inferences):

?- placement(11, 11, Vs, Sat), sat_count(Sat, Count).
Count = 15339263199580.

Lastly, and maybe most significantly, this approach works for any shape of tiles, and is not limited to squares or rectangles. For example, to handle 1x1 squares and a triangle shape as well as its vertical and horizontal reflections, use the following definition of tile_matrix/1:

tile_matrix([[1]]).
tile_matrix(T) :-
        T0 = [[1,1,1,1],
              [1,1,1,0],
              [1,1,0,0],
              [1,0,0,0]],
        (   T = T0
        ;   maplist(reverse, T0, T)
        ;   reverse(T0, T)
        ).

Allowing each of these shapes to be used between 0 and 7 times on a 9x7 board, I get, after a minute or so, Count = 58665048314 solutions.

Here is one of them, picked at random:

example with triangular shapes

Picking solutions in such a way that each of them is equally likely is also quite easy with CLP(B), even if the number of solutions is too large to enumerate them explicitly.

4
votes

I coded in SWI-Prolog

/*  File:    pack_squares.lp
    Author:  Carlo,,,
    Created: Nov 29 2012
    Purpose: http://stackoverflow.com/questions/13623775/prolog-constraint-processing-packing-squares
*/

:- module(pack_squares, [pack_squares/0]).
:- [library(clpfd)].

pack_squares :-
    maplist(square, [5,4,3,2], Squares),
    flatten(Squares, Coords),
    not_overlap(Squares),
    Coords ins 1..10,
    label(Coords),
    maplist(writeln, Squares),
    draw_squares(Squares).

draw_squares(Squares) :-
    forall(between(1, 10, Y),
           (   forall(between(1, 10, X),
              sumpts(X, Y, Squares, 0)),
           nl
           )).

sumpts(_, _, [], S) :- write(S).
sumpts(X, Y, [[X1,Y1, X2,Y2]|Qs], A) :-
    ( ( X >= X1, X =< X2, Y >= Y1, Y =< Y2 )
    ->  B is A+X2-X1+1
    ;   B is A
    ),
    sumpts(X, Y, Qs, B).

square(D, [X1,Y1, X2,Y2]) :-
    X1 + D - 1 #= X2,
    Y1 + D - 1 #= Y2.

not_overlap([_]).
not_overlap([A,B|L]) :-
    not_overlap(A, [B|L]),
    !, not_overlap([B|L]).

not_overlap(_, []).
not_overlap(Q, [R|Rs]) :-
    not_overlap_c(Q, R),
    not_overlap_c(R, Q),
    not_overlap(Q, Rs).

not_overlap_c([X1,Y1, X2,Y2], Q) :-
    not_inside(X1,Y1, Q),
    not_inside(X1,Y2, Q),
    not_inside(X2,Y1, Q),
    not_inside(X2,Y2, Q).

not_inside(X,Y, [X1,Y1, X2,Y2]) :-
    X #< X1 #\/ X #> X2 #\/ Y #< Y1 #\/ Y #> Y2.

here is the last lines displayed when running ?- aggregate_all(count,pack_squares,C)., notably C counts total placements

...
0002255555
0002255555
[6,6,10,10]
[7,2,10,5]
[4,3,6,5]
[5,1,6,2]
0000220000
0000224444
0003334444
0003334444
0003334444
0000055555
0000055555
0000055555
0000055555
0000055555
C = 169480.
1
votes

Here is a solution where disjoint only takes one line:

% disjoint(+Rectangle, +Rectangle)
disjoint([XA1,XA2,YA1,YA2],[XB1,XB2,YB1,YB2]) :-
   XB1 #>= XA2 #\/ XA1 #>= XB2 #\/
   YB1 #>= YA2 #\/ YA1 #>= YB2.

The model setup and labeling works as follows:

% squares(-List)
squares(L) :-
   maplist(square, [2,3,4,5], L),
   term_variables(L, V),
   place(L),
   label(V).

% square(+Integer, -Rectangle)
square(S, [X1,X2,Y1,Y2]) :-
   X1 in 0..8,
   X2 in 1..9,
   Y1 in 0..6,
   Y2 in 1..7,
   X2 #= X1+S,
   Y2 #= Y1+S.

% place(+List)
place([]).
place([A|L]) :-
   place(L, A),
   place(L).

% place(+List, +Rectangle)
place([], _).
place([A|L], B) :-
   disjoint(A, B),
   place(L, B).

Here is an example run:

Jekejeke Prolog 3, Runtime Library 1.3.7 (May 23, 2019)

?- squares(L), show(L).
555554444
555554444
555554444
555554444
55555333.
22...333.
22...333.
L = [[0,2,5,7],[5,8,4,7],[5,9,0,4],[0,5,0,5]]