A Generic Geometrical Constraint Kernel in Space and Time for Handling Polymorphic k-Dimensional Objects

(1)

A Generic Geometrical Constraint Kernel in Space and

Time for Handling Polymorphic k-Dimensional Objects

N. Beldiceanu

1

_{, M. Carlsson}

2

_{, E. Poder}

1

_{, R. Sadek}

1

_{, and C. Truchet}

3 1

_´

Ecole des Mines de Nantes, LINA FRE CNRS 2729, FR-44307 Nantes, France

{

Nicolas.Beldiceanu,Emmanuel.Poder,Rida.Sadek

}

@emn.fr

2

SICS, P.O. Box 1263, SE-164 29 Kista, Sweden

Mats.Carlsson@sics.se

3

Universit´e de Nantes, LINA FRE CNRS 2729, FR-44322 Nantes, France

Charlotte.Truchet@univ-nantes.fr

SICS Technical Report T2007:08

ISSN: 1100-3154

ISRN: SICS-T–2007/08-SE

Abstract: This report introduces a geometrical constraint kernel for handling the location in

space and time of polymorphic k-dimensional objects subject to various geometrical and time

constraints. The constraint kernel is generic in the sense that one of its parameters is a set of

constraints on subsets of the objects. These constraints are handled globally by the kernel. We first

illustrate how to model several placement problems with the constraint kernel. We then explain

how new constraints can be introduced and plugged into the kernel. Based on these interfaces,

we develop a generic k-dimensional lexicographic sweep algorithm for filtering the attributes of

an object (i.e., its shape and the coordinates of its origin as well as its start, duration and end in

time) according to all constraints where the object occurs. Experiments involving up to hundreds

of thousands of objects and

1 million integer variables are provided in 2, 3 and 4 dimensions,

both for simple shapes (i.e., rectangles, parallelepipeds) and for more complex shapes.

(2)

1 Introduction and Presentation of the Kernel

This report introduces a constraint kernel geost

(k, O, S, C) for handling in a generic

way a variety of geometrical constraints

C in space and time between polymorphic

k-dimensional objects

O (k ∈ N

+

_{), where each object takes a shape among a set of}

shapes described by

S during a given time interval and at a given position in space.

This line of research can be seen as a continuation and generalisation of previous work

on non-overlapping parallelepipeds [1–4].

Each shape is defined as a finite set of shifted boxes, where each shifted box is

described by a box in a k-dimensional space at a given offset (from the origin of the

shape) with given sizes. More precisely, a shifted box s

= sbox(sid , t[], l[]) ∈ S is

an entity defined by its shape id s.sid , shift offset s.t

[d], 0 ≤ d < k, and sizes s.l[d]

(s.l[d] > 0, 0 ≤ d < k). All attributes of a shifted box are integer values. Then, a shape

is defined as the union of shifted boxes sharing the same shape id. Each object o

=

object(id , sid , x[], start , duration, end ) ∈ O is an entity defined by its unique object

id o.id , shape id o.sid , origin o.x

[d], 0 ≤ d < k, start in time o.start , duration in time

o.duration

(o.duration ≥ 0) and end in time o.end .

1

_{All these attributes correspond}

to domain variables.

2

_{Typical constraints from the list of constraints}

_{C can express, for}

instance, the fact that a given subset of objects from

O do not pairwise overlap or that

they are all included within a given bounding box. Constraints always have two first

arguments

A

i

and

O

i

(followed by possibly some additional arguments) which resp.

specify:

– A list of distinct dimensions (integers in

[0, k − 1]) that the constraint considers.

– A list of identifiers of the objects to which the constraint applies.

Example 1. Assume we have a 3D placement problem (i.e., k

= 3) involving a set of

paral-lelepipeds

P and one subset P

′

of

P, where we want to express the fact that (1) no parallelepipeds

of

P should overlap, and (2) no parallelepipeds of P

′

should be piled. Constraints (1) and (2)

resp. correspond to non-overlapping

([0, 1, 2], P) and to non-overlapping([0, 1], P

′

_{). Within the}

first non-overlapping constraint, the argument

[0, 1, 2] expresses the fact that we consider a

non-overlapping constraint according to dimensions

0, 1 and 2 (i.e., given any pair of

paral-lelepipeds p

′

and p

′′

of

P there should exist at least one dimension d (d ∈ {0, 1, 2}) where

the projections of p

′

and p

′′

on d do not overlap). Similarly, the argument

[0, 1] of the second

non-overlapping constraint expresses the fact that, given any pair of parallelepipeds p

′

and p

′′

of

P

′

, there should exist at least one dimension d

(d ∈ {0, 1}) where p

′

and p

′′

do not overlap).

geost(k, O, S, C) is defined in the following way: given a constraint ctr

i

(A

i

,

O

i

)

from the list of constraints

C between a subset of objects O

i

⊆ O according to the

attributes

A

i

, let

MC

i

denote the sets of cliques stemming from the objects of

O

i

that

all overlap in time.

3

_{The constraints of geost}

_{(k, O, S, C) hold if and only if ∀ctr}

i

∈ C,

∀O

MCi

∈ MC

i

: ctr

i

(A

i

,

O

MCi

) holds.

1

The time dimension is treated specially since the duration attribute may not be fixed, which

is not the case for the sizes of a shifted box. Also, the geometrical constraints only apply on

objects that intersect in time.

2

A domain variable v is a variable ranging over a finite set of integers denoted by

dom(v); let

v and v resp. denote the minimum and maximum possible values for v.

3

In fact, these cliques (of an interval graph) are only used for defining the declarative semantics

of geost ’s constraints.

(3)

1 1 2 2 3 5 4 3 4 1 1 4 2 3 4 5 3 2 3 5 4 3 2 1 1 2 4 time objects 22 14 2 1 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 S9 S8 S7 S6 S5 S4 S3 S2 S1

object O3 is assigned shape S8 object O2 is assigned shape S5 object O1 is assigned shape S1

(E)

first clique={O1} second clique={O1,O2,O3} third clique={O2,O3,O4}

(D) object O4 (C) object O3

(B) object O2 (A) object O1

Shape of the fourth object

the third object

Potential shapes for

the second object

the first object

interval [14,21] O4 O3 O2 O1 O4 O2 O3 O4 O3 O1 O2 interval [2,9] O1 interval [10,13]

Fig. 1. Example with

4 objects, 9 shapes, one non-overlapping and one included constraints

Example 2. Fig. 1 presents a typical example of a dynamic 2D placement problem where one

has to place four objects, in time and within a given box, so that objects that overlap in time do

not overlap. Parts (A), (B), (C) and (D) resp. represent the potential shapes associated with the

four objects to place, where the origin of each object is represented by a black square

. Part (E)

shows the position of the four objects of the example as the time varies, where the first, second,

third and fourth objects were resp. assigned shapes S

1

, S

5

, S

8

and S

9

:

– During the first time interval

[2, 9] we have only object O

1

at position

(1, 2).

– Then, at instant

10 objects O

2

and O

3

both appear. Their origins are resp. placed at positions

(2, 1) and (4, 1).

– At instant

14 object O

1

disappears and is replaced by object O

4

. The origin of O

4

is fixed at

position

(1, 1). Finally, at instant 22 all three objects O

2

, O

3

and O

4

disappear.

The corresponding arguments are:

01 geost(2,

02 [object(1,1,[1,2], 2,12,14), object(2,5,[2,1],10,12,22), 03 object(3,8,[4,1],10,12,22), object(4,9,[1,1],14, 8,22)], 04 [sbox(1,[0,0],[2,1]), sbox(1,[0,1],[1,2]), sbox(1,[1,2],[3,1]), 05 sbox(2,[0,0],[3,1]), sbox(2,[0,1],[1,3]), sbox(2,[2,1],[1,1]), 06 sbox(3,[0,0],[2,1]), sbox(3,[1,1],[1,2]), sbox(3,[2,2],[3,1]), 07 sbox(4,[0,0],[3,1]), sbox(4,[0,1],[1,1]), sbox(4,[2,1],[1,3]), 08 sbox(5,[0,0],[2,1]), sbox(5,[1,1],[1,1]), sbox(5,[0,2],[2,1]), 09 sbox(6,[0,0],[3,1]), sbox(6,[0,1],[1,1]), sbox(6,[2,1],[1,1]), 10 sbox(7,[0,0],[3,2]), sbox(8,[0,0],[2,3]), sbox(9,[0,0],[1,4])],

11 [non-overlapping([0,1],[1,2,3,4]),included([0,1],[1,2,3,4],[1,1],[5,4])])

