On bounded Wang tilings

Marek Tyburec[Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 16000 Prague 6, Czech Republic] Jan Zeman

Wang tiles enable efficient pattern compression while avoiding the periodicity in tile distribution via programmable matching rules. However, most research in Wang tilings has considered tiling the infinite plane. Motivated by emerging applications in materials engineering, we consider the bounded version of the tiling problem and offer four integer programming formulations to construct valid or nearly-valid Wang tilings: a decision, maximum-rectangular tiling, maximum cover, and maximum adjacency constraint satisfaction formulations. To facilitate a finer control over the resulting tilings, we extend these programs with tile-based, color-based, packing, and variable-sized periodic constraints. Furthermore, we introduce an efficient heuristic algorithm for the maximum-cover variant based on the shortest path search in directed acyclic graphs and derive simple modifications to provide a $1/2$ approximation guarantee for arbitrary tile sets, and a $2/3$ guarantee for tile sets with cyclic transducers. Finally, we benchmark the performance of the integer programming formulations and of the heuristic algorithms showing that the heuristics provides very competitive outputs in a fraction of time. As a by-product, we reveal errors in two well-known aperiodic tile sets: the Knuth tile set contains a tile unusable in two-way infinite tilings, and the Lagae corner tile set is not aperiodic.

§ INTRODUCTION

Wang tiles, non-rotatable unit squares with colored edges, constitute a formalism introduced by Wang [1] to visualize the $\forall \exists \forall$ decidability problem of predicate calculus. Formulating an equivalent domino problem, Wang considered an infinite number of copies of an arbitrary set of Wang tiles and investigated whether there exists a simply-connected valid tiling of the infinite plane. Moreover, he conjectured in [2] that only the tile sets that form a torus, i.e., cover a periodic simply-connected rectangular domain with identical coloring of the opposite edges, generate infinite valid tilings. Berger [3] disproved the conjecture by finding a tile set that covers the infinite plane aperiodically by exploiting Kahr's reduction of the Turing halting problem [4, 5] to the origin-constrained domino problem [6]. Hence, the domino problem was proven to be undecidable and, consequently, no general finite algorithm for producing infinite valid tilings exists.

Far less attention has been paid to the finite version of the domino problem, bounded tiling, i.e., searching for a fixed-sized valid tiling generated by an arbitrary tile set. Although the problem is known to be -complete in general, e.g., [7] or <cit.>, most of the available approaches exploit specific properties of particular tile sets, e.g., [9, 10, 11, 12]. However, several closely related works address the tile packing problem for edge-matching puzzles, in which all tiles from the set are placed exactly once; see, e.g., [13, 14, 15] and [16] for an approach aiming at the famous Eternity II puzzle.

In this work, we investigate the bounded Wang tiling problem in its full generality. To this goal, we first survey the most significant aperiodic tile sets in Section <ref> and applications of Wang tiles in Section <ref>. In Section <ref>, we list available algorithms for generation of Wang tilings. Finally, our aims and contributions appear summarized in Section <ref>.

§.§ Aperiodic tile sets

[timeline color=black, line offset=2pt]
1964     $20{\small,}426$ and $104$ Wang tiles by Berger [17, 3]1966     $40$ Wang tiles over $16$ colors by Läuchli [18]1967     $104$ and $52$ Wang tiles over $8$ colors by Robinson [19, 20]1968     $92$ Wang tiles over $26$ colors by Knuth [21] (reducible to $86$ tiles [22])1971     $56$ Wang tiles over $12$ colors and $35$ tiles by Robinson [23]1973     $34$ Wang tiles by Penrose [24]1973     $32$ Wang tiles over $16$ colors by Robinson [24]1977     $24$ Wang tiles over $24$ colors by Robinson [24]1978     $16$ Wang tiles over $6$ colors by Ammann [25]1987     $24$ Wang tiles over $9$ colors by Grünbaum and Shephard [24]1995     $64$ Wang tiles by Senechal [26]1996     $14$ Wang tiles over $6$ colors by Kari [27]1996     $13$ Wang tiles over $5$ colors by Čulík [28]1999     $696$ Wang tiles by Kari and Papasoglu [29] (deterministic tile set)2008     $104$ Wang tiles by Ollinger [12]2018     $19$ Wang tiles over $16$ colors by Labbé [30] (self-similar tile set) 2021     $11$ Wang tiles over $5$ and $4$ colors by Jeandel and Rao [31]
List of aperiodic Wang tile sets.

The originally unexpected property of Wang tile sets—aperiodicity—resulted in a long-term competition among scientists in mathematical logic, computer science, discrete mathematics, and even recreational mathematicians to find the aperiodic tile set of the minimum cardinality <cit.>. Starting from the Berger tile set containing $20{\small,}426$ tiles in $1964$ [17, 3], it took almost $50$ years until the two sets of $11$ tiles were found and proved to be minimal [31]; see Fig. <ref> for a graphical overview of the selected historical developments.

In 1966, Läuchli sent to Wang an aperiodic set of $40$ tiles over $16$ colors, but it remained unpublished until 1975 [18]. Meanwhile, unaware of the Läuchli's result, Knuth [21] simplified Berger's set to $92$ tiles over $26$ colors; and Robinson developed sets of $104$ and $52$ tiles over $8$ colors in 1967 [20], of $56$ tiles over $12$ colors in 1971 [23], and anticipated an existence of a set of $35$ tiles [23].

In 1973, Penrose developed a new approach based on kite and dart tiling, leading to a set of $34$ tiles [24]. Robinson, being in contact with Penrose, modified Penrose's approach to reach a reduced set of $32$ tiles over $16$ colors [24]. Using the same technique together with Penrose rhombs tiling, Grünbaum and Shephard [24] obtained a set of $24$ tiles over $9$ colors in 1987.

Another two tile sets were discovered by Ammann. In 1978, he used the so-called Ammann bars to reach $16$ tiles over $6$ colors [25]. Building on the Ammann's A2 tiling, see, e.g., [24], Robinson obtained a set of $24$ tiles over $24$ colors in 1977 [24].

Subsequent size reduction of the smallest aperiodic set occurred in 1996, when Kari [27] developed a new method based on Mealy machines multiplying Beatty sequences and presented a set of $14$ tiles over $6$ colors. Čulík [28], using the same approach, reduced the set even further to $13$ tiles over $5$ colors.

The search for the minimal aperiodic set was concluded by Jeandel and Rao [31], who used an enumeration approach to find aperiodic sets of $11$ tiles over $4$ and $5$ colors and proved non-existence of an aperiodic set either containing $10$ or fewer tiles or labeled by less than $4$ colors.

In addition to the original Wang tiles, in 2006, Lagae and Dutré [32] described a subset of the Wang tiles, the corner tiles (we refer to Appendix <ref> for their relation to edge-based Wang tiles), with the adjacency rules stored in the colored corners instead of the edges. In the same year, they constructed multiple aperiodic sets of corner tiles [33], out of which the set of $44$ corner tiles over $6$ colors was the smallest one. The set was further simplified by Nurmi [34] to $30$ corner tiles over $6$ colors and both were claimed to be aperiodic.

§.§ Applications of Wang Tiles
Thanks to the property of particular tile sets to generate aperiodic tilings, Wang tiles gained interest among several disciplines. Building on the original purpose of Wang tiles, proofs in the first-order logic [2], they were also used in cellular automata theory [35], topology, group theory [36], and symbolic dynamical systems [37].

In computer graphics, Stam [38] adopted Wang tiles to generate aperiodic textures. After Cohen et al. [9] recognized that stochastic nonperiodic tilings are easier to handle, the Wang-tile-based approach to generating seamless textures became popular, also including the generation of point patterns [9, 39].

In science, Wang tiles and other related aperiodic tilings served as the key tool for understanding the $5$-fold symmetry of electron diffraction patterns of quasicrystals [40, 26]. Another application at the nanoscale involved molecular DNA-based realization of Wang tiles, introduced by Winfree et al. [41], which provided a self-assembly of biological nanostructures into aperiodic patterns. The self-assembly process of DNA Wang tiles also powered custom DNA-based computations [42], fueled by Turing completeness of Wang tiles [3, 18].

Beyond the nanoscale, Wang tiles have also been used for efficient compression [43] and reconstruction [44] of nonperiodic microstructures, speeding up numerical analyses of random heterogeneous materials [45, 46]. Furthermore, we have developed a bilevel optimization approach to design modular truss materials based on the corner Wang tiling formalism [47] and a clustering-based method for designing modular structures and mechanisms with continuum topology optimization [48]. Finally, Jílek et al. [49] developed a centimeter-scale self-assembly procedure building on the Wang tiling formalism.

§.§ Wang tiling generation algorithms

To the best of our knowledge, no general approaches to solving the bounded tiling problem have been reported in the literature; the only available results are specific to single families of tile sets [9, 10, 11, 12], or consider infinite thin strips [31]. In what follows, we describe the gist of three tiling algorithms: substitution-based, stochastic, and transducer-based.

Substitution-based tiling algorithm

Given a tile set $\mathcal{T}$, substitution is a map $S: \mathcal{T} \mapsto \mathfrak{T}$ that assigns a tiling $\mathfrak{T}_k$ to each tile $k \in \mathcal{T}$; we refer the reader to Section <ref> for the definitions. Consequently, arbitrary-sized tilings follow from a placed tile $k$ and a recursively applied substitution rule [12]. Hence, the tiling “grows” iteratively. Clearly, such a procedure has a low complexity, but only very specific tile sets allow for such substitution rule that generates valid tilings.

Stochastic tiling algorithms

In computer graphics, Wang tiles have mostly been used for generating visually appealing yet compressed textures. For this, it is essential to generate these nonperiodic patterns quickly, which is best achieved with stochastic tile sets—usually containing all combinations of edge labels for a given number of colors. For example, in the stochastic tiling algorithm [9], the tiling is generated row-wise, such that the neighbor of any tile that has already been placed can always be selected from at least two tiles at random. This approach was further extended towards the hash-based direct stochastic tiling algorithm [11]. Note that stochastic algorithms enable straightforward enforcement of several tile- or edge-based constraints.

Transducer-based tiling algorithm

The transducer-based tiling algorithm [31] builds on the fact that the $1$D domino problem is decidable and can be solved in a polynomial time because the bi-infinite path is formed by an arbitrary cycle in transducer graphs, see Section <ref> for clarification. To generate valid tilings of multiple rows, the product of several transducers must be computed. Hence, we must enumerate all feasible valid tilings for a requested height and unit width, and then find a path of the given length in the transducer graph of the just-formed tile set. Obviously, this approach works well for tiling thin strips; however, it is impractical for larger nearly-square domains.

§.§ Aims and novelty

In this contribution, we consider the bounded Wang tiling in its general form, thereby allowing arbitrary tile sets and control over the resulting tilings. As follows from the above state-of-the-art survey, no such method has been published yet.

We believe that development of such algorithms is important from multiple reasons. First, we have already investigated modeling and optimization of non-periodic and stochastic microstructures with Wang tilings, e.g., [43, 44, 45, 50, 46, 47]. We hope that the extension of our methods to more general tile sets would enable characterizing a broader class of non-periodic conventional materials and meta-materials [51, 52, 53] and thus improve upon the performance of optimized designs.

Apart from emerging applications in materials engineering, we believe that developing a unified methodology is of independent interest, e.g., for the verification of the results available in the literature. Here we justify this claim by finding two errors in well-established aperiodic tile sets.

To do this, we first provide the necessary definitions in Section <ref> to make the manuscript self-contained. The subsequent part is devoted to four integer programming formulations for generation of valid tilings: decision variant in Section <ref>, maximum rectangular valid tiling in Section <ref>, maximum-cover in Section <ref>, and maximum adjacency constraint satisfaction in Section <ref>. To allow for a finer control over the resulting tilings, we also include simple extensions to prescribe tile- and color-based boundary conditions, periodic constraints, and the tile-packing constraint in Section <ref>.

Due to the complexity of the proposed formulations, in Section <ref> we propose a heuristic graph-based algorithm to tackle the maximum-cover optimization variant from Section <ref>. The developed algorithm relies on solutions to shortest path problems in directed acyclic graphs, hence possesses a low asymptotic complexity. Further, we show that a slight modification maintains an approximation ratio of $2/3$ for the tile sets whose transducer graphs are cyclic.

Section <ref> collects results from the computational assessment of the integer programming formulations (Section <ref>) and heuristics (Section <ref>), and on the benchmarking of the periodic tile packing formulation against the algorithm of Lagae and Dutré [14] (Section <ref>). We close the section with two surprising observations found with integer programming for two well-known aperiodic tile sets: the Knuth [21] tile set of 92 tiles contains a tile unusable in infinite simply-connected valid tilings, Section <ref>, and the Lagae et al. [33] tile set of 44 corner tiles is not aperiodic, Section <ref>. We summarize our results in Section <ref>.

§ NOTATION AND PRELIMINARIES

Considering a finite set of color codes $\mathcal{C} = \{1, 2,\dots, n_\mathrm{c}\} \subset \mathbb{N}$, the (Wang) tile $k$ is a quadruple of the color codes $(c^\mathrm{n}_{k}, c^\mathrm{w}_{k}, c^\mathrm{s}_{k}, c^\mathrm{e}_{k})$, with $c^\mathrm{n}_{k}, c^\mathrm{w}_{k}, c^\mathrm{s}_{k}$, and $c^\mathrm{e}_{k} \in \mathcal{C}$ standing for the color codes of the north, west, south, and east edge of the tile $k$, respectively. Tiles can, therefore, be represented graphically as non-rotatable squares shown in Fig. <ref>. Without loss of generality, we further consider these squares to be of the unit size.

\begin{equation*}
(c_k^\mathrm{n}, c_k^\mathrm{w}, c_k^\mathrm{s}, c_k^\mathrm{e})\; \Leftrightarrow\; \WangTile{$c_k^\mathrm{n}$}{$c_k^\mathrm{w}$}{$c_k^\mathrm{s}$}{$c_k^\mathrm{e}$}
\end{equation*}

\begin{equation*}
\mathcal{T} =
\left\{
\WangTile{$c_1^\mathrm{n}$}{$c_1^\mathrm{w}$}{$c_1^\mathrm{s}$}{$c_1^\mathrm{e}$},\;
\WangTile{$c_2^\mathrm{n}$}{$c_2^\mathrm{w}$}{$c_2^\mathrm{s}$}{$c_2^\mathrm{e}$},\;
\dots,\;
\WangTile{$c_{n_\mathrm{t}}^\mathrm{n}$}{$c_{n_\mathrm{t}}^\mathrm{w}$}{$c_{n_\mathrm{t}}^\mathrm{s}$}{$c_{n_\mathrm{t}}^\mathrm{e}$}
\right\}
\end{equation*}

