libcompose.py

"""
High-level composition API for Truchet tiling.

Provides compose_* functions that place groups of tiles onto an existing Plane,
and make_* convenience wrappers that create a fresh Plane. Each returns a
Composition with the plane and tile outlines (for rendering).
"""

from __future__ import annotations

import os
import random as _random_module
from dataclasses import dataclass, field
from enum import Enum
from math import sqrt, hypot
from time import perf_counter

from z3 import Int, If, Solver, Optimize, sat, Or, And, Sum

from libgeom import (Point, rect, hexagon_in_rect,
                     oct_square_corners, corner_tri_vertices,
                     octagon_in_square)
from libtile import (TileType, edge_point_world, edge_points_with_normals,
                     oct_tile_type,
                     corner_tri_tile_type, enriched_square_tile_type,
                     hex_tile_type, tri30_tile_type,
                     equilateral_tri_tile_type)
from libtopo import Plane, ConstraintKind, Tile


class UnsatisfiableMatchingError(Exception):
    """Raised when solve_matching_assignment cannot find a valid matching assignment.

    This can occur when:
    1. No compatible matchings exist for a shared edge between tiles
    2. The Z3 solver determines the constraints are unsatisfiable

    Attributes:
        message: Human-readable error description
        tile_i: Index of first tile (for edge incompatibility)
        tile_j: Index of second tile (for edge incompatibility)
        edge_i: Edge index on tile_i (for edge incompatibility)
        edge_j: Edge index on tile_j (for edge incompatibility)
        positions: Full list of tile positions (for unsat case)
    """
    def __init__(self, message: str, *,
                 tile_i: int | None = None,
                 tile_j: int | None = None,
                 edge_i: int | None = None,
                 edge_j: int | None = None,
                 positions: list[TilePosition] | None = None):
        super().__init__(message)
        self.message = message
        self.tile_i = tile_i
        self.tile_j = tile_j
        self.edge_i = edge_i
        self.edge_j = edge_j
        self.positions = positions


@dataclass
class Composition:
    plane: Plane
    tile_outlines: list[list[Point]] = field(default_factory=list[list[Point]])


# ---------------------------------------------------------------------------
# compose_* functions — place tiles onto an existing Plane
# ---------------------------------------------------------------------------

def compose_square(plane: Plane, x: float, y: float, s: float,
                   tile_type: TileType | None = None) -> Composition:
    """Place a single square tile at (x, y) with side s."""
    if tile_type is None:
        tile_type = TileType.uniform("sq", 4, [0.5])
    verts = list(rect(x, y, s, s))
    plane.place_tile(tile_type, verts)
    return Composition(plane=plane, tile_outlines=[verts])


def compose_square_grid(plane: Plane, x: float, y: float, s: float,
                        rows: int, cols: int,
                        tile_type: TileType | None = None) -> Composition:
    """Place a rows x cols grid of square tiles."""
    if tile_type is None:
        tile_type = TileType.uniform("sq", 4, [0.5])
    outlines: list[list[Point]] = []
    for r in range(rows):
        for c in range(cols):
            verts = list(rect(x + c * s, y + r * s, s, s))
            plane.place_tile(tile_type, verts)
            outlines.append(verts)
    return Composition(plane=plane, tile_outlines=outlines)


def compose_488(plane: Plane, x: float, y: float, s: float) -> Composition:
    """Place a full 4.8.8 cell: octagon + 4 corner triangles + 4 enriched squares."""
    sq_corners, oct_v = oct_square_corners(x, y, s)
    tri_verts = corner_tri_vertices(sq_corners, oct_v)

    plane.place_tile(oct_tile_type(), oct_v)
    for tv in tri_verts:
        plane.place_tile(corner_tri_tile_type(), tv)

    sq_rects = [
        list(rect(x, y - s, s, s)),
        list(rect(x + s, y, s, s)),
        list(rect(x, y + s, s, s)),
        list(rect(x - s, y, s, s)),
    ]
    for sq_edge, sq_v in zip([2, 3, 0, 1], sq_rects):
        plane.place_tile(enriched_square_tile_type(sq_edge), sq_v)

    outlines: list[list[Point]] = [list(oct_v)] + [list(tv) for tv in tri_verts] + sq_rects
    return Composition(plane=plane, tile_outlines=outlines)


def compose_hex_corners(plane: Plane, x: float, y: float, L: float) -> Composition:
    """Place a hexagon + 4 corner 30-60-90 triangles in a 2L x L*sqrt(3) rect."""
    w = 2 * L
    h = L * sqrt(3)
    hex_v = hexagon_in_rect(x, y, L)
    r_tl: Point = (x, y)
    r_tr: Point = (x + w, y)
    r_br: Point = (x + w, y + h)
    r_bl: Point = (x, y + h)

    plane.place_tile(hex_tile_type(), hex_v)
    tri_verts_list = [
        [r_tl, hex_v[0], hex_v[5]],
        [r_tr, hex_v[2], hex_v[1]],
        [r_br, hex_v[3], hex_v[2]],
        [r_bl, hex_v[5], hex_v[4]],
    ]
    for tv in tri_verts_list:
        plane.place_tile(tri30_tile_type(), tv)

    outlines: list[list[Point]] = [list(hex_v)] + tri_verts_list
    return Composition(plane=plane, tile_outlines=outlines)


def compose_hex_stacked(plane: Plane, x: float, y: float, L: float,
                        n: int = 2) -> Composition:
    """Place n vertically stacked hexagons."""
    h = L * sqrt(3)
    outlines: list[list[Point]] = []
    for i in range(n):
        hex_v = hexagon_in_rect(x, y + i * h, L)
        plane.place_tile(hex_tile_type(), hex_v)
        outlines.append(list(hex_v))
    return Composition(plane=plane, tile_outlines=outlines)


def compose_equilateral_cell(plane: Plane, x: float, y: float,
                             L: float) -> Composition:
    """Place an equilateral triangle cell: 4 equilateral tris + 4 corner tris."""
    w = 2 * L
    h = L * sqrt(3)
    sl = L / 2
    ll = L * sqrt(3) / 2

    rc = list(rect(x, y, w, h))
    tl, tr, br, bl = rc

    t1 = [(x + sl, y), (x, y + ll), (x + L, y + ll)]
    t2 = [(x + sl, y), (x + L, y + ll), (x + w - sl, y)]
    t3 = [(x + sl, y + h), (x, y + ll), (x + L, y + ll)]
    t4 = [(x + w - sl, y + h), (x + sl, y + h), (x + L, y + ll)]

    teq = equilateral_tri_tile_type()
    t30 = tri30_tile_type()

    outlines: list[list[Point]] = []
    for tv in [t1, t2, t3, t4]:
        plane.place_tile(teq, tv)
        outlines.append(tv)

    corner_tris = [
        [tl, (x + sl, y), (x, y + ll)],
        [tr, (x + w, y + ll), (x + w - sl, y)],
        [br, (x + w - sl, y + h), (x + w, y + ll)],
        [bl, (x, y + ll), (x + sl, y + h)],
    ]
    for tv in corner_tris:
        plane.place_tile(t30, tv)
        outlines.append(tv)

    return Composition(plane=plane, tile_outlines=outlines)


def compose_hex_eq_grid(plane: Plane, x: float, y: float, L: float,
                        rows: int, cols: int,
                        hex_cols: int = 2) -> Composition:
    """Place a grid mixing hex cells and equilateral triangle cells.

    First hex_cols columns use hex+corners, remaining columns use equilateral cells.
    """
    w = 2 * L
    h = L * sqrt(3)
    all_outlines: list[list[Point]] = []

    for row in range(rows):
        for col in range(cols):
            gx = x + col * w
            gy = y + row * h
            if col < hex_cols:
                comp = compose_hex_corners(plane, gx, gy, L)
            else:
                comp = compose_equilateral_cell(plane, gx, gy, L)
            all_outlines.extend(comp.tile_outlines)

    return Composition(plane=plane, tile_outlines=all_outlines)


# ---------------------------------------------------------------------------
# Non-crossing perfect matchings
# ---------------------------------------------------------------------------

def noncrossing_matchings(n: int) -> list[list[tuple[int, int]]]:
    """All non-crossing perfect matchings of n points on a circle (0..n-1).

    n must be even. Returns list of matchings, each a sorted list of pairs.
    Counted by Catalan number C(n/2).
    """
    if n % 2 != 0:
        raise ValueError(f"n must be even, got {n}")
    if n == 0:
        return [[]]

    def _enum(pts: list[int]) -> list[list[tuple[int, int]]]:
        if len(pts) == 0:
            return [[]]
        if len(pts) == 2:
            return [[(pts[0], pts[1])]]
        result: list[list[tuple[int, int]]] = []
        first = pts[0]
        # first matches with pts[k] for odd k (1, 3, 5, ...)
        for k in range(1, len(pts), 2):
            partner = pts[k]
            left = pts[1:k]       # between first and partner
            right = pts[k + 1:]   # after partner
            for m_left in _enum(left):
                for m_right in _enum(right):
                    result.append(sorted([(first, partner)] + m_left + m_right))
        return result

    return _enum(list(range(n)))


def _catalan_table(n: int) -> list[int]:
    """Precompute Catalan numbers C(0)..C(n)."""
    c = [0] * (n + 1)
    c[0] = 1
    for i in range(1, n + 1):
        c[i] = c[i - 1] * 2 * (2 * i - 1) // (i + 1)
    return c