Its first argument

2 is the number of dimensions of the placement space we consider. Its second

(4)

we have. For instance, the

3 boxes of shape S

1

(depicted by

3 thick rectangles in Part (A) of

Fig. 1) respectively correspond to the

3 boxes declared at line 04 of the example. Finally, its last

argument gives the list of geometrical constraints imposed by geost : the first constraint expresses

a non-overlapping constraint between the four objects, while the second constraint imposes the

four objects to be located within the box containing all points

(x, y) such that 1 ≤ x ≤ 1 + 5 − 1

and

1 ≤ y ≤ 1+4−1. The constraints of geost hold since the four objects do not simultaneously

overlap in time and in space and since they are completely included within the previous box (i.e.,

see Part (E) of Fig. 1).

Within the scope of geost

(k, O, S, C), this report presents a filtering algorithm that

prunes the domain of each attribute of every object o

= object(id , sid , x[], start ,

duration, end

) ∈ O. All values found infeasible are deleted from the shape attribute

sid ; for the other attributes (i.e., the origin x[], the start, the duration and the end), the

minimum and maximum are adjusted. The approach presented in this report offers a

number of advantages:

– The main theoretical advantages are fourfold:

• First, the geometrical kernel makes it possible to integrate new geometrical

constraints as new applications and/or requirements show up. This is achieved

by providing for each geometrical constraint an API without knowing any

de-tails about the geometrical kernel. This contrasts with traditional approaches

where one has to come up with a rather involved filtering algorithm for each

global constraint.

• Second, while pruning the attribute of an object, the geometrical kernel takes

direct advantage of all geometrical constraints involving that object in order to

perform more deduction. This is a fundamental progress over the traditional

ap-proach where constraints only co-operate through the domains of their shared

variables.

• Even when we have three or four dimensions, the approach scales well since it

does not rely on building complex multi-dimensional data structure (e.g., like

quadtrees or octrees). It only stores a number of points in the order of O

(m · k)

where m is the total number of objects and k is the number of dimensions.

• Even if complex objects could be decomposed into boxes for which one links

the coordinates by external equality constraints this weakens a lot the

deduc-tion process as illustrated by the following example of Figure 2: if the shape s

(see Part (A)) is decomposed into two rectangles r

4 and r5 (see Part (C)) and

if the constraints linking the coordinates of the origins of r

4 and r5 are not

in-tegrated within the sweep process, infeasibility cannot be directly derived (see

Part (M)). In contrast our approach allows to detect infeasibility directly by

reasoning only on the coordinates of the origin of s.

– The main practical advantages are as follows:

• Having k dimensions allows to come up with a single constraint that can be

used for handling general non-overlapping constraints. This was originally

mo-tivated by a warehouse management problem where both two-dimensional and

three-dimensional sub-problems had to be solved. In the context of

three-di-mensional packing problems having an extra dimension also makes sense for

modelling the fact that we want to assign objects to a truck (in this context we

(5)

(M) values for (x1,y1)

Forbidden pairs of

ctr1: x1+2=x2 ctr2: y1+1=y2

ctr3: r1 and r4 do not overlap

ctr8: r3 and r5 do not overlap ctr7: r2 and r5 do not overlap ctr6: r1 and r5 do not overlap ctr4: r2 and r4 do not overlap ctr5: r3 and r4 do not overlap

(D) (C) (A) (B) = ctr2: ctr1: = +2 +1 2 1 3 1 = ctr2: ctr1: = +2 +1 4 3 5 3

and ctr2 are not integrated within the sweep algorithm external constraints ctr1 and ctr2 is not done if ctr1 to forbidden pairs of values for (x1,y1) through the Transmission of the forbidden pairs of values for (x2,y2)

for (x1,y1) pair of values Overall set of forbidden ctr6, ctr7 and ctr8 according to ctr8

according to ctr7 according to ctr6

Forbidden pairs of values for (x1,y1) according to ctr3, ctr4 and ctr5 Forbidden pairs of

values for (x1,y1) according to ctr5

Constraints of the problem where to put s

Placement space Decomposing s in two rectangles r4 and r5 Shape to place within (B)

so that it does not overlap rectangles r1, r2, r3

(N) Forbidden pairs of

values for (x2,y2)

Forbidden pairs of

values for (x2,y2) values for (x2,y2)

Forbidden pairs of Forbidden pairs of values for (x2,y2) according to (F)

(E) (G) (H)

(J)

(I) (K) (L)

according to ctr3 according to ctr4 values for (x1,y1) Forbidden pairs of 5 4 3 2 1 y1 x1 2 1 y2 x2

00000

11111

y1 x1 5

00000

11111

5 4 3 2 1 5 4 3 2 1 y2 x2 4

00000

11111

5 4 3 2 1 5 4 3 2 1 y1 x1 3 5 4 3 2 1 3 2 1

00000

11111

5 4 3 2 1 5 4

00000

11111

x2 x1 r5 r4 y1 y2 3 2 1 5

00000

11111

5 4 3 2 1 5 4 3 2 1 y1 x1 4 3 2 1 y2 x2

00000

11111

00000

11111

r3 r2 r1 s y1 x1 5 4 3 2 1 5

00000

11111

5 4 3 2 1 5 4

00000

11111

5 4 3 2 1

00000

11111

5 4 3 5 4 3 2 1 5 4 3 2 1 y2 x2 2 1 5 4 3 2 1 y1 x1 5 4 3 2 1 y2 x2 3 2 1 y1 x1 4

Fig. 2. Reasoning for detecting the infeasibility of the placement problem (dashed areas

corre-spond to initially forbidden pairs of values, while grey areas represent forbidden pairs of values

related to some non-overlapping constraints)

(6)

speak about an assignment dimension – see Part (H) of Figure 3) or the fact that

we do not want to place all the objects since there may simply be not enough

room (in this context we speak about a relaxation dimension).

• Factoring out the description of the shapes from the description of an object

makes sense in a lot of practical problems where a number of instances of the

same shape have to be considered (this is illustrated by Part (I) of Figure 3

where we have five objects but only three shapes: in fact the first, third and

fifth objects correspond to the first shape). This again occurs in the warehouse

management problem that originally motivated the constraint, where a major

car manufacturer has to pack within the same container the parts associated

with

24 instances of the same car model. By doing so we can decrease the

memory requirement (i.e., each complex shape is represented only once).

• Having a set of potential shapes for an object offers an extra modelling power

for representing directly the fact that objects may rotate, or for dealing with

tasks for which the duration depends on the machine where the task is actually

assigned.

• Having a temporal dimension allows to tackle dynamic placement problems

where objects are moving in time. Consider for instance a pick-up delivery

problem where objects are loaded or unloaded from a truck while visiting

dif-ferent locations. In this context the non-overlapping constraint applies only for

those objects which overlap in time. This is illustrated by Part (I) of Figure 3.

The report is organised as follows. Section 2 provides an overview of placement

problems that can be modelled with the constraints currently available in geost .

Sec-tion 3 presents the overall architecture of the geometrical kernel. It explains how to

define geometrical constraints in terms of a programming interface by the geometrical

kernel. Section 4 focusses on the main contribution of this report: a multi-dimensional

lexicographic sweep algorithm used for filtering the attributes of an object of geost .

Section 5 evaluates the scalability of the geost kernel as well as its ability to deal with a

variety of specific placement problems. Before we conclude, Section 6 compares geost

with related work and suggests future directions. Finally an annex provides the data sets

and test programs used in Section 5.

2 Modelling Problems with geost

As illustrated by Fig. 3 in the context of non-overlapping, geost allows to model directly

a large number of placement problems:

– Case (A) corresponds to a non-overlapping constraint among three segments.

– The second and third cases (B,C) correspond to a non-overlapping constraint