Graphical representation of (a) a Wang tile $k$, and of (b) a tile set $\mathcal{T}$.

A tile set $\mathcal{T}$ represents a finite collection of $n_\mathrm{t}$ tiles, see Fig. <ref>. When $\forall (c^\mathrm{n},c^\mathrm{w},c^\mathrm{s},c^\mathrm{e}) \in \mathcal{C}^4: (c^\mathrm{n},c^\mathrm{w},c^\mathrm{s},c^\mathrm{e}) \in \mathcal{T}$, we call the tile set complete.

Using the notation $\tilde{\bullet} = \bullet \cap \mathbb{Z}^2$ to denote an intersection of the set $\bullet$ with the integer lattice points, tiling $\mathfrak{T}^\mathcal{A}$ of a bounded domain $\mathcal{A} \in \mathbb{R}^2$ is an arrangement of copies of the tiles from $\mathcal{T}$ such that the tiles are placed at $\tilde{\mathcal{A}}$, and cover the entire domain $\mathcal{A}$, cf. Fig. <ref>. More formally, tiling is a mapping $\mathfrak{T}^\mathcal{A}: \tilde{\mathcal{A}} \rightarrow \mathcal{T}$ assigning a single tile $k \in \mathcal{T}$ to every coordinate $\left(i,j\right) \in \tilde{\mathcal{A}}$. Consequently, we call tilings $\mathfrak{T}^\mathcal{A}$ simply connected iff the domain $\mathcal{A}$ is so.

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0,0) grid (3.25,-2.25);
[step=1,black,xshift=0,yshift=0,very thin,dashed] (3.75,0) grid (4.0,-2.25);
[step=1,black,xshift=0,yshift=0,very thin,dashed] (0,-3) grid (3.25,-2.75);
[step=1,black,xshift=0,yshift=0,very thin,dashed] (3.75,-3) grid (4,-2.75);
[decoration=markings,mark=at position 1 with
[scale=1.5,>=stealth]>,postaction=decorate] (3.75,0)–(4.5,0) node[above right]$\mathcal{W}$;
[decoration=markings,mark=at position 1 with
[scale=1.5,>=stealth]>,postaction=decorate] (0,-2.75)–(0,-3.5) node[left]$\mathcal{H}$;
[label=$1$] at (0,0) ;
[label=$2$] at (1,0) ;
[label=$3$] at (2,0) ;
[label=$4$] at (3,0) ;
[label=$n_\mathrm{w}$] at (4,0) ;
[label=left:$1$] at (0,0) ;
[label=left:$2$] at (0,-1) ;
[label=left:$3$] at (0,-2) ;
[label=left:$n_\mathrm{h}$] at (0,-3) ;

Color matching among tiles $p$, $q$, and $r \in \mathcal{T}$.

The tiling $\mathfrak{T}^\mathcal{A}$ is rectangular if $\forall i \in \mathcal{H}, \mathcal{H} = \{1,\dots, n_\mathrm{h}\},$ and $\forall j \in \mathcal{W}, \mathcal{W} = \{1,\dots, n_\mathrm{w}\}$, it holds that $(i,j) \in \tilde{\mathcal{A}}$. Here, $\mathcal{H}$ and $\mathcal{W}$ are the sets of the height and width coordinates.

A valid tiling (Wang tiling) of $\mathcal{A}$, denoted by $\mathfrak{T}^\mathcal{A}_\mathrm{valid}$, is a tiling with equal color codes at the shared edges between all pairs of adjoining tiles. Therefore, the mapping $\mathfrak{T}^\mathcal{A}_\mathrm{valid}: \tilde{\mathcal{A}} \rightarrow \mathcal{T}$ satisfies, in addition to the requirements for $\mathfrak{T}^\mathcal{A}$, the additional constraints

\begin{equation}\label{eq:southNorthMapping}
c^\mathrm{s}_{\mathfrak{T}^\mathcal{A}_\mathrm{valid}(i,j)} = c^\mathrm{n}_{\mathfrak{T}^\mathcal{A}_\mathrm{valid}(i+1,j)}, \quad \forall (i,j) \in \tilde{\mathcal{A}}: (i+1,j) \in \tilde{\mathcal{A}},
\end{equation}

\begin{equation}\label{eq:eastWestMapping}
c^\mathrm{e}_{\mathfrak{T}^\mathcal{A}_\mathrm{valid}(i,j)} = c^\mathrm{w}_{\mathfrak{T}^\mathcal{A}_\mathrm{v}(i,j+1)}, \quad \forall (i,j) \in \tilde{\mathcal{A}}: (i,j+1) \in \tilde{\mathcal{A}},
\end{equation}

provided that the axes are oriented accordingly to Fig. <ref>. If such $\mathfrak{T}^\mathcal{A}_\mathrm{valid}$ exists, we say that the domain $\mathcal{A}$ admits a valid $\mathcal{T}$-tiling, or that it is tileable by $\mathcal{T}$.

Consider that $\mathcal{B}\subseteq \mathcal{A}$ and $\mathcal{B_{\max \mathrm{rect}}} \subseteq \mathcal{A}$ are simply connected, rectangular, and $\mathcal{T}$-tileable. Then, the maximum rectangular valid tiling $\mathfrak{T}^\mathcal{A}_{\mathrm{v},\max \mathrm{rect}}$ is a valid tiling of the domain $\mathcal{B}_{\max \mathrm{rect}}$, where $\{\mathcal{B}_{\max \mathrm{rect}} \subseteq \mathcal{A}, \forall \mathcal{B} \subseteq \mathcal{A}: \lvert\tilde{\mathcal{B}}_{\max \mathrm{rect}} \rvert \ge \lvert \tilde{\mathcal{B}} \rvert\}$. Here, the notation $\lvert \bullet \rvert$ denotes cardinality of the set $\bullet$.

The maximum cover $\mathfrak{T}^\mathcal{A}_\mathrm{v, {\max \mathrm{cov}}}$ is a valid tiling of $\mathcal{B_{\max \mathrm{cov}}}$, where $\mathcal{B}$ and $\mathcal{B}_{\max \mathrm{cov}}$ are arbitrary $\mathcal{T}$-tileable subdomains in $\mathcal{A}$ and $\{\mathcal{B}_{\max \mathrm{cov}} \subseteq \mathcal{A}, \forall \mathcal{B} \subseteq \mathcal{A}: \lvert\tilde{\mathcal{B}}_{\max \mathrm{cov}} \rvert \ge \lvert \tilde{\mathcal{B}} \rvert\}$.

A rectangular valid tiling is said to be periodic, if the color codes at the opposite sides of the rectangle match. If the valid tiling is not periodic, but the considered tile set allows for at least one periodic rectangular tiling, we call it nonperiodic. Finally, if no such periodic pattern exists and the tile set still allows for a valid tiling of the infinite plane, it is referred to as aperiodic. Similarly, the tile set $\mathcal{T}$ is periodic if it permits periodic valid tilings; and aperiodic if all feasible valid tilings are aperiodic.

Transducer graph[27] $G_\mathrm{t,h}$ of the tile set $\mathcal{T}$ is a directed (multi-)graph representation of a Mealy machine without any initial nor terminal state. It consists of $\lvert \mathcal{C}\rvert$ states (graph vertices) and $\lvert \mathcal{T} \rvert$ transitions (directed edges) $\mathcal{E}_\mathrm{h}$, where

\begin{equation}
\mathcal{E}_\mathrm{h} \coloneqq \bigcup_{k \in \mathcal{T}}\left( c_k^\mathrm{w} \xrightarrow{c_k^\mathrm{s}\vert c_k^\mathrm{n}} c_k^\mathrm{e}\right).
\end{equation}

For the dual transducer graph $G_\mathrm{t,v}$, composed of the dual Wang tiles [30] reflecting $\mathcal{T}$ along the major diagonal of the tiles, the edge set is defined as

\begin{equation}
\mathcal{E}_\mathrm{v} \coloneqq \bigcup_{k \in \mathcal{T}}\left( c_k^\mathrm{n} \xrightarrow{c_k^\mathrm{e}\vert c_k^\mathrm{w}} c_k^\mathrm{s}\right).
\end{equation}

To illustrate the construction, we include a visual example in Fig. <ref>.

6mm$\mathcal{T} = \left\{\WangTile{0}{1}{1}{0},
\WangTile{0}{1}{0}{1},
\WangTile{1}{0}{0}{1}
\right\}$

