One of my class at EPITA is about Constraint Programming. This is a technique to solve problems with
- A huge number of possible solutions. (2^1000 is not uncommon)
- Discrete variables (they can take a bounded number of values).
I wanted to share this method with you because it is able to solve an impressive number of problems without writing any line of code. How? Because it is declarative. You write the problem (CSP: Constraint Satisfaction Problem) and the solver is in charge of finding the solution.
Here are some basic examples. I used the trial version of IBM ILOG for all those examples.
Sudoku
The Sudoku is a perfect target for Constraint Programming.
// Variables range Size = 0..8; dvar int Sudoku[Size][Size] in 1..9; // Constraints subject to { forall(i in Size) { forall (j, k in Size: j != k) { // Lines & Columns Sudoku[i][j] != Sudoku[i][k]; Sudoku[j][i] != Sudoku[k][i]; // Squares Sudoku[(i % 3) * 3 + (j % 3)] [(i div 3) * 3 + (j div 3)] != Sudoku[(i % 3) * 3 + (k % 3)] [(i div 3) * 3 + (k div 3)]; } } // Initial Input forall (i, j in Size) { PreSolve[i][j] != 0 => Sudoku[i][j] == PreSolve[i][j]; } } |
Sadly, Constraint Programming allows to solve instantly all the sudokus. I chose the World Hardest Sudoku as an example.
PreSolve = [ [0 0 5 3 0 0 0 0 0] [8 0 0 0 0 0 0 2 0] [0 7 0 0 1 0 5 0 0] [4 0 0 0 0 5 3 0 0] [0 1 0 0 7 0 0 0 6] [0 0 3 2 0 0 0 8 0] [0 6 0 5 0 0 0 0 9] [0 0 4 0 0 0 0 3 0] [0 0 0 0 0 9 7 0 0] ]; // Solution [1 4 5 3 2 7 6 9 8] [8 3 9 6 5 4 1 2 7] [6 7 2 9 1 8 5 4 3] [4 9 6 1 8 5 3 7 2] [2 1 8 4 7 3 9 5 6] [7 5 3 2 9 6 4 8 1] [3 6 7 5 4 2 8 1 9] [9 8 4 7 6 1 2 3 5] [5 2 1 8 3 9 7 6 4] |
- Number of branches : 8
- Number of fails : 0
- Total memory usage : 1.1 Mb
- Time spent in solve : 0.02s
17x17 Challenge
I discovered the 17x17 Challenge. You have to assign one color per element of a matrix such as there is not any rectangle with all the corners of the same color.
// Variables dvar int Board[1..N][1..N] in 1..C; // Constraints subject to { forall (i, j in 1..(N - 1)) { forall (w in 1..(N - i), h in 1..(N - j)) { ! (Board[i][j] == Board[i + w][j] && Board[i][j + h] == Board[i + w][j + h] && Board[i][j + h] == Board[i + w][j]); } } } |
// Input int C = 14; int N = 4; // Solution 14x14 [1 1 2 4 3 3 2 1 4 4 1 3 4 4] [2 4 1 1 2 4 3 2 1 3 4 4 4 3] [2 3 4 3 1 2 2 3 2 3 4 2 1 4] [4 2 3 2 3 1 4 1 3 3 2 2 4 1] [4 2 3 3 1 3 2 4 1 2 1 4 2 4] [1 3 3 1 1 1 4 4 2 4 4 1 3 3] [3 2 2 4 2 4 4 1 2 3 3 1 1 2] [3 1 1 4 4 1 1 2 3 2 4 3 1 3] [3 3 2 1 3 4 3 4 4 2 1 2 1 1] [1 4 3 2 2 2 4 3 1 2 3 3 1 4] [3 1 4 1 4 3 4 3 4 2 2 4 3 2] [1 4 4 4 2 3 3 4 3 1 3 2 2 1] [4 4 2 3 4 2 3 2 1 4 2 3 2 1] [4 1 4 3 3 4 2 2 1 1 3 1 3 2] |
- Number of branches : 1,586,955
- Number of fails : 774,500
- Total memory usage : 13.6 Mb
- Time spent in solve : 89.40s
- Search speed (br. / s) : 17,750.3
For grids smaller than 14x14, the result is found instantly. A 14x14 takes 1min30. And for 15x15, I let it run during 24hours without any result 🙁 This naive approach will not give a solution any time soon for a 17x17 grid.
N-Queens
The 8-Queen problem and it's generalization to a NxN board is trivial to write in Constraint Programming.
There are three ways to code the problem. The unknown can be the board (boolean for the Queen presence) or the Queen position. A naive Queen position can be expressed with two coordinates X and Y. If we analyze the problem just a bit, we see that only one Queen can be per column. So we just store the column position of all the Queens.
The methods are sorted by speed of execution. The Board version is much slower than the Column one.
Board
dvar boolean Board[Size][Size]; subject to { forall (i in Size) { // One Queen per Line sum (j in Size) (Board[i][j]) == 1; // One Queen per Column sum (j in Size) (Board[j][i]) == 1; } forall (i, j, k, l in Size: i != k && j != l) { // One Queen per Diagonal abs(i - j) == abs(k - l) => Board[i][j] + Board[k][l] < = 1; } } |
Queen
dvar int QueenX[Size] in Size; dvar int QueenY[Size] in Size; subject to { forall (i, j in Size: i != j) { // One Queen per Line QueenX[i] != QueenX[j]; // One Queen per Column QueenY[i] != QueenY[j]; // Diagonals abs(QueenX[i] - QueenX[j]) != abs(QueenY[i] - QueenY[j]); } } |
Column
range Size = 1..N; dvar int Column[Size] in Size; subject to { forall (i, j in Size: i != j) { // One Queen per Line Column[i] != Column[j]; // Diagonal abs(Column[i] - Column[j]) != abs(i - j); } } |
Here is an example of result for a standard 8x8 board.
int N = 8; |---|---|---|---|---|---|---|---| | | | | | | X | | | |---|---|---|---|---|---|---|---| | | | | X | | | | | |---|---|---|---|---|---|---|---| | | X | | | | | | | |---|---|---|---|---|---|---|---| | | | | | | | | X | |---|---|---|---|---|---|---|---| | | | | | X | | | | |---|---|---|---|---|---|---|---| | | | | | | | X | | |---|---|---|---|---|---|---|---| | X | | | | | | | | |---|---|---|---|---|---|---|---| | | | X | | | | | | |---|---|---|---|---|---|---|---| |
Since the 8x8 example is instant to process, I made it run on a 100x100 board. You have to believe me to know that the result is correct 🙂
int N = 100; Column = [65 42 71 32 58 90 88 93 61 19 76 56 67 89 23 10 15 60 70 52 28 40 1 95 22 85 63 43 54 29 4 24 96 68 73 18 3 38 26 100 34 99 17 14 6 79 37 49 51 72 62 57 59 45 78 25 94 46 86 13 77 5 35 41 82 97 12 48 8 91 21 98 87 84 7 9 66 81 92 75 39 50 11 80 36 2 27 74 64 20 16 47 55 30 33 31 53 69 83 44]; |
- Number of branches : 30,958
- Number of fails : 12,535
- Total memory usage : 17.4 Mb
- Time spent in solve : 12.17s
Send More Money
Here we try to solve the Send More Money problem. We have to find what are the unique values for S, E, N, D, M ... such as this equation is true:
S E N D + M O R E --------- M O N E Y |
The implementation is fairly straightforward. We use 4 more variables that will be the carries.
{string} Letters = {"S", "E", "N", "D", "M", "O", "R", "Y"}; range Digit = 0..9; dvar int Values[Letters] in Digit; dvar int Carry[1..4] in 0..1; subject to { allDifferent (Values); Values["M"] != 0; Carry[4] == Values["M"]; Values["S"] + Values["M"] + Carry[3] == Values["O"] + 10 * Carry[4]; Values["E"] + Values["O"] + Carry[2] == Values["N"] + 10 * Carry[3]; Values["N"] + Values["R"] + Carry[1] == Values["E"] + 10 * Carry[2]; Values["D"] + Values["E"] == Values["Y"] + 10 * Carry[1]; } |
And we get instantly the following result:
9 5 6 7 + 1 0 8 5 --------- 1 0 6 5 2 |
Behind the scene, the algorithm used is Branch & Bound. It uses Look-ahead to reduce the domain definition of the variables and Backjumping in order to backtrack to the first instanced variable that caused a problem. There are heuristics to know the order of the variables to instanciate. A common approach is to try first the "hardest" variables in order to remove as many branches as possible.
If you want to know more, here are some links:
