PyCSP3
COP  difficult

# Problem Amaze

From Minizinc, Challenge 2012. Given a grid containing p pairs of numbers (ranging from 1 to p), connect the pairs (1 to 1, 2 to 2, …, $p$ to $p$) by drawing a line horizontally and vertically, but not diagonally. The lines must never cross. The number of lines must be minimized.

To build a COP (Constraint Optimization Problem) model, we need first to import the library PyCSP$^3$:

from pycsp3 import *


Then, we need some data. Here, we have a grid of size $5 \times 5$ with value 1 in cells at index (3,4) and (5,1), and value 2 in cells at index (2,2) and (4,2); here, p=2, and indexing is assumed to start at 1.

n, m = 5, 5

points = [
[(3,4), (5,1)],  # the two cells where value 1 is present
[(2,2), (5,5)]   # the two cells where value 2 is present
]


For representing a solution, we can fill up the grid with either value 0 (empty cell) or a line number (value from 1 to p).

We define a constant.

nPairs = len(points)  # number of pairs of points


We start our COP model by introducing an array $x$ of variables. Note the presence of a border. We use here a function that defines the specific domain of a variable at indexes (i,j) specifoed as parameters. For the border, the domain is singleton as it is {0}.

def domain_x(i, j):
return {0} if i in {0, n + 1} or j in {0, m + 1} else range(nPairs + 1)

# x[i][j] is the value at row i and column j (a boundary is put around the board).
x = VarArray(size=[n + 2, m + 2], dom=domain_x)


We can display the structure of the array, as well as the domains of the variables.

print("Array x: ", x)
for i in range(n+2):
print(f"Domain of variables on row {i}: {[x[i][j].dom for j in range(m+2)]}")

Array x:  [
[x[0][0], x[0][1], x[0][2], x[0][3], x[0][4], x[0][5], x[0][6]]
[x[1][0], x[1][1], x[1][2], x[1][3], x[1][4], x[1][5], x[1][6]]
[x[2][0], x[2][1], x[2][2], x[2][3], x[2][4], x[2][5], x[2][6]]
[x[3][0], x[3][1], x[3][2], x[3][3], x[3][4], x[3][5], x[3][6]]
[x[4][0], x[4][1], x[4][2], x[4][3], x[4][4], x[4][5], x[4][6]]
[x[5][0], x[5][1], x[5][2], x[5][3], x[5][4], x[5][5], x[5][6]]
[x[6][0], x[6][1], x[6][2], x[6][3], x[6][4], x[6][5], x[6][6]]
]
Domain of variables on row 0: [0, 0, 0, 0, 0, 0, 0]
Domain of variables on row 1: [0, 0..2, 0..2, 0..2, 0..2, 0..2, 0]
Domain of variables on row 2: [0, 0..2, 0..2, 0..2, 0..2, 0..2, 0]
Domain of variables on row 3: [0, 0..2, 0..2, 0..2, 0..2, 0..2, 0]
Domain of variables on row 4: [0, 0..2, 0..2, 0..2, 0..2, 0..2, 0]
Domain of variables on row 5: [0, 0..2, 0..2, 0..2, 0..2, 0..2, 0]
Domain of variables on row 6: [0, 0, 0, 0, 0, 0, 0]


Concerning the constraints, we start by setting the clues (fixed values) in the grid, by posting some constraints Intension.

satisfy(
# putting two occurrences of each value on the board
x[i][j] == v + 1 for v in range(nPairs) for i, j in points[v]
);


Interestingly, by calling the function solve(), we can check that the problem is satisfiable (SAT). We can also display the found solution. Here, we call the function values() that collects the values assigned to a specified list of variables. Because only some variables are involved in constraints, we obtain * in many cells (meaning that any value can be chosen).

if solve() is SAT:
print(values(x))

[
[*, *, *, *, *, *, *]
[*, *, *, *, *, *, *]
[*, *, 2, *, *, *, *]
[*, *, *, *, 1, *, *]
[*, *, *, *, *, *, *]
[*, 1, *, *, *, 2, *]
[*, *, *, *, *, *, *]
]