[>=stealth',shorten >=1pt,auto,node distance=2cm, semithick,every node/.style=inner sep=2pt, every state/.style=minimum size=0pt,inner sep=2pt]
(0) [state] $0$;
(1) [state, right=of 0] $1$;
[->] (0) edge [bend left=40, above, sloped, rotate=180] node(b) $0\vert1$ (1);
[->] (1) edge [bend left=40, above, sloped, rotate=180] node(b) $1\vert0$ (0);
[->] (1)  edge [loop above, sloped] node $0\vert0$ (1);

[>=stealth',shorten >=1pt,auto,node distance=2cm, semithick,every node/.style=inner sep=2pt, every state/.style=minimum size=0pt,inner sep=2pt]
(2) [state] $0$;
(3) [state, right=of 0] $1$;
[->] (0) edge [bend left=40, above, sloped, rotate=180] node(b) $0\vert1$ (1);
[->] (1) edge [bend left=40, above, sloped, rotate=180] node(b) $1\vert0$ (0);
[->] (0)  edge [loop above, sloped] node $1\vert1$ (0);

(b) Transducer and (c) dual transducer graphs of the tile set (a).

§ INTEGER PROGRAMMING FORMULATIONS

In this section, we introduce four integer programming formulations for the generation of valid tilings. The first one, in Section <ref>, develops a decision variant. In the later sections, we investigate the maximum rectangular tiling (Section <ref>), maximum cover (Section <ref>), and the maximum adjacency constraints satisfaction (Section <ref>). Finally, Section <ref> proposes several extensions to facilitate finer control over the resulting tilings.

§.§ Rectangular valid tiling

Let us now consider the fundamental problem of finding $\mathfrak{T}^\mathcal{A}_\mathrm{v}$ or proving it does not exist. From now on, we restrict $\mathcal{A}$ to be rectangular to simplify notation. However, the presented approach also extends to the general case.

To achieve this, we introduce $\forall (i,j,k) \in \mathcal{H}\times \mathcal{W}\times\mathcal{T}$ a binary decision variable $x_{i,j,k} \in \{0,1\}$ denoting the placement of the tile $k$ at the $(i,j)$ coordinate such that

\begin{equation}\label{eq:integrality}
x_{i,j,k}  = \left\{
\begin{array}{c l}
1 & \text{iff the tile $k$ lies at coordinate $(i,j)$,}\\
0 & \text{otherwise.}
\end{array}\right.
\end{equation}

Consequently, mapping $\mathfrak{T}^\mathcal{A}(i,j)$ is expressed as

\begin{equation}\label{eq:mapping}
\mathfrak{T}^\mathcal{A}(i,j) = \sum_{k \in \mathcal{T}} k x_{i,j,k},
\end{equation}

together with the requirement that every $(i,j)$ coordinate is occupied by one tile,

\begin{equation}\label{eq:sumTiles}
\sum_{k \in \mathcal{T}} x_{i,j,k} = 1,\quad \forall (i,j) \in \mathcal{H} \times \mathcal{W}.
\end{equation}

Similarly, the color codes of a tile placed at $(i,j)$ are expressed using the binary variables as

\begin{equation}\label{eq:northMap}
c_{\mathfrak{T}^\mathcal{A}(i,j)}^\mathrm{n} = \sum_{k \in \mathcal{T}} c_k^\mathrm{n} x_{i,j,k},
\end{equation}

\begin{equation}\label{eq:westMap}
c_{\mathfrak{T}^\mathcal{A}(i,j)}^\mathrm{w} = \sum_{k \in \mathcal{T}} c_k^\mathrm{w} x_{i,j,k},
\end{equation}

\begin{equation}\label{eq:southMap}
c_{\mathfrak{T}^\mathcal{A}(i,j)}^\mathrm{s} = \sum_{k \in \mathcal{T}} c_k^\mathrm{s} x_{i,j,k},
\end{equation}

\begin{equation}\label{eq:eastMap}
c_{\mathfrak{T}^\mathcal{A}(i,j)}^\mathrm{e} = \sum_{k \in \mathcal{T}} c_k^\mathrm{e} x_{i,j,k}.
\end{equation}

Inserting (<ref>) into (<ref>) and (<ref>) leads to the horizontal and vertical adjacency constraints expressed in terms of the decision variables, as

\begin{equation}
\sum_{k \in \mathcal{T}} c_k^\mathrm{s} x_{i,j,k} - \sum_{k \in \mathcal{T}} c_k^\mathrm{n} x_{i+1,j,k} = 0, \quad \forall (i,j) \in \mathcal{H}\setminus \{n_\mathrm{h}\}\times \mathcal{W},
\end{equation}

\begin{equation}
\sum_{k \in \mathcal{T}} c_k^\mathrm{e} x_{i,j,k} - \sum_{k \in \mathcal{T}} c_k^\mathrm{w} x_{i,j+1,k} = 0, \quad \forall (i,j) \in \mathcal{H}\times \mathcal{W}\setminus \{n_\mathrm{w}\}.
\end{equation}

Combining (<ref>), (<ref>), (<ref>), and (<ref>) then provides us with a complete binary linear programming representation of valid tiling $\mathfrak{T}^\mathcal{A}_\mathrm{valid}$.

For computational reasons, it proved to be advantageous to organize the constraints according to the color codes:

\begin{align}
\sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{s}=\ell] - \sum_{k \in \mathcal{T}} x_{i+1, j, k}[c_k^\mathrm{n}=\ell] = 0,\quad \forall (i,j,\ell) \in \mathcal{H}\setminus \{n_\mathrm{h}\} \times \mathcal{W}\times \mathcal{C},\label{eq:colorCona}\\
\sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{e}=\ell] - \sum_{k \in \mathcal{T}} x_{i, j+1, k}[c_k^\mathrm{w}=\ell] = 0,\quad \forall (i,j,\ell) \in \mathcal{H} \times \mathcal{W}\setminus \{n_\mathrm{w}\} \times \mathcal{C}\label{eq:colorConb},
\end{align}

where, in the Iverson notation [54], $\sum_{k \in \mathcal{T}} x_{i,j,k} \left[ c_k^\mathrm{s} = \ell \right]$ expresses that $x_{i,j,k}$ is added to the sum if and only if $c_k^\mathrm{s} = \ell$.

The constraint (<ref>) requires that the number of tiles at $(i,j)$ with the south edge colored by $\ell$ equals to the number of tiles at $(i+1,j)$ with the north edge marked by the same $\ell$, for all $\ell \in \mathcal{C}$. Because of (<ref>), there are either no tiles with the shared edge colored by $\ell$, or there is a single tile at each of the coordinates with its common edge labeled by $\ell$. Analogously to the vertical adjacency constraint, the horizontal constraint (<ref>) also enforces equality among the number of tiles at $(i,j)$ with the east edge colored by $\ell$ and the number of tiles at $(i,j+1)$ having the west edge colored by identical $\ell$.

Finally, combining (<ref>), (<ref>), and (<ref>), while noticing that the constraints (<ref>) naturally propagate with the adjacency constraints from the domain boundaries (compare (<ref>,<ref> with (<ref>)), leads to the binary programming formulation

find      𝐱

∑_k ∈𝒯 x_i,j,k[c_k^s=ℓ] - ∑_k ∈𝒯 x_i+1, j, k[c_k^n=ℓ] = 0,  ∀(i,j, ℓ) ∈ℋ∖{n_h}×𝒲×𝒞,

∑_k ∈𝒯 x_i,j,k[c_k^e=ℓ] - ∑_k ∈𝒯 x_i, j+1, k[c_k^w=ℓ] = 0,  ∀(i,j,ℓ) ∈ℋ ×𝒲∖{n_w} ×𝒞,

∑_k ∈𝒯 x_i,j,k = 1,
∀(i,j) ∈{1, n_h} ×𝒲,

∑_k ∈𝒯 x_i,j,k = 1,
∀(i,j) ∈ℋ ×{1, n_w},

x_i,j,k ∈{0,1},
∀(i,j,k) ∈ℋ ×𝒲×𝒯,

that provides a complete representation of the bounded tiling problem, i.e., all valid tilings solve the integer program, and conversely, all feasible solutions to (<ref>) are valid tilings. Moreover, observe that the problem consists of two totally unimodular constraints if considered independently: (<ref>,<ref>) representing row tilings, and (<ref>,<ref>) being column tilings. When considered simultaneously, the resulting problem becomes -complete [7, 8].

§.§ Maximum rectangular valid tiling

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (1.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (1.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (1.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (1.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-0.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-0.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-0.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-0.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;

Admissible tile placements (a)–(c) and (e)–(g), and forbidden placements (d) and (h) in the maximum rectangular valid tiling formulation.

When a solution to (<ref>) cannot be found in an acceptable time period or when no such solution exists, one can resort to relaxing the requirement of a valid tiling of $\mathcal{A}$ and search for a valid tiling of the largest rectangular subdomain.

Without loss of generality, let us assume that the maximum rectangular valid tiling always contains an anchor tile placed at $(1,1)$, i.e.,

\begin{equation}
\sum_{k \in \mathcal{T}}  x_{1,1,k} = 1.
\label{eq:rectOneEq}
\end{equation}

On the other hand, all the other coordinates may contain a tile or be empty, thus

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k} \le 1,
\quad \forall (i,j) \in \tilde{\mathcal{A}} \setminus (1,1).
\label{eq:rectAtMostOne}
\end{equation}

Let us now pick two vertically adjacent coordinates $(i,j)$ and $(i+1,j)$. If there is a tile $q$ placed at $(i+1,j)$, another tile $p$ has to be placed at $(i,j)$, as, otherwise, there is no simply-connected rectangular tiling containing both the tiles at $(1,1)$ and at $(i+1,j)$. Validity of the tiling also requires identical color codes at the shared edges. On the other hand, if no tile is placed at $(i+1,j)$, a coordinate $(i,j)$ may be either occupied or empty. The allowed and forbidden combinations are shown in Fig. <ref>–<ref>. Formally stated in terms of the decision variables, these considerations are expressed as

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{s}=\ell] - \sum_{k \in \mathcal{T}} x_{i+1, j, k}[c_k^\mathrm{n}=\ell] \ge 0,\quad \forall (i,j,\ell) \in \mathcal{H}\setminus \{n_\mathrm{h}\}\times \mathcal{W}\times \mathcal{C}.
\label{eq:rectAdj1}
\end{equation}

Similar arguments hold also for the coordinates $(i,j)$ and $(i,j+1)$, resulting in the constraints

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{e}=\ell] - \sum_{k \in \mathcal{T}} x_{i, j+1, k}[c_k^\mathrm{w}=\ell] \ge 0,\quad \forall (i,j,\ell) \in \mathcal{H} \times \mathcal{W}\setminus \{n_\mathrm{w}\} \times \mathcal{C}.
\label{eq:rectAdj2}
\end{equation}

The allowed and forbidden combinations for this case are shown in Fig. <ref>–<ref>.

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

[step=1,black,xshift=0,yshift=0,very thin,dashed] (0.25,0.80) grid (2.75,-1.75);
[label=$j$] at (1,0.75) ;
[label=$j+1$] at (2,0.75) ;
[label=left:$i$] at (0.3,0) ;
[label=left:$i+1$] at (0.3,-1) ;

Six possible placements of tiles $p$, $q$, $r$, and $s$. The combination (b) cannot appear in any rectangular tiling.

The developed constraints (<ref>)–(<ref>) enforce simple connectedness; however, they do not guarantee that the resultant tiling will be rectangular. For any $4$ adjacent tiles $p$, $q$, $r$, and $s$ placed at $(i,j)$, $(i+1,j)$, $(i,j+1)$, and $(i+1,j+1)$, respectively, these constraints allow for the assemblies shown in Fig. <ref>. Because the combination <ref> cannot appear in any simply-connected rectangular tiling, we must exclude it from the feasible set,

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i+1,j,k} + \sum_{k \in \mathcal{T}} x_{i,j+1,k} - \sum_{k \in \mathcal{T}} x_{i+1,j+1,k} \le 1, \quad \forall (i,j) \in \mathcal{H} \setminus\{n_\mathrm{h}\} \times \mathcal{W} \setminus\{n_\mathrm{w}\}.
\label{eq:rectRect}
\end{equation}

Finally, combining Eqs. (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) together with an objective function to maximize $\lvert \tilde{\mathcal{B}}_{\max \mathrm{rect}}\rvert$ provides us with the binary maximum rectangular valid tiling optimization program

max_𝐱      ∑_i ∈ℋ ∑_j ∈𝒲 ∑_k ∈𝒯 x_i,j,k

∑_k ∈𝒯 x_i,j,k[c_k^s=ℓ] - ∑_k ∈𝒯 x_i+1, j, k[c_k^n=ℓ] ≥0,  ∀(i,j,ℓ) ∈ℋ∖{n_h}×𝒲×𝒞,

∑_k ∈𝒯 x_i,j,k[c_k^e=ℓ] - ∑_k ∈𝒯 x_i, j+1, k[c_k^w=ℓ] ≥0,  ∀(i,j,ℓ) ∈ℋ×𝒲∖{n_w}×𝒞,

∑_k ∈𝒯 x_i+1,j,k + ∑_k ∈𝒯 x_i,j+1,k - ∑_k ∈𝒯 x_i+1,j+1,k ≤1,   ∀(i,j) ∈ℋ ∖{n_h}×𝒲 ∖{n_w},

∑_k ∈𝒯  x_1,1,k = 1,

∑_k ∈𝒯 x_i,j,k ≤1,
∀(i,j) ∈𝒜̃∖(1,1),

x_i,j,k ∈{0,1},
∀(i,j,k) ∈ℋ×𝒲×𝒯.

In contrast to (<ref>), a feasible solution to the optimization program (<ref>) can be found in a polynomial time, e.g., by tiling the first row or column of the $1$D bounded tiling problem. However, finding an optimal solution to (<ref>) is -hard, because the optimization problem (<ref>) is reducible to the decision version (<ref>) by fixing the value of the objective function to $\lvert \tilde{\mathcal{A}} \rvert$, which enforces equalities in (<ref>), (<ref>), and (<ref>), making the constraint (<ref>) redundant as a consequence.

§.§ Maximum cover

Another option for avoiding the infeasibility of (<ref>) rests in neglecting the requirement of (simple) connectedness, hence allowing for a placement of empty tiles (voids). In this section, we therefore search the maximum cover of $\mathcal{A}$, or equivalently a valid tiling of the (possibly disconnected) domain $\mathcal{B_{\max \mathrm{cov}}} \subseteq \mathcal{A}$. For the maximum cover formulation, we assume that any two adjacent tiles satisfy the edge-matching constraints of valid tilings, but these are also satisfied by any of the tile-void, void-tile, or void-void combination, where $\sum_{k \in \mathcal{T}} x_{i,j,k} = 0$ for a void located at $(i,j) \in \tilde{\mathcal{A}}$.

Thus, each coordinate $(i,j)$ is occupied either by a tile or a void, implying that

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k} \le 1,
\quad \forall (i,j) \in \mathcal{H}\times \mathcal{W},
\label{eq:coveringAtMostOne}
\end{equation}

and the vertical and horizontal edge matching conditions become

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{e}=\ell] + \sum_{k \in \mathcal{T}} x_{i, j+1, k}[c_k^\mathrm{w} \neq \ell] \le 1,\quad \forall (i,j,\ell) \in \mathcal{H}\times \mathcal{W} \setminus\{n_{\mathrm{w}}\}\times \mathcal{C},
\label{eq:coveringAdj1}
\end{equation}

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{s}=\ell] +  \sum_{k \in \mathcal{T}} x_{i+1, j, k}[c_k^\mathrm{n}\neq \ell] \le 1,\quad \forall (i,j,\ell) \in \mathcal{H} \setminus\{n_{\mathrm{h}}\}\times \mathcal{W}\times \mathcal{C}.
\label{eq:coveringAdj2}
\end{equation}

Finally, the combination of Eqs. (<ref>), (<ref>), (<ref>) with the objective function to maximize $\lvert \tilde{\mathcal{B}}_{\max \mathrm{cov}} \rvert$ leads to the binary optimization problem

max_𝐱      ∑_i ∈ℋ ∑_j ∈𝒲 ∑_k ∈𝒯 x_i,j,k

∑_k ∈𝒯 x_i,j,k[c_k^e=ℓ] + ∑_k ∈𝒯 x_i, j+1, k[c_k^w ≠ℓ] ≤1,  ∀(i,j,ℓ) ∈ℋ ×𝒲 ∖{n_w} ×𝒞,

∑_k ∈𝒯 x_i,j,k[c_k^s=ℓ] +  ∑_k ∈𝒯 x_i+1, j, k[c_k^n≠ℓ] ≤1,  ∀(i,j,ℓ) ∈ℋ ∖{n_h} ×𝒲 ×𝒞,

∑_k ∈𝒯 x_i,j,k ≤1,
∀(i,j) ∈ℋ ×𝒲,

x_i,j,k ∈{0,1},
∀(i,j,k) ∈ℋ ×𝒲 ×𝒯.

The program (<ref>) is trivially -hard: Requiring the objective function (<ref>) to be at least $\lvert \tilde{\mathcal{A}} \rvert$ implies that

\begin{equation}
\sum_{k \in \mathcal{T}} x_{i,j,k} = 1, \quad \forall (i,j) \in \mathcal{H} \times \mathcal{W},
\end{equation}

i.e., all positions are occupied by a Wang tile. Moreover, (<ref>) and (<ref>) require all adjacent tiles to share the color codes at their common edges. Consequently, the resulting tiling is void-free and valid, and solves the -complete bounded tiling problem.

§.§ Maximum adjacency constraints satisfaction

Because the decision problem (<ref>) also constitutes a specific instance of the constraint satisfaction problem (CSP), another optimization variant comes from the formulation of the max-CSP problem, maximizing the number of satisfied clauses—color matches in our case.

Therefore, for each vertical and horizontal edge we introduce a new variable	 $h_{i,j}^\mathrm{v} \in \mathbb{R}_{\ge 0}$, where $(i,j) \in \mathcal{H} \times \mathcal{W} \setminus n_\mathrm{w}$, and $h_{i,j}^\mathrm{h} \in \mathbb{R}_{\ge 0}$, with $(i,j) \in \mathcal{H} \setminus n_\mathrm{h} \times \mathcal{W}$, respectively. The adjacency constraints (<ref>) are then relaxed by considering

\begin{align}
\left| \sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{s}=\ell] - \sum_{k \in \mathcal{T}} x_{i+1, j, k}[c_k^\mathrm{n}=\ell] \right| \le h_{i,j}^\mathrm{h},\quad \forall (i,j, \ell) \in \mathcal{H}\setminus \{n_\mathrm{h}\} \times \mathcal{W} \times \mathcal{C},\label{eq:maxCSPcolorA}\\
\left| \sum_{k \in \mathcal{T}} x_{i,j,k}[c_k^\mathrm{e}=\ell] - \sum_{k \in \mathcal{T}} x_{i, j+1, k}[c_k^\mathrm{w}=\ell] \right| \le h_{i,j}^\mathrm{v},\quad \forall (i,j,\ell) \in \mathcal{H} \times \mathcal{W}\setminus \{n_\mathrm{w}\} \times \mathcal{C} \label{eq:maxCSPcolorB}
\end{align}

instead. Indeed, if $h_{i,j}^\mathrm{h} = 0$, the edge-matching requirement of the neighboring tiles at $(i,j)$ and $(i+1,j)$ is satisfied; and it is violated otherwise. Similarly, $h_{i,j}^\mathrm{v} = 0$ guarantees color matches among the tiles at $(i,j)$ and $(i,j+1)$.

Finally, rewriting absolute values in (<ref>) by two linear inequalities while supplying an objective function to maximize the number of color matches yields the binary optimization problem

max_𝐱      ∑_i∈ℋ ∑_j ∈𝒲 ∖n_w  (1 - h_i,j^v) + ∑_i∈ℋ ∖n_h ∑_j ∈𝒲 (1 - h_i,j^h)

∑_k ∈𝒯 x_i,j,k[c_k^s=ℓ] - ∑_k ∈𝒯 x_i+1, j, k[c_k^n=ℓ] ≤h_i,j^h,   ∀(i,j,ℓ) ∈ℋ∖{n_h} ×𝒲 ×𝒞,

∑_k ∈𝒯 x_i+1, j, k[c_k^n=ℓ] -∑_k ∈𝒯 x_i,j,k[c_k^s=ℓ] ≤h_i,j^h,   ∀(i,j,ℓ) ∈ℋ∖{n_h} ×𝒲 ×𝒞,

∑_k ∈𝒯 x_i,j,k[c_k^e=ℓ] - ∑_k ∈𝒯 x_i, j+1, k[c_k^w=ℓ] ≤h_i,j^v,   ∀(i,j,ℓ) ∈ℋ ×𝒲∖{n_w} ×𝒞,

∑_k ∈𝒯 x_i, j+1, k[c_k^w=ℓ] - ∑_k ∈𝒯 x_i,j,k[c_k^e=ℓ] ≤h_i,j^v,   ∀(i,j,ℓ) ∈ℋ ×𝒲∖{n_w} ×𝒞,

