A cellular automata rule placing a maximal number of dominoes in the square and diamond

The objective is to demonstrate that a probabilistic cellular automata rule can place reliably a maximal number of dominoes in different active area shapes, exemplarily evaluated for the square and diamond. The basic rule forms domino patterns, but the number of dominoes is not necessarily maximal and the patterns are not always stable. It works with templates derived from domino tiles. The first proposed enhancement (Rule Option 1) can form always stable patterns. The second enhancement (Rule Option 2) can maximize the number of dominoes, but the reached patterns are not always stable. All rules drive the evolution by specific noise injection.


Introduction
Pattern formation is an area of active research in various domains such as physics, chemistry, biology, computer science or natural and artificial life. Cellular automata (CA) are suitable and powerful tools for catching the influence of the microscopic scale onto the macroscopic behavior of such complex systems [1][2][3]. At the least, the 1-dimensional Wolfram's "Elementary" CA can be viewed as generating a large diversity of 2-dimensional patterns whenever the time evolution axis is considered as the vertical spatial axis, with patterns depending or not on the random initial configuration [4]. Regarding the agent-based Yamins-Nagpal "1D spatial computer" [5,6] the authors emphasize therein how the local-to-global CA paradigm can turn into the inverse global-to-local question, namely "given a pattern, which CA rules will robustly produce it?" Such CA rules can be found by (i) proper design, (ii) by exhaustive search or (iii) by heuristics like Genetic Algorithm (GA) or Simulated Annealing, methods which were applied to solve the Density Classification Problem [7], for instance.
The arrangement of dominoes in a grid of cells is a special case of pattern formation. Possible applications are: parcel packing encountered in different logistics settings, such as loading boxes on pallets, arrangements of pallets in trucks or cargo stowage [8]; the design of a sieve for rectangular particles with a maximum flow rate; or an optimal arrangement of nanoparticles; and so forth. Unlike the dimer in statistical mechanics [9,10]-a pure tiling problem-we do not allow dominoes to contact with one another. That means that space between dominoes is mandatory. As a result, the solution space is more complex than in the dimer case with tight compaction. Domino arrangement is closely related to short-range interaction couplings in spin systems [11].

Previous and related work
In further previous work [12][13][14], different patterns were generated by agents with embedded finite state control which was evolved by GA. Matching pattern templates were also applied during the training period, but are not part of the CA rule as in our current proposal. They were also defined in a different simple way in order to count the number of dominoes for the fitness function during the evolutionary process.
In [15,16], domino patterns were formed by moving agents. Agents' behavior was controlled by a finite state machine, evolved by GA. The effort to find such agents was quite high, especially to find agents that work on any field size. In order to avoid such a computational effort, a novel approach to construct directly the required CA rule is proposed that will be presented thereafter. It has also the potential to be applied to further pattern formations. In addition, the new Rule Option 1 is defined that drives the evolution to a stable patterns.
In [17], this approach was already applied for the domino problem in a square field n × n , n even. Herein the purpose is to prove the robustness of the model against field's shape changing. Now the square size is generalized to any odd-even n. Moreover, the rule is now extended to a rhombic diamond field's shape, a 4 -tilted square. Again in this case, no confusion must be made with the Aztec diamond, another pure dimer tiling problem [18].
Parallel Substitution Algorithm (PSA) [19] is a powerful generalization of CA, which was also inspiring this work. PSA allows to substitute small locally defined patterns P by other patterns Q in a non-conflicting way. Thereby, very complex computations and transformations can be performed in a decentralized and parallel way.
The problem of optimal domino placement is presented in Sect. 2. A basic probabilistic CA rule is designed in Sect. 3 that can form valid domino patterns. In Sect. 4, two rule enhancements (options) are presented. Rule Option 1 allows to stabilize the pattern, whereas Rule Option 2 drives the evolution to optimal patterns. Results of simulation, performance evaluation and robustness are discussed in Sect. 5 before Conclusion.