Next, we know that each fixed cell must have exactly one neighbour with the same value. We can enforce this with a group of constraints Count.

satisfy(
# each fixed cell has exactly one neighbour with the same value
Count([x[i - 1][j], x[i + 1][j], x[i][j - 1], x[i][j + 1]], value=v + 1) == 1
for v in range(nPairs) for i, j in points[v]
);


At this stage, we can display the internal representation of the posted constraints; this way, although a little bit technical, we can see that the constraints are correctly built. Note that the unary constraints have been collected together so as to form a constraint Instantiation, and that ‘eq’ stands for ‘equal to’.

print(posted())

instantiation(list:x[3][4] x[5][1] x[2][2] x[5][5], values:[1, 1, 2, 2])
count(list:[x[2][4], x[4][4], x[3][3], x[3][5]], values:[1], condition:(eq,1))
count(list:[x[4][1], x[6][1], x[5][0], x[5][2]], values:[1], condition:(eq,1))
count(list:[x[1][2], x[3][2], x[2][1], x[2][3]], values:[2], condition:(eq,1))
count(list:[x[4][5], x[6][5], x[5][4], x[5][6]], values:[2], condition:(eq,1))


We can run again the solver. Intersingly, we can see that building a valid solution is under progress.

if solve() is SAT:
print(values(x))

[
[*, *, *, *, *, *, *]
[*, *, 0, *, *, *, *]
[*, 0, 2, 0, 0, *, *]
[*, *, 2, 0, 1, 0, *]
[*, 0, *, *, 1, 0, *]
[0, 1, 1, *, 2, 2, 0]
[*, 0, *, *, *, 0, *]
]


When analysing the problem, one can find that any non-empty cell (i.e., any cell with a value different from 0) is such that if it is not an end-point then it has exactly two horizontal or vertical neighbours with the same value.

Interstingly, we can capture this rather complex constraint by a table that indicates the possible combinations of values for exactly 5 variables (forming a cross shape in the grid). When building the table, for simplicity, one can use the function ne() to stand for any value ’not equal’ to the specified parameter. The table is defined as follows:

T = ([(0, ANY, ANY, ANY, ANY)] +
[tuple(ne(v) if k in (i, j) else v for k in range(5))
for v in range(1, nPairs+1) for i, j in combinations(range(1, 5), 2)])


It permits to indicate that when the first value (for the cell situated at the center of the cross) of a tuple is not 0, exactly two other values in the tuple must be equal to this value. If we print the table, we can see the presence of unary conditions inside the so-called hybrid (smart) table. In the very near future, we shall let the user the opportunity to generate such tables when compiling.

print(T)

[(0, *, *, *, *), (1, ≠1, ≠1, 1, 1), (1, ≠1, 1, ≠1, 1), (1, ≠1, 1, 1, ≠1), (1, 1, ≠1, ≠1, 1), (1, 1, ≠1, 1, ≠1), (1, 1, 1, ≠1, ≠1), (2, ≠2, ≠2, 2, 2), (2, ≠2, 2, ≠2, 2), (2, ≠2, 2, 2, ≠2), (2, 2, ≠2, ≠2, 2), (2, 2, ≠2, 2, ≠2), (2, 2, 2, ≠2, ≠2)]


We can now post a group of constraints Extension after having identified the relevant free cells, those that are not intially fixed.

free_cells = [(i,j) for i in range(1, n + 1) for j in range(1, m + 1)
if (i, j) not in [p for pair in points for p in pair]]
satisfy(
# each free cell either contains 0 or has exactly two neighbours with its value
(x[i][j], x[i - 1][j], x[i + 1][j], x[i][j - 1], x[i][j + 1]) in T for i,j in free_cells
);


Let us have a look at the list of relevant cells.

print(free_cells)