be-tween rectangles where (B) is a special case where the sizes of all rectangles in the

second dimension are equal to

1; this can be interpreted as a machine assignment

problem.

– Case (D) corresponds to a non-overlapping constraint between rectangles where

each rectangle can have two orientations. This is achieved by associating with each

rectangle two shapes of respective sizes l

×h and h×l. Since their orientation is not

(7)

initially fixed, the included constraint enforces the three rectangles to be included

within the bounding box defined by the origin’s coordinates

1, 1 and sizes 8, 3.

– Case (E) corresponds to a non-overlapping constraint between more complex

ob-jects where each object is described by a given set of rectangles.

– Case (F) describes a placement problem where one has to first assign each rectangle

to a strip so that all rectangles that are assigned to the same strip do not overlap.

– Case (G) corresponds to a non-overlapping constraint between parallelepipeds.

– Case (H) can be interpreted as a non-overlapping constraint between parallelepipeds

that are assigned to the same container. The first dimension corresponds to the

identifier of the container, while the next three dimensions are associated with the

position of a parallelepiped inside a container.

– Case (I) describes a rectangle placement problem over three consecutive time-slots:

rectangles assigned to the same time-slot should not overlap in time. We initially

start with the three rectangles

1, 2 and 3. Rectangle 3 is no longer present at instant 2

(the triangle H within rectangle

3 at time 1 indicates that rectangle 3 will disappear

at the next time-point), while rectangle

4 appears at instant 2 (the triangle N within

rectangle

4 at time 2 denotes the fact that the rectangle 4 appears at instant 2).

Finally, rectangle

2 disappears at instant 3 and is replaced by rectangle 5.

2 1

(H)

(I)

(G)

(F)

(E)

(D)

(C)

(B)

(A)

2 3 1 2 3 1 2 1 2 3 4 1 2 3 4 5 1 2 34 5 5 2 3 4 3 2 1 3 2 1 6 7 8 5 4 3 2 1 3 2 1 6 7 8 5 4 3 2 1 5 6 7 8 1 2 3 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 1 1 8 7 5 4 3 2 1 1 2 3 4 1 2 3 1 2 3 1 3 2 6 23 1 1 2 3 1 2 3 1 2 4 4 2 3 5 2 1 1 1 3 2 1 3 2 1 3 1 3 2 1 3 2 1 2 2 3 3 1 3 2 2 3 1 time=3 time=2 time=1

Fig. 3. Nine typical examples of use of geost

3 Standard Representation of Geometrical Constraints

The key idea for handling multiple geometrical constraints in a common kernel is

the following. For each type of geometrical constraint found in

C (also called

exter-nal constraints), one has to provide a service that computes necessary conditions (also

called internal constraints) for a given object and shape. Given an external geometrical

constraint ectr

i

(A

i

,

O

i

) (A

i

⊆ {0, 1, . . . , k − 1}, O

i

⊆ O), one of its object o ∈ O

i

(8)

and one potential shape s of o, such a necessary condition generated by ectr

i

, o and s is

a unary

4

_{constraint ictr}

_{(o.x) such that: o.sid = s ∧ ectr}

i

(A

i

,

O

i

) ⇒ ictr (o.x). Now,

the key to being able to globally treat such necessary conditions in the kernel is to give

them a uniform representation. We have chosen the following one:

– A constraint

outbox(t, l) on o.x holds iff o.x is located outside the shifted box

defined by its origins point t

[d], 0 ≤ d < k, and sizes l[d], 0 ≤ d < k (i.e.,

∃d ∈ [0, k − 1] | o.x[d] < t[d] ∨ o.x[d] > t[d] + l[d] − 1).

Thus, an outbox corresponds to a box-shaped set of points that are infeasible for

o.x. The purpose of the introduction of outboxes is to have a common representation

for the kernel, suitable for the k-dimensional lexicographic sweep algorithm presented

in the next section, which considers all the outboxes, for a selected object and shape, in

one run.

Consequently, for each type of external geometrical constraint, found in

C a service

GenOutboxes(ectr

i

, o, s) : (ictrs), responsible for generating outboxes, must be

pro-vided. This service is assumed to generate outboxes that intersect the domains of the

origin coordinates of o. Also, if all attributes mentioned by ectr

i

belonging to objects

other than o are fixed, those outboxes are assumed to be necessary and sufficient

condi-tions, lest the kernel accept false solutions.

Example of External Geometrical Constraints. We now illustrate some external

geo-metrical constraints that are currently available within the constraint kernel. As we saw

in the introduction, an external constraint always has at least two arguments that resp.

correspond to a list of distinct dimensions and to a list of object identifiers to which the

constraint apply.

The included and non-overlapping external constraints. The included

(A

i

,

O

i

, t, l)

and the non-overlapping

(A

i

,

O

i

) external constraints take as input a list of distinct

dimensions

A

i

in

{0, 1, . . . , k − 1} and a list O

i

of distinct object identifiers of geost .

In addition, the included constraint considers a shifted box defined by its origin point

t[d], 0 ≤ d < k, and size l[d], 0 ≤ d < k.

The included constraint enforces for each object o (with o.id

∈ O

i

) and for any

corresponding shifted box s (with o.sid

= s.sid ) the condition ∀d ∈ A

i

| t[d] ≤

o.x[d] + s.t[d] ∧ o.x[d] + s.t[d] + s.l[d] − 1 ≤ t[d] + l[d] − 1 (i.e., s is included within

the shifted box attribute defined by the parameters t and l of the included constraint).

Depending on which shape of an object we actually consider, the included constraint

can be translated to

2k outbox constraints.

The non-overlapping constraint enforces the following condition: given two distinct

objects o and o

′

(with o.id , o

′

.id

∈ O

i

) that overlap in time, no shifted box s (with

o.sid

= s.sid ) should overlap any shifted box s

′

(with o

′

.sid

= s

′

.sid ); i.e. it should

hold that

∃d ∈ A

i

| o.x[d] + s.t[d] + s.l[d] ≤ o

′

.x[d] + s

′

.t[d] ∨ o

′

.x[d] + s

′

.t[d] +

s

′

.l[d] ≤ o.x[d] + s.t[d] (i.e., there exists a dimension where they do not intersect).

While focussing on an object o we can easily generate an

outbox constraint for each

object o

′

_{that should not overlap o by reusing the results of [2].}

4

_{Unary, since it involves the k coordinates of a single object.}

(9)

4 The Geometrical Kernel: a Generic k-dimensional

Lexicographic Sweep Algorithm

In this section, we first present the sweep algorithm used for filtering the coordinates of

the origin of an object o of geost when each object has one single shape. We initially

assume that time is treated exactly like the space dimensions, i.e. that the o.x array is

extended by one element. Toward the end of this section, we explain in detail how to

treat the time attributes of an object. We also assume for now that the shape attribute is

fixed, and explain later how to handle multiple potential shapes for an object (i.e.,

poly-morphism). We now introduce some notation used throughout this section.

Notation. Assume v and w are vectors of scalars of k components. Then v

← w

denotes the element-wise assignment of w to v, w

+ d (resp. w − d) denotes the

ele-ment-wise addition of d (resp.

−d) to w. Given a scalar d, 0 ≤ d ≤ k − 1, rot(v, d, k)

denotes the vector

(v[d], v[(d + 1) mod k], . . . , v[(d − 1) mod k]). That is, in the

ro-tated vector, v

[d] is the most significant element, which is what we need when running

the sweep algorithm on dimension d.

The Sweep Algorithm. This algorithm first considers all outboxes

IC

o

derived from

C where object o actually appears, and then performs a recursive traversal of the

place-ment space for each coordinate and direction (i.e.,

min or max). Without loss of

gen-erality, assume we want to adjust the minimum value of the d

th

coordinate o.x

[d],

0 ≤ d < k, of the origin of o. The algorithm starts its recursive traversal of the

place-ment space at point c

= rot(o.x, d, k) and could in principle explore all points of the

domains of o.x, one by one, in increasing lexicographic order, until a point is found

that is not inside any outbox, in which case c

[0] is the computed new minimum value.

To make the search efficient, instead of moving each time to the successor point, we