∑_k ∈𝒯 x_i,j,k = 1,
∀(i,j) ∈ℋ ×𝒲,

x_i,j,k ∈{0,1},
∀(i,j,k) ∈ℋ ×𝒲 ×𝒯,

that is -hard due to the reduction to (<ref>) after setting all $h_{i,j}^\mathrm{v}$ and $h_{i,j}^\mathrm{h}$ to zeros. A feasible solution can be found in a polynomial time by finding valid row/column tilings for each row/column, so that either term $\sum_{i\in \mathcal{H}} \sum_{j \in \mathcal{W} \setminus n_\mathrm{w}} h_{i,j}^\mathrm{v}$ or $\sum_{i\in \mathcal{H} \setminus n_\mathrm{h}} \sum_{j \in \mathcal{W}} h_{i,j}^\mathrm{h}$ equals zero.

§.§ Extensions

Up to now, we have focused solely on the (re)formulations of the bounded tiling problem, searching for arbitrary valid tilings. However, some potential applications may require finer control over the resulting tilings. Thus, in this section, we state some simple extensions to enforce tile- and color-based boundary conditions to solve the tile packing problem [14] and to enforce (variable-sized) periodic boundary conditions.

§.§.§ Tile-based boundary conditions

At first, we consider boundary conditions in the form of prescribed tiles. As the simplest one, we enforce the placement of a tile $k$ at $(i,j)$:

\begin{equation}
x_{i,j,k} = 1, \quad (i,j,k) \in \mathcal{H} \times \mathcal{W} \times \mathcal{T}.
\end{equation}

Similarly, we may prevent tile $k$ from being placed there:

\begin{equation}
x_{i,j,k} = 0, \quad (i,j,k) \in \mathcal{H} \times \mathcal{W} \times \mathcal{T}.
\end{equation}

Placement of an identical tile at the coordinates $(i,j) \in \tilde{\mathcal{A}}$ and $(p,q) \in \tilde{\mathcal{A}}$ requires

\begin{equation}
x_{i,j,k} -  x_{p,q,k} = 0, \quad \{i,p\} \in \mathcal{H}, \{j,q\} \in \mathcal{W}, \forall k \in \mathcal{T}.
\end{equation}

Conversely, different tiles at these coordinates are secured with

\begin{equation}
x_{i,j,k} + x_{p,q,k} \le 1,\quad\{i,p\} \in \mathcal{H}, \{j,q\} \in \mathcal{W}, \forall k \in \mathcal{T}.
\end{equation}

§.§.§ Color-based boundary conditions

In addition to the tile-based constraints, we may also enforce specific color codes for individual edges. To do this, the color of the north edge at $(i,j) \in \tilde{\mathcal{A}}$ is set to $\ell$ by

\begin{equation}
\sum_{k \in T} x_{i,j,k} [c^\mathrm{n}_k = \ell] = 1, \quad (i,j,\ell) \in \mathcal{H} \times \mathcal{W} \times \mathcal{C}.
\end{equation}

On the contrary, we may prevent this color by requiring

\begin{equation}
\sum_{k \in T} x_{i,j,k} [c^\mathrm{n}_k = \ell] = 0, \quad (i,j,\ell) \in \mathcal{H} \times \mathcal{W} \times \mathcal{C}.
\end{equation}

Further, the same color codes at the north edge of $(i,j) \in \tilde{\mathcal{A}}$ and at the west edge of $(p,q) \in\tilde{\mathcal{A}}$ are established with

\begin{equation}
\sum_{k \in T} x_{i,j,k} [c^\mathrm{n}_k = \ell] - \sum_{k \in T} x_{p,q,k} [c^\mathrm{w}_k = \ell] = 0, \quad \{i,p\} \in \mathcal{H}, \{j,q\} \in \mathcal{W}, \forall \ell \in \mathcal{C},
\end{equation}

and a different color with

\begin{equation}
\sum_{k \in T} x_{i,j,k} [c^\mathrm{n}_k = \ell] + \sum_{k \in T} x_{p,q,k} [c^\mathrm{w}_k = \ell] \le 1, \quad \{i,p\} \in \mathcal{H}, \{j,q\} \in \mathcal{W}, \forall \ell \in \mathcal{C}.
\end{equation}

§.§.§ Periodic tiling

In the domino problem, Wang [1] investigated the existence of tile sets admitting infinite aperiodic tilings. Here, we consider a similar setting for the finite domain $\mathcal{A}$: examining periodicity through periodic color-based boundary conditions.

We begin with requiring equal coloring at the fixed opposite domain boundaries,

\begin{align}
\sum_{k \in T} x_{1,j,k} [n_k = \ell] - \sum_{k \in T} x_{n_{\mathrm{t},h},j,k} [s_k = \ell] = 0, \quad \forall (j,\ell) \in \mathcal{W} \times \mathcal{C},\\
\sum_{k \in T} x_{i,1,k} [w_k = \ell] - \sum_{k \in T} x_{i,n_{\mathrm{t},w},k} [e_k = \ell] = 0, \quad \forall (i,\ell) \in \mathcal{H} \times \mathcal{C}.
\end{align}

When adding (<ref>) to the decision problem (<ref>), we thus ask for an existence of a fixed-sized periodic Wang tiling.

In a natural generalization, we ask for an existence of finite-sized periodic Wang tilings, thus relying on the maximum rectangular valid tiling formulation (<ref>). Naturally, the domain size is not known in this case. Therefore, we must consider $\forall (i,j,\ell) \in \mathcal{H} \times \mathcal{W}\times \mathcal{C}$ constraints of the form

∑_k ∈𝒯 x_i,j,k[e_k≠ℓ] + ∑_k ∈𝒯 x_i, 1, k[w_k=ℓ] -  ∑_k ∈𝒯 x_i,j+1,k [j<n_t,w] ≤1,

∑_k ∈𝒯 x_i,j,k[s_k ≠ℓ] + ∑_k ∈𝒯 x_1, j, k[n_k=ℓ] - ∑_k ∈𝒯 x_i+1,j,k[i<n_t,h] ≤1.

Here, (<ref>) prevents a color mismatch of the north edge of $(1,j) \in \tilde{\mathcal{A}}$ and the south edge of $(i,j) \in \mathcal{A}$ iff there is no tile placed at $(i,j+1) \in \tilde{\mathcal{A}}$. Similarly, in the case of (<ref>), we prevent a color mismatch of the west edge at $(i,1) \in \tilde{\mathcal{A}}$ and the east edge at $(i,j) \in \tilde{\mathcal{A}}$ iff the position $(i+1,j) \in \tilde{\mathcal{A}}$ is empty.

Finally, when adding the constraints (<ref>) to (<ref>), we usually search for the smallest periodic pattern rather than the largest,

\begin{align}\label{eq:smallestTiling}
\min_{\mathbf{x}}\, & \sum_{i \in \mathcal{H}} \sum_{j \in \mathcal{W}} \sum_{k \in \mathcal{T}} x_{i,j,k}.
\end{align}

§.§.§ Tile packing problem

Our last extension constitutes the setting of the tile-packing problem [14]: we require each tile to be placed exactly once yet form a fixed-sized valid tiling,

\begin{equation}\label{eq:packing}
\sum_{i \in \mathcal{H}} \sum_{j \in \mathcal{W}} x_{i,j,k} = 1, \quad \forall k \in \mathcal{T}.
\end{equation}

Note here that this extension requires that $\lvert \mathcal{T} \rvert = \lvert \tilde{\mathcal{A}} \rvert$ as, otherwise, no solution exists.

§ HEURISTIC ALGORITHM FOR THE MAXIMUM COVER TILING PROBLEM

In the previous sections, we have introduced several integer programming formulations for the bounded Wang tiling problem and their extensions. Because of the asymptotic complexity of the integer programming formulations, we further develop a simple heuristic algorithm for one of the optimization variants, the maximum cover.

§.§ Maximum row cover tilings

Let us start with revising the decision program (<ref>). In this formulation, neglecting any pair of the constraints (<ref>, <ref>) or (<ref>,<ref>) provides a totally unimodular constraint matrix, recall Section <ref>. Consequently, such simplified problems are deterministically solvable using the simplex method. Moreover, this setting agrees with the maximum flow problem [55], as (<ref>) and (<ref>) correspond to the flow balances in the source and sink, whereas (<ref>) and (<ref>) correspond to the Kirchhoff law equations. Further complexity reduction is possible by recognizing the (shortest) path problem structure, since the source and sink capacities are equal to one, allowing only a single source-to-sink path with positive flow to emerge. Omitting any of these constraint pairs produces valid tilings of (finite) stripes, i.e., of rows or columns. However, the edges shared by the neighboring stripes may not comply with the edge matching rules. Starting with this observation, we first focus on an efficient approach to generate valid tilings of the rows.

