This is the one hundred and first part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra
Today we’re going to solve Polyominoes. We have a set of polyominoes that needs to enclose an area of the greatest possible area. Let’s see the solution:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 |
var pieces = new []{ new []{ "XX", "X ", "X " }, new[]{ " XX", "XX " }, new[]{ "XX ", " XX" }, new[]{ "X ", "XX", "X " }, new[]{ "XX", "XX" }, new[]{ "XXXX" }, new[]{ "X ", "X ", "XX" } }; var totalLength = pieces.Select(p => (int)Math.Max(p.Length, p[0].Length)).Sum(); var maxLength = pieces.Max(p => (int)Math.Max(p.Length, p[0].Length)); var width = totalLength + 2; var height = width; var piecesCount = pieces.Length; // Right, down, left, up var directions = 4; var heads = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => Enumerable.Range(0, piecesCount).Select(p => Enumerable.Range(0, directions).Select(d => solver.CreateAnonymous(Domain.BinaryInteger) ).ToArray() ).ToArray() ).ToArray() ).ToArray(); var isInterior = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => solver.CreateAnonymous(Domain.BinaryInteger) ).ToArray() ).ToArray(); var isBoundary = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => solver.CreateAnonymous(Domain.BinaryInteger) ).ToArray() ).ToArray(); Func<string[], int, string[]> rotatePiece = (piece, direction) => { if(direction == 0){ return piece; }else if(direction == 2){ return piece.Select(l => string.Join("", l.Reverse())).Reverse().ToArray(); }else if(direction == 1){ return Enumerable.Range(0, piece[0].Length) .Select(c => string.Join("", Enumerable.Range(0, piece.Length).Select(w => piece[piece.Length - 1 - w][c]))) .ToArray(); }else { return Enumerable.Range(0, piece[0].Length) .Select(c => string.Join("", Enumerable.Range(0, piece.Length).Select(w => piece[w][piece[0].Length - 1 - c]))) .ToArray(); } }; // Put pieces on the board for(int p=0;p<piecesCount;++p){ solver.Set<Equal>( solver.FromConstant(1), solver.Operation<Addition>( Enumerable.Range(0, height).SelectMany(h => Enumerable.Range(0, width).SelectMany(w => Enumerable.Range(0, directions).Select(d => heads[h][w][p][d] ) ) ).ToArray() ) ); } // Make sure pieces don't fall off the board for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int p=0;p<piecesCount;++p){ for(int d = 0; d < directions;++d){ var piece = rotatePiece(pieces[p], d); if(d == 0){ if(w + piece[0].Length - 1 >= width){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } }else if(d == 1){ if(h + piece.Length - 1 >= height){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } }else if(d == 2){ if(w - piece[0].Length + 1 < 0){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } }else if(d == 3){ if(h - piece.Length + 1 < 0){ heads[h][w][p][d].Set<Equal>(solver.FromConstant(0)); } } } } } } // Set impacts IList<Tuple<int, int, int, int>>[][] allImpacts = Enumerable.Range(0, height).Select(h => Enumerable.Range(0, width).Select(w => new List<Tuple<int, int, int, int>>() ).ToArray() ).ToArray(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int p=0;p<piecesCount;++p){ for(int d = 0; d < directions;++d){ var piece = rotatePiece(pieces[p], d); var pieceWidth = piece[0].Length; var pieceHeight = piece.Length; for(int y=0;y<pieceHeight;++y){ for(int x=0;x<pieceWidth;++x){ var isPieceX = piece[y][x] == 'X'; if(!isPieceX) continue; var finalY = d == 0 || d == 1 ? h + y : y - pieceHeight + 1 + y; var finalX = d == 0 || d == 1 ? w + x : x - pieceWidth + 1 + x; if(finalY >= 0 && finalY < height && finalX >=0 && finalX < width){ allImpacts[finalY][finalX].Add(Tuple.Create(h, w, p, d)); } } } } } } } // Set the boundary for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ isBoundary[h][w].Set<Equal>( solver.Operation<Addition>( allImpacts[h][w].Select(t => heads[t.Item1][t.Item2][t.Item3][t.Item4]).ToArray() ) ); } } // Make sure pieces do not overlap for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ solver.Set<LessOrEqual>( solver.Operation<Addition>( allImpacts[h][w].Select(t => heads[t.Item1][t.Item2][t.Item3][t.Item4]).ToArray() ), solver.FromConstant(1) ); } } // Extend interior for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int y=-1;y<=1;++y){ for(int x=-1;x<=1;++x){ if(y == 0 && x == 0){ continue; } if(h + y < 0 || h + y >= height){ continue; } if(w + x < 0 || w + x >= width){ continue; } solver.Set<Equal>( solver.FromConstant(1), solver.Operation<MaterialImplication>( solver.Operation<Conjunction>( isInterior[h][w], isBoundary[h + y][w + x].Operation<BinaryNegation>() ), isInterior[h + y][w + x] ) ); } } } } // Make sure there is some exterior for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ if(h == 0 || h == height - 1 || w == 0 || w == width - 1){ solver.Set<Equal>( solver.FromConstant(0), isInterior[h][w] ); solver.Set<Equal>( solver.FromConstant(0), isBoundary[h][w] ); } } } // Make sure a field cannot be both exterior and boundary for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ solver.Set<LessOrEqual>( solver.Operation<Addition>( isInterior[h][w], isBoundary[h][w] ), solver.FromConstant(1) ); } } var goal = solver.Operation<Addition>(isInterior.SelectMany(x => x).ToArray()); solver.AddGoal("goal", goal); solver.Solve(); Console.WriteLine(solver.GetValue(goal)); var board = Enumerable.Range(0, height).Select(h => new String(' ', width).ToCharArray()).ToArray(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ for(int p=0;p<piecesCount;++p){ for(int d = 0; d < directions;++d){ var piece = rotatePiece(pieces[p], d); var pieceWidth = piece[0].Length; var pieceHeight = piece.Length; if(solver.GetValue(heads[h][w][p][d]) > 0.5){ for(int y=0;y<pieceHeight;++y){ for(int x=0;x<pieceWidth;++x){ var finalY = d == 0 || d == 1 ? h + y : y - pieceHeight + 1 + y; var finalX = d == 0 || d == 1 ? w + x : x - pieceWidth + 1 + x; if(finalY >= 0 && finalY < height && finalX >=0 && finalX < width){ board[finalY][finalX] = piece[y][x]; }else { Console.WriteLine("Error in " + h + " " + w + " " + p + " " + d + " " + y + " " + x); } } } } } } } } foreach(var l in board){ Console.WriteLine(string.Join("", l)); } Console.WriteLine(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ Console.Write(Math.Round(solver.GetValue(isInterior[h][w]))); } Console.WriteLine(); } Console.WriteLine(); for(int h = 0; h < height; ++h){ for(int w = 0; w < width; ++ w){ Console.Write(Math.Round(solver.GetValue(isBoundary[h][w]))); } Console.WriteLine(); } Console.WriteLine(); |
We start with the definition of pieces (lines 1-32). They are encoded as strings so we can rotate them easily later on. Next goes a couple of helper variables.
Our model has three sets of variables. First is heads
that indicates if a given piece was put on a board at a given place with a given rotation (lines 42-50).
Next, we have bit flags indicating if a given field belongs to the interior of the enclosed area (lines 52-56), and similar flags for the boundary (58-62).
Next, we have a helper function to rate a piece (lines 64-78).
And then we go with rules. First, we make sure that each piece is put somewhere. So we sum all of the flags for each piece, and set it to one (lines 80-94).
Then, we need to make sure that we don’t put pieces too close to the boundary. We set corresponding flags to zer (lines 96-122).
Next, we need to exclude some invalid solutions. First, for every field we calculate which pieces could “impact it”. So we check all the positions of all the pieces in the all possible rotations, and store the result in a helper array. This is in lines 124-154.
Now, we can make sure that there is a boundary. It needs to be only in the fields covered by the polyominoes, and nowhere else. So we take all the pieces that impact a given field, and make sure that the boundary is set to the sum of those. This is in lines 156-165.
Next, we need to make sure that pieces do not overlap. So for each field we make sure that there is at most one impacting piece selected (lines 167-177).
Next, we extend the interior. We encode the logic that if a given field is internal and there is a field right next to it that is not on the boundary, then the neighbouring field must be internal as well. This is in lines 180-210. We also need to make sure some fields are external (lines 212-226). We do it for a whole “frame” encompassing the board, otherwise the solver might decide to put everything “inside-out”, so the interior is actually the exterior that is implicitly bounded by the area outside of the board.
This works in the following way. When a field is truly in the interior, then it will expand as much as possible. However, fields outside of the shape cannot be marked as internal because then there would be a case of neighbours of the frame. We take any field which is right next to the frame, and since the frame is not on the boundary, it needs to be marked as internal as well. However, we explicitly block that.
Finally, we need to make sure there are no fields that are both interior and boundary (which otherwise would be allowed and might change the optimal solution). This is in lines 228-239).
Output:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 |
Tried aggregator 2 times. MIP Presolve eliminated 4659 rows and 4747 columns. MIP Presolve modified 688 coefficients. Aggregator did 1432 substitutions. Reduced MIP has 809 rows, 847 columns, and 4647 nonzeros. Reduced MIP has 847 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.03 sec. (10.92 ticks) Found incumbent of value 0.000000 after 0.06 sec. (33.44 ticks) Probing time = 0.00 sec. (2.70 ticks) Tried aggregator 1 time. Reduced MIP has 809 rows, 847 columns, and 4647 nonzeros. Reduced MIP has 847 binaries, 0 generals, 0 SOSs, and 0 indicators. Presolve time = 0.02 sec. (2.75 ticks) Probing time = 0.02 sec. (2.69 ticks) Clique table members: 1918. MIP emphasis: balance optimality and feasibility. MIP search method: dynamic search. Parallel mode: deterministic, using up to 8 threads. Root relaxation solution time = 0.02 sec. (17.35 ticks) Nodes Cuts/ Node Left Objective IInf Best Integer Best Bound ItCnt Gap * 0+ 0 0.0000 35.0000 565 --- 0 0 29.7678 243 0.0000 29.7678 565 --- * 0+ 0 2.0000 29.7678 657 --- 0 0 29.2337 242 2.0000 Cuts: 36 657 --- * 0+ 0 3.0000 29.2337 657 874.46% 0 0 28.7133 242 3.0000 Cuts: 99 816 857.11% 0 0 28.5944 248 3.0000 Cuts: 26 881 853.15% 0 0 28.3575 243 3.0000 Cuts: 21 958 845.25% 0 0 28.1942 256 3.0000 Cuts: 27 1036 839.81% 0 0 28.1076 244 3.0000 Cuts: 25 1120 836.92% * 0+ 0 4.0000 28.1076 1120 602.69% 0 0 28.1076 245 4.0000 Cuts: 12 1130 602.69% 0 0 28.0851 255 4.0000 Cuts: 9 1150 597.76% 0 0 27.9794 239 4.0000 Cuts: 28 1229 597.76% 0 0 27.8878 280 4.0000 Cuts: 25 1330 597.20% 0 0 27.8093 275 4.0000 Cuts: 23 1402 595.23% 0 0 27.6921 275 4.0000 Cuts: 40 1520 592.30% 0 0 27.6921 276 4.0000 Cuts: 7 1525 592.30% * 0+ 0 20.0000 27.6921 1525 38.46% 0 2 27.6921 276 20.0000 27.6060 1525 38.03% Elapsed time = 0.97 sec. (335.52 ticks, tree = 0.01 MB, solutions = 5) * 5+ 5 21.0000 27.5336 1944 31.11% * 56 40 integral 0 22.0000 27.4817 5646 24.92% * 63+ 38 25.0000 27.4817 5959 9.93% Clique cuts applied: 15 Cover cuts applied: 19 Implied bound cuts applied: 63 Mixed integer rounding cuts applied: 35 Zero-half cuts applied: 53 Gomory fractional cuts applied: 1 Root node processing (before b&c): Real time = 0.97 sec. (335.26 ticks) Parallel b&c, 8 threads: Real time = 0.52 sec. (250.13 ticks) Sync time (average) = 0.26 sec. Wait time (average) = 0.00 sec. ------------ Total (root+branch&cut) = 1.49 sec. (585.39 ticks) 25 XXXXXX XX X XX X X XX X X XX XX XXXXXXX 00000000000 00000000000 00001111000 00011111000 00111111000 00111111000 00011110000 00000000000 00000000000 00000000000 00011111100 00110000100 01100000100 01000000110 01000000100 01100001100 00111111100 00000000000 |