def sample_random_ncpm(n: int, rng: _random_module.Random,
                       min_span: int = 0) -> list[tuple[int, int]]:
    """Sample one random NCPM of n points (0..n-1) uniformly at random. O(n).

    Uses Catalan-weighted recursive decomposition: pick a partner for the
    first point in each sub-region with probability proportional to the
    product of Catalan numbers for the resulting sub-problems.

    min_span: if > 0, reject candidates whose circular distance to the
    first point is less than min_span (eliminates pips when min_span=2).
    """
    if n % 2 != 0:
        raise ValueError(f"n must be even, got {n}")
    if n == 0:
        return []

    cat = _catalan_table(n // 2)
    result: list[tuple[int, int]] = []

    def _sample(pts: list[int]):
        if len(pts) == 0:
            return
        if len(pts) == 2:
            result.append((pts[0], pts[1]))
            return
        first = pts[0]
        m = len(pts)
        # Candidates: pts[k] for odd k (1, 3, 5, ...)
        weights: list[int] = []
        candidates: list[int] = []
        for k in range(1, m, 2):
            left_size = k - 1       # points between first and partner
            right_size = m - k - 1  # points after partner
            w = cat[left_size // 2] * cat[right_size // 2]
            # Skip too-close candidates (circular distance on original ring)
            if min_span > 0 and w > 0:
                dist = pts[k] - first
                circ_dist = min(dist, n - dist)
                if circ_dist < min_span:
                    w = 0
            weights.append(w)
            candidates.append(k)

        # Weighted random choice
        total = sum(weights)
        if total == 0:
            raise UnsatisfiableMatchingError(
                f"No valid partner for point {first} with min_span={min_span}")
        r = rng.randrange(total)
        cumulative = 0
        chosen_k = candidates[-1]
        for k_idx, w in zip(candidates, weights):
            if w == 0:
                continue
            cumulative += w
            if r < cumulative:
                chosen_k = k_idx
                break

        partner = pts[chosen_k]
        result.append((first, partner))
        _sample(pts[1:chosen_k])
        _sample(pts[chosen_k + 1:])

    _sample(list(range(n)))
    return sorted(result)


def sample_constrained_ncpm(
    n: int,
    forced_pairs: list[tuple[int, int]],
    rng: _random_module.Random,
    min_span: int = 0,
) -> list[tuple[int, int]]:
    """Sample one random NCPM of n points with some pairs forced. O(n).

    forced_pairs: pairs that must appear in the matching (from neighbor
    edge constraints). These must be non-crossing among themselves.
    min_span: if > 0, skip free candidates whose circular distance to
    the first point is less than min_span.

    Raises UnsatisfiableMatchingError if the constraints can't be satisfied.
    """
    if n % 2 != 0:
        raise ValueError(f"n must be even, got {n}")
    if n == 0:
        return []
    if not forced_pairs:
        return sample_random_ncpm(n, rng, min_span=min_span)

    cat = _catalan_table(n // 2)

    # Build forced partner lookup
    forced_partner: dict[int, int] = {}
    for a, b in forced_pairs:
        forced_partner[a] = b
        forced_partner[b] = a

    result: list[tuple[int, int]] = []

    def _sample(pts: list[int]):
        if len(pts) == 0:
            return
        if len(pts) == 2:
            result.append((min(pts[0], pts[1]), max(pts[0], pts[1])))
            return
        if len(pts) % 2 != 0:
            raise UnsatisfiableMatchingError(
                f"Odd number of points ({len(pts)}) in sub-region")

        first = pts[0]
        m = len(pts)

        # If first has a forced partner, pair them directly
        if first in forced_partner:
            partner = forced_partner[first]
            k = pts.index(partner)
            if k % 2 == 0:
                raise UnsatisfiableMatchingError(
                    f"Forced pair ({first}, {partner}) encloses odd points")
            result.append((min(first, partner), max(first, partner)))
            _sample(pts[1:k])
            _sample(pts[k + 1:])
            return

        # First is free — pick a random partner from valid odd-indexed candidates
        weights: list[int] = []
        candidates: list[int] = []
        for k in range(1, m, 2):
            candidate = pts[k]

            # If candidate is forced to someone other than first, skip —
            # candidate is reserved for its forced partner
            if candidate in forced_partner:
                weights.append(0)
                candidates.append(k)
                continue

            left = pts[1:k]
            right = pts[k + 1:]
            left_set = set(left)
            right_set = set(right)

            # Check no forced pair is split across left/right
            valid = True
            for p in left:
                if p in forced_partner and forced_partner[p] not in left_set:
                    valid = False
                    break
            if valid:
                for p in right:
                    if p in forced_partner and forced_partner[p] not in right_set:
                        valid = False
                        break
            if not valid:
                weights.append(0)
                candidates.append(k)
                continue

            # Count free points in each sub-region
            left_free = sum(1 for p in left if p not in forced_partner)
            right_free = sum(1 for p in right if p not in forced_partner)

            if left_free % 2 != 0 or right_free % 2 != 0:
                weights.append(0)
                candidates.append(k)
                continue

            w = cat[left_free // 2] * cat[right_free // 2]
            # Skip too-close free candidates (circular distance on original ring)
            if min_span > 0 and w > 0:
                dist = candidate - first
                circ_dist = min(dist, n - dist)
                if circ_dist < min_span:
                    w = 0
            weights.append(w)
            candidates.append(k)

        total = sum(weights)
        if total == 0:
            raise UnsatisfiableMatchingError(
                f"No valid partner for point {first} in sub-region of {m} points")

        r = rng.randrange(total)
        cumulative = 0
        chosen_k = candidates[-1]
        for k_idx, w in zip(candidates, weights):
            if w == 0:
                continue
            cumulative += w
            if r < cumulative:
                chosen_k = k_idx
                break

        partner = pts[chosen_k]
        result.append((min(first, partner), max(first, partner)))
        _sample(pts[1:chosen_k])
        _sample(pts[chosen_k + 1:])

    _sample(list(range(n)))
    return sorted(result)


def tile_edge_points(tile_type: TileType, vertices: list[Point]
                     ) -> list[Point]:
    """Compute world-coordinate positions of all edge points for a tile."""
    pts: list[Point] = []
    for edge_idx, tvals in enumerate(tile_type.edge_points):
        for t in tvals:
            pts.append(edge_point_world(vertices, edge_idx, t))
    return pts


def tile_matchings(tile_type: TileType, vertices: list[Point]
                   ) -> list[list[tuple[Point, Point]]]:
    """For a tile, return all NCPMs as lists of world-coordinate endpoint pairs."""
    n = sum(len(ep) for ep in tile_type.edge_points)
    matchings = noncrossing_matchings(n)
    world_pts = tile_edge_points(tile_type, vertices)
    result: list[list[tuple[Point, Point]]] = []
    for matching in matchings:
        pairs = [(world_pts[a], world_pts[b]) for a, b in matching]
        result.append(pairs)
    return result


# ---------------------------------------------------------------------------
# make_* convenience functions — create a fresh Plane
# ---------------------------------------------------------------------------

def make_square(x: float = 0, y: float = 0, s: float = 1,
                tile_type: TileType | None = None) -> Composition:
    return compose_square(Plane(), x, y, s, tile_type)


def make_square_grid(x: float = 0, y: float = 0, s: float = 1,
                     rows: int = 2, cols: int = 2,
                     tile_type: TileType | None = None) -> Composition:
    return compose_square_grid(Plane(), x, y, s, rows, cols, tile_type)


def make_488(x: float = 0, y: float = 0, s: float = 1) -> Composition:
    return compose_488(Plane(), x, y, s)


def make_hex_corners(x: float = 0, y: float = 0,
                     L: float = 1) -> Composition:
    return compose_hex_corners(Plane(), x, y, L)


def make_hex_eq_grid(x: float = 0, y: float = 0, L: float = 1,
                     rows: int = 2, cols: int = 3,
                     hex_cols: int = 2) -> Composition:
    return compose_hex_eq_grid(Plane(), x, y, L, rows, cols, hex_cols)


# ---------------------------------------------------------------------------
# Edge signature functions for SAT-based matching assignment
# ---------------------------------------------------------------------------

def edge_local_pairing(matching: list[tuple[int, int]], tile_type: TileType
                       ) -> tuple[frozenset[tuple[int, int]], ...]:
    """Per-edge local pairing for a matching.

    Returns N frozensets, one per edge. Each frozenset contains
    (local_i, local_j) pairs of edge points paired within that edge.
    Points not in any local pair connect to other edges.
    """
    offsets: list[int] = []
    offset = 0
    for ep in tile_type.edge_points:
        offsets.append(offset)
        offset += len(ep)

    result: list[frozenset[tuple[int, int]]] = []
    for e in range(tile_type.n_edges):
        k = len(tile_type.edge_points[e])
        start = offsets[e]
        edge_indices = set(range(start, start + k))
        pairs: set[tuple[int, int]] = set()
        for a, b in matching:
            if a in edge_indices and b in edge_indices:
                la, lb = a - start, b - start
                pairs.add((min(la, lb), max(la, lb)))
        result.append(frozenset(pairs))
    return tuple(result)


def reversed_edge_pairing(pairing: frozenset[tuple[int, int]], k: int
                          ) -> frozenset[tuple[int, int]]:
    """Reverse a local pairing for edge traversal in opposite direction.
    k is the number of edge points on this edge."""
    reversed_pairs: set[tuple[int, int]] = set()
    for a, b in pairing:
        ra = k - 1 - a
        rb = k - 1 - b
        reversed_pairs.add((min(ra, rb), max(ra, rb)))
    return frozenset(reversed_pairs)


def matching_edge_pairings(tile_type: TileType
                           ) -> list[tuple[frozenset[tuple[int, int]], ...]]:
    """Precompute edge pairings for all NCPMs of a tile type."""
    n = sum(len(ep) for ep in tile_type.edge_points)
    matchings = noncrossing_matchings(n)
    return [edge_local_pairing(m, tile_type) for m in matchings]


EdgeSig = frozenset[tuple[int, int]]


@dataclass
class SignatureData:
    """Precomputed per-edge signature groupings for a tile type."""
    n_matchings: int
    # edge_groups[edge_idx] maps each unique signature -> list of matching indices with that signature
    edge_groups: list[dict[EdgeSig, list[int]]]
    # full_sig_groups maps (sig_e0, sig_e1, ...) -> list of matching indices
    full_sig_groups: dict[tuple[EdgeSig, ...], list[int]]


def compute_signature_data(tile_type: TileType) -> SignatureData:
    """Group matchings by their per-edge signature."""
    n = sum(len(ep) for ep in tile_type.edge_points)
    matchings = noncrossing_matchings(n)
    all_pairings = [edge_local_pairing(m, tile_type) for m in matchings]

    edge_groups: list[dict[EdgeSig, list[int]]] = []
    for edge_idx in range(tile_type.n_edges):
        groups: dict[EdgeSig, list[int]] = {}
        for m_idx, pairings in enumerate(all_pairings):
            sig = pairings[edge_idx]
            groups.setdefault(sig, []).append(m_idx)
        edge_groups.append(groups)

    full_sig_groups: dict[tuple[EdgeSig, ...], list[int]] = {}
    for m_idx, pairings in enumerate(all_pairings):
        full_sig_groups.setdefault(pairings, []).append(m_idx)

    return SignatureData(len(matchings), edge_groups, full_sig_groups)


class SolveStrategy(Enum):
    BASELINE = "baseline"   # Current: enumerate all matching pairs
    GROUPED = "grouped"     # Group by signature, no new Z3 vars
    TWOPHASE = "twophase"   # Phase 1: Z3 picks edge sigs. Phase 2: random matching


@dataclass
class TilePosition:
    tile_type: TileType
    vertices: list[Point]
    neighbors: list[tuple[int, int, int]] = field(default_factory=list[tuple[int, int, int]])
    # neighbors: list of (neighbor_idx, my_edge, their_edge)


def find_adjacency(positions: list[TilePosition], tol: float = 1e-6) -> None:
    """Find shared edges between tile positions by vertex matching.

    Populates each position's neighbors list in-place.
    """
    def quantize(p: Point) -> tuple[int, int]:
        return (round(p[0] / tol), round(p[1] / tol))

    # Map directed edge (va, vb) to (tile_idx, edge_idx)
    edge_map: dict[tuple[tuple[int, int], tuple[int, int]],
                   tuple[int, int]] = {}

    for i, pos in enumerate(positions):
        n = len(pos.vertices)
        for e in range(n):
            va = quantize(pos.vertices[e])
            vb = quantize(pos.vertices[(e + 1) % n])
            forward = (va, vb)
            reverse = (vb, va)
            if reverse in edge_map:
                j, ej = edge_map[reverse]
                pos.neighbors.append((j, e, ej))
                positions[j].neighbors.append((i, ej, e))
                del edge_map[reverse]
            else:
                edge_map[forward] = (i, e)


def validate_matching_assignment(positions: list[TilePosition],
                                  assignment: list[int]) -> None:
    """Verify all shared edges have compatible pairings.

    Raises AssertionError if any shared edge has incompatible pairings.
    """
    type_cache: dict[int, list[tuple[frozenset[tuple[int, int]], ...]]] = {}
    for pos in positions:
        key = id(pos.tile_type)
        if key not in type_cache:
            type_cache[key] = matching_edge_pairings(pos.tile_type)

    for i, pos in enumerate(positions):
        for j, my_edge, their_edge in pos.neighbors:
            if j <= i:
                continue
            pairings_i = type_cache[id(pos.tile_type)]
            pairings_j = type_cache[id(positions[j].tile_type)]
            k = len(pos.tile_type.edge_points[my_edge])
            sig_i = pairings_i[assignment[i]][my_edge]
            sig_j = pairings_j[assignment[j]][their_edge]
            assert sig_i == reversed_edge_pairing(sig_j, k), (
                f"Edge incompatibility: tile {i} edge {my_edge} "
                f"vs tile {j} edge {their_edge}")


def _add_baseline_constraints(
    s: Solver, variables: list[Int], positions: list[TilePosition],
    type_cache: dict[int, list[tuple[frozenset[tuple[int, int]], ...]]],
) -> None:
    """BASELINE strategy: enumerate all matching pairs per shared edge."""
    processed: set[frozenset[tuple[int, int, int]]] = set()

    for i, pos in enumerate(positions):
        for neighbor_idx, my_edge, their_edge in pos.neighbors:
            edge_key = frozenset({(i, my_edge), (neighbor_idx, their_edge)})
            if edge_key in processed:
                continue
            processed.add(edge_key)

            j = neighbor_idx
            pairings_i = type_cache[id(pos.tile_type)]
            pairings_j = type_cache[id(positions[j].tile_type)]
            k = len(pos.tile_type.edge_points[my_edge])
            k_j = len(positions[j].tile_type.edge_points[their_edge])

            if k != k_j:
                raise UnsatisfiableMatchingError(
                    f"Incompatible edge point counts: tile {i} edge {my_edge} has {k} points, "
                    f"but tile {j} edge {their_edge} has {k_j} points. "
                    f"Tile {i} type: {pos.tile_type.name}, "
                    f"Tile {j} type: {positions[j].tile_type.name}. "
                    f"Shared edges must have the same number of edge points.",
                    tile_i=i, tile_j=j, edge_i=my_edge, edge_j=their_edge)

            compatible: list = []
            for mi, pi in enumerate(pairings_i):
                for mj, pj in enumerate(pairings_j):
                    sig_i = pi[my_edge]
                    sig_j = pj[their_edge]
                    if sig_i == reversed_edge_pairing(sig_j, k):
                        compatible.append(And(variables[i] == mi,
                                              variables[j] == mj))

            if not compatible:
                raise UnsatisfiableMatchingError(
                    f"No compatible matchings found between tile {i} edge {my_edge} "
                    f"and tile {j} edge {their_edge}. "
                    f"Tile {i} type: {pos.tile_type.name}, "
                    f"Tile {j} type: {positions[j].tile_type.name}. "
                    f"Edge has {k} edge points but none of the {len(pairings_i)} x {len(pairings_j)} "
                    f"matching combinations produce compatible edge pairings.",
                    tile_i=i, tile_j=j, edge_i=my_edge, edge_j=their_edge)
            s.add(Or(*compatible))


def _add_grouped_constraints(
    s: Solver, variables: list[Int], positions: list[TilePosition],
    sig_cache: dict[int, SignatureData],
) -> None:
    """GROUPED strategy: group matchings by signature, O(1) lookup per signature."""
    processed: set[frozenset[tuple[int, int, int]]] = set()

    for i, pos in enumerate(positions):
        for neighbor_idx, my_edge, their_edge in pos.neighbors:
            edge_key = frozenset({(i, my_edge), (neighbor_idx, their_edge)})
            if edge_key in processed:
                continue
            processed.add(edge_key)

            j = neighbor_idx
            k = len(pos.tile_type.edge_points[my_edge])
            k_j = len(positions[j].tile_type.edge_points[their_edge])

            if k != k_j:
                raise UnsatisfiableMatchingError(
                    f"Incompatible edge point counts: tile {i} edge {my_edge} has {k} points, "
                    f"but tile {j} edge {their_edge} has {k_j} points. "
                    f"Tile {i} type: {pos.tile_type.name}, "
                    f"Tile {j} type: {positions[j].tile_type.name}. "
                    f"Shared edges must have the same number of edge points.",
                    tile_i=i, tile_j=j, edge_i=my_edge, edge_j=their_edge)

            groups_i = sig_cache[id(pos.tile_type)].edge_groups[my_edge]
            groups_j = sig_cache[id(positions[j].tile_type)].edge_groups[their_edge]

            clauses = []
            for sig_i, members_i in groups_i.items():
                sig_j_needed = reversed_edge_pairing(sig_i, k)
                if sig_j_needed in groups_j:
                    members_j = groups_j[sig_j_needed]
                    clause_i = Or([variables[i] == m for m in members_i])
                    clause_j = Or([variables[j] == m for m in members_j])
                    clauses.append(And(clause_i, clause_j))

            if not clauses:
                raise UnsatisfiableMatchingError(
                    f"No compatible matchings found between tile {i} edge {my_edge} "
                    f"and tile {j} edge {their_edge}. "
                    f"Tile {i} type: {pos.tile_type.name}, "
                    f"Tile {j} type: {positions[j].tile_type.name}. "
                    f"Edge has {k} edge points but no compatible signature groups.",
                    tile_i=i, tile_j=j, edge_i=my_edge, edge_j=their_edge)
            s.add(Or(*clauses))


def _solve_twophase(
    positions: list[TilePosition],
    sig_cache: dict[int, SignatureData],
    diverse: bool = True,
) -> list[int]:
    """Two-phase solver: Z3 picks edge signatures, Python picks matchings.

    Phase 1: One Z3 variable per shared edge choosing a compatible
    signature pair (domain ~6-20). Per-tile constraint ensures the
    combination of chosen edge sigs is jointly realizable.

    Phase 2: For each tile, randomly pick a matching consistent with
    the chosen signatures on all shared edges.
    """
    import random

    # Enumerate shared edges and their compatible signature pairs
    shared_edges: list[tuple[int, int, int, int]] = []  # (i, my_edge, j, their_edge)
    processed: set[frozenset[tuple[int, int]]] = set()

    for i, pos in enumerate(positions):
        for j, my_edge, their_edge in pos.neighbors:
            edge_key = frozenset({(i, my_edge), (j, their_edge)})
            if edge_key in processed:
                continue
            processed.add(edge_key)
            shared_edges.append((i, my_edge, j, their_edge))

    # For each shared edge, find compatible (sig_i, sig_j) pairs
    # and number the sigs for each (tile_type, edge) combination
    type_edge_sig_list: dict[tuple[int, int], list[EdgeSig]] = {}

    def get_sig_idx(tile_type_id: int, edge: int, sig: EdgeSig) -> int:
        key = (tile_type_id, edge)
        if key not in type_edge_sig_list:
            sd = sig_cache[tile_type_id]
            type_edge_sig_list[key] = list(sd.edge_groups[edge].keys())
        return type_edge_sig_list[key].index(sig)

    edge_compat: list[list[tuple[int, int]]] = []  # per shared edge: list of (sig_idx_i, sig_idx_j)
    for i, my_edge, j, their_edge in shared_edges:
        ti_id = id(positions[i].tile_type)
        tj_id = id(positions[j].tile_type)
        groups_i = sig_cache[ti_id].edge_groups[my_edge]
        groups_j = sig_cache[tj_id].edge_groups[their_edge]
        k = len(positions[i].tile_type.edge_points[my_edge])

        pairs = []
        for sig_i in groups_i:
            sig_j_needed = reversed_edge_pairing(sig_i, k)
            if sig_j_needed in groups_j:
                idx_i = get_sig_idx(ti_id, my_edge, sig_i)
                idx_j = get_sig_idx(tj_id, their_edge, sig_j_needed)
                pairs.append((idx_i, idx_j))
        edge_compat.append(pairs)

    # Build per-tile shared edge info: which shared_edge indices touch this tile,
    # and which side (i or j) is this tile on
    tile_shared: list[list[tuple[int, int, int]]] = [[] for _ in positions]
    # (shared_edge_idx, which_edge_of_tile, 0=i_side/1=j_side)
    for se_idx, (i, my_edge, j, their_edge) in enumerate(shared_edges):
        tile_shared[i].append((se_idx, my_edge, 0))
        tile_shared[j].append((se_idx, their_edge, 1))

    # Phase 1: Z3 picks one compatible pair per shared edge
    s = Solver()
    s.set("threads", os.cpu_count() or 1)
    edge_vars = [Int(f"se_{idx}") for idx in range(len(shared_edges))]
    for idx, pairs in enumerate(edge_compat):
        s.add(Or([edge_vars[idx] == p_idx for p_idx in range(len(pairs))]))

    # Per-tile realizability: the combination of chosen sigs must correspond
    # to at least one actual matching
    for tile_idx, pos in enumerate(positions):
        if not tile_shared[tile_idx]:
            continue

        tt_id = id(pos.tile_type)
        sd = sig_cache[tt_id]

        # Get the edge indices that are shared, and for each, the sig list
        shared_info = tile_shared[tile_idx]  # (se_idx, edge_of_tile, side)

        # Build a map: for each full sig tuple in this tile type,
        # determine which edge_var values would produce it
        # full_sig_groups maps (sig_e0, sig_e1, ...) -> matching indices
        # We need: for each valid full sig tuple, the conjunction of edge_var constraints

        # First, for shared edges, map: edge_of_tile -> (se_idx, side)
        edge_constraints: dict[int, tuple[int, int]] = {}
        for se_idx, edge_of_tile, side in shared_info:
            edge_constraints[edge_of_tile] = (se_idx, side)

        # For each full sig tuple that has matchings, check if it's consistent
        # with some edge_var assignment
        valid_clauses = []
        for full_sig, _matching_ids in sd.full_sig_groups.items():
            clause_parts = []
            feasible = True
            for edge_of_tile, (se_idx, side) in edge_constraints.items():
                target_sig = full_sig[edge_of_tile]
                target_idx = get_sig_idx(tt_id, edge_of_tile, target_sig)
                # Find which pair indices in edge_compat[se_idx] have this sig on the right side
                matching_pair_indices = []
                for p_idx, (pi, pj) in enumerate(edge_compat[se_idx]):
                    if (side == 0 and pi == target_idx) or (side == 1 and pj == target_idx):
                        matching_pair_indices.append(p_idx)
                if not matching_pair_indices:
                    feasible = False
                    break
                clause_parts.append(Or([edge_vars[se_idx] == p for p in matching_pair_indices]))
            if feasible:
                valid_clauses.append(And(clause_parts) if clause_parts else True)

        if valid_clauses:
            s.add(Or(valid_clauses))

    n_se = len(shared_edges)
    n_compat = sum(len(pairs) for pairs in edge_compat)
    print(f"  [twophase] {n_se} shared edges, {n_compat} compatible pairs, "
          f"{len(sig_cache)} tile types")

    t0 = perf_counter()
    if s.check() != sat:
        raise UnsatisfiableMatchingError("Two-phase: edge signature assignment unsatisfiable")
    t_phase1 = perf_counter() - t0

    model = s.model()

    # Decode chosen signatures per shared edge
    chosen_sigs: dict[tuple[int, int], EdgeSig] = {}  # (tile_idx, edge) -> sig
    for se_idx, (i, my_edge, j, their_edge) in enumerate(shared_edges):
        pair_idx = model[edge_vars[se_idx]].as_long()
        sig_idx_i, sig_idx_j = edge_compat[se_idx][pair_idx]
        ti_id = id(positions[i].tile_type)
        tj_id = id(positions[j].tile_type)
        chosen_sigs[(i, my_edge)] = type_edge_sig_list[(ti_id, my_edge)][sig_idx_i]
        chosen_sigs[(j, their_edge)] = type_edge_sig_list[(tj_id, their_edge)][sig_idx_j]

    # Phase 2: pick a random matching per tile consistent with chosen sigs
    result = []
    for tile_idx, pos in enumerate(positions):
        tt_id = id(pos.tile_type)
        sd = sig_cache[tt_id]

        # Filter full_sig_groups to those matching chosen sigs on shared edges
        candidates: list[int] = []
        for full_sig, matching_ids in sd.full_sig_groups.items():
            ok = True
            for edge_of_tile in range(pos.tile_type.n_edges):
                if (tile_idx, edge_of_tile) in chosen_sigs:
                    if full_sig[edge_of_tile] != chosen_sigs[(tile_idx, edge_of_tile)]:
                        ok = False
                        break
            if ok:
                candidates.extend(matching_ids)

        if not candidates:
            raise UnsatisfiableMatchingError(
                f"Two-phase: no matching for tile {tile_idx} with chosen edge sigs")
        result.append(random.choice(candidates))

    print(f"  [twophase] phase1={t_phase1*1000:.1f}ms")
    return result


# ---------------------------------------------------------------------------
# Bezier crossing detection + swap helpers
# ---------------------------------------------------------------------------

def _bezier_polyline(p1: Point, n1: Point, p2: Point, n2: Point,
                     n_segments: int = 16) -> list[Point]:
    """Flatten a cubic bezier into a polyline of n_segments+1 points.

    Control points computed inline: c = p + 0.4 * dist * n.
    """
    dist = hypot(p2[0] - p1[0], p2[1] - p1[1])
    alpha = 0.4 * dist
    c1 = (p1[0] + alpha * n1[0], p1[1] + alpha * n1[1])
    c2 = (p2[0] + alpha * n2[0], p2[1] + alpha * n2[1])
    pts: list[Point] = []
    for i in range(n_segments + 1):
        t = i / n_segments
        u = 1 - t
        x = (u*u*u * p1[0] + 3*u*u*t * c1[0] +
             3*u*t*t * c2[0] + t*t*t * p2[0])
        y = (u*u*u * p1[1] + 3*u*u*t * c1[1] +
             3*u*t*t * c2[1] + t*t*t * p2[1])
        pts.append((x, y))
    return pts


def _segments_intersect(p1: Point, p2: Point, p3: Point, p4: Point) -> bool:
    """Strict segment-segment intersection test (cross-product orientation).

    Returns False for shared endpoints and collinear cases.
    """
    def cross(o: Point, a: Point, b: Point) -> float:
        return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

    d1 = cross(p3, p4, p1)
    d2 = cross(p3, p4, p2)
    d3 = cross(p1, p2, p3)
    d4 = cross(p1, p2, p4)

    if ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and \
       ((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0)):
        return True
    return False


def _beziers_cross(poly1: list[Point], poly2: list[Point]) -> bool:
    """Check if two bezier polylines intersect (O(n^2) segment check)."""
    for i in range(len(poly1) - 1):
        for j in range(len(poly2) - 1):
            if _segments_intersect(poly1[i], poly1[i+1],
                                   poly2[j], poly2[j+1]):
                return True
    return False


def _chords_cross(a: int, b: int, c: int, d: int, n: int) -> bool:
    """Circle topology test: chord (a,b) crosses chord (c,d) on n points.

    Returns True iff exactly one of {c,d} is strictly inside arc a..b
    (going from a to b in increasing order mod n).
    """
    # Normalize so a < b
    if a > b:
        a, b = b, a
    if c > d:
        c, d = d, c
    # c is inside arc a..b iff a < c < b
    c_inside = a < c < b
    d_inside = a < d < b
    return c_inside != d_inside


def _try_chord_swap(
    matching: list[tuple[int, int]],
    forced_set: set[tuple[int, int]],
    idx_a: int, idx_b: int, n: int,
) -> list[tuple[int, int]] | None:
    """Try swapping two crossing pairs and return new matching if valid NCPM.

    Given pairs matching[idx_a]=(a1,a2) and matching[idx_b]=(b1,b2),
    try both alternative pairings:
      option1: (a1,b1)+(a2,b2)
      option2: (a1,b2)+(a2,b1)
    Validate NCPM (no chord crossings with any other pair) and that
    forced pairs are not disturbed.
    """
    a1, a2 = matching[idx_a]
    b1, b2 = matching[idx_b]

    # Don't touch forced pairs
    if (min(a1, a2), max(a1, a2)) in forced_set:
        return None
    if (min(b1, b2), max(b1, b2)) in forced_set:
        return None

    others = [matching[i] for i in range(len(matching))
              if i != idx_a and i != idx_b]

    for new_p1, new_p2 in [
        ((min(a1, b1), max(a1, b1)), (min(a2, b2), max(a2, b2))),
        ((min(a1, b2), max(a1, b2)), (min(a2, b1), max(a2, b1))),
    ]:
        # Check new pairs don't cross each other
        if _chords_cross(new_p1[0], new_p1[1], new_p2[0], new_p2[1], n):
            continue
        # Check new pairs don't cross any existing pair
        ok = True
        for op in others:
            if _chords_cross(new_p1[0], new_p1[1], op[0], op[1], n):
                ok = False
                break
            if _chords_cross(new_p2[0], new_p2[1], op[0], op[1], n):
                ok = False
                break
        if ok:
            result = others + [new_p1, new_p2]
            result.sort()
            return result
    return None


def _refine_no_bezier_crossings(
    matching: list[tuple[int, int]],
    tile_type: TileType,
    vertices: list[Point],
    forced_pairs: list[tuple[int, int]],
    max_passes: int = 5,
) -> list[tuple[int, int]]:
    """Iteratively detect and fix bezier crossings by swapping pairs.

    Returns the (possibly improved) matching.
    """
    forced_set = {(min(a, b), max(a, b)) for a, b in forced_pairs}
    total_eps = sum(len(ep) for ep in tile_type.edge_points)

    for _pass in range(max_passes):
        ep_data = edge_points_with_normals(tile_type, vertices)
        polylines = []
        for a, b in matching:
            p1, n1 = ep_data[a]
            p2, n2 = ep_data[b]
            polylines.append(_bezier_polyline(p1, n1, p2, n2))

        # Find first crossing
        found = False
        for i in range(len(matching)):
            for j in range(i + 1, len(matching)):
                if _beziers_cross(polylines[i], polylines[j]):
                    new_matching = _try_chord_swap(
                        matching, forced_set, i, j, total_eps)
                    if new_matching is not None:
                        matching = new_matching
                        found = True
                        break
            if found:
                break

        if not found:
            break  # No crossings or no fixable crossings

    return matching


def solve_matching_greedy(
    positions: list[TilePosition],
    seed: int | None = None,
    max_retries: int = 100,
    diversity_samples: int = 20,
    min_span: int = 0,
    avoid_crossings: bool = False,
) -> list[list[tuple[int, int]]]:
    """Greedy tile-by-tile matching solver. No Z3 — pure Python.

    Uses O(n) per-tile random NCPM sampling with neighbor constraints.
    No precomputation of all matchings — works for any EP count.

    diversity_samples: sample this many candidates per tile and pick the
    one whose edge signature differs most from already-assigned neighbors.

    min_span: reject candidates containing any pair with span < min_span.
    Use min_span=2 to eliminate pips (adjacent-point pairs).

    Returns list of matching pairs per tile (not indices).
    """
    rng = _random_module.Random(seed)

    n = len(positions)

    def _try_once(rng: _random_module.Random
                  ) -> tuple[list[list[tuple[int, int]]],
                             list[list[tuple[int, int]]]] | None:
        assignments: list[list[tuple[int, int]] | None] = [None] * n
        edge_pairings_cache: list[tuple[frozenset[tuple[int, int]], ...] | None] = [None] * n
        forced_per_tile: list[list[tuple[int, int]]] = [[] for _ in range(n)]

        for i in range(n):
            pos = positions[i]
            tt = pos.tile_type
            total_eps = sum(len(ep) for ep in tt.edge_points)

            # Collect forced pairs from already-assigned neighbors
            forced_pairs: list[tuple[int, int]] = []
            assigned_neighbor_sigs: list[tuple[frozenset[tuple[int, int]], ...]] = []
            for j, my_edge, their_edge in pos.neighbors:
                if assignments[j] is None:
                    continue
                # Get neighbor's edge pairing on the shared edge
                neighbor_sig = edge_pairings_cache[j][their_edge]
                k = len(tt.edge_points[my_edge])
                needed_sig = reversed_edge_pairing(neighbor_sig, k)
                # Convert local edge signature to global EP index pairs
                offset = sum(len(tt.edge_points[e]) for e in range(my_edge))
                for la, lb in needed_sig:
                    forced_pairs.append((offset + la, offset + lb))
                # Remember full sig for diversity scoring
                assigned_neighbor_sigs.append(edge_pairings_cache[j])

            # Sample candidates and pick best by diversity + span
            best_matching: list[tuple[int, int]] | None = None
            best_sig: tuple[frozenset[tuple[int, int]], ...] | None = None
            best_score = -1.0

            k_samples = diversity_samples if total_eps >= 4 else 1
            # Tiles with too few EPs can't avoid pips — skip min_span for sampling
            effective_min_span = min_span if total_eps >= 2 * min_span else 0

            for _ in range(k_samples):
                try:
                    matching = sample_constrained_ncpm(
                        total_eps, forced_pairs, rng,
                        min_span=effective_min_span)
                except UnsatisfiableMatchingError:
                    if effective_min_span > 0:
                        # min_span made it unsatisfiable — fall back without it
                        try:
                            matching = sample_constrained_ncpm(
                                total_eps, forced_pairs, rng, min_span=0)
                        except UnsatisfiableMatchingError:
                            return None  # true dead end — retry
                    else:
                        return None  # dead end — retry

                sig = edge_local_pairing(matching, tt)

                if k_samples == 1:
                    best_matching = matching
                    best_sig = sig
                    break

                # Score: prefer long-span pairs + neighbor diversity
                half = total_eps // 2
                span_score = sum(
                    min(b - a, total_eps - (b - a)) for a, b in matching
                ) / half

                div_score = sum(1 for ns in assigned_neighbor_sigs if ns != sig)

                score = span_score + div_score * 0.1
                if score > best_score:
                    best_score = score
                    best_matching = matching
                    best_sig = sig

            assignments[i] = best_matching
            edge_pairings_cache[i] = best_sig
            forced_per_tile[i] = forced_pairs

        return assignments, forced_per_tile

    for attempt in range(max_retries):
        result = _try_once(rng)
        if result is not None:
            if attempt > 0:
                print(f"  [greedy] succeeded on attempt {attempt + 1}")
            assignments, forced_per_tile = result
            if avoid_crossings:
                for i, pos in enumerate(positions):
                    assignments[i] = _refine_no_bezier_crossings(
                        assignments[i], pos.tile_type, pos.vertices,
                        forced_per_tile[i])
            return assignments
        # Reseed for next attempt
        rng = _random_module.Random(rng.randint(0, 2**32))

    raise UnsatisfiableMatchingError(
        f"Greedy solver: failed after {max_retries} attempts. "
        f"Try reducing max_eps or using a different seed.")


def solve_matching_assignment(
    positions: list[TilePosition],
    strategy: SolveStrategy = SolveStrategy.GROUPED,
    diverse: bool = True,
) -> list[int]:
    """SAT-solve which NCPM index each tile position gets.

    For each shared edge between positions[i].edge[e1] and positions[j].edge[e2]:
    precompute compatible (matching_i, matching_j) pairs by comparing
    edge pairings (with reversal), then encode as Z3 constraint.

    strategy controls how constraints are built:
    - BASELINE: enumerate all matching pairs (O(n²) per edge)
    - GROUPED: group by signature for O(1) lookup per signature

    Returns list of NCPM indices.
    Raises UnsatisfiableMatchingError if unsatisfiable (indicates a bug in tile construction).
    """
    t_start = perf_counter()

    # Two-phase: separate solver with edge-level Z3 variables
    if strategy == SolveStrategy.TWOPHASE:
        sig_cache: dict[int, SignatureData] = {}
        for pos in positions:
            key = id(pos.tile_type)
            if key not in sig_cache:
                sig_cache[key] = compute_signature_data(pos.tile_type)
        result = _solve_twophase(positions, sig_cache, diverse)
        t_total = perf_counter() - t_start
        print(f"  [solve] twophase | {len(positions)} tiles | "
              f"total: {t_total*1000:.1f}ms")
        return result

    # Cache by tile type identity
    if strategy == SolveStrategy.BASELINE:
        type_cache: dict[int, list[tuple[frozenset[tuple[int, int]], ...]]] = {}
        for pos in positions:
            key = id(pos.tile_type)
            if key not in type_cache:
                type_cache[key] = matching_edge_pairings(pos.tile_type)
        n_matchings_for = lambda pos: len(type_cache[id(pos.tile_type)])
    else:
        sig_cache: dict[int, SignatureData] = {}
        for pos in positions:
            key = id(pos.tile_type)
            if key not in sig_cache:
                sig_cache[key] = compute_signature_data(pos.tile_type)
        n_matchings_for = lambda pos: sig_cache[id(pos.tile_type)].n_matchings

    # Pre-filter: for each tile, only keep matchings compatible with at
    # least one neighbor matching on every shared edge.
    if strategy == SolveStrategy.GROUPED:
        valid_matchings: list[set[int] | None] = [None] * len(positions)
        for i, pos in enumerate(positions):
            if not pos.neighbors:
                continue
            sd_i = sig_cache[id(pos.tile_type)]
            valid = set(range(sd_i.n_matchings))
            for j, my_edge, their_edge in pos.neighbors:
                sd_j = sig_cache[id(positions[j].tile_type)]
                k = len(pos.tile_type.edge_points[my_edge])
                # Which of my matchings have a compatible partner on this edge?
                compatible_on_edge: set[int] = set()
                for sig_i, members_i in sd_i.edge_groups[my_edge].items():
                    sig_j_needed = reversed_edge_pairing(sig_i, k)
                    if sig_j_needed in sd_j.edge_groups[their_edge]:
                        compatible_on_edge.update(members_i)
                valid &= compatible_on_edge
            valid_matchings[i] = valid

        n_before = sum(n_matchings_for(pos) for pos in positions)
        n_after = sum(len(vm) if vm is not None else n_matchings_for(pos)
                      for pos, vm in zip(positions, valid_matchings))
        print(f"  [prefilter] {n_before} -> {n_after} total matchings "
              f"({100*n_after/n_before:.1f}%)")

    if diverse:
        s = Optimize()
    else:
        s = Solver()
        s.set("threads", os.cpu_count() or 1)
    variables: list[Int] = []
    for i, pos in enumerate(positions):
        v = Int(f"m_{i}")
        if strategy == SolveStrategy.GROUPED and valid_matchings[i] is not None:
            vm = sorted(valid_matchings[i])
            s.add(Or([v == m for m in vm]))
        else:
            s.add(v >= 0, v < n_matchings_for(pos))
        variables.append(v)

    if strategy == SolveStrategy.BASELINE:
        _add_baseline_constraints(s, variables, positions, type_cache)
    else:
        _add_grouped_constraints(s, variables, positions, sig_cache)

    if diverse:
        # Collect unique adjacent pairs
        adjacent_pairs: list[tuple[int, int]] = []
        seen: set[frozenset[int]] = set()
        for i, pos in enumerate(positions):
            for j, _, _ in pos.neighbors:
                pair_key = frozenset({i, j})
                if pair_key not in seen:
                    seen.add(pair_key)
                    adjacent_pairs.append((i, j))
        # Maximize number of adjacent pairs with different matchings
        if adjacent_pairs:
            diversity = sum(If(variables[i] != variables[j], 1, 0)
                            for i, j in adjacent_pairs)
            s.maximize(diversity)

    # Time the Z3 solver check
    t_solve_start = perf_counter()
    check_result = s.check()
    t_solve = perf_counter() - t_solve_start

    if check_result != sat:
        raise UnsatisfiableMatchingError(
            f"Z3 solver could not satisfy matching assignment constraints for "
            f"{len(positions)} tiles. This typically indicates incompatible tile "
            f"types or edge point configurations that cannot form a valid tiling.",
            positions=positions
        )

    model = s.model()
    result = [model[v].as_long() for v in variables]

    t_total = perf_counter() - t_start
    mode = f"{strategy.value}+diverse" if diverse else strategy.value
    print(f"  [solve] {mode} | {len(positions)} tiles | "
          f"Z3: {t_solve*1000:.1f}ms | total: {t_total*1000:.1f}ms")

    return result


# ---------------------------------------------------------------------------
# Grid builders — create list[TilePosition] with adjacency pre-computed
# ---------------------------------------------------------------------------

def square_grid(rows: int, cols: int, tile_type: TileType,
                s: float = 1.0) -> list[TilePosition]:
    """Grid of square tiles."""
    positions: list[TilePosition] = []
    for r in range(rows):
        for c in range(cols):
            verts = list(rect(c * s, r * s, s, s))
            positions.append(TilePosition(tile_type, verts))
    find_adjacency(positions)
    return positions


def equilateral_triangle_grid(n_rows: int, n_cols: int,
                              s: float = 1.0) -> list[TilePosition]:
    """Grid of equilateral triangles (up/down alternating).

    n_cols is triangles per row. (r+c) even = up, odd = down.
    """
    tt = TileType.uniform("tri-eq", 3, [1/3, 2/3])
    h = s * sqrt(3) / 2
    positions: list[TilePosition] = []
    for r in range(n_rows):
        for c in range(n_cols):
            x = c * s / 2
            y = r * h
            if (r + c) % 2 == 0:
                verts = [(x, y + h), (x + s, y + h), (x + s / 2, y)]
            else:
                verts = [(x + s / 2, y + h), (x + s, y), (x, y)]
            positions.append(TilePosition(tt, verts))
    find_adjacency(positions)
    return positions


def right_triangle_45_grid(rows: int, cols: int,
                           s: float = 1.0) -> list[TilePosition]:
    """Pairs of 45-45-90 triangles filling a square grid."""
    tt = TileType.uniform("tri45", 3, [1/3, 2/3])
    positions: list[TilePosition] = []
    for r in range(rows):
        for c in range(cols):
            x, y = c * s, r * s
            tl = (x, y)
            tr = (x + s, y)
            br = (x + s, y + s)
            bl = (x, y + s)
            positions.append(TilePosition(tt, [tl, tr, bl]))
            positions.append(TilePosition(tt, [br, bl, tr]))
    find_adjacency(positions)
    return positions


def hex_grid(rows: int, cols: int, L: float = 1.0,
             tile_type: TileType | None = None) -> list[TilePosition]:
    """Honeycomb grid of hexagons."""
    if tile_type is None:
        tile_type = TileType.uniform("hex", 6, [0.5])
    h_hex = L * sqrt(3)
    positions: list[TilePosition] = []
    for j in range(cols):
        for i in range(rows):
            x = j * 3 * L / 2
            y = i * h_hex - (j % 2) * h_hex / 2
            verts = hexagon_in_rect(x, y, L)
            positions.append(TilePosition(tile_type, verts))
    find_adjacency(positions)
    return positions


def hex_corners_grid(rows: int, cols: int,
                     L: float = 1.0) -> list[TilePosition]:
    """Grid of hexagons with 30-60-90 corner triangles filling gaps."""
    hex_tt = TileType.uniform("hex", 6, [1/3, 2/3])
    tri_tt = TileType.uniform("tri30", 3, [1/3, 2/3])
    w = 2 * L
    h = L * sqrt(3)
    positions: list[TilePosition] = []
    for r in range(rows):
        for c in range(cols):
            x = c * w
            y = r * h
            hex_v = hexagon_in_rect(x, y, L)
            positions.append(TilePosition(hex_tt, hex_v))
            r_tl = (x, y)
            r_tr = (x + w, y)
            r_br = (x + w, y + h)
            r_bl = (x, y + h)
            for tv in [
                [r_tl, hex_v[0], hex_v[5]],
                [r_tr, hex_v[2], hex_v[1]],
                [r_br, hex_v[3], hex_v[2]],
                [r_bl, hex_v[5], hex_v[4]],
            ]:
                positions.append(TilePosition(tri_tt, tv))
    find_adjacency(positions)
    return positions


def uniform_488_grid(n: int = 2, s: float = 1.0) -> list[TilePosition]:
    """4.8.8 grid with uniform octagon [0.5] + square [0.5] companions."""

    oct_tt = TileType.uniform("oct", 8, [0.5])
    sq_tt = TileType.uniform("sq", 4, [0.5])
    c = s / (2 + sqrt(2))
    positions: list[TilePosition] = []
    for r in range(n):
        for col in range(n):
            oct_v = octagon_in_square(col * s, r * s, s)
            positions.append(TilePosition(oct_tt, oct_v))
    for r in range(n - 1):
        for col in range(n - 1):
            cx = (col + 1) * s
            cy = (r + 1) * s
            sq_v = [(cx, cy - c), (cx + c, cy),
                    (cx, cy + c), (cx - c, cy)]
            positions.append(TilePosition(sq_tt, sq_v))
    find_adjacency(positions)
    return positions


def mixed_488_grid(n: int = 2, s: float = 1.0) -> list[TilePosition]:
    """4.8.8 grid with mixed octagon + corner triangles + enriched border squares."""

    oct_tt = oct_tile_type()
    tri_tt = corner_tri_tile_type()
    positions: list[TilePosition] = []

    for r in range(n):
        for col in range(n):
            x, y = col * s, r * s
            oct_v = octagon_in_square(x, y, s)
            sq_corners = list(rect(x, y, s, s))
            tl, tr, br, bl = sq_corners
            positions.append(TilePosition(oct_tt, oct_v))
            positions.append(TilePosition(tri_tt, [tl, oct_v[0], oct_v[7]]))
            positions.append(TilePosition(tri_tt, [tr, oct_v[2], oct_v[1]]))
            positions.append(TilePosition(tri_tt, [br, oct_v[4], oct_v[3]]))
            positions.append(TilePosition(tri_tt, [bl, oct_v[6], oct_v[5]]))

    for r in range(n):
        for col in range(n):
            x, y = col * s, r * s
            if r == 0:
                positions.append(TilePosition(
                    enriched_square_tile_type(2), list(rect(x, y - s, s, s))))
            if r == n - 1:
                positions.append(TilePosition(
                    enriched_square_tile_type(0), list(rect(x, y + s, s, s))))
            if col == 0:
                positions.append(TilePosition(
                    enriched_square_tile_type(1), list(rect(x - s, y, s, s))))
            if col == n - 1:
                positions.append(TilePosition(
                    enriched_square_tile_type(3), list(rect(x + s, y, s, s))))

    find_adjacency(positions)
    return positions


def solve_tile_types(positions: list[TilePosition],
                     min_eps: int = 1, max_eps: int = 3) -> None:
    """Use Z3 to solve for EP counts per edge, then replace tile_types in-place.

    Groups edges into equivalence classes (shared edges = same class),
    then solves one variable per class. This keeps the Z3 problem tiny
    regardless of tile count.

    Constraints:
    - Shared edges have the same EP count
    - Total EPs per tile is even
    - At least one tile has non-uniform EP counts (force irregularity)
    """
    # Build edge equivalence classes via union-find
    parent: dict[tuple[int, int], tuple[int, int]] = {}

    def find(k: tuple[int, int]) -> tuple[int, int]:
        while parent.get(k, k) != k:
            parent[k] = parent.get(parent[k], parent[k])
            k = parent[k]
        return k

    def union(a: tuple[int, int], b: tuple[int, int]):
        ra, rb = find(a), find(b)
        if ra != rb:
            parent[ra] = rb

    # Initialize all edges
    for i, pos in enumerate(positions):
        for e in range(pos.tile_type.n_edges):
            parent[(i, e)] = (i, e)

    # Merge shared edges
    for i, pos in enumerate(positions):
        for j_idx, my_edge, their_edge in pos.neighbors:
            union((i, my_edge), (j_idx, their_edge))

    # Map classes to Z3 variables
    class_vars: dict[tuple[int, int], Int] = {}
    class_id = 0
    for i, pos in enumerate(positions):
        for e in range(pos.tile_type.n_edges):
            root = find((i, e))
            if root not in class_vars:
                class_vars[root] = Int(f"ec_{class_id}")
                class_id += 1

    opt = Optimize()

    # Range constraints
    all_vars = list(class_vars.values())
    for v in all_vars:
        opt.add(v >= min_eps, v <= max_eps)

    # Evenness per tile
    for i, pos in enumerate(positions):
        edge_vars = [class_vars[find((i, e))]
                     for e in range(pos.tile_type.n_edges)]
        opt.add(Sum(edge_vars) % 2 == 0)

    # Minimize total EPs — keeps matchings tractable while naturally
    # producing irregular configs (boundary edges get lower counts)
    opt.minimize(Sum(all_vars))

    n_classes = len(class_vars)
    if opt.check() != sat:
        raise UnsatisfiableMatchingError("Cannot solve EP counts")

    model = opt.model()

    # Replace tile_types with solved EP configurations
    for i, pos in enumerate(positions):
        n = pos.tile_type.n_edges
        eps_per_edge = [model.eval(class_vars[find((i, e))]).as_long()
                        for e in range(n)]
        edge_points = [[k / (ep + 1) for k in range(1, ep + 1)]
                       for ep in eps_per_edge]
        name = f"{pos.tile_type.name}-{''.join(str(e) for e in eps_per_edge)}"
        pos.tile_type = TileType(name, n, edge_points)


def solve_tile_types_fast(positions: list[TilePosition],
                          min_eps: int = 1, max_eps: int = 4,
                          seed: int | None = None) -> None:
    """Solve EP counts per tile-type edge using GF(2) parity + random assignment.

    All tiles of the same type get identical edge points (Smith tile property).
    Uses union-find on (tile_type, edge_index) pairs to group edges that must
    match, solves parity over GF(2), then assigns random values.

    Replaces tile_types on positions in-place.
    """
    import random
    rng = random.Random(seed)

    # --- Union-find on (tile_type_id, edge_index) ---
    # Key insight: adjacency between tile instances constrains tile *types*.
    # If tile i (type A) edge e1 shares with tile j (type B) edge e2,
    # then type A's edge e1 must have the same EP count as type B's edge e2.
    parent: dict[tuple[int, int], tuple[int, int]] = {}

    def find(k: tuple[int, int]) -> tuple[int, int]:
        while parent.get(k, k) != k:
            parent[k] = parent.get(parent[k], parent[k])
            k = parent[k]
        return k

    def union(a: tuple[int, int], b: tuple[int, int]):
        ra, rb = find(a), find(b)
        if ra != rb:
            parent[ra] = rb

    # Collect unique tile types and initialize edges
    type_ids: dict[int, int] = {}  # id(tile_type) -> type_index
    type_objs: list[TileType] = []
    for pos in positions:
        key = id(pos.tile_type)
        if key not in type_ids:
            type_ids[key] = len(type_objs)
            type_objs.append(pos.tile_type)

    for ti, tt in enumerate(type_objs):
        for e in range(tt.n_edges):
            parent[(ti, e)] = (ti, e)

    # Merge edges from adjacency: type_of(i).edge[e1] == type_of(j).edge[e2]
    for pos in positions:
        ti = type_ids[id(pos.tile_type)]
        for j_idx, my_edge, their_edge in pos.neighbors:
            tj = type_ids[id(positions[j_idx].tile_type)]
            union((ti, my_edge), (tj, their_edge))

    # Collect unique edge classes
    class_ids: dict[tuple[int, int], int] = {}
    for ti, tt in enumerate(type_objs):
        for e in range(tt.n_edges):
            root = find((ti, e))
            if root not in class_ids:
                class_ids[root] = len(class_ids)

    n_classes = len(class_ids)

    # Per tile-type: list of edge class indices
    type_classes: list[list[int]] = []
    for ti, tt in enumerate(type_objs):
        type_classes.append([class_ids[find((ti, e))]
                             for e in range(tt.n_edges)])

    # --- GF(2) parity system ---
    # One equation per tile type: sum of edge parities ≡ 0 (mod 2)
    n_types = len(type_objs)
    rows: list[int] = []
    for classes in type_classes:
        row = 0
        for c in classes:
            row ^= (1 << c)
        rows.append(row)

    # Gaussian elimination
    pivot_col: list[int] = []
    rank = 0
    for col in range(n_classes):
        found = -1
        for r in range(rank, n_types):
            if rows[r] & (1 << col):
                found = r
                break
        if found == -1:
            continue
        rows[rank], rows[found] = rows[found], rows[rank]
        for r in range(n_types):
            if r != rank and rows[r] & (1 << col):
                rows[r] ^= rows[rank]
        pivot_col.append(col)
        rank += 1

    # Available values by parity
    even_vals = [v for v in range(min_eps, max_eps + 1) if v % 2 == 0]
    odd_vals = [v for v in range(min_eps, max_eps + 1) if v % 2 == 1]

    if not even_vals or not odd_vals:
        vals = even_vals or odd_vals
        values = [rng.choice(vals) for _ in range(n_classes)]
    else:
        parity = [0] * n_classes

        # Free variables get random parity
        pivot_set = set(pivot_col)
        for col in range(n_classes):
            if col not in pivot_set:
                parity[col] = rng.randint(0, 1)

        # Back-substitute
        for i in range(rank - 1, -1, -1):
            col = pivot_col[i]
            s = 0
            row = rows[i]
            for c in range(n_classes):
                if c != col and row & (1 << c):
                    s ^= parity[c]
            parity[col] = s

        values = [rng.choice(odd_vals if p == 1 else even_vals)
                  for p in parity]

    # Verify
    for classes in type_classes:
        assert sum(values[c] for c in classes) % 2 == 0

    # Build new TileTypes and replace on all positions
    new_types: dict[int, TileType] = {}  # type_index -> new TileType
    for ti, tt in enumerate(type_objs):
        eps_per_edge = [values[class_ids[find((ti, e))]]
                        for e in range(tt.n_edges)]
        edge_points = [[k / (ep + 1) for k in range(1, ep + 1)]
                       for ep in eps_per_edge]
        name = f"{tt.name}-{''.join(str(e) for e in eps_per_edge)}"
        new_types[ti] = TileType(name, tt.n_edges, edge_points)

    for pos in positions:
        ti = type_ids[id(pos.tile_type)]
        pos.tile_type = new_types[ti]


def build_plane_and_color(positions: list[TilePosition]
                          ) -> tuple[list[Tile], dict[int, int]]:
    """Build a Plane from tile positions and solve for minimum coloring.

    Returns (topo_tiles, coloring) where topo_tiles[i] is the Plane Tile
    for positions[i], and coloring maps face.id -> color index.
    """
    from libsat import min_coloring
    plane = Plane()
    topo_tiles = []
    for pos in positions:
        tile = plane.place_tile(pos.tile_type, pos.vertices)
        topo_tiles.append(tile)
    coloring = min_coloring(plane)
    return topo_tiles, coloring


# ---------------------------------------------------------------------------
# Tests
# ---------------------------------------------------------------------------

import unittest


class TestComposeSquare(unittest.TestCase):
    def test_single_square(self):
        comp = make_square()
        self.assertEqual(len(comp.plane.tiles), 1)
        self.assertEqual(len(comp.plane.faces), 8)
        self.assertEqual(len(comp.tile_outlines), 1)
        self._assert_no_conflicts(comp.plane)

    def test_custom_tile_type(self):
        tt = TileType("sq-irreg", 4,
                       [[0.5], [0.25, 0.5, 0.75], [0.5], [0.25, 0.5, 0.75]])
        comp = make_square(tile_type=tt)
        self.assertEqual(len(comp.plane.tiles), 1)
        self._assert_no_conflicts(comp.plane)

    def _assert_no_conflicts(self, plane: Plane):
        same = {frozenset({c.face_a.id, c.face_b.id})
                for c in plane.constraints if c.kind == ConstraintKind.SAME}
        diff = {frozenset({c.face_a.id, c.face_b.id})
                for c in plane.constraints if c.kind == ConstraintKind.DIFFERENT}
        self.assertEqual(len(same & diff), 0)


class TestComposeSquareGrid(unittest.TestCase):
    def test_2x2_grid(self):
        comp = make_square_grid(rows=2, cols=2)
        self.assertEqual(len(comp.plane.tiles), 4)
        self.assertEqual(len(comp.tile_outlines), 4)
        self._assert_no_conflicts(comp.plane)

    def test_custom_tile_type(self):
        tt = TileType("sq-irreg", 4,
                       [[0.5], [0.25, 0.5, 0.75], [0.5], [0.25, 0.5, 0.75]])
        comp = make_square_grid(rows=2, cols=2, tile_type=tt)
        self.assertEqual(len(comp.plane.tiles), 4)
        self._assert_no_conflicts(comp.plane)

    def _assert_no_conflicts(self, plane: Plane):
        same = {frozenset({c.face_a.id, c.face_b.id})
                for c in plane.constraints if c.kind == ConstraintKind.SAME}
        diff = {frozenset({c.face_a.id, c.face_b.id})
                for c in plane.constraints if c.kind == ConstraintKind.DIFFERENT}
        self.assertEqual(len(same & diff), 0)


class TestCompose488(unittest.TestCase):
    def test_tile_count(self):
        comp = make_488()
        self.assertEqual(len(comp.plane.tiles), 9)

    def test_face_count(self):
        comp = make_488()
        # All faces belong to a tile
        for f in comp.plane.faces:
            self.assertIsNotNone(f.tile)

    def test_outlines(self):
        comp = make_488()
        # 1 oct + 4 tri + 4 sq = 9
        self.assertEqual(len(comp.tile_outlines), 9)

    def test_no_conflicts(self):
        comp = make_488()
        same = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.SAME}
        diff = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.DIFFERENT}
        self.assertEqual(len(same & diff), 0)

    def test_compose_onto_existing_plane(self):
        plane = Plane()
        comp = compose_488(plane, 0, 0, 1)
        self.assertIs(comp.plane, plane)
        self.assertEqual(len(plane.tiles), 9)


class TestComposeHexCorners(unittest.TestCase):
    def test_tile_count(self):
        comp = make_hex_corners()
        self.assertEqual(len(comp.plane.tiles), 5)

    def test_outlines(self):
        comp = make_hex_corners()
        # 1 hex + 4 tri = 5
        self.assertEqual(len(comp.tile_outlines), 5)

    def test_shared_faces(self):
        comp = make_hex_corners()
        hex_face_ids = {id(f) for f in comp.plane.tiles[0].faces}
        shared = 0
        for t in comp.plane.tiles[1:]:
            for f in t.faces:
                if id(f) in hex_face_ids:
                    shared += 1
        self.assertEqual(shared, 12)  # 4 edges x 3 faces

    def test_no_conflicts(self):
        comp = make_hex_corners()
        same = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.SAME}
        diff = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.DIFFERENT}
        self.assertEqual(len(same & diff), 0)