[>=stealth',shorten >=1pt,auto,node distance=2cm, semithick,every node/.style=inner sep=2pt, every state/.style=minimum size=20pt,inner sep=2pt,on grid]

(s) [state] at (0,0) $s$;

(0) [state] at (2,2.25) $1$;
(1) [state] at (2,0.75) $2$;
(d) [] at (2,-0.75) $\vdots$;
(C) [state] at (2,-2.25) $\lvert\mathcal{C}\rvert$;

(0a) [state, right of=0]$1$;
(1a) [state, right of=1]$2$;
(da) [right of=d] $\vdots$;
(Ca) [state, right of=C]$\lvert\mathcal{C}\rvert$;

(0b) [right of=0a]$\dots$;
(1b) [right of=1a]$\dots$;
(db) [right of=da] $\ddots$;
(Cb) [right of=Ca]$\dots$;

(0c) [state, right of=0b]$1$;
(1c) [state, right of=1b]$2$;
(dc) [right of=db] $\vdots$;
(Cc) [state, right of=Cb]$\lvert\mathcal{C}\rvert$;

(t) [state] at (10,0)$t$;

(t1) [rotate=90] at (3,0) edges $\mathcal{E}_\mathrm{h}$;
(t2) [rotate=90] at (5,0) edges $\mathcal{E}_\mathrm{h}$;
(t3) [rotate=90] at (7,0) edges $\mathcal{E}_\mathrm{h}$;
[->] (s) edge (0);
[->] (s) edge (1);
[->] (s) edge (C);
[->] (0c) edge (t);
[->] (1c) edge (t);
[->] (Cc) edge (t);

(l1) [] at (3,3) Tile $1$;
(l2) [] at (5,3) Tile $2$;
(l3) [] at (7,3) Tile $\lvert \mathcal{W} \rvert$;

Transducer-based directed acyclic graph for generation of valid row tilings.

As follows from Section <ref>, any valid tiling of a row can be visualized as a $\lvert \mathcal{W} \rvert$-long path in the transducer graph $G_\mathrm{t,h}$, recall Section <ref>. To simplify subsequent developments, we represent the row-tiling problem by a transducer-based directed acyclic graph (DAG) composed of $\lvert \mathcal{W}\rvert +3$ vertex layers. While both the first and the last layer contain only a single vertex (the source $s$ and terminal $t$), the intermediate layers include $\lvert \mathcal{C} \rvert$ vertices to represent the vertical (east and west) color codes of the tiles, i.e., the states in the transducer graph. The source vertex is connected to all vertices in the second layer, facilitating an arbitrary coloring of the west edge of the first tile, and, similarly, all the vertices in the penultimate layer are linked to the terminal to allow for all colors in the last east edge. The intermediate layers are bridged with the transducer edges $\mathcal{E}_\mathrm{h}$; see Fig. <ref>. Consequently, any $s\rightarrow t$ path in the yet-established directed graph forms a valid tiling of the row, and conversely, any valid tiling builds a $s\rightarrow t$ path.

However, because such paths do not exist for tile sets that forbid a valid tiling of the row, we also need to incorporate voids. Clearly, we can add “void” tiles as edges that would interconnect the layers, i.e., any two consecutive layers would form a complete bipartite graph. However, such an approach requires adding at most $\lvert \mathcal{W} \rvert \lvert \mathcal{C} \rvert^2$ edges to the graph. Therefore, we add supplementary intermediate layers with a single vertex only, symbolizing the “void” tile type, and connect it to all vertices in the preceding and subsequent layer, see the dashed vertices and edges in Fig. <ref>. Consequently, we generate at most $2 \lvert \mathcal{W} \rvert \lvert \mathcal{C} \rvert$ new edges altogether.

[>=stealth',shorten >=1pt,auto,node distance=2cm, semithick,every node/.style=inner sep=2pt, every state/.style=minimum size=20pt,inner sep=2pt,on grid,node distance=4cm]

(s) [state] at (0,0) $s$;

(0) [state] at (2,2.25) $0$;
(1) [state] at (2,0.75) $1$;
(d) [] at (2,-0.75) $\vdots$;
(C) [state] at (2,-2.25) $\lvert\mathcal{C}\rvert$;

(0a) [state, right of=0]$0$;
(1a) [state, right of=1]$1$;
(da) [right of=d] $\vdots$;
(Ca) [state, right of=C]$\lvert\mathcal{C}\rvert$;

(0b) [right of=0a]$\dots$;
(1b) [right of=1a]$\dots$;
(db) [right of=da] $\ddots$;
(Cb) [right of=Ca]$\dots$;

(0c) [state, right of=0b]$0$;
(1c) [state, right of=1b]$1$;
(dc) [right of=db] $\vdots$;
(Cc) [state, right of=Cb]$\lvert\mathcal{C}\rvert$;

(t) [state] at (16,0)$t$;

(v0) [state, dashed] at (4,0)void;
(v1) [state, dashed] at (8,0)void;
(v2) [state, dashed] at (12,0)void;

(t1) [align=center] at (4,-2) + edges $\mathcal{E}_\mathrm{h}$;
(t2) [align=center] at (8,-2) + edges $\mathcal{E}_\mathrm{h}$;
(t3) [align=center] at (12,-2) + edges $\mathcal{E}_\mathrm{h}$;
[->] (s) edge [above,sloped] node 0.0 (0);
[->] (s) edge [above,sloped] node 0.0 (1);
[->] (s) edge [above,sloped] node 0.0 (C);
[->] (0c) edge [above,sloped] node 0.0 (t);
[->] (1c) edge [above,sloped] node 0.0 (t);
[->] (Cc) edge [above,sloped] node 0.0 (t);

[->, dashed] (0) edge [above,sloped] node 1.0 (v0);
[->, dashed] (1) edge [above,sloped] node 1.0 (v0);
[->, dashed] (C) edge [above,sloped] node 1.0 (v0);
[->, dashed] (v0) edge [above,sloped] node 0.0 (0a);
[->, dashed] (v0) edge [above,sloped] node 0.0 (1a);
[->, dashed] (v0) edge [above,sloped] node 0.0 (Ca);
[->, dashed] (0a) edge [above,sloped] node 1.0 (v1);
[->, dashed] (1a) edge [above,sloped] node 1.0 (v1);
[->, dashed] (Ca) edge [above,sloped] node 1.0 (v1);
[dashed] (v1) edge [above,sloped] node 0.0 ($(v1)!0.6!(0b)$);
[dashed] (v1) edge [above,sloped] node 0.0 ($(v1)!0.6!(1b)$);
[dashed] (v1) edge [above,sloped] node 0.0 ($(v1)!0.6!(Cb)$);
[->, dashed] ($(0b)!0.4!(v2)$) edge [above,sloped] node 1.0 (v2);
[->, dashed] ($(1b)!0.4!(v2)$) edge [above,sloped] node 1.0 (v2);
[->, dashed] ($(Cb)!0.4!(v2)$) edge [above,sloped] node 1.0 (v2);
[->, dashed] (v2) edge [above,sloped] node 0.0 (0c);
[->, dashed] (v2) edge [above,sloped] node 0.0 (1c);
[->, dashed] (v2) edge [above,sloped] node 0.0 (Cc);

(l1) [] at (4,3) Tile $1$;
(l2) [] at (8,3) Tile $2$;
(l3) [] at (12,3) Tile $\lvert \mathcal{W} \rvert$;

Transducer-based directed acyclic graph for computing the maximum row cover.

In addition, we assign unitary costs to the edges incoming to the void vertices and zero costs elsewhere. Hence, the $s\rightarrow t$ path cost is equivalent to the number of voids in the row tiling. Furthermore, because the emergent graph is acyclic and single-sourced, the maximum row-cover tiling is found in $\mathcal{O}(\lvert \mathcal{V} \rvert + \lvert \mathcal{E} \rvert)$ time using the DAG-shortest-path algorithm [55], where $\mathcal{V}$ denotes the set of the graph vertices and $\mathcal{E}$ the set of the graph edges. In our case, we have

\begin{align}
\lvert \mathcal{V} \rvert &= 2 + (\lvert \mathcal{W}\rvert + 1)\lvert \mathcal{C} \rvert + \lvert \mathcal{W}\rvert = 2 + \lvert \mathcal{W}\rvert + \lvert \mathcal{C}\rvert + \lvert \mathcal{W}\rvert\lvert \mathcal{C}\rvert,\\
\lvert \mathcal{E} \rvert &= 2\lvert \mathcal{C}\rvert + 2 \lvert \mathcal{W} \rvert \lvert \mathcal{C} \rvert + \lvert \mathcal{W} \rvert \lvert \mathcal{T} \rvert.
\end{align}

Thus, the overall asymptotic complexity to generate a maximum row cover tiling evaluates as

\begin{equation}
\mathcal{O}(\lvert \mathcal{V} \rvert + \lvert \mathcal{E} \rvert) = \mathcal{O}( 2 + \lvert \mathcal{W}\rvert + 3 \lvert \mathcal{C}\rvert + 3 \lvert \mathcal{W} \rvert \lvert \mathcal{C} \rvert + \lvert \mathcal{W} \rvert \lvert \mathcal{T} \rvert) = \mathcal{O}(\lvert \mathcal{W} \rvert \lvert \mathcal{C} \rvert + \lvert \mathcal{W} \rvert \lvert \mathcal{T} \rvert).
\end{equation}

Interestingly, the running time (but not the asymptotic complexity) of the DAG-shortest-path algorithm can be improved by recognizing that the topological order of the graph vertices—which is required for the DAG-shortest-path algorithm—is known from the graph construction method in advance.

Any path with total cost $c_\mathrm{t}$ contains exactly $c_\mathrm{t}$ voids in the row tiling. Because the shortest path algorithm therefore minimizes the number of voids, it generates the maximum row cover as its output. These considerations are summarized below.

The shortest path in the graph in Fig. <ref> is equivalent to the maximum row cover.

§.§ Tiling consecutive rows

Assuming already covered rows $i-1$ and $i+1$, e.g., initially by voids, we aim to generate the maximum cover of the $i$-th row. Interestingly, this only requires a minor modification of the graph in Fig. <ref>.

For this, we first check the north-east compatibility for each tile $k \in \mathcal{T}$ placed at $(i,j)$. Notice that the compatibility is never violated when the neighbors are voids. For color mismatch cases, we remove the edges denoting these incompatible tiles from the graph.

Assume that the rows $(i-1)$ and $(i+1)$ are voids. Then, clearly, inappropriate tiles at the $i$-th row may prevent the vertically-adjacent positions to be populated by tiles. To limit the appearance of such introduced voids, we include a small penalty of $\epsilon = 1/2 (\lvert \mathcal{W} \rvert+1)^{-1}$ to the tiles that admit a single vertical neighbor only, and $\epsilon = (\lvert \mathcal{W} \rvert+1)^{-1}$ to tiles not admitting any vertical neighbor. Notice that these costs are selected such that, in the worst case, the total penalty due to these void-preventing weights amounts to $\lvert\mathcal{W}\rvert/(\lvert \mathcal{W}\rvert+1)<1$, i.e., the maximum number of tiles is placed even if the void positions forbid any vertical neighbors. Hence, Proposition <ref> remains satisfied.

Consequently, we can build a simple heuristic algorithm,
Alg. <ref>, that requires $\lvert \mathcal{H} \rvert$ maximum row-cover iterations, rendering the overall complexity to be $\mathcal{O}(\lvert \tilde{\mathcal{A}} \rvert \lvert \mathcal{C} \rvert + \lvert \tilde{\mathcal{A}} \rvert \lvert \mathcal{T} \rvert)$.

Simple maximum cover heuristics
[1]
simpleMaximumCoverHeuristics$\mathcal{T}$, $\mathcal{A}$
$\mathfrak{T}$ $\gets$ ($\mathcal{A}$)
$G_{\mathrm{t},\mathrm{h}}$ $\gets$ ($\mathcal{T}$)
row $\gets$ $\{1,\dots,\lvert \mathcal{H} \rvert\}$
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $\mathcal{T}$, $\mathfrak{T}$, row)
shortestPath $\gets$ ($G_\mathrm{DAG}$)
$\mathfrak{T}$ $\gets$ ($\mathfrak{T}$, shortestPath, row)

Although Alg. <ref> usually generates relatively large ratio of the number of placed tiles $\lvert \tilde{\mathcal{B}}_\mathrm{cov}\rvert$ to $\lvert \tilde{\mathcal{A}}\rvert$, it probably lacks a guaranteed lower bound. Such bounds can, however, be provided by fairly straightforward modifications introduced next.

§.§ 1/2-approximation algorithm for general tile sets

In this section, we modify Alg. <ref> to maintain the $1/2$ approximation ratio. We start with the following observation:

Consider the maximum row-cover tiling of the odd rows of the initially void domain $\mathcal{A}$ given in Section <ref>. Then, $\lvert \tilde{\mathcal{B}}_\mathrm{cov}\rvert \ge 1/2 \lvert \tilde{\mathcal{B}}_{\max \mathrm{cov}} \rvert$.

Consider that the maximum row-cover problem alone terminates with $\lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}} \rvert$ tiles. Based on the maximum row-cover property in Proposition <ref>, none of the rows of $\tilde{\mathcal{A}}$ admit a tiling by more than $\lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}} \rvert$ tiles. Hence, we have $ \lvert\tilde{\mathcal{B}}_\mathrm{cov}\rvert \ge \lceil 1/2 \lvert\mathcal{H}\rvert \rceil \lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}}\rvert$ and $\lvert \tilde{\mathcal{B}}_{\max \mathrm{cov}}\rvert \le \lvert\mathcal{H}\rvert \lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}}\rvert$, so that $  \lvert \tilde{\mathcal{B}}_\mathrm{cov}\rvert \ge \lceil 1/2\lvert\mathcal{H}\rvert \rceil \lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}}\rvert \ge 1/2 \lvert\mathcal{H}\rvert \lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}}\rvert \ge 1/2 \lvert \tilde{\mathcal{B}}_{\max \mathrm{cov}} \rvert$, where $\lceil \bullet \rceil$ rounds $\bullet$ to the nearest greater or equal integer.

To exploit Proposition <ref> in Alg. <ref>, we modify the row processing order to $\{1,3,2,5,4\dots\}$. Indeed, then  each odd row contains exactly $\lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}}\rvert$ tiles. Nevertheless, covering the $i$-th (odd) row without acknowledging which tiles are placed in the $(i-2)$-th row may result in an unnecessarily empty $(i-1)$-th row. To avoid such situations, we do not check for compatibility with the $(i-1)$-th row voids, but rather we check using the dual transducer graph with the tiles in the $(i-2)$-th row. For each south color code in the $(i-2)$-th row, we find admissible colors (states) in the dual transducer graph as the states reachable by an edge-long path. Indeed, the reached states are exactly the admissible north colors of compatible tiles in the $i$-th row. For the special case of voids in the $(i-2)$-th row, all color codes are assumed to be compatible. Finally, we penalize the incompatibilities with the cost $\epsilon = 1/2 (\lvert \mathcal{W} \rvert+1)^{-1}$ as before. The final algorithm then reads as Alg. <ref>, allowing us to state the following, slightly stronger result:

$1/2$-approximation algorithm
[1]
maximumCoverApproximation050$\mathcal{T}$, $\mathcal{A}$
$\mathfrak{T}$ $\gets$ ($\mathcal{A}$)
$G_{\mathrm{t},\mathrm{h}}$, $G_{\mathrm{t},\mathrm{v}}$ $\gets$ ($\mathcal{T}$)
row $\gets$ $\{1,3,2,5,4,\dots\}$
row even
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $\mathcal{T}$, $\mathfrak{T}$, row)
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $G_{\mathrm{t},\mathrm{v}}$, $\mathfrak{T}$, row, $1$)
shortestPath $\gets$ ($G_\mathrm{DAG}$)
$\mathfrak{T}$ $\gets$ ($\mathfrak{T}$, shortestPath, row)

Assume a tile set $\mathcal{T}$ with the longest path in its transducer graph $G_{\mathrm{t},\mathrm{h}}$ of at least $2$. Then, Alg. <ref> terminates with $\lvert \tilde{\mathcal{B}}_\mathrm{cov}\rvert \ge 1/2 \lvert \tilde{\mathcal{A}}\rvert$.

When $\lvert\tilde{\mathcal{B}}_{\max \mathrm{rowcov}} \rvert = \lvert \mathcal{W}\rvert$, the proof follows directly from Proposition <ref>. For the other cases, the odd rows must contain $\lvert \tilde{\mathcal{B}}_{\max \mathrm{rowcov}} \rvert$ tiles due to Proposition <ref>. Because these row-covers are maximal, the sequence of consecutive voids in these rows cannot exceed two, as we could have placed an additional tile otherwise, contradicting with the maximum row-cover property. Moreover, without loss of generality, the cost of the shortest path in the $i$-th row is at most $\lvert\tilde{\mathcal{B}}_{\max \mathrm{rowcov}} \rvert + (\lvert \mathcal{W}\rvert - \lvert\tilde{\mathcal{B}}_{\max \mathrm{rowcov}} \rvert) \epsilon$, which occurs when the $(i-2)$-th and $i$-th row have the same tile-void patterns. Because the longest void sequence is at most two and the longest path in $G_{\mathrm{t},\mathrm{h}}$ is at least two, we can always place tiles to the north of the voids of the $i$-th row.

§.§ 2/3-approximation algorithm for tilesets with cyclic transducers

Another improvement in the approximation factor of Alg. <ref> is possible for tile sets with all the states in the transducer graphs $G_{\mathrm{t},\mathrm{h}}$ and $G_{\mathrm{t},\mathrm{v}}$ being in at least one graph cycle. Notice that this situation occurs for all tile sets that tile the infinite plane.

To this goal, we  modify the assignment of costs to graph, and the row processing order to $\{1,4,3,2,3,7,6,5,6,\dots\}$. We begin with (i) tiling the maximum row-cover of the first row. Then, we (ii) find the maximum row-cover of the $4$th row such that we penalize possible incompatibilities with the first row based on the dual transducer graph by $\epsilon$. The step (iii) encompasses finding a cover of the $3$rd row with penalized incompatibilities with the first row and enforced voids at even positions. Finally, we find the maximum covers for rows $2$ and $3$. We repeat the procedure for the row numbers iteratively increased by $3$, see Alg. <ref>.

$2/3$-approximation algorithm
[1]
maximumCoverApproximation067$\mathcal{T}$, $\mathcal{A}$
$\mathfrak{T}$ $\gets$ ($\mathcal{A}$)
$G_{\mathrm{t},\mathrm{h}}$, $G_{\mathrm{t},\mathrm{v}}$ $\gets$ ($\mathcal{T}$)
row $\gets$ $\{1,4,3,2,3,7,6,5,6\dots\}$
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $G_{\mathrm{t},\mathrm{v}}$, $\mathfrak{T}$, row, $2$)
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $G_{\mathrm{t},\mathrm{v}}$, $\mathfrak{T}$, row, $1$)
$G_\mathrm{DAG}$ $\gets$ ($G_\mathrm{DAG}$)
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $\mathcal{T}$, $\mathfrak{T}$, row)
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $\mathcal{T}$, $\mathfrak{T}$, row)
shortestPath $\gets$ ($G_\mathrm{DAG}$)
$\mathfrak{T}$ $\gets$ ($\mathfrak{T}$, shortestPath, row)

Then, we can make the following statement:

Consider that all states in the transducer graphs $G_{\mathrm{t},h}$ and $G_{\mathrm{t},\mathrm{v}}$ are in at least one graph cycle. Then, Alg. <ref> terminates with at least $\frac{2}{3} \lvert \tilde{\mathcal{A}}\rvert$ placed tiles.

Since the tile set allows for valid tiling of the row, the $\{1,4,\;\dots\}$ rows are occupied by exactly $\lvert \mathcal{W} \rvert$ tiles. The $\{3,6,\;\dots\}$ rows are then populated by at least $1/2\lvert \mathcal{W}\rvert$ tiles because each tile from rows $\{4,7,\dots\}$ admits a vertical neighbor. Finally, the $\{2,5,\;\dots\}$ rows contain at least the complement of the number of tiles used in the preceding row, because the tiles in the $\{1,4,6,\dots\}$ row admit a south neighbor. Depending on the number of rows, the algorithm places at least

\begin{equation}
\lvert \tilde{\mathcal{B}}_{\mathrm{cov}} \rvert \ge \min\{\lvert \tilde{\mathcal{A}} \rvert, \frac{3}{4}\lvert \tilde{\mathcal{A}} \rvert, \frac{2}{3}\lvert \tilde{\mathcal{A}} \rvert, \frac{3}{4}\lvert \tilde{\mathcal{A}} \rvert,\frac{7}{10}\lvert \tilde{\mathcal{A}} \rvert,\frac{2}{3}\lvert \tilde{\mathcal{A}} \rvert,\dots \} = \frac{2}{3} \lvert \tilde{\mathcal{A}} \rvert
\end{equation}