arrange the search so that it skips points that are known to be inside some outbox.

5

Thus, we compute the lexicographically smallest point c

′

such that:

1. c

′

is lexicographically greater than or equal to c,

2. every element of c

′

is in the domain of the corresponding element of o.x,

3. c

′

is not inside any outbox of

IC

o

.

If no such c

′

_{exists, the constraint fails. Otherwise, the minimum value of o.x}

_{[d] is}

adjusted to c

′

_{[0]. As we saw, the sweep algorithm moves in increasing lexicographic}

or-der a point c from its lexicographically smallest potential feasible position to its

lexico-graphically largest potential feasible position through all potential points. The algorithm

uses the following data structures:

– The current position c of the sweep.

– A vector n

[0..k − 1] that records knowledge about already encountered sets of

infeasible points while moving c from its first potential feasible position. The vector

n is always element-wise greater than c and maintained as follows. Let

inf, sup

denote the vectors

inf = rot(o.x, d, k) and sup = rot(o.x + 1, d, k):

(10)

• Initially, n = sup.

• Whenever an outbox f containing c is found, n is updated by taking the

ele-ment-wise minimal value of n and the upper boundary of

rot(f, d, k),

indicat-ing the fact that new candidate points can be found beyond that value.

• Whenever we skip to the next candidate point, we reset the elements of n that

were used to the corresponding values of

sup.

The following invariant holds for the vector n, and is used when advancing c to the

next candidate point. Let i be the smallest j such that n

[j + 1] = sup[j + 1] ∧ · · · ∧

n[k − 1] = sup[k − 1] and suppose c is known to be in some outbox. Then, the

next point, lexicographically greater than c and not yet known to be in any outbox,

is

(c[0], . . . , c[i − 1], n[i], inf[i + 1], . . . , inf[k − 1]).

Algorithm 1 implement this idea. The algorithm prunes the bounds of each

coordi-nate of every object wrt. its relevant outboxes, iterating to fix-point.

Efficiency. The main inefficiency in this sweep algorithm lies in searching the set

of outboxes (line 4 of

PruneMin). In order to make this search more efficient, we

can make the sweep algorithm more sophisticated by the following modifications to

PruneMin:

– We extend the state of the algorithm by an event point series, ordered in

lexi-cographically increasing order. These events correspond to the lexilexi-cographically

smallest (insert events) and largest (delete events) relevant infeasible point

associ-ated with each outbox ictr

o

∈ IC

o

. They are sorted in lexicographically increasing

order, and we maintain a pointer into the series in sync with point c.

– We maintain the set of active outboxes, corresponding to all outboxes ictr

o

∈ IC

o

such that c is between its lexicographically smallest and largest infeasible points.

This set is initially empty.

– When c is initialized in line 2 as well as when c is incremented in lines 7-17, the

relevant events up to point c from the event point series are processed, and the

corresponding outboxes are added to or deleted from the set of active outboxes.

– In line 4, only the active outboxes are considered.

Example 3. Fig. 4 illustrates the k-dimensional lexicographic sweep algorithm in the context of

k

= 2. Parts (A) and (B) provide the variables of the problem (i.e., the abscissa and ordinate

of each rectangle r

1

, r

2

, r

3

, r

4

and r

5

) as well as the non-overlapping constraint between the

five previous rectangles. On Part (D) we have represented the extreme possible feasible positions

of each rectangle (i.e., rectangles r

1

to r

4

): for instance the leftmost lower corner of rectangle

r

1

can only be fixed at positions

(1, 2), (1, 3), (1, 4), (2, 2), (2, 3), (2, 4), (3, 2), (3, 3), (3, 4),

(4, 2), (4, 3) and (4, 4). Parts (C) to (L) of Fig. 4 detail the different steps of the algorithm for

adjusting the minimum value of the abscissa of rectangle r

5

. Part (C) provides the outboxes

associated with the fact that we want to prune the coordinates of r

5

: constraints ctr

1

, ctr

2

, ctr

3

and ctr

4

resp. correspond to the fact that rectangle r

5

should not overlap rectangles r

1

, r

2

, r

3

and r

4

, while constraint ctr

5

represents the fact that the ordinate of r

5

should be different from

7. Part (D) represents the initialisation phase of the algorithm where we have all five outboxes

with their respective lexicographically smallest infeasible point (i.e.,

(1, 1) for ctr

1

,

(1, 3) for

ctr

2

,

(1, 7) for ctr

5

,

(1, 8) for ctr

3

and

(3, 1) for ctr

4

). Part (E) represents the first step of the

sweep algorithm where we start the traversal of the placement space at point c

= (1, 1). We first

(11)

00

11

00

11

00

11

00

11

FOR FILTERING THE ORIGIN OF THE INTERNAL CONSTRAINTS GENERATED

ctr1: outbox([1,1],[2,2]) ctr2: outbox([1,3],[6,4]) ctr3: outbox([1,8],[2,1]) ctr4: outbox([3,1],[5,3]) ctr5: outbox([1,7],[8,1]) DCTRS [object(1,1,[x1,y1],0,1,1),object(2,2,[x2,y2],0,1,1), geost( object(3,3,[x3,y3],0,1,1),object(4,4,[x4,y4]0,1,1), object(5,5,[x5,y5],0,1,1)], [shape(1,[0,2],[0,1]),shape(2,[0,3],[0,1]),shape(3,[0,1],[0,1]), shape(4,[0,1],[0,3]),shape(5,[0,5],[0,4])], [non−overlapping([0,1],[1,2,3,4,5])] )

FIFTH OBJECT, i.e. (x5,y5) (ICTRS)

10 10 8 8 9 7 6 5 4 3 2 1 7 6 5 4 3 2 1 10 8 8 9 7 6 5 4 3 2 1 7 6 5 4 3 2 1 1 2 3 4 5 6 7 1 2 3 4 5 6 7 9 8 8 10 10 8 8 9 7 6 5 4 3 2 1 7 6 5 4 3 2 1 ctr2 ctr5 ctr3 ctr4 ctr1 ctr5 ctr3 ctr4 ctr2 ctr1 ctr3 ctr4 ctr5 ctr2 ctr1 ctr4 ctr3 ctr5 ctr2 ctr1 ctr4 ctr5 ctr2 ctr4 ctr5 ctr2 ctr4 ctr5 ctr2 ctr4 ctr5 ctr2 constraint) ACTRS (active internal ctr1 ctr2 ctr5 ctr3 ctr4 constraint) DCTRS EXTERNAL CONSTRAINT (non−overlapping)

DCTRS DCTRS DCTRS DCTRS DCTRS DCTRS DCTRS ctr3 ctr5 ctr2 ctr4 ctr1 ctr3 ctr2 ctr4 ctr1 ctr3 ctr5 ctr4 ctr1 ctr3 ctr5 ctr2 ctr1 ctr5 ctr2 ctr4 ctr1 ctr4 ctr1 ctr2 ctr3 ctr3 ctr5 ctr4 ctr1 ctr4 ctr5 ctr3 ctr2 ctr5 ctr2 ctr4 ctr5 ctr2 ctr1 ctr3 x5 in 1..8, y5 in 1..8, y5<>7 x4 in 7..7, y4 in 1..1 x3 in 2..4, y3 in 8..9 x2 in 4..4, y2 in 6..6 x1 in 1..4, y1 in 2..4 VARIABLES ctr5 ctr2 ctr4 ctr3 ctr5 ctr2 ctr1 ACTRS CONFLICT CONFLICT ACTRS CONFLICT ACTRS CONFLICT ACTRS CONFLICT ACTRS CONFLICT ACTRS ACTRS CONFLICT ACTRS CONFLICT

SWEEP POINT: c=(1,1) SWEEP POINT: c=(1,3)

SWEEP POINT: c=(1,7) SWEEP POINT: c=(1,8) SWEEP POINT: c=(3,1)

SWEEP POINT: c=(3,4) SWEEP POINT: c=(3,7)

