Main

January 16, 2011

Some new Answer Set Programming Programs

Here are some new Answer Set Programming (ASP) programs not mentioned before. They are in about the order they where implemented. All of them are also at my Answer Set Programming page.

Other implementations: Almost all of these problem have been implemented earlier in a couple of CP systems, but instead of linking to all of these individually, I have linked to the problem entry at my Common constraint programming problems. That page contains all of the problems that have been implemented in at least two different CP systems (242 as of writing), and now also contains my ASP programs (83 right now). There are 59 problems which has been been implemented in at least 6 systems. (4 problems has been implemented in all 12 systems:(Survo Puzzle, Seseman, Quasigroup Completion, and Coins Grid).

Finding celebrities

Encoding: finding_celebrities.lp

From Uwe Hoffmann's 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.
(The paper includes an implementation in Scala.)

In summary: Find some distinct clique where everyone known a celebrity, and all celebrities only know all other celebrities.

Here's the code:

knows(adam, dan;alice;peter;eva).
knows(dan, adam;alice;peter).
knows(eva, alice;peter).
knows(alice, peter).
knows(peter, alice).

person(X) :- knows(X, _).
num_p(N) :- N = #sum [person(P) ].

% 1) a person is a celebrity if everyone
% knows P
celebrity(C) :-
person(C),
num_p(N),
N-1 { knows(P, C) : person(P) } N-1.

% 2) and the celebrities only know other
% celebrities, i.e.
% C is not a celebrity if he/she
% knows anyone that is not a celebrity)
:- celebrity(C), person(C), not celebrity(P), knows(C, P).


Solution:
celebrity(alice) celebrity(peter)
See other implementations

Building blocks

Encodings: building_blocks.lp
building_blocks2.lp: faster version

From Brown Buffalo logic puzzle collection (with implementations in Prolog): BuildingBlocksClues:
Each of four alphabet blocks has a single letter of the alphabet on each of its six sides. In all, the four blocks contain every letter but Q and Z. By arranging the blocks in various ways, you can spell all of the words listed below. Can you figure out how the letters are arranged on the four blocks?

BAKE ONYX ECHO OVAL

GIRD SMUG JUMP TORN

LUCK VINY LUSH WRAP
This was quite easy to implement. However, here - as well as in some other ASP encodings - I really miss the global constraint alldifferent that is so convenient in constraint programming systems. Here is my first version of ensuring that the letters of a word (L1, L2, L3, L4) should be on a different die. The words are represented as word(b,a,k,e)..
% the letters of each word must be on
% difference dice
diff(L1,L2,L3,L4) :-
    letters(L1;L2;L3;L4),
    d(D1;D2;D3;D4),
    letter(L1,D1), 
    letter(L2,D2), 
    letter(L3,D3), 
    letter(L4,D4),
    % explicitly state the diffs
    D1 != D2, 
    D1 != D3,
    D1 != D4,
    D2 != D3,
    D2 != D4,
    D3 != D4.
:- not diff(L1,L2,L3,L4), 
   word(L1,L2,L3,L4).
Performance: It takes 1 minute and 10 seconds for gringo/clasp to generate all 24 solutions. The grounding takes much of this time: 30 seconds.

The second version, building_blocks2.lp, use another representation, where the words instead are defined as

word(1, b).
word(1, a).
word(1, k).
word(1, e).
% ...


Or shorter as word(1, b;a;k;e)..

Now we can encode the problem as:

words(1..12).
letters(a;b;c;d;e;f;g;h;i;j;k;l;m;n;o;p;r;s;t;u;v;w;x;y).
d(1..4).

1 { letter(Letter, Dice) : d(Dice) } 1 :- letters(Letter).
6 { letter(Letter, Dice) : letters(Letter) } 6 :- d(Dice).

:- words(W),
word(W, L1),
word(W, L2),
L1 < L2,
letter(L1,D1),
letter(L2,D2),
D1==D2.


Performance: This version is much faster than the first approach. Generating all 24 solution now takes 10 seconds (which essentially is the grounding time).
But it's still slow: all the CP implementations generates all solutions well under one second.

See other implementations

Labeled dice

Encoding: labeled_dice.lp

From Jim Orlin's Colored letters, labeled dice: 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.
The encoding is very similar to the building blocks problem (see above)

See other implementations

SEND+MORE=MONEY, variant with fixed M

Encoding: send_more_money_b.lp

In A first look at Answer Set Programming I complained that Clingo was very slow on these kind of alphametic puzzles (the reason was the grounding time). However, when M is fixed to the value 1 it takes just 5 seconds to solve it, in comparison to 45 seconds when M is "free" (i.e. just > 0).

Project Euler problem 1

Encoding: euler1.lp

Project Euler problem 1:
If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23. Find the sum of all the multiples of 3 or 5 below 1000.
Here is the complete model.

#const n=1000.
number(1..n).
mod_test(N) :- number(N), N #mod 3 == 0.
mod_test(N) :- number(N), N #mod 5 == 0.

total(Sum) :- Sum = #sum [mod_test(N) : number(N) : N < n = N].

#hide.
#show total(Sum).


Note that the two mod_test predicates are to be seend as or, i.e. either multiple of 3 or multiple of 5.

Project Euler problem 2

Encoding: euler2.lp

Project Euler problem 2 involving Fibonacci numbers.
Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be: 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ... Find the sum of all the even-valued terms in the sequence which do not exceed four million.
#const n=45. % max value of N
#const limit = 4000000.

number(1..n).

fib(0, 1).
fib(1, 1).
fib(N, X1 + X2) :- number(N), N > 1, fib(N - 1, X1), fib(N - 2, X2).

total(Sum) :- Sum = #sum [fib(N, Res) : Res < limit : Res #mod 2 == 0 = Res].


The predicate total(Sum) contains the answer.

Prime numbers, small excursion in number theory

Encoding: prime_number.lp

This is small excursion in number theory. Unfortunately gringo don't support arbitrary precision, and the grounding takes long time so this is usable only for small(ish) numbers.

Some definition of primes:

A number N is a prime number if there is no number I: 2..N / 2 such that N % I == 0.

This translates quite well in ASP, where the set notation ({ ... } 0) mean that the set should be empty, i.e. the cardinality of the set is 0.