§.§ Iterative improvements

Similarly to finding the maximum row covers, we can search for the maximum cover of columns. When combining these two methods, we end up with our final algorithm that has the $\mathcal{O}(\lvert \tilde{\mathcal{A}} \rvert^2 \lvert \mathcal{C} \rvert + \lvert \tilde{\mathcal{A}} \rvert^2 \lvert \mathcal{T} \rvert + \lvert \mathcal{C} \rvert^2)$ complexity and provides the approximation ratios adjustable by algorithm choice (Algs. <ref>, <ref> or <ref>) at line <ref> of Alg. <ref>.

Final maximum cover heuristics
[1]
finalMaximumCoverHeuristics$\mathcal{T}$, $\mathcal{A}$
$\mathfrak{T}$ $\gets$ ($\mathcal{T}$,$\mathcal{A}$)
$G_{\mathrm{t},\mathrm{h}}$, $G_{\mathrm{t},\mathrm{v}}$ $\gets$ ($\mathcal{T}$)
numVoidsOld $\gets$ $\infty$
method $\gets$ “columns”
numVoidsOld$-$($\mathfrak{T}$) $> 0$
numVoidsOld $\gets$ ($\mathfrak{T}$)
row $\gets\{1,\dots,\lvert \mathcal{H}\rvert\}$
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{h}}$, $\mathfrak{T}$, row)
shortestPath $\gets$ ($G_\mathrm{DAG}$)
$\mathfrak{T}$ $\gets$ ($\mathfrak{T}$, shortestPath, row)
method $\gets$ “columns”
column $\gets\{1,\dots,\lvert \mathcal{W}\rvert\}$
$G_\mathrm{DAG}$ $\gets$ ($G_{\mathrm{t},\mathrm{v}}$, $\mathfrak{T}$, column)
shortestPath $\gets$ ($G_\mathrm{DAG}$)
$\mathfrak{T}$ $\gets$ ($\mathfrak{T}$, shortestPath, column)
method $\gets$ “rows”

Alg. <ref> runs in a polynomial time and terminates in a finite number of steps.

We have already shown that finding a maximum row-cover has $\mathcal{O}(\lvert \mathcal{W} \rvert \lvert \mathcal{C} \rvert + \lvert \mathcal{W} \rvert \lvert \mathcal{T} \rvert)$ complexity. Further, finding the $2$-long paths in the transducer graph possesses the $\lvert \mathcal{C}\rvert^2$ complexity and can be run only once prior to the algorithm main loop. Altogether, Alg. <ref> requires at most $4/3\lvert \mathcal{H}\rvert$ inner iterations so that we have the $\mathcal{O}(\lvert \tilde{\mathcal{A}} \rvert \lvert \mathcal{C} \rvert + \lvert \tilde{\mathcal{A}} \rvert \lvert \mathcal{T} \rvert +  \lvert \mathcal{C} \rvert^2)$ overall complexity.

Regardless of the method at line <ref> of Alg. <ref>, the improving loop runs at most $\lvert \tilde{\mathcal{A}} \rvert$ times. Consequently, the algorithm is finite and possesses the $\mathcal{O}(\lvert \tilde{\mathcal{A}} \rvert^2 \lvert \mathcal{C} \rvert + \lvert \tilde{\mathcal{A}} \rvert^2 \lvert \mathcal{T} \rvert + \lvert\mathcal{C}\rvert^2)$ complexity.

§ RESULTS

Having developed several exact and heuristic methods, this section is devoted to their numerical examination. We begin with assessing the performance of the integer programming formulations in Section <ref>. Then, in Section <ref>, we also relate these results to the outputs of the heuristic algorithms.

Extensions of the integer programs are investigated in subsequent sections. First, we demonstrate the usefulness of the packing constraint by comparing the efficiency of the solution to the tile-packing problem using our method with the times reported by Lagae and Dutré [14], Section <ref>. Subsequently, we also present two unexpected discoveries revealed when testing formulations: the Knuth [21] tile set contains a tile unusable in infinite tilings, Section <ref>, and the Lagae et al. [33] tile set of $44$ corner tiles lacks aperiodicity, Section <ref>.

We implemented all the methods described above in C++. As the integer programming solver, we used the state-of-the-art optimizer Gurobi 9.5.0 [56] dynamically linked to the compiled binary. Numerical tests were evaluated on a personal laptop running the Ubuntu $18.04$ operating system equipped with $24$ GB of RAM and Intel$^\text{\textregistered}$ Core$^\text{\textregistered}$ i5-8350U CPU clocked at 1.70GHz.

§.§ Integer programming formulations

In this section, we investigate the performance of all integer programming formulations from Section <ref>, i.e., the decision program (<ref>), the maximum rectangular tiling (<ref>), the maximum cover (<ref>), and the maximum adjacency constraint satisfaction problem (<ref>).

New tile sets (a) Finite1 of $7$ tiles over $4$ colors, and (b) Finite2 of $16$ tiles over $16$ colors used in our algorithmic tests.

26emTile set     22emSize     2c|Dec. prog. (<ref>)     2c|Max. rect. (<ref>)     2c|Max. cov. (<ref>)     2cMax. CSP (<ref>)

Time [s]     Objective     Time [s]     Objective     Time [s]     Objective     Time [s]     Objective

36emAperiodic1 ($11/4$) [31]     $20\times20$      $0.111$     $0$     $129.897$     $400$     $300.053$     $^*398$     $300.056$     $^*745$

$25\times25$     $90.810$     $0$     $300.070$     $^*150$     $300.083$     $^*606$     $300.071$     $^*1135$

$30\times30$     $300.069$     $^*$infeasible     $300.084$     $^*150$     $300.097$     $^*861$     $300.089$     $^*1628$

36emAperiodic2 ($13/5$) [28]     $20\times20$    $0.114$     $0$     $145.300$     $400$     $300.055$     $^*399$     $300.059$     $^*742$

$25\times25$     $178.337$     $0$     $300.070$     $^*125$     $300.082$     $^*612$     $300.078$     $^*1184$

$30\times30$     $300.069$     $^*$infeasible     $300.089$     $^*60$     $300.098$     $^*876$     $300.111$     $^*1655$

36emAperiodic3 ($14/6$) [27]     $20\times20$     $275.339$     $0$     $181.171$     $400$     $300.058$     $^*397$     $300.058$     $^*752$

$25\times25$     $300.057$     $^*$infeasible     $300.072$     $^*100$     $300.086$     $^*619$     $300.086$     $^*1178$

$30\times30$     $300.073$     $^*$infeasible     $300.092$     $^*90$     $300.107$     $^*863$     $300.104$     $^*1610$

36emAperiodic4 ($16/6$) [24]     $20\times20$     $0.142$     $0$     $171.136$     $400$     $176.584$     $400$     $71.141$     $760$

$25\times25$     $0.196$     $0$     $300.063$     $^*100$     $300.251$     $^*577$     $300.085$     $^*1030$

$30\times30$     $0.251$     $0$     $300.265$     $^*60$     $300.132$     $^*794$     $300.115$     $^*1616$

36emAperiodic5 ($56/12$) [23]     $20\times20$     $0.294$     $0$     $300.107$     $^*20$     $300.214$     $^*350$     $302.442$     $^*688$

$25\times25$     $0.440$     $0$     $300.155$     $^*25$     $300.354$     $^*553$     $300.197$     $^*1055$

$30\times30$     $0.648$     $0$     $300.228$     $^*30$     $300.434$     $^*795$     $301.102$     $^*1569$

36emStochastic1 ($8/2$) [9]     $20\times20$     $0.066$     $0$     $0.101$     $400$     $0.046$     $400$     $4.195$     $760$

$25\times25$     $0.091$     $0$     $0.125$     $625$     $0.089$     $625$     $5.598$     $1200$

$30\times30$     $0.116$     $0$     $0.225$     $900$     $0.110$     $900$     $10.021$     $1740$

36emStochastic2 ($16/4$) [32]     $20\times20$     $0.114$     $0$     $0.107$     $400$     $0.129$     $400$     $3.226$     $760$

$25\times25$     $0.141$     $0$     $0.175$     $625$     $0.210$     $625$     $6.118$     $1200$

$30\times30$     $0.183$     $0$     $0.217$     $900$     $0.283$     $900$     $6.846$     $1740$

36emPeriodic1 ($10/4$) [18]     $20\times20$     $0.121$     $0$     $107.475$     $400$     $111.982$     $400$     $54.696$     $760$

$25\times25$     $0.153$     $0$     $274.813$     $625$     $300.066$     $^*584$     $224.734$     $1200$

$30\times30$     $0.193$     $0$     $300.977$     $^*81$     $302.606$     $^*824$     $300.087$     $^*1628$

36emPeriodic2 ($30/17$) [34]     $20\times20$     $0.236$     $0$     $109.860$     $400$     $252.700$     $400$     $88.721$     $760$

$25\times25$     $0.325$     $0$     $300.103$     $^*25$     $300.222$     $^*545$     $300.150$     $^*1017$

$30\times30$     $0.473$     $0$     $300.158$     $^*30$     $300.300$     $^*786$     $300.204$     $^*1521$

($7/4$)     $20\times20$     $0.066$     infeasible     $300.025$     $^*120$     $300.051$     $^*378$     $300.076$     $^*725$

$25\times25$     $0.086$     infeasible     $300.038$     $^*125$     $300.054$     $^*585$     $300.061$     $^*1108$

$30\times30$     $0.105$     infeasible     $300.046$     $^*108$     $300.069$     $^*826$     $300.080$     $^*1628$

36emFinite2 ($16/16$)     $20\times20$     $0.100$     infeasible     $300.273$     $^*40$     $300.207$     $^*326$     $300.077$     $^*684$

$25\times25$     $0.133$     infeasible     $300.067$     $^*50$     $300.111$     $^*493$     $300.094$     $^*1029$

$30\times30$     $0.168$     infeasible     $300.084$     $^*30$     $300.131$     $^*690$     $300.128$     $^*1525$

Benchmark results. Values marked by an asterisk denote a premature termination of the integer programming solver.

We are unaware of any standard sets for bounded tiling problems except for the specific, mostly aperiodic tile sets listed in the literature, recall Section <ref>. Hence, we consider a set of benchmark problems consisting of five aperiodic tile sets ($11$ tiles over $4$ colors by Jeandel and Rao [31], $13$ tiles over $5$ colors  by Čulík [28], $14$ tiles over $6$ colors by Kari [27], $16$ tiles over $6$ colors by Ammann [24], and $56$ tiles over $12$ colors by Robinson [23]), two stochastic tile sets introduced in computer graphics ($8$ tiles over $2$ colors by Cohen et al. [9] and a set of $16$ tiles over $4$ edge colors by Lagae and Dutré  [32]), two periodic tile sets ($10$ tiles over $4$ colors by Wang [18] and the set of $30$ tiles over $17$ edge colors by Lagae et al. [33] and Nurmi [34]). In addition, in Fig. <ref>, we introduce two tile sets that do not allow for a valid tiling of the infinite domain.

For all these tile sets, we aimed at generating valid tilings sized, respectively, $20\times 20$, $25 \times 25$, and $30 \times 30$. The running time of the Gurobi solver was limited to $300$ seconds for the single-threaded mode.

The results shown in Tab. <ref> illustrate that the performance of the decision program (<ref>) surpasses any of the candidate variants. However, it failed to find an existent feasible solution in the time limit four times. In these cases, the output of the optimization problems (<ref>, <ref>, <ref>) provided at least some output. Interestingly, the decision problem (<ref>) also was more efficient in the case of proving that the domain $\lvert \mathcal{A}\rvert$ lacks $\mathcal{T}$-tilability.

Comparison of the optimization variants hints that the maximum cover (<ref>) and the maximum adjacency constraint satisfaction (<ref>) problems scale better than the maximum rectangular tiling (<ref>). Indeed, generating any smaller rectangular domain remains -complete, preventing any polynomial-time approximation algorithm to exist. On the other hand, both the formulations (<ref>) and (<ref>) admit simple heuristics, recall Section <ref>, allowing the solver to obtain higher-quality feasible solutions faster.

§.§ Heuristic algorithms

24emTile set     22emSize     4c|Alg. <ref> with Alg. <ref>     4c|Alg. <ref> with Alg. <ref>     4cAlg. <ref> with Alg. <ref>

$t$ [s]     min     avg     max     $t$ [s]     min     avg     max     $t$ [s]     min     avg     max

36emAperiodic1 ($11/4$) [31]     $20\times20$     $0.024$     $358$     $\mathbf{368.99}$     $380$     $0.056$     $342$     $360.22$     $372$     $0.029$     $334$     $347.60$     $360$

$25\times25$     $0.040$     $562$     $\mathbf{575.38}$     $586$     $0.101$     $543$     $563.67$     $582$     $0.046$     $524$     $537.68$     $547$

$30\times30$     $0.065$     $813$     $\mathbf{829.66}$     $843$     $0.160$     $792$     $812.33$     $835$     $0.080$     $758$     $779.91$     $798$

36emAperiodic2 ($13/5$) [28]     $20\times20$     $0.043$     $354$     $\mathbf{369.86}$     $381$     $0.054$     $326$     $359.86$     $373$     $0.051$     $326$     $353.21$     $373$

$25\times25$     $0.065$     $564$     $\mathbf{577.31}$     $590$     $0.099$     $504$     $564.92$     $581$     $0.118$     $520$     $557.49$     $586$

$30\times30$     $0.123$     $818$     $\mathbf{831.60}$     $847$     $0.179$     $800$     $817.11$     $837$     $0.186$     $760$     $806.65$     $838$

36emAperiodic3 ($14/6$) [27]     $20\times20$     $0.034$     $362$     $375.40$     $386$     $0.032$     $353$     $365.57$     $381$     $0.042$     $355$     $\mathbf{378.87}$     $388$

$25\times25$     $0.058$     $564$     $585.56$     $604$     $0.056$     $548$     $569.15$     $597$     $0.065$     $562$     $\mathbf{592.28}$     $604$

$30\times30$     $0.092$     $813$     $843.19$     $857$     $0.095$     $799$     $830.93$     $860$     $0.102$     $827$     $\mathbf{855.68}$     $871$

36emAperiodic4 ($16/6$) [24]     $20\times20$     $0.031$     $351$     $\mathbf{366.09}$     $381$     $0.077$     $293$     $333.98$     $351$     $0.092$     $281$     $339.30$     $368$

