The problem of out-competing the SAT solver with MultiwayDeletionsGraph apparently took me about 3-4 hours to knock out. Here are the results, starting with the SAT solver benchmark:

And here are the results from MultiwayDeletionsGraph:

Again, the MDG times are an order of magnitude faster and we get the same results in terms of cardinality for completed atlases $2\times 2$ or $3 \times 3$. However, in this case it was necessary to impose a fairly strict order which would cause the MDG to be a tree rather than a DAG. Especially for the $3 \times 3$ case it makes no sense to place a next tile unless it is adjacent to one of the previous tiles.

Let's generate tiles with SatisfiableQ SAT solving!

Those Wang tiles are aperiodic for fewer than 11 tiles or fewer than 4 colors.

Yes, are the parameters for the faster GenerateWangTiling like expanded capabilities for dealing with more than two colors?

We should test out the boundaries on the graph. And change edge matching rules into an exact cover matrix.

SAT solvers are so cool. In Python (brew install python@3.7), the help of some utility functions and R. Hettinger's SAT Solver we can solve Einstein's Puzzle in no time.

from pycosat import itersolve   # https://pypi.python.org/pypi/pycosat
from itertools import combinations
from sys import intern
from pprint import pprint


# SAT Utility Functions solve_one, from_dnf, one_of, https://github.com/redmage123/problem_solvers/blob/master/src/sat_utils.py.


class Q:
    'Quantifier for the number of elements that are true'

    def __init__(self, elements):
        self.elements = tuple(elements)

    def __lt__(self, n: int) -> 'cnf':
        return list(combinations(map(neg, self.elements), n))

    def __le__(self, n: int) -> 'cnf':
        return self < n + 1

    def __gt__(self, n: int) -> 'cnf':
        return list(combinations(self.elements, len(self.elements)-n))

    def __ge__(self, n: int) -> 'cnf':
        return self > n - 1

    def __eq__(self, n: int) -> 'cnf':
        return (self <= n) + (self >= n)

    def __ne__(self, n) -> 'cnf':
        raise NotImplementedError

    def __repr__(self) -> str:
        return f'{self.__class__.__name__}(elements={self.elements!r})'


def neg(element) -> 'element':
    'Negate a single element'
    return intern(element[1:] if element.startswith('~') else '~' + element)


def one_of(elements) -> 'cnf':
    'Exactly one of the elements is true'
    return Q(elements) == 1


def from_dnf(groups) -> 'cnf':
    'Convert from or-of-ands to and-of-ors'
    cnf = {frozenset()}
    for group in groups:
        nl = {frozenset([literal]): neg(literal) for literal in group}
        # The "clause | literal" prevents dup lits: {x, x, y} -> {x, y}
        # The nl check skips over identities: {x, ~x, y} -> True
        cnf = {clause | literal for literal in nl for clause in cnf
               if nl[literal] not in clause}
        # The sc check removes clauses with superfluous terms:
        #     {{x}, {x, z}, {y, z}} -> {{x}, {y, z}}
        # Should this be left until the end?
        sc = min(cnf, key=len)          # XXX not deterministic
        cnf -= {clause for clause in cnf if clause > sc}
    return list(map(tuple, cnf))


def make_translate(cnf):
    """Make translator from symbolic CNF to PycoSat's numbered clauses.
       Return a literal to number dictionary and reverse lookup dict
        >>> make_translate([['~a', 'b', '~c'], ['a', '~c']])
        ({'a': 1, 'c': 3, 'b': 2, '~a': -1, '~b': -2, '~c': -3},
         {1: 'a', 2: 'b', 3: 'c', -1: '~a', -3: '~c', -2: '~b'})
    """
    lit2num = {}
    for clause in cnf:
        for literal in clause:
            if literal not in lit2num:
                var = literal[1:] if literal[0] == '~' else literal
                num = len(lit2num) // 2 + 1
                lit2num[intern(var)] = num
                lit2num[intern('~' + var)] = -num
    num2var = {num: lit for lit, num in lit2num.items()}
    return lit2num, num2var


def translate(cnf, uniquify=False):
    'Translate a symbolic cnf to a numbered cnf and return a reverse mapping'
    # DIMACS CNF file format:
    # http://people.sc.fsu.edu/~jburkardt/data/cnf/cnf.html
    if uniquify:
        cnf = list(dict.fromkeys(cnf))
    lit2num, num2var = make_translate(cnf)
    numbered_cnf = [tuple([lit2num[lit] for lit in clause]) for clause in cnf]
    return numbered_cnf, num2var


def lookup_itersolve(symbolic_cnf, include_neg=False):
    numbered_cnf, num2var = translate(symbolic_cnf)
    for solution in itersolve(numbered_cnf):
        yield [num2var[n] for n in solution if include_neg or n > 0]


def solve_one(symcnf, include_neg=False):
    return next(lookup_itersolve(symcnf, include_neg))