The problem
Given is an array of N = (n + 2) × (n + 2) cells with state values s ∈ {0, 1} . The objective is to find a CA rule that can place a maximal number of dominoes within a given shape of active cells N active ≤ (n × n) surrounded by inactive border cells with value 0. A domino is given by two orthogonal adjacent cells of value 1. The constraint is that between dominoes there shall be at least one empty cell with value 0, i.e., dominoes are not allowed to touch each other.
The dominoes that we will use here in our solution are "domino tiles," and we propose to use them to cover the given shape. A domino tile consists of two pixels with value 1 (the kernel, in blue) and 10 surrounding pixels with value 0 (the hull, in green). In order to avoid confusion with the regular cells, we call the elements of a tile (or the later introduced templates) "pixels." Two types of dominoes are distinguished, the horizontal oriented domino ( D H ) and the vertical oriented ( D V ) (Fig. 1a). It is allowedand even necessary for a good solution-that green pixels from different domino tiles overlap. The possible levels of overlapping, from 2 to 4, are displayed in Fig. 1(b-d).
We call the level of overlapping cover level v.
Our intention is to show that the CA rule is robust (insensitive) against a change of the shape, and in addition, to improve the behavior of the CA rule (Option 1, in Sect. 4.1). Two shapes are considered, the square and the diamond. The diamond can be defined in the following way. Cells inside the diamond at position (x, y) (measured from a corner cell of the whole field) are given by the conditions The number of dominoes is denoted as where d H is the number of horizontal dominoes and d V is the number of vertical dominoes. A further requirement can be that the number of domino types shall be equal (or almost equal). We call is the maximal possible number of dominoes that can be placed into the field with overlapping.

Domino enumeration
We expose here a theoretical framework [20] to support the simulation results which will be presented in the sequel.

Dominoes in the square
Given a square array S n of (n + 2) × (n + 2) cells including a border with perimeter 4n + 4 enclosing a n × n field of order n 2 , the maximal domino number n (S n ) covering S n is given by the inductive formula n even: 0 = 0, 2 = 1, 4 = 4 and n = n−6 + 2(n − 2) n odd: 1 = 0, 3 = 2, 5 = 6 and n = n−6 + 2(n − 2) for n ≥ 6, where 2(n − 2) denotes the maximal number of dominoes in the crown surrounding the inner subgrid S n−6 . □ Setting m = ⌊n∕2⌋ and p = ⌊m∕3⌋ , this hierarchy of configurations can then be divided into the three equivalence classes S 0 , S 1 , S 2 according to m mod 3, as illustrated in Fig. 2.

Dominoes in the diamond
A diamond D n is given with perpendicular diagonals of length n and surrounded by a border with perpendicular diagonals of length n + 2.
• n odd . The order | D n | of D n fulfills including a border of length 2n + 2 and D n has an inner perimeter of length 2n − 2.
The maximal domino number n (D n ) is given by the inductive formula for n ≥ 6, where (n − 3) denotes the maximal number of dominoes in the crown surrounding the inner subgrid D n−6 . This hierarchy of configurations can be divided into the three equivalence classes D 0 , D 1 , D 2 according to m mod 3, as illustrated in Fig. 3.

3
A cellular automata rule placing a maximal number of dominoes… • n even . The order | D n | of D n fulfills including a border of length 2n + 4 and D n has an inner perimeter of length 2n. The maximal domino number n (D n ) is given by the inductive formula | D n | = 2 (2 + 4 + 6 + ⋯ + (n + 2)) = (n + 2)(n + 4)∕2 for n ≥ 6, where O(n) denotes the number of possible dominoes in the crown surrounding the inner subgrid D n−6 . The case n even is more intricate because the "crown" parameter is fluctuating and only fulfills a weak property. It is possible to define a generic family of diamonds as a complement of embedded squares, as illustrated in Fig. 4. The following three cases are considered (n ≥ 6).
In each case, the first term on the second side denotes the domino number in the square field × embedded into D n and the second term stands for the capacities in dominoes of the four remaining wedges. We should be aware that our construction can lead to a slight deficiency compared to the simulation. This can be explained by the fact that our theoretical construction is constrained by its own rule while the scenarios from the CA rule presented thereafter have more degrees of freedom. This deficit is due to the presence of isolated cells (Fig. 4) existing on some borders of our constrained system.