class TestComposeHexStacked(unittest.TestCase):
    def test_two_stacked(self):
        plane = Plane()
        comp = compose_hex_stacked(plane, 0, 0, 1, n=2)
        self.assertEqual(len(comp.plane.tiles), 2)
        self.assertEqual(len(comp.tile_outlines), 2)

    def test_shared_edge(self):
        plane = Plane()
        comp = compose_hex_stacked(plane, 0, 0, 1, n=2)
        h1_ids = {id(f) for f in comp.plane.tiles[0].faces}
        shared = [f for f in comp.plane.tiles[1].faces if id(f) in h1_ids]
        self.assertEqual(len(shared), 3)


class TestComposeEquilateralCell(unittest.TestCase):
    def test_tile_count(self):
        plane = Plane()
        comp = compose_equilateral_cell(plane, 0, 0, 1)
        # 4 equilateral + 4 corner = 8
        self.assertEqual(len(comp.plane.tiles), 8)

    def test_outlines(self):
        plane = Plane()
        comp = compose_equilateral_cell(plane, 0, 0, 1)
        self.assertEqual(len(comp.tile_outlines), 8)

    def test_no_conflicts(self):
        plane = Plane()
        comp = compose_equilateral_cell(plane, 0, 0, 1)
        same = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.SAME}
        diff = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.DIFFERENT}
        self.assertEqual(len(same & diff), 0)


