Bounded Wang tilings with integer programming and graph-based heuristics

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 provide 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 ∀∃∀ 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 N P-complete in general, e.g., [7] or [8,Theorem 7.2.1],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 1(a) and applications of Wang tiles in Section 1(b).In Section 1(c), we list available algorithms for generation of Wang tilings.Finally, our aims and contributions appear summarized in Section 1(d).

Aperiodic 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 [24,Chapter 11].Starting from the Berger tile set containing 20,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. 1 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 nonexistence 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 A 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-tilebased 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 T , substitution is a map S : T → T that assigns a tiling T k to each tile k ∈ T ; we refer the reader to Section 2 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 1D 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 2 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 2 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 3(a), maximum rectangular valid tiling in Section 3(b), maximum-cover in Section 3(c), and maximum adjacency constraint satisfaction in Section 3(d).To allow for a finer control over the resulting tilings, we also include simple Due to the complexity of the proposed formulations, in Section 4 we propose a heuristic graph-based algorithm to tackle the maximum-cover optimization variant from Section 3(c).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 5 collects results from the computational assessment of the integer programming formulations (Section 5(a)) and heuristics (Section 5(b)), and on the benchmarking of the periodic tile packing formulation against the algorithm of Lagae and Dutré [14] (Section 5(c)).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 simplyconnected valid tilings, Section 5(d), and the Lagae et al. [33] tile set of 44 corner tiles is not aperiodic, Section 5(e).We summarize our results in Section 6.

Notation and preliminaries
Considering a finite set of color codes C = {1, 2, . . ., n c } ⊂ N, the (Wang) tile k is a quadruple of the color codes (c n k , c w k , c s k , c e k ), with c n k , c w k , c s k , and c e k ∈ 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. 2a.Without loss of generality, we further consider these squares to be of the unit size.
A tile set T represents a finite collection of n t tiles, see Fig. 2b.When ∀(c n , c w , c s , c e ) ∈ C 4 : (c n , c w , c s , c e ) ∈ T , we call the tile set complete.
Using the notation • = • ∩ Z 2 to denote an intersection of the set • with the integer lattice points, tiling T A of a bounded domain A ∈ R 2 is an arrangement of copies of the tiles from T such that the tiles are placed at Ã, and cover the entire domain A, cf.Fig. 3.More formally, tiling is a mapping T A : Ã → T assigning a single tile k ∈ T to every coordinate (i, j) ∈ Ã.Consequently, we call tilings T A simply connected iff the domain A is so.
The tiling T A is rectangular if ∀i ∈ H, H = {1, . . ., n h }, and ∀j ∈ W, W = {1, . . ., n w }, it holds that (i, j) ∈ Ã.Here, H and W are the sets of the height and width coordinates.
A valid tiling (Wang tiling) of A, denoted by T A valid , is a tiling with equal color codes at the shared edges between all pairs of adjoining tiles.Therefore, the mapping T A valid : Ã → T satisfies, in addition to the requirements for T A , the additional constraints provided that the axes are oriented accordingly to Fig.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 T is periodic if it permits periodic valid tilings; and aperiodic if all feasible valid tilings are aperiodic. Transducer For the dual transducer graph G t,v , composed of the dual Wang tiles [30] reflecting T along the major diagonal of the tiles, the edge set is defined as To illustrate the construction, we include a visual example in Fig. 4.

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

