diff --git a/papers/medial_tire_cuts/experiments/medial_tire_dual_cut_experiment.py b/papers/medial_tire_cuts/experiments/medial_tire_dual_cut_experiment.py new file mode 100644 index 0000000..d2f6901 --- /dev/null +++ b/papers/medial_tire_cuts/experiments/medial_tire_dual_cut_experiment.py @@ -0,0 +1,344 @@ +"""Source-dual cut from a chained medial tire cut. + +Companion to ``run_medial_tire_cut_experiment.py``. Where that script reports +the cut graph of M(G), this one takes the same chained walk-depth labelling and +cut and reads it off as a *source-dual* cut: the planar dual of the source +triangulation G with the cut edges removed, as drawn for +``seed59_min5_dual_cut_1.png``. + +The dual of a plane triangulation G has one node per triangular face and one +edge per primal edge (joining the two faces that share it). Its faces are the +*vertices* of G, each bounded by ``deg(v)`` dual edges. A medial tire cut at an +annular medial vertex removes the dual edge of the corresponding primal edge; +the interesting quantity is how many of those removed (``missing``) dual edges +surround each dual face (vertex of G). For ``seed59`` at source 5 the maximum +is 3, around the degree-9 vertex 3. + +Four chained entry points (broad to narrow control): + + * ``random_dual_cut(n, ...)`` -- find a random maximal planar graph of a given + minimum degree, then defer to ``dual_cut_random_source``. + * ``dual_cut_random_source(G, ...)`` -- choose a random level source, then + defer to ``dual_cut_random_entry``. + * ``dual_cut_random_entry(G, source, ...)`` -- choose a random root entry + tooth, then defer to ``medial_tire_dual_cut``. + * ``medial_tire_dual_cut(G, source, entry_edge)`` -- the worker: chain the + walk-depth labelling/cut from the given root entry tooth and assemble the + source-dual cut. + +Run with the repo venv (networkx; matplotlib only for ``--png``): +``.venv/bin/python``. +""" + +from __future__ import annotations + +import argparse +import os +import random +import sys +from collections import defaultdict + +import networkx as nx + +_HERE = os.path.dirname(os.path.abspath(__file__)) +_MTD = os.path.normpath(os.path.join( + _HERE, "..", "..", + "medial_tire_decompositions_of_plane_triangulations", "experiments")) +sys.path.insert(0, _MTD) +sys.path.insert(0, _HERE) + +from tire_realization_analysis import ( # noqa: E402 + ekey, extract_tread, medial_graph, medial_tire_facemodel, + recognise, triangular_faces, +) +from run_medial_tire_cut_experiment import ( # noqa: E402 + _assemble_cut_graph, _cap_cut, _label_treads, + random_maximal_planar_min_degree, +) + + +# --------------------------------------------------------------------------- # +# Tread recognition and the source-dual graph. +# --------------------------------------------------------------------------- # + +def _build_treads(faces, levels): + """Recognise the full medial tire graph of every BFS-level tread. + + Returns ``(treads, skipped)`` where ``treads`` maps depth ``d`` to the + recognised ``(g, bij)`` and ``skipped`` lists ``(d, reason)`` for the rest. + """ + treads, skipped = {}, [] + for d in range(max(levels.values())): + tread = extract_tread(faces, levels, d) + if tread is None: + skipped.append((d, "no tread faces")) + continue + if len(tread["up"]) < 3: + skipped.append((d, f"only {len(tread['up'])} up teeth")) + continue + rec = recognise(medial_tire_facemodel(tread["tread_faces"]), tread) + if rec is None: + skipped.append((d, "not a valid full medial tire graph")) + continue + treads[d] = rec + return treads, skipped + + +def root_entry_choices(G, source): + """Edge indices of the root tread's up teeth -- the eligible entry teeth. + + Empty when ``source`` induces no recognised root tread. + """ + faces, _ = triangular_faces(G) + levels = nx.single_source_shortest_path_length(G, source) + treads, _ = _build_treads(faces, levels) + if not treads: + return [] + g, _bij = treads[min(treads)] + return sorted(g.up_edges) + + +def source_dual(G, faces): + """The planar dual of triangulation ``G``: one node per face, one edge per + primal edge (tagged ``primal``). Faces are indexed as in ``faces``.""" + edge_faces = defaultdict(list) + for fi, f in enumerate(faces): + for a, b in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])): + edge_faces[ekey(a, b)].append(fi) + D = nx.Graph() + D.add_nodes_from(range(len(faces))) + for e, fs in edge_faces.items(): + if len(fs) == 2: + D.add_edge(fs[0], fs[1], primal=e) + return D + + +def removed_dual_edges(results, cap_cuts): + """The set of primal edges whose dual edge a cut removes (cap cut plus every + tread cut).""" + removed = set() + for c in cap_cuts or []: + removed.add(c["medial_vertex"]) + for d in sorted(results): + g, bij = results[d]["g"], results[d]["bij"] + for c in results[d]["cuts"]: + if c.vertex is not None: + removed.add(bij[f"a{c.vertex}"]) + return removed + + +def dual_face_missing(G, removed): + """For each dual face (vertex ``v`` of ``G``), the number of bounding dual + edges removed by the cut.""" + return {v: sum(1 for w in G.neighbors(v) if ekey(v, w) in removed) + for v in G.nodes()} + + +# --------------------------------------------------------------------------- # +# The four chained entry points. +# --------------------------------------------------------------------------- # + +def medial_tire_dual_cut(G, source, entry_edge): + """Chain the walk-depth labelling/cut from root entry tooth ``entry_edge`` + and assemble the source-dual cut of ``G`` at level ``source``. + + ``entry_edge`` must be an up tooth of the root (lowest recognised) tread; + see ``root_entry_choices``. Returns a structured result dict. + """ + faces, emb = triangular_faces(G) + M = medial_graph(G) + levels = nx.single_source_shortest_path_length(G, source) + treads, skipped = _build_treads(faces, levels) + if not treads: + raise ValueError(f"level source {source} induces no recognised tread") + + g_root = treads[min(treads)][0] + if entry_edge not in g_root.up_edges: + raise ValueError( + f"entry edge {entry_edge} is not an up tooth of the root tread " + f"(choices: {sorted(g_root.up_edges)})") + + results = {} + _label_treads(treads, results, root_entry_edge=entry_edge) + cap_cuts = _cap_cut(G, emb, source, levels, results) + cut_graph, labels, warnings = _assemble_cut_graph(M, results, cap_cuts=cap_cuts) + + dual = source_dual(G, faces) + removed = removed_dual_edges(results, cap_cuts) + missing = dual_face_missing(G, removed) + + return { + "G": G, "M": M, "source": source, "entry_edge": entry_edge, + "faces": faces, "outer_face": 0, + "levels": levels, "treads": treads, "skipped": skipped, + "results": results, "cap_cuts": cap_cuts, "cut_graph": cut_graph, + "labels": labels, "warnings": warnings, + "dual": dual, "removed_dual_edges": removed, "dual_face_missing": missing, + "max_missing": max(missing.values()) if missing else 0, + } + + +def dual_cut_random_entry(G, source, rng=None): + """Pick a random root entry tooth at ``source``, then ``medial_tire_dual_cut``.""" + rng = rng or random.Random() + choices = root_entry_choices(G, source) + if not choices: + raise ValueError(f"level source {source} induces no recognised root tread") + return medial_tire_dual_cut(G, source, rng.choice(choices)) + + +def dual_cut_random_source(G, rng=None): + """Pick a random level source, then ``dual_cut_random_entry``. + + Sources are tried in random order; the first one inducing a recognised root + tread is used (a maximal planar graph always has at least one).""" + rng = rng or random.Random() + sources = sorted(G.nodes()) + rng.shuffle(sources) + for source in sources: + if root_entry_choices(G, source): + return dual_cut_random_entry(G, source, rng=rng) + raise ValueError("no level source induces a recognised root tread") + + +def random_dual_cut(n=20, seed=0, rng=None, min_degree=5, flips=400, attempts=1000): + """Find a random maximal planar graph of minimum degree ``min_degree``, then + ``dual_cut_random_source``. + + ``seed`` drives the graph sample; ``rng`` (defaulting to ``Random(seed)``) + drives the random source and entry choices, so the whole pipeline is + reproducible from ``(n, seed)``. + """ + rng = rng or random.Random(seed) + G, graph_seed = random_maximal_planar_min_degree( + n, seed, flips=flips, min_degree=min_degree, attempts=attempts) + result = dual_cut_random_source(G, rng=rng) + result["graph_seed"] = graph_seed + result["min_degree"] = min(dict(G.degree()).values()) + return result + + +# --------------------------------------------------------------------------- # +# Reporting and (optional) rendering. +# --------------------------------------------------------------------------- # + +def summary(result): + G, missing = result["G"], result["dual_face_missing"] + removed = result["removed_dual_edges"] + hist = defaultdict(int) + for k in missing.values(): + hist[k] += 1 + lines = [ + f"source-dual cut: n={G.number_of_nodes()} " + f"graph_seed={result.get('graph_seed', '?')} " + f"min_degree={result.get('min_degree', min(dict(G.degree()).values()))}", + f"level source: vertex {result['source']} " + f"root entry tooth: e{result['entry_edge']}", + f"recognised treads: {sorted(result['treads'])} " + f"skipped: {result['skipped']}", + f"removed source-dual edges ({len(removed)}): " + f"{sorted(removed)}", + f"dual-face missing-edge histogram (count by #removed around the dual " + f"face): {dict(sorted(hist.items()))} max={result['max_missing']}", + ] + for v in sorted(missing, key=lambda v: (-missing[v], v)): + if missing[v]: + inc = [ekey(v, w) for w in G.neighbors(v) if ekey(v, w) in removed] + lines.append(f" dual face v{v} (deg {G.degree(v)}): " + f"{missing[v]} missing -> {inc}") + return "\n".join(lines) + + +def draw_png(result, path, scale=6.0): + """Render the source-dual cut: dual nodes at face centroids, dual edges + drawn light gray where the cut removed them, labelled by missing count.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + from draw_medial_tire_cut import _source_layout # local import; needs numpy + + G, faces, dual = result["G"], result["faces"], result["dual"] + removed = result["removed_dual_edges"] + missing = result["dual_face_missing"] + pos_v = _source_layout(G) + + def centroid(fi): + xs = [pos_v[u][0] for u in faces[fi]] + ys = [pos_v[u][1] for u in faces[fi]] + return (sum(xs) / 3.0, sum(ys) / 3.0) + + pos = {fi: centroid(fi) for fi in dual.nodes()} + fig, ax = plt.subplots(figsize=(7, 7)) + # primal graph, faint, for orientation + for u, v in G.edges(): + ax.plot([pos_v[u][0], pos_v[v][0]], [pos_v[u][1], pos_v[v][1]], + color="0.85", lw=0.5, zorder=0) + for u, v, data in dual.edges(data=True): + cut = data["primal"] in removed + ax.plot([pos[u][0], pos[v][0]], [pos[u][1], pos[v][1]], + color="0.80" if cut else "0.25", + lw=1.0 if cut else 1.3, + linestyle=(0, (2, 2)) if cut else "solid", zorder=1) + for fi in dual.nodes(): + x, y = pos[fi] + # label each dual face's source vertex by its missing count instead: + ax.plot(x, y, "o", ms=4, color="#3a6ea5", zorder=2) + # annotate dual faces (vertices of G) with their missing count + for v in G.nodes(): + m = missing[v] + if m: + x, y = pos_v[v] + ax.text(x, y, str(m), color="#b03030", fontsize=8, + ha="center", va="center", zorder=3, + bbox=dict(boxstyle="circle,pad=0.1", fc="white", + ec="#b03030", lw=0.6)) + ax.set_title(f"source-dual cut (source {result['source']}, entry " + f"e{result['entry_edge']}); gray = edges missing after cuts\n" + f"red numbers = #missing dual edges around each dual face; " + f"max {result['max_missing']}", fontsize=9) + ax.set_aspect("equal") + ax.axis("off") + fig.tight_layout() + fig.savefig(path, dpi=150) + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("-n", type=int, default=20, help="number of vertices") + parser.add_argument("--seed", type=int, default=0, help="graph sample seed") + parser.add_argument("--min-degree", type=int, default=5) + parser.add_argument("--source", type=int, default=None, + help="fix the level source (default: random via rng)") + parser.add_argument("--entry", type=int, default=None, + help="fix the root entry tooth (requires --source)") + parser.add_argument("--png", metavar="PATH", help="render the dual cut to PNG") + args = parser.parse_args() + + rng = random.Random(args.seed) + if args.source is not None and args.entry is not None: + G, graph_seed = random_maximal_planar_min_degree( + args.n, args.seed, min_degree=args.min_degree) + result = medial_tire_dual_cut(G, args.source, args.entry) + result["graph_seed"] = graph_seed + result["min_degree"] = min(dict(G.degree()).values()) + elif args.source is not None: + G, graph_seed = random_maximal_planar_min_degree( + args.n, args.seed, min_degree=args.min_degree) + result = dual_cut_random_entry(G, args.source, rng=rng) + result["graph_seed"] = graph_seed + result["min_degree"] = min(dict(G.degree()).values()) + else: + result = random_dual_cut(n=args.n, seed=args.seed, + rng=rng, min_degree=args.min_degree) + + print(summary(result)) + if args.png: + draw_png(result, args.png) + print(f"wrote {args.png}") + + +if __name__ == "__main__": + main() diff --git a/papers/medial_tire_cuts/experiments/run_medial_tire_cut_experiment.py b/papers/medial_tire_cuts/experiments/run_medial_tire_cut_experiment.py index 8b28477..75e5e49 100644 --- a/papers/medial_tire_cuts/experiments/run_medial_tire_cut_experiment.py +++ b/papers/medial_tire_cuts/experiments/run_medial_tire_cut_experiment.py @@ -68,14 +68,23 @@ def _apex_vertex(g, bij, edge): return bij[g.apex_of_edge(edge)] -def _label_treads(treads, results): +def _label_treads(treads, results, root_entry_edge=None): """Fill ``results[d]`` with the walk-depth labelling and cuts for each - recognised tread ``d``, chaining child entries to parent down teeth.""" + recognised tread ``d``, chaining child entries to parent down teeth. + + The root tread (lowest recognised depth) is entered at ``root_entry_edge`` + when given -- it must be one of that tread's up teeth -- otherwise at an + arbitrary up tooth. + """ + root_d = min(treads) if treads else None for d in sorted(treads): g, bij = treads[d] parent = treads.get(d - 1) if parent is None: - entry_edge, start_depth = g.up_edges[0], 0 # arbitrary root entry + if d == root_d and root_entry_edge is not None: + entry_edge, start_depth = root_entry_edge, 0 + else: + entry_edge, start_depth = g.up_edges[0], 0 # arbitrary root entry else: pg, pbij = parent pdepth = results[d - 1]["depth"] diff --git a/papers/medial_tire_cuts/funcD_seed7_dual_cut.png b/papers/medial_tire_cuts/funcD_seed7_dual_cut.png new file mode 100644 index 0000000..faa519d Binary files /dev/null and b/papers/medial_tire_cuts/funcD_seed7_dual_cut.png differ