class TestComposeHexEqGrid(unittest.TestCase):
    def test_3x2_grid(self):
        comp = make_hex_eq_grid(rows=2, cols=3, hex_cols=2)
        # 2 hex cols x 2 rows x 5 tiles + 2 eq rows x 8 tiles = 20 + 16 = 36
        self.assertEqual(len(comp.plane.tiles), 36)

    def test_no_conflicts(self):
        comp = make_hex_eq_grid(rows=2, cols=3, hex_cols=2)
        same = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.SAME}
        diff = {frozenset({c.face_a.id, c.face_b.id})
                for c in comp.plane.constraints
                if c.kind == ConstraintKind.DIFFERENT}
        self.assertEqual(len(same & diff), 0)

    def test_no_different_at_polygon_vertices(self):
        comp = make_hex_eq_grid(rows=2, cols=3, hex_cols=2)
        from libtopo import VertexKind
        for c in comp.plane.constraints:
            if c.kind == ConstraintKind.DIFFERENT:
                fa, fb = c.face_a, c.face_b
                common = ({fa.v_start.id, fa.v_end.id}
                          & {fb.v_start.id, fb.v_end.id})
                if common:
                    vid = next(iter(common))
                    for v in comp.plane.vertices:
                        if v.id == vid:
                            self.assertEqual(v.kind, VertexKind.EDGE_POINT)
                            break