Space occupancy ratio
The space occupancy ratio k for a given domino k is defined as where k,i is the occupancy ratio (the inverse 1/v of the level of overlapping) of its pixel i. Since two dominoes cannot overlap, k,1 = k,2 = 1 always holds.
For illustration, the occupancy ratio of the surrounded tiles "k" in Fig. 5 yields: Evidently, the second and third configurations are optimal, but an overall configuration is constrained by the boundary conditions. For any field F of order | F | and containing d max dominoes, then Four typical arrangements of dominoes: (1) loosely coupled physical distancing configuration (2) tightly coupled orthotropic configuration (3) tightly coupled staggered configuration (4) tightly coupled isotropic configuration. Inset-The "domino tile" with numbered cells: the 2-cell kernel and the 10-cell hull for square S n and diamond D n , respectively. We observe that ratios | S n | ∕ n and | D n | ∕ n are decreasing Cauchy sequences which slowly converge towards the fixed limit of maximal occupancy as n approaches infinity and more generally for n even in any case. Our theoretical results are displayed in Table 1. Table 1 Domino enumeration for n × n fields with m = n∕2 and p = m∕3 Domino numbersn (S n ) in the squaren (D n ) in the diamond. Space occupancy ratio-| S n | ∕ n in the square-| D n | ∕ n in the diamond

The design idea
The first approach was to design a deterministic rule with synchronous updating.
After some experiments and experience from previous work, it showed to be very difficult if not even impossible to design such a rule that can converge always or with a high probability to the optimal or near-optimal aimed pattern. The second approach was to construct a probabilistic rule with synchronous updating. Indeed, such a rule was found for a field of size 6 × 6 by GA, where each cell is modeled as an agent that can turn in any direction. But the effort to find such rules is high and the good behavior cannot be guaranteed in general. The third and successful approach used here is the design of a probabilistic rule with asynchronous updating in a methodical way.

The basic rule
The basic idea is to modify the current configuration in a systematic way such that increasingly more dominoes appear and at last the CA evolves to a stable pattern. To do this, the CA configurations are searched for domino tile parts (specific local patterns), and if an almost correct tile part is found, it is corrected; otherwise, some random noise is injected.
The domino tile parts are called "templates" A i . They are systematically derived from the domino tiles (Fig. 6a). For each of the 12 pixels i of a domino (marked in red, carrying the domino pixel value dval(i)), a template A i is defined. A template can be seen as a copy of the tile, but shifted in space in a way that the pixel i corresponds to the center of the template.
In the computation, we represent a template as an array of size the reference pixel is called "reference value," refval(A i ) = val(A i , 0, 0) ∈ {0, 1} . Its value is equal to the red marked value of the corresponding tile pixel, refval(A i ) = dval(i) . The symbol # represents "Don't Care," meaning that a pixel with such a value is not used for matching (or does not exist, in another interpretation). Pixels with a value 0 or 1 are valid pixels; their values are equal to the values derived from the original tile. Some templates can be embedded into arrays smaller than (a � × b � ) when they have Don't Cares at their borders. Note that the valid pixels are asymmetrically distributed in a template because they are the result from shifting a tile.
Many of these templates are similar under mirroring, which can facilitate an implementation. For the vertical domino, a similar set of 12 templates is defined by 90 • rotation; altogether, we need 24 templates.
The templates A 7 − A 12 show white pixels that are not used because the template size (for the later described matching process) was restricted to (5 × 5) . As an example, the reduced template A 9 is marked in Fig. 6b by the blue square. The implementation with these incomplete templates worked very well, but further investigations are necessary to prove to which extent templates can be incomplete.
We need also to define the term "neighborhood template" that is later used in the matching procedure. The neighborhood template A * i is the template A i in which the reference value is set to #, in order to exclude the reference pixel from the matching process. The cell processing scheme is:  ). B * (x, y) denotes the states of the neighbors of (x, y) within a (5 × 5)-window, excluding the center. A new generation at time-step t + 1 is declared after N active cell updates (sub-steps) during the compute interval between t and t + 1 . We will update each cell once in a time interval in a random order which is re-computed for each new time-step.
The following rule is applied: The neighborhood templates A * i are tested against the corresponding cell neighbors B * (x, y) in the (5 × 5)-window at current position (x, y). Thereby, the marked reference position ( x, y) = (0, 0) of a neighborhood template is aligned with the center of the window. If all values match, then the state of cell (x, y) is set to the value refval(A i ) . Otherwise, with probability 0 , the cell is set randomly to either 0 or 1, or remains unchanged with probability 1 − 0 . The rule corrects the state of a cell to the reference value if a matching neighborhood is detected; otherwise, cell's state is changed randomly. It is possible that several neighborhood templates match (then tiles are overlapping), but there can be no conflicts because then all templates have the same reference value as derived from the tile. As no conflicts can arise, the sequence of testing the neighborhood templates does not matter, and one could skip further tests after a first match. It is important to note that the rule obeys the criterion of stability, which means that a field filled with dominoes without gaps (uncovered cells) is stable because we have matching at every site. Otherwise some random noise is injected in order to drive the evolution to an aimed pattern.