Rectangular valid tiling
Let us now consider the fundamental problem of finding T A v or proving it does not exist.From now on, we restrict A to be rectangular to simplify notation.However, the presented approach also extends to the general case.
To achieve this, we introduce ∀(i, j, k) ∈ H×W ×T a binary decision variable x i,j,k ∈ {0, 1} denoting the placement of the tile k at the (i, j) coordinate such that Consequently, mapping T A (i, j) is expressed as together with the requirement that every (i, j) coordinate is occupied by one tile, Similarly, the color codes of a tile placed at (i, j) are expressed using the binary variables as Inserting ( 7) into (1a) and (1b) leads to the horizontal and vertical adjacency constraints expressed in terms of the decision variables, as Combining ( 4), ( 5), (6), and ( 8) then provides us with a complete binary linear programming representation of valid tiling T A valid .
For computational reasons, it proved to be advantageous to organize the constraints according to the color codes: where, in the Iverson notation [54], k∈T x i,j,k [c s k = ] expresses that x i,j,k is added to the sum if and only if c s k = .The constraint (9a) requires that the number of tiles at (i, j) with the south edge colored by equals to the number of tiles at (i + 1, j) with the north edge marked by the same , for all ∈ C. Because of ( 6), there are either no tiles with the shared edge colored by , or there is a single tile at each of the coordinates with its common edge labeled by .Analogously to the vertical adjacency constraint, the horizontal constraint (9b) also enforces equality among the number of tiles at (i, j) with the east edge colored by and the number of tiles at (i, j + 1) having the west edge colored by identical .
Finally, combining ( 4), ( 6), and ( 9), while noticing that the constraints (6) naturally propagate with the adjacency constraints from the domain boundaries (compare (10d,10e with ( 6)), leads to the binary programming formulation find x (10a) 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 (10) are valid tilings.Moreover, observe that the problem consists of two totally unimodular constraints if considered independently: (10c,10e) representing row tilings, and (10b,10d) being column tilings.When considered simultaneously, the resulting problem becomes N P-complete [7,8].

Maximum rectangular valid tiling
When a solution to (10) 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 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., On the other hand, all the other coordinates may contain a tile or be empty, thus 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 simplyconnected 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. 5a-5d.Formally stated in terms of the decision variables, these considerations are expressed as Similar arguments hold also for the coordinates (i, j) and (i, j +1), resulting in the constraints The allowed and forbidden combinations for this case are shown in Fig. 5e-5h.
The developed constraints ( 11)-( 14) 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. 6.Because the combination 6b cannot appear in any simply-connected rectangular tiling, we must exclude it from the feasible set, Finally, combining Eqs. ( 4), ( 11), ( 12), ( 13), (14), and ( 15) together with an objective function to maximize | Bmax rect | provides us with the binary maximum rectangular valid tiling optimization program In contrast to (10), a feasible solution to the optimization program ( 16) can be found in a polynomial time, e.g., by tiling the first row or column of the 1D bounded tiling problem.However, finding an optimal solution to ( 16) is N P-hard, because the optimization problem ( 16) is reducible to the decision version (10) by fixing the value of the objective function to | Ã|, which enforces equalities in (16b), (16c), and (16f), making the constraint (16d) redundant as a consequence.

Maximum cover
Another option for avoiding the infeasibility of (10) 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 A, or equivalently a valid tiling of the (possibly disconnected) domain B max cov ⊆ 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 k∈T x i,j,k = 0 for a void located at (i, j) ∈ Ã.Thus, each coordinate (i, j) is occupied either by a tile or a void, implying that and the vertical and horizontal edge matching conditions become Finally, the combination of Eqs. ( 17), (18a), (18b) with the objective function to maximize | Bmax cov | leads to the binary optimization problem The program ( 19) is trivially N P-hard: Requiring the objective function (19a) to be at least i.e., all positions are occupied by a Wang tile.Moreover, (19b) and (19c) 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 N P-complete bounded tiling problem.

Maximum adjacency constraints satisfaction
Because the decision problem (10) 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 v i,j ∈ R ≥0 , where (i, j) ∈ H × W \ n w , and h h i,j ∈ R ≥0 , with (i, j) ∈ H \ n h × W, respectively.The adjacency constraints (9) are then relaxed by considering instead.Indeed, if h h i,j = 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 v i,j = 0 guarantees color matches among the tiles at (i, j) and (i, j + 1).
Finally, rewriting absolute values in ( 21) by two linear inequalities while supplying an objective function to maximize the number of color matches yields the binary optimization problem max that is N P-hard due to the reduction to (10) after setting all h v i,j and h h i,j 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 i∈H j∈W\nw h v i,j or i∈H\n h j∈W h h i,j 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 (variablesized) 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): Similarly, we may prevent tile k from being placed there: Placement of an identical tile at the coordinates (i, j) ∈ Ã and (p, q) ∈ Ã requires Conversely, different tiles at these coordinates are secured with

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) ∈ Ã is set to by On the contrary, we may prevent this color by requiring Further, the same color codes at the north edge of (i, j) ∈ Ã and at the west edge of (p, q) ∈ Ã are established with and a different color with

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 A: examining periodicity through periodic color-based boundary conditions.We begin with requiring equal coloring at the fixed opposite domain boundaries, When adding (31) to the decision problem (10), 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 (16).Naturally, the domain size is not known in this case.Therefore, we must consider ∀(i, j, ) ∈ H × W × C constraints of the form Here, (32a) prevents a color mismatch of the north edge of (1, j) ∈ Ã and the south edge of (i, j) ∈ A iff there is no tile placed at (i, j + 1) ∈ Ã.Similarly, in the case of (32b), we prevent a color mismatch of the west edge at (i, 1) ∈ Ã and the east edge at (i, j) ∈ Ã iff the position Finally, when adding the constraints (32) to ( 16), we usually search for the smallest periodic pattern rather than the largest,

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, Note here that this extension requires that |T | = | Ã| 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 (10).In this formulation, neglecting any pair of the constraints (10b, 10d) or (10c,10e) provides a totally unimodular constraint matrix, recall Section 3(a).Consequently, such simplified problems are deterministically solvable using the simplex method.Moreover, this setting agrees with the maximum flow problem [55], as (10d) and (10e) correspond to the flow balances in the source and sink, whereas (10b) and (10c) 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.As follows from Section 2, any valid tiling of a row can be visualized as a |W|-long path in the transducer graph G t,h , recall Section 2. To simplify subsequent developments, we represent the row-tiling problem by a transducer-based directed acyclic graph (DAG) composed of |W| + 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 |C| 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 E h ; see Fig. 7. Consequently, any s → t path in the yet-established directed graph forms a valid tiling of the row, and conversely, any valid tiling builds a s → 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 |W||C| 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. 8. Consequently, we generate at most 2|W||C| new edges altogether.
In addition, we assign unitary costs to the edges incoming to the void vertices and zero costs elsewhere.Hence, the s → 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 O(|V| + |E|) time using the DAG-shortest-path algorithm [55], where V denotes the set of the graph vertices and E the set of the graph edges.In our case, we have Thus, the overall asymptotic complexity to generate a maximum row cover tiling evaluates as 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 t contains exactly c 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.
Proposition 4.1.The shortest path in the graph in Fig. 8 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. 8.
For this, we first check the north-east compatibility for each tile k ∈ 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 = 1/2(|W| + 1) −1 to the tiles that admit a single vertical neighbor only, and = (|W| + 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 |W|/(|W| + 1) < 1, i.e., the maximum number of tiles is placed even if the void positions forbid any vertical neighbors.Hence, Proposition 4.1 remains satisfied.

1/2-approximation algorithm for general tile sets
In this section, we modify Alg. 1 to maintain the 1/2 approximation ratio.We start with the following observation: , where • rounds • to the nearest greater or equal integer.
To exploit Proposition 4.2 in Alg. 1, we modify the row processing order to {1, 3, 2, 5, 4 . . .}. Indeed, then each odd row contains exactly | Bmax rowcov | 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 = 1/2(|W|+1) −1 as before.The final algorithm then reads as Alg. 2, allowing us to state the following, slightly stronger result:

2/3-approximation algorithm for tilesets with cyclic transducers
Another improvement in the approximation factor of Alg. 2 is possible for tile sets with all the states in the transducer graphs G t,h and G t,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, . . .}.We begin with (i) tiling the maximum row-cover of the first row.Then, we (ii) find the maximum row-cover of the 4th row such that we penalize possible incompatibilities with the first row based on the dual transducer graph by .The step (iii) encompasses finding a cover of the 3rd 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. 3.Then, we can make the following statement: Lemma 4.1.Consider that all states in the transducer graphs G t,h and G t,v are in at least one graph cycle.Then, Alg. 3 terminates with at least 2  3 | Ã| placed tiles.Proof.Since the tile set allows for valid tiling of the row, the {1, 4, . . .} rows are occupied by exactly |W| tiles.The {3, 6, . . .} rows are then populated by at least 1/2|W| tiles because each tile from rows {4, 7, . . .} admits a vertical neighbor.Finally, the {2, 5, . . .} rows contain at least the complement of the number of tiles used in the preceding row, because the tiles in the {1, 4, 6, . . .} row admit a south neighbor.Depending on the number of rows, the algorithm places at least tiles.