prime(N) :- number(N), N > 1, { number(I) : I > 1 : I < 1+(N #div 2) : N #mod I == 0} 0.

Another way of defining prime numbers: A number N ( > 1) is a prime number if it has no divisors.

divisor(N, I) :- number(I), I > 1, I < 1+(N #div 2), N #mod I == 0, number(N), N > 1.
prime2(N) :- number(N), N > 1, {divisor(N,I) : number(I) } 0.

A Coin problem

Encoding: coins3.lp

From "The ECLiPSe Book" pages 99f and 234 ff (the ECLiPSe encoding is at page 236.).
What is the minimum number of coins that allows one to pay _exactly_ any amount smaller than one Euro? Recall that there are six different euro cents, of denomination 1, 2, 5, 10, 20, 50
This was harder than I first thought, and was on the wrong track at first. The solution I settled on was to define x(Coin, N) so that N is the maximum value for the N in y (the "matrix" for all combinations of values).

#const max_val = 99.
#const max_amount = 5.
coins(1;2;5;10;25;50).
val(1..max_val). % total
amount(0..max_amount). % max amount per coin

% ensure that all number 1..max_val
% are summed in some way
Value #sum [y(Value, Coin, N) : coins(Coin) : N*Coin <= Value : amount(N) = N*Coin] Value :- val(Value).

% x(Coin, N):
% for each coins Coin, N is the maximum value
% in all y's for this Coin
x(Coin, N) :-
coins(Coin),
amount(N),
N = #max[y(Value, Coin, N2) : val(Value) : amount(N2) = N2].

total(Total) :- Total = #sum [x(Coin, N) = N].

% Object: Minimize the total number of coins.
#minimize [x(Value, N) = N].


See other implementations

Set Covering Deployment

Encoding: set_covering_deployment.lp

From Math World SetCoveringDeployment:
Set covering deployment (sometimes written "set-covering deployment" and abbreviated SCDP for "set covering deployment problem") seeks an optimal stationing of troops in a set of regions so that a relatively small number of troop units can control a large geographic region. ReVelle and Rosing (2000) first described this in a study of Emperor Constantine the Great's mobile field army placements to secure the Roman Empire.
Here is my encoding, except for the matrix of neighbour regions (in matrix(Region1, Region2):

#const n=8. % number of cities
cities(1..n).
army(0..1).

% unique indices of the cities
1 { x(City, Army1, Army2) : army(Army1;Army2) } 1 :- cities(City).

% Constraint 1: There is always an army in a city (+ maybe a backup)
% Or rather: Is there is a backup, there must be an an army
:- x(City, Army1, Army2), Army1 < Army2.

% Constraint 2: There should always be an backup army near every city
:- x(City, Army11, Army12), S = #sum[ x(City2, Army21, Army22) : matrix(City, City2) = Army22], Army11 + S < 1.

% show the selected cities with name
selected(City, Name, Army1, Army2) :- x(City, Army1, Army2), Army1 > 0, name(City, Name).

#minimize [x(City, Army1, Army2) = Army1 + Army2].


The solution found is the following:
selected(3,britain,1,1) selected(1,alexandria,1,1)
which means that we should place one army and one reserv in each britain and alexandria.

(There are 6 optimal solutions of this problem.)

See other implementations

Sicherman Dice

Encoding: sicherman_dice.lp

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.
#const n = 6. % number of side per die
#const m = 10. % max integer on a die
#const start = 1. % start value (set to 0 for the extra two solutions)

% standard distribution of 2 dice
dist(2, 1). dist(3, 2). dist(4, 3).
dist(5, 4). dist(6, 5). dist(7, 6).
dist(8, 5). dist(9, 4). dist(10,3).
dist(11,2). dist(12,1).

sides(1..n).
values(start..m).
dists(2..12).
dist_sums(S) :- dist(_, S).

1 { d1(Side, Value) : values(Value) } 1 :- sides(Side).
1 { d2(Side, Value) : values(Value) } 1 :- sides(Side).

% combine the dice so we can #sum over them
comb(S1,V1,S2,V2) :- d1(S1,V1), d2(S2,V2).
1 { comb(S1,V1,S2,V2) : values(V1;V2) } 1 :- sides(S1;S2).

% symmetry breaking
:- d1(S1a,V1a), d1(S1b,V1b), S1a < S1b, V1a > V1b.
:- d2(S2a,V2a), d2(S2b,V2b), S2a < S2b, V2a > V2b.
:- d1(S,V1),d2(S,V2), V1 > V2.

Sum #sum [ comb(S1, V1, S2, V2) : sides(S1;S2) : values(V1;V2) : V1+V2 == K] Sum :- dist(K, Sum).

#hide.
#show d1(Side, Value).
#show d2(Side, Value).


And it shows the two solutions mentioned above.

See other implementations

Set Partition

Encoding: set_partition.lp

Problem formulation taken from Koalog (could not find the exact page now)
This is a partition problem.
Given the set S = {1, 2, ..., n},
it consists in finding two sets A and B such that:
- A U B = S
- |A| = |B|
- sum(A) = sum(B)
- sum_squares(A) = sum_squares(B)
See other implementations

Averbach problem 1.4

Encoding: averbach_1.4.lp

Problem 1.4 from the book "Problem Solving Through Recreational Mathematics" by Averbach and Chein.
Messr Baker, Dyer, Farmer, Glover, and Hosier are seated around a circular table, playing poker. Each gentleman is the namesake of the profession of one of the others.

The dyer is seated two places to the left of Mr Hosier.
The baker sits two places to Mr Baker's right.
The farmer is seated to the left of Mr Farmer.
Mr Dyer is on the glover's right.

What is the name of the dyer?
This is one of my standard problem for testing the modulo operator (#mod in Clingo). The constraint The dyer is seated two places to the left of Mr Hosier. (which may be stated as dyer = (Hosier - 2) mod 5 in, say, MiniZinc) is translated into:

:- x(_, p_dyer, P1), x(n_hosier, _, P2), P1 != (P2 - 2) #mod 5.

The integrity constraint means that we remove all the solutions that contains the stated condition (hence the !=). Sometimes it's confusing to state the constraint in this manner; on the other hand sometimes it's very convenient.

See other implementations

Scheduling speakers

Encoding: scheduling_speakers.lp

Scheduling of 6 speakers in 6 slots. From Rina Dechter, Constraint Processing, page 72.

See other implementations

Traffic lights problem, CSPLib #16

Encoding: traffic_lights.lp

This is a small problem from the CSPLib where the object is to find different combinations of lights for cars and pedestrians.

See other implementations

Daily schedule, small scheduling problem

Encoding: daily_schedule.lp

This is a small problem I found the other day at the site for the logic programming system Hansei (in ML):
Solving the logical puzzle (scheduling problem) posed by Mikael More on October 25, 2010:

For the daily schedule for Monday to Wednesday:
One of the days I'll shop.
One of the days I'll take a walk.
One of the days I'll go to the barber.
One of the days I'll go to the supermarket.
The same day as I go to the supermarket, I'll shop.
The same day as I take a walk I'll go to the barber.
I'll go to the supermarket the day before the day I'll take a walk.
I'll take a walk Tuesday.

Which are all possible daily schedules, regarding those four events?
...

That is, on Monday I go to the supermarket and shop, on Tuesday I walk and take a haircut. There is only one schedule satisfying the constraints.
It is as direct as possible:

days(monday;tuesday;wednesday).
actions(shop;walk;barber;supermarket).
% an action is done exactly once
1 { x(Day, Action) : days(Day) } 1 :- actions(Action).
% The same day as I go to the supermarket, I'll shop.
:- x(Day1, supermarket), x(Day2, shop), Day1 != Day2.
% The same day as I take a walk I'll go to the barber.
:- x(Day1, walk), x(Day2, barber), Day1 != Day2.
% I'll go to the supermarket the day before the day I'll take a walk.
:- x(Day1, supermarket), x(Day2, walk), Day1 >= Day2.
% I'll take a walk Tuesday.
x(tuesday, walk).

Zebra problem

Encoding: zebra.lp

From Wikipedia: Zebra_Puzzle
1. There are five houses.
2. The Englishman lives in the red house.
3. The Spaniard owns the dog.
4. Coffee is drunk in the green house.
5. The Ukrainian drinks tea.
6. The green house is immediately to the right of the ivory house.
7. The Old Gold smoker owns snails.
8. Kools are smoked in the yellow house.
9. Milk is drunk in the middle house.
10. The Norwegian lives in the first house.
11. The man who smokes Chesterfields lives in the house next to the man with the fox.
12. Kools are smoked in the house next to the house where the horse
is kept. (should be ".. a house ...", see Discussion section)
13. The Lucky Strike smoker drinks orange juice.
14. The Japanese smokes Parliaments.
15. The Norwegian lives next to the blue house.

Now, who drinks water? Who owns the zebra? In the interest of clarity, it must be added that each of the five houses is painted a different color, and their inhabitants are of different national extractions, own different pets, drink different beverages and smoke different brands of American cigarets [sic]. One other thing: in statement 6, right means your right.
The complete code:

% domains
colors(red;green;ivory;yellow;blue).
nationalities(english;spaniard;ukrainian;norwegian;japanese).
animals(dog;fox;horse;zebra;snails).
drinks(coffee;tea;milk;orange_juice;water).
cigarettes(old_gold;kools;chesterfields;lucky_strike;parliaments).

% 1. There are five houses.
houses(1..5).

% alldifferents
1 { color(House, Color) : colors(Color) } 1 :- houses(House).
1 { color(House, Color) : houses(House) } 1 :- colors(Color).

1 { nationality(House, Nationality) : nationalities(Nationality) } 1 :- houses(House).
1 { nationality(House, Nationality) : houses(House) } 1 :- nationalities(Nationality).

1 { animal(House, Animal) : animals(Animal) } 1 :- houses(House).
1 { animal(House, Animal) : houses(House) } 1 :- animals(Animal).

1 { drink(House, Drink) : drinks(Drink) } 1 :- houses(House).
1 { drink(House, Drink) : houses(House) } 1 :- drinks(Drink).

1 { smoke(House, Cigarette) : cigarettes(Cigarette) } 1 :- houses(House).
1 { smoke(House, Cigarette) : houses(House) } 1 :- cigarettes(Cigarette).

next_to(H1,H2) :- houses(H1;H2), |H1-H2| == 1.

% 2. The Englishman lives in the red house.
:- color(H1, red), nationality(H2, english), H1 != H2.

% 3. The Spaniard owns the dog.
:- nationality(H1, spaniard), animal(H2,dog), H1 != H2.

% 4. Coffee is drunk in the green house.
:- color(H1, green), drink(H2, coffee), H1 != H2.

% 5. The Ukrainian drinks tea.
:- nationality(H1, ukrainian), drink(H2, tea), H1 != H2.

% 6. The green house is immediately to the right of the ivory house.
:- color(H1, green), color(H2, ivory), H1 != H2 + 1.

% 7. The Old Gold smoker owns snails.
:- smoke(H1,old_gold), animal(H2, snails), H1 != H2.

% 8. Kools are smoked in the yellow house.
:- smoke(H1, kools), color(H2, yellow), H1 != H2
.
% 9. Milk is drunk in the middle house.
:- not drink(3, milk).

% 10. The Norwegian lives in the first house.
:- not nationality(1,norwegian).

% 11. The man who smokes Chesterfields lives in the house next to the man with the fox.
:- smoke(H1, chesterfields), animal(H2, fox), not next_to(H1,H2).

% 12. Kools are smoked in the house next to the house where the horse
% is kept. (should be ".. a house ...", see Discussion section)
:- smoke(H1, kools), animal(H2, horse), not next_to(H1,H2).

% 13. The Lucky Strike smoker drinks orange juice.
:- smoke(H1, lucky_strike), drink(H2, orange_juice), H1 != H2.

% 14. The Japanese smokes Parliaments.
:- nationality(H1, japanese), smoke(H2, parliaments), H1 != H2.

% 15. The Norwegian lives next to the blue house.
:- nationality(H1, norwegian), color(H2, blue), not next_to(H1, H2).

has_zebra(Nationality) :- nationality(House, Nationality), animal(House, zebra).
drinks_water(Nationality) :- nationality(House, Nationality), drink(House,water).

% for output:
house(House, Color, Nationality, Animal, Drink, Cigarette) :-
houses(House),
color(House, Color),
nationality(House, Nationality),
animal(House, Animal),
drink(House, Drink),
smoke(House, Cigarette).


See other implementations

KenKen

Encoding: kenken2.lp

Wikipedia: KenKen
KenKen or KEN-KEN is a style of arithmetic and logical puzzle sharing several characteristics with sudoku. The name comes from Japanese and is translated as "square wisdom" or "cleverness squared".
...
The objective is to fill the grid in with the digits 1 through 6 such that:

* Each row contains exactly one of each digit
* Each column contains exactly one of each digit
* Each bold-outlined group of cells is a cage containing digits which
achieve the specified result using the specified mathematical operation:
addition (+),
subtraction (-),
multiplication (x),
and division (/).
(Unlike in Killer sudoku, digits may repeat within a group.)

...
More complex KenKen problems are formed using the principles described above but omitting the symbols +, -, x and/, thus leaving them as yet another unknown to be determined.
It took a while to get both principal representation and then the conditions correct. The hints are represented with the operators explicitly stated:

p("+", 11, 1,1, 2,1).
p("/", 2, 1,2, 1,3).
...


Here are some of constraints for handling 2 parameter hints. % 2 arguments
:- p("+", Res, X1,X2, Y1,Y2), x(X1,X2, A), x(Y1,Y2, B), A+B != Res.
:- p("*", Res, X1,X2, Y1,Y2), x(X1,X2, A), x(Y1,Y2, B), A*B != Res.
...


See other implementations

17x17 grid problem

Encoding: 17.lp

I wrote about this grid problem last August in Nontransitive dice, Ormat game, 17x17 challenge.

This ASP encoding solves 14x14 quite fast (< 1 second) using the following parameters:

clingo -c num_rows=14 -c num_cols=14 --stat --heuristic=Vsids --restarts=2000,1.5 --shuffle=1,1 17.lp

with these statistics:
Time        : 0.470
  Prepare   : 0.220
  Prepro.   : 0.050
  Solving   : 0.200
Choices     : 5276
Conflicts   : 3881
Restarts    : 1
However, for larger values (n = 15..16) and experimenting with different solve parameters, I have had no luck.

Kakuro

Encoding: kakuro.lp

Kakuro is yet another grid puzzle. From Wikipedia: Kakuro:
The object of the puzzle is to insert a digit from 1 to 9 inclusive into each white cell such that the sum of the numbers in each entry matches the clue associated with it and that no digit is duplicated in any entry. It is that lack of duplication that makes creating Kakuro puzzles with unique solutions possible, and which means solving a Kakuro puzzle involves investigating combinations more, compared to Sudoku in which the focus is on permutations. There is an unwritten rule for making Kakuro puzzles that each clue must have at least two numbers that add up to it. This is because including one number is mathematically trivial when solving Kakuro puzzles; one can simply disregard the number entirely and subtract it from the clue it indicates.
The encoding is quite fast, solves the problem instance in 0.090s. However, it isn't very general since it only supports clues of length 2..5. Also, it's a tedious job stating the alldifferent constraint for the different arguments.

See other implementations

Young Tableaux and partition

Encoding: young_tableaux.lp

See Wikipedia: Young_tableau:
A Young diagram (also called Ferrers diagram, particularly when represented using dots) is a finite collection of boxes, or cells, arranged in left-justified rows, with the row lengths weakly decreasing (each row has the same or shorter length than its predecessor).
Also see MathWorld: YoungTableau.

The five partitions of 4 are

{4}, {3,1}, {2,2}, {2,1,1}, {1,1,1,1}

The corresponding standard Young tableaux are:
1.   1 2 3 4

2.   1 2 3         1 2 4    1 3 4
     4             3        2

3.   1 2           1 3
     3 4           2 4

4    1 2           1 3      1 4 
     3             2        2 
     4             4        3

5.   1
     2
     3
     4
I thought this should be quite complicated to implement my standard approach of this in ASP. But actually it was quite easy to translate from the MiniZinc model:

See other implementations

Airline crew allocation

Encoding: crew.lp

From the Gecode example crew:
Example: Airline crew allocation Assign 20 flight attendants to 10 flights. Each flight needs a certain number of cabin crew, and they have to speak certain languages. Every cabin crew member has two flights off after an attended flight.
This encoding solves the problem very fast for a schedule all 20 flight attendants, which is the object. However, if we regard it as an optimization problem where the object is to minimize the number of person filling the schedule (i.e. 19), it takes "forever". It finds 19 very fast, but to prove it takes very long (as do many CP system).

See other implementations

Pi Day Sudoku 2009

Encodings: sudoku_pi.lp
sudoku_pi_gcc.lp: with "occurrence count"

Via 360's: Pi Day Sudoku 2009:
Each row, column, and region contains the digits 1-9 exactly once plus three [Pi] symbols. There's a printable .pdf file here
Also see BrainFrees Puzzles: Pi Day 2009.
This was quite easy to encode but very time consuming to convert all the 12 regions to predicates. The cardinality constraints ("exact 3 occurrences of Pi") is encoded quite directly in ASP.

#const n=12.

% domains
rows(1..n).
cols(1..n).
val(1;2;3;4;5;6;7;8;9;p).
val1_9(1..9).
regions(Region) :- region(Region, _, _).

% unique index
1 { x(Row, Col, Val) : val(Val) } 1 :- rows(Row), cols(Col).

% assign hints
x(Row, Col, Val) :- hint(Row, Col, Val).

% all different 1..9
1 { x(Row, Col, Val) : rows(Row) } 1 :- val1_9(Val), cols(Col).
1 { x(Row, Col, Val) : cols(Col) } 1 :- val1_9(Val), rows(Row).

% 3 occurrences of pi per row and column
3 { x(Row, Col, p) : cols(Col) } 3 :- rows(Row).
3 { x(Row, Col, p) : rows(Row) } 3 :- cols(Col).

% handle the regions
% all different for 1..9 in regions
1 { x(Row, Col, Val) : region(Region, Row, Col) } 1 :- regions(Region), val1_9(Val).

% 3 occurrences of pi in the regions
3 { x(Row, Col, p) : region(Region, Row, Col) } 3 :- regions(Region).


However, I miss the generality of the global constraint global cardinality count, and did implement (in sudoku_pi_gcc.lp) some sort of ersatz which handles both the 1 and 3 counts (by the value Occ for each value occ(Val, Occ) instead of stating them separately:

% occ(Number, Occurrences).
occ(1;2;3;4;5;6;7;8;9, 1).
occ(p, 3).

% ..
.
% rows and columns
Occ { x(Row, Col, Val) : rows(Row) } Occ :- val1_9(Val), cols(Col), occ(Val, Occ).
Occ { x(Row, Col, Val) : cols(Col) } Occ :- val1_9(Val), rows(Row), occ(Val, Occ).
% handle the regions
Occ { x(Row, Col, Val) : region(Region, Row, Col) } Occ :- regions(Region), val1_9(Val), occ(Val, Occ).


Both versions solve the problem in about 0.04 seconds but have slightly different number of choices and conflicts:

sudoku_pi.lp: 7 choices, 2 conflicts, 0 restarts
sudoku_pi_gcc.lp: 39 choices, 12 conflicts, 0 restarts

Steiner triplets

Encoding: steiner.lp

This implements the Steiner triplet problem. Formulation from BProlog
The ternary Steiner problem of order n is to find n(n-1)/6 sets of elements in {1,2,...,n} such that each set contains three elements and any two sets have at most one element in common.

For example, the following shows a solution for size n=7:

{1,2,3}, {1,4,5}, {1,6,7}, {2,4,6}, {2,5,7}, {3,4,7}, {3,5,6}

Problem taken from:
C. Gervet: Interval Propagation to Reason about Sets: Definition and Implementation of a PracticalLanguage, Constraints, An International Journal, vol.1, pp.191-246, 1997.
Also see MathWorld SteinerTripleSystem.

I'm not really happy with this encoding. It could probably be done in a more neater way since modeling set constraints is quite natural in ASP. Also, I have not found a way of encode a good symmetry breaking to ensure that the sets are in increasing order.

#const n = 7.
#const nb = n * (n-1) #div 6.
% domain
val(1..n).
sets(1..nb).

0 { x(Set, Val) } 1 :- sets(Set), val(Val).

% 3 values of each set
3 { x(Set, Val) : val(Val) } 3 :- sets(Set).

% count the number of common occurrences of each pair of Sets
check(Set1, Set2, Val, N) :-
sets(Set1;Set2),
Set1 < Set2,
val(Val),
N = #count { x(Set1,Val), x(Set2,Val) }.

% ensure that there are at most 1 occurrence
% with the same value (i.e. N=2) for any two sets.
:- sets(Set1;Set2), Set1 < Set2, 2 #sum[check(Set1, Set2, Val, N) : N == 2 ].


Here are some statistics for clingo for generating 1 solution with default parameters:
n    time   Choices Conflicts Restarts
--------------------------------------
 7    0.01s    183     116      1
 9    0.02s    412     136      1
13    0.33s   6269    2385      6
15    2.35s  22801   13755     10
19    > 2h
With --restarts=10,2.5 it's a bit faster, but not much.
 7    0.01s     57     20      1
 9    0.03s    833    329      4
13    0.32s   6516    2387     6
15    0.46s   6744    1134     5
19    ?


See other implementations

January 09, 2011

Rogo grid puzzle in Answer Set Programming (Clingo) and MiniZinc

Encoding/models:
ASP (Clingo): rogo.lp, rogo2.lp
MiniZinc: rogo.mzn, rogo2.mzn
(See below for some problem instances.)

Later update: I have now also implemented versions with symmetry breaking constraints in the two encodings: rogo2.lp and rogo2.mzn. See below for more into

In Operations Research, Sudoko, Rogo, and Puzzles Mike Trick presented the Rogo puzzle, a grid puzzle where the object is to find a loop of a number of steps and get as many points as possible. He writes
The loop must be a real loop: it must return where it started and canít cross itself or double back. Steps can be either horizontal or vertical: they cannot be diagonal. The loop cannot include any of the black squares. ... Rogo was created by two faculty members (Nicola Petty and Shane Dye)  at the University of Canterbury in Christchurch, New Zealand.  Not surprisingly, Nicola and Shane teach management science there:  the problem has a very strong operations research flavor.
From Creative Heuristics Ltd (the Rogo puzzle site): Rogo is an entirely new type of puzzle. The object is to collect the biggest score possible using a given number of steps in a loop around a grid. The best possible score for a puzzle is given with it, so you can easily check that you have solved the puzzle. Rogo puzzles can also include forbidden squares, which must be avoided in your loop..

Below I have assumed that the path must be in exactly the given number of steps and the programs are written with this assumption in mind (and further reading of the problem descriptions support this interpretation). However, my first approach - probably caused by sloppy reading - was that the optimal path could possible be in lesser steps. It took several hours trying to implement a program supporting this, but was abandoned after re-reading of the problem descriptions. This seems to be a harder problem; maybe it could be used as a variant of the original problem?

I was inspired to solve these puzzle in part because of Mike's last words: Creating a solver would make a nice undergraduate project (and I suspect there are at least a few master's theses and perhaps a doctoral dissertation on algorithmic aspects of creating and solving these). One other reason was to see how to do this with Answer Set Programming - here using Clingo (a Potassco tool) and to compare it with a Constraint Programming system, MiniZinc.

Some more links about Rogo:
  • Instructions
  • Rogo blog
  • YouTube clip.
  • Nicola Petty, Shane Dye: Determining Degree Of Difficulty In Rogo, A TSP-based Paper Puzzle (PDF)
    From the Conclusions: The Rogo puzzle format has a number of aspects that can be controlled to potentially affect degree of difficulty of solving. As a pilot, this study showed that there are many aspects of puzzle-solving related to the nature of the puzzle that can be explored, and there appear to be some general effects, though there are still marked individual differences between people solving the puzzles. This research has the potential to provide interesting insights into both human behaviour, and the nature of puzzles.

    Note: I didn't noticed this until I was almost finished with this blog post (and have just glanced through it).

Example

Here is the Rogo example from Mike Trick's site (pictures borrowed from his site; click on them for larger versions).

Problem:

Rogo puzzle, problem.

One solution:

Rogo puzzle, solution

Note that there is not an unique solution to these puzzles. All three problem instances I tested have more than one solution. For example, the Mike Trick problem has 48 solutions including path symmetries. Since there are 12 steps there are (removing the path symmetry) 48 / 12 = 4 distinct different paths. These different paths are shown below as MiniZinc solutions, where the first step have been fixed to (2,2), i.e. x[1]=2 and y[1]=2<, and it has 3 points (points[1]):

points = array1d(1..12, [3, 0, 0, 2, 0, 0, 2, 0, 0, 1, 0, 0]);
x = array1d(1..12, [2, 3, 4, 4, 5, 5, 5, 4, 3, 2, 2, 2]);
y = array1d(1..12, [2, 2, 2, 3, 3, 4, 5, 5, 5, 5, 4, 3]);
----------
points = array1d(1..12, [3, 0, 0, 2, 0, 0, 2, 0, 0, 1, 0, 0]);
sum_points = 8;
x = array1d(1..12, [2, 3, 3, 4, 5, 5, 5, 4, 3, 2, 2, 2]);
y = array1d(1..12, [2, 2, 3, 3, 3, 4, 5, 5, 5, 5, 4, 3]);
----------
points = array1d(1..12, [3, 0, 0, 1, 0, 0, 2, 0, 0, 2, 0, 0]);
sum_points = 8;
x = array1d(1..12, [2, 2, 2, 2, 3, 4, 5, 5, 5, 4, 3, 3]);
y = array1d(1..12, [2, 3, 4, 5, 5, 5, 5, 4, 3, 3, 3, 2]);
----------
points = array1d(1..12, [3, 0, 0, 1, 0, 0, 2, 0, 0, 2, 0, 0]);
sum_points = 8;
x = array1d(1..12, [2, 2, 2, 2, 3, 4, 5, 5, 5, 4, 4, 3]);
y = array1d(1..12, [2, 3, 4, 5, 5, 5, 5, 4, 3, 3, 2, 2]);
----------

Answer Set Programming, Clingo

Here is the full ASP (Clingo) encoding, without the data:

% domains
rows(1..r).
cols(1..c).

% max number of steps
steps(1..max_steps).

% define adjacency between cells
adj(R,C, R1,C1) :- rows(R;R1), cols(C;C1), |R-R1| + |C-C1|==1.

% the path: unique index
0 { path(I, Row, Col) : steps(I) } 1 :- rows(Row), cols(Col).
1 { path(I, Row, Col) : rows(Row) : cols(Col) } 1 :- steps(I).

% close the circuit: ensure that the first and last cells
% in the path are connected.
:- path(1, R1, C1), path(max_steps, R2, C2), not adj(R1,C1,R2,C2).

% remove bad paths
:- path(I-1,R1,C1), path(I,R2,C2), not adj(R1,C1, R2,C2).

% no black cells in the path
:- path(I, R,C), black(R,C).

% total points, needed since
% "Optimization:" don't show the proper value.
total(Total) :- Total = #sum[got_points(R,C,Value) = Value].

% list the cells in path with points
got_points(R,C, Value) :- point(R,C,Value), path(I, R, C).

% maximize the number of points
% #maximize [ path(I,R,C) : steps(I) : point(R,C,P) = P ].

% alternative: we can add an second objective to
% start with the cell with lowest indices
#maximize [ path(I,R,C) : steps(I) : point(R,C,P) = P@2 ].
#minimize [ path(1,R,C) = R*c+C@1].

#hide.
#show path(I, Row, Col).
#show total(Total).
#show got_points(R,C,Value).

Here is the encoding for Mike Trick's problem instance:
#const max_steps = 12.
#const r = 5.
#const c = 9.

% the point cells
point(1,1,2).
point(2,2,3).
point(2,5,1).
point(2,8,2).
point(3,9,3).
point(4,3,2).
point(5,5,2).
point(5,8,1).

% black cells (to avoid)
black(3,7).
black(4,4).
The solution of Mike Trick's problem (edited), using the following command line:
clingo --heuristic=Vmtf --stat rogo_data_mike_trick.lp rogo.lp
total: 8
got_points:
  5,5,2
  4,3,2
  2,5,1
  2,2,3

path:
  1,2,2
  2,2,3
  3,2,4
  4,2,5
  5,3,5
  6,4,5
  7,5,5
  8,5,4
  9,5,3
  10,4,3
  11,4,2
  12,3,2
Statistics for this solution:
Models      : 1     
  Enumerated: 6
  Optimum   : yes
Optimization: 184 20 
Time        : 0.960
  Prepare   : 0.060
  Prepro.   : 0.020
  Solving   : 0.880
Choices     : 19826
Conflicts   : 16539
Restarts    : 1

Atoms       : 912   
Bodies      : 22839 
Tight       : Yes

  Deleted   : 10406 
Update With the following symmetry breaking added, the problem is solved in 0.58 seconds.

% symmetry breaking: the cell with the lowest coordinates
% should be in the first step
:- path(1, R, C), steps(Step), Step > 1, path(Step, R2, C2),
R*c+C > R2*c+C2.


The statistics for this variant:

Time        : 0.580
  Prepare   : 0.080
  Prepro.   : 0.030
  Solving   : 0.470
Choices     : 8727
Conflicts   : 6914
Restarts    : 2
End of update

Some notes:
One nice feature in Clingo (and lparse) is that it is possible to have many optimization objective. Here we we first maximize the number of points (#maximize [ path(I,R,C) : steps(I) : point(R,C,P) = P@2 ] , and as a second objective (with lower priority @1) we minimize the start cell to start with the cell with the lowest coordinate: #minimize [ path(1,R,C) = R*c+C@1].. Sometimes this is faster than the plain #maximize objective, sometimes not.

The size of "core" of the encoding is quite small. Here is the code with the comments and the helper predicates (for outputs) removed.

rows(1..r).
cols(1..c).
steps(1..max_steps).
adj(R,C, R1,C1) :- rows(R;R1), cols(C;C1), |R-R1| + |C-C1|==1.
0 { path(I, Row, Col) : steps(I) } 1 :- rows(Row), cols(Col).
1 { path(I, Row, Col) : rows(Row) : cols(Col) } 1 :- steps(I).
:- path(1, R1, C1), path(max_steps, R2, C2), not adj(R1,C1,R2,C2).
:- path(I-1,R1,C1), path(I,R2,C2), not adj(R1,C1, R2,C2).
:- path(I, R,C), black(R,C).
#maximize [ path(I,R,C) : steps(I) : point(R,C,P) = P ].

The corresponding "core" of the MiniZinc program (see below for the full code) is larger.

Constraint Programming, MiniZinc

Here is the full MiniZinc code (without data). Compared with the ASP approach, the decision variables are represented in another way: the paths is represented by two arrays x and y, and the points are collected in a separate array (points) so we can simply sum over it for the optimization.

include "globals.mzn";

int: W = 0; % white (empty) cells
int: B = -1; % black cells
int: max_val = max([problem[i,j] | i in 1..rows, j in 1..cols]);

% define the problem
int: rows;
int: cols;
int: max_steps; % max length of the loop
array[1..rows, 1..cols] of int: problem;

% the coordinates in the path
array[1..max_steps] of var 1..rows: x :: is_output;
array[1..max_steps] of var 1..cols: y :: is_output;

% the collected points
int: max_point = max([problem[i,j] | i in 1..rows, j in 1..cols]);
array[1..max_steps] of var 0..max_point : points :: is_output;

% objective: sum of points in the path
int: max_sum = sum([problem[i,j] | i in 1..rows, j in 1..cols where problem[i,j] > 0]);
var 0..max_sum: sum_points :: is_output;

% solve satisfy;
solve maximize sum_points;
% solve :: int_search(x ++ y, first_fail, indomain_min, complete) maximize sum_points;

% all coordinates must be unique
constraint forall(s in 1..max_steps, t in s+1..max_steps) (
x[s] != x[t] \/ y[s] != y[t]
);

% calculate the points (to maximize)
constraint forall(s in 1..max_steps) (
points[s] = problem[x[s], y[s]]
)
/\
sum_points = sum(points);

% ensure that there are no black cells
% in the path
constraint forall(s in 1..max_steps) (
problem[x[s],y[s]] != B
);

% get the path
constraint forall(s in 1..max_steps-1) (
abs(x[s] - x[s+1]) + abs(y[s] - y[s+1]) = 1
)
/\ % close the path around the corner
abs(x[max_steps] - x[1]) + abs(y[max_steps] - y[1]) = 1;


Except for more declaration of the arrays and decision variables, this code don't have more logic than the ASP encoding. However it is more verbose.

The solution for Mike Trick's problem, using LazyFD, takes 1.1 second, slightly slower than using Clingo (see below for more comparison of times):

points = array1d(1..12, [0, 0, 3, 0, 0, 2, 0, 0, 2, 0, 0, 1]);
sum_points = 8;
x = array1d(1..12, [2, 2, 2, 3, 3, 4, 5, 5, 5, 4, 3, 2]);
y = array1d(1..12, [4, 3, 2, 2, 3, 3, 3, 4, 5, 5, 5, 5]);


Update
After some thoughts I decided to try the same symmetry breaking that was an option the ASP encoding. It is implemented in rogo2.mzn and use the following extra constraint that ensure that the cell with the lowest coordinate is in the first step. % symmetry breaking: the cell with lowest coordinates % should be in the first step. constraint forall(i in 2..max_steps) ( x[1]*cols+y[1] < x[i]*cols+y[i] );

With this model, LazyFD solves Mike Trick's problem in 0.627 seconds. Also see under "Comparison" below.

End of update

Comparison

I was curious how well the systems should do so here is a "recreational" comparison. Please don't read too much into it:
  • it is just 3 problems.
  • there is probably better heuristics for both Clingo and Gecode/fz
The following 3 problem instances was used in the test. Unfortunately, I have not found any direct links for the two latter instances (see below for links to my encodings). Instead a variant of the coding used in MiniZinc is shown, where "_" ("W" in the MiniZinc code) is white/blank, "B" is a black cell to avoid, and a number represent a point of the cell.
  • Mike Trick's example. 5 rows, 9 column, 12 steps; good: 6, best: 8.
    %1 2 3 4 5 6 7 8 9  
     2,_,_,_,_,_,_,_,_, % 1
     _,3,_,_,1,_,_,2,_, % 2
     _,_,_,_,_,_,B,_,2, % 3
     _,_,2,B,_,_,_,_,_, % 4
     _,_,_,_,2,_,_,1,_, % 5
    
  • The Paper Rogo puzzle from Creative Heuristics Ltd for 20110106. 9 rows, 7 columns, 16 steps; good: 28, best: 31 points.
     %1 2 3 4 5 6 7
      B,_,6,_,_,3,B, % 1
      2,_,3,_,_,6,_, % 2
      6,_,_,2,_,_,2, % 3
      _,3,_,_,B,B,B, % 4
      _,_,_,2,_,2,B, % 5
      _,_,_,3,_,_,_, % 6
      6,_,6,B,_,_,3, % 7
      3,_,_,_,_,_,6, % 8
      B,2,_,6,_,2,B, % 9
    
  • The Paper Rogo puzzle from Creative Heuristics Ltd for 20110107. 12 rows, 7 columns, 16 steps; good: 34 points, best: 36 points.
     %1 2 3 4 5 6 7
      4,7,_,_,_,_,3, % 1
      _,_,_,_,3,_,4, % 2
      _,_,4,_,7,_,_, % 3
      7,_,3,_,_,_,_, % 4
      B,B,B,_,3,_,_, % 5
      B,B,_,7,_,_,7, % 6
      B,B,_,_,_,4,B, % 7
      B,4,4,_,_,_,B, % 8
      B,_,_,_,_,3,B, % 9
      _,_,3,_,4,B,B, % 10
      3,_,_,_,_,B,B, % 11
      7,_,7,4,B,B,B  % 12
    
For the ASP encoding I used clingo (a combined grounder and solver) and the parameter --heuristi=Vmtf after some minimal testing with different parameters. For MiniZinc, both solvers LazyFD and Gecode/fz where used and I settled with the search heuristic solve minimize sum_points which seems to be fast enough for this experiment. Note that LazyFD tend to behave best without any explicit search heuristics . (Also: there is always a better search heuristics that the one you settle with.)

The time reported is the total time, i.e. including the grounding time/convert to FlatZinc.

Update I have added the second version of the MiniZinc model with the added symmetry breaking constraint as a separate entry: End of update

Mike Trick problem

Clingo: 0.96 seconds, 19826 choices, 16539 conflicts, 1 restart.
Clingo with symmetry breaking: 58 seconds, 8727 choices, 6914 conflicts, 2 restarts.
LazyFD: 1.1 second. (failure is not reported)
LazyFD with symmetry breaking: 0.6 seconds (failure is not reported)
Gecode/fz: 2.62s, 92113 failures2577853 failures
Gecode/fz: with symmetry breaking: 0.4 seconds, 9418 failures

20110106 problem

Clingo: 1:57.07 minutes, 1155290 choices, 1044814 conflicts, 1 restart
Clingo with symmetry breaking: 20.4 seconds, 157146 choices, 135178 conflicts, 3 restarts
LazyFD: 2:58 minutes
LazyFD with symmetry breaking: 19.9 seconds (failure is not reported)
Gecode/fz: 58.6 seconds, 1380512 failures
Gecode/fz with symmetry breaking: 7.8 seconds, failures

20110107 problem

Clingo: 3:13.72 1541808 choices, 1389396 conflicts, 1 restart
Clingo with symmetry breaking: 31.6 seconds 178301 choices, 151439 conflicts, 1 restart
LazyFD: 2:55.18 minutes
LazyFD with symmetry breaking: 44.5 seconds (failure is not reported)
Gecode/fz: 1:54.50 minutes 2577853 failures
Gecode/fz with symmetry breaking: 11.3 seconds, failures

Here we see that Gecode/fz without (symmetry breaking) is the fastest for the two larger problems (but the slowest on the first), and both Clingo and LazyFD was each placed second in the harder problems. So, I'm not really sure we can draw any real conclusions from this small test.

Update With symmetry breaking added the result is more unanimous. All three solvers benefit much from this. Gecode/fz is faster on all three problems, and the other two still one second place each. We also see how much symmetry breaking can do. End of update

Some comments/conclusions

Here are some general comments about the problem and the coding.

ASP

As mentioned above, my initial assumption of the problem was that the given number of steps was not always the best path length, and trying to code this was not easy. After a long time sitting with this approach (say 6 hours coding time?), I re-read the description and couldn't find any real support for this assumption, so I skipped it in favor for the "fixed length" approach.

To get to the final version it took several hours, say 9-10 hours total coding (in smaller sessions over several days). This includes time for coding the initial assumption/interpretation. It also includes the time trying to get the first assumption approach to work with the incremental solver iclingo which - I hoped - should first try the longest length and if that gave no solution it should try at the the lengths in decrement of 1; but I didn't get this to work.

As usual I realized that several of the predicates I first thought was needed could just be thrown away. An example is this snipped to ensure "positive" that there is a path between the cell R1,C1 and the cell R2,C2.

{ x(I-1, R1, C1) : adj(R1,C1, R2,C2) : rows(R1;R2) : cols(C1;C2) : not black(R1,C1) : not black(R2,C2) } max_steps :- x(I, R2, C2), steps(I).

Instead it could be replaced with :- x(I, R,C), black(R,C).
:- x(I-1,R1,C1), x(I,R2,C2), not adj(R1,C1, R2,C2).


Well, I hope that with time I'll recognize these things faster.

MiniZinc

It took about 30 minutes to code the MiniZinc model (after a unnecessarily long detour of debugging bad thinking and a bug in the data). The reason for this much shorter time is two-fold:
  • It's much easier to code something when the general strategy (approach) of the problem is known and has been coded in another system. All the (de)tours when writing the ASP version made much more comfortable with the given problem.
  • I'm much more comfortable coding in MiniZinc than in ASP for two reasons: 1) I have programmed in MiniZinc much longer than in ASP. 2) I am also more used with the high-level CP approach of stating the requirements/constraints with traditional loops and using arrays/matrices.
Programs and data files Here are all the files used, both program and data files for the two systems. ASP (Clingo)
MiniZinc

December 31, 2010

A first look at Answer Set Programming

Almost exactly 2 years ago I started this blog with the following tag line: This is my blog about constraint programming and related paradigms. One of the "related paradigms" I then thought of was Answer Set Programming (ASP). From the Wikipedia article Answer Set Programming:
Answer set programming (ASP) is a form of declarative programming oriented towards difficult (primarily NP-hard) search problems. It is based on the stable model (answer set) semantics of logic programming. In ASP, search problems are reduced to computing stable models, and answer set solvers -- programs for generating stable models are used to perform search. The computational process employed in the design of many answer set solvers is an enhancement of the DPLL algorithm and, in principle, it always terminates (unlike Prolog query evaluation, which may lead to an infinite loop).
In 2001 I first read about Answer Set Programming as a byproduct from an interest in planning and automated theorem proving. I then implemented the ASP program who_killed_agatha.lp, but didn't pursue this interest much more. After that I have - in different times - been very fascinated by this paradigm, especially how succint and clear many of these programs are.

The last weeks I did a more systematic take on ASP and have implemented about 50 different ASP programs (see my Answer Set Programming Page), among them many of the learning problems I use as a stepping stone when learning a new CP system. A full list of the implemented programs and problem instances is shown below.

Answer Set Programming is a huge subject and I will barely scratch the surface here. The object of this blog post is to show how some standard CP problems (and not so standard problems) can be implemented in ASP. I realize that having this approach don't show off some very nice features/applications of ASP, such as the connection to planning, database query, etc. See this as ground for future projects. Also, I have been more concerned about understand how to program in ASP than finding the best heuristics (solver parameters) for each problem. That said, there are some benchmark comparisons below which may not be representative for what ASP systems can do.

My system

All programs was run on this computer: Linux Ubuntu 10.4, 64-bit with 8 cores (Intel i7 930, 2.80GHz) and 12 Gb RAM.

Grounders and solvers

The ASP systems I have worked with solves a problem in two steps. Given a file in apropriate format:
  • first a grounder "grounds" the problem file(s) to an intermittent format. In principle this means that all possible values of the variables are generated including all the combinations that are forbidden by integrity constraints. Note: This can be very large files.
  • then a solver reads this grounded file and - hopefully - solves the problem by generating answers sets (solutions).
After some considerations I decided to go with Potassco (the Potsdam Answer Set Solving Collection ) tools which includes the grounder gringo, the solver clasp and a program that combines these two: clingo, and it is Clingo that I tend to use. The versions I have been using are clasp version 1.3.6, gringo version 3.0.3, clingo version 3.0.3. For some tests with lparse/smodels: lparse version 1.1.2 and smodels 2.34.

Another option was to use the grounder/solver lparse/smodels but I decided to use Potassco since it seems to be much more active and it include a lot of other tools, for example clingcon (using constraint programming domains etc), which I will study more later.

Another reason to use the Potassco tools was the impressive results in The Second Answer Set Programming Competition (for year 2010). For more about this, see the report The Second Answer Set Programming Competition (PDF).

However, with some few problem I also compared with with lparse/smodels, e.g. send_more_money.lp.

Note that I have not tweaked very much with the large amount of flags for the clasp solver. In some problems (e.g. the alphanumeric) I tested some of them but didn't find any that speeded it up considerably.

More information about ASP

Here are some useful links for Answer Set Programming:

ASP Systems

There are many ASP systems. Here are those I have tested most.

Some programs

Here is some examples how one can encode some problems with ASP (or rather in lparse/gringo format). Let's start with Sudoku.

Sudoku

Program: sudoku.lp

Here is one encoding of the problem (inspired by this). The hints part is described below.
% domains
val(1..9).
border(1;4;7).
% alldifferent rows, columns, values
1 { x(X,Y,N) : val(N) } 1 :- val(X;Y).
1 { x(X,Y,N) : val(X) } 1 :- val(N;Y).
1 { x(X,Y,N) : val(Y) } 1 :- val(N;X). 
% alldifferent: boxes
1 { x(X,Y,N) : val(X;Y):
    X1<=X:X<=X1+2:Y1<=Y:Y<=Y1+2 } 1 :- val(N), border(X1;Y1).
One really big difference with coding ASP and CP is that in CP the data is often in matrix/array form (except maybe for the Prolog based CP systems). In ASP there are no arrays or matrices, so instead one have to state the order and "coordinates" by enumerating the combinations. E.g. the following problem instance is normally coded in CP as something like (where "_" must replaced with "0" in some systems).
 _,_,6, _,_,_, _,9,_,
 _,_,_, 5,_,1, 7,_,_,
 2,_,_, 9,_,_, 3,_,_,

 _,7,_, _,3,_, _,5,_,
 _,2,_, _,9,_, _,6,_,
 _,4,_, _,8,_, _,2,_,

 _,_,1, _,_,3, _,_,4,
 _,_,5, 2,_,7, _,_,_,
 _,3,_, _,_,_, 8,_,_
In ASP encodings this is stated as predicates x(Row, Column, Value), e.g.:
x(1, 3, 6).
x(1, 8, 9).
x(2, 4, 5).
x(2, 6, 1).
...
x(8, 6, 7).
x(9, 2, 3).
x(9, 7, 8).
This lack of direct matrix representation (and no direct representation of arrays) has been one of the biggest hurdle when translating the MiniZinc models to ASP. (It takes a lot of time to manually convert the matrix representation to this predicate variant, so I've created a simple Perl program for this: matrix2predicate.pl)

OK, now some explanations of the ASP encoding above.
val and border is "domain predicates" which contains the different values (1..9, and 1,4,7 respectively) to use in the program. They are used to generate different values later on.

Next we have the constraints that all rows, columns, and boxes should contain different values. In Constraint Programming this is (normally) done with a bunch of alldifferent constraints, but in ASP we (normally) have to use another approach. (however there are some systems, e.g. clingcon mentioned above, that is a hybrid of ASP and CP which includes some global constraints.)

Instead one have to use some other features which can be described as.
  • "generators"
  • cardinality constraint
  • integrity constraint
The first line: 1 { x(X,Y,N) : val(N) } 1 :- val(X;Y). is to ensure that the indices in the "matrix" (X and Y) are unique, i.e. that the we have exactly 1 occurrence of x(1,1,Val) (for some value Val), 1 occurence of x(1,2, Val) etc.

The left "1" around the curly braces is the lower bound, and the right "1" is the upper bound of the expression. This is the cardinality constraints (think atleast and atmost in CP). The part after the colon (:) - : val(N) - is a condition filter which here means that N must take the values of the domain val/1 (i.e. the values 1..9).

The right part (the body) val(X;Y) is used as a "generator" of the expression; I tend to think of this as a "for loop" over X and Y that "drives" the left part.

Continuing with the next two lines: 1 { x(X,Y,N) : val(X) } 1 :- val(N;Y). 1 { x(X,Y,N) : val(Y) } 1 :- val(N;X).

These act as alldifferent constraints: for each combination of val(N) and val(Y) (val(X;Y) is a shorthand for val(X) and val(Y)) there must be exactly one occurrence of X (the rows). And the next line is for Y.

The last constraint
1 { x(X,Y,N) : val(X;Y):
    X1<=X:X<=X1+2:Y1<=Y:Y<=Y1+2 } 1 :- val(N), border(X1;Y1).
is for the boxes, and we can here see how the filter conditions (the ":" parts) can be combined. Note that the "loop" is over X1 and Y1.

Note: There is an alternative version of the constraints for the rows, columns, values which seems to be slighly faster:
% alternative:
:- 2 { x(X,Y,N) : val(N) }, val(X;Y).
:- 2 { x(X,Y,N) : val(X) }, val(N;Y).
:- 2 { x(X,Y,N) : val(Y) }, val(N;X). 
These "headless" expressions (starts with :-) are integrity constraints which are used to remove certain results: if the expression is true then it don't belong to a solution. They all states that the occurrences of each value (X;Y) cannot be 2 or larger.

nqueens

Program: nqueens.lp

Let's continue with another standard CP problems and compare the performances between encoding in ASP and some CP systems. Here is a complete encoding for n-queens problem in ASP: nqueens.lp.
#const n = 8.
number(1..n).
% alldifferent
1 { q(X,Y) : number(Y) } 1 :- number(X).
1 { q(X,Y) : number(X) } 1 :- number(Y).
% remove conflicting answers
:- number(X1;X2;Y1;Y2), q(X1,Y1), q(X2,Y2), X1 < X2, Y1 == Y2.
:- number(X1;X2;Y1;Y2), q(X1,Y1), q(X2,Y2), X1 < X2, Y1 + X1 == Y2 + X2.
:- number(X1;X2;Y1;Y2), q(X1,Y1), q(X2,Y2), X1 < X2, Y1 - X1 == Y2 - X2.
The first line (#const n = 8.) is a constant, and can be redefined from the command line with different grounders/solvers:
$ clingo -c n=12 nqueens.lp
$ gringo -c n=12 nqueens.lp | clasp
$ lparse -c n=12 nqueens.lp | smodels
$ lparse -c n=12 nqueens.lp | clasp
The program use the same principles as for the Sudoku encoding. We represent the queens as a predicate q(X, Y) where X is the position in the "array" q, and Y is the value. (Or maybe as a 2-dimensional matrix of row X and column Y, but I tend to think about arrays with a unique index.)
% alldifferent
1 { q(X,Y) : number(Y) } 1 :- number(X).
1 { q(X,Y) : number(X) } 1 :- number(Y).
The first line ensures unicity of the X:s, and the second line ensure the unicity of the Y:s.

After that is the constraints (and integrity constraints) stating that the queens can not be placed on the same line, same row or column.

Since it is a standard benchmark problem in CP, I thought it might be interesting to see how well this ASP program do.

Below are some statistics for generating the first solution for different values of n, with two different grounders (gringo and lparse) and two solvers (clasp and smodels). I had a tiny time limit for these problems, 2 minutes (except for N=200). Here smodels didn't do very well so I skipped it for N larger than 50.

In this table the times for the grounders and the solvers are separared. As we can see there is quite a difference of the times for the two grounders (gringo, lparse) and the solvers (clasp, smodels).
 N   Grounder Time    Solver  Time   Choices Conflicts Restarts
 ----------------------------------------------------------------
  8  gringo   0.006s  clasp    0.004s     5      2      0
  8  gringo   0.006s  smodels  0.008s     5      -      -
  8  lparse   0.023s  clasp    0.02s      4      1      0
  8  lparse   0.023s  smodels  0.02s      3      -      -

 50  gringo   1.14s   clasp    0.424s    366    252     2
 50  gringo   1.14s   smodels  > 2min      -      -     -
 50  lparse   3.34s   clasp    0.37s     211    137     1 
 50  lparse   3.34s   smodels  > 2min      -      -     -

100  gringo  15.5s    clasp    3.81s     839    505     3
100  lparse  33.1s    clasp    3.4s      599    332     2

200! gringo 3:38min   clasp    45.3s    9782   7747     9
200  lparse 6:30min   clasp    38.5s    3669   2347     6
Note: The results above was obtained by the following:
 $ time gringo -c n=200 nqueens.lp > xxx
 $ time clasp xxx
However, it seems that piping the results directly to the solver is faster. The following $ time gringo -c n=200 nqueens.lp | clasp takes less time than running the programs separately: 4:01 minutes (compare with 3:38min + 45s = 4.23min).

Using lparse, $ time lparse -c n=200 nqueens.lp | clasp: 6:28min (compare with 6:30min + 38.5s = 7:08min).

Compare these numbers with the MiniZinc model using the same modelling approach: queens.mzn with Gecode/fz as solver. These times are much better except maybe for small N. Or rather, the ASP solvers seems to be quite competitive, it is the ASP grounders that take long time. MiniZinc uses the same principle: first the MiniZinc model is converted to a FlatZinc file, and then the CP solvers works on the FlatZinc file. The times below includes both these steps.
  N    Time Failures
  ------------------
   8   0.08s     3
  50   0.2s      2
 100   0.5s     60
 200   2.9s      8
 500  21.0s     25 
1000 1:44min     0
  
With the Google CP Solver program nqueens2.py (Python interface) we have even better results (for larger N times also without output):
  N    Time Failures Time without output
  -----------------------------------
  8    0.04s     3
  50   0.08s    16
 100   0.2s     12   0.1s
 200   0.6s      1   0.5s
 500   4.0s     25   2.9s 
1000  16.1s     0   12.1s 
 
Clingcon
As mentioned above, there is a hybrid of ASP and CP (clingcon which seems to be much faster. However, the example of n-queens in the distribution (queens.lp) don't give correct result. For n=8 it should be 92 solutions, but it generates 1066 solutions. There is an alldifferent predicate in the program but it is commented which probably is the cause of this; removing the comment generates a segmentation fault.

magic square

Program: magic_square.lp

Magic squares is another standard CP problem.
#const n = 3.               % the size
#const s = n*(n*n + 1) / 2. % the num
% domains
size(1..n).
val(1..n*n).
% unique index of x
1 { x(Row, Col, N) : val(N) } 1 :- size(Row;Col).
% alldifferent values of x
1 { x(Row, Col, N) : size(Row;Col) } 1 :- val(N).
% sum rows
:- not s #sum[ x(Row, Col, Val) : size(Col) : val(Val) = Val ] s, size(Row).
% sum columns
:- not s #sum[ x(Row, Col, Val) : size(Row) : val(Val) = Val ] s, size(Col).
% Sum diagonals 1 and 2
:- not s #sum[ x(I, I, Val) : size(I) : val(Val) = Val ] s.
:- not s #sum[ x(I, n-I+1, Val) : size(I) : val(Val) = Val ] s.
#hide.
#show x(Row, Col, N).
Note here how the constraints of sums are encoded. We use an integrity constraint :- not s #sum[ x(Row, Col, Val) : size(Col) : val(Val) = Val ] s, size(Row). stating that if the sum is not s, then it should be removed. This is kind of awkward sometimes but it is often required.

Here are some statistics for 1 solution of different n:s. Running with > $ time gringo -c n=N magic_square.lp | clasp --stat
 N  gringo     clasp   Choices   Conflicts Restarts
 --------------------------------------------------
 3   0.004s    0.009s     244       230       1
 4   0.007s    0.041s    1764      1545       5
 5   0.011s    0.061s    1946      1618       5
 6   0.014s    8.509s  184347    172478      16
 7   0.035s    6.007s   78131     72744      14
 8   0.042s   1:56min  925330    858613      20
 9   0.046s  2:32:37h  36674602 32898767     29
In contrast to the n-queens problem the grounding step is quite fast, it is the solver step that takes time. This also means that it may be room for a lot of improvements via different parameter settings of clasp.

Compare this with MiniZinc model (magic_square.mzn) that use about the same approach using Gecode/fz as the solver:
 N  time   failures
 ------------------
 3  0.077s         3
 4  0.076s         7
 5  0.077s       678
 6  0.083s       810
 7  14:52m  99401166
Here the gringo/clasp implementation manage to handle larger N than the MiniZinc+Gecode/fz approach. However the plain Gecode version (magic-square) using standard settings is faster on some of the problem instances. It's quite interesting that the ASP program is faster than Gecode for N=8 (and assuminly also for N=9).
 N  time   failures
 ------------------
 3  0.008s        6
 4  0.013s       892
 5  0.427s     72227
 6  0.009s        27
 7  2.628s     481301
 8  > 1hour
Note: the Gecode mode has some symmetry breaking: x(1,1) > x(1,n) and x(1,1) > x(n,1). When I tried these in the ASP program it seems to be slower, not faster.

At last the mandatory declaimer: It's perhaps misleading to compare times on different systems/approaches like this. But it might give a good feeling of these systems/models.

Papers on Comparisons ASP vs CP

There have been some more systematic comparisons of ASP and CP, which - as I understand - is written by ASP researchers: These two papers are interesting since they show ASP in many areas where it do quite well: graph related problems. Many of the grid puzzles in the second paper use ASP's ability to encode transitive closure (see below) rather easy. This is not so easy (or "natural") to code in most CP systems.

Other programs

Here is another samples of ASP programs together with some code and comments.

Who killed Agatha

Program: who_killed_agatha.lp As mentioned above, this was actually my first program in ASP, written when I first read about ASP, about 2002. (And that's why this problem has been a standard problem since.)

Here we see more of the declarative aspect of ASP and its use of defined predicates (hates/2, richer/2) which are then used in integrity constraints.
% "Someone who lives in Dreadsbury Mansion, killed aunt Agatha. Agatha, 
% the butler, and Charles live in Dreadsbury Mansion, and are the only 
% people who live therein."

% The people in this drama
person(agatha;butler;charles).

% Exactly one person killed agatha
1 {killed(agatha,agatha), killed(butler,agatha), killed(charles,agatha)} 1.

% A killer always hates her victim, and is never richer than her victim
hates(Killer,Victim) :- person(Killer), person(Victim), killed(Killer,Victim).
:- richer(Killer,Victim), person(Killer), person(Victim), killed(Killer,Victim). 

% a person cannot be richer than him/herself
richer(X,Y) :- person(X), person(Y), X != Y.

% Charles hates no one that aunt Agatha hates
0 { hates(charles,X) } 0 :- person(X), hates(agatha,X).

% Agatha hates everyone except the butler, ...
hates(agatha,X) :- person(X), X != butler.
% ... and no one the butler does not hate.
:- hates(agatha,X), person(X), hates(butler,X).

% the butler hates everyone not richer than aunt Agatha
:- hates(butler, X),  person(X), richer(X,agatha).

% No one [in Dreadsbury Mansion] hates everyone
{hates(X,Y):person(Y)} 2 :- person(X).

% Who killed aunt Agatha?
{killed(X,agatha):person(X)}.

#hide.
#show killed(X,Y).

Diet

Program: diet1.lp

This is a simple optmization problem where the object is to minimize the cost of the products given that we must exceed some limits of the ingredients
#const n = 10.
amount(0..n). % max amount of each product

%    food                  calories  chocolate   sugar     fat    price
%                          (ounces)    (ounces) (ounces) (ounces) 
food_a(chocolate_cake,       400,      3,          2,        2,   50).
food_a(chocolate_ice_cream,  200,      2,          2,        4,   20).
food_a(cola,                 150,      0,          4,        1,   30).
food_a(pineapple_cheesecake, 500,      0,          4,        5,   80).

% Minimum limits
limits(calories, 500).
limits(chocolate, 6).
limits(sugar, 10).
limits(fat, 8).

% extract the different nutrition types
calories(Food, Amount)  :- food_a(Food, Amount, _, _, _, _).
chocolate(Food, Amount) :- food_a(Food, _, Amount, _, _, _).
sugar(Food, Amount)     :- food_a(Food, _, _, Amount, _, _).
fat(Food, Amount)       :- food_a(Food, _, _, _, Amount, _).

% extrace the price
price(Food, Price) :- food_a(Food, _, _, _, _, Price).
prices(Price) :- price(Food, Price).

% food(chocolate_cake;chocolate_ice_cream;cola;pineapple_cheesecake).
food(F) :- food_a(F, _,_, _, _, _). 

% each food has exactly one Price and one Amount
1 { food_price_amount(Food,Price,Amount) : amount(Amount) } 1 :- 
                                                          food(Food),
                                                          price(Food, Price).

:- #sum[food_price_amount(F, _, A)  : calories(F, C)  = C*A] L-1, limits(calories, L).
:- #sum[food_price_amount(F, _, A)  : chocolate(F, C) = C*A] L-1, limits(chocolate, L).
:- #sum[food_price_amount(F, _, A)  : sugar(F, C)     = C*A] L-1, limits(sugar, L).
:- #sum[food_price_amount(F, _, A)  : fat(F, C)       = C*A] L-1, limits(fat, L).

#minimize [food_price_amount(F,Price,Amount) = Price*Amount ].
The way I have encoded it is rather complex but most of it is to extract different values from the predicate food_a. I have a feeling that there is a much better way...

A comment on the #sum construct. This code
:- #sum[food_price_amount(F, _, A)  : calories(F, C)  = C*A] L-1, limits(calories, L).
means that the sum of the products of the amount of food (A) and the calories in the food (C) should be larger than the limit (L) of this food. Since it's a integrity constrain we must negate the expression.

The objective is to minimize (#minimize) the cost of all selected products.

For me this problem is easier to encode with a traditional "array/matrix" approach in CP, and it is also easier to generalize (by just adding a dimension or two to the arrays). Compare for example with the MiniZinc model diet1.mzn

However, I like that it is possible to use symbolic domains directly instead of using integers (and enums).

SEND+MORE=MONEY and other alphametic problems

Program: send_more_money.lp
Program: send_more_money_any_base.lp
Program: send_most_money.lp
Program: least_diff.lp

Here is my encoding of SEND+MORE=MONEY (send_more_money.lp). Note the use of the domain predicate letter) to define the "slots" in the "array" x(Letter, Value which are then extracted in the main predicate smm. Also note that we have to use :- not smm in order to remove the combinations we don't want (i.e. remove everything that is not satified by the rules in smm).
letter(s;e;n;d;m;o;r;y).
values(0..9).
% exact 1 occurrence of each letter
1 { x(L,Val) : values(Val) } 1 :- letter(L).
% 0..1 occurrences of each value
{ x(L,Val) : letter(L) } 1 :- values(Val).
% no digit can be given to two different symbols
% :- letter(L), letter(L1), x(L,V1), x(L1,V1), L != L1.
smm :- 
   values(S;E;N;D;M;O;R;Y),
   x(s,S), x(e,E), x(n,N), x(d,D),
   x(m,M), x(o,O), x(r,R), x(y,Y), 
   M > 0,  S > 0,
   S*1000+E*100+N*10+D + M*1000+O*100+R*10+E ==
   M*10000+O*1000+N*100+E*10+Y.
:- not smm.
I was surprised - and somewhat disappointed - that it took so long to solve these kind of alphametic problems with standard ASP tools. It took about 45 seconds to solve this with clingo (gringo + clasp). lparse accept this encoding with just some warnings and it takes 25s.

After a while I realized that it is the grounding (gringo) that takes so long, almost all of the 45 seconds. Given this grounding as a file, then the solvers - clasp and smodel - solves this in no time, about 0.004s.

clingcon:
There is also an example of this problem in clingcon's distribution (sendmoremoney1.lp). It works and is very fast (0.001s). However, I have not fully understood how to encode other similar problems with optimizations, different domains, etc (e.g. SEND + MOST + MONEY) etc, so understanding this is a further project. Also, it seems to have some bugs in the current version, especially regarding global constraints.

send_most_money.lp
This solves the problem SEND+MOST=MONEY and maximizes MONEY. It seems to works since it show the correct maximum value of MONEY (10768). Almost, since the value of Optimization is not this value but a much larger number (4949944124). I read in the Potassco documentation that maximization is converting to minimization with the weights inversed. Maybe the weird Optimization value is a consequence of that?
Answer: 1
x(e,6) x(n,7) x(d,3) x(m,1) x(o,0) x(t,5) x(y,8) x(s,9) y(money,10768) 
Optimization: 4949944232
Answer: 2
x(e,7) x(n,8) x(d,4) x(m,1) x(o,0) x(t,2) x(y,6) x(s,9) y(money,10876) 
Optimization: 4949944124
OPTIMUM FOUND
As with send_more_money.pl it takes a long time to solve it: 2:24 minutes with clingo (gringo/clasp). With lparse/smodels it takes much longer.
gringo/clasp: 2:24 minutes
gringo/smodels: over 7 minutes and then some "strange" result is shown from
                smodels (e.g. not y(money,98979)) and I stopped.
lparse/smodels: over 10 minutes for lparse
lparse/clasp:  over 10 minutes for lparse
send_more_money_any_base.lp: SEND+MORE=MONEY for "any" base
For n=10 it takes the same time as for send_more_money.lp. For n=11 it takes 1:37 minutes to show all 3 solutions.

least_diff.lp
I have also implemented the least diff problem, but I stopped the program after 30 minutes without any answer.

Well, maybe my approach in these problems is not a good example of how to use ASP. However, they are standard examples in CP.

Other standard CP/OR problems

Here are some other standard CP or OR problems.

all_interval.lp: All interval.

langford.lp: Langford numbers
This follows quite nicely the MiniZinc version in langford.mzn
#const k = 4.
val(1..k).   % for solution
pos(1..2*k). % for position
% alldifferent position
1 { position(I, N) : pos(N) } 1 :- pos(I).
1 { position(I, N) : pos(I) } 1 :- pos(N).

% The difference in position between the two I's is I.
:- position(I+k, N1), position(I, N2), N1 != N2 + I+1, val(I).

% solution: unique index
1 { solution(I, N) : val(N) } 1 :- pos(I).

% exactly two occurrences of 1..k
2 { solution(I, N) : pos(I) } 2 :- val(N).

% channel solution <-> position
:- position(I, P1), solution(P1, I2), I != I2, val(I;I2).
:- position(I+k, P2), solution(P2, I2), I != I2, val(I;I2).
alldifferent_except_0.lp: All different except 0
Since there is no concept of global constraints in ASP, one has to code these kind of constraints directly:
#const n = 6.
#const m = 9.
values(0..m).
ix(1..n).
% unique indices of x, 1..n
1 { x(I, Val) : values(Val) } 1 :- ix(I).
% alldifferent except 0
% If Val > 0 then there must be 0..1 
% occurrences of Val in x.
{ x(I, Val) : ix(I) } 1 :- values(Val), Val > 0.

% Additional constraint: There must be exactly 2 zeros.
% 2 #sum [ x(I, 0) : ix(I) = 1 ] 2.
% Alternative:
:- not 2 #sum [ x(I, 0) : ix(I) = 1 ] 2.
bus_scheduling.lp: Simple scheduling problem

organize_day.lp: Another simple scheduling problem: How to organize a day
The no overlap requirement are done with the following. A task is defined as the predicate task(TaskId, StartTime, EndTime). The predicate no_overlap is very much Prolog, apart from that the order of rules and the order of predicates don't matter in ASP.
% No overlap of tasks
no_overlap(Task1, Task2) :-
        task(Task1, Start1, End1),
        task(Task2, Start2, End2),
        End1 <= Start2.

no_overlap(Task1, Task2) :-
        task(Task1, Start1, End1),
        task(Task2, Start2, End2),
        End2 <= Start1.

:- tasks(Task1;Task2), Task1 != Task2, not no_overlap(Task1, Task2).
photo.lp: Photo problem.
This is an example from Mozart/Oz tutorial on FD. The object is to maximize the number of satisfied preferences. I like this direct approach.
persons(betty;chris;donald;fred;gary;mary;paul).
positions(1..7).
% Preferences:
pref(betty, gary;mary).
pref(chris, betty;gary).
pref(fred, mary;donald).
pref(paul, fred;donald).
% alldifferent positions
1 { position(Person, Pos) : positions(Pos) } 1 :- persons(Person).
1 { position(Person, Pos) : persons(Person) } 1 :- positions(Pos).
next_to(P1,P2) :-
        pref(P1,P2),
        position(P1,Pos1),
        position(P2,Pos2),
        |Pos1-Pos2| == 1.
% maximize the number of satisfied preferences
#maximize [ next_to(P1,P2) ].
post_office_problem2.lp: Post office problem

coloring.lp: Simple map coloring problem
countries(belgium;denmark;france;germany;netherlands;luxembourg).
colors(red;green;blue;white).
arc(france,belgium;luxembourg;germany).
arc(luxembourg,germany;belgium).
arc(netherlands,belgium).
arc(germany,belgium;netherlands;denmark).
neighbour(X,Y) :- arc(X,Y).
neighbour(Y,X) :- arc(X,Y).
1 {color(X, C) : colors(C) } 1 :- countries(X). 
:- color(X1, C), color(X2, C), neighbour(X1,X2).
:- color(germany, red).
A variant to the generator line (line 9) is to explicit state the different alternative with disjunction using '|' (bar):
color(X, red) | color(X, green) | color(X, blue) | color(X, white) :- countries(X).
However, then one has to use the parameter --shift to gringo/clingo for this to work.

xkcd.lp: xkcd problem
This is a subset sum (or knapsack depending on the definition) problem from xkcd. The object is to select dishes so they give a total of 15.05 (as 1505 in the encoding).
#const total = 1505.
#const n = 10.
amount(0..n).
food(mixed_fruit;french_fries;side_salad;host_wings;mozzarella_sticks;samples_place).
price(mixed_fruit,215).
price(french_fries,275).
price(side_salad,335).
price(host_wings,355).
price(mozzarella_sticks,420).
price(samples_place,580).
prices(P) :- price(_, P).
% each food has exactly one amount
1 { food_amount(Food, Amount) : amount(Amount) } 1 :- food(Food).
% sum to the exact amount
total [ food_amount(F, Amount) : price(F, Price) : prices(Price) : amount(Amount) = Price*Amount ] total.
One can easily change the total (e.g. to 1506) with clingo -c total=1506 xkcd.lp. Also, here is an (not uncommon) example that the data part is larger than the actual program.

subset_sum.lp: Another subset sum problem

3_jugs.lp: 3 jugs problem (as a graph problem)
Here is a graph approach to the 3 jugs problem which shows another use of an adjacency predicate (adj). Here we also add the edge visited with an index to get the order in the graph traversal. The graph is of the states of the different ways of pouring from the jugs. The object is to get the shortest path from node 1 to node 15. Note that we select the edges and then with the predicate selected_nodes calculates which node was visited at what time.
% g(EdgeId, From, To).
g(1, 1, 2).
g(2, 1, 9).
g(3, 2, 3).
g(4, 3, 4).
g(5, 3, 9).
g(6, 4, 5).
g(7, 5, 6).
g(8, 5, 9).
g(9, 6, 7).
g(10, 7, 8).
g(11, 7, 9).
g(12, 8, 15).
g(13, 9, 10).
g(14, 10, 2).
g(15, 10, 11).
g(16, 11, 12).
g(17, 12, 2).
g(18, 12, 13).
g(19, 13, 14).
g(20, 14, 2).
g(21, 14, 15).

#const start = 1.
#const end = 15.
edges(1..21).
ix(1..21).

% ensure that the index is unique
{ selected(Ix, Edge) : edges(Edge) } 1 :- ix(Ix).
% ensure that we visit an edge atmost once
{ selected(Ix, Edge) : ix(Ix) } 1 :- edges(Edge).

% define adjacency and add the Edge to selected
adj(X, Y, I) :- g(Edge, X, Y), selected(I, Edge).
adj(X, Y, I) :- g(Edge, X, Z), adj(Z, Y, I+1), selected(I, Edge).

% init the problem: from start to end 
% with index 1)
:- not adj(start, end, 1).

% Here we check which nodes that was involved
% in the selected edges.
selected_nodes(Ix, X, Y) :- selected(Ix, Edge), g(Edge, X, Y).

% minimize the number of edges
#minimize [ selected(Ix, Edge) : edges(Edge) : ix(Ix) ].
The solution is (edited). The first number is the order index.
selected edges:
	1,2
	2,13
	3,15
	4,16
	5,18
	6,19
	7,21

selected_nodes
	1,1,9
	2,9,10
	3,10,11
	4,11,12
	5,12,13
	6,13,14
	7,14,15
I.e. the nodes are visited in the following order (as it should): 1,9,10,11,12,13,14,15 via the edges 2,13,15,16,18,19,21.

The second definition of adj:
adj(X, Y, I) :- g(Edge, X, Z), adj(Z, Y, I+1), selected(I, Edge).
could as well be written with g/3 and adj/3 swapped:
adj(X, Y, I) :- adj(Z, Y, I+1), g(Edge, X, Z), selected(I, Edge).
ASP don't care about the order of rules in the predicates, in contrast to Prolog where order can matter very much.

Grid puzzles

There are quite a lot of grid puzzles in my MiniZinc collection and I have implemented some of them in ASP. The problem instances are shown below in the full list of files.

minesweeper.lp: Minesweeper. This encoding was inspired by the "Phase1" encoding from Nitisha Desai's (?) .
nums(0..8).
rows(1..r).
cols(1..c).
%a cell can be a number or a mine
1 {number(R,C,Z) : nums(Z), mine(R, C)} 1 :- rows(R), cols(C).
% defining adjacency
adj(R,C,R1,C1) :- rows(R;R1), cols(C;C1), |R-R1| + |C-C1|==1.
adj(R,C,R1,C1) :- rows(R;R1), cols(C;C1), |R-R1|==1, |C-C1|==1.
% N mines around a number N
N {mine(R2, C2) : adj(R2,C2,R1,C1)} N :- number(R1,C1,N).
quasigroup_completion.lp: Quasigroup completion
This is pure Latin square conditions that we saw from the Sudoku encoding:
values(1..n).
% all different rows, columns, and ensure that we use 
% distinct values
1 { cell(Row, Col, Val) : values(Row) } 1 :- values(Col;Val).
1 { cell(Row, Col, Val) : values(Col) } 1 :- values(Row;Val).
1 { cell(Row, Col, Val) : values(Val) } 1 :- values(Row;Col).
survo puzzle.lp: Survo puzzle.
The lines that handle the row and column sums are the following. Note how we extract the lower/upper bounds of the sum from the predicates rowsum (and colsum).
% Row sums
:- not RowSum #sum [matrix(Row, Col, Val) : values(Val) : cols(Col) = Val] RowSum,
                                                               rows(Row), rowsums(Row, RowSum).
% Col sums
:- not ColSum #sum [matrix(Row, Col, Val) : values(Val) : rows(Row) = Val] ColSum,
                                                              cols(Col), colsums(Col, ColSum).
discrete_tomography.lp: Discrete tomography

hidato.lp: Hidato puzzle

strimko2.lp: Strimko puzzle

futoshiki.lp: Futoshiki puzzle

fill_a_pix.lp: Fill-a-pix puzzle
This is a Minesweeper variant where the object is to fill a picture. It's almost as the Minesweeper encoding (see above) except that "this" cell is also counted in the hints (which I tend to forget) so the sum around and including a cell can be 9 instead of 8 (as in Minesweeper). Also, the third variant of adj/4 must be added to regard that each cell is its own adjacent.
nums(0..9).
size(1..n).
adj(R,C,R1,C1) :- size(R;C;R1;C1), |R-R1| + |C-C1|==1.
adj(R,C,R1,C1) :- size(R;C;R1;C1), size(C1), |R-R1|==1, |C-C1|==1.
adj(R,C,R,C) :- size(R;C).
N { x(R2,C2) : adj(R2,C2,R1,C1) } N :- hint(R1,C1,N).
seseman.lp: Seseman convent problem

coins_grid.lp: Coins grid problem
The object is to minimize the sum of the quadratic distances from the main diagonal. As noted before this is a problem that is very easy for MIP solvers, but harder for traditional CP solvers.
#const n = 10. % 31.
#const c = 4. % 14.
values(0..1).
size(1..n).
1 { x(Row, Col, Val) : values(Val) } 1 :- size(Row;Col).
c [ x(Row, Col, Val) : size(Row) : values(Val) = Val ] c :- size(Col).
c [ x(Row, Col, Val) : size(Col) : values(Val) = Val ] c :- size(Row).
#minimize [ x(Row, Col, Val) : values(Val) = Val*|Row-Col|*|Row-Col| ].
It's interesting to see that for a small problem (n=10, c=4) it is much faster than for example the MiniZinc model coins_grid.mzn and Gecode/fz as solver (I have tried to get a better heuristics but not succeeded). It takes 0.93 seconds to solve this with gringo/clasp, and more than 5 minutes with Gecode/fz. Even with the solution is assigned (z = 98) Gecode/fz takes 7 second.

A MIP solver solves the full problem (n=31, c=14) in 0.1 second, but it takes much longer for the ASP solver (I haven't waited for it to end).

Set covering/set parition

For some reason I like set covering and set partition problems. Most are from OR books or other sources (Winston, Taha, etc). They share a core of constraints but most have some specific twist.

set_covering3.lp: Assigning senators to committee (Katta G Murty)
Here is a version where the senator is represented as integers from 1..10. As we have seen before, the problem instance is represented not as a matrix but as the predicate senator/2.
senator(1, southern;conservative;republicans).
senator(2, southern;liberals;republicans).
senator(3, southern;liberals;democrats).
senator(4, southern;democrats).
senator(5, southern;conservative;democrats).
senator(6, northern;conservative;democrats).
senator(7, northern;conservative;democrats).
senator(8, northern;liberals;republicans).
senator(9, northern;liberals;democrats).
senator(10, northern;liberals;republicans).
senators(1..num_senators).
groups(southern;northern;liberals;conservative;democrats;republicans).
1 { selected(Senator) : senators(Senator) : senator(Senator, Group) } :- groups(Group).
#minimize [selected(Senator) : senators(Senator) ].
The ";"'s in the predicates are enumerating shortcuts, so first line for senator 1 is really a shortcut for these three lines:
senator(1, southern).
senator(1, conservative).
senator(1, republicans).
set_covering3_b.lp: Assigning senators to committee (Katta G Murty)
This is the "symbolic" version of the above.

set_covering.lp: Placing firestations (Winston: "Operations Research")

set_covering_b.lp: Placing firestations, alternative approach

set_covering2.lp: Number of security telephones on campus (Taha: "Operations Research")

set_covering4.lp: Set convering (and set partition) problem (Lundgren, Rönnqvist, Värbrand "Optimeringslära")

set_covering_opl.lp: Selecting workers to build a house (from an OPL model)

set_covering_skiena.lp: Set covering problem from Skiena
Here I use a constant set_covering which is default 1 which mean to run it as a set covering problem. Setting this to 0, (clingo -c set_covering=0 set_covering_skiena.lp) it will be a set partition problem.

combinatorial_auction.lp: Combinatorial auction


bus_scheduling_csplib.lp: Bus driver scheduling from CSPLib. The problem instances (see below) are converted from CSPLib #22. I have not really tried more than the simplest problem instance (t1) with my ASP encoding.

assignment.lp: Assignment Problems from Winston "Operations Research"

ski_assignment.lp: Ski assignment problem

Asorted puzzles

And last some problems from the worlds of recreational mathematics/logic .

safe_cracking.lp
This is an alphametic problem, but in contrast to send_more_money etc it is solved fast. The first version also used x(5,C5) in the safe_cracking predicate and it took almost 2 seconds to solve. After commenting out this variables, it took 0.2 second.
values(1..9).
% all different
1 { x(I, Val) : values(Val) } 1 :- values(I).
1 { x(I, Val) : values(I) } 1 :- values(Val).
safe_cracking :-
    x(1,C1), x(2,C2),x(3,C3),x(4,C4), % x(5,C5),
    x(6,C6),x(7,C7),x(8,C8), x(9,C9),  
    C4 - C6 == C7,
    C1 * C2 * C3 == C8 + C9,
    C2 + C3 + C6 < C8,
    C9 < C8.
:- not safe_cracking.        
% no fix points
:- x(I,I), values(I).
marathon2.lp: Marathon puzzle (XPress)

mr_smith.lp: Smith family problem (IF Prolog), a logic problem
Here we can see how the ASP constructs can be used to express different boolean expressions.
persons(mr_smith;mrs_smith;matt;john;tim).
value(go;stay).
% unique indices of person
1 { action(P, T) : value(T) } 1 :- persons(P).

% If Mr Smith comes, his wife will come too.
action(mrs_smith, go) :- action(mr_smith, go).

% At least one of their two sons Matt and John will come.
1 { action(matt, go), action(john, go) }.

% Either Mrs Smith or Tim will come but not both.
1 { action(mrs_smith, go), action(tim, go) } 1.

% Either Tim and John will come or neither will come.
:- action(tim, T1), action(john, T2), T1 != T2.

% If Matt comes then John and his father will also come.
2 { action(john, go), action(mr_smith, go) } 2 :- action(matt, go).
just_forgotten.lp: Just forgotten puzzle (Enigma 1517)

place_number.lp: Place number problem

a_round_of_golf.lp: A round of gold (Dell Logic Puzzles)
The hardest part to get right in this encoding was the requirement that in MiniZinc is quite simple to state:
(
 (score[Frank] = score[Sands] + 4 /\ score[caddy] = score[Sands] + 7 )
 \/
 (score[Frank] = score[Sands] + 7 /\ score[caddy] = score[Sands] + 4 )
)
I have translated it as the following. Note the integrity constraint to ensure that we don't want both of them to be false.
p3d1 :-
     score(frank, FrankScore),
     last_name(Sands, sands),
     first_names(Sands),first_names(Caddy),
     score(Sands, SandsScore),
     job(Caddy, caddy),
     score(Caddy, CaddyScore),
     FrankScore == SandsScore + 4,
     CaddyScore == SandsScore + 7.

p3d2 :-
     score(frank, FrankScore),
     last_name(Sands, sands),
     first_names(Sands),first_names(Caddy),
     score(Sands, SandsScore),
     job(Caddy, caddy),
     score(Caddy, CaddyScore),
     FrankScore = SandsScore + 7,
     CaddyScore = SandsScore + 4.

:- not p3d1, not p3d2. % not both false

Transitive closure

One thing that should also be mentioned is transitive closure which is some kind of trademark for Answer Set Programming (as well as planning). I have not used this much I have implemented mostly my traditional CP "learning problems" that does not use this type of structure. Here are two examples of transitive closures, and it is the succintness of these encodings that I made me so fascinated by ASP.

Hamiltonian cycle

For example Hamiltonian cycle can be encoded like this in ASP. The predicate e(X,Y) is the edge between two nodes, and v(X) is the vertices. The cycle is hardcoded to start from node 0.
{in(X,Y)} :- e(X,Y).

:- 2 {in(X,Y) : e(X,Y)}, v(X).
:- 2 {in(X,Y) : e(X,Y)}, v(Y).

r(X) :- in(0,X), v(X). % start from node 0
r(Y) :- r(X), in(X,Y), e(X,Y).

:- not r(X), v(X).

Clique

Here is an encoding of the problem finding maximum cliques, from Wikipedia Answer_set_programming#Large_clique:
v(X) :- e(X,Y).

n {in(X) : v(X)}.
% n {in(X) : v(X)} n. % exact size of clique

% Clique criteria
:- in(X), in(Y), v(X), v(Y), X!=Y, not e(X,Y), not e(Y,X).

% Find maximum clique
#maximize [in(X) : v(X)].
I have implemented an example of the maximum clique problem in the ASP program clique.lp with the problem instances clique_data1.lp, clique_data2.lp, clique_data2.lp.

Future projects

As mentioned above, Answer Set Programming is a large subject and this blog post has only showed a small part. Here are some of my possible future projects regarding Answer Set Programming:
  • learn how to tweak clasp etc.
  • look more into clingcon, the hybrid of ASP and CP
  • look more into iclingo (another Potassco tool), an incremental grounder/solver which can be used to incrementally step different variables for example in planning problems
  • planning: I first read about ASP when I was interested in automated planning so it might be time to continue this track
  • transitive closures: There are a lot of fun problems that use transitive closures. One example that I have in mind is the Rogo puzzle that Mike Trick wrote about some weeks ago in Operations Research, Sudoko, Rogo, and Puzzles. ASP might be good and natural approach to this. (I did a first try with ASP on Rogo but didn't get it right...)
  • DLV and other ASP grounders/solvers. There is a lot of more or less experimental tools that use different techniques and encoding schemes.

My Answer Set Programming programs

Here are my ASP programs collected for easy reference together with their problem instances. Almost all of of them have been implemented in some Constraint Programming system before, for example MiniZinc. See Common constraint programming problems for a list of different CP implementations.

All are written to comply with the Potassco tools (gringo/clasp/clingo), but should be quite easy to use with to other ASP systems, e.g. Lparse/Smodels.