jix · November 20, 2022 01:21
diff --git a/enum_pentasquare.sage b/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)
	#!/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)