(A) (B) (C) (D) (E) (F) (G) (H) (I) (J) (K) (L) 1 2 3 4 5 6 7 1 2 3 4 5 6 7 9 8 8 10 (delayed internal 10 8 8 9 7 6 5 4 3 2 1 7 6 5 4 3 2 1 1 2 3 4 5 6 7 1 2 3 4 5 6 7 9 8 8 10 10 8 8 9 7 6 5 4 3 2 1 7 6 5 4 3 2 1 1 2 3 4 5 6 7 1 2 3 4 5 6 7 9 8 8

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

00000000

11111111

r3 r3 c=(3,8) r₁ 4 5 r 3 r r₃ 3 r 3 r r₂ 1 r r₁ r₁ r

Fig. 4. Illustration of the lexicographic sweep algorithm for adjusting the minimum value of the

(12)

PROCEDURE

FilterCtrs(k, O, S, C) : bool

1: nonfix

← true

// fixpoint not yet reached

2: while nonfix do

3:

nonfix

← false

// assumes no filtering will be done

4:

for all o

∈ O do

5:

I

←

S

e∈C

GenOutboxes(e, o, o.sid)

// build the set of outboxes on o

6:

I

← I ∪

S

0≤d<k

possible outboxes corresponding to holes in o.x

[d]

7:

for d

← 0 to k − 1 do

8:

if

¬PruneMin(o, d, k, I) ∨ ¬PruneMax(o, d, k, I) then

9:

return false

// no feasible origin

10:

else if o.x was pruned then

11:

nonfix

← true

// fixpoint not yet reached

12:

end if

13:

end for

14:

end for

15: end while

16: return true

// feasible origin

PROCEDURE

PruneMin(o, d, k, I) : bool

1: b

← true

// b

= true while we have not failed

2: c

← o.x

// initial position of the point

3: n

← o.x + 1

// upper limits+1 in the different dimensions

4: while b

∧ ∃f ∈ I | c ∈ f do

5:

n

← min(n, f.t + f.l)

// update vector n according to an outbox f containing c

6:

b

← false

// no new point to jump to yet

7:

for j

← k − 1 downto 0 do

8:

j

′

_{← (j + d) mod k}

_{// rotation wrt. d, k}

9:

c[j

′

_{] ← n[j}

′

_]

// use vector n to jump

10:

n[j

′

_{] ← o.x[j}

′

_{] + 1}

_{// reset component of n to maximum value}

11:

if c

[j

′

_{] ≤ o.x[j}

′

_{] then}

12:

b

← true

// jump target found

13:

j

← 0

// exit for loop

14:

else

15:

c[j

′

_{] ← o.x[j}

′

_]

_{// reset component of c, for exhausted a dimension}

16:

end if

17:

end for

18: end while

19: if b then

20:

o.x[d] ← max(o.x[d], c[d])

21: end if

22: return b

Algorithm

1:

FilterCtrs is the main filtering algorithm associated with

geost(k, O, S, C), where k, O, S and C resp. correspond to the number of

di-mensions, to the objects, to the shapes and to the external geometrical constraints.

PruneMin adjusts the lower bound of the d

th

coordinate of the origin of object o

where I is the set of outboxes associated with object o (since

PruneMax is similar

to

PruneMin it is omitted). The given fixpoint loop is an over-simplification. The

implementation maintains a set of objects that need filtering. Whenever an object o is

pruned, all non-fixed objects connected to o by an external constraint are added to this

set. When the set becomes empty, the fixpoint is reached.

(13)

transfer to the list of active outboxes all outboxes for which the first lexicographically smallest

infeasible point is lexicographically greater than or equal to the current position of the sweep

c

= (1, 1) (i.e., constraint ctr

1

= outbox([1, 1], [2, 2])). We then search through the list of

active constraints (represented on the figure by a box with the legend ACTRS on top of it) the

first constraint for which c

= (1, 1) is infeasible. In fact, since ctr

1

is infeasible (represented on

the figure by a box with the legend CONFLICT on top of it) we compute the vector f

= (3, 3)

that tells how to get the next potentially feasible point in the different dimensions. Consequently

the sweep moves to the next position

(1, 3) (see Part (F)) and the process is repeated until we

finally find a feasible point for all outboxes (i.e., point

(3, 8) in Part (L)). Note that, when the

lexicographically largest infeasible point associated with an active outbox is lexicographically

less than the current position of the sweep, we remove that constraint from the list of active

outboxes. This is for instance the case in Part (I), where we remove constraint ctr

3

from the list

of active outboxes (since its lexicographically largest infeasible point

(2, 8) is lexicographically

less than the position of the sweep c

= (3, 1)).

Complexity. Rather than analysing the complexity of the geost kernel for a fixed k,

which depends both on the type of each external constraint (i.e., the complexity of a

given external constraint for generating all its corresponding outboxes as well as their

number), we rather focus on

PruneMin for adjusting the minimum value of the d

th

coordinate of the origin of an object. Assuming that the maximum number of outboxes

is equal to n we give an upper bound on the maximum number of jumps of

PruneMin

(i.e., the maximum number of times the sweep is moved).

First note that we always jump to an upper border (+1) of an outbox (i.e., see line 5

of

PruneMin) or that we reset some coordinates of the sweep to its minimum value (i.e.,

see line 15 of

PruneMin). Consequently, all coordinates of the sweep are always equal

to an upper border (+1) of some outboxes or to a minimum possible value. Since we

want to evaluate the maximum number of jumps, let us assume that for every dimension

d (0 ≤ d < k) the upper limits of all the n outboxes are distinct. Having this in

mind we can construct a maximum of

(n + 1)

k

points. Even if we found a systematic

construction where this number of jumps is reached, the performance evaluation of

Section 5 indicates that we can handle a reasonable number of objects for k

= 2, 3, 4.

From a memory consumption point of view, the algorithm only records the

coor-dinates of the sweep from one invocation to the next, in order not to restart the search

from scratch (i.e.,

2k points for each object).

Handling Time. Given an object o

∈ O of geost, the sweep algorithm that we have

introduced in the previous section can be easily adapted to handle the start in time

o.start , duration in time o.duration and end in time o.end . Beside maintaining bound

consistency for the constraint o.end

= o.start + o.duration, we add an extra time

dimension to the geometric coordinates of object o. Roughly, this new time coordinate

corresponds to o.start resp. o.end depending on whether we are adjusting the minimum

or maximum.

Handling Polymorphism. In order to handle the fact that objects can have several

potential shapes we modify the previous algorithm in the following way. For adjusting

the minimum value of the coordinate of the origin of an object that has more than one

shape we call the sweep algorithm for each potential shape of the object (i.e., for each

(14)

value of its shape variable). Then we take the smallest minimum value obtained (i.e.,

we use constructive disjunction) and prune the shape variable of an object if we did not

find any feasible point for a given potential shape of that object.

Other Internal Constraints. The standard representation of geometrical constraints

given in Section 3 is an over-simplification. For some constraints, e.g. distance

con-straints, outboxes are not a suitable representation, as the set of forbidden coordinates

cannot be covered by a small number of boxes. Therefore, the constraint kernel

inter-nally rather handles internal constraints, with an appropriate internal API. In order to

have a compact representation which can be used efficiently by the geometrical kernel a

set of infeasible points is in fact defined implicitly by providing the following functions:

–

LexInfeasible(ictr , minlex , d, k, o) : (found , p) when minlex = true

(respec-tively false), returns the smallest (respec(respec-tively largest) infeasible lexicographical

point p associated with the internal geometrical constraint ictr (according to the

fact that we prune the d

th

coordinate of the origin of object o, i.e., the ordering

among the different dimensions is d,

(d + 1) mod k, . . . , (d − 1) mod k)

com-patible with the domains of the coordinates of the origin of o. If no such point

exists, found is set to false (otherwise found is set to true).

–

IsFeasible(ictr , min, d, k, o, c) : (feasible , f ) sets feasible to true if point c is

feasible according to the internal constraint ictr ; if this is not the case, sets feasible

to false, and computes the forbidden region f according to the fact that we prune

the minimum (min

= true) or the maximum (min = false) value of the d

th

coordinate of o: we first maximise the size of f in dimension

(d − 1) mod k, then

maximise the size of f in dimension

(d − 2) mod k and so on until we reach the

most significant dimension d. Part (A) (respectively Part (B)) of Figure 5 illustrates

the computation of the forbidden region f in the context of k

= 2 and d = 0

(respectively d

= 1).

–

CardInfeasible(ictr , k, o) : (n) returns an estimation of the number n of infeasible

points for the origin of object o under the assumption that constraint ictr holds. This

information is used as a heuristics for ordering the internal constraints checked by

the geometrical kernel.

Algorithm 1 can be adapted to this API: (1) The variable I initialised at lines 5 and 6

of

FilterCtrs would be a set of internal constraints instead of outboxes. (2) Line 4 of

PruneMin would use IsFeasible.

5 Performance Evaluation

We evaluate the implementation

6

_{of the geost kernel from three perspectives:}

Wanting to measure the speed and the scalability of the sweep algorithm for

find-ing a first solution on loosely constrained placement problems (i.e., 20% spare space),

we generated one set of random problem instances of m k-dimensional boxes for

6

The experiments were run in SICStus Prolog 4 compiled with gcc -02 version 4.0.2 on a 3GHz

Pentium IV with 1MB of cache.

(15)

most significant

dimension

(A)

(B)

most significant dimension

f.max[1] f.min[1])=c[1] f.max[0] f.min[0]=c[0] f.min[0]=c[0] f.max[0] f.min[1]=c[1] f.max[1]

dimension

least significant

o.x[0] o.x[0] o.x[1] o.x[1] o.x[0] o.x[0] o.x[1] o.x[1]

infeasible points of the largest lexicographical infeasible point infeasible largest box f

infeasible largest box f (forbidden region)

infeasible points of the largest lexicographical infeasible point (forbidden region) infeasible point smallest lexicographical internal constraint smallest lexicographical infeasible point internal constraint

Fig. 5. Illustration of an internal geometrical constraint and of how to compute the forbidden

(16)

k

∈ {2, 3, 4} involving t ∈ {1, 16, 256, 1024} distinct types of boxes, and for m ∈

{1024, 2048, . . . , 262144}. The results for k = 2 are shown in Fig. 6 (top left) and

indicates that the approach is sensible to the number of distinct types of boxes. It can

typically pack

1024 2D, 3D and 4D distinct boxes in at most 200 msec. The longest

time,

13694 seconds (close to 4 hours), was obtained for packing 262144 4D

paral-lelepipeds (over

1 million domain variables) with a memory consumption of 351MB.

Wanting to get an idea of the performance of the geost kernel on very tight

place-ment problems (i.e., 0% spare space), we considered the perfect squared squares

prob-lem [1, 5] as well as the 3D pentominoes probprob-lem [6]:

– A perfect squared square of order n is a square that can be tiled with n smaller

squares where each of the smaller squares has a different integer size. We used the

data available (i.e., the size of the small squares to pack) from the catalogue [7] and

tested the corresponding

207 instances. The labelling strategy is roughly to repeat

the following, first for the x dimension, then for the y dimension:

1. Find the smallest position where some square can be placed.

2. Find a square to place in that position.

– Pentominoes are pieces made of

5 connected unit cubes laid on a plane surface.

Their shapes look like the

12 letters F , I, L, P , N , T , U , V , W , X, Y and Z. We

considered the problem of finding the different ways of putting

12 distinct shapes

that can be reflected and rotated in a box having a volume of

60 unit cubes. Our

labelling strategy is roughly to repeat the following:

1. Find a slot in the space that has not yet been filled by some piece.

2. Find a piece that can fill that slot.

Fig. 6 (top right) and Table 1 respectively report, for the squared squares and the

pentominoes problems, the time and number of backtracks for exploring all the search

space

7

without breaking any symmetry. For the squared squares problems the

maxi-mum time of

1585 seconds was spent on problem 48; on the other hand, 148 problems

were completely solved within

60 seconds. For the 3D pentomino packing instances,

performance results for comparison can be found in [6]. However, they stop the search

when the first 100 solutions have been found, so the results are only partly comparable.

Finally, wanting to compare the geost kernel with a recent exact state of the art

method for the 2D orthogonal packing problem [4], we reused the benchmarks

pro-posed by Clautiaux et al. [8]. This is a feasibility problem which consists in determining

whether a set of rectangles that cannot be rotated, can be packed or not into a rectangle

of fixed size. In these instances the discrepancy between the sum of the areas of the

rectangles to pack and the area of the big rectangle vary from 0% to 20%. We have

41 instances involving between

10 and 23 rectangles. Moreover, from these 41 instances,

26 instances are not feasible. In order to break symmetries between multiple rectangles

of the same shape we added lexicographic ordering constraints. All x coordinates were

labelled followed by all y coordinates, by decreasing rectangle size. Values were tried

by increasing value. Fig. 6 (bottom) compares our results with the ones reported in [4].

Note that the sequence order for the curves differs, since the instances of each curve

(17)

are ordered by increasing y value. We solved all instances and are comparable with [4],

although

8 instances are much easier for [4] and 10 instances are much easier for us.

Note that for the last three problems (i.e., Squared Squares, Pentominoes and 2D

orthogonal packing) extra filtering algorithms mostly based on cumulative relaxation

were integrated within our kernel. Since this report focusses on the constraint kernel

and because of space limitations these methods were not detailed.

6 Related Work and Future Directions

The rectangles packing problem has been studied by Clautiaux et al. [8, 4] where

schedul-ing-based reasoning is used [1]. The use of sweep algorithms in constraint filtering

algorithms was introduced in [3] and applied to the non-overlapping 2D rectangles

con-straints. This report generalizes and extends that work in several ways.

– The 2D sweep is generalized to a lexicographic sweep, independent of the number

of dimensions.

– The notion of forbidden regions for non-overlapping rectangles is generalized to

necessary conditions for general geometric constraints.

The idea of generating necessary conditions is reminiscent of indexicals [9], a.k.a.

projection constraints [10]. An indexical for a constraint c

(x

1

, . . . , x

n

) computes a

unary constraint on a single variable x

i

, i.e. a set S of values such that c

⇒ x

i

∈ S,

in reaction to domain changes in x

1

, . . . , x

n

. The constraint kernel then immediately

enforces x

i

∈ S. Our kernel generalizes this in two ways:

– We compute necessary conditions in the form of k-dimensional forbidden regions.

– We treat all such forbidden regions, for a selected object and shape, in one run of

the sweep algorithm. Projecting a single forbidden region on one coordinate often

does not yield any pruning, whereas considering the union of forbidden regions is

much more effective.

Dal Pal`u et al. in [11] proposed a constraint solver specialized for 3D discrete

do-mains. Their solver was targeted to the study of problems in molecular, chemical and

crystal structures. Our work, however, remains in the setting of mainstream finite

do-main constraint systems, whereas our kernel internally handles k-dimensional objects.

Even though the geost kernel has been designed over discrete domains, it could

rather easily be extended to continuous domains with the coordinates of the objects

approximated by the floating-point numbers F . Since switching from N to F may cause

rounding errors at this level, the sweep algorithm needs to handle these rounding errors

when moving the sweep out of an outbox constraint. If the projections of the forbidden

regions on all dimensions are intervals of real bounds we can proceed as follows. On

continuous domains, an outbox will have an very thin strip at the border where the

feasibility of the corresponding internal constraint is unknown. The region inside this

strip is strictly forbidden, and outside, the constraints is certainly satisfied. The outbox

must be computed including this strip, by taking lower and upper approximations of

the region’s coordinates. In that case, the solutions are guaranteed to be valid, but the

(18)

solver may not be complete, because it may (rarely) happen that the real forbidden

region allows positions that are forbidden by its approximation.

This research was conducted under the European Union project “Net-WMS”, a

ma-jor task of which is to study packing problems in warehouse management. In this

con-text, our constraint kernel is a step towards being able to capture a large set of packing

rules in a constraint programming setting. Future work involves extending our set of

external geometric constraints to include such packing rules.

7 Conclusion

The main contribution of this report is a geometrical constraint kernel for handling the

location in space and time of polymorphic k-dimensional objects subject to various

geometrical and time constraints. The constraint kernel is generic in the sense that one

of its parameters is a set of constraints on subsets of the objects. These constraints are

handled globally by the kernel.

We have presented a sweep algorithm for filtering the attributes of the objects.

Thank to its architecture, new geometric constraints can be plugged into this sweep

algorithm without modifying it. The strong point of this sweep algorithm is that it

con-siders all the geometrical constraints for a selected object and shape in one run. As

a first result, more deduction can be performed by combining sets of forbidden points

coming from multiple geometrical constraints. Secondly, it can handle within one single

constraint problems involving up several tens of thousands of objects without memory

consumption problems, which is often a weak point for constraint programming

envi-ronment. We have also shown that we could handle tight 2D or 3D placement problems,

which were traditionally solved by specific approaches.

Acknowledgements

This research was conducted under European Union Sixth Framework Programme

Con-tract FP6-034691 “Net-WMS”.

configuration backtracks (1st) time (1st) backtracks (all) time (all) solutions

20 × 3 × 1 1434 1840 47381 49740 8

15 × 4 × 1 290 560 888060 939060 1472

12 × 5 × 1 1594 1850 3994455 4112870 4040

10 × 6 × 1 111 260 9688985 10726810 9356

10 × 3 × 2 1267 2370 1203511 1778980 96

6 × 5 × 2 157 730 n/a n/a n/a

5 × 4 × 3 3567 14930 n/a n/a n/a

Table 1. Performance evaluation. 3D pentomino packing instances. Time in milliseconds. “n/a”

(19)

10 100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09 1000 10000 100000 1e+06 time m t=1 t=16 t=256 t=1024 100 1000 10000 100000 1e+06 1e+07 0 50 100 150 200 250 backtracks time 1 10 100 1000 10000 100000 0 5 10 15 20 25 30 35 40 45 time (we) time [4] 1 10 100 1000 10000 100000 0 5 10 15 20 25 30 35 40 45 backtracks (we) backtracks [4]

Fig. 6. Performance evaluation. Top left: scalability, t

∈ {1, 16, 256, 1024}. Top right: Perfect

Squared Squares, runtime and backtracks. Bottom: 2D Orthogonal Packing, runtime (left) and

backtracks (right). Time in milliseconds. In each curve, the instances are ordered by increasing y

value.

(20)

8 Annex: Benchmark Data and Code

8.1 Random Instances

:- use_module(library(lists)). :- use_module(library(random)). :- use_module(library(clpfd)). limit_rand(2, 10). % d=2: li=rand(1..33) limit_rand(3, 4). % d=3: li=rand(1..11) limit_rand(4, 3). % d=4: li=rand(1.. 6)

% T is the number of types of objects (i.e., the number of shapes) % K is the number of dimensions

% N is the number of objects

% generates a geost constraint with one non-overlapping constraint and search for a first solution % a run with: top(1, 2, 1024). top(1, 2, 2048). ... ...

% a run with: top(16, 2, 1024). top(16, 2, 2048). ... ... % a run with: top(256, 2, 1024). top(256, 2, 2048). ... ... % a run with: top(1024, 2, 1024). top(1024, 2, 2048). ... ... % a run with: top(1, 3, 1024). top(1, 3, 2048). ... ... % a run with: top(16, 3, 1024). top(16, 3, 2048). ... ... % a run with: top(256, 3, 1024). top(256, 3, 2048). ... ... % a run with: top(1024, 3, 1024). top(1024, 3, 2048). ... ... % a run with: top(1, 4, 1024). top(1, 4, 2048). ... ... % a run with: top(16, 4, 1024). top(16, 4, 2048). ... ... % a run with: top(256, 4, 1024). top(256, 4, 2048). ... ... % a run with: top(1024, 4, 1024). top(1024, 4, 2048). ... ... paper :-pow2(1024, N), t_param(T), k_param(K), top(T, K, N), fail. paper. t_param(1). t_param(16). t_param(256). t_param(1024). k_param(2). k_param(3). k_param(4). pow2(P, P). pow2(P, R) :-Q is P<<1, pow2(Q, R). top(T, K, N) :-Goal = top(T, K, N), T1 is T+1, N1 is N+1,

M1 is (N // T)+1, % number of objects plus one of a given type gen_ints(0, K, _Dimensions), % generates the list of dimensions gen_ints(1, N1, _ObjectsId), % generates the list of objects id length(Zeros, K), % generates a list with K ’0’ for the shifts domain(Zeros, 0, 0),

gen_shps(1, T1, K, Zeros, Shapes1, Vol), % generates T random shapes sort_shapes(Shapes1, Shapes2),

(21)

% place the origins of the objects in a k-dimensionsal box of size Limit % Limit is computed by majoring the needed volume by 20 percent and % by taking the 1/K root

Limit is integer(floor(exp(NeededVol+((NeededVol*20)//100),1/K)))+5, gen_objects(1, T1, 1, M1, K, Limit, Objects, Variables),

statistics(runtime, _), statistics_memory(Membase),

( diffn(Objects, Shapes2, [fixall(F)]), search_fixall(Variables, F),

statistics(runtime, [_,T2]), statistics_memory(Mem),

format(’goal=˜q CPU=˜d memory=˜d\n’, [Goal,T2,Mem-Membase]) -> true ; format(’goal=˜q failed\n’, [Goal])

). statistics_memory(Mem) :-garbage_collect, statistics(program, [P|_]), statistics(global_stack, [G|_]), statistics(local_stack, [L|_]), statistics(trail, [T|_]), statistics(choice, [C|_]), Mem is P+G+L+T+C. sort_shapes(Shapes1, Shapes4) :-tag_shapes(Shapes1, Shapes2), keysort(Shapes2, Shapes3), rebuild_shapes(Shapes3, Shapes4, 0). tag_shapes([], []). tag_shapes([shape(_,Off,Size1)|S1], [Size2-Off|S2]) :-negate_shape(Size1, Size2), tag_shapes(S1, S2). rebuild_shapes([], [], _). rebuild_shapes([Size1-Off|S1], [shape(J,Off,Size2)|S2], I) :-J is I+1, negate_shape(Size1, Size2), rebuild_shapes(S1, S2, J). negate_shape([], []). negate_shape([X|Xs], [Y|Ys]) :-Y is -X, negate_shape(Xs, Ys). search_fixall(_, 1). search_fixall([], _). search_fixall([X|Xs], F) :-indomain(X), search_fixall(Xs, F). gen_objects(T, T, _, _, _, _, [], []) :- !.

gen_objects(J, T, Oid, M1, K, Limit, Objects, Variables) :-J < T,

gen_objs(1, M1, Oid, J/*shape*/, K, Limit, Objs1, Vars1), J1 is J+1,

NextOid is Oid+M1-1,

gen_objects(J1, T, NextOid, M1, K, Limit, Objs2, Vars2), append(Objs1, Objs2, Objects),

append(Vars1, Vars2, Variables). gen_objs(M, M, _, _, _, _, [], []) :- !.

gen_objs(J, M, Oid, S, K, Limit, [object(Oid,S,Origins)|R], Variables) :-J < M,

gen_origins(0, K, Limit, Origins), J1 is J+1,

Oid1 is Oid+1,

(22)

append(Origins, Vars, Variables). gen_origins(K, K, _, []) :- !. gen_origins(J, K, Limit, [O|R])

:-J < K, O in 1..Limit, J1 is J+1,

gen_origins(J1, K, Limit, R). gen_shps(M, M, _, _, [], 0) :- !.

gen_shps(J, M, K, Zeros, [shape(J,Zeros,Sizes)|R], Volum) :-J1 is J+1,

limit_rand(K, Limit),

gen_sizes(0, K, Limit, Sizes, Vol), gen_shps(J1, M, K, Zeros, R, Vol1), Volum is Vol+Vol1. gen_sizes(K, K, _, [], 1) :- !. gen_sizes(J, K, L, [S|R], V) :-J < K, random(1, L, S), J1 is J+1, gen_sizes(J1, K, L, R, V1), V is S*V1. gen_ints(M, M, []) :- !. gen_ints(J, M, [J|R]) :-J < M, J1 is J+1, gen_ints(J1, M, R).

(23)

8.2 Perfect Squared Squares

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

:-findall(Size, posted_domains(_,Size), Sizes), sumlist(Sizes, Total),

format(’TOTAL domain size=˜d\n’, [Total]). posted_domains(ID, Total)

:-data(ID, _, Limit, Rev), reverse(Rev, Sizes),

constraints(ID, diffn, Xs, _Ys, _Sizes, _Limit, [-c,-k], _), print_square_domains(Xs, Sizes, DomSizes),

sumlist(DomSizes, Total),

format(’square instance=˜w space=˜d*˜d total domain size=˜d\n’, [ID,Limit,Limit,Total]). print_square_domains([], [], []).

print_square_domains([X|Xs], [S|Ss], [Sz|Szs]) :-fd_dom(X, Xd),

fd_size(X, Sz),

format(’square ˜dx˜d, domain=˜q, domain size=˜d\n’, [S,S,Xd,Sz]), print_square_domains(Xs, Ss, Szs).

runfirst(ID,Srch) :-Opt=[-c,-k],

statistics(runtime, _),

solve(ID, diffn, Srch, _, _, Opt), statistics(runtime, [_,Time]), fd_statistics(backtracks, B),

logit(’instance=˜w opt=˜w solutions=˜d backtracks=˜d time=˜d\n’, [ID,Opt,1,B,Time]). runall(ID,Srch)

:-Opt=[-c,-k],

statistics(runtime, _),

count(ID, diffn, Srch, N, Opt), statistics(runtime, [_,Time]), fd_statistics(backtracks, B),

logit(’instance=˜w opt=˜w solutions=˜d backtracks=˜d time=˜d\n’, [ID,Opt,N,B,Time]). runall(Srch) :-retractall(btr_time/2), data(ID, _, _, _), statistics(runtime, _), count(ID, diffn, Srch, N, [-c,-k]), statistics(runtime, [_,Time]), fd_statistics(backtracks, Btr),

logit(’instance=˜w opt=˜w solutions=˜d backtracks=˜d time=˜d\n’, [ID,[-c,-k],N,Btr,Time]), assertz(btr_time(Btr,Time)),

fail. runall(_)

:-findall(B-T, btr_time(B,T), Btr1), keysort(Btr1, Btr2),

findall(T-B, btr_time(B,T), Time1), keysort(Time1, Time2), tell(’squares.dat’), write_plot_data(Btr2, Time2, 0), told, keys_and_values(Time2, Times, Btrs), sumlist(Btrs, SumBtr), sumlist(Times, SumTime),

format(’TOTAL backtracks=˜d time=˜d\n’, [SumBtr,SumTime]). write_plot_data([], [], _).

write_plot_data([B-_|Bs], [T-_|Ts], I) :-J is I+1,

(24)

write_plot_data(Bs, Ts, J). runfirst(Srch)

:-Opt = [-c,-k], data(ID, _, _, _),

(solve(ID, diffn, Srch, _, _, Opt) -> true), fd_statistics(backtracks, B),

fd_statistics(resumptions, R),

logit(’instance=˜w opt=˜w backtracks=˜d resumptions=˜d\n’, [ID,Opt,B,R]), fail. runfirst(_). logit(Fmt, Args) :-prolog_flag(argv, [Logfile]), !, open(Logfile, append, S), format(S, Fmt, Args), close(S). logit(Fmt, Args) :-format(Fmt, Args).

count(ID, Type, Search, N, Opt)

:-findall(1, solve(ID,Type,Search,_,_,Opt), L), length(L, N).

squares(ID, Type, Search, Opt) :-solve(ID, Type, Search, Xs, Ys, Opt), writeq(Xs), nl,

writeq(Ys), nl, nl.

solve(ID, Type, Search, Xs, Ys, Opt)

:-constraints(ID, Type, Xs, Ys, Sizes, Limit, Opt, 0), search(Search, Xs, Ys, Sizes, Limit, 0).

search(primal, Xs, Ys, _, _Limit, _Fixall) :-labeling([bisect], Xs),

labeling([bisect], Ys).

search(dual, Xs, Ys, _, Limit, Fixall) :-dual_labeling(Xs, 1, Limit, Fixall), dual_labeling(Ys, 1, Limit, Fixall). search(lex, Xs, Ys, Ss, Limit, _)

:-transpose([Xs,Ys,Ss], Sqs), search_lex(Sqs, [], [1,1], Limit). search_lex([], _, _, _) :- !.

search_lex(Sqs, Done, P, Limit) :-member(Pl, Done),

in_square(P, Pl, P1, Limit), !, search_lex(Sqs, Done, P1, Limit). search_lex(Sqs, Done, P, Limit)

:-search_lex2(Sqs, Sqs1, Done, Done1, P, P1, Limit), search_lex(Sqs1, Done1, P1, Limit).

search_lex2([Sq|Sqs], Sqs, Done, [Sq|Done], P, P1, Limit) :-assign2(Sq, P), in_square(P, Sq, P1, Limit). search_lex2([Sq|Sqs], [Sq|Sqs1], Pl, Pl1, P, P1, Limit) :-Sq = [X,Y,_], lex_chain([P,[X,Y]], [op(#<)]), search_lex2(Sqs, Sqs1, Pl, Pl1, P, P1, Limit). assign2([X,Y,_], [U,V]) :-clpfd:’$fd_in_interval’(X, U, U, 1), clpfd:’$fd_in_interval’(Y, V, V, 0), clpfd:’$fd_evaluate_indexical’(RC, Global), clpfd:evaluate(RC, Global).

in_square(P, [X,Y,S], P1, Limit) :-P = [:-Px,:-Py],

(25)

X1 is X+S-1, Y1 is Y+S-1, Px in X..X1, Py in Y..Y1, ( Y1=:=Limit -> Px1 is Px+1, P1 = [Px1,1] ; Y2 is Y1+1, P1 = [Px,Y2] ). dual_labeling(_, _, _, 1) :- !. dual_labeling([], _, _, _) :- !. dual_labeling(L, I, Limit, Fixall)

:-dual_labeling(L, L1, I, Limit, J, Fixall), dual_labeling(L1, J, Limit, Fixall). dual_labeling([], [], _, J, J, _) :- !. dual_labeling([X|L1], L2, I, J0, J, Fixall)

:-( integer(X) -> dual_labeling(L1, L2, I, J0, J, Fixall) ; X #= I, dual_labeling(L1, L2, I, J0, J, Fixall) ; X #> I, fd_min(X, J1), J2 is min(J0,J1), L2 = [X|L3], dual_labeling(L1, L3, I, J2, J, Fixall) ).

constraints(ID, Type, Xs, Ys, Sizes, Limit, Opt, Fixall) :-generate_squares(ID, Xs, Ys, Sizes, Limit),

order_squares(Xs, Ys, Sizes),

% state_asymmetry(Xs, Ys, Sizes, Limit), state_no_overlap(Type, Xs, Ys, Sizes, Fixall), ( memberchk(-c, Opt) -> true

; cumulative(Xs, Sizes, Sizes, Limit), cumulative(Ys, Sizes, Sizes, Limit) ),

( memberchk(-k, Opt) -> true ; knapsack(0, Limit, Xs, Sizes),

knapsack(0, Limit, Ys, Sizes) ), true. knapsack(L, L, _, _) :- !. knapsack(I, L, Xs, Sizes) :-J is I+1, crossings(Xs, Sizes, Bs, J), scalar_product(Sizes, Bs, #=, L, [consistency(domain)]), knapsack(J, L, Xs, Sizes). crossings([], [], [], _). crossings([X|Xs], [S|Ss], [B|Bs], J) :-J0 is J-S+1, X in J0..J #<=> B, crossings(Xs, Ss, Bs, J). % for compatibility cumulative(Os, Ds, Hs, L) :-mktasks(Os, Ds, Hs, Tasks), cumulative(Tasks, [limit(L),global(true)]). mktasks([], [], [], []). mktasks([O|Os], [D|Ds], [H|Hs], [task(O,D,E,H,0)|Tasks]) :-E in 0..1000, mktasks(Os, Ds, Hs, Tasks).

generate_squares(ID, Xs, Ys, Sizes, Size) :-data(ID, _, Size, Rev),