Iterative improvements
Similarly to finding the maximum row covers, we can search for the maximum cover of columns.

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 5(a).Then, in Section 5(b), we also relate these results to the outputs of the heuristic algorithms.
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 ® Core ® i5-8350U CPU clocked at 1.70GHz.

Integer programming formulations
In this section, we investigate the performance of all integer programming formulations from Section 3, i.e., the decision program (10), the maximum rectangular tiling ( 16), the maximum cover (19), and the maximum adjacency constraint satisfaction problem (22).
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 1(a).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 Algorithm 4 Final maximum cover heuristics 1: function FINALMAXIMUMCOVERHEURISTICS(T , A) T ← generateInitialCover(T ,A) 3: method ← "columns" if method=="rows" then 9: for row ← {1, . . ., |H|} do 10: T ← updateTiling(T, shortestPath, row) T ← updateTiling(T, shortestPath, column) return T 25: end function 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. 9, 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×20, 25×25, and 30 × 30.The running time of the Gurobi solver was limited to 300 seconds for the singlethreaded mode.
The results shown in Tab. 1 illustrate that the performance of the decision program (10) 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 (16,19,22)  provided at least some output.Interestingly, the decision problem (10) also was more efficient in the case of proving that the domain |A| lacks T -tilability.
Comparison of the optimization variants hints that the maximum cover (19) and the maximum adjacency constraint satisfaction (22) problems scale better than the maximum rectangular tiling (16).Indeed, generating any smaller rectangular domain remains N P-complete, preventing any polynomial-time approximation algorithm to exist.On the other hand, both the formulations (19) and ( 22) admit simple heuristics, recall Section 3, allowing the solver to obtain higher-quality feasible solutions faster.

