Main

September 11, 2010

Results of MiniZinc Challenge 2010

The results of MiniZinc Challenge 2010 has been published: MiniZinc Challenge 2010 Results:
The entrants for this year (with their descriptions, when provided):

In addition the challenge organisers entered the following FlatZinc implementations:

  • Chuffed (description). A C++ FD solver using Lazy clause generation.
  • Fzntini. Translates to SAT, uses TiniSAT.
  • G12/FD. A Mercury FD solver (the G12 FlatZinc interpreter's default solver).
  • G12/CPLEX. Translates to MIP, uses CPLEX12.1.

As per the challenge rules, these entries are not eligible for prizes, but do modify the scoring results. Furthermore, entries in the FD search category (Gecode, JaCoP, Chuffed, G12/FD) were automatically included in the free search category, while entries in the free search category (Fzn2smt, Fzntini, SCIP, CPLEX and promoted FD entries) were automatically included in the parallel search category.

I'm really curious about Chuffed (which I have not tested). As far as I know, it is not public available.

Results

This year the results are not in fixed result lists. Instead there is a Javascript application where one can choose different combinations of solvers and problems.

Update some hours later: I have been informed that my summaries was not correct so I have removed them to not confuse anyone. Sorry for any confusions. Please see the Results page.

See also

The results of the earlier MiniZinc Challenges:
* MiniZinc Challenge 2009 Results. Results of the MiniZinc challenge 2009.
* MiniZinc Challenge 2008 Results. Results of the MiniZinc challenge 2008

August 31, 2010

Nontransitive dice, Ormat game, 17x17 challenge

Here are some new models done the last weeks and some older not mentioned before. (I have also spent some time with Numberjack, a Python based solver, and will blog about this later on.)

Nontransitive dice

MiniZinc: nontransitive_dice.mzn
Comet: nontransitive_dice.co

(Sidenote: This was triggered by a problem in Numberjack's Tutorial.)

Nontransitive dice is presented at Wikipedia as:
A set of nontransitive dice is a set of dice for which the relation "is more likely to roll a higher number" is not transitive. See also intransitivity. This situation is similar to that in the game Rock, Paper, Scissors, in which each element has an advantage over one choice and a disadvantage to the other.
And then gives the following example:
Consider a set of three dice, A, B and C such that
* die A has sides {2,2,4,4,9,9},
* die B has sides {1,1,6,6,8,8}, and
* die C has sides {3,3,5,5,7,7}.

Then:
* the probability that A rolls a higher number than B is 5/9 (55.55 %),
* the probability that B rolls a higher number than C is 5/9, and
* the probability that C rolls a higher number than A is 5/9.

Thus A is more likely to roll a higher number than B, B is more likely to roll a higher number than C, and C is more likely to roll a higher number than A. This shows that the relation "is more likely to roll a higher number" is not transitive with these dice, and so we say this is a set of nontransitive dice.
How easy is it to model such nontransitive dice in a high level constraint programming system such as MiniZinc or Comet? The basic MiniZinc model nontransitive_dice.mzn - with no bells and whistles - is this (somewhat edited):
int: m = 3; % number of dice
int: n = 6; % number of sides of each die

% the dice
array[1..m, 1..n] of var 0..n*2: dice :: is_output;

% the competitions: 
% The last wrap around is the one that breaks 
% the transitivity.
array[0..m-1, 1..2] of var 0..n*n: comp :: is_output;

solve satisfy;
constraint
   % order the number of each die
   forall(d in 1..m) (
       increasing([dice[d,i] | i in 1..n])
   )
;
 
% and now we roll...
% Number of wins for [A vs B, B vs A]
constraint
   forall(d in 0..m-1) (
      comp[d,1] = sum(r1, r2 in 1..n) (
          bool2int(dice[1+(d mod m), r1] >      
                   dice[1+((d + 1) mod m), r2]))
      /\
      comp[d,2] = sum(r1, r2 in 1..n) (
          bool2int(dice[1+((d+1) mod m), r1] >
                   dice[1+((d) mod m), r2]))
   )
;

% non-transitivity
% All dice 1..m-1 must beat the follower, 
% and die m must beat die 1
constraint 
  forall(d in 0..m-1) (
    comp[d,1] > comp[d,2]
   )
;
The result of the m competitions is placed in the (m x 2) matrix comp where the winner is in comp[i,1] and the loser in comp[i,2] for the match i. The last constraint section in the code above ensures that the winner is always in comp[i,1].

In the full model I have added the following:
  • max_val: maximum value of the dice, to be minimized
  • max_win: maximum number of winnings, to be maximized
  • gap and gap_sum: the difference of wins of a specific competition, to be maximized or minimized
  • the example setups from the Wikipedia page
The Comet version nontransitive_dice.co, includes about the same functions as the MiniZinc model. Here is the basic model (also edited for presentational purposes):
import cotfd;
int m = 3; // number of dice
int n = 6; // number of sides of each die
Solver cp();

// the dice
var{int} dice[1..m, 1..n](cp, 1..n*2);

// The competitions: 
var{int} comp[0..m-1, 1..2](cp, 0..n*n);

explore  {
  // symmetry breaking: order the number of each die
  forall(d in 1..m) {
    forall(i in 2..n) {
      cp.post(dice[d,i-1] <= dice[d,i]);
    }
  }

  // and now we roll...
  // Number of wins for [A vs B, B vs A]
  forall(d in 0..m-1) {
    cp.post(comp[d,1] == 
        sum(r1 in 1..n, r2 in 1..n) (
            dice[1+(d % m), r1] >
            dice[1+((d + 1) % m), r2]));

    cp.post(comp[d,2] == 
         sum(r1 in 1..n, r2 in 1..n) (
            dice[1+((d+1) % m), r1] >
            dice[1+((d) % m), r2]));
  }

  forall(d in 0..m-1) {
    cp.post(comp[d,1] > comp[d,2]);
  }


} using {

  labelFF(dice);
  labelFF(comp);      

  cout << "dice:" << endl;
  forall(i in 1..m) {
    forall(j in 1..n) {
      cout << dice[i,j] << " ";
    }
    cout << endl;
  }
  cout << endl;
}
The full Comet model when minimzing max_value for 3 dice with 6 sides of each die (with domain 1..12), gives this result, i.e. the maximum value used in the dice is 4:
solution #1
dice:
1 1 1 2 4 4
1 1 1 3 3 3 
1 2 2 2 2 2 
comp:
1 vs 2: 15 (p:0.416667) 12 (p:0.333333) 
2 vs 3: 18 (p:0.500000) 15 (p:0.416667) 
3 vs 1: 15 (p:0.416667) 13 (p:0.361111) 
gap[3,3,2]
max_val: 4
max_win: 18
gap_sum: 8
Somewhat related is another dice problem: sicherman_dice.mzn.

17x17 challenge (not solved, though)

Last November, William Gasarch published a challenge in The 17x17 challenge. Worth $289.00. This is not a joke.. Also, see bit-player's The 17×17 challenge, and O.R. by the Beach Let’s Join O.R. Forces to Crack the 17×17 Challenge.

The objective is to create a 17x17 matrix where there are no colored sub-rectangles, i.e. no (sub)rectangles in the matrix can have the same color (number). Here's a solution of the 5x5x3 problem, i.e. 5 x 5 matrix and with 3 colors (1..3):
1 3 3 2 1 
3 1 2 2 1 
3 2 1 2 1 
2 2 2 1 1 
1 1 1 1 2 
When I first saw this problem in December, I tried some models but didn't pursued it further. However, when Karsten Konrad published some OPL models, I picked it up again. We also discussed this problem to see what a join force could do.

Here is my take in MiniZinc: 17_b.mzn. It is inspired (or rather a translation) of Karsten Konrad's OPL model from his blog post Let’s do real CP: forbiddenAssignment. Karsten earlier posted two OPL models: However, 17x17x4 is very hard, and - as I understand it - no one has yet cracked it.

Here are three MiniZinc models with some different approaches. A Comet version is mentioned below.

17_b.mzn
This is a translation of Karsten Konrad's OPL model in Let’s do real CP: forbiddenAssignment, which uses a forbiddenAssignments constraint. However, MiniZinc don't has this constraint as a built-in so I had to roll my own, here called existential_conflict. The forbidden is an table of forbidden combinations in the matrix. For 4 colors it look like this:
1 1 1 1
2 2 2 2
3 3 3 3
4 4 4 4
i.e. there can be no rectangle with has all the same values. (If you now start to think about global cardinality constraint, please read on.) Here is the complete model:
int: NbColumns = 10;
int: NbRows = 10;
int: NbColors = 4;
set of int: Columns = 1..NbColumns;
set of int: Rows = 1..NbRows;
set of int: Colors = 1..NbColors;

array[Rows,Columns] of var Colors: space :: is_output;
array[Colors, 1..4] of int: forbidden = 
   array2d(Colors, 1..4,[ i | i in Colors, j in 1..4]);

% The pattern x must not be in table table.
predicate extensional_conflict(array[int, int] of var int: table, 
                               array[int] of var int: x) =
   not exists(pattern in index_set_1of2(table)) (
      forall(j in index_set_2of2(table)) (
           x[j] = table[pattern, j]
      )
   )
;

solve :: int_search(
   [space[i,j] | i in Rows, j in Columns], 
   first_fail, 
   indomain_min, 
  complete) 
satisfy;

constraint
  space[1,1] = 1 /\
  space[NbRows,NbColumns] = 2 /\
   
  forall(r in Rows, r2 in 1..r-1, c in Columns, c2 in 1..c-1) (
    extensional_conflict(forbidden, 
       [space[r,c], space[r2,c], space[r,c2], space[r2,c2]])
  )
17_b2.mzn
This is almost exactly the same as the one above, except that it has another implementation of existential_conflict, a kind of MIP version:
predicate extensional_conflict(array[int, int] of var int: table, 
                               array[int] of var int: x) =
 forall(pattern in index_set_1of2(table)) (
     sum(j in index_set_2of2(table)) ( 
           bool2int(x[j] = table[pattern, j])
     ) < 4
  )
;
17_b3.mzn
This third version use the constraint global cardinality constraint instead of forbidden assignments/existential_conflict. We create a temporary array (gcc), and then check that there is no value of 4 (i.e. the number of occurrences of each number must be less than 4). Also, I skipped the symmetry breaking of assigning the corner values in the matrix space (commented below).
constraint
  forall(r in Rows, r2 in 1..r-1, c in Columns, c2 in 1..c-1) (
    let {
       array[1..4] of var 0..4: gcc
    } in
    global_cardinality(
        [space[r,c], space[r2,c], space[r,c2], space[r2,c2]], 
        gcc) 

    /\ 
    forall(i in 1..4) (  gcc[i] < 4  )
  )
;
As you can see these are rather simple (a.k.a. naive) models. I suspect that there are better search strategies than the one I tested, for example Gecode's size_afc_max with indomain_random which is one of the better:
solve :: int_search(
         [space[i,j] | i in Rows, j in Columns], 
         size_afc_max,
         indomain_random,
         fail) 
     satisfy;
However, neither of these models could find any solution of the original 17x17x4 problem. The longest take was 36 hours on my new 8 core Linux (Ubuntu) and 12 Gb RAM, but no luck.

Here is a solution of the 17x17x5 (i.e. 5 colors), using the LazyFD solver (it took 1:19 minutes), c.f. Karsten's Improved Model for Coloring Problem.
1 5 1 5 1 5 3 2 5 5 4 4 2 3 1 3 2 
5 1 4 3 5 5 2 3 2 2 3 2 4 1 4 4 4 
3 2 2 2 1 3 2 4 1 3 5 3 4 4 2 5 5 
2 2 3 3 5 4 4 2 5 4 4 1 3 1 1 1 5 
4 4 5 3 1 4 1 1 3 5 2 2 5 3 2 5 3 
5 1 3 1 1 3 3 2 2 1 4 1 5 4 2 4 5 
4 2 4 3 5 3 5 1 2 4 2 5 2 4 1 2 1 
2 4 2 4 1 5 4 3 4 2 5 3 4 1 5 2 3 
4 3 5 2 1 3 4 3 5 2 5 5 3 3 4 1 2 
1 5 3 1 2 1 4 1 3 2 3 5 2 4 2 5 4 
2 1 3 5 2 1 5 3 2 3 5 4 5 2 4 1 1 
5 3 5 4 3 4 3 2 1 3 3 2 1 2 1 2 4 
4 3 3 4 4 1 5 5 1 5 4 1 2 2 5 3 3 
5 4 4 1 3 2 1 4 2 5 4 5 3 5 3 1 3 
2 4 1 3 4 2 1 5 1 3 1 4 5 5 2 4 1 
3 3 5 2 4 2 5 4 3 1 1 2 1 5 3 1 4 
1 5 2 3 4 2 1 5 4 4 5 1 1 2 3 3 2 
I also wrote a Comet model, 17_b.co, and experimented with different cardinality constraints (cardinality, atmost, a decomposition of forbiddenAssignments) and different labelings.

Ormat game

Ormat game is yet another grid problem. From bit-player: The Ormat Game (where the nice pictures has been replace with 0/1 matrices below):
Here’s the deal. I’m going to give you a square grid, with some of the cells colored and others possibly left blank. We’ll call this a template. Perhaps the grid will be one of these 3×3 templates:
 
   1 0 0      1 1 1     1 1 1
0 1 1 1 1 1 1 1 1
0 1 1 1 1 1 1 1 0
(1) (2) (3)
You have a supply of transparent plastic overlays that match the grid in size and shape and that also bear patterns of black dots:
   1 0 0     1 0 0     0 1 0
0 1 0 0 0 1 1 0 0
0 0 1 0 1 0 0 0 1
(a) (b) (c)
0 0 1 0 1 0 0 0 1
1 0 0 0 0 1 0 1 0
0 1 0 1 0 0 1 0 0
(d) (e) (f)
Note that each of these patterns has exactly three dots, with one dot in each row and each column. The six overlays shown are the only 3×3 grids that have this property.

Your task is to assemble a subset of the overlays and lay them on the template in such a way that dots cover all the colored squares but none of the blank squares. You are welcome to superimpose multiple dots on any colored square, but overall you want to use as few overlays as possible. To make things interesting, I’ll suggest a wager. I’ll pay you $3 for a correct covering of a 3×3 template, but you have to pay me $1 for each overlay you use. Is this a good bet?
This is the first version in MiniZinc: ormat_game.mzn where the problem instances are included in the model file. As with the 17x17 problem above, this is - in concept, at least - quite simple. The model, sans the overlays, is shown below:
include "globals.mzn"; 
% the number of overlays, n!
int: f :: is_output = product([i | i in 1..n]);

% Which overlay to use.
array[1..f] of var 0..1: x :: is_output;

% number of overlays used (to minimize)
var 1..f: num_overlays :: is_output = sum(x);

% There are n! possible overlays
% They are in a separate .dzn file
%    ormat_game_n3.dzn etc
array[1..f, 1..n,1..n] of 0..1: overlays;

solve :: int_search(
   x, 
   smallest, 
   indomain_min, 
   complete) 
  minimize num_overlays;

constraint 
    % if problem has a black cell (=1) then there 
    % must be a selected overlay that has a black cell
    forall(i,j in 1..n) (
        (
        problem[i,j] = 1 -> 
            exists(o in 1..f) (
                               x[o]*overlays[o,i,j] = 1 
                               )
        )                                            
    )
    /\ % the inverse: wherever a selected overlay 
       %  has a black cell, problem must have 
       %  a black cell
    forall(o in 1..f) (
        x[o] = 1 -> 
           forall(i,j in 1..n) (
              overlays[o,i,j] = 1 -> problem[i,j] = 1
           ) 
    ) 
;

% Problem grid 3
include "ormat_game_n3.dzn";
array[1..n, 1..n] of 0..1: problem = array2d(1..n, 1..n,
[
   1,1,1,
   1,1,1,
   1,1,0
]);
The only constraints are:
  • if the problem matrix has a black cell (1) then there at least one of the selected overlay matrices must have a 1 in this cell
  • the converse: for all the black cells in the selected overlays the problem must has a black cell.
The output and solution:
Problem:
111
111
110

Solution:
x: [0, 1, 0, 1, 1, 1]
num_overlays: 4

Overlay #2
100
001
010

Overlay #4
001
100
010

Overlay #5
010
001
100

Overlay #6
001
010
100
Ideally, the overlays should be calculated in real-time, but it was much easier to generate overlay files for each size (using ormat_game_generate.mzn and some post-processing in Perl). These must be included in the model as well. Note: there are n! overlay matrices for problems of size n x n. However, changing the model file each time a new problem is made or tested is not very convenient. It's better to have problem instances in separate files which include the model and the overlay files of appropriate size. These following problem instances includes the general model ormat_game_nodel.mzn and an overlay file: Note: As of writing, there is problems running problems of size 7 on my 64-bit Ubuntu; it works on the 32-bit machine (also Linux), though. This has been reported to the G12 team.

Other new models

Here are some other new models (or not mentioned earlier):

Contiguity using regular

contiguity_regular.mzn: decomposition of the global constraint contiguity using regular constraint.

Combinatorial auction

These two MiniZinc models implements a simple combinatorial auction problem (Wikipedia): combinatorial_auction.mzn, and a "talkative" MIP version combinatorial_auction_lp.mzn (which is translated from Numberjack's Tutorial)

Compare with the Gecode version combinatorial_auction.cpp.

August 27, 2010

MiniZinc version 1.1.6 released

MiniZinc release 1.1.6 has been released. It can be downloaded here.

From the NEWS

G12 MiniZinc Distribution 1.1.6
-------------------------------

Changes in this release:

* We have modified the decomposition of the global constraint lex_lesseq
in order to avoid the introduction of an auxiliary Boolean variable.
(Thanks to Chris Mears and Todd Niven for pointing this out.)


Bugs fixed in this release:

* mzn2fzn now correctly computes the set of output variables when the
output item contains let expressions. [Bug #141]

* A bug that caused mzn2fzn to infer incorrect bounds for integer and
float var array elements has been fixed. [Bug #149]

* mzn2fzn now prints the source locations of all solve (output) items when
the are multiple such items. [Bug #143]

* mzn2fzn now flattens par expressions containing the built-in operation
pow/2 correctly.

* mzn2fzn now flattens arrayNd expressions containing arrays of strings
correctly.

* The mzn script no longer aborts if the model contains an array of
decision variables. [Bug #140]

July 22, 2010

Minizinc version 1.1.5 released

Version 1.1.5 of MiniZinc has been released. It can be downloaded here.

From the NEWS:

Bugs fixed in this release:

* We have fixed a number of problems that caused stack overflows in mzn2fzn.

* The FlatZinc interpreter's MIP backend no longer reports "unknown" for
satisfaction problems.

July 20, 2010

Minizinc version 1.1.4 released

MiniZinc version 1.1.4 has been released. It can be downloaded here.

From the NEWS:

Changes in this release:

* We have added a library of global constraint definitions and FlatZinc built-in redefinitions suitable for use with LP/MIP solvers; both are in the "linear" directory of the MiniZinc library.

Bugs fixed in this release:

* Some performance issues with mzn2fzn that occurred when flattening models that generate and traverse large arrays have been fixed.

* An omission that caused mzn2fzn not to recognise the MiniZinc built-in function round/1 has been corrected.

* A bug in flatzinc that caused the MIP backend to abort when the model instance contained an unused set parameter has been fixed. [Bug #134]

* A bug in mzn2fzn that caused it not to place domain constraints on the FlatZinc variables generated for an array of variables introduced via a let expression has been fixed. [Bug #133]

* The implementation of the div propagator in flatzinc's FD backend has been modified to avoid potentially long fixpoint computations.

July 18, 2010

Some new MiniZinc models

Here are some new MiniZinc models. Some are completely new - or previously unannounced - but there are also some which has been implemented in some other constraint programming system before.

New models

Here are some new - or at least previously not announced - models:

Latin square card puzzle

latin_square_card_puzzle.mzn: Latin square card puzzle

I found this problem in Mario Livio's nice PopSci book about the development of group theory The Equation that couldn't be solved, page 22:
... Incidentally, you may get a kick out of solving this eighteenth century card puzzle: Arrange all the jacks, queens, kings, and aces from a deck of cards in a square so that no suit or value would appear twice in any row, column, or the two main diagonals.
My general approach is to use integers of the following form. n is the size of the matrix (here 4) and m is the number we use for modulo operation (here 10). These values are calculated automatically by the model depending on n.
 % values: i mod 10 
  0, 1, 2, 3,  % suite 1: 0 div 10
 10,11,12,13,  % suite 2: 1 div 10
 20,21,22,23,  % suite 3: 2 div 10
 30,31,32,33   % suite 4: 3 div 10
What then want that the values divided by 10 (the suites) should be different in each row, column, and diagonals, and also the values by modulo 10 (the values) is different in each row, column, and diagonals.

With the symmetry breaking constraint that the value 0 must be in the upper leftmost cell, there are 72 solutions for n = 4. Here is one of them:
 0 33 22 11
21 12  3 30
13 20 31  2
32  1 10 23
Note: There are solutions for n = 4 and n = 5 but not for n = 6. The n = 6 problem is the same as Euler's 36 officer's problem, which thus is not solvable. Also see MathWorld.

Investment problem

This problem is from ORMM (Operations Research Models and Methods):
A portfolio manager with a fixed budget of $100 million is considering the eight investment opportunities shown in Table 1. The manager must choose an investment level for each alternative ranging from $0 to $40 million. Although an acceptable investment may assume any value within the range, we discretize the permissible allocations to intervals of $10 million to facilitate the modeling. This restriction is important to what follows. For convenience we define a unit of investment to be $10 million. In these terms, the budget is 10 and the amounts to invest are the integers in the range from 0 to 4.
Here are two implementations:

Arg max

argmax.mzn

argmax/argmin predicate
  • argmax_ge(pos, x)
    pos is the index x for the maximum value(s). There can be many maximum values.
  • argmax_gt(pos, x)
    pos is the index x for the maximum value. There can be only one maximum value.
  • argmin_le(pos, x)
    pos is the index x for the minimum value(s). There can be many minimum values.
  • argmin_lt(pos, x)
    pos is the index x for the minimum value. There can be only one maximum value.

Permutation number

permutation_number.mzn

A permutation number is the number of transpositions in a permutation, see Wikipedia Permutation.

Sicherman dice

sicherman_dice.mzn

From Wikipedia Sicherman_dice:
Sicherman dice are the only pair of 6-sided dice which are not normal dice, bear only positive integers, and have the same probability distribution for the sum as normal dice.

The faces on the dice are numbered 1, 2, 2, 3, 3, 4 and 1, 3, 4, 5, 6, 8.
I read about this problem in a book/column by Martin Gardner long time ago, and got inspired to model it now by the recent WolframBlog post Sicherman Dice

Here is the vital part of the code:
array[2..12] of int: standard_dist = 
       array1d(2..12, [1,2,3,4,5,6,5,4,3,2,1]);
% the two dice
array[1..n] of var 1..m: x1 :: is_output;
array[1..n] of var 1..m: x2 :: is_output;
constraint
forall(k in 2..12) (
  standard_dist[k] = 
    sum(i,j in 1..n) ( 
       bool2int(x1[i]+x2[j] == k)
    )
)
% symmetry breaking
/\ increasing(x1) 
/\ increasing(x2)

/\ % x1 is less or equal to x2
forall(i in 1..n) (
   x1[i] <= x2[i]
)
% /\ lex_lesseq(x1, x2)
This model gets the two different solutions, first the standard dice and then the Sicherman dice:
[1, 2, 3, 4, 5, 6]
[1, 2, 3, 4, 5, 6]

--

[1, 2, 2, 3, 3, 4]
[1, 3, 4, 5, 6, 8]
Extra: If we also allow that 0 (zero) is a valid value then the following two solutions are also shown. The only thing we have to do is to change the domains of x1 and x2:
% the two dice
array[1..n] of var 0..m: x1 :: is_output;
array[1..n] of var 0..m: x2 :: is_output;
Here here are the two new solutions:
[0, 1, 1, 2, 2, 3]);
[2, 4, 5, 6, 7, 9]);

--

[0, 1, 2, 3, 4, 5]);
[2, 3, 4, 5, 6, 7]);
These two extra cases are mentioned in MathWorld SichermanDice.

Translations of other models

The following MiniZinc models is ports from models written in at least one other constraint programming system before, mostly Comet:

June 11, 2010

MiniZinc version 1.1.3 released

MiniZinc version 1.1.3 has been released. It was be downloaded here.

From the NEWS file:

Changes in this release:

* We have added a new script, mzn, that allows output items to work with two-pass MiniZinc evaluation. (The script requires a Unix-like system -- we hope to lift restriction in later versions.)
* The files alldifferent.mzn, atmost1.mzn, atmost.mzn and atleast.mzn have been added to the MiniZinc globals library. At the moment these files merely cause all_different.mzn, at_most1.mzn etc to be included. Eventually the latter will be replaced by the former.

Bugs fixed in this release:

* A bug in mzn2fzn that caused an internal error when flattening predicates with a reified form has been fixed. [Bug #131]
* The MiniZinc type checker now correctly reports an error for attempts to use the built-in function index_set/1 with arrays that have more than one dimension. [Bug #68]
* The broken definition for refied all_different for the lazy clause generation solver has been fixed.
* A bug where mzn2fzn was mishandling arrays of strings has been fixed.

May 23, 2010

Two new tools for MiniZinc, and a paper

Some days ago, the G12 group released two new tools for MiniZinc that I haven't used that much, but will hopefully do in the not to far future. Also, a recent paper is mentioned.

G12 IDE

The G12 IDE is an application for writing, running, visualizing, and debugging MiniZinc models. It is based on (the Java IDE) Eclipse. It can be downloaded here.

fzn2xcsp

fzn2xcsp is a tool for converting a subset of FlatZinc to XCSP 2.1. As of writing this tools is only available in the development version.

Paper: Philosophy of the MiniZinc challenge

Peter J. Stuckey, Ralph Becket, and Julien Fischer: Philosophy of the MiniZinc challenge (Springer Link). It is published in the Constraints Journal, but the paper is not available there (yet).
Abstract
MiniZinc arose as a response to the extended discussion at CP2006 of the need for a standard modelling language for CP. This is a challenging problem, and we believe MiniZinc makes a good attempt to handle the most obvious obstacle: there are hundreds of potential global constraints, most handled by few or no systems. A standard input language for solvers gives us the capability to compare different solvers. Hence, every year since 2008 we have run the MiniZinc Challenge comparing different solvers that support MiniZinc. In this report we discuss the philosophy behind the challenge, why we do it, how we do it, and why we do it that way.

Beside being an interesting paper about the MiniZinc challenge (see the MiniZinc homepage for links to the last two year's challenges, and this year's), it is also the first constraint programming paper where I'm mentioned (in the Acknowledgment). Thanks for this, I'm honored and appreciate it very much.

May 14, 2010

Optimizing shopping baskets: The development of a MiniZinc model

Yesterday (Thursday) was a holiday in Sweden and I had a busy day. Apart from writing a modulo propagator for ECLiPSe/ic, I also spend a lot of time with the Stack Overflow question How to use constraint programming for optimizing shopping baskets?. It asked for some tips (in a Java constraint programming system) for solving a problem of optimizing a shopping basket where different shops has different prices for items, as well as a delivery cost for orders below a certain limit.

You can read what I wrote and the developments of the models in my answer (I am hakank).

Update 2010-05-16: Added a better version. See below, Version 6.

First version

Even though the question was about a Java constraint system (I mentioned JaCoP, though), I rather quickly throw together a MiniZinc model instead as some kind of a prototype for a Java model. I'm a very lazy person and MiniZinc is excellent for prototyping.

This first version is shopping_basket.mzn (which actually include the next requirement: to handle shops which don't have all items). It solved the very simple example described in the problem section without any problem.

This model is quite straightforward: x is a matrix of 0/1 decision variables which shop to buy which item (the "knapsack" part). The more tricky part is to calculate the delivery costs if the total is below a certain limit (delivery_costs_limits). Here is the constraint part of the model:
constraint
   % we must buy all items (from some shop)
   forall(i in 1..num_items) (
       sum([x[i,j] | j in 1..num_shops]) = 1
   )
   /\ 
   total = sum(i in 1..num_items) (
       % costs of the products
       sum(j in 1..num_shops) (
          x[i,j]*costs[i,j]
       )
   ) + sum(j in 1..num_shops) (
       %  and delivery costs if total for a shop < limit
       let {
         var int: this_cost = sum([costs[i,j]*x[i,j] | i in 1..num_items])
       } in
       delivery_costs[j]*bool2int(this_cost > 0 /\ this_cost < delivery_costs_limits[j])
   );

The N/A problem: Not all shops sells all items

The next question was how to handle the "N/A" cases, i.e. for shops that don't have some item. It was solved by setting the price to a large number (999999 or 99999). However this is not very nice to the constraint solver. I really wanted to set these N/A to 0 (zero) but didn't get this correct (see below for a wrong turn on this).

As mentioned above the model shopping_basket.mzn includes the N/A problem.

Another representation

Almost always there are different approached how to model a problem. In the next version, shopping_basket2.mzn, I represented x - instead of a 0/1 matrix - as an array of length 1..num_items with the domains 1..num_shops to represent which shop to buy an item. The idea was that the constraint solvers now had much less work to do.

The constraint section is now simpler in some part, but calculating the delivery costs required a little more code:
constraint
 total = sum(i in 1..num_items) (
     costs[i,x[i]]
 ) + 
 sum(j in 1..num_shops) (
   let {
   var int: this_cost = 
     sum([costs[i,j]*bool2int(x[i]==j) | i in 1..num_items])
   } in
   delivery_costs[j]*bool2int(
      this_cost > 0 /\ 
      this_cost < delivery_costs_limits[j])
 );

Real data

The simple example data was easily solved for these two models, but the real challenge began when handling real data: 29 items and 37 shops. The two simplistic models was, however, not strong enough to handle this, so I started to think harder about the problem.

A wrong turn

The first attempt was continuing with the matrix approach and trying to be clever to represent N/A by 0 (zero), see above. This was done in shopping_basket3.mzn, but this model is plain wrong! The reason for this rather embarrissing bug was that the test case I used was not correct (completely mea culpa of course).

Final(?) version, part I: limiting the domains

OK, back to the drawing board. After some deep thoughts (i.e. an afternoon nap) I realized that the first models was too much of a proof-of-concept, too much prototype. In these I had broken one of the cardinality rules of constraint programming: not thinking about the domains of all decision variables.

In the first models the total cost (total) was defined as
var int: total :: is_output;
and the temporary variables this_cost was also defined as var int. This means that there is no limits of these variables (not even that they should be positive) which gives little hints for the constraint solvers.

The remedy is shown in the final(?) version shopping_basket5.mzn (shopping_basket4.mzn is basically the same thing for the matrix approach, and was not published).

Here the maximum total for any shop is first calculated, which can be quite large given the 99999 for N/A, but it is still limited and positive.
% total price
int: max_total :: is_output = 
      sum(i in 1..num_items) (
         max([costs[i,j] | j in 1..num_shops] )
      );
total is then defined with this limit:
var 0..max_total: total :: is_output; 
The temporary variable this_cost is handled accordingly.

Final(?) version, part II: labeling

This model was now in better shape, but still too slow. For all these versions I tried a lot of different labelings (search heuristics) and with many solvers (including the two MIP solvers, MiniZinc/mip and ECLiPSe/eplex, but they didn't accept the models).

I first saw a real improvement with the combination of Gecode/fz and the largest,indomain_max labeling. It solved the problem in 25 seconds (and with 69964 failures) :
total = 42013
x = [13, 20, 17, 18, 18, 13, 17, 17, 20, 13, 17, 17, 13, 13, 
     18, 17, 13, 13, 17, 8, 13, 36, 10, 17, 13, 13, 17, 20, 13] 
Just a minute or so after publishing this result (which seemed to please the person asking the question), I tested a different labeling, largest,indomain_split, and it solved the problem in 12 seconds (51013 failures). It gave a slightly different solution of where to buy item 12 (shown as bold):
total = 42013
x = [13, 20, 17, 18, 18, 13, 17, 17, 20, 13, 17, 13, 13, 13, 
     18, 17, 13, 13, 17, 8, 13, 36, 10, 17, 13, 13, 17, 20, 13]. 
These two solutions are the only two solutions for the (optimal) total of 42013. The test for all solutions takes about 1 second.

Some comments

This was a fun problem, quite simple but instructive what will happen when ignoring declaration of proper domain variables in the first versions of the models.

Even though the "real problem" was of the magnitude of 29 items and 37 shops (29 * 37 =1073), I am somewhat surprised that the problem was so hard to solve. This reminds me of the coins_grid problem which MIP solvers solve very fast, but the constraint solver have problems with it.

Update 2010-05-16: Version 6: Now faster

Well, version 5 was not the final version. Maybe this version 6 (shopping_basket6.mzn) is?

By two simple changes the problem is now solved in about 2 seconds and with 4460 failures (compared with 12 seconds and 51013). Here are the changes:
* The first change was the representation of N/A:s from 99999 to 0 (zero) . As noted above the former representation is not a good since it make all calculations, and the domains, unnecessary large.
* And finally, the following constraint states that all prices to consider must be larger than 0:
   forall(i in 1..num_items) (
      costs[i,x[i]] > 0
   ) 
Sigh. Of course! No defense of my silliness there is.

This also lessens the surprise considerable in Some comment above.

May 11, 2010

MiniZinc version 1.1.2 released

Version 1.1.2 of MiniZinc has been released. Download.

Changes from version 1.1.1 (from NEWS):

G12 MiniZinc Distribution 1.1.2
--------------------------------
Bugs fixed in this release:

* The file diffn.mzn is now included in globals.mzn.

* A bug in mzn2fzn that caused it to abort when flattening int2float expressions has been fixed.

* An error in the FlatZinc specification has been fixed. All var types may have assignments, not just arrays.

* A bug in mzn2fzn that caused it generate an array_int_element constraint where an array_var_int_element constraint was required has been fixed. [Bug #122]

* A bug that caused mzn2fzn to generate invalid FlatZinc rather than emit an error message when the bound of an unbounded variable is taken has been fixed.


May 06, 2010

Some 40 new ECLiPSe models

During the last days I have published about 40 new ECLiPSe models on my ECLiPSe page. They are all listed below.

A major part are ports from my SICStus Prolog models that I wrote some months ago, and I have not changed much from the SICStus models. Also, most of them - in turn - was first written in MiniZinc (see my MiniZinc page).

However, some small puzzles are completely new. By a freaky coincidence two are some kind of weighting problems (for comparison, the corresponding MiniZinc model are also linked):

The new models

Here is all the new ECLiPSe models.

Related

I think I have forgot to mention Helmut Simonis' great ECLiPSE ELearning Website, which includes 20 interesting chapters, most with video lectures. Even if the code examples are in ECLiPSe, the material is definitely of general interest for anyone who want to learn more about Constraint Programming.

March 29, 2010

MiniZinc Challenge 2010

MiniZinc Challenge 2010 has been announced:

The Challenge

The aim of the challenge is to start to compare various constraint solving technology on the same problems sets. The focus is on finite domain propagation solvers. An auxiliary aim is to build up a library of interesting problem models, which can be used to compare solvers and solving technologies.

Entrants to the challenge provide a FlatZinc solver and global constraint definitions specialized for their solver. Each solver is run on 100 MiniZinc model instances. We run the translator mzn2fzn on the MiniZinc model and instance using the provided global constraint definitions to create a FlatZinc file. The FlatZinc file is input to the provided solver. Points are awarded for solving problems, speed of solution, and goodness of solutions (for optimization problems).

....

For the rules, see: MiniZinc Challenge 2010 -- Rules.

Also, see the earlier challenges:

March 26, 2010

MiniZinc version 1.1.1 released

MiniZinc version 1.1.1 has been released. Download.

From the NEWSG12 MiniZinc Distribution version 1.1.1
---------------------------------------

Bugs fixed in this release:

* A bug that caused predicate arguments to be incorrectly flattened in
reifying contexts has been fixed. [Bug #109]

* A bug in mzn2fzn that caused incorrect bounds to be calculated for the
result of a mod operation has been fixed. [Bug #107]

* A bug in mzn2fzn that caused out of range array accesses to be generated in
reified contexts, instead of just making the result of the reification
false. [Bug #110]

* The omission of the null annotation from the Zinc / MiniZinc specification
has been fixed.

* The rostering problem in the MiniZinc benchmark suite (benchmarks/roster),
has been reformulated. The old formulation was always unsatisfiable under
the change to the semantics of the mod operation introduced in MiniZinc 1.1.
[Bug #108]

* A bug in mzn2fzn that caused it to emit the null/0 annotation in the
generated FlatZinc. [Bug #111]

March 17, 2010

MiniZinc version 1.1 released

MiniZinc version 1.1 has been released. See below how my existing (and maybe others') MiniZinc models are affected by the changes.

From the NEWS:

G12 MiniZinc Distribution version 1.1
-------------------------------------

Changes to the MiniZinc language:

* The following built-in operations have been introduced:

int: lb_array(array[int] of var int)
float: lb_array(array[int] of var float)
set of int: lb_array(array[int] of var set of int)

int: ub_array(array[int] of var int)
float: ub_array(array[int] of var float)
set of int: ub_array(array[int] of var set of int)

set of int: dom_array(array[int] of var int)

These new operations are synonyms for the following existing built-in
MiniZinc operations:

int: lb(array[$T] of var int)
float: lb(array[$T] of var float)
set of int: lb(array[$T] of var set of int)

int: ub(array[$T] of var int)
float: ub(array[$T] of var float)
set of int: ub(array[$T] of var set of int)

set of int: dom(array[$T] of var int)

These latter operations are now deprecated. Support for them will
be removed in the next release. This change is being made in order
to preserve compatibility with the full Zinc language.

Note: that only the versions of lb, ub and dom that take an array
as a an argument are deprecated. The MiniZinc lb, ub and dom operations
on non-array values are *not* deprecated.


Changes to the FlatZinc language:

* Boolean variable expressions as constraints are no longer supported.
All constraints in FlatZinc must now be predicate applications.

* String parameters are no longer supported. String literals are restricted
to appearing as the arguments of annotations.

* Set of bool and set of float parameters and literals are no longer
supported.

* The int_float_lin/4 objective expression is no longer supported.

* FlatZinc now has two additional evaluation outcomes: "unknown"
for when search terminates without having explored the whole search
space and "unbounded", for when the objective of an optimization
problem is unbounded.

* The semantics of the int_div/3 and int_mod/3 built-ins has been changed.
See the ``Specification of FlatZinc'' for further details.


Other Changes:

* The single pass MiniZinc interpreter, minizinc, has been deprecated.
It will be removed in a future release.

* The MiniZinc-to-FlatZinc converter, mzn2fzn, has been rewritten.
The new implementation is smaller and more efficient.
Computation of variable bounds has also been improved.

* mzn2fzn now outputs singleton sets as ranges. [Bug #94]

* A bug that caused expressions containing abs/1 to be incorrectly
flattened has been fixed. [Bug #91]

* The FlatZinc interpreter's finite-domain backend now implements
global_cardinality_low_up as a built-in.

* The FlatZinc interpreter's lazy clause generation solver now supports
the int_mod/3 built-in.

* Two additional modes of operation have been added to the FlatZinc
solution processing tools, solns2dzn, that allow it to extract the first
or last solution from a FlatZinc output stream. Also, there is no longer
a default mode of operation for solns2dzn, it must now be specified by
the user or an error will occur.

* The following new global constraints have been added to the MiniZinc
library:

all_equal
decreasing
diffn
lex2
lex_greater (Synonym for lex_less with arguments swapped.)
lex_greatereq (Synonym for lex_lesseq with arguments swapped.)
sliding_sum
strict_lex2

* The following synonyms for existing global constraints have been added
to the MiniZinc library (the existing name is given in parentheses):

alldifferent (all_different)
atleast (at_least)
atmost (at_most)
atmost1 (at_most1)

* The sequence constraint is deprecated. Models should use the new
sliding_sum constraint instead.

* The 'table' constraint decompositions in the MiniZinc library have been
modified so as to fit better with the G12 MiniZinc-to-FlatZinc conversion:
now no scaling constraints are created.

* The decompositions of the constraints in the 'lex' family have been
tweaked to enable a little more propagation.

Documents

Here are the three new specification documents for MiniZinc version 1.1 (and Zinc version 0.11):

Comments

The next few days I will change my MiniZinc models so they comply to version 1.1, and I have already started this work. Update 2010-03-20: These changes has now been done, including updating the SVN repository.

The following will be changed:

  • lb/ub for arrays

    lb(array) and ub(array) will be changed to lb_array(array) and ub_array(array) respectively.
  • Comparing/copying arrays

    One thing that should not work even in MiniZinc version 1.0.3 - but for some reason did - was copying/equality/comparison of arrays in the constraint section or in predicates. This don't work in MiniZinc version 1.1. E.g., the following no longer works:

    int: n = 4;
    array[1..n] of var 1..n: x;
    constraint
    x = [1,2,3,4] % no longer works
    ;

    Instead, the arrays must now be handled element-wise. Since many of my models use the above construct, especially for testing the global constraints, the models use a new predicate family cp<n>d (where <n> is the dimension, 1, 2, etc), e.g. cp1d and cp2d. Example of one version of cp1d:

    int: n = 4;
    array[1..n] of var 1..n: x;
    solve satisfy;
    % arrays of 1d where both arguments are var int
    predicate cp1d(array[int] of var int: x, array[int] of var int: y) =
    assert(index_set(x) = index_set(y),
    "cp1d: x and y have different sizes",
    forall(i in index_set(x)) ( x[i] = y[i] ));
    constraint
    cp1d(x, [1,2,3,4]) % this works
    ;

    Some examples are collected in the model copy_arrays.mzn.

    I estimate that over 200 of my models have to be fixed in this way.As mentioned above, some of models are now already changed.

  • Renamed models

    Some of my MiniZinc models has been renamed since they clash with new built-in predicates:

After the changes are done, I will also update the G12's MiniZinc SVN repository, the hakank directory.

Two more things

Also see, my MiniZinc Page.

March 12, 2010

Pi Day Sudoku 2009 - the models (MiniZinc and Comet)

The 14 March is Pi Day (π Day) and is a celebration of Pi (3.14 in the mm.dd date notation). There are a lot of activities for this day. One of these activities is the Pi Day Sudoku.

Almost exactly one year ago, I blogged about the Pi Day Sudoku 2009 problem in these two posts: The problem is an extended version of a Sudoku problem:
Rules: Fill in the grid so that each row, column, and jigsaw region contains 1-9 exactly once and π [Pi] three times.

(Click to enlarge the picture)

Since it was a competition (closed June 1, 2009), I didn't published any models of this problem when blogging about the problem. Now, a year after, it seems to be a good time to publish them. I then implemented two version, one in MiniZinc, and one in Comet:: Both models use the same approach (details differs however). It is the same as for the plain Sudoku, with two exceptions:
  • The common approach in plain Sudoku is to use the global constraint alldifferent for stating that the values in the rows, columns, and regions should be different. Since there should be 3 occurrences of π in each row, column and region, this approach don't work. As mentioned in Solving Pi Day Sudoku 2009 with the global cardinality constraint, my first approach was a home made version of the constraint alldifferent except 0 and Pi but it was too slow. Instead (and via a suggestion of Mikael Zayenz Lagerkvist) I changed to the global constraint global_cardinality, which was much faster.
  • The other difference to the standard approach is that the regions are not 3x3 (or MxM), but "jigsaw regions". It was not especially hard to make this representation (though boring to type them in). The constraint for checking the regions are (in MiniZinc):
      /\ % the regions
      forall(i in 1..num_regions) (
        check([x[regions[j,1],regions[j,2]] | j in 1+(i-1)*n..1+((i-1)*n)+11 ])
      )
    
Here I will not publish the solution to the puzzle since I gather that there are a lot of people out there that will solve it by there own. And there is at least one valid solution out there.

Summary

This was a fun problem, especially since I learned some new things by implement the models. As a constraint programming challenge is was quite harder than this year puzzle: Pi Day 2010:
Rules: Fill in the grid so that each row, column, and block contains 1-9 exactly once. This puzzle only has 18 clues! That is conjectured to be the least number of clues that a unique-solution rotationally symmetric puzzle can have. To celebrate Pi Day, the given clues are the first 18 digits of π = 3.14159265358979323...
[Yes, I have implemented a MiniZinc model for this problem as well; it is a standard Sudoku problem. No, I will not publish the model or a solution until the deadline, June 1, 2010.]

For more about Pi Day Sudoku 2010, see the blog 360: Pi Day Sudoku is back.

Also, see the following models that implements Killer Sudoku, which use the same approach as Pi Day Sudoku 2009:

The MiniZinc code

Here is the MiniZinc model (sudoku_pi.mzn), slightly edited.
include "globals.mzn";
int: n = 12;
int: X = 0;  % the unknown
int: P = -1; % π (Pi)

predicate check(array[int] of var int: x) =
   global_cardinality(x, occ) % :: domain
;

array[1..n, 1..n] of var -1..9: x :: is_output;
array[-1..9] of 0..3: occ = array1d(-1..9, [3, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1]);
array[1..11] of 0..3: occ2 = [3, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1];

% solve satisfy;
solve :: int_search([x[i,j] | i,j in 1..n], first_fail, indomain_min, complete) satisfy;

constraint

  % copy the hints
  forall(i in 1..n, j in 1..n) (
      x[i,j] != 0
      /\
      if dat[i,j] != X then
        x[i,j] = dat[i,j]
      else 
        true
      endif
  )

  /\ % rows
  forall(i in 1..n) (
    check([x[i,j] | j in 1..n]) 
  )

  /\ % columns
  forall(j in 1..n) (
    check([x[i,j] | i in 1..n])
  )

  /\ % the regions
  forall(i in 1..num_regions) (
    check([x[regions[j,1],regions[j,2]] | j in 1+(i-1)*n..1+((i-1)*n)+11 ])
  )
;

output 
[ show(occ) ++ "\n"] ++
[
  if j = 1 then "\n" else " " endif ++
   show(x[i,j])
  | i in 1..n, j in 1..n

] ++ ["\n"];


% data
array[1..n,1..n] of int: dat = array2d(1..n, 1..n,
[
 4,9,7, P,5,X,X,X,X, X,X,X,
 X,P,X, 8,X,X,9,6,1, 5,2,X,
 X,8,X, 1,X,X,X,P,X, 7,X,X,
 X,X,X, X,X,X,X,P,X, 4,X,X,
 5,3,9, 6,X,X,X,X,X, X,X,X,

 9,4,X, P,P,P,7,X,X, X,X,X,
 X,X,X, X,X,6,2,5,P, X,7,4,
 X,X,X, X,X,X,X,X,P, P,3,8,
 X,7,8, 4,6,9,X,X,X, X,X,X,

 X,X,3, X,P,X,X,4,7, 1,6,9,
 X,X,4, X,1,X,X,X,6, X,P,X,
 X,X,X, X,X,X,X,X,4, X,5,X
 ]);


% The regions
int: num_regions = 12;
array[1..num_regions*12, 1..2] of int: regions  = array2d(1..num_regions*12, 1..2,
[
  % Upper left dark green
  1,1  , 1,2  , 1,3  , 
  2,1  , 2,2  , 2,3  , 
  3,1  , 3,2  , 
  4,1  , 4,2  ,  
  5,1  , 5,2  , 
 
  % Upper mid light dark green
  1,4  ,  1,5  ,  1,6  ,  1,7  ,  1,8  ,  1,9  , 
  2,4  ,  2,5  ,  2,6  ,  2,7  ,  2,8  ,  2,9  , 

  % Upper right green
  1,10  ,  1,11  ,  1,12  , 
  2,10  ,  2,11  ,  2,12  , 
  3,11  ,  3,12  , 
  4,11  ,  4,12  , 
  5,11  ,  5,12   , 

  % Mid upper left "blue"
  3,3  ,  3,4  , 3,5  ,  3,6  , 
  4,3  ,  4,4  , 4,5  ,  4,6  , 
  5,3  ,  5,4  , 5,5  ,  5,6  , 

  % Mid Upper right blue
  3,7  ,  3,8  ,  3,9  ,  3,10  , 
  4,7  ,  4,8  ,  4,9  ,  4,10  , 
  5,7  ,  5,8  ,  5,9  ,  5,10  , 

  % Mid left green
  6,1  ,  6,2  , 6,3  , 
  7,1  ,  7,2  , 7,3  , 
  8,1  ,  8,2  , 8,3  , 
  9,1  ,  9,2  , 9,3  , 

  % Mid left blue
  6,4  , 6,5  , 
  7,4  , 7,5  , 
  8,4  , 8,5  , 
  9,4  , 9,5  , 
  10,4 , 10,5  , 
  11,4 , 11,5  , 

  % Mid mid green
  6,6  , 6,7  , 
  7,6  , 7,7  , 
  8,6  , 8,7  , 
  9,6  , 9,7  , 
  10,6 , 10,7  , 
  11,6 , 11,7  , 

  % Mid right blue
  6,8  ,  6,9  , 
  7,8  ,  7,9  , 
  8,8  ,  8,9  , 
  9,8  ,  9,9  , 
  10,8 ,  10,9  , 
  11,8 ,  11,9  , 

  % Mid right green
  6,10  ,  6,11  ,  6,12  , 
  7,10  ,  7,11  ,  7,12  , 
  8,10  ,  8,11  ,  8,12  , 
  9,10  ,  9,11  ,  9,12  , 

  % Lower left dark green
  10,1  , 10,2  ,  10,3  , 
  11,1  , 11,2  ,  11,3  , 
  12,1  , 12,2  ,  12,3  , 12,4  , 12,5  ,  12,6  , 

  % Lower right dark green
  10,10  ,  10,11  , 10,12  , 
  11,10  ,  11,11  , 11,12  , 
  12,7   ,  12,8   ,  12,9  , 12,10  , 12,11  ,  12,12  
]);

The Comet code

The Comet model (sudoku_pi.co) use the same principle as the MiniZinc model. However, the representation of the regions are different where I instead of a matrix use a more object oriented approach with two tuple for the structures. For some reasons that I have forgot now, I didn't create a function check in this Comet model, instead stated the cardinality constraints directly.
import cotfd;
int t0 = System.getCPUTime();

int n = 12;
int P = -1; // Pi
int X = 0; // unknown 
range R = -1..9; 

set{int} V = {-1,1,2,3,4,5,6,7,8,9};

// regions where 1..9 is alldiff + 3 Pi
tuple Pos {
  int row;
  int col;
}

tuple Region {
  set{Pos} p;
}

int num_regions = 12;
Region regions[1..num_regions] = 
[
 // Upper left dark green
 Region({Pos(1,1),Pos(1,2),Pos(1,3),
         Pos(2,1),Pos(2,2),Pos(2,3),
         Pos(3,1),Pos(3,2),
         Pos(4,1),Pos(4,2), 
         Pos(5,1), Pos(5,2)}),
 
 // Upper mid light dark green
 Region({Pos(1,4), Pos(1,5), Pos(1,6), Pos(1,7), Pos(1,8), Pos(1,9),
         Pos(2,4), Pos(2,5), Pos(2,6), Pos(2,7), Pos(2,8), Pos(2,9)}),

 // Upper right green
 Region({Pos(1,10), Pos(1,11), Pos(1,12),
         Pos(2,10), Pos(2,11), Pos(2,12),
         Pos(3,11), Pos(3,12),
         Pos(4,11), Pos(4,12),
         Pos(5,11), Pos(5,12) }),

 // Mid upper left "blue"
 Region({Pos(3,3), Pos(3,4),Pos(3,5), Pos(3,6),
         Pos(4,3), Pos(4,4),Pos(4,5), Pos(4,6),
         Pos(5,3), Pos(5,4),Pos(5,5), Pos(5,6)}),

 // Mid Upper right blue
 Region({Pos(3,7), Pos(3,8), Pos(3,9), Pos(3,10),
        Pos(4,7), Pos(4,8), Pos(4,9), Pos(4,10),
        Pos(5,7), Pos(5,8), Pos(5,9), Pos(5,10)}),

 // Mid left green
 Region({Pos(6,1), Pos(6,2),Pos(6,3),
         Pos(7,1), Pos(7,2),Pos(7,3),
         Pos(8,1), Pos(8,2),Pos(8,3),
         Pos(9,1), Pos(9,2),Pos(9,3)}),

 // Mid left blue
 Region({Pos(6,4),Pos(6,5),
         Pos(7,4),Pos(7,5),
         Pos(8,4),Pos(8,5),
         Pos(9,4),Pos(9,5),
         Pos(10,4),Pos(10,5),
         Pos(11,4),Pos(11,5)}),

 // Mid mid green
 Region({Pos(6,6),Pos(6,7),
         Pos(7,6),Pos(7,7),
         Pos(8,6),Pos(8,7),
         Pos(9,6),Pos(9,7),
         Pos(10,6),Pos(10,7),
         Pos(11,6),Pos(11,7)}),

 // Mid right blue
 Region({Pos(6,8), Pos(6,9),
         Pos(7,8), Pos(7,9),
         Pos(8,8), Pos(8,9),
         Pos(9,8), Pos(9,9),
         Pos(10,8), Pos(10,9),
         Pos(11,8), Pos(11,9)}),

 // Mid right green
 Region({Pos(6,10), Pos(6,11), Pos(6,12),
         Pos(7,10), Pos(7,11), Pos(7,12),
         Pos(8,10), Pos(8,11), Pos(8,12),
         Pos(9,10), Pos(9,11), Pos(9,12)}),

 // Lower left dark green
 Region({Pos(10,1),Pos(10,2), Pos(10,3),
         Pos(11,1),Pos(11,2), Pos(11,3),
         Pos(12,1),Pos(12,2), Pos(12,3),Pos(12,4),Pos(12,5), Pos(12,6)}),

 // Lower right dark green
 Region({Pos(10,10), Pos(10,11),Pos(10,12),
         Pos(11,10), Pos(11,11),Pos(11,12),
         Pos(12,7),Pos(12,8), Pos(12,9),Pos(12,10),Pos(12,11), Pos(12,12)})

 ];

// the hints
int data[1..n,1..n] = 
[
 [4,9,7, P,5,X,X,X,X, X,X,X],
 [X,P,X, 8,X,X,9,6,1, 5,2,X],
 [X,8,X, 1,X,X,X,P,X, 7,X,X],
 [X,X,X, X,X,X,X,P,X, 4,X,X],
 [5,3,9, 6,X,X,X,X,X, X,X,X],

 [9,4,X, P,P,P,7,X,X, X,X,X],
 [X,X,X, X,X,6,2,5,P, X,7,4],
 [X,X,X, X,X,X,X,X,P, P,3,8],
 [X,7,8, 4,6,9,X,X,X, X,X,X],

 [X,X,3, X,P,X,X,4,7, 1,6,9],
 [X,X,4, X,1,X,X,X,6, X,P,X],
 [X,X,X, X,X,X,X,X,4, X,5,X]
 ];

Integer num_solutions(0);

Solver m();
var{int} x[1..n,1..n](m, R);

// int occ[-1..9] = [3,0,1,1,1,1,1,1,1,1,1];
var{int} occ[-1..9](m,0..3);
int occ_count[-1..9] = [3,0,1,1,1,1,1,1,1,1,1];

exploreall {

  // get the hints
  forall(i in 1..n) {
    forall(j in 1..n) {
      int c = data[i,j];
      if (c == P) {
        cout << "P";
      } else {
        cout << data[i,j];
      }
      if (c != 0) {
        m.post(x[i,j] == data[i,j]);
      }
      cout << " ";
    }
    cout << endl;
   }

  forall(i in 1..n, j in 1..n) {
    m.post(x[i,j] != 0);
  }

  forall(i in -1..9) {
    m.post(occ[i] == occ_count[i]);
  }  
  
  // rows
  forall(i in 1..n) {
    m.post(cardinality(occ, all(j in 1..n) x[i,j])); 
  }
    
  // columns
  forall(j in 1..n) {
    m.post(cardinality(occ, all(i in 1..n) x[i,j]));
  }

  // regions
  forall(r in 1..num_regions) {
    m.post(cardinality(occ, all(i in regions[r].p) x[i.row,i.col])); 

  }

} using {

  // reversing i and j gives faster solution
  forall(i in 1..n, j in 1..n: !x[i,j].bound()) {
    tryall(v in V : x[i,j].memberOf(v)) by(v) 
      m.label(x[i,j], v);
    onFailure
      m.diff(x[i,j], v);
  }

  int t1 = System.getCPUTime();
  cout << "time:      " << (t1-t0) << endl;
  cout << "#choices = " << m.getNChoice() << endl;
  cout << "#fail    = " << m.getNFail() << endl;
  cout << "#propag  = " << m.getNPropag() << endl;
  

  num_solutions++;

   forall(i in 1..n) {
     forall(j in 1..n) {
       int c = x[i,j];
       if (c == P) {
         cout << "P";
       } else {
         cout << x[i,j];
       }
       cout << " ";
     }
     cout << endl;
   }
}

cout << "\nnum_solutions: " << num_solutions << endl;
cout << endl << endl;

int t1 = System.getCPUTime();
cout << "time:      " << (t1-t0) << endl;
cout << "#choices = " << m.getNChoice() << endl;
cout << "#fail    = " << m.getNFail() << endl;
cout << "#propag  = " << m.getNPropag() << endl;

January 14, 2010

Some new MiniZinc models, mostly Enigma problems

The last week I have solved some (more) of New Scientist's Enigma problems. Note: You may have to be a subscriber to read all the articles, however there is a grace period with some free peeks (7 it seems to be).

Enigma problems

Here are the new models of Enigma problems. All are written in MiniZinc.
  • enigma_counting_pennies.mzn: Enigma Counting pennies (from 2005)
    Here all the 24 ways of permuting a list of 4 items are needed. Instead of calculating them, I simply listed all the variants and used them later as index.
  • enigma_843.mzn: Enigma How many are whole? (#843)
    In this model there was problem using the following type of reifications:
      square(ONE) /\ ( square(TWO) \/ square(THREE) \/ square(FOUR)) % /\ ..
      % or
      (square(ONE) /\ ( bool2int(square(TWO)) + bool2int(square(THREE)) + bool2int(square(FOUR)) = 1) % /\ ...
    
    where square is defined as:
    predicate square(var int: y) =
      let {
         var 1..ceil(sqrt(int2float(ub(y)))): z
      } in
       z*z = y;
    
    In some of the later versions of MiniZinc (probably 1.0.2), the handling of reification was changed. My fix was to add a boolean matrix which handled the boolean constraints. It is not beatiful:
       % In the list ONE TWO THREE FOUR just the first and one other are perfect squares.
       square(ONE) /\ 
       [bools[1,j] | j in 1..4] = [square(ONE),square(TWO),square(THREE),square(FOUR)] /\
       sum([bool2int(bools[1,j]) | j in 1..4]) = 2 /\
       % ...
    
  • enigma_1530.mzn: Enigma Tom Daley (#1530)
    A simple alphametic puzzle.
  • enigma_1535.mzn: Enigma Back to front (#1535)
    Magic square with some more constraints.
  • enigma_1555.mzn: Enigma Not a square (#1555)
    Another alphametic puzzle on a grid.
  • enigma_1557.mzn: Enigma Reverse Division (#1557)
    Division and reversing some numbers.
  • enigma_1568.mzn: Enigma Odd puzzle (#1568)
    Long multipication alphametic puzzle.
  • enigma_1570.mzn: Enigma Set of cubes (#1570)
    Pandigital problem over at set of cubes. Here I simply listed the cubes between 0^3 and 1000^3 instead of calculating them. This may be consider cheating...
  • enigma_1575.mzn: Enigma All our days (#1575)
    Another alphametic puzzle with some more constraints.
Some of these problem may have been better solved using other means, e.g. pure (or not so pure) mathematics, but this is after all a blog about constraint programming, and modellng them is good exercise. And fun.

I have also solved - or at least modeled - some of the open problems: 1573, 1574, 1576, and 1577, but will not publish the models until they are closed. I "accidentally" published problem #1574 the other day before I realized it was still open, but so be it.

Some other models/data files

Here are two simple new problems taken from Choco's (forthcoming site ) example distribution: Also, a new data file for the Bridge and Torch model: bridge_and_torch_problem8.dzn, which is, as I understand it, the same as Choco's U2planning problem.

Birthdays 2010 puzzle

Last, here is a very simple puzzle of my own based on a coincidence of birth years of me and my brother and the current year:
This year (2010) my older brother Anders will be the same age as the year I was born in the last century, and I will be the same age as the year he was born in the last century. My brother is four years older than me. How old are we this year?
A (unnecessary complex) MiniZinc model of this problem is birthdays_2010.mzn. Of course, it can be solved much easier with the simple equation: {H+A=110,A-H=4} (link to Wolfram Alpha), or in MiniZinc:
constraint H+A = 110 /\ A-H = 4;
But this latter version is not very declarative, and I - in general - prefer the declarative way.

Update: Added link to Wolfram Alpha for the simple equation above.

January 04, 2010

Finding celebrities at a party

This problem is from Uwe Hoffmann Finding celebrities at a party (PDF):
Problem: Given a list of people at a party and for each person the list of people they know at the party, we want to find the celebrities at the party. A celebrity is a person that everybody at the party knows but that only knows other celebrities. At least one celebrity is present at the party.
(This paper contains an implementation of the problem in Scala.)

Note: The original of this problem is Richard Bird and Sharon Curtis: "Functional pearls: Finding celebrities: A lesson in functional programming" (J. Funct. Program., 16(1):13–20, 2006), but I (nor Hoffmann) have not been able to access this paper. Update: I have now got hold of this paper. Thank you SW!

The example in Hoffmann's paper is to find of who are the celebrity/celebrities in this party graph:
Adam  knows {Dan,Alice,Peter,Eva},
Dan   knows {Adam,Alice,Peter},
Eva   knows {Alice,Peter},
Alice knows {Peter},
Peter knows {Alice}
Solution: the celebrities are Peter and Alice.

MiniZinc model

The MiniZinc model of this problem is finding_celebrities.mzn.

Following are some discussion of the two constraints
  • All must knows a celebrity
  • Celebrities only knows a celebrity
But first a comment how the party graph is represented.

Party graph

Here I have chosen to represent the party as an array of sets:
  Adam  (1) knows {Dan,Eva,Alice,Peter}  {2,3,4,5}
  Dan   (2) knows {Adam,Alice,Peter}     {1,4,5}
  Eva   (3) knows {Alice,Peter}          {4,5}
  Alice (4) knows {Peter}                {5}
  Peter (5) knows {Alice}                {4}
This is coded in MiniZinc as:
party = [
          {2,3,4,5}, % 1, Adam
          {1,4,5},   % 2, Dan 
          {4,5},     % 3, Eva
          {5},       % 4, Alice
          {4}        % 5, Peter
         ];
This corresponds to the following incidence matrix, which is calculated in the model. For simplifying the calculations, we assume that a person know him/herself (this is also handled in the model).
% 1 2 3 4 5
  1,1,1,1,1, % 1
  1,1,0,1,1, % 2
  0,0,1,1,1, % 3
  0,0,0,1,1, % 4
  0,0,0,1,1  % 5
This conversion from incidence sets (party) to incidence matrix (graph) is done by the set2matrix predicate:
predicate set2matrix(array[int] of var set of int: s,
                     array[int,int] of var int: mat) =
     forall(i in index_set(s)) ( graph[i,i] = 1)
     /\
     forall(i,j in index_set(s) where i != j) (
      j in party[i] <-> graph[i,j] = 1
    )

All must know a celebrity

I started this model by consider the celebrities as a clique, and therefore using the constraint clique (which is still included in the model). However, is not enough to identify the clique since there may be other cliques but are not celebrities. In the above example there is another clique: {Dan,Adam}.

In fact, the clique constraint is not needed at all (and it actually makes the model slower). Instead we can just look for person(s) that everybody knows, i.e. where there are all 1's in the column of a celebrity in the party graph matrix (graph in the model). This is covered by the constraint:
forall(i in 1..n) (
   (i in celebrities -> sum([graph[j,i] |j in 1..n]) = n)
)  
This is a necessary condition but not sufficient for identifying celebrities.

Celebrities only know each other

We must also add the constraint that all people in the celebrity clique don't know anyone outside this clique. This is handled by the constraint the all celebrities must knows the same amount of persons, i.e. the size (cardinality) of the clique:
forall(i in 1..n) (
   i in celebrities -> sum([graph[i,j] | j in 1..n]) = num_celebrities
)
As noted above, a person is assumed to know him/herself.

Running the model

If we run the MiniZinc model finding_celebrities.mzn we found the following solution:
celebrities = 4..5;
num_celebrities = 2;
which means that the celebrities are {4,5}, i.e. {Alice, Peter}.

Another, slightly different party

We now change the party matrix slightly by assuming that Alice (4) also know Adam (1), i.e. the following incidence matrix:
 % 1 2 3 4 5
   1,1,1,1,1, % 1
   1,1,0,1,1, % 2
   0,0,1,1,1, % 3
   1,0,0,1,1, % 4
   0,0,0,1,1  % 5
]);
This makes Alice a non celebrity since she knows the non celebrity Adam, and this in turn also makes Peter a non celebrity since he knows the non celebrity Alice. In short, in this party there are no celebrities.

The model also contains a somewhat larger party graph with 10 persons.

Update about 8 hours later There is now a version which uses just set constraints, i.e. no conversion to an incidence matrix: finding_celebrities2.mzn. The constraint sections is now:
constraint
  num_celebrities >= 1

  /\ % all persons know the celebrities,
     % and the celebrities only know celebrities
  forall(i in 1..n) (
     (i in celebrities -> 
             forall(j in 1..n where j != i) (i in party[j]))
     /\
     (i in celebrities -> 
             card(party[i]) = num_celebrities-1)
  )
I have kept the same representation of the party (the array of sets of who knows who) as the earlier model which means that a person is now not assuming to know him/herself. The code reflects this change by using where j != i and num_celebrities-1.

December 29, 2009

1 year anniversary and Secret Santa problem II

Exactly one year ago, I started this blog. Little did I know what to expect. Through this blog and its accompanying pages - I have met many very interesting persons that has learned me much in my pursuit of learning constraint programming. I am very grateful for your help and inspiration! And thanks to all my (other) readers.

I hope the following year will be as rewarding as this last.

Secret Santa problem II

As an anniversary gift, here is another Secret Santa problem (compare with Merry Christmas: Secret Santas Problem) with a slightly different touch.

The problem formulation is from Maple Primes forum Secret Santa Graph Theory:
Every year my extended family does a "secret santa" gift exchange. Each person draws another person at random and then gets a gift for them. At first, none of my siblings were married, and so the draw was completely random. Then, as people got married, we added the restriction that spouses should not draw each others names. This restriction meant that we moved from using slips of paper on a hat to using a simple computer program to choose names. Then people began to complain when they would get the same person two years in a row, so the program was modified to keep some history and avoid giving anyone a name in their recent history. This year, not everyone was participating, and so after removing names, and limiting the number of exclusions to four per person, I had data something like this:

Name: Spouse, Recent Picks

Noah: Ava. Ella, Evan, Ryan, John
Ava: Noah, Evan, Mia, John, Ryan
Ryan: Mia, Ella, Ava, Lily, Evan
Mia: Ryan, Ava, Ella, Lily, Evan
Ella: John, Lily, Evan, Mia, Ava
John: Ella, Noah, Lily, Ryan, Ava
Lily: Evan, John, Mia, Ava, Ella
Evan: Lily, Mia, John, Ryan, Noah
I have interpreted the problem as follows:
  • one cannot be a Secret Santa of one's spouse nor of oneself
  • one cannot be a Secret Santa for somebody two years in a row
  • objective: maximize the "Secret Santa distance", i.e. the number of years since the last assignment of the same person
My MiniZinc model for this problem is secret_santa2.mzn.

This is a more traditional linear programming problem compared to Secret Santa I, using a distance matrix for maximizing the "Secret Santa Distance". M is a "large" number (number of persons + 1) for coding that there have been no previous Secret Santa assignment.
rounds = array2d(1..n, 1..n, [
%N  A  R  M  El J  L  Ev 
 0, M, 3, M, 1, 4, M, 2, % Noah
 M, 0, 4, 2, M, 3, M, 1, % Ava 
 M, 2, 0, M, 1, M, 3, 4, % Ryan
 M, 1, M, 0, 2, M, 3, 4, % Mia 
 M, 4, M, 3, 0, M, 1, 2, % Ella
 1, 4, 3, M, M, 0, 2, M, % John
 M, 3, M, 2, 4, 1, 0, M, % Lily
 4, M, 3, 1, M, 2, M, 0  % Evan
]);

The original problem don't say anything about single persons, i.e. those without spouses. In this model, singleness (no-spouseness) is coded as spouse = 0, and the no-spouse-Santa constraint has been adjusted to takes care of this.

The constraint part is the following, where n is the number of persons:
constraint
   all_different(santas)

   /\ % no Santa for one self or the spouse
   forall(i in 1..n) (
      santas[i] != i /\
      if spouses[i] > 0 then santas[i] != spouses[i] else true endif
   )

   /\ % the "santa distance" (
   forall(i in 1..n) ( santa_distance[i] = rounds[i,santas[i]] )

   /\ % cannot be a Secret Santa for the same person two years in a row.
   forall(i in 1..n) (
      let { var 1..n: j } in
       rounds[i,j] = 1 /\ santas[i] != j
   )

   /\
   z = sum(santa_distance)

Solution

This model gives - when using solve satisfy and the constraint z >= 67 - the following 8 solutions with a total of Secret Santa distance of 67 (sum(santa_distance)). If all Secret Santa assignments where new it would have been a total of n*(n+1) = 8*9 = 72. As we can see there is always one Santa with a previous existing assignment. (With one Single person and the data I faked, we can get all brand new Secret Santas. See the model for this.)
santa_distance: 9 9 4 9 9 9 9 9
santas        : 7 5 8 6 1 4 3 2

santa_distance: 9 9 4 9 9 9 9 9
santas        : 7 5 8 6 3 4 1 2

santa_distance: 9 9 9 4 9 9 9 9
santas        : 7 5 6 8 1 4 3 2

santa_distance: 9 9 9 4 9 9 9 9
santas        : 7 5 6 8 3 4 1 2

santa_distance: 9 9 9 9 4 9 9 9
santas        : 4 7 1 6 2 8 3 5

santa_distance: 9 9 9 9 4 9 9 9
santas        : 4 7 6 1 2 8 3 5

santa_distance: 9 9 9 9 9 9 4 9
santas        : 4 7 1 6 3 8 5 2

santa_distance: 9 9 9 9 9 9 4 9
santas        : 4 7 6 1 3 8 5 2

December 25, 2009

Merry Christmas: Secret Santas Problem

Here is a fun little problem related to the holiday. Merry Christmas, everyone! (For the Swedish readers: Sorry for the one day off greeting.)

This problem is from the Ruby Quiz#2 Secret Santas
Honoring a long standing tradition started by my wife's dad, my friends all play a Secret Santa game around Christmas time. We draw names and spend a week sneaking that person gifts and clues to our identity. On the last night of the game, we get together, have dinner, share stories, and, most importantly, try to guess who our Secret Santa was. It's a crazily fun way to enjoy each other's company during the holidays.

To choose Santas, we use to draw names out of a hat. This system was tedious, prone to many "Wait, I got myself..." problems. This year, we made a change to the rules that further complicated picking and we knew the hat draw would not stand up to the challenge. Naturally, to solve this problem, I scripted the process. Since that turned out to be more interesting than I had expected, I decided to share.

This weeks Ruby Quiz is to implement a Secret Santa selection script.

Your script will be fed a list of names on STDIN.

...

Your script should then choose a Secret Santa for every name in the list. Obviously, a person cannot be their own Secret Santa. In addition, my friends no longer allow people in the same family to be Santas for each other and your script should take this into account.
The MiniZinc model secret_santa.mzn skips the parts of input and mailing. Instead, we assume that the friends are identified with a unique number from 1..n, and the families are identified with a number 1..num_families.

We use two arrays:
  • the array x represents whom a person should be a Santa of. x[1] = 10 means that person 1 is a Secret Santa of person 10, etc.
  • the family array consists of the family identifier of each person.
Now, the three constraints can easily be stated in a constraint programming system like MiniZinc:
  • "everyone gives and received a Secret Santa gift": this is handled with a permutation of the values 1..n using all_different(x).
  • "one cannot be one own's Secret Santa". This is captured in the no_fix_point predicate, stating that there can be no i for which x[i] = i (i.e. no "fix point").
  • "no Secret Santa to a person in the same family". Here we use the family array and checks that for each person (i), the family of i (family[i]) must not be the same as the family of the person that receives the gift (family[x[i]]).
Here is the complete MiniZinc model (in a slightly compact form):
include "globals.mzn"; 
int: n = 12;
int: num_families = 4;
array[1..n] of 1..num_families: family = [1,1,1,1, 2, 3,3,3,3,3, 4,4];
array[1..n] of var 1..n: x :: is_output;

% Ensure that there are no fix points in the array.
predicate no_fix_points(array[int] of var int: x) = 
      forall(i in index_set(x)) ( x[i] != i  );

solve satisfy;

constraint
  % Everyone gives and receives a Secret Santa
  all_different(x) /\      

  % Can't be one own's Secret Santa
  no_fix_points(x) /\   

  % No Secret Santa to a person in the same family
  forall(i in index_set(x)) (   family[i] != family[x[i]]  )
; 

% output (just for the minizinc solver)
output [
   "Person " ++ show(i) ++ 
   " (family: " ++ show(family[i]) ++ ") is a Secret Santa of " ++ 
    show(x[i]) ++ 
   " (family: " ++ show(family[x[i]]) ++ ")\n"
   | i in 1..n
] ++ 
["\n"];

Solution

Here is the first solution (of many):
[10, 9, 8, 5, 12, 4, 3, 2, 1, 11, 7, 6]
This means that person 1 should be a Secret Santa of person 10, etc.

The minizinc solver gives the following, using the output code (slightly edited):
Person  1 (family: 1) is a Secret Santa of 10 (family: 3)
Person  2 (family: 1) is a Secret Santa of  9 (family: 3)
Person  3 (family: 1) is a Secret Santa of  8 (family: 3)
Person  4 (family: 1) is a Secret Santa of  5 (family: 2)
Person  5 (family: 2) is a Secret Santa of 12 (family: 4)
Person  6 (family: 3) is a Secret Santa of  4 (family: 1)
Person  7 (family: 3) is a Secret Santa of  3 (family: 1)
Person  8 (family: 3) is a Secret Santa of  2 (family: 1)
Person  9 (family: 3) is a Secret Santa of  1 (family: 1)
Person 10 (family: 3) is a Secret Santa of 11 (family: 4)
Person 11 (family: 4) is a Secret Santa of  7 (family: 3)
Person 12 (family: 4) is a Secret Santa of  6 (family: 3)

Bales of Hay

As an extra, here is another MiniZinc model: bales_of_hay.mzn which solves the following problem (from The Math Less Traveled the other day):
You have five bales of hay.

For some reason, instead of being weighed individually, they were weighed in all possible combinations of two. The weights of each of these combinations were written down and arranged in numerical order, without keeping track of which weight matched which pair of bales. The weights, in kilograms, were 80, 82, 83, 84, 85, 86, 87, 88, 90, and 91.

How much does each bale weigh? Is there a solution? Are there multiple possible solutions?
The answer? There is a unique solution (when bales are ordered on weight): 39, 41, 43, 44, 47.

November 08, 2009

Update on Nonogram: Jan Wolter's Survey and my own new benchmark

Survey of Paint-by-Number Puzzle Solvers

In Some new models, etc I mentioned the great Survey of Paint-by-Number Puzzle Solvers, created by Jan Wolter (also author of the Nonogram solver pbnsolve).

In this survey he included both Gecode's Nonogram solver written by Mikael Lagerkvist as well as my own Nonogram model (with Gecode/FlatZinc).

Since the last update on the former blog post the following has happeded:
  • both our solver has now the assessment: "An amazingly good solver, especially for a simple demo program", and are placed 4, 5, and 6 of 10 tested systems
  • my Gecode/FlatZinc model has been tested for "Full results"; it came 4 out of 5.
  • my Nonogram model with the Lazy FD solver is now included in the "Sample result", at 6'th place
It seems than Wolter has appreciated constraint programming as a general tool for solving these kind of combinatorial problems, much for its ease of experimentation, e.g. with labeling strategies and (for the MiniZinc models) changing solvers:

From the analysis of Lagerkvist's Gecode model:
This is especially impressive because the faster solvers are large, complex programs custom designed to solve paint-by-number puzzles. This one is a general purpose solver with just a couple hundred custom lines of code to generate the constraints, run the solver, and print the result. Considering that this is a simple application of a general purpose solving tool rather than hugely complex and specialized special purpose solving tool, this is an astonishingly good result.

Getting really first class search performance usually requires a lot of experimentation with different search strategies. This is awkward and slow to do if you have to implement each new strategies from scratch. I suspect that a tool like Gecode lets you try out lots of different strategies with relatively little coding to implement each one. That probably contributes a lot to getting to better solvers faster.
From the analysis of my MiniZinc model:
If you tried turning this into a fully useful tool rather than a technology demo, with input file parsers and such, it would get a lot bigger, but clearly the constraint programming approach has big advantages, achieving good search results with low development cost.

...

These two results [Gecode/FlatZinc and LazyFD] highlight the advantage of a solver-independent modeling language like MiniZinc. You can describe your problem once, and then try out a wide variety of different solvers and heuristics without having to code every single one from scratch. You can benefit from the work of the best and the brightest solver designers. It's hard to imagine that this isn't where the future lies for the development of solvers for applications like this.
And later in the Conclusions:
The two constraint-based systems by Lagerkvist and Kjellerstrand come quite close in performance to the dedicated solvers, although both are more in the category of demonstrations of constraint programming than fully developed solving applications. The power of the underlaying search libraries and the ease of experimentation with alternative search heuristics obviously serves them well. I think it very likely that approaches based on these kinds of methods will ultimately prove the most effective.
I think this is an important lesson: Before starting to write very special tools, first try a general tool like a constraint programming system and see how well it perform.

The Lazy FD solver and the Lion problem

Most of the problems in the Sample Results where solved by some solver within the time limit of 30 minutes. However, one problem stand out as extra hard: The Lion problem. When I tested the MiniZinc's Lazy FD solver on my machine I was very excited that it took just about 2 minutes, and mentioned this to Wolter. He also tested this but on his 64-bit machine it took 80 minutes to solve (and since it was above the time limit this is not in the result table). This is how he describes the Lazy FD solver:
But the remarkable thing is that [the Lazy FD solver] solves almost everything. Actually, it even solved the Lion puzzle that no other solver was able to solve, though, on my computer, it took 80 minutes to do it. Now, I didn't allow a lot of other solvers 80 minutes to run, but that's still pretty impressive. (Note that Kjellerstrand got much faster solving times for the Lion than I did. Lagerkvist also reported that his solver could solve it, but I wasn't able to reproduce that result even after 30 CPU hours. I don't know why.)
After some discussion, we come to the conclusion that the differences was probably due to the fact that I use a 32-bit machine (and the 32-bit version of MiniZinc) with 2 Gb memory, and Wolter use a 64-bit machine with 1 Gb memory.

One should also note that the all other solvers was compiled without optimization, including Gecode/FlatZinc; however the LazyFD does not come with source so it is running optimized. This may be an unfair advantage to the LazyFD solver.

My own benchmark of the Sample Results

The times in the Sample Results is, as mentioned above, for solvers compiled with no optimization. I have now run the same problems on my machine (Linux Mandriva, Intel Dual 3.40GHz, with 2Gb memory), but the solvers uses the standard optimization. All problems was run with a time limit of 10 minutes (compared to Wolters 30 minutes) and searched for 2 solutions, which checks for unique solutions. The last three problems (Karate, Flag, Lion) has multiple solutions, and it is considered a failure if not two where found in the time limit. I should also note that during the benchmark I am using the machine for other things, such as surfing etc.

The problems
I downloaded the problems from Wolter's Webbpbn: Puzzle Export. For copyright reasons I cannot republish these models, but it is easy to download each problem. Select ".DZN" for the MiniZinc files, and "A compiled in C++ format" for Gecode. There is no support for Comet's format, but it's quite easy to convert a .dzn file to Comet's.

The solvers + labeling strategies
Here is a description of each solver and its labeling strategy:
  • fz, "normal" (column_row)
    MiniZinc model with Gecode/FlatZinc. The usual labeling in nonogram_create_automaton2.mzn, i.e. where the columns are labeled before rows:
    solve :: int_search(
          [x[i,j] | j in 1..cols, i in 1..rows], 
          first_fail, 
          indomain_min, 
          complete) 
    satisfy;
    
  • fz, "row_column"
    MiniZinc model with Gecode/FlatZinc. Here the order of labeling is reversed, rows are labeled before columns. Model is nonogram_create_automaton2_row_column.mzn
    solve :: int_search(
          [x[i,j] | i in 1..rows, j in 1..cols], 
          first_fail, 
          indomain_min, 
          complete) 
    satisfy;
    
  • fz, "mixed"
    MiniZinc model with Gecode/FlatZinc: nonogram_create_automaton2_mixed.mzn.
    I have long been satisfied with the "normal" labeling in the MiniZinc model because P200 (the hardest problem I until now has tested) was solved so fast. However, the labeling used in the Comet Nonogram model described in Comet: Nonogram improved: solving problem P200 from 1:30 minutes to about 1 second, and which is also used in the Gecode model, is somewhat more complicated since it base the exacl labeling by comparing the hints for the rows and the column.

    I decided to try this labeling in MiniZinc as well. However, labeling in MiniZinc is not so flexible as in Comet and Gecode. Instead we have to add a dedicated array for the labeling (called labeling):
    array[1..rows*cols] of var 1..2: labeling;
    
    and then copy the element in the grid to that array based on the relation between rows and column hints:
    constraint
          % prepare for the labeling
          if rows*row_rule_len < cols*col_rule_len then
               % label on rows first
               labeling = [x[i,j] | i in 1..rows, j in 1..cols]
          else 
               % label on columns first
               labeling = [x[i,j] | j in 1..cols, i in 1..rows]
          endif
          /\
          % .... 
    
    and last, the search is now based just on this labeling array:
    solve :: int_search(
            labeling,
            first_fail, 
            indomain_min, 
            complete)
    satisfy;
    
  • jacop, "normal"
    MiniZinc normal model with JaCoP/FlatZinc using the same model as for fz "normal".
  • lazy, satisfy
    Model: nonogram_create_automaton2_mixed.mzn. This use the MiniZinc LazyFD solver with the search strategy:
    solve satisfy;
    
    This labeling is recommended by the authors of LazyFD. See MiniZinc: the lazy clause generation solver for more information about this solver.

    Note: The solver in MiniZinc latest official version (1.0.3) don't support set vars. Instead I (and also Jan Wolter) used the latest "Release Of The Day" version (as per 2009-11-02).
  • Comet, normal
    Model: nonogram_regular.co. This is the Comet model I described in Comet: Nonogram improved: solving problem P200 from 1:30 minutes to about 1 second. No changes has been done.
  • Gecode, normal
    This is the Nonogram model distributed with Gecode version 3.2.1. The labeling is much like the one used in the Comet model, as well as fz, "mixed". (In fact the labeling in the Gecode model was inspired by the labeling in the Comet model).


Here is the results. For each model (+ labeling strategy) two values are presented:
  • time (in seconds)
  • number of failures if applicable (the LazyFD solver always returns 0 here).
The result
Model fz
normal
fz
row_column
fz
mixed
jacop
normal
lazy
satisfy
Comet
normal
Gecode
normal
Dancer
(#1)
0.48s
(0)
0.31s
(0)
1.00s
(0)
3.64s
(0)
0.91s
(0)
0.691s
(0)
0.199s
(0)
Cat
(#6)
0.24s
(0)
0.24s
(0)
0.25s
(0)
1.20s
(0)
1.13s
(0)
0.6s
(0)
0.025s
(0)
Skid
(#21)
0.24s
(13)
0.23s
(3)
0.28s
(13)
0.78s
(13)
1.37s
(0)
0.586s
(0)
0.217s
(0)
Bucks
(#27)
0.32s
(3)
0.32s
(9)
0.37s
(3)
0.96s
(3)
2.37s
(0)
1.366s
(5)
0.026s
(2)
Edge
(#23)
0.16s
(25)
0.17s
(25)
0.18s
(25)
0.59s
(25)
0.31s
(0)
0.521s
(43)
0.175s
(25)
Smoke
(#2413)
0.27s
(5)
0.27s
(8)
0.28s
(8)
0.83s
(5)
1.44s
(0)
0.616s
(14)
0.275s
(5)
Knot
(#16)
0.42s
(0)
0.43s
(0)
0.48s
(0)
1.19s
(0)
8.15s
(0)
1.307s
(0)
0.329s
(0)
Swing
(#529)
0.95s
(0)
0.94s
(0)
0.96s
(0)
2.19s
(0)
21.94s
(0)
1.782s
(0)
0.431s
(0)
Mum
(#65)
0.64s
(20)
0.64s
(22)
0.66s
(22)
1.68s
(20)
16.34s
(0)
1.268s
(39)
0.491s
(22)
Tragic
(#1694)
340.32s
(394841)
1.02s
(255)
436.92s
(394841)
--
(198329)
45.97s
(0)
477.39s
(702525)
1.139s
(255)
Merka
(#1611)
--
(361587)
1.44s
(79)
--
(294260)
--
(136351)
80.92s
(0)
1.654s
(46)
0.645s
(13)
Petro
(#436)
2.97s
(1738)
143.09s
(106919)
3.42s
(1738)
7.27s
(1738)
9.86s
(0)
3.103s
(3183)
151.09s
(106919)
M_and_m
(#4645)
1.41s
(89)
601.27s
(122090)
1.59s
(89)
3.43s
(89)
66.98s
(0)
2.215s
(162)
2.797s
(428)
Signed
(#3541)
1.87s
(929)
23.12s
(6484)
28.23s
(6484)
5.75s
(929)
73.02s
(0)
20.369s
(12231)
1.648s
(929)
Light
(#803)
600.47s--
(400660)
--
(621547)
--
(485056)
601.53s--
(171305)
8.64s
(0)
--
(0)
--
(538711)
Forever
(#6574)
4.14s
(17143)
7.86s
(30900)
6.22s
(17143)
12.21s
(17143)
3.27s
(0)
7.5s
(27199)
8.077s
(30900)
Hot
(#2040)
--
(303306)
--
(330461)
--
(248307)
--
(119817)
165.72s
(0)
--
(0)
--
(312532)
Karate
(#6739)
95.78s
(215541)
67.27s
(130934)
133.43s
(215541)
373.02s
(215541)
19.32s
(0)
120.02s
(272706)
80.56s
(170355)
Flag
(#2556)
--
(1686545)
5.69s
(14915)
7.93s
(14915)
--
(243222)
9.29s
(0)
7.28s
(24678)
3.998s
(16531)
Lion
(#2712)
--
(542373)
--
(1124697)
--
(420216)
--
(187215)
115.56s
(0)
--
(0)
--
(869513)

Some conclusions, or rather notes

Here are some conclusions (or notes) about the benchmark.
  • The same solver Gecode/FlatZinc is here compared with three different labelings. There is no single labeling that is better than the other. I initially has some hopes that the "mixed" labeling should take the best labeling from the two simpler row/columns labelings, but this is not really the case. For example for Tragic the row_column strategy is better than "normal" and "mixed". I am, however, somewhat tempted, to use the "row_column" labeling, but the drawback is that "P200" problem (not included in Wolfter's sample problems) takes much longer with this labeling.
  • The same model and labeling but with different solvers is compared: Gecode/FlatZinc is faster than JaCoP/FlatZinc on all the problems. For the easier problems this could be explained by the extra startup time of Java for JaCoP, but that is not the complete explanation for the harder problems. Note: Both Gecode/FlatZinc and JaCoP/FlatZinc has dedicated and fast regular constraints (whereas the LazyFD, and the Comet solvers use a decomposition).
  • The LazyFD solver is the only one that solves all problems (including Lion), but is somewhat slower on the middle problems than most of the others. It emphasizes that this is a very interesting solver.
  • It is also interesting to compare the result of the Comet model and Gecode/FlatZinc "mixed", since they use the same principle of labeling. However there are some differences. First, the MiniZinc model with Gecode/FlatZinc use a dedicated regular constraint, and Comet use my own decomposition of the constraint. For the Merka problem the Comet version outperforms the Gecode/FlatZinc version, otherwise it's about the same time (and number of failures).
  • The Light problem: It is weird that this problem was solved in almost exactly 10 minutes (the timeout is 10 minutes) for Gecode/FlatZinc and JaCoP/FlatZinc. The solutions seems correct but I'm suspicious of this. Update: Christian Schulte got me on the right track. Here is was happened: The first (unique) solution was found pretty quick and was printed, but the solvers could not prove a unique solution so it timed out. JaCoP/FlatZinc actually printed "TIME-OUT" but I didn't observe that. Case closed: They both FAILED on this test. Thanks, Christian.End update
As said above, I can only agree with Jan Wolter in his comment that the ease of experimenting, for example changing labeling, and solver for the FlatZinc solvers, is a very nice feature.

Last word

No benchmark or comparison of (constraint programming) models is really complete without the reference to the article On Benchmarking Constraint Logic Programming Platforms. Response to Fernandez and Hill's "A Comparative Study of Eight Constraint Programming Languages over the Boolean and Finite Domains" by Mark Wallace, Joachim Schimpf, Kish Shen, Warwick Harvey. (Link is to ACM.)

From the Abstract:
... The article analyses some pitfalls in benchmarking, recalling previous published results from benchmarking different kinds of software, and explores some issues in comparative benchmarking of CLP systems.

October 31, 2009

Some new models, etc

The readers that follows me on Twitter (hakankj) have probably already has seen this, but here follows a list of models, etc, not yet blogged about.

MiniZinc: Different models

The following four models are translation of my Comet models: And here are some new models:

Nonogram related

A large 40x50 Nonogram problem instance in MiniZinc: nonogram_stack_overflow.dzn to be used with nonogram_create_automaton2.mzn model. The problem was mentioned in Stack Overflow thread Solving Nonograms (Picross). It is solved in under 1 second and 0 failures.

Today my Nonogram model nonogram_create_automaton2.mzn was included in the great Survey of Paint-by-Number Puzzle Solvers (created by Jan Wolter).

My MiniZinc model is described and analyzed here. I'm not at all surprised that it's slower compared to the other solvers; it was quite expected.

Some comments:
Assessment: Slow on more complex puzzles.

...

Results: Run times are not really all that impressive, especially since it is only looking for one solution, not for two like most of the other programs reviewed here. I don't know what the differences are between this and Lagerkvist's system, but this seems much slower in all cases, even though both are ultimately being run in the Gecode library.
Update 2009-11-01
I later realized two things:

1) That the mzn2fzn translator did not use the -G gecode flag, which means that Gecode/FlatZinc uses a decomposition of the regular constraint instead of Gecode's built in, which is really the heart in this model. The model is basically two things: build an automata for a pattern and then run regular on it.
2) When Jan compiled Gecode, he set off all optmization for comparison reason. This is quite unfortunate since Gecode is crafted with knowledge of the specific optimizations.

I there have run all problems by myself and see how well it would done (at least the ballpart figure) when using an "normal" optimized Gecode and with the -G gecode flag for mzn2fzn.

Explanation of the values:
Problem: Name of the problem instance
Runtime: The value of runtime from Gecode/FlatZinc solver
Solvetime: The value of solvetime from Gecode/FlatZinc solver
Failures: Number of failures:
Total time: The Unix time for running the complete problem, including the time of mzn2fzn (which was not included in the benchmark).
A "--" means that a solution was not reached in 10 minutes.
Problem Runtime Solvetime Failures Total time
Dancer 0.002 0.000 0 1.327s
Cat 0.009 0.002 0 0.965s
Skid 0.012 0.006 13 0.660s
Bucks 0.015 0.004 3 0.866s
Edge 0.008 0.005 25 0.447s
Smoke 0.011 0.004 5 0.963s
Knot 0.026 0.006 0 1.450s
Swing 0.059 0.012 0 1.028s
Mum 0.120 0.093 20 1.811s
Tragic 6:24.273 6:23.607 394841 6:28,10
Merka -- -- -- --
Petro 2.571 2.545 1738 4.071s
M&M 0.591 0.510 89 1.961s
Signed 1.074 1.004 929 2.461s
Hot -- -- -- --
Flag -- [lazy solver: 2 solutions in 10 seconds] -- -- --
Lion -- -- -- --

It's interesting to note that the Lazy solver finds some solutions quite fast for the "Flag" problem. However, there where no other big differences compared to Gecode/FlatZinc. I also tested the problems with JaCoP's FlatZinc solver which solved the problems in about the same time as Gecode/FlatZinc with no dramatic differences.

As mentioned above, the exact values is not really comparable to the benchmark values. But it should give an indication of the result when using -G gecode and a "normal" optimized Gecode.

Unfortunately, I cannot link to the specific models due to copyright issues, but they can all be downloaded from the page Web Paint-by-Number Puzzle Export.

(End update)

Update 2009-11-02
Jan Wolter now have rerun the tests of my solver with the -G gecode option, and the time is much more like mine in the table above. The analysis is quite different with an assessment of Pretty decent, and the following under Result:
...

When comparing this to other solvers, it's important to note that nonogram_create_automaton2.mzn contains only about 100 lines of code. From Kjellerstrand's blog, it is obvious that these 100 lines are the result of a lot of thought and experimentation, and a lot of experience with MiniZinc, but the Olšák solver is something like 2,500 lines and pbnsolve approaches 7,000. If you tried turning this into a fully useful tool rather than a technology demo, with input file parsers and such, it would get a lot bigger, but clearly the CSP approach has big advantages.

(End update 2)

Also, the Gecode nonogram solver is included in the survey: called Lagerkvist. I'm not sure when it was added to the survey. It use the latest version of Gecode, so it must have been quite recent.

Some comments:
Assessment: Pretty decent.

...

Puzzles with blank lines seem to cause the program to crash with a segmentation fault.

Otherwise it performs quite well. There seems to be about 0.02 seconds of overhead in all runs, so even simple puzzles take at least that long to run. Aside from that, it generally outperforms the Wilk solver. It's not quite in the first rank, especially considering that it was only finding one solution, not checking for uniqueness, but it's still pretty darn good.

...
I mailed Jan today about the -solutions n options in the Gecode related solvers (he tested my MiniZinc model with Gecode/FlatZinc), as well as some other comments about my model. Also, I will download the problems tested and play with them.

Tailor/Essence': Zebra puzzle

Suddenly I realized that there where no Zebra problem, neither at Andrea's or here. So here it is: zebra.eprime: Zebra puzzle.

Gecode: Words square problem

See Gecode: Modeling with Element for matrices -- revisited for an earlier discussion of this problem.

Thanks to Christian Schulte my Word square model is now much faster: word_square2.cpp.

From the comment in the model:
Christian Schulte suggested 
using branch strategy INT_VAL_SPLIT_MIN 
instead of INT_VAL_MAX .
This made an improvement for size 7 from
322 failures to 42 and from 1:16 minutes
to 10 seconds (for 1 solution).
 
Now it manage to solve a size 8 in a reasonable 
time (1:33 minutes and 1018 failures):
  abastral
  bichrome
  achiotes
  shippers
  troponin
  rotenone
  amerinds
  lessness
 
But wait, there's more!

In the SVN trunk version of Gecode, there is now a version of this model: examples/word-square.cpp where Christian Schulte and Mikael Zayenz Lagerkvist has done a great job of improving it (using a slighly smaller word list, though). It solves a size 8 problem in about 14 seconds. There is also two different approaches with different strengths, etc. I have great hopes that it will improved even further...

October 14, 2009

MiniZinc version 1.0.3 released

MiniZinc version 1.0.3 has been released. It can be downloaded here.

From the NEWS:


G12 MiniZinc Distribution version 1.0.3
---------------------------------------

Bugs fixed in this release:

* A fencepost error that was being introduced into flattened array access
reifications has been fixed.

* Common subexpression elimination has been improved in order to eliminate
redundant int and float linear equations during flattening.

* A bug that caused flattening to abort if array_*_element built-ins were
redefined has been fixed. [Bug #82]

* A bug in the implementation of the FlatZinc set_lt and set_gt built-ins
has been fixed. Note that the expected outputs for the corresponding
tests in the FCTS were also previously incomplete.

* The omission of the string_lit tag from the XML-FlatZinc DTD has been
corrected.

October 08, 2009

MiniZinc: All my public MiniZinc models are now at G12 Subversion repository

All my public MiniZinc models (and data files) are now in the G12 SVN MiniZinc examples repository: http://www.g12.cs.mu.oz.au/mzn, subdirectory hakank/.

The file hakank/README states the following:
In this directory I have collected all the MiniZinc
models, data files, and tools that are available from
   http://www.hakank.org/minizinc/

Any new models on the site will be transferred to this
SVN directory as soon as possible.

Hakan Kjellerstrand, hakank@bonetmail.com
http://www.hakank.org/minizinc/
This means that the models (and other files) which is published at My MiniZinc page is the master, and I will put the files to the G12 SVN repository as soon as possible, hopefully directly after.

The structure in the repository is exactly the same as for the web site, which means that all models are directly in the hakank directory, and then two specific data collections in nonogram_examples and sudoku_problems

Some statistics of the collection (as of writing):
* .mzn files: 742
* .dzn files: 234
* .pl files (Perl program): 2
* README files: 1
* index.html files: 3
* .zip files: 2
*.java files: 1

For a total of 985 files.

I hope there are some use of them.

Also, see The MiniZinc Wiki.

October 01, 2009

MiniZinc: 151 new Nonogram problem instances (from JaCoP)

The latest versions of JaCoP (download here) includes 151 Nonogram instances (ExamplesJaCoP/nonogramRepository/). I wrote about these and the included Nonogram solver in JaCoP: a request from the developers (Knapsack and Geost) and Nonogram labeling.

When I (beta) tested the new FlatZinc support for this version, I converted these instances to MiniZinc's data file (.dzn) for some tests. It is a fun way of testing a system.

Now all these problem instances have been published. They are to be used with the MiniZinc Nonogram model nonogram_create_automaton2.mzn. See At last 2: A Nonogram solver using regular written in "all MiniZinc" for more about this model.

All problem instances - in MiniZinc's .dzn format - are available in www.hakank.org/minizinc/nonogram_examples, and also packed in the Zip file nonogram_examples.zip. The name of each file corresponds to the files in JaCoP's ExamplesJaCoP/nonogramRepository/; the file data000.nin is here called nonogram_jacop_data000.dzn, etc. (A larger file, nonogram_examples_with_fzn.zip also includes the generated .fzn files. Note: I used mzn2fzn for MiniZinc version ROTD/2009-09-13 for this. See below how to generate these files.)

Batch running with JaCoP/Fz2jacop

Running many FlatZinc models with JaCoP's Fz2jacop is not optimal because of the overhead of the Java startup. Instead I used a Java program, JaCoPFznTest.java (based on a program from Krzysztof Kuchcinski, one of JaCoP's main developers. Thanks, Kris.).

The main call is

fz.main(new String[] {"-s", "-a", m});

which means that statistics should be shown and also require to return all solutions of the problem.

The FlatZinc (.fzn) files was generated from the .dzn files with the following command:

foreach_file 'mzn2fzn -G jacop --data $f -o $f.fzn nonogram_create_automaton2.mzn' '*.dzn'

ExamplesJaCoP.Nonogram vs. JaCoP/Fz2jacop

This time I just compared the run time of the two different approaches for solving Nonograms with JaCoP:
* The Java program ExamplesJaCoP.Nonogram
* Running JaCoP's MiniZinc/FlatZinc solver Fz2jacop

The program instances is, with one exception the same as in the distributed ExamplesJaCoP.Nonogram. Since I wanted both program to solve for all solutions, I excluded instance #83 (see below) since it has too many solutions. However, the P200 problem in included since it is hardwired in ExamplesJaCoP.Nonogram. This means that it it still 151 problems running.

The result:
* ExampleJaCoP.Nonogram took 14.8 seconds
* The Fz2jacop version took 17.8 seconds

I'm quite impressed with both of these results, especially Fz2jacop's As shown in JaCoP: a request from the developers (Knapsack and Geost) and Nonogram labeling, many problems are solved in "0" milliseconds, but still.

Instance #83

As mentioned above, problem instance #83 was excluded since it has too many solutions. Here are two of them. It looks like a person (alien?) standing on a spaceship waving.
                              #           #         
                            #       #               
                      # # # #       # # #           
                          #   # # #   # # #         
                  #                   #     #       
                #         #     # # # # #           
              #     #   #           # # #     #     
              #   #   #             # # #   #       
              # #   #               # # #           
              # #                   #   #           
                #                   #   #           
          # # # # # # # # # # # # # # # # #         
      # # # # # # # # # # # # # # # # # # # # #     
    # # # #     # #     # # #     # #     # # # #   
  # # # # #     # #     # # #     # #     # # # # # 
    # # # #     # #     # # #     # #     # # # #   
      # # # # # # # # # # # # # # # # # # # # #     
          # # # # # # # # # # # # # # # # #         
              # # # # # # # # # # # # #             
                  #       #       #                 
                # #       #       # #               
                #         #         #               
              # #         #         # #             
              #           #           #             
            # # #       # # #       # # #           
----------
                              #           #         
                            #       #               
                      # # # #       # # #           
                          #   # # #   # # #         
                  #                   #     #       
                #         #     # # # # #           
              #     #   #           # # #   #       
              #   #   #             # # #     #     
              # #   #               # # #           
              # #                   #   #           
                #                   #   #           
          # # # # # # # # # # # # # # # # #         
      # # # # # # # # # # # # # # # # # # # # #     
    # # # #     # #     # # #     # #     # # # #   
  # # # # #     # #     # # #     # #     # # # # # 
    # # # #     # #     # # #     # #     # # # #   
      # # # # # # # # # # # # # # # # # # # # #     
          # # # # # # # # # # # # # # # # #         
              # # # # # # # # # # # # #             
                  #       #       #                 
                # #       #       # #               
                #         #         #               
              # #         #         # #             
              #           #           #             
            # # #       # # #       # # #           
----------
This output was generated with the (very new) option {noresult} in mzn_show.pl. (This option don't display the real numeric results.)

September 30, 2009

This weeks news

This is mostly a consolidation (with some additional information) of the tweets on twitter.com/hakankj since last time.

mzn_show.pl

The Perl program mzn_show.pl is used for making little better output from the FlatZinc solvers.

The parameters to the program are now:
  • {tr:from:to:}: translate digit to another digit/character
  • {trvar:var:from:to:}: translate digit to another digit/character for a specific variable var
  • {trtr:from_string:replacement_string:}: translates all digits in from_string to the corresponding digit/character in replacement_string (inspired from the Unix tr command)
  • {trtrvar:var:from_string:replacement_string:}: as trtr but for a specific variable var
  • {nospaces}: don't show space after character
  • {noprintf}: don't try to be clever with printf stuff
Example: for showing a nicer picture of a Nonogram:

flatzinc nonogram_create_automaton2.fzn | mzn_show.pl "{tr:1: :} {tr:2:#:} {nospace}"

Where the {tr:from:to:} means translate digit from to character to, and {nospace} means that no spaces are shown after each digit/character.

This is now probably better written with trtr:

flatzinc nonogram_create_automaton2.fzn | mzn_show.pl "{trtr:12: #:} {nospace}"

Also, see below for some more examples.

MiniZinc/Comet: Nonogram model fixed for "zero length" clues

I discovered that my MiniZinc and Comet models had a little bug in the generation of the automaton for "zero length" clues (i.e. the row/column is completely empty, just 0's).

The fixed versions is here:
* MiniZinc: nonogram_create_automaton2.mzn
* Comet: nonogram_regular.co

It was quite tricky to debug the generated automata, but I found a fix for it surprisingly soon. I won't publish the fixes here, but it suffice to say that the MiniZinc model is - shall we say - not more beautiful than before.

A test instance for this: T2 (source unknown)
* Comet: nonogram_t2.co. * MiniZinc: nonogram_t2.dzn

T2 has the following two solutions:
   ##   ##  ###
  #  # #  # #  #
  #      #  #  #
  #     #   ###
  #  # #  # #
   ##   ##  #

----------

   ##   ##  ###
  #  # #  # #  #
  #     #   #  #
  #      #  ###
  #  # #  # #
   ##   ##  #

One more Nonogram example

Besides the T2 problem instance mentioned above, there is a new Nonogram example: Gondola for both Comet and MiniZinc:
* Comet: nonogram_gondola.co
* MiniZinc: nonogram_gondola.dzn.
For the curious, here is the solution, using the following command (see above for the mzn_show.pl part):
$ mzn2fzn -G jacop --data nonogram_gondola.dzn nonogram_create_automaton2.mzn
$ java JaCoP.fz.Fz2jacop -a -s  nonogram_create_automaton2.fzn | mzn_show.pl "{trtr:12: #:} {nospaces}"
  #####      ######
  ######     #    #        #
     ###  ###########     ###
  ######   # #    #       # #
  #######  # # # ##   #   ###
     ####    #    #  ##   # ####
  #######    #  # # ##    ### #
  #######     #  # ###    # #
     ####     # #  # # #########
     ####   ######## # #
     ####  #     ####  #   ###
     #### # #######    #  #####
     #### # #  ## # ####  #    #
     #### #########      ## # ##
     #### # ###   #      ##    #
     #### # ######       #   # #
     #### ########      ###    #
     ########## ###    ##### ###
      #### # ## ###    #####  ##
      ### #####  ##     ########
      ## ######  ###      #    #
      # ############     #
  ####################   # #
    ## #########################
   ##   ### ####################
  ##     ##### ### ## ## ## ## #
  #   ##  ######################
            ####################
    ###       ##################
  #        ##

Updated Comet's SONET model

In About 20 more constraint programming models in Comet Pierre Schaus (on of Comet's developers) commented that my use of tryall outside the using block in the SONET problem (sonet_problem.co) was not to recommend.

I have updated the model with his suggestion that uses or instead:
m.post(or(ring in 1..r) (rings[ring,client1] + rings[ring, client2] >= 2));
Thanks for the correction (suggestion), Pierre.

New MiniZinc model: Letter Square problem

letter_square.mzn is (yet another) grid problem. From my comment in the model file:

This problem is from the swedish book Paul Vaderlind: Vaderlinds nya hjärngympa ('Vaderlind's new brain gymnastics'), page 63ff. Unfortunately, I don't know the origin of this problem.

The objective is to create a matrix where all values on each row/column is different, except for the blanks.

A set of hints is given: the (or some) first seen non blank of each
- upper row
- lower row
- left column
- right column
seen from that view.

Example:
    B B
  -------
 |A   B C| 
B|  B C A| 
B|B C A  | 
 |C A   B| 
 --------
  C
This model codes the hints as follows:
 blank -> 0
 A     -> 1
 B     -> 2
 C     -> 3

row_upper = [0,2,2,0];
row_lower = [3,0,0,0];
col_left  = [0,2,0,0];
col_right = [0,0,0,0];
Note: there are no hints for the right column.

Here is yet another opportunity to use mzn_show.pl (0 is here translated to blank, etc):

solver model.fzn | mzn_show.pl "{trtr:0123456: ABCDEF:}"

With the following output:
  A B   C
      A B C
  B C   A
    A C   B
  C   B   A
The following problem instances from the book has been created:

September 23, 2009

MiniZinc Challenge 2009 Results

The result of MiniZinc Challenge 2009 is presented in MiniZinc Challenge 2009 Results:
There were two entrants this year:

* Gecode
* SICStus

In addition, the challenge organisers entered the following three FlatZinc implementations:

* ECLiPSe
* G12/FD
* G12/LazyFD

As per the challenge rules, these entries are not eligible for prizes, but do modify the scoring results.

Summary of Results

fd_search

sicstus 1651.8
eclipse_ic 322.1
gecode 4008.8
g12_fd 2040.6
g12_lazyfd 1376.6

free_search

sicstus 1841.0
gecode 4535.5
g12_fd 1112.4
g12_lazyfd 2511.1


Congratulations to the Gecode team!

September 21, 2009

A few new MiniZinc models, and a lot of improved

Some news about my MiniZinc models.

New MiniZinc models

This last weeks I have implemented the following new MiniZinc models:

Corrected some models

When testing the MiniZinc/FlatZinc support for the new version of JaCoP , I found problems in some models. These are now corrected:
  • strech_path.mzn: The former implementation was not correct.
  • min_index.mzn and max_index.mzn:
    The minimum(x[i], x) and maximum(x[i], x) don't work with the current MiniZinc ROTD version. Substituted to x[i] = min(x) and x[i] = max(x).

Improved all global constraints models

The global constraints section of My MiniZinc Page contains about 160 decompositions of global constraints from Global Constraint Catalog (and some not in the Catalog). The following improvements has been done on all models, especially for the older models:
  • Corrected the links to Global Constraint Catalog in the presentation of the constraint (only older models)
  • Removed some strange characters in the quoted text from Global Constraint Catalog (I hope all these has been removed now).
  • Made older models more general by using index_set, ub, lb, instead of assuming that all arrays start with index 1 etc. Some examples of this generality
           ...
           let {
             int: lbx = min(index_set(x)),
             int: ubx = max(index_set(x))
           } in
             forall(i in lbx+1..ubx) (
               forall(j in i+1..ubx-1) (
                  % ...
               )
             )
    
            forall(i in index_set(x)) (
              all_different([x[i,j] | j in index_set(x)) 
            )
          
            
            

September 16, 2009

JaCoP version 2.4 released

From JaCoP's news JaCoP version 2.4 is released:
Dear all,

We are happy to announce the release of a new version of our
Java-based solver JaCoP. This new version 2.4 has a number of new
features in addition to some bug fixes. The most important additions
in this version are:

1. The flatzinc interface that makes it possible to run minizinc
programs using JaCoP. The distribution contains number of different
minizinc examples.

2. Geometrical constraint, geost, based on pruning algorithms
originally proposed by Nicolas Beldiceanu et al. This constraint makes
it possible to define placement problems of non-convex objects in
k-dimensional space.

3. Knapsack constraint, which is based on the work published by Irit
Katriel et al. We extend the original work in number of ways, for example by
making it possible to use non-boolean quantity variables.

4. Set constraints defining typical operation on sets using set
interval variables.

This work would not be possible without help of several people. We
would like to thank Hakan Kjellerstrand for his help in testing
flatzinc parser as well as providing a number of examples in minizinc
format. We would also like to thank Meinolf Sellmann for his comments
on the initial implementation of knapsack constraint which have helped to
improve it further. Marc-Olivier Fleury has implemented most of the
functionality of the geost constraint. Robert Åkemalm has implemented the
first version of set constraint package. Wadeck Follonier has implemented
the first version of the Knapsack constraint. We would like to thank them
for their contributions.

As always feel free to send us ( radoslaw [dot] szymanek [at] gmail [dot] com )
feedback. We are always looking for cooperation to improve JaCoP. If you miss
some functionality of JaCoP we can help you to develop it so it can be done
efficiently and fast.


best regards,
Radoslaw Szymanek and Kris Kuchcinski
The latest version of JaCoP can be downloaded here.

For more information, see:
Also see My JaCoP page.

MiniZinc/FlatZinc support

I have especially tested the FlatZinc solver (Fz2jacop) and it is fast. For example, here is the statistics for nonogram_create_automaton2.mzn with the P200 problem (I wrote about this some days ago here):
Model variables : 629
Model constraints : 50

Search CPU time : 640ms
Search nodes : 1040
Search decisions : 520
Wrong search decisions : 520
Search backtracks : 520
Max search depth : 22
Number solutions : 1

Total CPU time : 1010ms
Note: JaCoP uses a special optimized regular constraint.

The FlatZinc solver has the following options:
$ java JaCoP.fz.Fz2jacop --help
Usage: java Fz2jacop [] .fzn
Options:
    -h, --help
        Print this message.
    -a, --all-solutions
    -t , --time-out 
         - time in second.
    -s, --statistics
    -n , --num-solutions 
         - limit on solution number.
Great work!

September 13, 2009

At last 2: A Nonogram solver using regular written in "all MiniZinc"

The model: nonogram_create_automaton2.mzn.

nonogram_regular.mzn

In At last, a Nonogram solver using regular constraint in MiniZinc I wrote about a Nonogram solver in MiniZinc using the regular solver, nonogram_regular.mzn. The drawback of this version is that an external program (make_nonogram_automata.pl) was needed to convert the Nonogram patterns (clues) to an automaton (finite states) for use in the regular constraint.

nonogram_create_automaton.mzn

In an update some hours later in the same blog post, the variant nonogram_create_automaton.mzn was mentioned. This is an "all MiniZinc" solution where the model also calculating the automata. The drawback of this was that the states is var int (decision variables), so it couldn't use the optimized regular constraint that some solvers (e.g. Gecode/FlatZinc) has implemented (this optimized version of regular is in fact the reason that Gecode/FlatZinc is so fast for the nonogram_regular.mzn model.) Instead a tweaked version of MiniZinc's standard regular constrant was used.

nonogram_create_automaton2.mzn

After a short discussion about this with Mikael Lagerkvist (of the Gecode team) the other day, we agreed that it would be nice thing to have the states calculated as "par" variables (i.e. not decision variables), thus the optimized regular could be used.

So I went got back to the drawing board and wondered how to do this. The solution I've found is not pretty, in fact it is hairy, very hairy: nonogram_create_automaton.mzn.

It is now as fast as nonogram_regular.mzn for solvers that use optimized version of regular, especially on the P200 problem.

Problem instances

Here are the problem instances I have tested (they are the same as used before.)

Comparison of calculating of the states

Now, let's compare the two different versions of calculating the automaton states.

First method: using decision variables

This is the version that is used in nonogram_create_automaton.mzn. The states are represented by the 2-dimensional array states.
states[1,1] = 1
/\ 
states[1,2] = 2
/\
forall(i in 2..len-1) (
   if i in zero_positions then
      states[i,1] = i+1 /\
      states[i,2] = 0 /\
      states[i+1,1] = i+1 /\
      states[i+1,2] = i+2
   else 
      if not(i - 1 in zero_positions) then
        states[i,1] = 0 /\
        states[i,2] = i+1
      else 
        true
      endif
   endif
)
/\
states[len,1] = len
/\
states[len,2] = 0
Quite neat, and relatively easy to understand.

Second method: no decision variables

And this is the non-decision variable version of calculating the finite states used in nonogram_create_automaton2.mzn. Note that states are represented by a 1-dimensional array, which is - as I have understood it - a requirement for this kind of initialisation of an array. It is also the cause of this hairyness.
array[1..2*len] of 0..len*2: states = 
[1, 2] ++
[
   if i div 2 in zero_positions then
       if i mod 2 = 0 then
        0
       else
        (i div 2) + 1
       endif
   elseif (i-1) div 2 in zero_positions then
       if i mod 2 = 0 then
        (i div 2)+1
       else
        (i div 2)+2
       endif
   else
     if not( (((i-1) div 2) - 1) in zero_positions) then
        if i mod 2 = 0 then
           (i div 2) + 1
        else 
          if (i div 2) + 1 in zero_positions then
              (i div 2) + 2
          else 
              0
          endif
        endif
      else
         if i mod 2 = 0 then
             (i div 2) + 1
         else 
            if not((i div 2) + 1 in zero_positions) then
               0
            else 
               (i div 2) + 2 
            endif
         endif
      endif
   endif
| i in 3..2*(len-1)]
++
[len, 0]
It could most probably be done in a more elegant fashion. And maybe I will think more about this later on.

September 09, 2009

At last, a Nonogram solver using regular constraint in MiniZinc

Here it is: nonogram_regular.mzn, a MiniZinc solver for Nonogram problems, using regular constraint.

In Comet version 2.0 released I wrote about the rewritten automata handler for the new Comet model nonogram_automaton.co (using the built-in constraint regular and helper function Automaton, both new in Comet version 2.0). This inspired me to finish the project of a "real" Nonogram solver for MiniZinc which use the regular constraint instead of the old (very slow) model nonogram.mzn.

Since MiniZinc has a built-in regular constraint, the hardest part was to create an automaton (finite state machine) given a Nonogram pattern. To be honest, I didn't write it purely in MiniZinc; instead a Perl program was written for this conversion (see below for more information about this program). Update: Well, I do have a version fully written in MiniZinc, nonogram_create_automaton.mzn, but it is not fast enough to be really interesting (must faster than the old version, nonogram.mzn, though). End of update this time.

The conversion pattern -> automaton is the same as described in Comet: regular constraint, a much faster Nonogram with the regular constraint, some OPL models, and more. To quote verbatim:
For the Nonogram clue [3,2,1] - which represents the regular expression 
"0*1110*110*10*" - the following automata (transition matrix) is 
generated:
1 2
0 3
0 4
5 0
5 6
0 7
8 0
8 9
9 0

Note that the regular function uses 0 (zero) as the failing state, so the states 
must start with 1..

Results: P200

In fact, the only model I was really curious about was the P200 problem since it has been the one that was the challenge for the Comet Nonogram solver. See the following posts for more about this quest: First, here is a picture of the solution of the P200 problem (nonogram_p200_aut.dzn, generated by the minizinc solver (currently the only solver that use output as formatting output):
 ##   ##     ###
####  # #    # ####
 #### # ##   #     #
####  # # #  #  #  #
 ##   # # ##  ###  #####
      #   # #   #  ##  #
     ###  # #####  #  ##
   ### ##  ##    # ##  ##
  ## #  ####   # # #    #
 ##      ## # ## #     ##
 #  #     # ### ##   ###
 #  # ##  #######  ###
 # ##     ## # #####
 ###  ## ##  #  ##
  ###   ##   #   ##
    #####    #    ##
  ##   ##    #     ##
 ####   ##   #    ##
######   ##  ### ##
#######   #### ###  ##
 #######    ####   ####
#######      #      ####
######       #     ####
 ####       ##      ##
  ##         #

P200: How does different solver do?

Here is a small benchmark for solving the P200 problem with different solvers. It was run on a Linux machine (Mandriva), Dual CPU 3.40GHz, 2Gb memory. "Runtime" is the total runtime, including converting from MiniZinc model to FlatZinc, startup, and also small time for running a wrapper program (where I can choose solver etc). Also, the solvers searched through the whole search tree for a solution (which happens to be exactly one).

All versions of the solvers was "the latest" as of 2009-09-08, i.e. the versions from respective CVS, SVN, or release of the day.

The search strategy used was first_fail where the columns (j) is labeled before rows (i):
solve :: int_search(
        [x[i,j] | j in 1..col_max, i in 1..row_max],
        first_fail,
        indomain_min,
        complete
        ) 
    satisfy;
For the MiniZinc/lazy solver, I also tested with solve satisfy since it often do very well without any specific search heuristics.
  • Gecode/FlatZinc: runtime 1.0s, solvetime 0.2s, failures 520 (also see below)
  • MiniZinc/minizinc: runtime 9s, 7697 choice points
  • MiniZinc/lazy: runtime 6s, choice points 2572
  • MiniZinc/lazy (with solve satisfy): runtime 13s, choice points 0
  • MiniZinc/fd: runtime 14s, choice points 11615
  • MiniZinc/fdmip: runtime 14s, choice points 11615
  • ECLiPSe/ic: runtime 109s
  • ECLiPSe/fd: runtime 34s
The Gecode/FlatZinc solver was by far the fastest. Here is the full statistics (for one random instance):
runtime:       0.207 (207.693000 ms)
solvetime:     0.198 (198.988000 ms)
solutions:     1
variables:     625
propagators:   50
propagations:  22940
nodes:         1041
failures:      520
peak depth:    22
peak memory:   1220 KB
mzx.pl nonogram_regular.mzn --fz --data nonogram_p200_aut.dzn  0,96s user 0,10s system 98% cpu 1,078 total
For comparison, here is the statistics for running the Comet model nonogram_automaton.co:
time:      459
#choices = 520
#fail    = 794
#propag  = 693993
comet nonogram_automaton.co  1,32s user 0,09s system 90% cpu 1,558 total
And, finally, the statistics for the Gecode (C++) program nonogram solving the P200 problem:
$ time nonogram -solutions 0 9
# ....
Initial
        propagators:  50
        branchings:   25

Summary
        runtime:      1.342 (1342.926000 ms)
        solutions:    1
        propagations: 35728
        nodes:        1409
        failures:     704
        peak depth:   25
        peak memory:  1027 KB
nonogram -solutions 0 9  1,45s user 0,06s system 99% cpu 1,517 total
It seems that the Gecode/FlatZinc version compares quite good.

For completion, here is the time for generating the automata of the P200 problem using the Perl program:
$ time perl make_nonogram_automata.pl nonogram_p200.dzn 1 > nonogram_p200_aut.dzn
0,20s user 0,02s system 96% cpu 0,233 total

Model and problem instances

Below are the problem instances as MiniZinc data file (.dzn), all the Nonogram problems listed on My Comet page, and some more. For each instance there are two variants: the "normal" version (which is the input file of the transformation) called "nonogram_name.dzn", and the generated automata version, called "nonogram_name_aut.mzn". It is the latter version ("_aut") that is used with nonogram_regular.mzn.

The program for converting to automata

As noted above, the Perl program make_nonogram_automata.pl converts a Nonogram problem instannce in "normal" pattern (clue) format to a format using automata.

The requirement of the indata file is semi-strict: the names must be the one in the example below, and each patterns must be on a separate line. The program use regular expressions to extract the information and can handle some variants in the format. The result is printed to standard output.
%% This is the problem instance of 'Hen', in "normal" pattern format.
%% Comments is keeped as they are
rows = 9;
row_rule_len = 2;
row_rules = array2d(1..rows, 1..row_rule_len,
       [0,3, 
        2,1, 
        3,2, 
        2,2, 
        0,6, 
        1,5, 
        0,6, 
        0,1, 
        0,2]);

   cols = 8;
   col_rule_len = 2;
   col_rules = array2d(1..cols, 1..col_rule_len,
       [1,2,
        3,1, 
        1,5, 
        7,1, 
        0,5, 
        0,3, 
        0,4, 
        0,3]);
To run the program:
$ perl make_nonogram_automata.pl nonogram_hen.dzn 1 > nonogram_hen_aut.dzn
The second parameter parameter 1 gives some more debugging info. The result is nonogram_hen_aut.dzn.

Minor note: For easy of debugging and further developments I decided to keep the 1-based version of arrays from the Comet model (Perl is 0-based) which made the code somewhat uglier.

September 04, 2009

The MiniZinc Wiki

The MiniZinc Wiki in brand new and includes (right now) for example: I have great hopes for the MiniZinc tutorial ("a work in progress") which now includes some commented models. Hopefully this will expand into a full tutorial of MiniZinc.

Worth of notice is also the SVN repository of MiniZinc examples: http://www.g12.cs.mu.oz.au/mzn. From the README file:
This subversion repository can be read and written by anyone in the constraint programming community who wishes to contribute MiniZinc models to the public domain.

August 11, 2009

Strimko - Latin squares puzzle with "streams"

Via the post A New Twist on Latin Squares from the interesting math blog 360, I learned the other day about a new grid puzzle Strimko, based on Latin squares (which for example Sudoku also is).

The rules of Strimko are quite simple:
Rule #1: Each row must contain different numbers.
Rule #2: Each column must contain different numbers.
Rule #3: Each stream must contain different numbers.
The stream is a third dimension of the puzzle: an connected "stream" of numbers which also must be distinct.

Here is an example of a Strimko puzzle (weekly #068, same as in the 360 post):
Strimko puzzle
(the link goes to a playable version).

MiniZinc model

This problem was simple to solve in MiniZinc: strimko.mzn, and the data file strimko_068.dzn.

Here are the constraints of the model: three calls to all_different, one for each rows, cols (these two are the latin square requirement), and also for each stream. The places is for handling the "hints" (the given numbers) in the puzzle.
constraint 
  % latin square
  forall(i in 1..n) (
      all_different([ x[i, j] | j in 1..n]) /\
      all_different([ x[j, i] | j in 1..n])
  )
  /\ % streams
  forall(i in 1..n) (
     all_different([x[streams[i,2*j+1],streams[i,2*j+2]] | j in 0..n-1])
  )
  /\ % placed
  forall(i in 1..num_placed) (
      x[placed[i,1], placed[i,2]] = placed[i,3]
  )
;
The data file strimko_068.dzn is shown below:. Note the representation of the puzzle: Each cell is represented with two numbers, row,column.
% Strimko Set 068
n = 4;
streams = array2d(1..n, 1..n*2, [
                    1,1, 2,2, 3,3, 4,4,
                    2,1, 1,2, 1,3, 2,4,
                    3,1, 4,2, 4,3, 3,4,
                    4,1, 3,2, 2,3, 1,4
                   ]);
num_placed = 3;
placed = array2d(1..num_placed, 1..3, [
                             2,2,3,
                             2,3,2,
                             3,3,1
                           ]);
Solution:
  4 2 3 1
  1 3 2 4
  2 4 1 3
  3 1 4 2

A better version

But it was quite boring and error prone to code the problem instance in that representation. There is - of course - a simpler way by represent the streams themselves, i.e. give each stream a unique "stream number", and all cells with the same numbers must be distinct. Data file: strimko2_068.dzn
% Strimko Weekly Set 068
streams = array2d(1..n, 1..n, [
                    1,2,2,4,
                    2,1,4,2,
                    3,4,1,3,
                    4,3,3,1
                  ]);
% ...
The model also needed a simple change to reflect this representation: strimko2.mzn:
  % ...
  /\ 
  % streams
  forall(s in 1..n) (
     all_different([x[i,j] | i,j in 1..n where streams[i,j] = s])
  )
  % ....
The main advantage of the second model is that it makes it easier to state the problem instances. Also it is somewhat smaller to represent the grid: n x n*2 versus n x n. However, I haven't noticed any difference between the performance of these two models, but the problems are all very small so it may not be visible.

The problem instances

Here are all the implemented models and the problem instances: For more information about Strimko, and some problems, see: For more about MiniZinc and some of my models, see My MiniZinc page.

August 02, 2009

MiniZinc: Some new implemented global constraints (decompositions)

For some reason I have implemented some more (decompositions of) global constraints from the great Global constraint catalog; all models use the stated example from that site.

I have tried as much as possible - and mostly succeeded - to make the predicates fully multi-directional, i.e. so that each and every parameter of a predicate can either be fix or free (decision variable).

See my other implementations of global constraint in MiniZinc here. Right now there are now about 150 of them, about of half of the 313 listed in the Global constraint catalog. Some other new MiniZinc models:

MiniZinc: the lazy clause generation solver

In the following two papers (for CP2009), a new Zinc/MiniZinc solver using both finite domain and SAT is discussed: the lazy clause generation solver:

  • T. Feydy and P.J. Stuckey: Lazy clause generation reengineered:
    Abstract Lazy clause generation is a powerful hybrid approach to com-
    binatorial optimization that combines features from SAT solving and fi-
    nite domain (FD) propagation. In lazy clause generation finite domain
    propagators are considered as clause generators that create a SAT de-
    scription of their behaviour for a SAT solver. The ability of the SAT
    solver to explain and record failure and perform conflict directed back-
    jumping are then applicable to FD problems. The original implemen-
    tation of lazy clause generation was constructed as a cut down finite
    domain propagation engine inside a SAT solver. In this paper we show
    how to engineer a lazy clause generation solver by embedding a SAT
    solver inside an FD solver. The resulting solver is flexible, efficient and
    easy to use. We give experiments illustrating the effect of different design
    choices in engineering the solver.

  • A. Schutt, T. Feydy, P.J. Stuckey, and M. Wallace: Why cumulative decomposition is not as bad as it sounds.
    AbstractThe global cumulative constraint was proposed for mod-
    elling cumulative resources in scheduling problems for finite domain (FD)
    propagation. Since that time a great deal of research has investigated new
    stronger and faster filtering techniques for cumulative, but still most of
    these techniques only pay off in limited cases or are not scalable. Re-
    cently, the “lazy clause generation” hybrid solving approach has been
    devised which allows a finite domain propagation engine possible to take
    advantage of advanced SAT technology, by “lazily” creating a SAT model
    of an FD problem as computation progresses. This allows the solver to
    make use of SAT nogood learning and autonomous search capabilities.
    In this paper we show that using lazy clause generation where we model
    cumulative constraint by decomposition gives a very competitive im-
    plementation of cumulative resource problems. We are able to close a
    number of open problems from the well-established PSPlib benchmark
    library of resource-constrained project scheduling problems.

This solver (the lazy solver) has been included in the MiniZinc distribution since version 1.0, and has been improved each version, e.g. by adding primitives it supports. See the NEWS file for more information.


I have tested the lazy solver with some problems, for example nonogram.mzn (which is implemented inefficient without the regular constraint), and am very impressed by it. It is very fast, though it cannot solve the P200 problem in a reasonable time.

Some comments about the lazy solver (some limitation is hopefully lifted in the future):

  • it cannot handle set vars.
  • all decision variables must be explicit bounded, e.g. var int: x; will not work, instead use something like var 1..10: x;.
  • Labeling: The papers state that the default labeling (e.g. solve satisfy) is often very good.

Conclusion: The lazy solver is definitely worth more testing.

July 08, 2009

MiniZinc version 1.0.1 released

From the MiniZinc mailing list:

Version 1.0.1 of the G12 MiniZinc distribution has been released.

Version 1.0.1 of the distribution can be downloaded from Download.
Snapshots of the development version of the G12 MiniZinc distribution
are also available from this page.

Further information about MiniZinc and FlatZinc is available from MiniZinc and FlatZinc

Bugs may be reported via the G12 bug tracking system, which can be accessed at http://bugs.g12.csse.unimelb.edu.au.


From the NEWS file:

G12 MiniZinc Distribution version 1.0.1

---------------------------------------

There have been no changes to the definitions of the MiniZinc and FlatZinc
languages in this release. There have been no changes to the definition
of XML-FlatZinc.


Changes in this release:

* MiniZinc tools command line changes

The MiniZinc interpreter and MiniZinc-to-FlatZinc converter now
recognise files with the .dzn extension as MiniZinc data files, i.e.
you no longer need to use the --data option to use such files as data
files.

* FlatZinc output processing tool

We have added a new tool, solns2dzn, that can be used to process the
output of FlatZinc implementations in various ways. For example, it
can extract each individual solution and write it to a separate file.

* The FlatZinc interpreter's lazy clause generation solver now supports
the int_times/3 built-in.

* The global_cardinality_low_up global constraint has been added to the
MiniZinc library.

* The MiniZinc-to-FlatZinc converter now propagates annotations through
assertions during flattening, For example, the following fragment of
MiniZinc:

predicate foo(int: x, array[int] of var int y);

predicate bar(int: x, array[int] of var int y) =
assert(x > 3, "value of x must be greater than 3", foo(x, y));

constraint bar(4, ys) :: baz;

will be flattened into the following fragment of FlatZinc:

constraint foo(4, ys) :: baz;


Bugs fixed in this release:

* A bug in the implementation of the FlatZinc built-in int_mod/3 in the
FlatZinc interpreter's finite-domain backend has been fixed.

* The MiniZinc-to-FlatZinc converter now does a better job of extracting
output variables from output items.

* The solver-specific global constraint definitions for the G12 lazy
clause generation solver are now documented. They are in the
``g12_lazyfd'' directory of the MiniZinc library.

* A bug that caused a segmentation fault in the type checker when
checking large FlatZinc instances has been fixed.

* A bug where a cycle of equalities caused mzn2fzn to go into
a loop has been fixed. [Bug #65]

* mzn2fzn now imposes constraints on arguments induced by predicate
parameter types. [Bug #69]

* Flattening of set2array coercions is now supported. [Bug #70]

* A bug that caused mzn2fzn to abort on min/max expressions over empty arrays
has been fixed. [Bug #71]

* mzn2fzn now computes bounds on set cardinality variables correctly.

June 24, 2009

My old swedish blog posts about constraint programming translated (via Google)

Before starting this blog (i.e. My Constraint Programming blog) late December 2008, I blogged about constraint programming in my Swedish blog hakank.blogg). Translations of those Swedish blog post has not been collected before, and now it is time.

So, here are links to most of the blog posts from the category Constraint Programming, translated via Google's Language Tools. Most of the translation is intelligible, but if you have some questions about some particular, please feel to mail me: hakank@bonetmail.com for more information.


November 4, 2005
Swedish: Sesemans matematiska klosterproblem samt lite Constraint Logic Programming
Translated: Sesemans kloster mathematical problems and little Constraint Logic Programming ("kloster" means "convent").


April 18, 2006
Swedish: Choco: Constraint Programming i Java
Translated: Choco: Constraint Programming in Java


The two post above was written when I had just a cursory interest in constraint programming. From about February 2008 and onward, it became my main interest.


April 5, 2008
Swedish: Constraint Programming: Minizinc, Gecode/flatzinc och ECLiPSe/minizinc,
Translated: Constraint Programming: Minizinc, Gecode / flatzinc and ECLIPSE / minizinc.


April 14, 2008
Swedish: MiniZinc-sida samt fler MiniZinc-exempel
Translated: MiniZinc-page and multi-MiniZinc example.


April 21, 2008
Swedish: Applikationer med constraint programming, lite om operations research samt nya MiniZinc-modeller
Translated: Applications of constraint programming, a little of operations research and new MiniZinc models


April 26, 2008
Swedish: Ett litet april-pyssel
Translated: A small-craft in April (a better translation would be "A small April puzzle")


April 27, 2008
Swedish: Mitt OpenOffice Calc-/Excel-skov
Translated: My Open Office Calc-/Excel-skov (better translation: My Open Office Calc-/Excel period).


May 26, 2008
Swedish: Några fler MiniZinc-modeller, t.ex. Smullyans Knights and Knaves-problem
Translated: Some more MiniZinc models, eg Smullyans Knights and Knaves-problem Smullyans Knights and Knaves problem


June 2, 2008
Swedish: MiniZinc/FlatZinc version 0.8 släppt
Translated: MiniZinc / FlatZinc version 0.8 released


June 2, 2008
Swedish: Nya MiniZinc-modeller, flera global constraints, däribland clique
Translated: New MiniZinc models, several global constraints, including Clique


June 5, 2008
Swedish: Mats Anderssons tävling kring fotbolls-EM 2008 - ett MiniZinc-genererat tips
Translated: Mats Andersson racing around championship in 2008 - a MiniZinc-generated tips (it is about a competition how to predict Soccer World Cup 2008 using MiniZinc)


June 24, 2008
Swedish: Tre matematiska / logiska pyssel med constraint programming-lösningar: n-puzzle, SETS, M12 (i MiniZinc)
Translated: Three mathematical / logical craft with constraint programming solutions: n-puzzle, SETS, M12 (in MiniZinc) ("craft" should be translated to "puzzles")


June 24, 2008
Swedish: MiniZinc-nyheter
Translated: MiniZinc news


June 29, 2008
Swedish: Gruppteoretisk lösning av M12 puzzle i GAP
Translated: Group Theoretical solution of the M12 puzzle in GAP (well, this is not really a constraint programming solution, but it is another way of solving the M12 puzzle blogged about in June 24)


June 30, 2008
Swedish: Gruppteoretisk lösning av M12 puzzle i GAP - take 2
Translated: Group Theoretical solution of the M12 puzzle in GAP - take 2


July 4, 2008
Swedish: Martin Chlond's Integer Programming Puzzles i MiniZinc
Translated: Martin's Chlond Integer Programming Puzzles in MiniZinc


July 7, 2008
Swedish: Fler MiniZinc modeller kring recreational mathematics
Translated: More MiniZinc models on recreational mathematics


July 20, 2008
Swedish: Fler constraint programming-modeller i MiniZinc, t.ex. Minesweeper och Game of Life
Translated: More constraint programming models in MiniZinc, eg Minesweeper och Game of Life Minesweeper and the Game of Life


August 17, 2008
Swedish: JaCoP - Java Constraint Programming solver
Translated: JaCoP - Java Constraint Programming solver


September 14, 2008
Swedish: Constraint programming i Java: Choco version 2 släppt - samt jämförelse med JaCoP
Translated: Constraint programming in Java: Choco version 2 released - and comparison with JaCoP


September 28, 2008
Swedish: Constraint programming: Fler MiniZinc-modeller, t.ex. Martin Gardner och nonogram
Translated: Constraint programming: More MiniZinc models, eg Martin Gardner och nonogram Martin Gardner and nonogram


December 27, 2008
Swedish: Constraint programming-nyheter samt nya MiniZinc-modeller
Translated: Constraint programming, news and new MiniZinc models


December 29, 2008
Swedish: My Constraint Programming Blog
Translated: My Constraint Programming Blog


After that, entries about constraint programming at my swedish blog posts where just summaries of the stuff written here at My Constraint Programming blog.

June 13, 2009

Miscellaneous news

Here are some miscellaneous news.

Choco

Version 2.1

I have forgot to mention that Choco have released version 2.1 . Download this version via the download page. The Sourceforge download page is here.

Changed models for version 2.1

I have changed some of my Choco models for version 2.1. The biggest change was the use of cumulative in FurnitureMoving.java where the call to the cumulative constraint have changed and now use TaskVariable.

I added the following
// Create TaskVariables
TaskVariable[] t = new TaskVariable[starts.length];
   for(int i = 0; i < starts.length; i++){
      t[i] = makeTaskVar("", starts[i], ends[i], durations[i]);
   }
and changed the call to cumulative to:
 
m.addConstraint(cumulative(null, t, constantArray(_heights), constant(NumPersons)));
Also, the line
durations[i] = makeConstantVar("durations" + i, durationsInts[i]);
was changed to
durations[i] = makeIntVar("durations" + i, durationsInts[i], durationsInts[i]);
For some of the other models (compiled for version 2.0.*), I just recompiled the source and it worked without any change in the code.

MiniZinc

List of global constraints

At the MiniZinc site, there is now a list of the supported global constraints in MiniZinc: MiniZinc: Global Constraints.

MiniZinc Library

MiniZinc Library is a great collections of problems, examples, and tests. (For some reason is not linked from the main MiniZinc site.)

Output in Minizinc

I wrote here about the change of output in MiniZinc version 1.0 . Only the distributed solver minizinc supports the output statement (for some kind of "pretty printing" of the output), but the solvers using the FlatZinc format has no such support. Since I want to see matrices presented in a normal way for all solvers, I wrote a simple Perl program to show matrices more pretty than just a single line. Also, single arrays are shown without array1d(...).

The Perl program is mzn_show.pl and it simply takes the result from a solver and reformats the result.

An example: The output from the debruijn_binary.mzn from a solver using FlatZinc file, such as Gecode/FlatZinc and MiniZinc's fdmip is:
bin_code = array1d(1..8, [0, 0, 0, 0, 1, 0, 1, 1]);
binary = array2d(1..8, 1..4, [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0]);
x = array1d(1..8, [0, 1, 2, 5, 11, 6, 12, 8]);
When filtering this with mzn_show.pl it will be shown as:
bin_code: 0 0 0 0 1 0 1 1
% binary = array2d(1..8, 1..4, [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0]);
binary:
  0 0 0 0
  0 0 0 1
  0 0 1 0
  0 1 0 1
  1 0 1 1
  0 1 1 0
  1 1 0 0
  1 0 0 0
x: 0 1 2 5 11 6 12 8
Very simple, but quite useful.

Somewhat related: The NEWS file of the latest ROTD (Release of The Day) states the following:
FlatZinc output processing tool

We have added a new tool, solns2dzn, that can be used to process the output of FlatZinc implementations in various ways. For example, it can extract each individual solution and write it to a separate file.
This program is distributed as source (Mercury), not as an executable, so I have not been able to test it.

Hidato and exists

In the hidato.mzn model I have changed the exists construct
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]
     )
  )
)
to a a version using just var ints:
forall(k in 1..r*c-1) (
    let {
       var 1..r: i,
       var 1..c: j,
       var {-1,0,1}: a,
       var {-1,0,1}: b
    }
    in
    k = x[i, j] % fix this k
    /\
    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]
)
This replacement of exists with int vars, if possible, seems to always be more effective.

However, there is one use of exists where it is harder to replace in this way. As an example, take Pandigital numbers in any base (the model includes a presentation of the problem) where - among other things - we want to find some array indices of the integer array x in order to find the positions of three numbers num1, num2 and res (the result of num1 * num2).
% num1. 
exists(j in 1..x_len) (
   j = len1 /\
   toNum([x[i] | i in 1..j], num1, base)
)

/\  % num2
exists(j, k in 1..x_len) (
   j = len1+1 /\ 
   k = len1+len2 /\ k > j  /\
   toNum([x[i] | i in j..k], num2, base)
)

/\ % the product
exists(k in 1..x_len) (
   k = len1+len2+1 /\
   toNum([x[i] | i in k..x_len], res, base)
)
Using this approach, we have to use exists since indices in a range (e.g. j in 1..j) must not be an int var.

Also

Both Choco and MiniZinc sites now have links to my site, which is quite fun:
* Choco, Users (last)
* MiniZinc and FlatZinc (also last)

Work related page in english

The page Håkan Kjellerstrand, CV and work related interests is a english version of my work related page. (The original, swedish, version is here.)

May 29, 2009

MiniZinc Challenge 2009

From MiniZinc Challenge 2009:
The aim of the challenge is to start to compare various constraint solving technology on the same problems sets. The focus is on finite domain propagation solvers. An auxiliary aim is to build up a library of interesting problem models, which can be used to compare solvers and solving technologies.
Announcement of the results will be at CP2009: 20 September 2009.

Also see
* MiniZinc Challenge 2009 -- Rules
* MiniZinc Challenge 2008 and Results.

May 20, 2009

MiniZinc version 1.0 released!

MiniZinc version is now official, i.e. version 1.0. Congratulations to the G12 team!

This version can be downloaded from here.

Comparing to the beta version 0.9 (released Christmas last year) the following has changed:


G12 MiniZinc Distribution version 1.0
-------------------------------------

* Licence change

The source code in the G12 MiniZinc distribution has now been released
under a BSD-style licence. See the files README and COPYING in the
distribution for details.

The MiniZinc examples, global constraint definitions and libraries
have been placed in the public domain.

* XML-FlatZinc

We have defined an XML representation for FlatZinc called XML-FlatZinc.
Two new tools, fzn2xml and xml2fzn, can be used to convert between FlatZinc
and XML-FlatZinc.

* FlatZinc Conformance Test Suite

We have added a suite of conformance tests for FlatZinc implementations.
It includes tests for built-in constraints, output and the behaviour of the
standard search annotations.

Changes to the MiniZinc language:

* Reification has been fixed. It is now a top-down process that correctly
handles partial functions such as integer division. Users can now also
supply alternative definitions for reified forms of predicates (this is
useful if a backend does not provide reified forms of all predicates).

* Users can supply alternative definitions for FlatZinc built-in constraints
(e.g., one can force the generated FlatZinc to use just int_lt rather than
int_lt and int_gt).

* A new variable annotation has been added: is_output is used to indicate
variables to be printed as part of the solution if no output item is
supplied. This annotation is converted to output_var or output_array as
appropriate by mzn2fzn.

Changes to the FlatZinc language:

* The outdomain_min and outdomain_max value choice methods are now supported
in the finite-domain solver backend.

* A new search annotation, seq_search, allows a sequential ordering to
be imposed on search annotations.

* The standard solve annotations now use nested annotations instead of
strings to describe variable selection strategies, value choice methods,
and exploration strategies.

* FlatZinc model instances may now contain bodyless predicate declarations.
This is to allow tools to type check FlatZinc that contains non-standard
built-in predicates.

* Two new annotations have been added that allow functional relationships
between variables to be defined: is_view_var on a variable declaration
states that this var is defined as a function of some other variables by
a constraint; defines_view_var(x) on a constraint states that the
constraint provides a definition for the view variable x.

* The FlatZinc specification now specifies how multiple solutions should be
output.

* Two new variable annotations have been added to indicate which variables
should be printed if a solution is found: output_var is used
for non-array variables; output_array([IndexSet1, IndexSet2, ...]) is used
for array variables.

* Output items are no longer supported in FlatZinc. The built-in string
operations, show/1 and show_cond/3 have also been removed from the
language.

Changes to the G12 MiniZinc-to-FlatZinc converter:

* The converter now outputs array_xxx_element constraints instead of
array_var_xxx_element constraints when the array argument is a
constant.

* An error is now reported if a variable is defined in a let expression
in a negated or reified context.

* The ZINC_STD_SEARCH_DIRS environment variable is no longer supported.
The new environment variable MZN_STDLIB_DIR or the command line option
``--stdlib-dir'' can be used to set the MiniZinc library directory.

* String array lookups are now supported.

* Comparison of fixed string expressions is now supported.

* Bodyless predicates in MiniZinc are now emitted at the head of the
generated FlatZinc. For backwards compatibility this behaviour
can be disabled using the ``--no-output-pred-decls'' command line
option.


Changes to the G12 MiniZinc and FlatZinc interpreters:

* There is a new solver backend for the FlatZinc interpreter based upon
the G12 lazy clause generation solver. This backend is selected with
the ``lazy'' or ``lazyfd'' argument to the interpreter's ``--backend''
option.

* The implementation of the int_negate/2 builtin constraint has been
fixed.

* The interpreters now take a flag (--solver-statistics or --solver-stats)
causing any statistical information gathered by the solver to be appended to
the output in the form of a Zinc comment.

* The interpreters now take a flag (-a or --all-solutions) that causes
then to return all solutions.

* The interpreters now take a flag (-n or --num-solutions) taking an
integer argument giving the maximum number of solutions to display.

* The MiniZinc interpreter behaviour has changed for models with no
output item: now only the values of variables annotated with
is_output are printed.

Other Changes:

* The following global constraints have been added to the MiniZinc library:

at_most1
nvalue
precedence
table

For developers of FlatZinc solvers there is a transition guide to the changes.

Some comments

One of the greatest news is that it is now possible to state the number of solutions: -a for all solutions, or -n or --num-solutions for some fixed number of solutions. I have missed that feature for the distributed solvers since day one. Other solvers, e.g. Gecode/FlatZinc, ECLiPSe's, and, SICStus Prolog supported this feature already.

I'm also excited that string array lookups is supported. And there are other stuff which I will now test more in detail.

My MiniZinc page

On My MiniZinc page there are over 540 MiniZinc models which right now are for MiniZinc version 0.9. I have just started to convert these to the new version and it will take some time. Some models don't have to be changed at all, though.


Update some hours later
All models has now been converted to version 1.0. Here are some findings while doing this conversion:

* If there are any file in the current directory with the same name as in lib/minizinc/std/ (etc) then there will be a conflict with message structure error: more than one solve item

* solve: the search strategies should not have quotes. Also the search strategy occurrences is not available for the minizinc/flatzinc solvers.

* if there is no output clause, then define some (or all) decision variables as :: is_output and they will be printed.

* it seems that it is just the minizinc solver that handles the output section. For other solvers in the MiniZinc distribution (e.g. mip) the output is done via :: is_output annotation.

* always add "\n" last in the output clause.


May 09, 2009

Learning Constraint Programming IV: Logical constraints: Who killed Agatha? revisited

Here is a comparison of how different constraint programming systems implements the logical constraints in the problem Who Killed Agatha?, also known as The Dreadsbury Mansion Murder Mystery.. In Learning constraint programming - part II: Modeling with the Element constraint the problem was presented and showed how the systems implements the Element constraints.

Problem formulation from The h1 Tool Suite
Someone in Dreadsbury Mansion killed Aunt Agatha. Agatha, the butler, and Charles live in Dreadsbury Mansion, and are the only ones to live there. A killer always hates, and is no richer than his victim. Charles hates noone that Agatha hates. Agatha hates everybody except the butler. The butler hates everyone not richer than Aunt Agatha. The butler hates everyone whom Agatha hates. Noone hates everyone. Who killed Agatha?
Originally from F. J. Pelletier: Seventy-five problems for testing automatic theorem provers., Journal of Automated Reasoning, 2: 191-216, 1986.

Here we see compare how different constraint programming systems implement the three emphasized conditions in the problem formulation above:
  • the concept of richer
  • Charles hates noone that Agatha hates
  • No one hates everyone
All models defines the concepts of hates and richer as two matrices. The matrix declarations are omitted in the code snippets below.

Models

Here are the different models for the Who killed Agatha? problem. JaCoP and Choco has two version for how to implement the Element constraint, see the link above. Also, there is no Gecode/R version of this problem.

Defining the concept of richer

First we define the concept of richer, which consists of two parts:
  • No one is richer than him-/herself
  • if i is richer than j then j is not richer than i
This is an antisymmetric relation which is explored somewhat more (with an alternative predicate of the concept) in the MiniZinc model antisymmetric.mzn.

The logical concept used here is equivalence (if and only of).

MiniZinc

% define the concept of richer: no one is richer than him-/herself
forall(i in r) (
   richer[i,i] = 0
)

/\ % if i is richer than j then j is not richer than 
forall(i, j in r where i != j) (
   richer[i,j] = 1 <-> richer[j,i] = 0
)

Comet

Note that Comet don't have support for equivalence directly, instead we have to use two implications. Update: equivalence in Comet is written as == (I tested <=> which didn't work). Thanks to Pascal Van Hentenryck for pointing this out.
// no one is richer than him-/herself
forall(i in r)
  m.post(richer[i,i] == 0);

// if i is richer than j then j is not richer than i
forall(i in r, j in r : i != j) {
  /* earlier version: two implications
  m.post(richer[i,j] == 1 => richer[j,i] == 0);
  m.post(richer[j,i] == 0 => richer[i,j] == 1);
  */
  // equivalence
  m.post((richer[i,j] == 1) == (richer[j,i] == 0));
}

JaCoP

//  no one is richer than him-/herself
for(int i = 0; i < n; i++) {
    store.impose(new XeqC(richer[i][i], 0));
}

//  if i is richer than j then j is not richer than i
for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
        if (i != j) {
             store.impose(
                          new Eq(
                                 new XeqC(richer[i][j], 1),
                                 new XeqC(richer[j][i], 0)
                                 )
                          );

        }
    }
}

Choco

//   a) no one is richer than him-/herself
for(int i = 0; i < n; i++) {
    m.addConstraint(eq(richer[i][i], 0));
}

//   b) if i is richer than j then j is not richer than i
for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
        if (i != j) {
             m.addConstraint(
                             ifOnlyIf(
                                    eq(richer[i][j], 1),
                                    eq(richer[j][i], 0)
                                      )
                             );
        }
    }
}

Gecode

// no one is richer than him-/herself
for(int i = 0; i < n; i++) {
  rel(*this, richer_m(i,i), IRT_EQ, 0, opt.icl());
}

// if i is richer than j then j is not richer than i
for(int i = 0; i < n; i++) {
  for(int j = 0; j < n; j++) {
    if (i != j) {
      post(*this, tt(
              eqv(
                 richer_m(j,i) == 1, // <=>
                 richer_m(i,j) == 0)), 
                 opt.icl());
    }
  }
}

Charles hates noone that Agatha hates

Here is the definitions of the condition Charles hates noone that Agatha hates, which simply mean the implication:
  if Agatha hates X then Charles don't hate X

MiniZinc

When starting to modeling these kind of problems, I tend to follow the order of the conditions, which here means that the Charles part is before the Agatha part. When remodeling in another system the order tends to be fixed (cf the Comet version).
forall(i in r) (
   hates[charles, i] = 0 <- hates[agatha, i] = 1
)

Comet

forall(i in r)
  m.post(hates[agatha, i] == 1 => hates[charles, i] == 0);

JaCoP

for(int i = 0; i < n; i++) {
    store.impose(
                 new IfThen(
                            new XeqC(hates[agatha][i], 1),
                            new XeqC(hates[charles][i], 0)
                            )
                 );
}

Choco

I tend to copy/paste the models for Choco and JaCoP and just change the functions that are different. A consequence of this is that some special feature in one of these two systems is not used.
for(int i = 0; i < n; i++) {
    m.addConstraint(
                    implies(
                            eq(hates[agatha][i], 1),
                            eq(hates[charles][i], 0)
                            )
                    );
}

Gecode

for(int i = 0; i < n; i++) {
   post(*this, tt(
                  imp(
                        hates_m(i,agatha) == 1, 
                        // =>
                        hates_m(i,charles) == 0)), 
                        opt.icl());
}

No one hates everyone

This is the last condition to compare: No one hates everyone. It is implemented by a sum of the number of persons that each person hates, and this sum must be 2 or less. Note that it is possible to hate oneself.

MiniZinc

forall(i in r) (
  sum(j in r) (hates[i,j]) <= 2  
)

Comet

  forall(i in r) 
    m.post(sum(j in r) (hates[i,j]) <= 2);

JaCoP

Note: We could save the XlteqC constraint by restrict the domain of a_sum to 0..2 instead of 0..n (=3) but this explicit use of XlteqC is more clear.
for(int i = 0; i < n; i++) {
    FDV a[] = new FDV[n];
    for (int j = 0; j < n; j++) {
        a[j] = new FDV(store, "a"+j, 0, 1);
        a[j] = hates[i][j];
    }
    FDV a_sum = new FDV(store, "a_sum"+i, 0, n);
    store.impose(new Sum(a, a_sum));
    store.impose(new XlteqC(a_sum, 2));
}

Choco

Note: sum is an operator, which makes the condition somewhat easier to state than in JaCoP.
for(int i = 0; i < n; i++) {
    IntegerVariable a[] = makeIntVarArray("a", n, 0, 1);
    for (int j = 0; j < n; j++) {
        a[j] = hates[i][j];
    }
    m.addConstraint(leq(sum(a), 2));
}

Gecode

In Gecode this condition is quite easy to state by using linear. In order to use this there is a Matrix "view" of the hates matrix hates_m.
Matrix hates_m(hates, n, n);
// ...
for(int i = 0; i < n; i++) {
  linear(*this, hates_m.row(i), IRT_LQ, 2, opt.icl());
}

End comment

The mandatory end comment: There are probably better ways of modeling the problem than shown above, either by changing some details or by model the problem completely different. Maybe this will be done sometime...

May 08, 2009

Learning Constraint Programming III: decomposition of a global constraint: alldifferent_except_0

The series Learning Constraint Programming is a preparation for My talk at SweConsNet Workshop 2009: "Learning Constraint Programming (MiniZinc, JaCoP, Choco, Gecode/R, Comet, Gecode): Some Lessons Learned". Confusingly, the entries is not numbered in any logical order. Sorry about that.

Here are the previous entries:

Introduction

The global constraint alldifferent_except_0 (or the more general variant alldifferent_except_c) is one of my favorite global constraints. It is very handy to use when 0 (or any constant c) is coded as an unknown or irrelevant value. Then we can constraint the rest to be all distinct.

The great Global Constraint Catalog entry alldifferent_except_0) explains this constraint as:
Purpose
Enforce all variables of the collection VARIABLES to take distinct values, except those variables that are assigned to 0.

Example
​(<5,​0,​1,​9,​0,​3>)​
The alldifferent_except_0 constraint holds since all the values (that are different from 0) 5, 1, 9 and 3 are distinct.

Models

I have modeled a decomposition of alldifferent_except_0 in the following models, where the constraint is just tested perhaps combined with some other constraint, e.g. sorted or that there must be at least some zeros:

- MiniZinc alldifferent_except_0.mzn
- Comet: alldifferent_except_0.co
- Gecode/R: all_different_except_0.rb
- Choco: AllDifferentExcept0_test.java
- JaCoP: AllDifferentExcept0_test.java
- Gecode: alldifferent_except_0.cpp

Some models using alldifferent_except_0

And here is some real use of the constraint:

- Nonogram (Comet): nonogram.co. (A faster model using the regular constraint, nonogram_regular.co, is described here and here).
- I wrote about alldifferent_except_0 in Pi Day Sudoku 2009. However, as faster way of solving the problem was found and is described in Solving Pi Day Sudoku 2009 with the global cardinality constraint). Note: the competition is still on, so there is no link to any model.
- Sudoku generate (Comet): sudoku_generate.co
- all paths graph (MiniZinc): all_paths_graph.mzn
- Cube sum (MiniZinc): cube_sum.mzn
- Message sending (MiniZinc): message_sending.mzn

As the first two entries indicates, there may be faster solutions than using (a decomposition) of alldifferent_except_0, but even as a decomposition is a great conceptual tool when modeling a problem.

Implementations

In the implementations below we also see how to define a function (predicate) in the constraint programming systems.

For the Gecode/R model there are different approaches:
- "standard" ("direct") approach where we loop over all different pairs of elements and ensures that if both values are different from 0 then they should be different
- using count
- using global cardinality ("simulated" in Gecode/R, see below)

Also, in some models we use the slighly more general version alldifferent_except_c where c is any constant (e.g. "Pi" in the Pi Day Sudoku puzzle mentioned above.

MiniZinc:

Model: alldifferent_except_0.mzn.
predicate all_different_except_0(array[int] of var int: x) =
   let {
      int: n = length(x)
   }
   in
   forall(i,j in 1..n where i != j) (
        (x[i] > 0 /\ x[j] > 0) 
        -> 
        x[i] != x[j]
   )
;

// usage:
constraint all_different_except_0(x);

Comet:

Model: alldifferent_except_0.co.
function void alldifferent_except_0(Solver m, var{int}[] x) {
  int n = x.getSize();
  forall(i in 1..n, j in i+1..n) {
    m.post(
           x[i] > 0 && x[j] > 0 
           => 
           x[i] != x[j]
           );
  }
}

// usage
exploreall {
  // ...
  alldifferent_except_0(m, x);
}

Gecode/R

Model:all_different_except_0.rb.

When modeling the constraint in Gecode/R, I experimented with different approaches. The reification variant all_different_except_0_reif is actually quite fast.
# The simplest and the fastest implementation 
# using count for 1..max (poor man's global cardinality)
def all_different_except_0
  (1..self[0].domain.max).each{|i|
    self.count(i).must <= 1
  }
end

# global cardinality version using an extra array with the counts
def global_cardinality(xgcc)
  (self[0].domain.max+1).times{|i|
    xgcc[i].must == self.count(i)
  }
end

#
# The standard approach using reification.
#
def all_different_except_0_reif(x)
  n = x.length
  b1_is_an bool_var_matrix(n,n)
  b2_is_an bool_var_matrix(n,n)
  b3_is_an bool_var_matrix(n,n)
  n.times{|i|
    n.times{|j|
      if i != j then 
        x[i].must_not.equal(0, :reify => b1[i,j]) 
        x[i].must_not.equal(0, :reify => b2[i,j]) 
        x[i].must_not.equal(x[j], :reify => b3[i,j])
        (b1[i,j] & b2[i,j]).must.imply(b3[i,j])
      else 
        b1[i,j].must.true
        b2[i,j].must.true
        b3[i,j].must.true
      end
    }
  }
end

    # ...
    # usage:
    x.all_different_except_0
    # all_different_except_0_gcc(x)
    # all_different_except_0_reif(x)

Choco

Model: AllDifferentExcept0_test.java.

Note that here alldifferent_except_0 is derived from the more general version alldifferent_except_c.
//
// decomposition of alldifferent except 0
//
public void allDifferentExcept0(CPModel m, IntegerVariable[] v) {
    allDifferentExceptC(m, v, 0);
}

//
// slightly more general: alldifferent except c
//
public void allDifferentExceptC(CPModel m, IntegerVariable[] v, int c) {
    int len = v.length;
    for(int i = 0; i < v.length; i++) {
        for(int j = i+1; j < v.length; j++) {
            m.addConstraint(ifThenElse(
                                       and(
                                           gt(v[i], c), 
                                           gt(v[j], c)
                                           ), 
                                       neq(v[i], v[j]),
                                       TRUE
                                       )
                            );
        }
    }    
}

// ...
// usage: 

    m.addConstraint(allDifferent(queens));

JaCoP

Model: AllDifferentExcept0_test.java

This is exactly the same approach as the Choco version.
//
// decomposition of alldifferent except 0
//
public void allDifferentExcept0(FDstore m, FDV[] v) {
    allDifferentExceptC(m, v, 0);
} // end allDifferentExcept0


//
// slightly more general: alldifferent except c
//
public void allDifferentExceptC(FDstore m, FDV[] v, int c) {
    int len = v.length;
    for(int i = 0; i < v.length; i++) {
        for(int j = i+1; j < v.length; j++) {
            m.impose(new IfThen(
                                       new And(
                                           new XneqC(v[i], c), 
                                           new XneqC(v[j], c)
                                           ), 
                                       new XneqY(v[i], v[j])
                                       )
                            );
        }
    }        
} // end allDifferentExceptC


        // ...
        // usage:
        allDifferentExcept0(store, x);

Gecode

alldifferent_except_0.cpp

The Gecode version is very succint since it use overloaded boolean operators. Very nice.
void alldifferent_except_0(Space& space, IntVarArray x, 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] != 0 && x[j] != 0, 
           // =>
           x[i] != x[j])),
           icl
           );
    }
  }
} // alldifferent_except_0


// ...
// usage:
    alldifferent_except_0(*this, x, opt.icl());

May 05, 2009

Learning constraint programming - part II: Modeling with the Element constraint

As indicated last in Learning constraint programming (languages) - part I here are some findings when implementing Crossword, Word square, and Who killed Agatha?. See links below for the implementations.

Spoiled
The first constraint programming system i learned after constraint logic programming in Prolog was MiniZinc. When implemented the problems below I realized that I have been quite spoiled by using MiniZinc. The way MiniZinc (and also Comet) supports the Element constraint, i.e. access of variable arrays/matrices, is very straightforward in these systems and it doesn't matter whether the array to access is an array of integers or of non-variable variable integers. In the Java (Choco, JaCoP) and C++ systems (Gecode), however, this is another matter. Due to different circumstances I have not implemented these models in Gecode/R.

Element in MiniZinc and Comet
Accessing arrays and matrices in MiniZinc and Comet is simply done by using the [] construct, no matter what the type of the array or the index are (I assume integers and variable integers here). For the other systems we must explicitly use the Element constraint (called nth in Choco).

Crossword

This is a standard constraint programming problem, used as a running example in Apt's great book Principles of Constraint Programming. Here is a formulation of the problem (stated differently than in the book):
Place the words listed below in the following crossword. The '#' means a blocked cell, and the numbers indicate the overlappings of the words.

      1   2   3   4   5
    +---+---+---+---+---+
  1 | 1 |   | 2 |   | 3 |
    +---+---+---+---+---+
  2 | # | # |   | # |   |
    +---+---+---+---+---+
  3 | # | 4 |   | 5 |   |
    +---+---+---+---+---+
  4 | 6 | # | 7 |   |   |
    +---+---+---+---+---+
  5 | 8 |   |   |   |   |
    +---+---+---+---+---+
  6 |   | # | # |   | # |
    +---+---+---+---+---+
                         
We can use the following words
* AFT * ALE * EEL * HEEL * HIKE * HOSES * KEEL * KNOT * LASER * LEE * LINE * SAILS * SHEET * STEER * TIE

Models


MiniZinc: crossword.mzn
Choco: CrossWord.java
JaCoP: CrossWord.java
Gecode/R: (Not implemented)
Comet: crossword.co
Gecode: crossword.cpp

Note: I have seen more general models for solving crossword problems in Choco, JaCoP, Gecode/R, and Gecode with constructions other that the simple Elements used here. Since I wanted to compare the same way of solving the problem using the same Element-construct this may be an unfair comparison between the systems. Well, this is at least a finding how to implement this problem by Element...

Explanation of variables
The matrix A is the individual chars of the words (Comet variant):
int A[1..num_words,1..word_len] = 
  [
   [h, o, s, e, s], //  HOSES
   [l, a, s, e, r], //  LASER
   [s, a, i, l, s], //  SAILS
   [s, h, e, e, t], //  SHEET
   [s, t, e, e, r], //  STEER
   [h, e, e, l, 0], //  HEEL
   [h, i, k, e, 0], //  HIKE
   [k, e, e, l, 0], //  KEEL
   [k, n, o, t, 0], //  KNOT
   [l, i, n, e, 0], //  LINE
   [a, f, t, 0, 0], //  AFT
   [a, l, e, 0, 0], //  ALE
   [e, e, l, 0, 0], //  EEL
   [l, e, e, 0, 0], //  LEE
   [t, i, e, 0, 0]  //  TIE
];
overlapping is the matrix of the overlapping cells)
This is the Comet version:
[
 [1, 3, 2, 1],   //  s
 [1, 5, 3, 1],   //  s 
  
 [4, 2, 2, 3],   //  i
 [4, 3, 5, 1],   //  k
 [4, 4, 3, 3],   //  e
  
 [7, 1, 2, 4],   //  l
 [7, 2, 5, 2],   //  e
 [7, 3, 3, 4],   //  e
  
 [8, 1, 6, 2],   //  l
 [8, 3, 2, 5],   //  s
 [8, 4, 5, 3],   //  e
 [8, 5, 3, 5]    //  r
 ];
E is the variable array of which the word to use for the different overlappings. This is in fact the only variable (array) that is needed in the problem, apart from the utility/convenience variables.

The main constraint for the crossword example in each system is stated thus:

MiniZinc:
forall(i in 1..num_overlapping) (
   A[E[overlapping[i,1]], overlapping[i,2]] =  A[E[overlapping[i,3]], overlapping[i,4]]
)
Comet:
  forall(i in 1..num_overlapping) {
    cp.post(A[E[overlapping[i,1]], overlapping[i,2]] ==  A[E[overlapping[i,3]], overlapping[i,4]], onDomains);
  }
Choco:
Note that Choco has a special Element which support two dimensional arrays (matrix), which we use.
for(int I = 0; I < num_overlapping; I++) {
  IntegerVariable tmp = makeIntVar("tmp" + I, 1, 26);
  M.addConstraint(nth(E[overlapping[I][0]], W[overlapping[I][1]], AX, tmp));
  M.addConstraint(nth(E[overlapping[I][2]], W[overlapping[I][3]], AX, tmp));
}
JaCoP:
Here we had used some trickery by using a transposed version of the word matrix since JaCoP has no special Element constraint for two dimensional arrays.
for (int I = 0; I < num_overlapping; I++) {
   FDV tmp = new FDV(store, "TMP" + I, 0, num_words*word_len);
   store.impose(new Element(E[overlapping[I][0]], words_t[overlapping[I][1]], tmp));
   store.impose(new Element(E[overlapping[I][2]], words_t[overlapping[I][3]], tmp));
}
Gecode:
This is more complicated compared to the two Java systems since in Gecode we use an array (of length rows*cols) to simulate the matrix. (There is a Matrix "view" in Gecode but the indices must be of type integer, not IntVar, so it can not be used.) Also, the constraints plus and mult takes IntVar as argument.

The first overlapped crossing is "expanded" like this (Gecode is 0-based):
   A[E[overlapping[i,0]], overlapping[i,1]] // MiniZinc/Comet
   =>
   a1 = A[ E[I*4+0] * word_len + overlapping[I*4+1]] // Gecode
Here is the complete code. The comments hopefully explains what is going on.

First we define an utility function for accessing the element according to the above principle.
/**
 *
 * Special version of element for an array version of a "matrix" words,
 * E is an integer variable array, C is an array of IntVars for
 * the offset j in the words "matrix".
 *
 * The call 
 *    element_offset(*this, words, E[i], word_len_v, C[j], res, opt.icl());
 *
 * corresponds to:
 *    res = words[E[i], j] --> words[E[i]*word_len+J]
 *
 */
void element_offset(Space& space,
                   IntArgs words,
                   IntVar e,
                   IntVar word_len,
                   IntVar c,
                   IntVar res,
                   IntConLevel icl = ICL_DOM) {

      element(space, words, 
              plus(space, 
                   mult(space, 
                        e, 
                        word_len, icl), 
                   c, icl), 
              res, icl);
}

The function is then used as follows:
for(int I = 0; I < num_overlapping; I++) {
   IntVar e1(*this, 0, num_overlapping*4);
   IntVar e2(*this, 0, num_overlapping*4);

   IntVarArray o(*this, 4, 0, num_overlapping*4);
   for(int J = 0; J < 4; J++) {
     post(*this, o[J] == overlapping[I*4+J], opt.icl());
   }

   element(*this, E, o[0], e1, opt.icl());      // e1 = E[I*4+0]
   element(*this, E, o[2], e2, opt.icl());      // e2 = E[I*4+2]

   IntVar a1(*this, 0, num_words*word_len);
      
   element_offset(*this, A, e1, word_len_v, o[1], a1, opt.icl());
   element_offset(*this, A, e2, word_len_v, o[3], a1, opt.icl());
}
(The same element_offset function is also used in the word_square problem below.) It took quite a time to get the function and temporary variables (and their domains) right. With training (and the element_offset as a skeleton) similiar problems should be easier to implement.

Note: this is not a bashing of Gecode. Gecode is a great system and it happens that for this specific problem, Gecode does not has the appropriate support. I should also mention that it was a long time since I programmed in C++ and am little rusty.

As mentioned earlier, I have been very spoiled by the MiniZinc (and Comet) constructs. Also: I'm a very 'lazy' person (in the Perl sense of the word lazy) and likes the agile programming languages - Perl, Ruby, Python etc - much for their high level constructs.

Word square problem

The word problem is a cousin to the crossword problem, and is described in Wikipedia's Word_square:
A word square is a special case of acrostic. It consists of a set of words, all having the same number of letters as the total number of words (the "order" of the square); when the words are written out in a square grid horizontally, the same set of words can be read vertically.
An example of order 7 found by the Comet model where we see that the first row word (aalborg) is also the first column word.

aalborg
abaiser
lantaca
bittman
osamine
recanes
granese

Models

Here are the models for solving the Word square problem:
MiniZinc: word_square.mzn
Choco: WordSquare.java
JaCoP: WordSquare.java
Gecode/R: (Not implemented it Gecode/R)
Comet: word_square.co
Gecode: word_square.cpp

It is somewhat easier than the crossword problem. As before, E is the array of the index of the words to use, and words is an matrix of the words. Also, these models is an experiment of how to read a file, the word list /usr/dict/words (standard on Unix/Linux systems).

MiniZinc
forall(I, J in 1..word_len) (
  words[E[I], J] = words[E[J],I]
)
Comet
  forall(i in 1..word_len) {
    forall(j in 1..word_len) {    
      cp.post(words[E[i], j] == words[E[j],i], onDomains);
    }
  }
JaCoP
// 
// The overlappings (crossings).
// Note that we use a transposed word matrix for the Element.
//
for(int i = 0; i < word_length ; i++) {
    for(int j = 0; j < word_length ; j++) {
        // Comet: words[E[i], j] ==  words[E[j],i]
        FDV tmp = new FDV(store, "tmp" + i + " " + j, 0, dict_size);
        store.impose(new Element(E[i], words[j], tmp));
        store.impose(new Element(E[j], words[i], tmp));
    }
}
Choco
// Constants for the nth constraint below
IntegerVariable [] C = new IntegerVariable[dict_size];
for (int I = 0; I < word_length; I++) {
    C[I] = makeIntVar("C"+I, I,I);
}

// The overlappings (crossings)
for(int I = 0; I < word_length ; I++) {
    for(int J = 0; J < word_length ; J++) {
        // Comet: words[E[i], j] ==  words[E[j],i]
        IntegerVariable tmp = makeIntVar("tmp" + I + " " + J, 0, dict_size);
        M.addConstraint(nth(E[I], C[J], words, tmp));
        M.addConstraint(nth(E[J], C[I], words, tmp));
    }
}
Gecode
Note that this model use the same function element_offset that was used in the Crossword problem. It took some time to realize that it also could be used here.
// convenience variables for the element constraints below
// since element, plus, and mult wants IntVars.
IntVar word_len_v(*this, word_len, word_len);
IntVarArray C(*this, word_len, 0, word_len-1);
for(int i = 0; i < word_len; i++) {
  rel(*this, C[i], IRT_EQ, i, opt.icl());
}

for(int i = 0; i < word_len; i++) {
  for(int j = 0; j < word_len; j++) {
    // words[E[i], j] ==  words[E[j],i]

    IntVar tmp(*this, 0, num_words);

    // tmp == words[E[i], j] --> words[E[i]*word_len+j]
    element_offset(*this, words, E[i], word_len_v, C[j], tmp, opt.icl());

    // tmp == words[E[j], i]  --> words[E[j]*word_len+i]
    element_offset(*this, words, E[j], word_len_v, C[i], tmp, opt.icl());
  }
}

Who killed Agatha?

This is a standard benchmark for theorem proving, also known as The Dreadsbury Mansion Murder Mystery.

Problem formulation from The h1 Tool Suite
Someone in Dreadsbury Mansion killed Aunt Agatha. Agatha, the butler, and Charles live in Dreadsbury Mansion, and are the only ones to live there. A killer always hates, and is no richer than his victim. Charles hates noone that Agatha hates. Agatha hates everybody except the butler. The butler hates everyone not richer than Aunt Agatha. The butler hates everyone whom Agatha hates. Noone hates everyone. Who killed Agatha?
Originally from F. J. Pelletier: Seventy-five problems for testing automatic theorem provers., Journal of Automated Reasoning, 2: 191-216, 1986.

Models


MiniZinc: who_killed_agatha.mzn
JaCoP : WhoKilledAgatha.java
JaCoP : WhoKilledAgatha_element.java
Choco: WhoKilledAgatha.java
Choco: WhoKilledAgatha_element.java
Comet: who_killed_agatha.co
Gecode: who_killed_agatha.cpp

In Some new Gecode models I wrote about the findings in implemented this problem in Gecode and compared to Comet/MiniZinc.

The models use two 3x3 matrices for representing the two relations hates and richer with 0..1 as domain (i.e. boolean). The Element constraints is used for implementing the condition A killer always hates, and is no richer than his victim. where the_killer is an integer variable; the_victim is in some models replaced by agatha (the integer 0). The interesting thing here is that at least one of the indices are integer variables, which caused the difficulties in the two problems above.

These models also use a lot of boolean constructs. A comparison of how these are implemented in the different CP systems may be described in a future blog post.

MiniZinc
hates[the_killer, the_victim] = 1 /\
richer[the_killer, the_victim] = 0
Comet:
m.post(hates[the_killer, the_victim] == 1);
m.post(richer[the_killer, the_victim] == 0);
Note: In the models below I have simplified the problem by using agatha (defined as the integer 0) instead of the integer variable the_victim. This is not a problem since we know that Agatha is the victim, and is the reason why the Element is easier to use than for Crossword and Word square.

JaCoP variant 1 (no Element):
JaCoP don't have a direct support for the case when the index i (in matrix[i][j]) is an integer variable so the first variant of modeling the condition A killer always hates, and is no richer than his victim. does not use Element at all. In WhoKilledAgatha.java we simply loop over all integers (0..2), check if "this" i equals the_killer and then we can use two integers for accessing the matrices. Also, note the IfThen construct.
for(int i = 0; i < n; i++) {
    store.impose(
                 new IfThen(
                            new XeqC(the_killer, i),
                            new XeqC(hates[i][agatha], 1)
                            )
                 );
    store.impose(
                 new IfThen(
                            new XeqC(the_killer, i),
                            new XeqC(richer[i][agatha], 0)
                            )
                 );
}
This was the first variant I implemented, but then I recall the "trickery" used in Crossword and Word square where the matrices where transposed and Element could be used. The problem with this approach is that all constraints must be rewritten in a way that may be confusing. Come to think of it, maybe the names of the matrices should have been changed to is_hated_by and poorer.

JaCoP variant 2 (transposed matrices, Element)
This method of transposition and using Element is implemented in WhoKilledAgatha_element.java. The constraint is now much simpler:
int shift = -1;
for(int i = 0; i < n; i++) {
    store.impose(new Element(the_killer, hates[agatha], one, shift));
    store.impose(new Element(the_killer, richer[agatha], zero, shift));
}
Note: the Element in JaCoP defaults to start index as 1, but has support for shifting it to 0, by using -1 as the shift parameter. Choco variant 1 (no Element)
I implemented exact the same principle that was used in the two JaCoP model in the two Choco models. The first - no Element - is WhoKilledAgatha.java:

for(int i = 0; i < n; i++) {
    m.addConstraint(implies(
                            eq(the_killer,i),
                            eq(hates[i][agatha], 1)
                            )
                    );
    m.addConstraint(implies(
                            eq(the_killer, i),
                            eq(richer[i][agatha], 0)
                            )
                    );
}
Choco variant 2 (transposed matrices, nth)
Note: one and zero are integer variables since nth cannot handle plain integers as the last argument.
for(int i = 0; i < n; i++) {
   m.addConstraint(nth(the_killer, hates[agatha], one));
   m.addConstraint(nth(the_killer, richer[agatha], zero)


Comments

Here we have seen - not surprisingly - that using the Element constraint is quite different depending of which CP system we use and it can be easy or not so easy. It was my explicit intension to see how to solve the same problem as similiar as possible. We should also note that most (if not all) problems can be modeled in many ways, some not using Element at all.

One last comment: The two Java models of Who killed Agatha? took quite a long time to implement. The main reason for that was not the the handling of Element but was due to a bug of confusing the two matrices in one of the conditions. Sigh.

April 04, 2009

MiniZinc/FlatZinc support in SICStus Prolog version 4.0.5

The other day I downloaded a demo of SICStus Prolog version 4.0.5, with the sole intention of testing the new support for MiniZinc/FlatZinc. The version I downloaded was for Linux 32 bit glibc 2.7 (my machine use glibc 2.6 but it seems to work well anyway).

The support for MiniZinc/FlatZinc in SICStus Prolog is described more in 10.35 Zinc interface - library(zinc). Some restrictions are described in Notes:
  • Domain variables
    Only variables with finite integer domains are supported. This includes boolean variables which are considered finite integer domain variables with the domain 0..1.
  • Optimization problems
    Only variables with finite integer domains can be optimized in minimize and maximize solve items. The int_float_lin/4 expression as described in the FlatZinc specification is thus not supported.
  • Solve annotations
    • The solve annotations currently supported are int_search/4, bool_search/4, and labelling_ff/0.
    • The FlatZinc specification describes several exploration strategies. Currently, the only supported exploration strategy is "complete".
    • When no solve annotation is given, a most constrained heuristic is used on all problem variables (excluding those that have a var_is_introduced annotation, see below). This corresponds to labeling/2 of library(clpfd) with the option ffc.
    • The choice method "indomain_random" as described in the FlatZinc specification uses random_member/2 of library(random). The random generator of SICStus is initialized using the same seed on each start up, meaning that the same sequence will be tried for "indomain_random" on each start up. This behavior can be changed by setting a different random seed using setrand/1 of library(random).
  • Constraint annotations
    Constraint annotations are currently ignored.
  • Variable annotations
    Variable annotations are currently ignored, except var_is_introduced, which means that the corresponding variable is not considered in any default labeling (such as when no search annotation is given or when the labelling_ff search annotation is given).

Testing

For testing the MiniZinc solver I used exactly the same principle as I do for the ECLiPSe solver, so hooking it up into in my system was very easy. All this is done via a Perl script of my own. It consists of generating a Prolog file, here called prolog_file.pl with the content below . model.mzn is the MiniZinc file to run, and number_of_solutions is the number of solutions to generate (an integer, or all for all solutions).

%
% Generated by mzx.pl
%
:- use_module(library(zinc)).

go :-
   mzn_run_file('model.mzn',
                       [solutions(number_of_solutions)]),
   halt.
And then running the following command (from the Perl program):
sicstus -l prolog_file.pl --goal go.

Findings

Most things works very well and with about the same performance as the ECLiPSe solver. I will investigate some more before (perhaps) buying a private license of SICStus Prolog (or upgrading from my old version 3.9.1 if that is possible). However, I did find some problems.
  • global_cardinality/2
    The support for the builtin global_cardinality/2 is broken. The following error occurs:
    ! Existence error
    ! `global_cardinality/2' is not defined
    
    Example: sudoku_gcc.mzn.

    There is an easy fix which works but is slower that using an builtin support: in the file globals.mzn (in the SICStus Prolog distribution), just use the decomposition variant (the commented) instead.
  • cumulative
    For the model furniture_moving.mzn the following error occurs:
    ! Instantiation error in argument 2 of user:cumulative/2
    ! goal:  cumulative([task(_4769,30,_4771,3,1),task(_4777,10,_4779,1,2),task(_4785,15,_4787,3,3),task(_4793,15,_4795,2,4)],[limit(_4801)])
    
    Note: the model cumulative_test_mats_carlsson.mzn works without problem (this is a simple example from one of Mats Carlsson's lecures).
  • integer overflow
    For the Grocery example (grocery.mzn) an integer overflow error is thrown:
    ! Item ending on line 3:
    ! Representation error in argument 1 of user:range_to_fdset/2
    ! CLPFD integer overflow
    ! goal:  range_to_fdset(1..359425431,_4771)
    
    Notes: MiniZinc's solver also overflows, so this is probably not a big thing. The solvers for Gecode/FlatZinc and ECLiPSe ic handles this problems correctly, though.
  • Statistics
    I have not seen any relevant statistics (e.g. number of failures, nodes, propagations etc) for the SICStus MiniZinc solver. The standard SISCtus Prolog predicate statistics/0 is somewhat useful, but is not what is needed when writing MiniZinc models and comparing with other versions and/or solvers.

    What I have in mind is something like the statistics from Gecode (and Gecode/FlatZinc):
    runtime:       50
    solutions:     1
    propagators:   8
    propagations:  3995
    nodes:         249
    failures:      124
    peak depth:    12
    peak memory:   27 KB
    

Final note

I have contacted the developers of SICStus Prolog about these things. They responsed that the bugs are now fixed and will be included in the next version (4.0.6). They also indicated that more detailed statistics may make it in a later version. That is great!

I have now also added the SICStus Prolog solver on my MiniZinc page.

March 24, 2009

Gecode version 3.0.1 and Gecode/FlatZinc 1.5 released

Version 3.0.1 of Gecode is relased. It contains mostly bug fixes:

This is a bug fix release fixing an embarassing bug in reified Boolean linear constraints.
  • Finite domain integers

    • Bug fixes
      • IntSetArgs no longer inherit from PrimArgArray, which was wrong as IntSet is no primitive type and hence does not support vararg initializers. (minor)
      • Fixed bug in reified Boolean linear constraints (an optimization is currently disabled and will be active in the next release: the optimization was incorrect and was never tested). (major, thanks to Alberto Delgado)


  • Example scripts
    • Additions
    • Bug fixes
      • The examples now pass the c-d and a-d command line options correctly to Gist. (minor)
      • The Steel Mill Slab Design example had two bugs in the search heuristic and a missing redundant constraint. (minor, bugzilla entry, thanks to Chris Mears)




Gecode/FlatZinc has been updated to version 1.5. The bug fix here is very interesting (and exiting) since Gist now also works with Gecode/FlatZinc.

Gist in Gecode and Gecode/FlatZinc
Gist for Gecode has been around for time, and was officially in the distribution in version 3.0.0. In Modeling with Gecode (PDF) there is a section how to program and use Gist.

Stable Marriage Problem.
As a first example of Gist with Gecode/FlatZinc (i.e. MiniZinc models): stable_marriage.pdf is a PDF file which shows the full search tree for the MiniZinc model stable_marriage.mzn. The green nodes are solutions, the blue nodes are choice points, the red squares are failures, and the big red triangles are a sub tree of failures.

Running the model with fz -mode stats stable_marriage.fzn gives the following result.

wife : [7, 5, 9, 8, 3, 6, 1, 4, 2]
husband: [7, 9, 5, 8, 2, 6, 1, 4, 3]
----------
wife : [6, 5, 9, 8, 3, 7, 1, 4, 2]
husband: [7, 9, 5, 8, 2, 1, 6, 4, 3]
----------
wife : [6, 4, 9, 8, 3, 7, 1, 5, 2]
husband: [7, 9, 5, 2, 8, 1, 6, 4, 3]
----------
wife : [6, 1, 4, 8, 5, 9, 3, 2, 7]
husband: [2, 8, 7, 3, 5, 1, 9, 4, 6]
----------
wife : [6, 4, 1, 8, 5, 7, 3, 2, 9]
husband: [3, 8, 7, 2, 5, 1, 6, 4, 9]
----------
wife : [6, 1, 4, 8, 5, 7, 3, 2, 9]
husband: [2, 8, 7, 3, 5, 1, 6, 4, 9]
----------
==========

runtime: 10
solutions: 6
propagators: 346
propagations: 12426
nodes: 67
failures: 28
peak depth: 8
peak memory: 132 KB


Also see My MiniZinc Page for MiniZinc models which may be used with Gecode/FlatZinc.

March 12, 2009

MiniZinc Challenge 2008 Results

The MiniZinc Challenge 2008 was held in the summer 2008. I don't know exactly when the result was published (probably in the last week), but here it is; MiniZinc Challenge 2008 Results.

The summary states (higher result better):


The results of the challenge are available here; congratulations go to GeCode on a very convincing win! The summary of results is given in this table:

Contestant Total Score
eclipse_fd 787.3
eclipse_ic 938.8
g12_fd 1655.1
gecode 3418.8

Congratulations to the Gecode team!

This result confirms my impression about the MiniZinc solvers: Gecode/flatzinc is my favorite solver, since it is often very fast, and also that it shows important statistics. See Pi Day Sudoku 2009 for some examples of the latter.

As of writing (2009-03-12) the links to the models and data files in the result page don't work but all the files is included in the last ROTD (release of the day).

By the way, two of the models are from my MiniZinc collection:
* debruijn_binary.mzn: de Bruijn sequences
* quasiGroup7.mzn: quasigroup problem 7

March 11, 2009

Pi Day Sudoku 2009

Brainfreeze Puzzles has a Pi day Sudoku competition:


Pi Day is a celebration of the number π that occurs every March 14 (3.14...). Math geeks all over the world celebrate their math-geekiness on this day with pie-eating contests, recitations of the digits of pi, and occasionally fundraisers where math faculty get hit in the face with pies. At Brainfreeze Puzzles we celebrate Pi Day - how else? - with a Sudoku puzzle.

...

Rules: Fill in the grid so that each row, column, and jigsaw region contains 1-9 exactly once and π [Pi] three times.

The puzzle is also here (as a PDF file)

Models

I programmed this puzzle in both Comet (using the constraint programming module) and MiniZinc using different solvers. The models use the same principle of solving the problem: a (homebrewn) all different except 0 constraint and counting exactly 3 occurrences of Pi for each row/column and region.

All the solvers give the same unique answer, which of course may mean that there is some systematic error on my part.


Comet


Comet solves the problem in about 9 seconds with the following statistics, using exploreall to make sure that it is a unique solution.

time: 8821
#choices = 9999
#fail = 16231
#propag = 23168886

Somewhat surprisingly selecting j before i in the forall loop, gives an improvment from 15 seconds to 9 seconds. This was the same labeling strategy that was used for improving Nonogram (see Comet: Nonogram improved: solving problem P200 from 1:30 minutes to about 1 second). Here is the labeling:

// reversing i and j gives faster solution
forall(j in 1..n, i in 1..n: !x[i,j].bound()) {
tryall(v in V : x[i,j].memberOf(v)) by(v)
m.label(x[i,j], v);
onFailure
m.diff(x[i,j], v);
}

MiniZinc and Gecode/flatzinc

Interestingly, the same improvement was observed with MiniZinc and the Gecode/flatzinc solver (the interface to Gecode). Here is the labeling function, with the same swapping of i and j (I didn't know that this mattered in MiniZinc).


solve :: int_search([x[i,j] | j, i in 1..n], "input_order", "indomain", "complete") satisfy;

With the swap, the time went from 16 seconds to 10 seconds with the following statistics:

propagators: 9460
propagations: 13461717
failures: 11221
clones: 11221
commits: 31478
peak memory: 5638 KB
mzx.pl sudoku_pi.mzn --num 2 9,41s user 0,72s system 99% cpu 10,139 total

Fzntini, FlatZinc, ECLiPSe and searching for the first solution

MiniZinc's builtin solver flatzinc and fzntini (a SAT solver) only shows one solution and no statistics.

* fzntini solve problem in about 13 seconds
* flatzinc takes 3 seconds.

The ECLiPSe's MiniZinc solver ic takes 6.5 seconds using the same labeling as Gecode/flatzinc. No statistics is shown for this solver either.

For comparison both Comet and Gecode/flatzinc shows the first solutions in 2.5 seconds.

Solution

Since it is a competition I won't show my solution or the models on/via the blog. Sorry about that.

Reference

I learned about Pi Days Sudoku from the 360 blog: Pi Day Sudoku 2009.

February 28, 2009

Comet: regular constraint, a much faster Nonogram with the regular constraint, some OPL models, and more

Since the last time some more Comet model has been written.

Much faster Nonogram model using the regular constraint

In More Comet models, e.g. Nonogram, Steiner triplets, and different set covering problems I presented nonogram.co, a solver for Nonogram puzzles, and also noted that is was quite slow: Well, it is nice to have some more optimization to do (or more probably, a complete remodelling)....

Inspired by the announcement this week of the ECLiPSe example of Nonogram solver using the regular constraint (see nono_regular.ecl.txt and regular.ecl.txt) - and also my earlier Gecode/R model nonogram.rb which used "regular expressions" - I created a Nonogram model in Comet using the regular constraint: nonogram_regular.co.

Let us first look at the regular constraint.


Regular constraint
The regular constraint (see the Comet model regular.co for my implementation) is a global constraint using a DFA (deterministic finite automata) which accepts (or not accepts) an input sequence given a "transition matrix". The constraint was presented by Gilles Pesant in "A Regular Language Membership Constraint for Finite Sequences of Variables" (2004).

My implementation of the regular constraint is heavily borrowed from MiniZinc's builtin regular predicate (from lib/zinc/globals.mzn); the translation to Comet was very straightforward (omitting just some asserts).

An exemple of the usage of the constraint: We want to match the regular expression "123*21" (i.e. first "1", then "2", then zero or more "3", then "2", and last a "1"). Note: The length of the sequence is 10, so there must be 6 "3":s since "123" and "21" are "anchored" at the beginning and the end.

int len = 10;
int n_states = 5; // number of states
int input_max = 3; // the states are 1,2, and 3
int initial_state = 1; // we start with the 1 state
set{int} accepting_states = {n_states}; // This is the last state
// The transition matrix
int transition_fn[1..n_states, 1..input_max] =
[[2, 0, 0], // transitions from state 1: "1" -> state 2
[0, 3, 0], // transitions from state 2: "2" -> state 3
[0, 4, 3], // transitions from state 3: "2" -> state 4 | "3" -> state 3 (loop)
[5, 0, 0], // transitions from state 4: "1" -> state 5
[0, 0, 0]]; // transitions from state 5: END state

...

exploreall {
regular(reg_input, n_states, input_max, transition_fn, initial_state, accepting_states);
} using {
// ....
}

The unique sequence resulting from this automata is thus 1 2 3 3 3 3 3 3 2 1.

For using regular in the Nonogram problem, the automata for each row/column clue, must be built, preferably by a function. In regular.co this is done with the function make_transition_matrix, which by far was the hardest part of this problem (and surely could be written in a smoother way).

For the Nonogram clue [3,2,1] - which represents the regular expression "0*1110*110*10*" - the following automata (transition matrix) is generated:
1 2
0 3
0 4
5 0
5 6
0 7
8 0
8 9
9 0


Note that the regular function uses 0 (zero) as the failing state, so the states must start with 1. This is taken care in the model as the resulting value is subtracted by 1.

As usual in my models, the regular constraint is just a "convenience function", i.e. not using special written propagation methods etc.

The regular constraint has - of course - more applications than for Nonogram solving. I plan to look more into this.


The Nonogram solver: nonogram_regular.co
After the regular constraint and the automata generator was written, it was quite easy to change the old Nonogram solver to use these new tools. The result is in nonogram_regular.co. I was quite curious how fast this version should be compared to the former slow model. In short: it is much faster.

As comparison the Lambda instace took about 12 seconds with the old model; with the new model it takes 0.5 seconds, which is a nice improvement. I never finished a run for Nonunique problem with the older model; the new model takes 0.5 seconds (including the startup of the Comet program). Etc.

The P200 problem nonogram_p200.co now takes 2:05 minutes, which can be compared with 57 seconds for the Gecode/R model. Thus, the Comet model is still slow compared to Gecode version and the ECLiPSe version, which solves the P200 problem in just a second or two. Maybe some creative labeling or a "proper" regular constraint can fasten things up...

Update some hour later: One way to gain about 30 seconds to 1:30 minutes on the P200 problem was to explicitly state the consistency level to onDomain, e.g. cp.post(a[m] == q0, onDomains), and to use another labeling strategy:
forall(i in 1..rows, j in 1..cols : !board[i,j].bound()) {
// label(board[i,j]); // the former strategy
tryall(v in 1..2) m.post(board[i,j] == v, onDomains);
}


Some more Nonogram problems have been coded:

An aside about regular expressions
I have been very interested in regular expressions (especially the more expressive Perl type) for a long time. In 1997 I wrote a simple Perl program MakeRegex which returns a regular expression given a list of words. It was later Appletized in MakeRegexApplet. Now there are better programs/modules for this.

OPL Models

One of my projects is to translate the OPL models from Pascal Van Hentenryck "The OPL Optimization Programming Language" into Comet. Even if OPL and Comet are different languages the reading of the book has been very awarding. Thanks again, Pascal!

Some OPL models already have been published, but now I've been more systematic and started from the beginning. More models to come.

Finding: arrays in a tuple
In Comet: New models, e.g. Einstein puzzle, KenKen, Kakuro, Killer Sudoku, Stigler's Diet problem, translations of OPL models I wrote

One big difference here: In Comet it is not allowed to have an array in a tuple, so the use data must be a separate matrix.

I've found a way of using arrays in a tuple, by first initialize the values in a matrix and then - by using the all function - "slice" the values into the tuple array. This has been done in the models fixed_charge2.co and production3.co.

Example from fixed_charge2.co:

tuple productType {
int profit;
set{Machines} machines;
int[] use;
}
int use[Products, Resources] = [[3,4],[2,3], [6,4]];
productType Product[Products] =
[productType(6, {shirtM}, all(i in Resources) use[shirt,i]),
productType(4, {shortM}, all(i in Resources) use[shorts,i]),
productType(7, {pantM}, all(i in Resources) use[pants,i])];


Combining different models

The alphametic problem SEND + MOST = MONEY has the additional requirement to maximize the value of MONEY. The older model send_most_money.co does that and nothing more.

One natural extension of the problem is the following:
* first find the maximum value of MONEY
* then find all the solutions with this value.

The model send_most_money2.co has two functions:
* smm which is just the constraint for the alphametic problem
* send_most_money which has a parameter money.

send_most_money is first called with 0, indicating that is should maximize the value of MONEY, and then returns that value. The next call to send_most_money is with the calculated MONEY value, which indicates that all solutions should be generated.

The answer is

check all solutions for MONEY = 10876
x[9,7,8,2,1,0,4,6] MONEY: 10876
x[9,7,8,4,1,0,2,6] MONEY: 10876


Jim Orlin's Logic Puzzle

In Colored letters, labeled dice: a logic puzzle Jim Orlin stated a logic puzzle:
My daughter Jenn bough a puzzle book, and showed me a cute puzzle. There are 13 words as follows: BUOY, CAVE, CELT, FLUB, FORK, HEMP, JUDY, JUNK, LIMN, QUIP, SWAG, VISA, WISH.

There are 24 different letters that appear in the 13 words. The question is: can one assign the 24 letters to 4 different cubes so that the four letters of each word appears on different cubes. (There is one letter from each word on each cube.) It might be fun for you to try it. I’ll give a small hint at the end of this post. The puzzle was created by Humphrey Dudley.

...

If anyone wants to write an Excel spreadsheet and solve it via integer programming, please let me know. I’d be happy to post the Excel spreadsheet if you send it to me, or link to it if you post it and send me the URL.

This was a fun puzzle, so I modeled the problem in labeled_dice.co, and mailed Jim the link.

Some days later he wrote Update on Logic Puzzle where the contributed models were presented. There where another Comet model WordPuzzleSchaus.co by Pierre Schaus, one of Comet's developers. Pierres model use an different and more elegant approach than mine.

Jim also linked to some other logical puzzles from the great collection brownbuffalo.sourceforge.net (which have solutions in ECLiPSe/Prolog). One of these puzzles was Building Blocks, in the same vein as his labeled dice puzzle. Hence I had to make a Comet model of this problem: building_blocks.co.

He concludes the post with the following:

Incidentally, the ease for solving it using Constraint Programming makes me think that Constraint Programming should be considered a fundamental tool in the OR toolkit.


Other Comet models this week

Here are the other Comet models created/published this week. Some of them where very easy to do since they are translations of my MiniZinc models.

January 18, 2009

Some other Gecode/R models, mostly recreational mathematics

Here are some new Gecode/R models. They are mostly in recreation mathematics since I still learning the limits of Gecode/R and want to use simple problems.


SEND + MOST = MONEY


This is a simple alphametic puzzle: SEND + MOST = MONEY and the model is based on the example in the Gecode/R distribution (send_most_money.rb). The difference with its cousin SEND + MORE = MONEY is that there are many solutions, and we want to maximize the value of MONEY (10876).

The model is send_most_money2.rb.

Now, there are actually 4 different solutions for MONEY = 10876 which we want to find.


s e n d m o s t y
3 7 8 2 1 0 9 4 6
5 7 8 2 1 0 9 4 6
3 7 8 4 1 0 9 2 6
5 7 8 4 1 0 9 2 6

In order to do that, the model first solves the maximation problem and assigns the value (10876) to max_value; the second steps finds all 4 solutions. These two steps use the same code with the difference that the second step activates the following constraint

if (max_money > 0) then
money.must == max_money
end

Nothing fancy, but it quite easy to to it in Ruby and Geocde/R.


Steiner triplets


This is another standard problem in the constraint programming (and recreational mathematics) community: Find a set of triplet of numbers from 1 to n such that any two (different) triplets have at most one element in common.

For more about this, see Mathworld: SteinerTripleSystem, and Wikipedia Steiner_system.

The Gecode/R model is steiner.rb.

The problem can be simply stated as:

nb = n * (n-1) / 6 # size of the array
sets_is_an set_var_array(nb, [], 1..n)
sets.must.at_most_share_one_element(:size => 3)
branch_on sets, :variable => :smallest_unknown, :value => :min

This, however, is very slow (and I didn't care to wait that long for a solution). I tried some other branch strategies but found none that made some real difference.

When the following constraint was added, it really fasten things up. This in effect the same as the constraint sets.must.at_most_share_one_element(:size => 3) above:

n.times{|i|
sets[i].size.must == 3
(i+1..n-1).each{|j|
(sets[i].intersection(sets[j])).size.must <= 1
}
}

The first 10 solutions for n = 7 took about 1 seconds. The first solution is:

Solution #1
{1, 2, 3}
{1, 4, 5}
{1, 6, 7}
{2, 4, 6}
{2, 5, 7}
{3, 4, 7}
{3, 5, 6}
memory: 25688
propagations: 302
failures: 0
clones: 2
commits: 16

For n = 9, 10 solutions took slightly longer, 1.3 seconds.


Devil's Words


Gecode/R model: devils_word.rb.

Devil's word is a "coincidence game" where a the ASCII code version of a name, often a famous person's, is calculated to sum up to 666 and some points is made about that fact (which of course is nonsense).

There are 189 different ways my own name HakanKjellerstrand (where the umlaut "å" in my first name is replaced with "a") can be "devilized" to 666. With output it took about 2 seconds to generate the solutions, without output it took 0.5 seconds.

The first solution is:

+72 +97 +107 +97 +110 -75 +106 +101 +108 -108 -101 +114 +115 +116 +114 -97 -110 -100


Also, see
* MiniZinc model: devils_word.mzn
* my CGI program: Devil's words.

And see Skeptic's law of truly large numbers (coincidence) for more about coincidences. The CGI program mentioned above was presented in the swedish blog post Statistisk data snooping - att leta efter sammanträffanden ("Statistical data snooping - to look for coincidences") which contains some more references to these kind of coincidences.


Pandigital Numbers, "any" base


Pandigital numbers are a recreational mathemathics construction. From MathWorld Pandigital Number

A number is said to be pandigital if it contains each of the digits
from 0 to 9 (and whose leading digit must be nonzero). However,
"zeroless" pandigital quantities contain the digits 1 through 9.
Sometimes exclusivity is also required so that each digit is
restricted to appear exactly once.

The Gecode/R model pandigital_numbers.rb extends this to handle "all" bases. Or rather bases from 2 to 10, since larger bases cannot be handled by the Gecode solver.

For base 10 using the digits 1..9 there are 9 solutions:

4 * 1738 = 6952 (base 10)
4 * 1963 = 7852 (base 10)
18 * 297 = 5346 (base 10)
12 * 483 = 5796 (base 10)
28 * 157 = 4396 (base 10)
27 * 198 = 5346 (base 10)
39 * 186 = 7254 (base 10)
42 * 138 = 5796 (base 10)
48 * 159 = 7632 (base 10)

For base 10, using digits 0..9, there are 22 solutions.

Here is the number of solutions for base from 5 to 10 and start either 0 or 1 (there are no solutions for base 2..4):
* base 2, start 0: 0
* base 2, start 1: 0
* base 3, start 0: 0
* base 3, start 1: 0
* base 4, start 0: 0
* base 4, start 1: 0
* base 5, start 0: 0
* base 5, start 1: 1
* base 6, start 0: 0
* base 6, start 1: 1
* base 7, start 0: 2
* base 7, start 1: 2
* base 8, start 0: 4
* base 8, start 1: 4
* base 9, start 0: 10
* base 9, start 1: 6
* base 10, start 0: 22
* base 10, start 1: 9

See also the MiniZinc model pandigital_numbers.mzn with a slightly different approach using the very handy MiniZinc construct exists instead of looping through the model for different len1 and len2.


SEND + MORE = MONEY (any base)


And talking of "all base" problems, send_more_money_any_base.rb is a model for solving SEND+MORE=MONEY in any base from 10 and onward.


Pattern
The number of solutions from base 10 and onward are the triangle numbers, i.e. 1,3,6,10,15,21,... I.e.

* Base 10: 1 solution
* Base 11: 3 solutions
* Base 12: 6
* Base 13: 10
* Base 14: 15
* Base 15: 21
* etc

For more about triangle numbers, see Wikipedia Triangular numbers.

There seems to be a pattern of the solutions given a base b:

 
           S  E  N  D  M  O  R  Y
           ----------------------
 Base 10:  9  5  6  7  1  0  8  2
 Base 11: 10  7  8  6  1  0  9  2
          10  6  7  8  1  0  9  3
          10  5  6  8  1  0  9  2
 Base 12:
          11  8  9  7  1  0 10  3
          11  8  9  6  1  0 10  2
          11  7  8  9  1  0 10  4
          11  6  7  9  1  0 10  3
          11  6  7  8  1  0 10  2
          11  5  6  9  1  0 10  2
 ...
 Base 23:
          22 19 20 18  1  0 21 14
          22 19 20 17  1  0 21 13
 ...

Some patterns:

S: always base-1 e.g. 9 for base 10
M: always 1 e.g. 1 any base
0: always 0 e.g. 0 any base
R: always base-2 e.g. 8 for base 10
E, N, D: from base-3 down to 5 e.g. {5,6,7,8,9} for base 10
Y: between 2 and ??? e.g. {2,3,4} for base 12


I haven't read any mathematical analysis of these patterns. There is an article in math Horizons April 2006: "SENDing MORE MONEY in any base" by Christopher Kribs-Zaleta with the summary: Dudeney's classic "Send More Money" cryptarithmetic puzzle has a unique solution in base ten. This article generalizes explores solutions in bases other than ten. I haven't read this article, though.


Also, see my MiniZinc model send_more_money_any_base.mzn.

January 15, 2009

Some models in Gecode/R (Ruby interface to Gecode)

Gecode/R is a great Ruby interface to Gecode (implemented in C++). At last I now have done some modeling in Gecode/R and it's nice.

The models and some other information about Gecode/R can be seen at My Gecode/R page.

Since Ruby is a very high level language, the modelling can be done "high levelish", more so than for example with the Java solvers Choco and JaCoP. And that is a feature I really like.

I'm still learning Gecode/R and have surely missed some stuff. There are not many examples in the package (these are also commented at Examples). Some things have been of great help
* The Sitemap
* The Documentation, especially the RDocs
* The test files in the specs directory.


An example: Survo Puzzle
From Wikipedia's Survo Puzzle


In a Survo puzzle the task is to fill an m * n table by integers 1,2,...,m*n so
that each of these numbers appears only once and their row and column sums are
equal to integers given on the bottom and the right side of the table.
Often some of the integers are given readily in the table in order to
guarantee uniqueness of the solution and/or for making the task easier.

E.g. the puzzle 128/2008 is presented with the following clues:


* * * * * * 30
* * 18 * * * 86
* * * * * * 55
22 11 42 32 27 37

where * marks a digit to find. The number to the right is the row sums which the row must sum to, and the last row is the column sums.

The unique solution of the problem is (by running the program ruby survo_puzzle.rb survo_puzzle_128_2008.txt)


Solution #1
4 1 10 5 3 7 = 30
12 8 18 16 15 17 = 86
6 2 14 11 9 13 = 55
= = = = = =
22 11 42 32 27 37

Statistics:
propagations: 3784
failures: 174
clones: 179
memory: 25740
commits: 461
Number of solutions: 1

The relevant constraint programming code is below (slightly edited). I think it's quite nice.


....
def initialize(clues, rowsums, colsums)
r = rowsums.length # number of rows
c = colsums.length # number of columns
x_is_an int_var_matrix(r, c, 1..r*c) # the solution matrix
x.must_be.distinct # all values in x must be distinct
# values in clues with values > 0 are copied straight off to x
r.times{|i|
c.times{|j|
x[i,j].must == clues[i][j] if clues[i][j] > 0
}
}
r.times{|i| x[i,0..c-1].sum.must == rowsums[i] } # check row sums
c.times{|j| x.transpose[j,0..r-1].sum.must == colsums[j] } # check column sums
branch_on x, :variable => :smallest_size, :value => :min
end

The full program is here: survo_puzzle.rb, and three data files:
survo_puzzle_126_2008.txt
survo_puzzle_127_2008.txt
survo_puzzle_128_2008.txt


The Gecode/R models
Below are the models shown at My Gecode/R page. The selection is mostly for comparison with models implemented in other constraint programming languages, see respectively:
My MiniZinc page (which has a lot more models)
My JaCoP page
My Choco page

The models first written (e.g. diet, least_diff) are very "un-Rubyish" since I wanted to model straight after the MiniZinc models. Maybe I Rubyfi these later.

After a while the modeling went quite easy, and both de Bruijn and Minesweeper was done surprisingly fast. I do, however, miss MiniZinc's sum construct, since it which make some things easier (e.g. the neighbour summing in minesweeper.rb).

The execution time of the models os approximately the same time as the corresponding MiniZinc model with the Gecode/flatzinc solver which is normally quite fast. The big exception of these examples is coins_grid which seems to be slow for the constraint programming systems, but fast with linear programming systems, e.g. with the MiniZinc solvers ECLiPSe/ic and MiniZinc/mip.

References to the problems etc are in the header of the model.

January 06, 2009

Map coloring problem: Lichtenstein

In The Chromatic Number of Liechtenstein bit-player asked (2008-10-28)) the following about coloring the map of Lichtenstein:

It seems that Liechtenstein is divided into 11 communes, which emphatically do not satisfy the connectivity requirement of the four color map theorem. Just four of the communes consist of a single connected area (Ruggell, Schellenberg and Mauren in the north, and Triesen in the south). All the rest of the communes have enclaves and/or exclaves.

...

In the map above, each commune is assigned its own color, and so we have an 11-coloring. It’s easy to see we could make do with fewer colors, but how many fewer? I have found a five-clique within the map; that is, there are five communes that all share a segment of border with one another. It follows that a four-coloring is impossible. Is there a five-coloring? What is the chromatic number of Liechtenstein?

I wrote a MiniZinc model for this minimizing problem: lichtenstein_coloring.mzn.

The model has two variable arrays:
* color_communes: the color of the 11 communes
* color: the color of the 29 en-/exclaves

Objective: minimize the number of colors used.

The Gecode/flatzinc solver gives the following solutions in less than 1 second, which states that 5 different colors (n_colors) is sufficient. The model allows for up to 11 different colors, hence the large color numbers.

n_colors: 5
color_communes: [1, 1, 10, 8, 8, 1, 9, 9, 8, 10, 11]
color: [1, 1, 1, 1, 1, 1, 10, 10, 8, 8, 8, 8, 8, 1, 9, 9, 9, 9, 9, 9, 8, 10, 10, 11, 11, 11, 11, 11, 11]

Optimal solution found.

runtime: 290
solutions: 4
propagators: 1235
propagations: 1992711
failures: 1045
clones: 1046
commits: 2757
peak memory: 1414 KB

Times for other MiniZinc solvers:
* Minizinc's flatzinc: 2 seconds,
* Minizinc's fdmip: 2 seconds,
* ECLiPSe's ic: 4 seconds
* tinifz: 5.5 seconds


Also see
The chromatic number of Lichtenstein by Michi (from whom I borrowed the edges).

January 05, 2009

Tom Schrijvers: "Monadic Constraint Programming" (in Haskell)

Tom Schrijvers presents in the blog post Monadic Constraint Programming a draft version of the paper Monadic Constraint Programming written by him, Peter Stuckey, and Philip Wadler:


A constraint programming system combines two essential components: a constraint solver and a search engine. The constraint solver reasons about satisfiability of conjunctions of constraints, and the search engine controls the search for solutions by iteratively exploring a disjunctive search tree defined by the constraint program. In this paper we give a monadic definition of constraint programming where the solver is defined as a monad threaded through the monadic search tree. We are then able to define search and search strategies as first class objects that can themselves be built or extended by composable search transformers. Search transformers give a powerful and unifying approach to viewing search in constraint programming, and the resulting constraint programming system is first class and extremely flexible.


Prototype code in Haskell can be downloaded here.

December 29, 2008

Temporal reasoning model in MiniZinc

temporal_reasoning.mzn is a MiniZinc model of temporal reasoning. The example is from From Krzysztof R. Apt Principle of Constraint Programming, page 23ff:

The meeting ran non-stop the whole day.
Each person stayed at the meeting for a continous period of time.
The meeting began while Mr Jones was present and finished
while Ms White was present.
Ms_White arrived after the meeting has began.
In turn, Director Smith, was also present but he arrived after Jones had left.
Mr Brown talked to Ms White in presence of Smith.
Could possibly Jones and White have talked during this meeting?

The coding was inspired by the ECLiPSe (Prolog) model in Apt's presentation of chapter 2, ch2-sli.pdf.gz (gzipped PDF file), slides 15ff.


Also see My MiniZinc page for other MiniZinc models.

Welcome to my My Constraint Programming Blog

Welcome to my My Constraint Programming Blog!

This is an extension of my "normal" swedish blog hakank.blogg, and will contain news etc about constraint programming and related paradigms. It will also link to my newly written constraint programming models.

As stated (in swedish) in Constraint programming-nyheter samt nya MiniZinc-modeller (~ Constraint programming news and some new MiniZinc models) the target group for this kind of things (especially in swedish) is quite small. Hence this new blog, and in english.

Some links as introduction to what I have done so far:
* My Constraint Programming page

* My MiniZinc page
* My JaCoP page
* My Choco page

The latter three pages contains information about the systems and some models. I regard the MiniZinc page as my main constraint programming page, since MiniZinc is - as time of writing - my favorite system.

If you know swedish, you may also read the Constraint programming category at hakank.blogg.


News


And we start with some new MiniZinc models written this weekend.

Three Rosetta Code programs, just to test the limits of MiniZinc.

* 99_bottles_of_beer.mzn: 99 bottles of beer
* knapsack_problem.mzn: Knapsack problem
* pyramid_of_numbers.mzn: Pyramid of numbers

And then an operations research model.
* sportsScheduling.mzn Sport scheduling, using channeling for symmetry breaking. I didn't found out how to generate the channeling matrix automatically, so a Perl one-line is used instead (contained in the model). It was inspired by the Essence' model sportsScheduling.eprime from the Tailor (Minion) distribution.