Test on the square
The rule was tested for N run = 1000 runs on ( 5 × 5)-fields with random initial configuration, T Limit = 50 (simulation time-step limit), with 0 = 0.5 . We know that there exist valid solutions with 6, 5, or 4 dominoes. The system converges after a few iterations into one of two solution classes. In the class I (stable), the 4 dominoes are covering the square totally without gaps and the reached pattern is stable. In the class II (partially stable), the 4 dominoes are covering the diamond not totally with at least one gap. The gap cells are toggling their state values ( 0 ↔ 1) due to the injected noise that never ends. Nevertheless, the class II patterns contains 4 dominoes which do not change (neither position, nor orientation or number).
The number of stable patterns (class I) with d = 6, 5, 4 dominoes was 41, 514, and 225, respectively. The average number of time-steps to reach stability for the class I patterns was t avrg = 4.0 . The number of partially stable class II patterns with an isolated toggling state was 194. The evolution of a class I and a class II pattern is depicted in Fig. 7b, c.

Test on the diamond
The rule was tested on ( 7 × 7 ) initially random fields for N run = 1000 , with T Limit = 50 and 0 = 0.5 . This diamond has the same number of active cells (25) as the (5 × 5)-square. We know that there exist only solutions with 4 dominoes. The system converges after a few iterations into one of the two solution classes; 166/1000 runs evolved into a stable class I pattern, and the rest evolved into class II patterns. The average number of time-steps to reach stability for the class I patterns was t avrg = 5.1 . A simulation of a class I and a class II pattern is depicted in Fig. 7d, e.
These tests showed that always either class I (stable) or class II (partially stable with gaps (isolated uncovered cells)) patterns evolved rapidly. The number of dominoes was ranging from the minimum to the maximum, also validated by further experiments not outlined here.

3 4 Rule enhancements
The basic rule was enhanced by adding optional rules on top: • (Option 1) The aim was to reach only stable class I patterns. • (Option 2) The aim was to reach preferably patterns with a maximal number of dominoes (max patterns).
A variable hit(x, y) was added to the cell's state then becoming (s, hit). All hits can be seen as an additional layer, also called hit matrix. The hit matrix stores the number of template hits for every site (x, y) that was selected for computation and updating.  • 100, if one neighborhood template with the reference value 1 matches. Note that 1-valued pixels are not allowed to overlap. The number 100 was arbitrarily chosen in order to differentiate such a hit from the other.
Recall that we are using asynchronous updating, where a time-step interval consists of N active sequential substeps. The value of a hit can be up-to-date, computed already during the current time interval, or it can be the old value from the previous time interval. In the case of a stable pattern the hit is equal to the cover level. In the case of a transient pattern a hit can be equal to the cover level in regions where the pattern shows already stable dominoes. So the hit matrix approximates the cover level values.
The enhanced execution mode is now for each cell: (1) compute the new state s ′ according to the basic rule, (2) compute the hit value, (3) modify the new state to s ′′ by Option 1 and then (4) modify the new state to s ′′′ by Option 2. Because of the asynchronous updating scheme, the new state replaces the state s then immediately.