# SAT Solver Solution to the Einstein Puzzle, https://rhettinger.github.io/einstein.html.


houses = ['1',          '2',      '3',           '4',         '5']

groups = [
    ['Red',    'Green',    'Ivory',       'Yellow',     'Blue'],
    ['Englishman',      'Spaniard',   'Ukrainian',       'Japanese', 'Norwegian'],
    ['Coffee',     'Tea',    'Milk',        'Orange juice',    'Water'],
    ['Old Gold', 'Kools', 'Lucky Strike', 'Parliament',   'Chesterfield'],
    ['Dog',     'Snails',    'Zebra',        'Horse',      'Fox'],
]

values = [value for group in groups for value in group]


def comb(value, house):
    'Format how a value is shown at a given house'
    return intern(f'{value} {house}')


def found_at(value, house):
    'Value known to be at a specific house'
    return [(comb(value, house),)]


def same_house(value1, value2):
    'The two values occur in the same house:  Englishman1 & red1 | Englishman2 & red2 ...'
    return from_dnf([(comb(value1, i), comb(value2, i)) for i in houses])


def consecutive(value1, value2):
    'The values are in consecutive houses:  green1 & white2 | green2 & white3 ...'
    return from_dnf([(comb(value1, i), comb(value2, j))
                     for i, j in zip(houses, houses[1:])])


def beside(value1, value2):
    'The values occur side-by-side: blends1 & fox2 | blends2 & fox1 ...'
    return from_dnf([(comb(value1, i), comb(value2, j))
                     for i, j in zip(houses, houses[1:])] +
                    [(comb(value2, i), comb(value1, j))
                     for i, j in zip(houses, houses[1:])]
                    )


cnf = []

# each house gets exactly one value from each attribute group
for house in houses:
    for group in groups:
        cnf += one_of(comb(value, house) for value in group)

# each value gets assigned to exactly one house
for value in values:
    cnf += one_of(comb(value, house) for house in houses)


# The Englishman lives in the Red house.
cnf += same_house('Englishman', 'Red')
cnf += same_house('Spaniard', 'Dog')  # The Spaniard owns the Dog.
cnf += same_house('Coffee', 'Green')  # Coffee is drunk in the Green house.
cnf += same_house('Ukrainian', 'Tea')  # The Ukrainian drinks Tea.
# The Green house is immediately to the right of the Ivory house.
cnf += consecutive('Ivory', 'Green')
cnf += same_house('Old Gold', 'Snails')  # The Old Gold smoker owns Snails.
cnf += same_house('Kools', 'Yellow')  # Kools are smoked in the Yellow house.
cnf += found_at('Milk', 3)  # Milk is drunk in the middle house.
cnf += found_at('Norwegian', 1)  # The Norwegian lives in the first house.
# The man who smokes Chesterfields lives in the house next to the man with the Fox.
cnf += beside('Chesterfield', 'Fox')
# Kools are smoked in the house next to the house where the Horse is kept.
cnf += beside('Kools', 'Horse')
# The Lucky Strike smoker drinks Orange juice.
cnf += same_house('Lucky Strike', 'Orange juice')
# The Japanese smokes Parliaments.
cnf += same_house('Japanese', 'Parliament')
# The Norwegian lives next to the Blue house.
cnf += beside('Norwegian', 'Blue')

pprint(solve_one(cnf))

['Yellow 1',  'Norwegian 1',  'Water 1',  'Kools 1',  'Fox 1',  'Blue 2',  'Ukrainian 2',  'Tea 2',  'Chesterfield 2',  'Horse 2',  'Red 3', 'Englishman 3',  'Milk 3',  'Old Gold 3',  'Snails 3',  'Ivory 4',  'Spaniard 4',  'Orange juice 4',  'Lucky Strike 4',  'Dog 4',  'Green 5',  'Japanese 5',  'Coffee 5',  'Parliament 5',  'Zebra 5']

in what looks like machine code form.

Table[
  WolframAlpha[
   "penrose tiles " <> n <> " steps", {{"Result", 1}, 
    "Content"}], {n, {"4", "5"}}]

Entity["NonperiodicTiling", "KitesAndDartsTiling"]["Dataset"]

ResourceData["Demonstrations Project: Making Patterns with Wang Tiles"]

Exact Cover and Boolean Satisfiability problems are proven to be NP-complete, right?

patterns = {
   {1, 2, 3, 4},
   {4, 3, 2, 1},
   {2, 1, 4, 3},
   {3, 4, 1, 2}
   };