class TestNoncrossingMatchings(unittest.TestCase):
    def test_zero(self):
        self.assertEqual(noncrossing_matchings(0), [[]])

    def test_two(self):
        self.assertEqual(noncrossing_matchings(2), [[(0, 1)]])

    def test_four_catalan_2(self):
        ms = noncrossing_matchings(4)
        self.assertEqual(len(ms), 2)  # C(2) = 2
        self.assertIn([(0, 1), (2, 3)], ms)
        self.assertIn([(0, 3), (1, 2)], ms)

    def test_six_catalan_5(self):
        ms = noncrossing_matchings(6)
        self.assertEqual(len(ms), 5)  # C(3) = 5

    def test_eight_catalan_14(self):
        ms = noncrossing_matchings(8)
        self.assertEqual(len(ms), 14)  # C(4) = 14

    def test_twelve_catalan_132(self):
        ms = noncrossing_matchings(12)
        self.assertEqual(len(ms), 132)  # C(6) = 132

    def test_odd_raises(self):
        with self.assertRaises(ValueError):
            noncrossing_matchings(3)

    def test_all_matchings_are_perfect(self):
        for matching in noncrossing_matchings(8):
            pts = set()
            for a, b in matching:
                pts.add(a)
                pts.add(b)
            self.assertEqual(pts, set(range(8)))