[(1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (2, 1), (2, 3), (2, 4), (2, 5), (3, 1), (3, 2), (3, 3), (3, 5), (4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (5, 2), (5, 3), (5, 4)]


And let us print the very last posted constraint. Here, we can see that the hybrid table has been automatically compiled into a starred one.

print(posted()[-1])

extension(list:[x[5][4], x[4][4], x[6][4], x[5][3], x[5][5]], supports:(0,*,*,*,*)(1,0,0,1,1)(1,0,1,0,1)(1,0,1,1,0)(1,0,1,1,2)(1,0,1,2,1)(1,1,0,0,1)(1,1,0,1,0)(1,1,0,1,2)(1,1,0,2,1)(1,1,1,0,0)(1,1,1,0,2)(1,1,1,2,0)(1,1,1,2,2)(1,2,0,1,1)(1,2,1,0,1)(1,2,1,1,0)(1,2,1,1,2)(1,2,1,2,1)(2,0,0,2,2)(2,0,2,0,2)(2,0,2,1,2)(2,0,2,2,0)(2,0,2,2,1)(2,1,0,2,2)(2,1,2,0,2)(2,1,2,1,2)(2,1,2,2,0)(2,1,2,2,1)(2,2,0,0,2)(2,2,0,1,2)(2,2,0,2,0)(2,2,0,2,1)(2,2,2,0,0)(2,2,2,0,1)(2,2,2,1,0)(2,2,2,1,1))


We can run the solver.

if solve() is SAT:
print(values(x))

[
[*, 0, 0, 0, 0, 0, *]
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 2, 2, 2, 2, 0]
[0, 0, 0, 1, 1, 2, 0]
[0, 0, 0, 1, 0, 2, 0]
[0, 1, 1, 1, 0, 2, 0]
[*, 0, 0, 0, 0, 0, *]
]


Concerning optimization, we want to minimize the sum of values assigned to $x$.

minimize(
Sum(x)
);


We can run again the solver, with this optimization task. Note that we need to check that the status returned by the solver is now OPTIMUM.

if solve() is OPTIMUM:
print(values(x))
print("Bound: ", bound())

[
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 2, 2, 2, 2, 0]
[0, 0, 0, 1, 1, 2, 0]
[0, 0, 0, 1, 0, 2, 0]
[0, 1, 1, 1, 0, 2, 0]
[0, 0, 0, 0, 0, 0, 0]
]
Bound:  20


Finally, we give below the model in one piece, with two minor changes: we use ExactlyOne as a specialized version of Count, and we use x.beside(i,j) to denote the set of cells that are aside the cell with coordinates (i,j). Here the data is expected to be given by the user (in a command line).

from pycsp3 import *

n, m, points = data        # points[v] gives the pair of points for value v+1
nPairs = len(points)  # number of pairs of points
free_cells = [(i, j) for i in range(1, n + 1) for j in range(1, m + 1)
if  (i, j) not in [tuple(p) for pair in points for p in pair]]

# x[i][j] is the value at row i and column j (a boundary is put around the board).
x = VarArray(size=[n + 2, m + 2], dom=lambda i, j: {0} if i in {0, n + 1} or j in {0, m + 1} else range(nPairs + 1))

T = ({(0, ANY, ANY, ANY, ANY)}
|{tuple(ne(v) if k in (i, j) else v for k in range(5)) for v in range(1, nPairs + 1)
for i, j in combinations(range(1, 5), 2)})

satisfy(
# putting two occurrences of each value on the board
[x[i][j] == v + 1 for v in range(nPairs) for i, j in points[v]],

# each cell with a fixed value has exactly one neighbour with the same value
[ExactlyOne(x.beside(i,j), value=v + 1) for v in range(nPairs) for i, j in points[v]],

# each empty cell either contains 0 or has exactly two neighbours with the same value
[(x[i][j], x[i - 1][j], x[i + 1][j], x[i][j - 1], x[i][j + 1]) in T for i, j in free_cells]
)

minimize(
Sum(x)
)