Tiling a square

I wondered if there was a way to automate the process of finding a way to tile a tile into a square.

The idea is to represent the tile with a matrix of $0$s for blank space and $1$s for filled spaces like in $\{\{0,1,1\}, \{1,0,0\}\}$.

The square (or a general rectangular area) is represented by an array as well.

The function should then have as arguments two lists (representing the tile and the area) and should find the ways (if they exist) to do the tiling.

This is quite an hard problem for me, I have tried playing with Plus, so that if the result matrix has $2$ it means that the tile overlaps while if it is $0$ there is a gap but I can't find a way to tell Mathematica that they should fit the area.

-
Man thanks so much!I'll read it now but at a first glance seems very useful. Thanks again! – Gianpiero Cea Jun 14 '12 at 18:31
@belisarius I think this is nice and should be posted as an asnwer – Vitaliy Kaurov Jun 14 '12 at 19:48
@VitaliyKaurov Ok, done. I usually don't post as an answer something that came up from a Google search ... but that also depends on how long it took me to find the answer. – belisarius has settled Jun 14 '12 at 21:50

In this article the author solves the problem of tiling a rectangle by using pieces taken from a set of polyominoes, which are plane geometric figures formed by joining one or more equal squares edge to edge. For example, these are the pentaminoes, polyominoes formed by joining 5 squares:

Of course this problem is more difficult than the one you asked for, but it is also more interesting and ... there is Mathematica code in the article!

-
+1 Thank you for posting. I think that Mathematica Journal is great and people can find many gems there. I wish more people would read it. – Vitaliy Kaurov Jun 14 '12 at 22:03
Some time ago i had found that great resource,the Mathematica Journal,but my skills in Mathematica were really poor(and still they are :D)and so i eventually forgot it.Anyway i'm getting better and better and i think i'll start giving a deeper look to it.Thank you all for remembering me this great source! – Gianpiero Cea Jun 15 '12 at 8:20
+1 For letting me know about the existance of the Mathematica Journal. – Voyska Oct 20 '12 at 3:58
@belisarius Tried the Mathematica Journal code, and it won't run. Going through the issues with it show that the first line is unterminated, or...there is one undeclared variable "allPieces" that might have been meant to actually be "pieces" The function "tile" can not be included inside the module of the function "tess" because they both declare the same inputs. Love these Journal articles, but I do wish they would post working code. Does anyone have the actual code that runs? From all the errors in the article it's not clear exactly what is going to make it run. – R Hall Dec 14 '13 at 13:23
@RHall I don't remember testing the code, and probably I didn't because I would posted a warning if it wasn't working. If you really want others looking into it I think a "what is missing in this Mma journal code?" question could attract a lot of attention in this site – belisarius has settled Dec 14 '13 at 15:31

With some diffidence (because there appears to be a Mathematica bug: see below), I would like to offer an answer in the spirit of the OP's original attempt to solve the problem algebraically.

Solution

This problem can be formulated as a binary integer linear program. The reformulation represents the square (or more generally, a rectangle as implemented below) as a vector. Its coefficients represent occupancy: $0$ for empty, $1$ for full. In this representation, every possible tile placement corresponds to a unique vector indicating the cells occupied by the tile.

Let's begin by generating the vectors corresponding to all possible placements of tiles in all possible orientations. It will be convenient to specify as little as possible, beginning with a set of unique tile shapes.

As a running example I use the L, C, T, and Z tiles:

tileSet = {{{1, 1, 1}, {1, 0, 0}}, {{1, 1, 1}, {1, 0, 1}}, {{1, 1, 1}, {0, 1, 0}}, {{0, 1, 1}, {1, 1, 0}}};
ArrayPlot[#, Frame -> False, ImageSize -> 10 Max[Length[First[#]]]] & /@ tileSet


(There is no restriction that tiles have a common bounding box, but--to make sure they can be properly positioned--it is important that their matrix descriptions not have any zero padding along outside rows or columns.)

We need to rotate (and, if we wish, flip) all these tiles. An elegant way (perhaps not the simplest) is to generate the group action abstractly in terms of generators (alpha and beta).

apply[s_List, alpha] := Reverse /@ s;       (* Horizontal reflection *)
apply[s_List, beta] := Transpose[s];        (* Diagonal reflection *)
apply[s_List, g_List] := Fold[apply, s, g]; (* Group composition *)
group = FoldList[Append, {}, Riffle[ConstantArray[alpha, 4], beta]];


(Riffle and ConstantArray are just a trick, specific to Coxeter groups, to generate a unique set of group elements. If you don't want to allow reflections, then eliminate the four orientation-reversing group elements, which are the ones using odd numbers of generators.)

Apply the group action to the tiles and eliminate any duplicates:

tiles = Union[Flatten[Outer[apply[#1, #2] &, tileSet, group, 1], 1]];
ArrayPlot[#, Frame -> False, ImageSize -> 10 Max[Length[First[#]]]] & /@ tiles


These tiles need to be converted from arrays (suitable for graphical display) into vectors (for input to the linear program). Because these vectors represent all possible positions (translations) of the oriented tiles into the square, here is where we specify the square's dimensions, $m$ rows by $n$ columns. I use $6$ by $9$ for this example (which is a significant size). makeAllTiles creates all possible positions of a single tile within a rectangle; mapping it over the set of all oriented tiles and flattening the result gives what we want. (No need any more to remove duplicates, because there is no possibility of an accidental match.)

makeAllTiles[tile_, m_Integer, n_Integer] :=
With[{m0 = Length[tile], n0 = Length[First[tile]]},
Flatten[Table[ArrayPad[tile, {{i, m - m0 - i}, {j, n - n0 - j}}],
{i, 0, m - m0}, {j, 0, n - n0}], 1]];
allTiles = Flatten[makeAllTiles[#, 6, 9] & /@ tiles, 1];
ArrayPlot[options = Transpose[Flatten /@ allTiles], ImageSize -> 800]


Each possible position of each possible orientation of each tile shape is represented by a column vector in this array. The rows correspond to cells in the rectangle.

It's now surprisingly straightforward to complete the setup: we seek a subset of columns summing to the vector representing a filled square: all the coefficients in this sum must be $1$'s. The subset can be picked out by an indicator vector $x$. We therefore require its entries to be integral, bounded below by $0$, and bounded above by $1$. A feasible vector, let's say, will be one where the corresponding sum fills the square, but sometimes fills its cells multiple times (due to overlap of tiles): that is, a lower bound $b$ for the result should be a $1$ in each cell of the square. Finally, we minimize the sum of the cell values: if a solution exists, this sum will equal the number of squares. The sum is a linear combination with unit coefficients, represented as a vector $c$.

b = ConstantArray[1, Length[options]]; (* Lower bounds on cell values *)
c = ConstantArray[1, Length[First[options]]]; (* Sum of cell values *)
lu = ConstantArray[{0, 1}, Length[First[options]]]; (* 0-1 constraints *)
x = Quiet[LinearProgramming[c, options, b, lu, Integers, Tolerance -> 0.0005],
{LinearProgramming::lpip, LinearProgramming::lpsnf}];
If[!ListQ[x] || Max[options.x] > 1, x = {}];
solution = allTiles[[Select[x Range[Length[x]], # > 0 &]]];
ArrayPlot[z = Sum[i solution[[i]], {i, 1, Length[solution]}],
PlotRange -> {0, Max[z]}, Mesh -> True, ColorFunction -> "Rainbow"]


Solution is a list of tiles (in their orientations and positions).

Limitations

LinearProgramming (apparently) can return only one solution even when many are optimal. To explore all solutions, or a range of solutions, may require different methods. (There are still some tricks we can play: e.g., by permuting the order of columns in Options, we might wind up with a different solution.)

As with the systematic searching methods, when this problem gets big, the computing time grows rapidly. You can see this is so, just by considering the quadratic growth in the list of possible positions of the tiles as the square's dimensions increase. Computing time for the example is around one second (increasing as Tolerance decreases).

Finally, I had to fudge something to get this solution: LinearProgramming, with its default setting for Tolerance, returned a non-optimal solution! (options.x had some 2's in it, indicating some double overlaps.) By starting with Tolerance set to $0.5$ and gradually decreasing it, within a few steps a valid solution was obtained. (With the value set to 10^-1, though, I raised an Assert in the underlying C code :-(. That is a bug, at least in version 8.0.0.0.)

-
Wow! I love answers like this. It's humbling and enlightening all at once. – Mr.Wizard Jun 16 '12 at 15:29
Great answer! Perhaps your dihedral groups mastering may help me in this question on math.se related to this question in this site – belisarius has settled Jun 16 '12 at 15:35
@belisarius Actually, I had been using these methods to study your peg solitaire game (and I have even simplified it using exactly the same LinearProgramming treatment described here). – whuber Jun 16 '12 at 18:43
I think the problem is that you are trying to minimize the number of pieces used to cover the square when you use c={1,1,1,1...}. If instead you tried to minimize the number of total squares by using c=Total@options it seems to work fine without the Tolerance specification. What do you think? – Rojolalalalalalalalalalalalala Aug 8 '12 at 1:51
@Rojo That's a good idea. What I was really seeking was any feasible solution. I could not find an option in LinearProgramming to do that and so, as a work-around, was forced to find a feasible solution by optimizing something. – whuber Aug 8 '12 at 13:01

So here's a clarification to the Mathematica Journal post that David Carraher has corrected and updated so it works in Mathematica 9 posted here. Now we all have working code for this cool Mathematica solution!

-