class TestTileMatchings(unittest.TestCase):
    def test_square_basic(self):
        tt = TileType.uniform("sq", 4, [0.5])
        verts = list(rect(0, 0, 1, 1))
        ms = tile_matchings(tt, verts)
        self.assertEqual(len(ms), 2)  # C(2) = 2
        for matching in ms:
            self.assertEqual(len(matching), 2)

    def test_triangle_2ep(self):
        tt = TileType.uniform("tri", 3, [1/3, 2/3])
        verts = [(0.0, 0.0), (1.0, 0.0), (0.5, 0.866)]
        ms = tile_matchings(tt, verts)
        self.assertEqual(len(ms), 5)  # C(3) = 5


class TestMakeConvenience(unittest.TestCase):
    """Verify make_* functions are thin wrappers creating fresh Planes."""

    def test_make_square_creates_fresh_plane(self):
        c1 = make_square()
        c2 = make_square()
        self.assertIsNot(c1.plane, c2.plane)

    def test_make_square_grid_creates_fresh_plane(self):
        c1 = make_square_grid()
        c2 = make_square_grid()
        self.assertIsNot(c1.plane, c2.plane)

    def test_make_488_creates_fresh_plane(self):
        c1 = make_488()
        c2 = make_488()
        self.assertIsNot(c1.plane, c2.plane)

    def test_make_hex_corners_creates_fresh_plane(self):
        c1 = make_hex_corners()
        c2 = make_hex_corners()
        self.assertIsNot(c1.plane, c2.plane)

    def test_make_hex_eq_grid_creates_fresh_plane(self):
        c1 = make_hex_eq_grid()
        c2 = make_hex_eq_grid()
        self.assertIsNot(c1.plane, c2.plane)


