Polyomino puzzle
The general challenge posed is to tile a given region with a given set of polyominoes.
The puzzle in the picture comes from a Stackoverflow's question and is used extensively in this tutorial.
Generalities
A polyomino is a simply connected tile obtained by gluing together rookwise connected unit squares.
A tiling of a region by a set of polyominoes is a partition of the region into images of the tiles by isometries.
A tiling by translation is a tiling where isometries are restricted to translations.
First idea
Our reference puzzle has 13 polyominoes (which can be rotated and reflected) and a 8x8 board.
In this case, the classical binary encoding is quite straightforward. A piece can be:
- the original way or flipped left-right (1 bit);
- rotated counterclockwise by 0 (i.e. none), 90, 180 or 270 degrees (2 bits);
- at position
(x, y), wherexandygo from0to7(3 bits for each coordinate).
This sums up to 9 bits for a piece. Because we have 13 pieces, a total of 117 bits is required.
*** ***
* * *
1st piece 2nd piece
F R Y X F R Y X
| 0 00 010 010 | 1 01 101 100 | ...
^ ^ ^ ^ ^ ^ ^ ^
| | | | | | | +---- x position 4 **
| | | | | | +-------- y position 5 *
| | | | | +------------ rotated by 90 deg **
| | | | +-------------- flipped
| | | |
| | | +------------------- x position 2
| | +----------------------- y position 2 ***
| +--------------------------- not rotated *
+----------------------------- not flipped
The fitness can be calculated by placing each piece in the frame, ignoring any parts that lie out of the frame, and then adding up the number of empty squares. When that hits zero, we have a solution.
The model is simple but there are some questionable aspects:
-
The user should select the smallest alphabet that permits a natural expression of the problem.
(Genetic Algorithms - David E. Goldberg - chapter about Codings)
Unfortunately, the same configuration can be expressed in multiple ways. E.g.
PIECE EQUIVALENT CODINGS *** 0 00 ... ... * 1 00 ... ... --- ** 0 00 ... ... * 1 10 ... ... ** -
Considering the full range of coordinates (
[0;7] x [0;7]) for every piece is too much (and somewhat misleading for fitness evaluation).
We can address these issues by:
- enumerating the valid configurations of a single piece;
- giving up the binary encoding (which doesn't offer any advantage) and using a vector of shapes.
Start coding
A piece, a piece-on-board, or a full configuration of the board can be represented using integer matrices:
using shape = ultra::matrix<int>;
For instance, a piece can be represented as:
// `0` marks an empty square
shape t_tetromino = { {1, 1, 1},
{0, 1, 0} };
or, even better:
shape t_tetromino = { {'T', 'T', 'T'},
{ 0 , 'T', 0 } };
(using different characters/values allows piece recognition when pieces are placed side by side on the board).
The same piece at (2, 2) can be represented as a piece-on-board:
{ {0, 0, 0 , 0 , 0 , 0 0, 0},
{0, 0, 0 , 0 , 0 , 0 0, 0},
{0, 0, 'T', 'T', 'T', 0, 0, 0},
{0, 0, 0 , 'T', 0 , 0, 0, 0},
{0, 0, 0 , 0 , 0 , 0 0, 0},
{0, 0, 0 , 0 , 0 , 0 0, 0},
{0, 0, 0 , 0 , 0 , 0 0, 0},
{0, 0, 0 , 0 , 0 , 0 0, 0} }
The configuration in the first picture could be expressed as a full configuration:
shape solution = { {'T', 'T', 'T', 'C', 'C', 'C', 'Q', 'Q'},
{'H', 'T', 'H', 'C', '*', 'C', 'Q', 'Q'},
{'H', 'H', 'H', '*', '*', '*', 'S', 'S'},
{'H', '~', 'H', '*', '=', '=', 'S', '!'},
{'~', '~', '=', '=', '=', 'S', 'S', '!'},
{'~', '>', '>', '>', 'L', 'L', 'L', '!'},
{'R', 'R', 'R', '>', 'L', 't', 'L', 'L'},
{'R', 'R', 'R', '>', 't', 't', 't', 't'} };
The code for enumerating valid configurations is:
const std::size_t board_height = 8;
const std::size_t board_width = 8;
std::vector<std::vector<shape>> piece_masks;
std::size_t add_piece_variants(const shape &piece)
{
std::set<shape> ms;
const shape empty(board_height, board_width); // filled with `0`s
for (unsigned reflection(0); reflection <= 1; ++reflection)
for (unsigned rotation(0); rotation <= 3; ++rotation)
for (unsigned y(0); y < board_height; ++y)
for (unsigned x(0); x < board_width; ++x)
{
shape flipped(reflection ? ultra::fliplr(piece) : piece);
shape flip_rot(ultra::rot90(flipped, rotation));
shape piece_on_board(put(flip_rot, y, x));
if (piece_on_board != empty)
ms.insert(piece_on_board);
}
piece_masks.emplace_back(ms.begin(), ms.end());
return ms.size();
}
There are two interesting points:
- configurations are initially inserted in a set (
std::set<shape> ms) and then in the final data structure (piece_masks). The set filters out duplicates; - the
putfunction tries to place a piece on the board at position(y, x). If parts of the piece lie out of the frame (putreturnsempty), the placement is not considered.
piece_masks contains the objects available in our combinatorial problem (we don't directly store them in the genome, preferring a simple integer index):
piece_masks[i][j] = /* j-th configuration of the i-th piece */
This setup reduces the search space from \(2^{117}\) to about \(6.19E28 ≈ 2^{96}\) elements.
Genetic Algoritms (GA) related coding
The for loop defines the chromosome format:
ga::problem prob;
prob.params.population.individuals = 500;
prob.params.evolution.generations = 1000;
// The chromosome is a sequence of bounded integers (indices) used to access
// the `piece_masks` data structure.
for (const auto &piece : piece_masks)
prob.insert( interval(0, piece.size()) );
Graphically:
Piece 1 Piece 2 Piece 13
↕ ↕ ↕
[ integer in range R0 | integer in range R1 | ... | integer in range R12 ]
The other essential element of a GA is the fitness function, which simply counts the non-empty locations:
auto f = [](const ga::individual &ind)
{
shape board(board_height, board_width);
for (std::size_t i(0); i < ind.parameters(); ++i)
board += piece_masks[i][ind[i]];
// Number of non-empty squares.
double filled(std::ranges::count_if(board,
[](auto v) { return v != 0; }));
return filled;
};
Finally, to start the search, use the following code:
ga::search search(prob, f);
auto result = search.run(10);
How does it work?
...
0.000 0: 50
0.001 1: 52
0.002 2: 53
0.003 3: 54
0.005 5: 55
0.006 7: 56
0.009 10: 57
0.015 15: 58
0.019 18: 59
0.025 21: 60
0.086 43: 61
0.139 56: 62
[INFO] Evolution completed at generation: 1002. Elapsed time: 6.494
Run 9 TRAINING. Fitness: 62
Best result:
150 48 21 61 43 0 84 8 181 131 26 83 123
fitness 62
L L L J A D D D
L L L J A A D .
I + J J A D D D
I G M M M M C C
I G G B M B C C
I E G B B B K F
E E E H H F + F
E H H H . F K K
Seems good, but
- how to be sure?
- can it be improved?
(the full source code of the example is contained in the examples/polyomino01.cc file)
Random search
A random search is a good reference test and it isn't difficult to code:
void random_put(const shape &base)
{
int max(0);
while (max < 64)
{
shape board(base);
for (const auto &piece : piece_masks)
board += ultra::random::element(piece);
// Number of non-empty squares.
int filled(std::ranges::count_if(board,
[](auto v) { return v != 0; }));
if (filled > max)
{
max = filled;
std::cout << filled << '\n';
print_board(board);
}
}
}
This is a completely uninformed random search and the infrastructure (shape, piece_masks, put, add_piece_variant) is shared with the previous example.
Only a small improvement has been introduced. Consider the following situation:
It's a legal but illogical placement. Given our set of pieces, every surrounded single square cannot be filled. Every time we place a piece in such a way, we're giving away our winning chances.
Excluding this kind of placement, we obtain a smaller piece_masks data structure and a slightly smaller search space of \(2.86E28 ≈ 2^{95}\) elements.
Even so, after many hours, the best result obtained is a board with 57 locations filled (a blind brute force search is even worse). The basic GA-based example reaches a better configuration (62 filled locations) in a few seconds.
The source code of this test is contained in the examples/polyomino00.cc file.
There are two new functions: circled checks if a location/square is surrounded (locations at North, South, East and West are occupied or out of the frame), while the circled_zero counts how many surrounded empty locations are present on the board.
Improving the search
Of course the "circled location" trick can be embedded in the GA code. Clearly it's a just a small improvement.
We could look for circled locations not only during the initial enumeration, but even during the evolution. The "circled location" can be used as a effective heuristics:
auto n_circled(static_cast<double>(circled_zero(board)));
// Number of non-empty squares.
double filled(std::ranges::count_if(board,
[](auto v) { return v != 0; }));
return {-n_circled, filled};
changing the fitness from scalar to multi value (and falling back upon lexicographical comparison), the search is restricted / forced toward circled-free-boards (-n_circled tends to 0).
Other equally valid options are:
-
Invert the terms
return {filled, -n_circled};(strangely works well);
-
Combine the values
return filled - n_circled;
Indeed we're improving but a lot of (long) runs are required to find a solution. This kind of combinatorial problem can be very hard to solve with GA (because of the many local optima).
All right, we still have a few tricks up our sleeve.
Consider the following situation:
Current board Available pieces
----------
|SSLLL | IIII RRRR
|SSL | RRRR
----------
The fitness function prefers an illegal variant (placing the rectangle with a partial overlap) to a legal one:
---------- ----------
|SSLL+RRR| Fitness = 15 |SSLLL | Fitness = 12
|SSL RRRR| |SSLIIII |
---------- ----------
The fitness function is misleading: an illegal position shouldn't be scored more than a legal one. It's not a difficult change to code:
auto f = [](const ga::individual &ind)
{
shape board(board_height, board_width);
for (std::size_t i(0); i < ind.size(); ++i)
{
const auto mask(piece_masks[i][ind[i]]);
if (!crash(board, mask))
board += mask;
}
auto n_circled(static_cast<double>(circled_zero(board)));
double filled(std::ranges::count_if(board,
[](auto v) { return v != 0; }));
return filled - n_circled;
};
The crash function signals an illegal placement (overlapping piece) and the score due to the conflicting polyomino is entirely skipped.
With this change the fitness function becomes more informative and almost correct.
Last but not least we observe that identifying when a restart is appropriate is a difficult task. It happens that the evolution get caught in a local optimum (near the absolute optimum) and a sequence of mutations could be enough to move on. Unfortunately also local optima far from the absolute one(s) are frequent (and they're without hope).
In this situation ALPS meta-heuristics helps to avoid premature convergence (see 9).
(the source code of the improved example is contained in the examples/polyomino02.cc file)
Conclusion
GA aren't the most appropriate tool for this kind of puzzle. GA tend to produce good but sub-optimal solutions and this behaviour is acceptable only for some combinatorial optimization problems.
There are better approaches and further improvements (e.g. see the original question on Stackoverflow and @TodorBalabanov's answer).
Anyway the last example proposed finds a solution almost always (often in a short time) and GA prove themselves a general, viable path especially when previous knowledge isn't available.