$25\times25$     $0.052$     $555$     $\mathbf{573.19}$     $591$     $0.149$     $469$     $524.82$     $549$     $0.187$     $442$     $533.27$     $562$

$30\times30$     $0.096$     $795$     $\mathbf{825.44}$     $860$     $0.234$     $666$     $758.52$     $785$     $0.276$     $743$     $773.66$     $802$

36emAperiodic5 ($56/12$) [23]     $20\times20$     $0.054$     $344$     $\mathbf{360.55}$     $381$     $0.149$     $256$     $341.48$     $364$     $0.171$     $290$     $332.83$     $349$

$25\times25$     $0.110$     $540$     $\mathbf{563.41}$     $607$     $0.289$     $402$     $527.68$     $563$     $0.472$     $484$     $529.63$     $552$

$30\times30$     $0.147$     $782$     $\mathbf{811.17}$     $856$     $0.473$     $598$     $782.99$     $809$     $0.601$     $706$     $759.48$     $786$

36emStochastic1 ($8/2$) [9]     $20\times20$     $0.014$     $400$     $\mathbf{400.00}$     $400$     $0.012$     $400$     $\mathbf{400.00}$     $400$     $0.014$     $400$     $\mathbf{400.00}$     $400$

$25\times25$     $0.013$     $625$     $\mathbf{625.00}$     $625$     $0.014$     $625$     $\mathbf{625.00}$     $625$     $0.021$     $625$     $\mathbf{625.00}$     $625$

$30\times30$     $0.016$     $900$     $\mathbf{900.00}$     $900$     $0.016$     $900$     $\mathbf{900.00}$     $900$     $0.019$     $900$     $\mathbf{900.00}$     $900$

36emStochastic2 ($16/4$) [32]     $20\times20$     $0.013$     $400$     $\mathbf{400.00}$     $400$     $0.015$     $400$     $\mathbf{400.00}$     $400$     $0.015$     $400$     $\mathbf{400.00}$     $400$

$25\times25$     $0.017$     $625$     $\mathbf{625.00}$     $625$     $0.017$     $625$     $\mathbf{625.00}$     $625$     $0.019$     $625$     $\mathbf{625.00}$     $625$

$30\times30$     $0.025$     $900$     $\mathbf{900.00}$     $900$     $0.022$     $900$     $\mathbf{900.00}$     $900$     $0.026$     $900$     $\mathbf{900.00}$     $900$

36emPeriodic1 ($10/4$) [18]     $20\times20$     $0.026$     $342$     $\mathbf{354.82}$     $374$     $0.059$     $325$     $339.42$     $353$     $0.044$     $323$     $337.66$     $346$

$25\times25$     $0.043$     $533$     $\mathbf{553.13}$     $573$     $0.103$     $507$     $528.02$     $546$     $0.065$     $512$     $527.11$     $537$

$30\times30$     $0.065$     $770$     $\mathbf{797.26}$     $830$     $0.165$     $745$     $764.69$     $784$     $0.097$     $729$     $758.36$     $770$

36emPeriodic2 ($30/17$) [34]     $20\times20$     $0.092$     $366$     $\mathbf{382.23}$     $400$     $0.419$     $259$     $337.18$     $368$     $0.346$     $281$     $345.72$     $369$

$25\times25$     $0.153$     $559$     $\mathbf{595.04}$     $620$     $0.808$     $503$     $533.30$     $568$     $0.649$     $537$     $563.63$     $575$

$30\times30$     $0.237$     $825$     $\mathbf{854.86}$     $887$     $1.259$     $724$     $760.08$     $802$     $1.131$     $768$     $800.74$     $822$

($7/4$)     $20\times20$     $0.023$     $353$     $360.71$     $369$     $0.047$     $352$     $\mathbf{362.98}$     $376$     $0.054$     $348$     $357.72$     $366$

$25\times25$     $0.039$     $548$     $\mathbf{561.15}$     $570$     $0.073$     $348$     $559.37$     $580$     $0.091$     $539$     $554.12$     $566$

$30\times30$     $0.056$     $789$     $807.79$     $819$     $0.111$     $789$     $\mathbf{809.48}$     $830$     $0.128$     $791$     $804.21$     $826$

36emFinite2 ($16/16$)     $20\times20$     $0.065$     $334$     $\mathbf{344.36}$     $363$     $0.113$     $286$     $331.80$     $353$     $0.104$     $291$     $328.83$     $355$

$25\times25$     $0.101$     $518$     $536.20$     $552$     $0.198$     $493$     $529.24$     $545$     $0.140$     $511$     $\mathbf{544.78}$     $559$

$30\times30$     $0.191$     $745$     $\mathbf{770.34}$     $789$     $0.393$     $695$     $761.45$     $780$     $0.382$     $678$     $750.89$     $790$

Numerical tests of the maximum-cover heuristics, Alg. <ref>, initialized based on Algs. <ref>, <ref>, and <ref>. Best mean runs are highlighted in bold.

Second, we compare the performance of the maximum cover formulation (<ref>) solved with the heuristic Alg. <ref> supplied with three different initial coverings, i.e., based on Algs. <ref>, <ref> and <ref>.

Alg. <ref> ran sequentially. In order to limit the dependence of the heuristic algorithm on the ordering of tiles, we randomized the edge order in the directed acyclic graphs. Thus, we evaluated Alg. <ref> $100$ times for each of the tested option, and listed the best, worst, and mean results in Tab. <ref>.

From Tab. <ref>, it follows that the initialization with the cover from Alg. <ref> is the most efficient for the tested tile sets, both in terms of speed and performance. The remaining two initializations seem to be fairly comparable on average. While for Alg. <ref>, at least $82\%$ of tiles were always placed, only more than $60\%$ followed from Alg. <ref>. Using Alg. <ref>, we obtained at least $70\%$ tile placement.

When comparing Tab. <ref> with Tab. <ref>, a few patterns emerge. First, the heuristic algorithm always generates valid tilings if (any of) the stochastic tile sets are used. For aperiodic and periodic tile sets, Gurobi required a considerably longer time to reach feasible solutions of a similar quality, but usually surpassed the developed algorithms in the time limit of $300$ s. In the case of Alg. <ref>, it can be seen that the resulting covers are very competitive to the outputs of (<ref>) and also obtained in much shorter times.

Periodic packing of a complete set of $625$ tiles over $5$ colors.

§.§ Periodic tile packing problem

Tile set     Time     Time

(<ref>,<ref>,<ref>)     Lagae and Dutré [14]

Stochastic edge ($16/2$)     $<1$ sec.     $<1$ sec.

Stochastic edge ($81/3$)     $<1$ sec.     $<1$ sec.

Stochastic edge ($256/4$)     $9$ sec.     $140$ days

Stochastic edge ($625/5$)     $4$ days     -

Periodic tile packing problem: comparison of core times needed to find a single feasible solution by integer programming (second column) and by the backtracking method (third column) proposed in Lagae and Dutré [14] to find a feasible solution.

As the second numerical example, we consider the periodic tile packing problem investigated in computer graphics applications [14]. Considering a complete edge tile set, Lagae and Dutré searched for a periodic square valid tiling with each tile from the tile set used exactly once. Clearly, such tilings not only contain the entire (textural) information stored in individual tiles but also maintain compatibility with the traditional periodic arrangement.

While Lagae and Dutré [14] proposed a backtracking-based algorithm to generate periodic packings, we rely here on a solution to the decision program (<ref>) supplemented with the packing (<ref>) and fixed periodicity (<ref>) constraints. The resulting core times spent in the search for a single feasible solution (Tab. <ref>) illustrates the higher effectiveness of our method.

§.§ Unusable tile in the Knuth tile set

One of the oldest aperiodic tile sets, containing $92$ tiles over $26$ colors, is from Knuth <cit.>. Generating valid tilings from the Knuth tile set using the decision program (<ref>) together with the tile-based boundary conditions, recall Section <ref>, led to an unexpected observation that enforced placement of the $\beta U S$ tile makes the program (<ref>) infeasible under certain circumstances.

After a careful investigation, it indeed turned out that there is not any $2 \times 2$ valid tiling with the $\beta U S$ tile placed at $(2,2)$. Moreover, there is also not any $4 \times 3$ valid tiling with the $\beta U S$ tile placed at $(3,1)$. Thus, using the maximum-cover optimization variant (<ref>) and the $\beta U S$ tile enforced at the respective coordinate, there are exactly $31$ optimal solutions with the objective function equal to $3$, and $498$ optimal solutions with the objective function equal to $11$.

Consequently, the $\beta U S$ tile can appear only in the strip of at most $3$ consecutive infinite columns and does not allow for simply-connected valid tilings of the infinite plane. In a private communication, Knuth confirmed the issue and discovered another $5$ tiles that are unnecessary but usable in infinite valid tilings, allowing for a possible reduction of the tile set to $86$ tiles. For more information, we refer the interested reader to Knuth's discussion about the reduced tile set <cit.>.

§.§ Periodicity of the Lagae corner tile set

Analogously to the Wang tiles, with the connectivity information stored in the edges, Lagae and Dutré [32] introduced corner tiles with  colored corners. As Wang [18] noted in 1975, these formalisms are interchangeable if the (infinite) domino problem is considered, because every set of Wang tiles can be represented by sets of corner tiles with greater or equal cardinality [33]. However, corner tiles avoid the so-called corner problem of Wang tiles in computer graphics [32], motivating Lagae et al. [33] to develop conversion methods for transforming Wang tiles to corner tiles, and vice versa. A direct product of these conversions are aperiodic tile sets of corner tiles [33].

Two of these methods, called horizontal and vertical translations, were used to convert the Ammann set of $16$ Wang tiles over $6$ colors [24] to the set of $44$ corner tiles over $6$ colors, and the resulting isomorphic corner tile sets were claimed aperiodic [33]. In 2016, Nurmi [34] noticed that, in this set, $14$ tiles are unusable in infinite valid tilings, and reduced the corner tile set to $30$ tiles over $6$ colors. Quite surprisingly, neither Lagae et. al. nor Nurmi recognized that the tile set forms a torus, and is therefore periodic, as we show next.

Instead of developing a new formulation for another type of tiles, we first notice that corner tiles are actually a subset of Wang tiles, and therefore every set of corner tiles can be represented by a set of Wang tiles with the same cardinality, see Appendix <ref>. For these tiles, we solve the rectangular tiling formulation (<ref>) with periodic boundary conditions (<ref>) and an objective function to find the smallest tiling (<ref>). As its output, we receive the optimal value of $6$ and $12$ optimal periodic rectangular tilings of the size $2 \times 3$. Not surprisingly, all these solutions follow from only two periodic patterns shown in Fig. <ref> by translations over the infinite plane.

Rectangular periodic valid tilings. Translating a $2\times3$ rectangle over the infinite valid tiling generated from (a) or (b) leads to $6$ different periodic patterns of the same size. Consequently, the tile set allows for $12$ periodic rectangles of the size $2\times 3$.

Having revealed the smallest periodic solutions, it remains to be shown why the Lagae conversion methods failed. Lagae et al. [33] mentioned that their methods lack bijectiveness in general but they assumed it was not the case here. Therefore, we believe it is useful to state the conditions under which the methods are bijective and show that they are not satisfied for the Ammann tile set.

The horizontal translation method from Lagae et al. [33] is bijective iff the dual transducer graph $G_\mathrm{T,v}$ of the input tile set $\mathcal{T}$ does not contain any parallel arcs.

The horizontal translation method is formally a mapping $\mathcal{T} \times \mathcal{T} \mapsto \mathcal{T}_\mathrm{corner}$ that generates $\forall (p,q) \in \mathcal{T}^2: c_p^\mathrm{e} = c_q^\mathrm{w}$ a corner tile $(c_p^\mathrm{n}, c_p^\mathrm{s}, c_q^\mathrm{s}, c_q^\mathrm{n})$. To be bijective, the cardinality of the output needs to be equal to the cardinality of the input, and the mapping has to produce unique output for each input. Consequently, all the tiles $p \in \mathcal{T}$ in the original tile set must be uniquely determined by $c_p^\mathrm{n}$ and $c_p^\mathrm{s}$, as the color codes of the vertical edges of $\mathcal{T}$ are avoided in the construction of $\mathcal{T}_\mathrm{corner}$.

Let us now consider that the dual transducer graph contains a parallel arc connecting the state $c^\mathrm{n}$ with $c^\mathrm{s}$. Then, there may exist two tiles colored by $(c^\mathrm{n}, c_p^\mathrm{w}, c^\mathrm{s}, c_p^\mathrm{e})$ and $(c^\mathrm{n}, c_q^\mathrm{w}, c^\mathrm{s}, c_q^\mathrm{e})$ that are indistinguishable in $\mathcal{T}_\mathrm{corner}$, which contradicts the bijection. For the other option, if the transducer graph does not contain any parallel arcs, then each $c_q^\mathrm{n}$, $c_q^\mathrm{s}$ identifies with a single arc labeled by $c_q^\mathrm{w} \vert c_q^\mathrm{e}$, i.e., with a single tile, which completes the proof.

Rotating the tile set by $90$ degrees, the arguments in Lemma <ref> provide us with the conditions for the bijectiveness of the vertical translation method:

The vertical translation method of Lagae et al. [33] is bijective iff the transducer graph $G_\mathrm{T,h}$ of the input tile set $\mathcal{T}$ does not contain any parallel arcs.

