CSP easy

Problem BIBD

From CSPLib (Problem 28): ``Balanced Incomplete Block Design (BIBD) generation is a standard combinatorial problem from design theory, originally used in the design of statistical experiments but since finding other applications such as cryptography. It is a special case of Block Design, which also includes Latin Square problems. BIBD generation is described in most standard textbooks on combinatorics. A BIBD is defined as an arrangement of $v$ distinct objects into $b$ blocks such that each block contains exactly $k$ distinct objects, each object occurs in exactly $r$ different blocks, and every two distinct objects occur together in exactly $\lambda$ blocks. Another way of defining a BIBD is in terms of its incidence matrix, which is a $v$ by $b$ binary matrix with exactly $r$ ones per row, $k$ ones per column, and with a scalar product of $\lambda$ between any pair of distinct rows. A BIBD is therefore specified by its parameters $(v,b,r,k,\lambda)$.’’

An example of a solution for $(7,7,3,3,1)$ is:

1 1 0 0 1 0 
0 1 0 1 0 0 
0 1 1 0 0 1 
1 0 0 0 0 1 
0 0 0 1 1 1 
0 0 1 0 1 0 
1 0 1 1 0 0 

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

from pycsp3 import *

Then, we need some data. Here, we need five integers $v$, $b$, $r$, $k$, $l$ (for $\lambda$) for specifying a unique instance; possibly, $b$ and $r$ can be set to 0, so that these values are automatically computed according to a classical BIBD template.

v, b, r, k, l = 7, 7, 3, 3, 1

We start our CSP model by introducing a two-dimensional array $x$ of variables.

# x[i][j] is the value of the matrix at row i and column j
x = VarArray(size=[v, b], dom={0, 1})

We can display the structure of the array, as well as the domain of the first variable (note that all variables have the same domain).

print("Array x: ", x)
print("Domain of any variable: ", x[0][0].dom)

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 any variable:  0 1

Concerning the constraints, we first impose that the number of 1 per row is $r$. We use a constraint Sum as follows:

satisfy(
   # constraints on rows
   Sum(row) == r for row in x
);

By calling the function solve(), we can check that the problem is satisfiable (SAT). We can also display the values assigned to the variables in the solution that has been found.

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

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

Now, we impose that the number of 1 per column is $k$:

satisfy(
   # constraints on columns
   Sum(col) == k for col in columns(x)
);

We can check that we can find a solution with these two groups of constraints.

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

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

Finally, we have to post some kind of scalar constraints with respect to the value of $l$.

satisfy(
   # scalar constraints with respect to lambda
   row1 * row2 == l for row1, row2 in combinations(x, 2)
);

We do obtain a BIDB solution, when we execute:

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

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

Note that the number of solutions is $151,200$ (we do not try to find them here, because it is rather inefficient).

As a matter of fact, this problem has many symmetries. It is known that we can break variable symmetries by posting a lexicographic constraint between any two successive rows and any two successive columns. The reader is invited to reflect about the proof. For posting lexicographic constraints, we use the constraint LexIncreasing, with the named parameter matrix set to True since the constraint must be applied on each row and each column of the specified two-dimensional array. On the other hand, it is relevant to tag this constraint because it clearly informs us that it is inserted for breaking symmetries (and besides, tags may be exploited by solvers): tagging is made possible by putting in a comment line an expression of the form tag(), with a token (or a sequence of tokens separated by a white-space) between parentheses. We write:

satisfy(
    # Increasingly ordering both rows and columns  tag(symmetry-breaking)
    LexIncreasing(x, matrix=True)
);

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

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

Now, we can compute the number of solutions.

if solve(sols=ALL) is SAT:
    print("Number of solutions: ", n_solutions())

Number of solutions:  1

Finally, we give below the model in one piece. Here the data is expected to be given by the user (in a command line).

from pycsp3 import *

v, b, r, k, l = data 
b = (l * v * (v - 1)) // (k * (k - 1)) if b == 0 else b  
r = (l * (v - 1)) // (k - 1) if r == 0 else r  

# x[i][j] is the value of the matrix at row i and column j
x = VarArray(size=[v, b], dom={0, 1})

satisfy(
   # constraints on rows
   [Sum(row) == r for row in x],

   # constraints on columns
   [Sum(col) == k for col in columns(x)],

   # scalar constraints with respect to lambda
   [row1 * row2 == l for row1, row2 in combinations(x, 2)],
    
   # Increasingly ordering both rows and columns  tag(symmetry-breaking)
   LexIncreasing(x, matrix=True)
)