class TestEdgeLocalPairing(unittest.TestCase):
    def test_square_basic_ncpm0(self):
        """Square [0.5]: NCPM [(0,1),(2,3)] — all points pair across edges."""
        tt = TileType.uniform("sq", 4, [0.5])
        matching = [(0, 1), (2, 3)]
        pairings = edge_local_pairing(matching, tt)
        self.assertEqual(len(pairings), 4)
        for p in pairings:
            self.assertEqual(p, frozenset())

    def test_square_2ep_local_pair(self):
        """Square [1/3, 2/3]: NCPM that pairs within an edge."""
        tt = TileType.uniform("sq", 4, [1/3, 2/3])
        # NCPM [(0,1), (2,3), (4,5), (6,7)] — all local pairs
        matching = [(0, 1), (2, 3), (4, 5), (6, 7)]
        pairings = edge_local_pairing(matching, tt)
        for p in pairings:
            self.assertEqual(p, frozenset({(0, 1)}))

    def test_square_2ep_cross_edge(self):
        """Square [1/3, 2/3]: NCPM that pairs across edges."""
        tt = TileType.uniform("sq", 4, [1/3, 2/3])
        # NCPM [(0,7), (1,2), (3,4), (5,6)]
        matching = [(0, 7), (1, 2), (3, 4), (5, 6)]
        pairings = edge_local_pairing(matching, tt)
        for p in pairings:
            self.assertEqual(p, frozenset())