Rule option 1: stabilizing the pattern
We have seen in Sect. 3.1 that the basic rule can evolve class I and class II patterns. We define now an additional rule that will turn class II patterns into class I patterns.
Analyzing the class II patterns, we can see isolated toggling cells with hit = 0 . The idea is to disseminate this information to the cells in the von Neumann neighborhood. If a cell in the neighborhood detects a hit-zero cell, it will produce additional noise in order to drive the evolution to a stable pattern without gaps. The optional rule is Two tests were performed with 100 runs and T limit = 300 , one on a (5 × 5)-square and another on a (7 × 7)-diamond.

Test on diamond
All 100 reached patterns were stable. The number of dominoes was d = 4 , the only possibility as expected, with t avrg = 37.0 (min 2 -max 197).

Rule option 2: maximizing the number of dominoes
The idea is to maximize the overlap between tiles by destroying non-overlapping situations ( hit = 0 ) through additional noise, allowing a reordering with higher hit s �� (x, y) = random ∈ {0, 1} with probability 1 ∃hit(x ± 1, y ± 1) = 0 s � (x, y) rates. First, the new state s ′′ (or s ′ if Option 1 is not applied) is computed and then the hit matrix. Then, the new state is modified to s ′′′ : When this option is applied it is not clear whether a reached domino pattern remains stable. In fact, it turned out that stability can only be reached if there exists a tiling, where every tile overlaps with at least another one or the border (called totally overlapping tiling), i.e., the cover level is v ≥ 2 everywhere inside the active area, e.g., a totally overlapping tiling exists for (10 × 10) but not for (8 × 8) square fields. Therefore, the number of dominoes will reach the maximum and remain stable in a (10 × 10) field, whereas the number of dominoes in a (8 × 8) field is reaching a maximum, and then it is decreasing and fluctuating and is driving again towards another maximum, and so forth.
The max patterns with d = 6 were stable. The patterns with d = 5 were changing between different solutions but did not reach a max pattern because of the limited time. These tests suggest that Option 1 is not really helpful when Option 2 is used in order to find a max pattern. One can see Option 1 only as another unnecessary source of noise. Therefore, in the following we will use only Option 2 because our objective is to find max patterns.

3
A cellular automata rule placing a maximal number of dominoes… 5 Performance and robustness

