My Constraint Programming Blog: 12 more Gecode models. e.g. circuit, clique, and Hidato puzzle

Here are 12 more Gecode models and some comments. Note that more information and references to the problem is always presented is the model. These and other Geode models are collected at My Gecode page.

alldifferent_except_0.cpp: Global constraint of alldifferent_except_0 (decomposition), also the slightly more general alldifferent_except_c

This is one of my first test models how to create (a decomposition of) a global constraint as well as how to call a function. Note the use of imp (implication):
```
void alldifferent_except_c(Space& space, IntVarArray x, int c, IntConLevel icl = ICL_BND) {
  for(int i = 0; i < x.size(); i++) {
    for(int j = i+1; j < x.size(); j++) {
      post(space,
           tt(imp(x[i] != c && x[j] != c, // =>
                  x[i] != x[j])),
           icl
           );
    }
  } 
}
```
Just for fun, the model also has the constraint that there should be at least one 0.

discrete_tomography.cpp: Discrete tomography

For more about discrete tomography, see for example Christoph Dürr's Discrete Tomography and the comments in the ECLiPSe code tomo.ecl.txt.

The coding of the model was quite straightforward. The main lesson learned here is to be able to state (from the command line) which problem to solve. I borrowed the principle used in Gecode's nonogram example.

First in the model is the global variables for accessing the specific problem:

namespace {
  extern const int* problems[];         // List of problems
  extern const unsigned int n_examples; // Number of specifications
}

And last in the file is the specifications of the problems:

namespace {
  // ...
  // Specification for problem 1
  const int p1[] =
    { 11, 12,
      // Rows sums
      0,0,8,2,6,4,5,3,7,0,0,
      // Column sums
      0,0,7,1,6,3,4,5,2,7,0,0
    };
  
  // ...

To access the problem the following is used (slightly compressed)

class DiscreteTomography : public Example {
protected:
  const int* prob;  // Problem to be used
  IntVarArray x;     // Fields of matrix
  int rows(void) const { return prob[0]; } // Return rows of matrix
  int cols(void) const { return prob[1]; } // Return columns of matrix
public:
  DiscreteTomography(const SizeOptions& opt) 
    : prob(problems[opt.size()]),
      x(*this, rows()*cols(), 0, 1)
  {
    // ...

labeled_dice.cpp: Labeled dice problem

Jim Orlin wrote about this puzzle in Colored letters, labeled dice: a logic puzzle and I presented a Comet solution (labeled_dice.co) in Comet: regular constraint, a much faster Nonogram with the regular constraint, some OPL models, and more.

The Gecode model use one element constraint and a couple of count constraints. It was simpler than I first thought.
building_blocks.cpp: Building Blocks (Dell Logic Puzzle)

This problem almost the same as labeled_dice (see above) and should perhaps be incorporated into that model.
young_tableaux.cpp: Young tableaux and partitions

Young tableaux are described in the Wikipedia article Young tableau and MathWorld article YoungTableau, and is one of my standard comparison models.

The tricky part here was to calculate the partition (from the Young tableau), especially the following (MiniZinc/Comet) where p is the partition element, and x is the Young tableau matrix. p[i] == sum(j in 1..n) (x[i,j] <= n) The corresponding Gecode was:
```
for(int i = 0; i < n; i++) {
  // p[i] == sum(j in 1..n) (x[i,j] <= n)
  IntVar t(*this, 0, n+1);
  count(*this, m.row(i), n+1, IRT_EQ, t, opt.icl());
  post(*this, p[i] == n-t, opt.icl());
}
```
Not very advanced, but it took some time.
set_covering_deployment.cpp: Set covering deployment

Set covering deployment is presented in the MathWorld Set Covering Deployment and uses a small example of placing Roman armies.

Well, the hardest part was the second requirement: "There should always be an backup army near every city". See the model how I implemented that.
subset_sum.cpp: Subset sum problem

This problem is from Katta G. Murty Optimization Models for Decision Making (PDF) page 340:
Example 7.8.1

A bank van had several bags of coins, each containing either 16, 17, 23, 24, 39, or 40 coins. While the van was parked on the street, thieves stole some bags. A total of 100 coins were lost. It is required to find how many bags were stolen.
Comments The subset sum constraint is simply this constraint where x is the values to choose from and tot is the sum:
```
 linear(space, values, x, IRT_EQ, tot,  icl);
```
The model is slightly more general, as the total (to sum) can be set from the command line.

sonet.cpp: SONET problem

This is a standard problem in constraint programming. My solution is based on the ESSENCE' (Andrea Rendl's Tailor) model sonet_problem.eprime.

The biggest problem was to translate the following requirement if there is a demand between 2 nodes, then there has to exist a ring, on which they are both installed which use an exists construct in the ESSENCE' model (as well as my MiniZinc model sonet_problem.mzn. The correspondent Gecode code is:

for(int client1 = 0; client1 < n; client1++) {
  for(int client2 = client1 + 1; client2 < n; client2++) {
    if(demand[client1*n+client2] == 1) {
      // ESSENCE' code:
      //   exists ring : int(1..r) .
      //      rings[ring,client1] + rings[ring, client2] >= 2)

      // here we use reification: 
      BoolVarArray bs(*this, r, 0, 1);
      for(int ring = 0; ring < r; ring++) {
        // does this ring satisfy the condition?
        // note: the "~" is needed
        bs[ring] = post(*this, ~(rings_m(client1, ring) + rings_m(client2,ring) >= 2), opt.icl());
      }
      // at least one success is needed
      linear(*this, bs, IRT_GQ, 1, opt.icl());
    }
  }
}

I really like that Gecode supports this higher level of stating logical operations and reifications.

hidato.cpp: Hidato puzzle

Hidato is yet another grid puzzle. For more about Hidato see the following sources (and games to play):

www.shockwave.com (Shockwave)
www.hidato.com

The first source explains the problem as:

Puzzles start semi-filled with numbered tiles. The first and last numbers are circled. Connect the numbers together to win. Consecutive number must touch horizontally, vertically, or diagonally.

Conceptually the problem can be modeled as follows (which was done in the MiniZinc model hidato.mzn), e.g. one forall loop and four exists loops. It works but is inefficient:

  % ...
  forall(k in 1..r*c-1) (
     exists(i in 1..r, j in 1..c) (
        k = x[i, j] % fix this k /\
        exists(a, b in  {-1, 0, 1}
          where 
          i+a >= 1 /\ j+b >=  1 /\
          i+a <= r /\ j+b <= c
          /\ not(a = 0 /\ b = 0) 
        ) 
        (
          % find the next k
          k + 1 = x[i+a, j+b]  
        )
     )
  )

I started to follow this logic in Gecode with a lot of boolean arrays and reifications, but after some thoughts here is what was used instead.

    for(int k = 1; k < n*n; k++) {
      IntVar i(*this, 0, n-1);
      IntVar j(*this, 0, n-1);
      IntVar a(*this, -1, 1);
      IntVar b(*this, -1, 1);

      // 1) First: fix this k, i.e.
      //        k == board[i*n+j]
      IntVar inj(*this, 0,n*n-1);
      post(*this, inj == i*n+j, opt.icl());
      element(*this, board, inj, k, opt.icl());

      // 2) Then, find the position of the next value, i.e.
      //     k + 1 == board[(i+a)*n+(j+b)]
      IntVar ai(*this, 0, n-1); // this takes care of the borders of the matrix
      IntVar jb(*this, 0, n-1); // ibid.
      post(*this, a+i == ai,opt.icl());
      post(*this, j+b == jb,opt.icl());
      IntVar ai_plus_jb(*this, 0,n*n-1);
      post(*this, ai_plus_jb == ai*n+jb, opt.icl());
      element(*this, board, ai_plus_jb, k+1, opt.icl());
    }

Note that the borders of the matrix is handled by the domains of the IntVar's ai and jb.

circuit_orbit.cpp: Global constraint circuit (decomposition) using permutation orbits

circuit is a global constraint I have been fascinated by for a time, and I first implemented a decomposition in MiniZinc (circuit_test.mzn). Gecode has a builtin version of circuit, but I thought it could be fun to implement this decomposition as well.

It uses some observation I noted about permutation orbits when playing around with another problem.

Also, it sets the parameter c_d (recomputation commit distance) to size*2 in order to save memory. For example, with the default value of requires 600Mb memory for the problem size=500, but with c-d 1000 requires only 12Mb (for more about this, see the discussion and comments in de Bruijn sequences in Gecode (and other systems)).
nadel.cpp: Nadel's construction problem ("soft" version)

This a construction problem attributed by Rina Dechter "Constraint Processing" (page 5) to B.A. Nadel's article "Constraint satisfaction algorithms" (1989). I not sure about the solution, since I have not been able to find Nadel's article.

I haven't found any solution that satisfies all the constraints presented in the problem (which may indicate that my model is wrong).

However this can be used as a study of reifications: this "soft" version minimized the broken constraints, and generates 28 different solutions with 1 broken constraint. The broken constraints are either
- steep_slopes constraints, or
- near_dump constraints.
clique.cpp: Global constraint clique (decomposition), maximize the cardinality of cliques in a graph

This decomposition of the global constraint clique for a graph g uses the (obvious) observation that if two different nodes is in the same clique (s) then there must be a connection between the two nodes in the boolean (channeled) representation of s. This actually seems to be enough for use as a constraint (which is an example of the declarative nature of constraint programming). In MiniZinc parlance this is formulated in the model clique.mzn as:
```
link_set_to_booleans(x, x_bool) /\
forall(i,j in 1..n where i != j) (
   (x_bool[i] /\ x_bool[j]) -> ( g[i,j] > 0)
) /\
c = card(s)
```
The Gecode variant is slightly more complicated, mainly because the imp function don't accept IntArgs as a second argument.
```
// x_bool is a boolean representation of the clique set x
BoolVarArray x_bool(*this, n, 0, 1);
channel(*this, x_bool, x); // channel the clique to a boolean array
for(int i = 0; i < n; i++) {
  for(int j = 0; j < n; j++) {
    if (i != j) {
      BoolVar b(*this, 0, 1);
      if (g[i*n+j] > 0) {
        post(*this, b == 1, opt.icl());
      } else {
        post(*this, b == 0, opt.icl());
      }
      post(*this, tt(imp(
                        x_bool[i] == 1 && x_bool[j] == 1, // ->
                        b)),
          opt.icl());
     }
  }
}
```
The model contains the small example from Global constraint catalog: clique, a graph of 10 nodes, a bigger random graph of 100 nodes, and a graph randomizer. These examples is choosen by re-commenting the source code.

The included random graph of size 100 finds a maximum clique (size 6) in about 2 seconds. and generates the four size 6 cliques in about the same time.