class TestReversedEdgePairing(unittest.TestCase):
    def test_empty(self):
        self.assertEqual(reversed_edge_pairing(frozenset(), 2), frozenset())

    def test_single_pair_k2(self):
        # (0,1) with k=2: reversed to (0,1) — symmetric
        self.assertEqual(
            reversed_edge_pairing(frozenset({(0, 1)}), 2),
            frozenset({(0, 1)}))

    def test_single_pair_k3(self):
        # (0,1) with k=3: reversed to (1,2)
        self.assertEqual(
            reversed_edge_pairing(frozenset({(0, 1)}), 3),
            frozenset({(1, 2)}))


class TestSolveMatchingAssignment(unittest.TestCase):
    def test_square_grid_2x2(self):
        """2x2 grid of square [0.5] tiles — always satisfiable."""
        tt = TileType.uniform("sq", 4, [0.5])
        positions: list[TilePosition] = []
        for r in range(2):
            for c in range(2):
                verts = list(rect(c, r, 1, 1))
                positions.append(TilePosition(tt, verts))
        find_adjacency(positions)
        result = solve_matching_assignment(positions)
        self.assertEqual(len(result), 4)

    def test_square_2ep_grid_3x3(self):
        """3x3 grid of square [1/3, 2/3] — satisfiable with edge constraints."""
        tt = TileType.uniform("sq", 4, [1/3, 2/3])
        positions: list[TilePosition] = []
        for r in range(3):
            for c in range(3):
                verts = list(rect(c, r, 1, 1))
                positions.append(TilePosition(tt, verts))
        find_adjacency(positions)
        result = solve_matching_assignment(positions)
        self.assertEqual(len(result), 9)
        # Verify edge compatibility
        self._check_edge_compatibility(positions, result)

    def test_triangle_grid(self):
        """Small triangle grid — satisfiable."""
        tt = TileType.uniform("tri", 3, [1/3, 2/3])
        h = sqrt(3) / 2
        s = 1.0
        # 4 triangles: 2 up + 2 down
        positions = [
            TilePosition(tt, [(0, h), (s, h), (s/2, 0)]),       # up 0
            TilePosition(tt, [(s/2, h), (s, 0), (3*s/2, h)]),   # down 0 (reversed winding)
        ]
        find_adjacency(positions)
        result = solve_matching_assignment(positions)
        self._check_edge_compatibility(positions, result)

    def test_unsatisfiable_raises_exception(self):
        """Incompatible tile edge configurations should raise UnsatisfiableMatchingError."""
        # Create tiles with incompatible edge point structures
        # tt1: edges with different point counts - [1 pt, 2 pts, 1 pt, 2 pts]
        # tt2: edges all with 1 point
        # When they share an edge, tt1's 2-point edge can't match tt2's 1-point edge
        tt1 = TileType("sq1", 4, [[0.5], [1/3, 2/3], [0.5], [1/3, 2/3]])
        tt2 = TileType.uniform("sq2", 4, [0.5])

        # Place them so tt1's edge 1 (2 points) shares with tt2's edge 3 (1 point)
        positions = [
            TilePosition(tt1, list(rect(0, 0, 1, 1))),
            TilePosition(tt2, list(rect(1, 0, 1, 1))),
        ]
        find_adjacency(positions)

        with self.assertRaises(UnsatisfiableMatchingError) as ctx:
            solve_matching_assignment(positions)

        # Verify exception contains useful debugging info
        exc = ctx.exception
        self.assertIn("Incompatible edge point counts", exc.message)
        self.assertIsNotNone(exc.tile_i)
        self.assertIsNotNone(exc.tile_j)
        self.assertIsNotNone(exc.edge_i)
        self.assertIsNotNone(exc.edge_j)

    def test_square_grid_2x2_grouped(self):
        """2x2 grid with GROUPED strategy."""
        tt = TileType.uniform("sq", 4, [0.5])
        positions: list[TilePosition] = []
        for r in range(2):
            for c in range(2):
                verts = list(rect(c, r, 1, 1))
                positions.append(TilePosition(tt, verts))
        find_adjacency(positions)
        result = solve_matching_assignment(positions, strategy=SolveStrategy.GROUPED)
        self.assertEqual(len(result), 4)
        validate_matching_assignment(positions, result)

    def test_strategies_agree(self):
        """All strategies produce valid (not necessarily equal) assignments."""
        tt = TileType.uniform("sq", 4, [1/3, 2/3])
        positions = [TilePosition(tt, list(rect(c, r, 1, 1)))
                     for r in range(3) for c in range(3)]
        find_adjacency(positions)

        for strategy in SolveStrategy:
            result = solve_matching_assignment(positions, strategy=strategy)
            validate_matching_assignment(positions, result)

    def test_diverse_produces_valid_assignment(self):
        """diverse=True produces a valid assignment."""
        tt = TileType.uniform("sq", 4, [1/3, 2/3])
        positions = [TilePosition(tt, list(rect(c, r, 1, 1)))
                     for r in range(3) for c in range(3)]
        find_adjacency(positions)
        result = solve_matching_assignment(positions, diverse=True)
        validate_matching_assignment(positions, result)

    def _check_edge_compatibility(self, positions: list[TilePosition], assignment: list[int]):
        """Verify all shared edges have compatible pairings."""
        validate_matching_assignment(positions, assignment)


class TestFindAdjacency(unittest.TestCase):
    def test_two_squares(self):
        tt = TileType.uniform("sq", 4, [0.5])
        positions = [
            TilePosition(tt, list(rect(0, 0, 1, 1))),
            TilePosition(tt, list(rect(1, 0, 1, 1))),
        ]
        find_adjacency(positions)
        self.assertEqual(len(positions[0].neighbors), 1)
        self.assertEqual(len(positions[1].neighbors), 1)
        # Square 0 edge 1 (right) shares with square 1 edge 3 (left)
        n_idx, my_e, their_e = positions[0].neighbors[0]
        self.assertEqual(n_idx, 1)
        self.assertEqual(my_e, 1)
        self.assertEqual(their_e, 3)

    def test_no_shared_edges(self):
        tt = TileType.uniform("sq", 4, [0.5])
        positions = [
            TilePosition(tt, list(rect(0, 0, 1, 1))),
            TilePosition(tt, list(rect(5, 5, 1, 1))),
        ]
        find_adjacency(positions)
        self.assertEqual(len(positions[0].neighbors), 0)
        self.assertEqual(len(positions[1].neighbors), 0)


if __name__ == "__main__":
    unittest.main()