Heuristic algorithms
Second, we compare the performance of the maximum cover formulation (19) solved with the heuristic Alg. 4 supplied with three different initial coverings, i.e., based on Algs. 1, 2 and 3.
Alg. 4 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. 4 100 times for each of the tested option, and listed the best, worst, and mean results in Tab. 2.
From Tab. 2, it follows that the initialization with the cover from Alg. 1 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. 1, at least 82% of tiles were always placed, only more than 60% followed from Alg. 2. Using Alg.placement.When comparing Tab. 1 with Tab. 2, 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. 1, it can be seen that the resulting covers are very competitive to the outputs of (19) and also obtained in much shorter times.

Periodic tile packing problem
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.proposed in Lagae and Dutré [14] to find a feasible solution.
While Lagae and Dutré [14] proposed a backtracking-based algorithm to generate periodic packings, we rely here on a solution to the decision program (10) supplemented with the packing (34) and fixed periodicity (31) constraints.The resulting core times spent in the search for a single feasible solution (Tab.3) 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 [21, Exercise 5 in Section 2.3.4.3].Generating valid tilings from the Knuth tile set using the decision program (10) together with the tile-based boundary conditions, recall Section 3(e)(i), led to an unexpected observation that enforced placement of the βU S tile makes the program (10) infeasible under certain circumstances.
After a careful investigation, it indeed turned out that there is not any 2 × 2 valid tiling with the βU S tile placed at (2, 2).Moreover, there is also not any 4 × 3 valid tiling with the βU S tile placed at (3, 1).Thus, using the maximum-cover optimization variant (19) and the β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 β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 [22, Exercise 221 in Section 7.2.2.1].

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 A. For these tiles, we solve the rectangular tiling formulation (16) with periodic boundary conditions (32) and an objective function to find the smallest tiling (33).As its output, we receive the optimal value of 6 and 12 optimal periodic rectangular tilings of the size 2×3.Not surprisingly, all these solutions follow from only two periodic patterns shown in Fig. 11 by translations over the infinite plane.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.
Lemma 5.1.The horizontal translation method from Lagae et al. [33] is bijective iff the dual transducer graph G T,v of the input tile set T does not contain any parallel arcs.
Proof.The horizontal translation method is formally a mapping T × T → T corner that generates ∀(p, q) ∈ T 2 : c e p = c w q a corner tile (c n p , c s p , c s q , c n q ).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 ∈ T in the original tile set must be uniquely determined by c n p and c s p , as the color codes of the vertical edges of T are avoided in the construction of T corner .
Let us now consider that the dual transducer graph contains a parallel arc connecting the state c n with c s .Then, there may exist two tiles colored by (c n , c w p , c s , c e p ) and (c n , c w q , c s , c e q ) that are indistinguishable in T corner , which contradicts the bijection.For the other option, if the transducer graph does not contain any parallel arcs, then each c n q , c s q identifies with a single arc labeled by c w q |c e q , i.e., with a single tile, which completes the proof.
Rotating the tile set by 90 degrees, the arguments in Lemma 5.1 provide us with the conditions for the bijectiveness of the vertical translation method: Lemma 5.2.The vertical translation method of Lagae et al. [33] is bijective iff the transducer graph G T,h of the input tile set T does not contain any parallel arcs.
For the Ammann tile set, we obtain the transducer graph G T,h = G T,v shown in Fig. 12.Clearly, there exist parallel arcs 1 → 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.as well as its N P-hard optimization variants relaxing some of the initial assumptions: tilability of the entire rectangular domain leading to the maximum rectangular tiling formulation (16), simple-connectedness to the maximum-cover program (19), and tiling validity to the maximum constraint satisfaction problem (22).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 N P-hardness of all optimization formulations, we developed simple yet efficient heuristic algorithms for the maximum-cover variant (19).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 × 20, 25 × 25, and 30 × 30 revealed that the decision program (10) 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 (19) and maximum adjacency constraint satisfaction problems (22) appear to be similarly efficient.The remaining formulation-maximum rectangular tiling ( 16)-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 (19) for 300 s.Testing revealed that the simplest setup of Alg. 4 initialized with the cover generated by Alg. 1 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.This is, however, also maintained if we denote their edges by labels in the form of tuples of the corner codes, (c nw k , c ne k ), (c nw k , c sw k ), (c sw k , c se k ), and (c ne k , c se k ), 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 where n 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. 13.