0.55![>=stealth',shorten >=1pt,auto,node distance=2cm, semithick,every node/.style=inner sep=2pt, every state/.style=minimum size=0pt,inner sep=2pt]
(0) [state] $0$;
(1) [state, below=of 0] $1$;
[->] (0) edge [bend left=40, above, sloped, rotate=180] node(b) $1\vert0$ (1);
[->] (0)  edge [loop above, sloped] node $3\vert4$ (0);
[->] (1)  edge [loop below, sloped] node $5\vert2$ (1);
[->] (1) edge [bend left=40, above, sloped]  node $2\vert4$ (0);
[->] (1) edge node [above, sloped] $5\vert3$ (0);

(2) [state, right=of b.north] $2$;
(5) [state, below right=of 2] $5$;
(4) [state, above right=of 2] $4$;
(3) [state, below right=of 4] $3$;
[->] (2) edge node [above, sloped] $3\vert2$ (3);
[->] (2) edge [loop left] node $3\vert5$ (2);
[->] (2) edge [bend left, sloped, above] node $3\vert3$ (4);
[->] (2) edge [bend left, sloped, above] node $1\vert1$ (5);
[->] (5) edge [bend left, sloped, above] node $2\vert5$ (2);
[->] (3) edge [bend left, sloped, above] node $4\vert3$ (4);
[->] (3) edge [bend left, sloped, above] node $0\vert1$ (5);
[->] (5) edge [bend left, sloped, above] node $2\vert2$ (3);
[->] (3) edge [loop right] node $4\vert2$ (3);
[->] (4) edge [bend left, sloped, above] node $0\vert1$ (2);
[->] (4) edge [bend left, sloped, above] node $0\vert0$ (3);

Transducer graph of the Ammann set of $16$ Wang tiles over $6$ colors.

For the Ammann tile set, we obtain the transducer graph $G_\mathrm{T,h} = G_\mathrm{T,v}$ shown in Fig. <ref>. Clearly, there exist parallel arcs $1 \rightarrow 0$. Moreover, using the same approach, we can show that the horizontal translation method also fails for the Robinson tile set of $24$ tiles over $24$ colors [24], contrary to the claims in [33], and the corresponding corner tile set is also periodic.

§ CONCLUSIONS

In this contribution, we investigated methods generating bounded Wang tilings for arbitrary tile sets. To this goal, we developed an -complete binary linear programming formulation (<ref>), as well as its -hard optimization variants relaxing some of the initial assumptions: tilability of the entire rectangular domain leading to the maximum rectangular tiling formulation (<ref>), simple-connectedness to the maximum-cover program (<ref>), and tiling validity to the maximum constraint satisfaction problem (<ref>). In addition, we supplemented these formulations with extensions enabling control over individual tiles and their colors, including a tile-packing constraint to enable generation of periodic tile packings, and introduced the variable-sized periodic constraints to facilitate computation of the smallest periodic patterns.

Motivated by the -hardness of all optimization formulations, we developed simple yet efficient heuristic algorithms for the maximum-cover variant (<ref>). These algorithms rely on the fact that generating the maximum row cover is equivalent to finding the shortest path in a directed acyclic graph. Moreover, well-chosen costs of the graph edges also maintain color matches with the neighboring rows. Thus, in the simplest case, a heuristic solution follows from a sequential generation of row-cover tilings. Moreover, with simple modifications to the row processing order, we showed how to provide a $1/2$ approximation factor for general tile sets and a $2/3$ guarantee for tile sets whose transducer graphs are cyclic.

We illustrated the effectiveness of these methods on a collection of $11$ tile sets. Generating tilings sized $20\times20$, $25\times25$, and $30\times30$ revealed that the decision program (<ref>) is the most efficient for our test problems. However, when a time limit is imposed or if the tile set does not allow for valid tiling of the entire domain, then the maximum cover (<ref>) and maximum adjacency constraint satisfaction problems (<ref>) appear to be similarly efficient. The remaining formulation—maximum rectangular tiling (<ref>)—exhibits the worst scalability.

The usefulness of integer programming extensions was demonstrated by means of three illustrative problems: showing a better solution efficiency to the tile packing problem than the Lagae and Dutré [14] backtracking approach, revealing an unusable tile in the Knuth [21] tile set, and proving that the Lagae et al. [33] tile set of corner tiles lacks aperiodicity. For the latter case, we also included an explanation for why the tile set construction method failed.

Furthermore, we tested the performance of the heuristic algorithms against the Gurobi solver tackling integer problem (<ref>) for $300$ s. Testing revealed that the simplest setup of Alg. <ref> initialized with the cover generated by Alg. <ref> produces the best results, obtained faster and (in most cases) competitive with the outputs of Gurobi. Somewhat surprisingly, the variants with guaranteed lower bounds exhibited slightly worse performance on average.

Having summarized our contributions, we believe that this work has not only introduced new methods that can possibly be applied to materials engineering, but also a simple and quite extensible framework to verify theoretical results on Wang tilings.

§ CORNER TILES REPRESENTED AS WANG TILES

Each corner tile in the corner tile set is defined by a quadruple of color codes $(c_k^\mathrm{nw},c_k^\mathrm{sw},c_k^\mathrm{se},c_k^\mathrm{ne})$, with $c_k^\mathrm{nw}$, $c_k^\mathrm{sw}$, $c_k^\mathrm{se}$, and $c_k^\mathrm{ne}$ denoting the colors of the northwest, southwest, southeast, and northeast corner of the $k$-th tile. Similarly to Wang tiles, corner tiles are assembled such that the color codes at the adjoining corners match.

This is, however, also maintained if we denote their edges by labels in the form of tuples of the corner codes, $(c_k^\mathrm{nw}, c_k^\mathrm{ne})$, $(c_k^\mathrm{nw}, c_k^\mathrm{sw})$, $(c_k^\mathrm{sw}, c_k^\mathrm{se})$, and $(c_k^\mathrm{ne}, c_k^\mathrm{se})$, each of which denotes a single edge label of the north, west, south, and east edge, respectively. Consequently, we can compute unique color codes as

\begin{equation}
c_k^\mathrm{n} = c_k^\mathrm{nw} + c_k^\mathrm{ne} n_\mathrm{vc},
\end{equation}

\begin{equation}
c_k^\mathrm{w} = c_k^\mathrm{nw} + c_k^\mathrm{sw} n_\mathrm{vc},
\end{equation}

\begin{equation}
c_k^\mathrm{s} = c_k^\mathrm{sw} + c_k^\mathrm{se} n_\mathrm{vc},
\end{equation}

\begin{equation}
c_k^\mathrm{e} = c_k^\mathrm{ne} + c_k^\mathrm{se} n_\mathrm{vc},
\end{equation}

where $n_\mathrm{vc}$ stands for the number of colors used in the corner tile set. Graphical illustration of the corner-edge tile equivalence is shown in Fig. <ref>.

$(c_k^\mathrm{nw}, c_k^\mathrm{sw})\qquad\qquad\qquad\quad$

$(c_k^\mathrm{sw},c_k^\mathrm{se})$$\qquad\qquad\qquad\quad(c_k^\mathrm{ne}, c_k^\mathrm{se})$
A corner tile expressed using edge formalism.

Data accessibility Source code available at: <https://gitlab.com/tyburec/tilopt>.

Authors' Contributions M.T.: conceptualization, methodology, software, writing – original draft, J.Z.: methodology,  writing – review & editing, supervision, funding acquisition.

Competing Interests The authors declare there is no competing interest.

Funding This research was funded by the Czech Science Foundation, project No. 19-26143X.

Acknowledgments We thank Stephanie Krueger for proofreading the initial draft of this manuscript.

[1]
Wang H. 1963  Dominoes and the AEA case of the decision problem. Symposium on the Mathematical Theory of Automata pp. 23–55.

[2]
Wang H. 1961  Proving theorems by pattern recognition–II. Bell System
Technical Journal 40, 1–41.

[3]
Berger R. 1966  The undecidability of the domino problem. Memoirs of the
American Mathematical Society.

[4]
Turing AM. 1937  On computable numbers, with an application to the
Entscheidungsproblem. Proceedings of the London Mathematical Society
s2–42, 230–265.

[5]
Davis M. 1958 Computability & Unsolvability.
Dover Books on Computer Science Series. Dover Publications.

[6]
Kahr AS, Moore EF, Wang H. 1962  Entscheidungsproblem reduced to the $\forall
\exists \forall$ case. Proceedings of the National Academy of Sciences
48, 365–377.

[7]
Lewis HR. 1978 pp. 35–47.
In Complexity of solvable cases of the decision problem for the predicate
calculus, pp. 35–47. IEEE.

[8]
Lewis H, Papimitriou C. 1981 Elements of the theory of computation.
Prentice-Hall Software Series. Pearson Education Canada.

[9]
Cohen MF, Shade J, Hiller S, Deussen O. 2003  Wang tiles for image and texture
generation. ACM Transactions on Graphics 22, 287–294.

[10]
Derouet-Jourdan A, Kaji S, Mizoguchi Y. 2019  A linear algorithm for brick
Wang tiling. Japan Journal on Industrial and Applied Mathematics
36, 749–761.

[11]
Lagae A, Dutré P. 2005  A procedural object distribution function. ACM Transactions on Graphics 24, 1442–1461.

[12]
Ollinger N. 2008  Two-by-two substitution systems and the undecidability of the
domino problem. In Logic and Theory of Algorithms pp. 476–485.
Springer Berlin Heidelberg.

[13]
Kovalsky SZ, Glasner D, Basri R. 2015  A global approach for solving
edge-matching puzzles. SIAM Journal on Imaging Sciences 8,

[14]
Lagae A, Dutré P. 2007  The Tile Packing Problem. Geombinatorics
17, 8–18.

[15]
Rui Yu CR, Agapito L. 2016  Solving jigsaw puzzles with linear programming.
In Richard C. Wilson ERH, Smith WAP, editors, Proceedings of the British
Machine Vision Conference (BMVC) pp. 139.1–139.12. BMVA Press.

[16]
Salassa F, Vancroonenburg W, Wauters T, Della Croce F, Berghe GV. 2017  MILP
and max-clique based heuristics for the Eternity II puzzle. arXiv:

[17]
Berger R. 1964 The Undecidability of the Domino Problem.
PhD thesis Harvard University.

[18]
Wang H. 1975  Notes on a class of tiling problems. Fundamenta
Mathematicae 82, 295–305.

[19]
Robinson RM. 1967  Seven polygons which permit only nonperiodic tilings of the
plane. Notices of American Mathematical Society p. 835.

[20]
Poizat B. 1980  Une théorie finiement axiomatisable et superstable. Groupe d'étude des théories stables 3, 1–9.

[21]
Knuth DE. 1968 The Art of Computer Programming, Volume 1: Fundamental
Addison-Wesley Educational Publishers Inc.

[22]
Knuth DE. 2018 The Art of Computer Programming, Volume 4B, Fascicle 5:
The: Mathematical Preliminaries Redux; Backtracking; Dancing Links.
Addison-Wesley Professional.

[23]
Robinson RM. 1971  Undecidability and nonperiodicity for tilings of the plane.
Inventiones Mathematicae 12, 177–209.

[24]
Grünbaum B, Shephard GC. 2016 Tilings and patterns.
Mineola, New York: Dover Publications.

[25]
Robinson RM. 1978  Undecidable tiling problems in the hyperbolic plane. Inventiones Mathematicae 44, 259–264.

[26]
Senechal M. 1996 Quasicrystals and Geometry.
Cambridge University Press.

[27]
Kari J. 1996  A small aperiodic set of Wang tiles. Discrete Mathematics
160, 259–264.

[28]
Čulík K. 1996  An aperiodic set of 13 Wang tiles. Discrete
Mathematics 160, 245–251.

[29]
Kari J, Papasoglu P. 1999  Deterministic aperiodic tile sets. Geometric
and Functional Analysis 9, 353–369.

[30]
Labbé S. 2019  A self-similar aperiodic set of 19 Wang tiles. Geometriae Dedicata 201, 81–109.

[31]
Jeandel E, Rao M. 2021  An aperiodic set of 11 Wang tiles. Advances in
Combinatorics pp. 1–37.

[32]
Lagae A, Dutré P. 2006  An alternative for Wang tiles: Colored edges
versus colored corners. ACM Transactions on Graphics 25,

[33]
Lagae A, Kari J, Dutré P. 2006  Aperiodic sets of square tiles with
colored corners. Report CW.

[34]
Nurmi T. 2016  From checkerboard to cloverfield: Using Wang tiles in seamless
non-periodic patterns. In Bridges Finland Conference Proceedings.

[35]
Kari J. 1990  Reversibility of 2D cellular automata is undecidable. Physica D: Nonlinear Phenomena 45, 379–385.

[36]
Conway J, Lagarias J. 1990  Tiling with polyominoes and combinatorial group
theory. Journal of Combinatorial Theory, Series A 53,

[37]
Mozes S. 1989  Tilings, substitution systems and dynamical systems generated by
them. Journal d'Analyse Mathématique 53, 139–186.

[38]
Stam J. 1997 Aperiodic texture mapping.
European Research Consortium for Informatics and Mathematics.
Tech. rep., R046.

[39]
Hiller S, Deussen O, Keller A. 2001  Tiled blue noise samples. In Proceedings of the Vision Modeling and Visualization Conference pp.
265–272 Stuttgart, Germany.

[40]
Radin C. 1987  Low temperature and the origin of crystalline symmetry. International Journal of Modern Physics B 1, 1157–1191.

[41]
Winfree E, Liu F, Wenzler LA, Seeman NC. 1998  Design and self-assembly of
two-dimensional DNA crystals. Nature 394, 539–544.

[42]
Seeman NC, Mao C, LaBean TH, Reif JH. 2000  Logical computation using
algorithmic self-assembly of DNA triple-crossover molecules. Nature
407, 493–496.

[43]
Novák J, Kučerová A, Zeman J. 2012  Compressing random
microstructures via stochastic Wang tilings. Physical Review E
86, 4–7.

[44]
Doškář M, Novák J, Zeman J. 2014  Aperiodic
compression and reconstruction of real-world material systems based on Wang
tiles. Physical Review E 90, 062118.

[45]
Doškář M, Novák J. 2016  A jigsaw puzzle framework for
homogenization of high porosity foams. Computers & Structures
166, 33–41.

[46]
Doškář M, Zeman J, Rypl D, Novák J. 2020  Level-set
based design of Wang tiles for modelling complex microstructures. Computer-Aided Design 123, 102827.

[47]
Tyburec M, Zeman J, Doškář M, Kružík M,
Lepš M. 2020  Modular-topology optimization with Wang tilings: An
application to truss structures. Structural and Multidisciplinary
Optimization 63, 1099–1117.

[48]
Tyburec M, Doškář M, Zeman J, Kružík M. 2022
Modular-topology optimization of structures and mechanisms with free material
design and clustering. Computer Methods in Applied Mechanics and

[49]
Jílek M, Somr M, Kulich M, Zeman J, Přeučil L. 2021  Towards a
passive self-assembling macroscale multi-robot system. IEEE Robotics and
Automation Letters 6, 7293–7300.

[50]
Doškář M, Zeman J, Jarušková D, Novák J.
2018  Wang tiling aided statistical determination of the Representative
Volume Element size of random heterogeneous materials. European Journal
of Mechanics – A/Solids 70, 280–295.

[51]
Coulais C, Teomy E, de Reus K, Shokef Y, van Hecke M. 2016  Combinatorial
design of textured mechanical metamaterials. Nature 535,

[52]
Yang W, Liu Q, Gao Z, Yue Z, Xu B. 2018  Theoretical search for heterogeneously
architected 2D structures. Proceedings of the National Academy of
Sciences 115, E7245–E7254.

[53]
Nežerka V, Somr M, Janda T, Vorel J, Doškář M,
Antoš J, Zeman J, Novák J. 2018  A jigsaw puzzle metamaterial
concept. Composite Structures 202, 1275–1279.

[54]
Knuth DE. 1992  Two notes on notation. The American Mathematical Monthly
99, 403–422.

[55]
Korte B, Vygen J. 2006 Combinatorial Optimization.

[56]
Gurobi Optimization, LLC. 2022  Gurobi Optimizer Reference Manual.