enum_pentasquare.sage

#!/usr/bin/env sage
import itertools

# cycle decomposition of a permutation
def cycles(perm):
    seen = [False] * len(perm)
    for i in range(len(perm)):
        if seen[i]:
            continue
        cycle = [i]
        seen[i] = True
        while not seen[i := perm[i]]:
            seen[i] = True
            cycle.append(i)
        yield cycle

for n in (5..8):
    print("===", n, "===", flush=True)

    counter = 0

    # The tiling we're looking for has fewer faces than vertices so it's faster to
    # enumerate candidates for the dual graph. There's a catch though, because in our
    # tiling the corners of the square can be degree-2 vertices which turn into parallel
    # edges of the dual graph. The enumeration doesn't produce parallel edges, but as
    # the square corners are the only place where we can have a degree 2 vertex, it's
    # faster to ignore them for the initial enumeration and to insert them manually
    # after selecting the outside face.
    #
    # Unfortunately this does mean we can only restrict the minimum degree to 5 - c
    # where c is the maximal number of square corners in one pentagon. This makes the
    # planar graph enumeration quite a bit slower but it's still fast enough.
    for dual_g in graphs.planar_graphs(n + 1, minimum_degree=2, minimum_connectivity=2):
        hist = dual_g.degree_histogram() + [0, 0, 0]

        # dual_g vertices with degree < 5 are either a) the outside face of our tiling
        # or b) pentagons that use square corners we still need to add. Let's call them
        # bad dual_g vertices
        d2 = hist[2]
        d3 = hist[3]
        d4 = hist[4]

        # count missing number of half edges
        missing = d4 + 2 * d3 + 3 * d2

        # we can add at most 4 dual_g edges from the dual_g vertex for the outside face
        # to up to 4 corners:
        if missing > 8 or (d2 + d3 + d4) > 5:
            continue

        # find the set of bad dual_g vertices
        bad = set(v for v in dual_g if dual_g.degree(v) < 5)

        # The closed neighborhood of the dual_g vertex corresponding to the tiling's
        # outside face must contain all bad dual_g vertices, but isn't necessarily a bad
        # vertex itself. This finds all candidates:
        outside_candidates = set(bad) if len(bad) > 4 else None

        for b in bad:
            bn = set(dual_g.neighbors(b))
            bn.add(b)
            if outside_candidates is None:
                outside_candidates = bn
            else:
                outside_candidates &= bn

        # With no bad vertices, any vertex could correspond to the outside face. This
        # never happens for small enough values of n. As is, this is would most likely
        # be quite slow and would really need some way to only make non-isomorphic
        # choices.
        if not bad:
            outside_candidates = set(dual_g)

        # If no candidates are left, dual_g cannot be the dual graph of the kind of
        # tiling we're looking for.
        if not outside_candidates:
            continue

        # Go through the candidates, most of the time there is only one
        for outside in outside_candidates:
            corners = []

            # The bad vertices besides outside correspond to the pentagons that need
            # square corners
            for needs_corners in bad:
                if needs_corners == outside:
                    continue

                # Add as many corners as required to make it at least a pentagon
                corners.extend([needs_corners] * (5 - dual_g.degree(needs_corners)))

            # If we needed more than 4 corners, that's not a valid tiling
            if len(corners) > 4:
                continue

            # If we have corners left, we can add them to any pentagon touching the
            # outside, but we also need to consider the case where we don't add them, as
            # split corners shared by two or more pentagons are already considered by
            # the enumartion of dual_g.
            for extra_corner_count in (0..(4 - len(corners))):
                for extra_corners in itertools.combinations(
                            dual_g.neighbors(outside), extra_corner_count):

                    # To add the parallel edges to dual_g we first convert it from the
                    # `get_embedding()` representation which doesn't handle parallel
                    # edges, to a combinatorial embedding of dual_g's dual, i.e. the
                    # pentagon tiling of the square.
                    #
                    # See also https://en.wikipedia.org/wiki/Combinatorial_map
                    #
                    # The edge-flip involution is fixed and exchanges adjacent odd and
                    # even pairs.
                    idmap = {}
                    face_ccw = [None] * dual_g.num_edges() * 2
                    face_cw = [None] * dual_g.num_edges() * 2
                    for v, ns in dual_g.get_embedding().items():
                        for n in ns:
                            if (v, n) not in idmap:
                                idmap[v,n] = len(idmap)
                                idmap[n,v] = len(idmap)

                        for n, nn in zip(ns, ns[1:] + ns[:1]):
                            face_ccw[idmap[v, n]] = idmap[v, nn]
                            face_cw[idmap[v, nn]] = idmap[v, n]

                    # Then we can add the square corners one by one, with some linked
                    # list manipulation
                    for rb in [*corners, *extra_corners]:
                        new = len(face_ccw)

                        face_ccw.extend([None] * 2)
                        face_cw.extend([None] * 2)

                        orig_r_next = face_ccw[idmap[outside, rb]]
                        face_ccw[idmap[outside, rb]] = new
                        face_cw[new] = idmap[outside, rb]
                        face_ccw[new] = orig_r_next
                        face_cw[orig_r_next] = new

                        orig_rb_next = face_cw[idmap[rb, outside]]
                        face_cw[idmap[rb, outside]] = new + 1
                        face_ccw[new + 1] = idmap[rb, outside]
                        face_cw[new + 1] = orig_rb_next
                        face_ccw[orig_rb_next] = new + 1

                    # From the face links and we can compute the vertex links
                    vertex_ccw = [face_cw[i ^^ 1] for i in range(len(face_cw))]

                    # And from those, we can build a graph
                    edges = {}
                    tiling = Graph()
                    for i, half_edges in enumerate(cycles(vertex_ccw)):
                        tiling.add_vertex(i)
                        for half_edge in half_edges:
                            edges.setdefault(half_edge // 2, []).append(i)

                    for edge in edges.values():
                        tiling.add_edge(*edge)

                    counter += 1

                    # There's quite a few things we haven't checked yet. We would need
                    # to make sure that the combinatorial embedding corresponds to a
                    # geometric embedding where exactly 5 vertices of every pentagon is
                    # a corner with <180deg and every additional vertex is on a side.
                    # But we're lucky and this excludes enough candidates to manually
                    # check this.
                    print(f"==> Found candidate #{counter} with {n} pentagons!")

                    # We don't actually use our embedding for plotting this, so we plot
                    # a few random layouts in the hope that one obviously looks like the
                    # tiling we're looking for.
                    for plot_variation in range(4):
                        fname = f"penta_{n}_{plot_variation}.png"
                        print(f"  writing {fname!r}")
                        plot(tiling).save_image(fname)