math-research/papers/level_switching/experiments/stress_test_termination.py

"""Stress-test the preprocessing algorithm on random maximal-outerplanar
graphs of varying size. Track whether the algorithm always reaches
all-depth-0 and how many steps it takes."""
import random
import sys
import networkx as nx
from collections import Counter


def face_edges(f):
    return {frozenset((f[0], f[1])), frozenset((f[1], f[2])),
            frozenset((f[0], f[2]))}


def random_triangulation(num_vertices, seed=None):
    """Random triangulation of a convex polygon on `num_vertices` vertices.
    Returns (outer_edges, chords, faces).

    Uses a recursive splitting algorithm. The polygon vertices are
    labelled 0..n-1 around the outer cycle."""
    rng = random.Random(seed)
    n = num_vertices
    outer = [(i, (i + 1) % n) for i in range(n)]
    chords = []
    faces = []

    def split(start, end):
        """Triangulate the sub-polygon from `start` to `end` (inclusive),
        going through vertices start, start+1, ..., end on the outer cycle.
        The "chord" closing this sub-polygon is (start, end)."""
        if (end - start) % n <= 1:
            return
        if (end - start) % n == 2:
            mid = (start + 1) % n
            faces.append(tuple(sorted((start, mid, end))))
            return
        # Pick a random vertex strictly between start and end as the apex
        length = (end - start) % n
        offset = rng.randint(1, length - 1)
        mid = (start + offset) % n
        faces.append(tuple(sorted((start, mid, end))))
        if mid != (start + 1) % n:
            chords.append(tuple(sorted((start, mid))))
            split(start, mid)
        if end != (mid + 1) % n:
            chords.append(tuple(sorted((mid, end))))
            split(mid, end)

    # Start by picking an outer-cycle vertex as the "root" for splitting;
    # we triangulate the polygon as if the chord 0--(n-1) were closing
    # the sub-polygon on the other side. To get a proper triangulation
    # of the polygon, we use vertex 0 as the splitting root.
    split(0, n - 1)
    # Note: this doesn't add a chord from 0 to n-1 since that's an
    # outer edge.

    return outer, chords, faces


def compute_depths(faces, outer_set):
    D = nx.Graph()
    D.add_nodes_from(range(len(faces)))
    for i, fi in enumerate(faces):
        for j, fj in enumerate(faces):
            if i < j and face_edges(fi) & face_edges(fj):
                D.add_edge(i, j)
    B = [i for i, f in enumerate(faces)
         if len(face_edges(f) & outer_set) >= 1]
    if not B:
        return {i: float('inf') for i in range(len(faces))}
    return {i: min(nx.shortest_path_length(D, i, b) for b in B)
            for i in range(len(faces))}


def check_balanced(F_idx, faces, depth_, outer_set):
    F = faces[F_idx]
    fe = face_edges(F)
    for e in fe:
        if e in outer_set:
            continue
        cands = [j for j in range(len(faces))
                 if j != F_idx and e in face_edges(faces[j])]
        if not cands:
            continue
        Fp_idx = cands[0]
        if depth_[Fp_idx] != depth_[F_idx] - 1:
            continue
        Fp = faces[Fp_idx]
        d = depth_[F_idx]
        ok = True
        for e2 in face_edges(Fp):
            if e2 == e:
                continue
            if e2 in outer_set:
                continue
            others = [j for j in range(len(faces))
                      if j != Fp_idx and e2 in face_edges(faces[j])]
            if not others or depth_[others[0]] != d - 2:
                ok = False
                break
        if ok:
            return True, F_idx, Fp_idx, e
    return False, None, None, None


def apply_switch(faces, uv, wx):
    u, v = uv
    w, x = wx
    new_faces = [f for f in faces
                 if set(f) != {u, v, w} and set(f) != {u, v, x}]
    new_faces.append(tuple(sorted((u, w, x))))
    new_faces.append(tuple(sorted((v, w, x))))
    return new_faces


def run_algorithm(faces, outer_set, max_steps=2000, seed=None):
    rng = random.Random(seed)
    balanced_steps = 0
    preprocess_steps = 0
    depth_history = []
    for step in range(max_steps):
        depth = compute_depths(faces, outer_set)
        d_max = max(depth.values())
        depth_history.append(d_max)
        if d_max == 0:
            return {'terminated': True, 'steps': step,
                    'balanced': balanced_steps,
                    'preprocess': preprocess_steps,
                    'depth_history': depth_history}
        # Pick any maximum-depth face
        max_d_faces = [i for i, d in depth.items() if d == d_max]
        F_idx = rng.choice(max_d_faces)
        F = faces[F_idx]

        ok, _, fp_idx, e = check_balanced(F_idx, faces, depth, outer_set)
        if ok:
            Fp = faces[fp_idx]
            u, v = tuple(e)
            w = [vert for vert in F if vert not in (u, v)][0]
            x = [vert for vert in Fp if vert not in (u, v)][0]
            faces = apply_switch(faces, (u, v), (w, x))
            balanced_steps += 1
            continue

        # Preprocess: pick a random depth-(d-1) neighbour
        choices = []
        for e_test in [frozenset((F[0], F[1])), frozenset((F[1], F[2])),
                       frozenset((F[0], F[2]))]:
            if e_test in outer_set:
                continue
            cands = [j for j in range(len(faces))
                     if j != F_idx and e_test in face_edges(faces[j])]
            if cands and depth[cands[0]] == d_max - 1:
                choices.append((e_test, cands[0]))
        if not choices:
            return {'terminated': False, 'reason': 'no depth-(d-1) neighbour',
                    'depth_history': depth_history,
                    'final_depth': d_max}
        e, fp_idx = rng.choice(choices)
        Fp = faces[fp_idx]
        u, v = tuple(e)
        w = [vert for vert in F if vert not in (u, v)][0]
        x = [vert for vert in Fp if vert not in (u, v)][0]
        faces = apply_switch(faces, (u, v), (w, x))
        preprocess_steps += 1

    return {'terminated': False, 'reason': 'max_steps reached',
            'depth_history': depth_history,
            'final_depth': max(compute_depths(faces, outer_set).values())}


if __name__ == '__main__':
    results = []
    sizes = [10, 14, 18, 24, 30, 40, 60, 80]
    trials_per_size = 20

    for n in sizes:
        n_terminated = 0
        max_steps = 0
        max_init_depth = 0
        failed = []
        for trial in range(trials_per_size):
            seed = hash((n, trial)) & 0xffffffff
            outer, chords, faces = random_triangulation(n, seed=seed)
            outer_set = {frozenset(e) for e in outer}
            init_depth = max(compute_depths(faces, outer_set).values())
            max_init_depth = max(max_init_depth, init_depth)
            result = run_algorithm(faces, outer_set,
                                   max_steps=10 * len(faces),
                                   seed=seed + 1)
            if result['terminated']:
                n_terminated += 1
                max_steps = max(max_steps, result['steps'])
            else:
                failed.append((seed, result))
        print(f'n={n}: {n_terminated}/{trials_per_size} terminated; '
              f'max init depth {max_init_depth}, max steps {max_steps}')
        if failed:
            for seed, res in failed[:3]:
                print(f'  FAIL seed={seed}: {res}')
            results.append((n, failed))

    if not results:
        print('\nAll trials terminated successfully.')
    else:
        print(f'\n{sum(len(f) for _, f in results)} failures across {len(results)} sizes')