polyomino_tutorial
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.
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 transposition is a tiling where isometries are restricted to translations.
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 straight. 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).
It 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)
but, 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 make up for these issues by:
- enumerating the valid configurations of a single place;
- giving up the binary encoding (which doesn't offer any advantage) and using a vector of shapes.
A piece, a single-piece-on-the-board or a full configuration of the board can be represented using integer matrices:
using shape = vita::matrix<int>;For instance:
// `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 for distinct pieces allows piece recognition when pieces are placed side by side on the board).
The configuration in the first picture could be expressed as:
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 the enumeration of 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 ? vita::fliplr(piece) : piece);
shape piece1(vita::rot90(flipped, rotation));
shape piece_on_board(put(piece1, empty, 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, the placement is not considered (putreturnsempty).
piece_masks contains the objects available in our combinatorial problem (we don't directly store them in the chromosome, preferring a simple integer index):
piece_masks[i][j] = /* matrix representing the j-th configuration of the i-th piece */ This setup reduces the search space from 2^117 to about 6.19E28 ≈ 2^96 elements.
Now the GAs related coding:
problem prob;
prob.env.individuals = 500;
prob.env.generations = 1000;
for (const auto &piece : piece_masks)
prob.sset.insert<ga::integer>( range(0, piece.size()) );The for loop defines the chromosome format:
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:
auto f = [](const i_ga &ind) -> fitness_t
{
shape board(board_height, board_width);
for (std::size_t i(0); i < ind.size(); ++i)
board += piece_masks[i][ind[i].as<int>()];
double filled(std::count_if(board.begin(), board.end(),
[](unsigned v) { return v != 0; }));
return {filled};
};which simply counts the non-empty locations.
Finally, to start the search:
ga_search<i_ga, std_es, decltype(f)> search(prob, f);
auto result = search.run(10);How does it work?
...
Run 9. 0 ( 0%): fitness (39)
Run 9. 0 ( 0%): fitness (40)
Run 9. 0 ( 1%): fitness (41)
Run 9. 0 ( 3%): fitness (43)
Run 9. 0 ( 5%): fitness (47)
Run 9. 0 ( 46%): fitness (48)
Run 9. 0 ( 85%): fitness (50)
Run 9. 1 ( 18%): fitness (51)
Run 9. 3 ( 34%): fitness (52)
Run 9. 4 ( 98%): fitness (53)
Run 9. 5 ( 88%): fitness (54)
Run 9. 7 ( 51%): fitness (55)
Run 9. 9 ( 81%): fitness (56)
Run 9. 13 ( 36%): fitness (57)
Run 9. 20 ( 99%): fitness (58)
Run 9. 38 ( 11%): fitness (59)
Run 9. 42 ( 4%): fitness (61)
Run 9. 58 ( 23%): fitness (62)
[INFO] Elapsed time: 2.674s
[INFO] Training fitness: (62)
Best result:
154 149 47 22 147 65 89 103 3 141 78 71 15
fitness (62)
J C C A B B L L
J C C A A B L L
J J J A B B + L
G G F F H E E E
D G + F H H K E
D D D F F H K .
D M D . I H K K
M M M M I I I I
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)
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 += vita::random::element(piece);
int filled(std::count_if(board.begin(), board.end(),
[](unsigned 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.
We have introduced just a small improvement. 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.
The circled function checks if a location is surrounded (locations at North, South, East and West are occupied or out of the frame) while the circled_zero function counts how many surrounded empty locations are present on the board.
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 (Xeon E3-1230v3), the best result obtained is a board with 57 locations filled (a blind brute force search is even worse). The basic GAs-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)
Of course the "circled location" trick can be embedded in the GAs code, clearly it's a just a small improvement, yet it's an interesting opportunity.
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 searching heuristics:
auto n_circled(static_cast<double>(circled_zero(board)));
double filled(std::count_if(board.begin(), board.end(),
[](unsigned 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 GAs (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) toa 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 i_ga &ind) -> fitness_t
{
shape board(board_height, board_width);
for (std::size_t i(0); i < ind.size(); ++i)
{
const auto mask(piece_masks[i][ind[i].as<int>()]);
if (!crash(board, mask))
board += mask;
}
auto n_circled(static_cast<double>(circled_zero(board)));
double filled(std::count_if(board.begin(), board.end(),
[](unsigned v) { return v != 0; }));
return {filled - n_circled};
};The crash function signals an illegal placement (overlapping piece) and 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 are frequent (and they're without hope).
ALPS helps in this kind of problems (see [2]).
(the source code of the improved example is contained in the examples/polyomino02.cc file)