init = {1, 2, 3, 4};
size = 5;
tiles = GenerateWangTiling[patterns, init, size];
Show@MapIndexed[
  Graphics@WangTile[
     Reverse[#2 {-1, 1}],
     #1] &,
  tiles,
  {2}]

I don't know how else to put it, that GenerateWangTiling handles the Wang tiling problems in such a way that we can find a deterministic form of the SAT, Boolean Satisfiability Problem. And we could rewrite all of these functions here, except then we wouldn't just be able to get it over with in the beginning...so all of those Jeandel-Rao functions for the eleven tiles, I didn't know they were that significant in aperiodic tiling theory. I kind of wish we would change the generation of these 3x3 arrangements, the validity of these tiles, yet we changed nothing.

patterns = {
   {1, 2, 3, 4},
   {4, 3, 2, 1},
   {2, 1, 4, 3},
   {3, 4, 1, 2}
   };

init = {1, 2, 3, 4};
sizes = Range[5, 50, 5];
times = Table[
  AbsoluteTiming[
    GenerateWangTiling[
      patterns,
      init,
      size];][[1]],
  {size, sizes}]
ListLinePlot[times,
 PlotLabels -> sizes,
 AxesLabel -> {"Size", "Time"}]

I think the biggest thing for us is how we got timestamps for the Wang tiling benchmarks which just shows how efficient the template-based approach can be, what with all these von Neumann neighborhood templates that we use, derived from the Jeandel-Rao tiles, that we compare with the original Wang tiling...yes, I have the results from the SAT solver and that's how we ride off into the sunset..we have this MultiwayDeletionsGraph approach that significantly outperforms the SAT solver in terms of speed. I mean sure, the SAT solver demonstrates reasonable efficiency for Wang tiling which makes it a viable option for such combinatorial problems. The trick here is that when compared with MDG, this result gets inverted; the latter shows superior performance, particularly when handling complex tiling problems efficiently.

Graphics[WangTile[{0, 0}, #], ImageSize -> 50] & /@ patterns

While there are a lot of ways that we can generate computational overhead, the GenerateWangTiling function exists. From the periodic boundary conditions to a variety of tiling problems, then there is nothing other than the versatility and robust nature of the GenerateWangTiling function which highlights how flexible these tiling problems can be, either through direct SAT solving or through constructing and analyzing templates, which is really quite idiomatic in its insights but we can see that the SAT solver versus MDG, both methods provided valuable insights but differed in complexity and execution time.

The role that Jeandel-Rao tiles play in aperiodic tiling is the hallmark of the pitfalls of the complexity of the nuances that are involved in such tiling problems. The SAT solver could be effectively utilized and adapted accordingly to study the aperiodic tilings. We could do a walkthrough of the impact crater of these complex combinatorial problems like the Queen's board game of Isolation, wherein the insights cultivate a beneficial environment for future research and applications in tiling theory and combinatorial optimization.

tileCountPeriodic = Length@GenerateWangTiling[
   patterns, {}, {3, 3}, All, "Boundary" -> "Periodic"]
tileCountAny = Length@GenerateWangTiling[
   patterns, {}, {3, 3}, All, "Boundary" -> "Any"]
Print["Set of 3x3 tilings with periodic boundaries cardinality, is ", \
tileCountPeriodic]
Print["Set of 3x3 tilings with any boundaries cardinality, is ", \
tileCountAny]

I think this is essentially defining how we can display tiling using Show@MapIndexed, in the sense of starting out with a 5x5 tiling pattern where each Wang tile is colored with red, green, blue, or gray, and arranged in such a way that the edges of the adjacent tiles match in color. It's not like we're doing this in some single-state problem where the agent knows the exact state. I can't quite put my finger on how we get certain custom and abstract functions like GenerateWangTiling and WangTile. They just aren't what I'm used to when I think of visual representations of the SAT (Boolean Satisfiability Problem) solver applied to Wang tiling. Because if there's one thing to note, it's that Wang tiles are a type of mathematical tiling that use square tiles with colored edges, it's our way of consolidating the rules; there are so many ways to connect a node and draw things that we normally think of as just adjacent tiles that have the same color on their touching edges.

AbsoluteTiming[
Show@MapIndexed[
Graphics@WangTile[
Reverse[#2 {-1, 1}], #1] &,
GenerateWangTiling[JeandelRaoTiles, {}, 20],
{2}]]

You'd think that we're demonstrating the use of Mathematica-specific syntax to solve complex concepts like Wang tiles, but it's actually the other way around. With the conditional reduction that we explore in fields related to computational geometry, tiling theory, and combinatorics, we can find how some spatial tiling patterns can generalize, generate and visualize the Wang tilings...at the intersection of mathematics and computer science, where again, as a best practice, we utilize the AbsoluteTiming function to measure how long each tiling generation takes and in plain language, it takes a ridiculous amount of computational complexity and scalability to reconstruct what happened when we pressed these keys and derived meaningful insights atomically, I can't take my eyes off the & /@ syntax that vectorizes maps the WangTile function to each element in a list, for what it's worth. It's the difference between an introduction and a question that allows us to see that period boundary conditions show how the tiling pattern repeats itself at the boundaries, and that, is how the output 2 is trying to say that there are two such tilings, really. It's easier to train the humans than it is to have the human user interface adapt and that's the kind of stuff that we talk about, which periodic boundary has which constraint and it's something that we're a little scared to admit or acknowledge.