Figure 2 :
Figure 2: Graphical representation of (a) a Wang tile k, and of (b) a tile set T .

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

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

Proposition 4 . 2 .
Consider the maximum row-cover tiling of the odd rows of the initially void domain A given in Section 4(a).Then, | Bcov | ≥ 1/2| Bmax cov |.Proof.Consider that the maximum row-cover problem alone terminates with | Bmax rowcov | tiles.Based on the maximum row-cover property in Proposition 4.1, none of the rows of Ã admit a tiling by more than | Bmax rowcov | tiles.Hence, we have

Proposition 4 . 3 .
Assume a tile set T with the longest path in its transducer graph G t,h of at least 2.Then, Alg. 2 terminates with | Bcov | ≥ 1/2| Ã|.Proof.When | Bmax rowcov | = |W|, the proof follows directly from Proposition 4.2.For the other cases, the odd rows must contain | Bmax rowcov | tiles due to Proposition 4.1.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 | Bmax rowcov | + (|W| − | Bmax rowcov |) , 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 t,h is at least two, we can always place tiles to the north of the voids of the i-th row.
end function | Ã| 2 |T | + |C| 2 ) complexity and provides the approximation ratios adjustable by algorithm choice (Algs.1, 2 or 3) at line 2 of Alg. 4. Proposition 4.4.Alg. 4 runs in a polynomial time and terminates in a finite number of steps.Proof.We have already shown that finding a maximum row-cover has O(|W||C| + |W||T |) complexity.Further, finding the 2-long paths in the transducer graph possesses the |C| 2 complexity and can be run only once prior to the algorithm main loop.Altogether, Alg. 3 requires at most 4/3|H| inner iterations so that we have the O(| Ã||C| + | Ã||T | + |C| 2 ) overall complexity.Regardless of the method at line 2 of Alg. 4, the improving loop runs at most | Ã| times.Consequently, the algorithm is finite and possesses the O(| Ã| 2 |C| + | Ã| 2 |T | + |C| 2 ) complexity. )

Figure 9 :
Figure 9: 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.

Figure 10 :
Figure 10: Periodic packing of a complete set of 625 tiles over 5 colors.

Figure 11 :
Figure 11: Rectangular periodic valid tilings.Translating a 2 × 3 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 × 3.

Figure 12 :
Figure 12: Transducer graph of the Ammann set of 16 Wang tiles over 6 colors.

Figure 13 :
Figure 13: A corner tile expressed using edge formalism.A Corner tiles represented as Wang tilesEach corner tile in the corner tile set is defined by a quadruple of color codes(c nw k , c sw k , c se k , c ne k ), with c nw k , c sw k , c se k ,and cne k 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 nw k , c ne k ), (c nw k , c sw k ), (c sw k , c se k ), and (c ne k , c se k ), 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 3. If such T A valid exists, we say that the domain A admits a valid T -tiling, or that it is tileable by T .Consider that B ⊆ A and B max rect ⊆ A are simply connected, rectangular, and T -tileable.Then, the maximum rectangular valid tiling T A v,max rect is a valid tiling of the domain B max rect , where {B max rect ⊆ A, ∀B ⊆ A : | Bmax rect | ≥ | B|}.Here, the notation |•| denotes cardinality of the set •.The maximum cover T A v,max cov is a valid tiling of B max cov , where B and B max cov are arbitrary T -tileable subdomains in A and {B max cov ⊆ A, ∀B ⊆ A : | Bmax cov | ≥ | B|}.
graph[27] G t,h of the tile set T is a directed (multi-)graph representation of a Mealy machine without any initial nor terminal state.It consists of |C| states (graph vertices) and |T | transitions (directed edges) E h , where

Table 1 :
Benchmark results.Values marked by an asterisk denote a premature termination of the integer programming solver.
3, we obtained at least 70% tile

Table 2 :
Numerical tests of the maximum-cover heuristics, Alg. 4, initialized based on Algs. 1, 2, and 3. Best mean runs are highlighted in bold.

Table 3 :
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)