Performance for different field sizes
The basic rule with Option 2 (maximizing dominoes) was simulated for the square and diamond with 0 = 0.5 and 2 = 0.05 . The simulation time limit T limit was set large enough to yield max patterns for every simulation run in each test set of runs, e.g., T limit = 40 000 for n = 12, 16 , two slow converging cases. The number of runs in a test set for averaging was N run = 100 for n ≤ 10 and 20 for n > 10 . The average time t avrg is shown in Table 2. It corresponds to the Work per Cell if each cell is considered as a processor, because the total (parallel) Work is t avrg × N active . We divide further the work per cell by the problem size N active . This quotient E = t avrg ∕N active is a constant if the work per cell would increase linear with N active . So this measure can show us a superlinear increase in work per cell if it increases with problem size. The minimal and maximal times were also recorded for each run and related to the average time. On average, the minimal/ maximal time recorded was 0.085/4.02 times t avrg . These extensive simulations confirmed that the basic CA rule with Option 2 converges to patterns with a maximal number of dominoes if the time limit is large enough, no matter whether the square shape or the diamond shape was used. It was surprising that E shows some remarkable peaks and drops, in addition of a weak increase with problem size. In the investigated range E was highest for n = 12 in the diamond and for n = 18 in the square. An explanation is that max patterns have a certain frequency among all valid patterns, and that there a harder tasks (field sizes) where the frequency of max patterns is very low and therefore more difficult to reach by the walk through the pattern space. For example, there are only 4 highly symmetric optimal solutions in the diamond D 12 but a lot of non-symmetric optimal solutions in D 14 . In Fig. 4, only symmetric solutions were shown. Further research is necessary to make a sound statement about the time complexity for large n. Figure 8 shows some patterns evolved by the CA rule with Option 2 for the square and the diamond for n = 12, 13 . All are max patterns except (d). (a, c, e, g) are balanced. (a, b, c, f) are totally overlapping (no v = 1 is visible). Some larger balanced max patterns evolved with Option 2 are shown in Fig. 9.
Our theoretical results displayed in Table 1 should be compared with simulation results in Table 2. | S n−2 | and | D n−2 | in Table 1 correspond with N active in Table 2 while n and n tie in d max . Let us observe the slight deficit for 16 , 17 . Looking at Fig. 9, let us observe a slight deficit for 24 , 25 and a big deficit for 26 . Looking back to Fig. 4, let us observe again the presence of isolated cells on the diamond's border, that is likely to explain this phenomenon of domino's deficit.

Robustness
We have seen that the rule works reliable for two different shapes, and it works also for the circle and rectangle when being tested. Further experiments have shown that the CA rule is also robust (insensitive) against the initial configuration and the order in which the cells are updated.  a, c, e, g are balanced. a, b, c, f are totally overlapping (cover level v > 1)

3
A cellular automata rule placing a maximal number of dominoes…

Initial configuration
In all the simulations conducted here, the initial cell states were random. As already shown in [17], the rule works similarly well if the initial configuration is black or white everywhere.

Self-stabilization
This term was introduced in 1973 by Dijkstra [21,22]. It means that a system will always converge to a desired system state even if it is disturbed. This is a very important feature for systems to be reliable. When in our system failures appear from time to time, they will not influence the overall dynamics, because if the system is in an early stage it will result in just another random initial state, and if the system is in a late stage, then errors are corrected or the system will drive to another valid solution.

Updating sequence
We used here for each time interval a different random order in such a way that every cell is updated once. In the former work [17], there was no ordering, cells were picked up at random and updated. Thereby, a cell can be updated never or multiple times in a time interval (updating N active times). By this method, the domino patterns evolved similarly.
In addition we tested also the row-by-row sweep sequence and the "checkerboard sequence." For the checkerboard sequence, first the cells are updated for (x + y) even, then for (x + y) odd. Also these tests suggest that the cell-by-cell updating order does not significantly influence the overall evolution. An explanation is that the rule is probabilistic, and any deterministic order is destroyed by the rule's inherent randomness. This is useful when the computation is executed on a parallel supercomputer with several computing nodes / cores. Then all cores can compute in parallel and the order of processing the border cells between the subareas (including data exchange between neighboring cores) is not critical and can be scheduled in a way that the slowdown is minimal.

Conclusion
A probabilistic CA rule with two options was designed that can form high quality domino patterns. The basic rule uses 24 matching templates derived from the two (3 × 4) domino tiles. Each selected cell is tested against the templates and is adjusted in the case that a template matches in the neighborhood. The basic rule is sub-optimal with respect to the number of placed dominoes. The evolved patterns are not always stable if there exist isolated gaps. Rule Option 1 distributes the gap information to the neighborhood, then the additional noise drives the evolution to stable patterns. Rule Option 2 injects noise where there is no overlap (overlap level 1) which drives the evolution to a maximal number of dominoes. A reached optimal pattern remains stable, if it is totally overlapping (no overlap level 1 exists). The basic rule and its options are robust against the active shape (square, diamond) and the order in which the cells are updated.
Funding Open Access funding enabled and organized by Projekt DEAL..

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.