pcb/scripts/close_remaining.py

#!/usr/bin/env python3
"""Close every connection the staged flow leaves open, with one router.

The earlier residual scripts each hard-coded a shape: drop a via beside pad A,
A* on one nominated layer, drop a via beside pad B.  That shape cannot escape a
0.8-mm BGA ball whose four dogbone sites are taken, and it cannot bridge a
power pour that signal traces have cut into islands.  Those are exactly the
failures left on this board, so this stage replaces the shape with a real maze:

* the search runs over (x, y, layer) on all five routable layers at a 0.1-mm
  pitch and changes layer through a via wherever one legally fits;
* start and goal are *every* piece of copper in the two groups being joined -
  pad, track or via - so the router leaves from whichever is cheapest instead
  of one nominated pad;
* clearance uses the exact pad/track/via geometry at the real net-class values,
  including the board.kicad_dru rule that relaxes clearance against the Ti375
  lands to 0.10 mm inside the escape field;
* copper is necked to 0.10 mm inside the U1 courtyard, which is the only width
  that fits between two 0.8-mm-pitch lands, and runs at class width outside;
* cells inside a foreign non-ground pour are refused, so a new trace can never
  cut a power island in half.

Everything added is locked.  Run the KiCad DRC afterwards.
"""
from __future__ import annotations

import argparse
import heapq
import math
import os
import shutil
import time
from collections import defaultdict
from pathlib import Path

import numpy as np
from matplotlib.path import Path as MplPath

import pcbnew

import route_local_power as rlp

ROOT = Path(__file__).resolve().parents[1]

ROUTABLE = ("F.Cu", "L3_SIG", "L5_SIG", "L6_SIG", "B.Cu")
STEP = 0.1
VIA_PENALTY = 1.5               # mm of path cost charged for one layer change
VIA_D, VIA_DRILL = 0.45, 0.225
HOLE_TO_HOLE = 0.25             # board design rule
EDGE_MARGIN = 0.30

ESCAPE_REFS = frozenset({"U1"})
ESCAPE_CLEARANCE = 0.10         # "Ti375 0.8-mm BGA escape clearance" rule
FIELD_WIDTH = 0.10              # only width that threads two 0.8-mm lands

# Routing widths.  POWER's 0.4-mm class width is a plane-drop default, not a
# width the residual connections can carry through congested copper.
ROUTE_WIDTH = {"Default": 0.12, "BGA_ESCAPE": 0.10, "CLOCK": 0.15,
               "POWER": 0.25}

# Pours that a new trace must not cross, because splitting them strands copper.
# The board-wide GND pours simply re-pour around a trace, so they are not
# blockers.
POUR_SAFE_NETS = frozenset({"GND"})
POUR_MARGIN = 0.05

# Copper the router may rip up and hand back to itself to re-route.  The
# deterministic stages (power loops, RF/USB geometry, BGA dogbones) are all
# locked, so "unlocked and in a signal class" is exactly the autorouted copper.
RIPPABLE_CLASSES = frozenset({"Default", "BGA_ESCAPE", "CLOCK"})


def mm(v):
    return pcbnew.ToMM(v)


def MM(v):
    return pcbnew.FromMM(v)


# --------------------------------------------------------------------------
# board model
# --------------------------------------------------------------------------
def _project_clearances(board: pcbnew.BOARD) -> dict[str, float]:
    """Net-class name -> clearance in mm, read from the sibling .kicad_pro."""
    import json
    path = Path(board.GetFileName()).with_suffix(".kicad_pro")
    out = {"Default": 0.1}
    if not path.exists():
        return out
    data = json.loads(path.read_text())
    for entry in data.get("net_settings", {}).get("classes", []):
        if "clearance" in entry:
            out[entry["name"]] = float(entry["clearance"])
    return out


class Model:
    def __init__(self, board: pcbnew.BOARD) -> None:
        self.board = board
        self.layer_num = {board.GetLayerName(l): l
                          for l in board.GetEnabledLayers().CuStack()}
        self.routable = [self.layer_num[n] for n in ROUTABLE
                         if n in self.layer_num]
        self.layer_index = {l: i for i, l in enumerate(self.routable)}

        # Net-class clearances come from the project file: the pcbnew Python
        # binding hands back an opaque object for GetNetClass().
        self.class_clearance = _project_clearances(board)
        self.clearance: dict[str, float] = {}
        self.classname: dict[str, str] = {}
        for _code, info in board.GetNetsByNetcode().items():
            name = info.GetNetname()
            cls = info.GetNetClassName()
            self.classname[name] = cls
            self.clearance[name] = self.class_clearance.get(cls, 0.1)

        self.pads = []          # (shape, layerset, net, ref)
        self.pad_holes = []     # (pos, radius, drill, net, uuid, rippable)
        self.tracks = []        # (a, b, half, layer, net, uuid, rippable)
        self.vias = []          # pad holes + board vias, same shape
        for fp in board.GetFootprints():
            ref = fp.GetReference()
            for pad in fp.Pads():
                self.pads.append((rlp.pad_shape(pad), pad.GetLayerSet(),
                                  pad.GetNetname(), ref))
                drill = max(mm(pad.GetDrillSizeX()), mm(pad.GetDrillSizeY()))
                if drill > 0:
                    p = pad.GetPosition()
                    self.pad_holes.append(((mm(p.x), mm(p.y)), 0.0, drill / 2,
                                           pad.GetNetname(), None, False))
        self.load_copper()

        # filled pours, per copper layer, as polygons
        self.fills = []         # (net, layer, MplPath, bbox)
        for zone in board.Zones():
            if zone.GetIsRuleArea():
                continue
            net = zone.GetNetname()
            for layer in zone.GetLayerSet().CuStack():
                polys = zone.GetFilledPolysList(layer)
                if polys is None:
                    continue
                for i in range(polys.OutlineCount()):
                    self.fills.append((net, layer) + _outline_path(polys.Outline(i)))

        # rule areas that forbid copper outright
        self.keepouts = []      # (MplPath, bbox, layerset, no_track, no_via)
        self.rf_areas = []      # RF_NO_DIGITAL: blocks non-POWER track and via
        for zone in board.Zones():
            if not zone.GetIsRuleArea():
                continue
            outline = zone.Outline()
            for i in range(outline.OutlineCount()):
                path, bbox = _outline_path(outline.Outline(i))
                entry = (path, bbox, zone.GetLayerSet())
                if zone.GetZoneName() == "RF_NO_DIGITAL":
                    self.rf_areas.append(entry)
                if zone.GetDoNotAllowTracks() or zone.GetDoNotAllowVias():
                    self.keepouts.append(entry + (zone.GetDoNotAllowTracks(),
                                                  zone.GetDoNotAllowVias()))

        u1 = board.FindFootprintByReference("U1")
        self.u1_courtyard = None
        if u1 is not None:
            poly = u1.GetCourtyard(pcbnew.F_Cu)
            if poly and poly.OutlineCount():
                self.u1_courtyard = _outline_path(poly.Outline(0))

        box = self.board.GetBoardEdgesBoundingBox()
        self.edge = (mm(box.GetX()) + EDGE_MARGIN, mm(box.GetY()) + EDGE_MARGIN,
                     mm(box.GetRight()) - EDGE_MARGIN,
                     mm(box.GetBottom()) - EDGE_MARGIN)

    def load_copper(self) -> None:
        """(Re)read every track and via from the board.

        Copper is identified by KIID, never by the address of a SWIG proxy:
        the same board item can be handed out as several different Python
        objects, and a proxy for a removed item may be freed at any time.
        """
        self.tracks = []
        vias = []
        for t in self.board.GetTracks():
            net = t.GetNetname()
            uuid = t.m_Uuid.AsString()
            rippable = (not t.IsLocked() and bool(net) and
                        self.classname.get(net) in RIPPABLE_CLASSES)
            if t.GetClass() == "PCB_VIA":
                p = t.GetPosition()
                vias.append(((mm(p.x), mm(p.y)),
                             mm(t.GetWidth(pcbnew.F_Cu)) / 2,
                             mm(t.GetDrill()) / 2, net, uuid, rippable))
            else:
                self.tracks.append(((mm(t.GetStart().x), mm(t.GetStart().y)),
                                    (mm(t.GetEnd().x), mm(t.GetEnd().y)),
                                    mm(t.GetWidth()) / 2, t.GetLayer(),
                                    net, uuid, rippable))
        self.vias = self.pad_holes + vias

    def route_width(self, net: str) -> float:
        return ROUTE_WIDTH.get(self.classname.get(net, "Default"), 0.12)

    def clearance_to(self, net: str, other: str) -> float:
        return max(self.clearance.get(net, 0.1), self.clearance.get(other, 0.1))


def _outline_path(outline):
    pts = [(mm(outline.CPoint(i).x), mm(outline.CPoint(i).y))
           for i in range(outline.PointCount())]
    arr = np.asarray(pts)
    bbox = (arr[:, 0].min(), arr[:, 1].min(), arr[:, 0].max(), arr[:, 1].max())
    return MplPath(arr), bbox


# --------------------------------------------------------------------------
# connectivity groups
# --------------------------------------------------------------------------
def net_groups(board: pcbnew.BOARD, net: str):
    conn = board.GetConnectivity()
    items = [p for fp in board.GetFootprints() for p in fp.Pads()
             if p.GetNetname() == net]
    items += [t for t in board.GetTracks() if t.GetNetname() == net]
    index = {it.m_Uuid.AsString(): n for n, it in enumerate(items)}
    parent = list(range(len(items)))

    def find(a):
        while parent[a] != a:
            parent[a] = parent[parent[a]]
            a = parent[a]
        return a

    for i, it in enumerate(items):
        for other in conn.GetConnectedItems(it, 0):
            j = index.get(other.m_Uuid.AsString())
            if j is not None:
                ra, rb = find(i), find(j)
                if ra != rb:
                    parent[rb] = ra

    out = defaultdict(list)
    for i, it in enumerate(items):
        out[find(i)].append(it)
    return list(out.values())


# --------------------------------------------------------------------------
# rasterised routing field
# --------------------------------------------------------------------------
def _seg_distance(X, Y, a, b):
    ax, ay = a
    bx, by = b
    dx, dy = bx - ax, by - ay
    length2 = dx * dx + dy * dy
    if length2 <= 1e-12:
        return np.hypot(X - ax, Y - ay)
    t = np.clip(((X - ax) * dx + (Y - ay) * dy) / length2, 0.0, 1.0)
    return np.hypot(X - (ax + t * dx), Y - (ay + t * dy))


class Field:
    """Blocked bitmaps for one net over one rectangular region.

    `blocked[w][layer]` is True where a trace of half-width w may not sit;
    `via_blocked` is True where a through via may not sit.
    """

    def __init__(self, model: Model, net: str, region, widths,
                 ignore_rippable: bool = False,
                 protected: frozenset[str] = frozenset()) -> None:
        self.model = model
        self.net = net
        self.ignore_rippable = ignore_rippable
        self.protected = protected
        x0, y0, x1, y1 = region
        # Snap the grid to the absolute STEP lattice.  Ball centres, dogbone
        # sites and the mid-lines between two 0.8-mm lands are all whole
        # multiples of 0.1 mm; on an unaligned grid the legal corridor between
        # two lands (0.05 mm either side of its centre) falls between samples
        # and the escape field looks sealed.
        x0 = math.floor(x0 / STEP) * STEP
        y0 = math.floor(y0 / STEP) * STEP
        self.x0, self.y0 = x0, y0
        self.nx = int(math.ceil((x1 - x0) / STEP)) + 1
        self.ny = int(math.ceil((y1 - y0) / STEP)) + 1
        self.xs = x0 + STEP * np.arange(self.nx)
        self.ys = y0 + STEP * np.arange(self.ny)
        self.X, self.Y = np.meshgrid(self.xs, self.ys)
        self.nl = len(model.routable)
        self.widths = widths

        self.blocked = {w: np.zeros((self.nl, self.ny, self.nx), bool)
                        for w in widths}
        self.via_blocked = np.zeros((self.ny, self.nx), bool)

        self.in_u1 = self._mask_path(model.u1_courtyard) if model.u1_courtyard \
            else np.zeros((self.ny, self.nx), bool)
        self._mark_edges()
        self._mark_obstacles()
        self._mark_pours()
        self._mark_rule_areas()

    # -- helpers ----------------------------------------------------------
    def _window(self, bbox, pad):
        i0 = max(int(math.floor((bbox[0] - pad - self.x0) / STEP)), 0)
        j0 = max(int(math.floor((bbox[1] - pad - self.y0) / STEP)), 0)
        i1 = min(int(math.ceil((bbox[2] + pad - self.x0) / STEP)) + 1, self.nx)
        j1 = min(int(math.ceil((bbox[3] + pad - self.y0) / STEP)) + 1, self.ny)
        if i0 >= i1 or j0 >= j1:
            return None
        return slice(j0, j1), slice(i0, i1)

    def _mask_path(self, path_bbox):
        path, bbox = path_bbox
        mask = np.zeros((self.ny, self.nx), bool)
        win = self._window(bbox, 0.0)
        if win is None:
            return mask
        sy, sx = win
        pts = np.column_stack([self.X[sy, sx].ravel(), self.Y[sy, sx].ravel()])
        mask[sy, sx] = path.contains_points(pts).reshape(self.X[sy, sx].shape)
        return mask

    @staticmethod
    def _dilate(mask, radius: float):
        """Grow a boolean mask by `radius` mm.

        The keep-out rules are written with intersectsArea(), so a via or track
        violates as soon as its *body* touches the area - testing the centre
        cell alone lets copper hang over the edge.
        """
        cells = int(math.ceil(radius / STEP))
        if cells <= 0:
            return mask
        ny, nx = mask.shape
        out = mask.copy()
        for dj in range(-cells, cells + 1):
            for di in range(-cells, cells + 1):
                if (dj or di) and dj * dj + di * di <= cells * cells:
                    dst = (slice(max(0, -dj), ny - max(0, dj)),
                           slice(max(0, -di), nx - max(0, di)))
                    src = (slice(max(0, dj), ny - max(0, -dj)),
                           slice(max(0, di), nx - max(0, -di)))
                    out[dst] |= mask[src]
        return out

    def _in_field(self, points) -> bool:
        """True when every point lies inside the Ti375 escape field."""
        if self.model.u1_courtyard is None:
            return False
        path, _bbox = self.model.u1_courtyard
        return bool(path.contains_points(np.asarray(points)).all())

    def _threshold(self, sy, sx, clearance: float, inside: bool):
        """Clearance owed per cell: the escape-field minimum where both this
        cell and the obstacle sit inside U1's courtyard, the net-class value
        everywhere else.  Mirrors the "Ti375 escape field copper clearance"
        rule in board.kicad_dru."""
        if not inside or clearance <= ESCAPE_CLEARANCE:
            return clearance
        return np.where(self.in_u1[sy, sx], ESCAPE_CLEARANCE, clearance)

    def _mark_edges(self):
        ex0, ey0, ex1, ey1 = self.model.edge
        for w in self.widths:
            half = w / 2
            bad = ((self.X < ex0 + half) | (self.X > ex1 - half) |
                   (self.Y < ey0 + half) | (self.Y > ey1 - half))
            self.blocked[w][:, bad] = True
        half = VIA_D / 2
        self.via_blocked |= ((self.X < ex0 + half) | (self.X > ex1 - half) |
                             (self.Y < ey0 + half) | (self.Y > ey1 - half))

    def _mark_obstacles(self):
        model = self.model
        net = self.net
        for shape, layers, other, ref in model.pads:
            if other == net:
                continue
            cx, cy, hx, hy, r = shape
            bbox = (cx - hx - r, cy - hy - r, cx + hx + r, cy + hy + r)
            win = self._window(bbox, 0.6)
            if win is None:
                continue
            sy, sx = win
            if ref in ESCAPE_REFS:
                clearance = ESCAPE_CLEARANCE
            else:
                clearance = self._threshold(
                    sy, sx, model.clearance_to(net, other),
                    self._in_field([(cx - hx - r, cy - hy - r),
                                    (cx + hx + r, cy + hy + r),
                                    (cx - hx - r, cy + hy + r),
                                    (cx + hx + r, cy - hy - r)]))
            dx = np.maximum(np.abs(self.X[sy, sx] - cx) - hx, 0.0)
            dy = np.maximum(np.abs(self.Y[sy, sx] - cy) - hy, 0.0)
            dist = np.maximum(np.hypot(dx, dy) - r, 0.0)
            for li, layer in enumerate(model.routable):
                if not layers.Contains(layer):
                    continue
                for w in self.widths:
                    self.blocked[w][li, sy, sx] |= dist < clearance + w / 2
            self.via_blocked[sy, sx] |= dist < clearance + VIA_D / 2

        for a, b, half, layer, other, _uuid, rippable in model.tracks:
            if other == net or layer not in model.layer_index:
                continue
            if self.ignore_rippable and rippable and other not in self.protected:
                continue
            li = model.layer_index[layer]
            bbox = (min(a[0], b[0]), min(a[1], b[1]),
                    max(a[0], b[0]), max(a[1], b[1]))
            win = self._window(bbox, half + 0.6)
            if win is None:
                continue
            sy, sx = win
            clearance = self._threshold(sy, sx, model.clearance_to(net, other),
                                        self._in_field([a, b]))
            dist = _seg_distance(self.X[sy, sx], self.Y[sy, sx], a, b) - half
            for w in self.widths:
                self.blocked[w][li, sy, sx] |= dist < clearance + w / 2
            self.via_blocked[sy, sx] |= dist < clearance + VIA_D / 2

        for pos, radius, drill, other, _uuid, rippable in model.vias:
            soft = (self.ignore_rippable and rippable and
                    other not in self.protected)
            bbox = (pos[0] - radius, pos[1] - radius,
                    pos[0] + radius, pos[1] + radius)
            win = self._window(bbox, max(radius, drill) + 0.6)
            if win is None:
                continue
            sy, sx = win
            centre = np.hypot(self.X[sy, sx] - pos[0], self.Y[sy, sx] - pos[1])
            if not soft:
                # a hole always owes hole-to-hole, even to its own net
                self.via_blocked[sy, sx] |= (centre <
                                             drill + VIA_DRILL / 2 + HOLE_TO_HOLE)
            if other == net or soft:
                continue
            clearance = self._threshold(
                sy, sx, model.clearance_to(net, other),
                self._in_field([(pos[0] - radius, pos[1] - radius),
                                (pos[0] + radius, pos[1] + radius),
                                (pos[0] - radius, pos[1] + radius),
                                (pos[0] + radius, pos[1] - radius)]))
            dist = np.maximum(centre - radius, 0.0)
            for w in self.widths:
                self.blocked[w][:, sy, sx] |= dist < clearance + w / 2
            self.via_blocked[sy, sx] |= dist < clearance + VIA_D / 2

    def _mark_pours(self):
        """Refuse to route inside a foreign pour that a trace would split."""
        model = self.model
        for other, layer, path, bbox in model.fills:
            if other == self.net or other in POUR_SAFE_NETS:
                continue
            if layer not in model.layer_index:
                continue
            li = model.layer_index[layer]
            win = self._window(bbox, POUR_MARGIN)
            if win is None:
                continue
            sy, sx = win
            pts = np.column_stack([self.X[sy, sx].ravel(),
                                   self.Y[sy, sx].ravel()])
            inside = path.contains_points(pts, radius=POUR_MARGIN)
            inside = inside.reshape(self.X[sy, sx].shape)
            for w in self.widths:
                self.blocked[w][li, sy, sx] |= inside

    def _mark_rule_areas(self):
        model = self.model
        is_power = model.classname.get(self.net) == "POWER"
        for path, bbox, layers, no_track, no_via in model.keepouts:
            mask = self._mask_path((path, bbox))
            grown = {w: self._dilate(mask, w / 2) for w in self.widths}
            for li, layer in enumerate(model.routable):
                if no_track and layers.Contains(layer):
                    for w in self.widths:
                        self.blocked[w][li] |= grown[w]
            if no_via:
                self.via_blocked |= self._dilate(mask, VIA_D / 2)
        if is_power:
            return
        # board.kicad_dru: no track and no via for non-POWER classes inside
        # the GNSS front-end area, on every layer the area covers.
        for path, bbox, layers in model.rf_areas:
            mask = self._mask_path((path, bbox))
            grown = {w: self._dilate(mask, w / 2) for w in self.widths}
            covered = False
            for li, layer in enumerate(model.routable):
                if layers.Contains(layer):
                    covered = True
                    for w in self.widths:
                        self.blocked[w][li] |= grown[w]
            if covered:
                self.via_blocked |= self._dilate(mask, VIA_D / 2)

    # -- indexing ---------------------------------------------------------
    def cell(self, x, y):
        i = int(round((x - self.x0) / STEP))
        j = int(round((y - self.y0) / STEP))
        if 0 <= i < self.nx and 0 <= j < self.ny:
            return j, i
        return None

    def xy(self, j, i):
        # plain floats: pcbnew's FromMM rejects numpy scalars
        return (round(float(self.x0) + int(i) * STEP, 4),
                round(float(self.y0) + int(j) * STEP, 4))


# --------------------------------------------------------------------------
# search
# --------------------------------------------------------------------------
def seed_cells(field: Field, model: Model, items):
    """Grid cells that lie inside a group's existing copper."""
    cells = set()
    for it in items:
        cls = it.GetClass()
        if cls == "PAD":
            shape = rlp.pad_shape(it)
            cx, cy, hx, hy, r = shape
            bbox = (cx - hx - r, cy - hy - r, cx + hx + r, cy + hy + r)
            win = field._window(bbox, 0.0)
            if win is None:
                continue
            sy, sx = win
            dx = np.maximum(np.abs(field.X[sy, sx] - cx) - hx, 0.0)
            dy = np.maximum(np.abs(field.Y[sy, sx] - cy) - hy, 0.0)
            inside = np.maximum(np.hypot(dx, dy) - r, 0.0) <= 0.0
            layers = it.GetLayerSet()
            for li, layer in enumerate(model.routable):
                if not layers.Contains(layer):
                    continue
                for j, i in zip(*np.nonzero(inside)):
                    cells.add((li, int(j) + sy.start, int(i) + sx.start))
        elif cls == "PCB_VIA":
            p = it.GetPosition()
            radius = mm(it.GetWidth(pcbnew.F_Cu)) / 2
            pos = (mm(p.x), mm(p.y))
            win = field._window((pos[0] - radius, pos[1] - radius,
                                 pos[0] + radius, pos[1] + radius), 0.0)
            if win is None:
                continue
            sy, sx = win
            inside = np.hypot(field.X[sy, sx] - pos[0],
                              field.Y[sy, sx] - pos[1]) <= radius
            for li in range(field.nl):
                for j, i in zip(*np.nonzero(inside)):
                    cells.add((li, int(j) + sy.start, int(i) + sx.start))
        else:
            layer = it.GetLayer()
            if layer not in model.layer_index:
                continue
            li = model.layer_index[layer]
            a = (mm(it.GetStart().x), mm(it.GetStart().y))
            b = (mm(it.GetEnd().x), mm(it.GetEnd().y))
            half = mm(it.GetWidth()) / 2
            bbox = (min(a[0], b[0]), min(a[1], b[1]),
                    max(a[0], b[0]), max(a[1], b[1]))
            win = field._window(bbox, half)
            if win is None:
                continue
            sy, sx = win
            inside = _seg_distance(field.X[sy, sx], field.Y[sy, sx], a, b) <= half
            for j, i in zip(*np.nonzero(inside)):
                cells.add((li, int(j) + sy.start, int(i) + sx.start))
    return cells


NEIGHBOURS = [(-1, 0), (1, 0), (0, -1), (0, 1),
              (-1, -1), (-1, 1), (1, -1), (1, 1)]


def astar(field: Field, narrow_w: float, wide_w: float,
          starts, goals, deadline: float):
    nl, ny, nx = field.nl, field.ny, field.nx
    if not starts or not goals:
        return None
    narrow = field.blocked[narrow_w]
    wide = field.blocked[wide_w]
    via_blocked = field.via_blocked
    in_u1 = field.in_u1

    goal_set = set(goals)
    gj = np.fromiter((g[1] for g in goals), int)
    gi = np.fromiter((g[2] for g in goals), int)
    gy0, gy1 = gj.min(), gj.max()
    gx0, gx1 = gi.min(), gi.max()

    def heuristic(j, i):
        dj = 0 if gy0 <= j <= gy1 else min(abs(j - gy0), abs(j - gy1))
        di = 0 if gx0 <= i <= gx1 else min(abs(i - gx0), abs(i - gx1))
        return math.hypot(dj, di) * STEP

    INF = math.inf
    best = {}
    prev = {}
    heap = []
    for (li, j, i) in starts:
        if wide[li, j, i] and narrow[li, j, i]:
            continue
        key = (li, j, i)
        best[key] = 0.0
        heapq.heappush(heap, (heuristic(j, i), 0.0, key))

    diag = STEP * math.sqrt(2.0)
    stop = time.time() + deadline
    counter = 0
    while heap:
        counter += 1
        if (counter & 0x3FF) == 0 and time.time() > stop:
            return None
        _f, cost, key = heapq.heappop(heap)
        if cost > best.get(key, INF):
            continue
        if key in goal_set:
            path = [key]
            while path[-1] in prev:
                path.append(prev[path[-1]])
            return list(reversed(path))
        li, j, i = key
        here_narrow = in_u1[j, i]
        for dj, di in NEIGHBOURS:
            nj, ni = j + dj, i + di
            if not (0 <= nj < ny and 0 <= ni < nx):
                continue
            use_narrow = here_narrow and in_u1[nj, ni]
            grid = narrow if use_narrow else wide
            # Both ends must clear at the width this edge is emitted with.
            # Distance to each obstacle (rounded rect, disc, capsule) is a
            # convex function, so its minimum along the segment is attained at
            # an endpoint - checking the two cells is exact, not a sample.
            if grid[li, nj, ni] or grid[li, j, i]:
                continue
            step = diag if dj and di else STEP
            nkey = (li, nj, ni)
            ncost = cost + step
            if ncost < best.get(nkey, INF):
                best[nkey] = ncost
                prev[nkey] = key
                heapq.heappush(heap, (ncost + heuristic(nj, ni), ncost, nkey))
        if via_blocked[j, i]:
            continue
        for nli in range(nl):
            if nli == li:
                continue
            grid = narrow if here_narrow else wide
            if grid[nli, j, i] or grid[li, j, i]:
                continue
            nkey = (nli, j, i)
            ncost = cost + VIA_PENALTY
            if ncost < best.get(nkey, INF):
                best[nkey] = ncost
                prev[nkey] = key
                heapq.heappush(heap, (ncost + heuristic(j, i), ncost, nkey))
    return None


# --------------------------------------------------------------------------
# emission
# --------------------------------------------------------------------------
def collapse(points):
    if len(points) < 3:
        return points
    out = [points[0]]
    for k in range(1, len(points) - 1):
        ax, ay = out[-1]
        bx, by = points[k]
        cx, cy = points[k + 1]
        if abs((bx - ax) * (cy - by) - (by - ay) * (cx - bx)) < 1e-9:
            continue
        out.append(points[k])
    out.append(points[-1])
    return out


def decompose(field: Field, model: Model, path, narrow_w: float,
              wide_w: float):
    """Cell path -> [(layer, width, polyline)] plus via positions.

    The width is decided per edge, exactly as the search decided it: the
    escape-field width only where both ends of that edge are inside U1's
    courtyard.  Deciding it per run instead would emit class-width copper
    through a corridor that was only ever checked at the narrow width.
    """
    runs = []
    vias = []
    for u, v in zip(path, path[1:]):
        if u[0] != v[0]:
            vias.append((u[1], u[2]))
            continue
        width = (narrow_w if field.in_u1[u[1], u[2]] and field.in_u1[v[1], v[2]]
                 else wide_w)
        if runs and runs[-1][0] == u[0] and runs[-1][1] == width \
                and runs[-1][2][-1] == u:
            runs[-1][2].append(v)
        else:
            runs.append([u[0], width, [u, v]])
    out = []
    for li, width, cells in runs:
        pts = collapse([field.xy(j, i) for _l, j, i in cells])
        out.append((model.routable[li], width, pts))
    return out, [field.xy(j, i) for j, i in vias]


def emit(board, model: Model, field: Field, net: str, path,
         narrow_w: float, wide_w: float):
    """Turn a cell path into locked tracks and vias.  Returns (tracks, vias)."""
    netinfo = board.FindNet(net)
    geometry, vias = decompose(field, model, path, narrow_w, wide_w)
    created = []
    added_t = added_v = 0
    for layer, width, pts in geometry:
        for a, b in zip(pts, pts[1:]):
            if abs(a[0] - b[0]) < 1e-6 and abs(a[1] - b[1]) < 1e-6:
                continue
            t = pcbnew.PCB_TRACK(board)
            t.SetStart(pcbnew.VECTOR2I_MM(*a))
            t.SetEnd(pcbnew.VECTOR2I_MM(*b))
            t.SetWidth(MM(width))
            t.SetLayer(layer)
            t.SetNet(netinfo)
            t.SetLocked(True)
            board.Add(t)
            t.thisown = False
            model.tracks.append((a, b, width / 2, layer, net,
                                 t.m_Uuid.AsString(), False))
            created.append(t)
            added_t += 1
    for pos in vias:
        v = pcbnew.PCB_VIA(board)
        v.SetPosition(pcbnew.VECTOR2I_MM(*pos))
        v.SetWidth(MM(VIA_D))
        v.SetDrill(MM(VIA_DRILL))
        v.SetLayerPair(pcbnew.F_Cu, pcbnew.B_Cu)
        v.SetNet(netinfo)
        v.SetLocked(True)
        board.Add(v)
        v.thisown = False
        model.vias.append((pos, VIA_D / 2, VIA_DRILL / 2, net,
                           v.m_Uuid.AsString(), False))
        created.append(v)
        added_v += 1
    return added_t, added_v, created


def group_bbox(items):
    xs, ys = [], []
    for it in items:
        box = it.GetBoundingBox()
        xs += [mm(box.GetX()), mm(box.GetRight())]
        ys += [mm(box.GetY()), mm(box.GetBottom())]
    return min(xs), min(ys), max(xs), max(ys)


def group_centre(items):
    x0, y0, x1, y1 = group_bbox(items)
    return ((x0 + x1) / 2, (y0 + y1) / 2)


def conflicts(model: Model, net: str, geometry, vias, protected):
    """Rippable board items the proposed route cannot coexist with."""
    victims = set()
    for a, b, half, layer, other, uuid, rippable in model.tracks:
        if not rippable or uuid is None or other == net or other in protected:
            continue
        clearance = model.clearance_to(net, other)
        hit = False
        for glayer, gwidth, pts in geometry:
            if glayer != layer:
                continue
            gh = gwidth / 2
            for p, q in zip(pts, pts[1:]):
                if rlp.segments_distance(p, q, a, b) < gh + half + clearance:
                    hit = True
                    break
            if hit:
                break
        if not hit:
            for v in vias:
                if rlp.segment_distance(v, a, b) < VIA_D / 2 + half + clearance:
                    hit = True
                    break
        if hit:
            victims.add(uuid)
    for pos, radius, drill, other, uuid, rippable in model.vias:
        if not rippable or uuid is None or other == net or other in protected:
            continue
        clearance = model.clearance_to(net, other)
        hit = False
        for _glayer, gwidth, pts in geometry:
            gh = gwidth / 2
            for p, q in zip(pts, pts[1:]):
                if rlp.segment_distance(pos, p, q) < gh + radius + clearance:
                    hit = True
                    break
            if hit:
                break
        if not hit:
            for v in vias:
                gap = math.hypot(v[0] - pos[0], v[1] - pos[1])
                if gap < VIA_D / 2 + radius + clearance or \
                        gap < VIA_DRILL / 2 + drill + HOLE_TO_HOLE:
                    hit = True
                    break
        if hit:
            victims.add(uuid)
    return victims


def rip(board, model: Model, victims: set):
    """Delete the named copper, returning recipes that can rebuild it.

    Rip-up is only ever worth doing if the board ends up better for it, so
    every removal is reversible: if the nets we displaced cannot be put back,
    the whole attempt is rolled back and the board is left exactly as it was.
    """
    doomed = [t for t in board.GetTracks() if t.m_Uuid.AsString() in victims]
    recipes = []
    for item in doomed:
        if item.GetClass() == "PCB_VIA":
            p = item.GetPosition()
            recipes.append(("via", (mm(p.x), mm(p.y)),
                            mm(item.GetWidth(pcbnew.F_Cu)), mm(item.GetDrill()),
                            item.GetNetname(), item.IsLocked(), None))
        else:
            recipes.append(("track",
                            (mm(item.GetStart().x), mm(item.GetStart().y)),
                            (mm(item.GetEnd().x), mm(item.GetEnd().y)),
                            mm(item.GetWidth()), item.GetNetname(),
                            item.IsLocked(), item.GetLayer()))
        board.Remove(item)
    model.load_copper()
    return recipes


def restore(board, model: Model, recipes) -> None:
    """Rebuild copper removed by rip()."""
    for recipe in recipes:
        if recipe[0] == "via":
            _kind, pos, width, drill, net, locked, _layer = recipe
            item = pcbnew.PCB_VIA(board)
            item.SetPosition(pcbnew.VECTOR2I_MM(*pos))
            item.SetWidth(MM(width))
            item.SetDrill(MM(drill))
            item.SetLayerPair(pcbnew.F_Cu, pcbnew.B_Cu)
        else:
            _kind, start, end, width, net, locked, layer = recipe
            item = pcbnew.PCB_TRACK(board)
            item.SetStart(pcbnew.VECTOR2I_MM(*start))
            item.SetEnd(pcbnew.VECTOR2I_MM(*end))
            item.SetWidth(MM(width))
            item.SetLayer(layer)
        item.SetNet(board.FindNet(net))
        item.SetLocked(locked)
        board.Add(item)
        item.thisown = False
    model.load_copper()


def undo(board, model: Model, items) -> None:
    """Remove copper this run added."""
    for item in items:
        board.Remove(item)
    model.load_copper()


def _attempt(board, model, net, source, target, region, widths,
             narrow_w, wide_w, deadline, ignore_rippable, protected):
    field = Field(model, net, region, widths,
                  ignore_rippable=ignore_rippable, protected=protected)
    starts = seed_cells(field, model, source)
    goals = seed_cells(field, model, target)
    path = astar(field, narrow_w, wide_w, starts, goals, deadline)
    return field, path


def _regions(source, target, margins):
    sb, tb = group_bbox(source), group_bbox(target)
    base = (min(sb[0], tb[0]), min(sb[1], tb[1]),
            max(sb[2], tb[2]), max(sb[3], tb[3]))
    return [(base[0] - m, base[1] - m, base[2] + m, base[3] + m)
            for m in margins]


def route_once(board, model: Model, net: str, source, target, protected,
               margins, deadline, ignore_rippable=False):
    """One connection attempt.  Returns (field, path, region) or None."""
    width = model.route_width(net)
    widths = sorted({FIELD_WIDTH, width})
    for region in _regions(source, target, margins):
        field, path = _attempt(board, model, net, source, target, region,
                               widths, widths[0], widths[-1], deadline,
                               ignore_rippable, protected)
        if path:
            return field, path, region, widths
    return None


def close_net(board, model: Model, net: str, protected, margins, deadline,
              depth=0, max_depth=0):
    """Join every open group of `net`.  Returns (ok, created, tracks, vias)."""
    created, tracks, vias = [], 0, 0
    while True:
        board.BuildConnectivity()
        groups = net_groups(board, net)
        if len(groups) < 2:
            return True, created, tracks, vias
        groups.sort(key=len, reverse=True)
        found = route_once(board, model, net, groups[-1], groups[0],
                           protected, margins, deadline)
        if found is None:
            if depth < max_depth:
                ok, extra, t, v, _r = ripup_transaction(
                    board, model, net, groups[-1], groups[0], protected,
                    margins, deadline, depth + 1, max_depth)
                if ok:
                    created += extra
                    tracks += t
                    vias += v
                    continue
            return False, created, tracks, vias
        field, path, _region, widths = found
        t, v, items = emit(board, model, field, net, path,
                           widths[0], widths[-1])
        created += items
        tracks += t
        vias += v


def ripup_transaction(board, model: Model, net: str, source, target,
                      protected, margins, deadline, depth, max_depth):
    """Route `net` through autorouted copper, then put the displaced nets back.

    Either the whole thing succeeds or the board is restored exactly as it
    was, so a rip-up can never leave the routing worse than it found it.
    Returns (ok, created, tracks, vias, ripped).
    """
    found = route_once(board, model, net, source, target, protected,
                       margins, deadline, ignore_rippable=True)
    if found is None:
        return False, [], 0, 0, 0
    field, path, _region, widths = found
    geometry, vias_xy = decompose(field, model, path, widths[0], widths[-1])
    victims = conflicts(model, net, geometry, vias_xy, protected)
    if not victims:
        return False, [], 0, 0, 0
    victim_nets = {t[4] for t in model.tracks if t[5] in victims}
    victim_nets |= {v[3] for v in model.vias if v[4] in victims}
    victim_nets.discard(net)

    recipes = rip(board, model, victims)
    found = route_once(board, model, net, source, target, protected,
                       margins, deadline)
    if found is None:
        restore(board, model, recipes)
        return False, [], 0, 0, 0
    field, path, _region, widths = found
    tracks, vias, created = emit(board, model, field, net, path,
                                 widths[0], widths[-1])

    guard = frozenset(set(protected) | {net})
    for victim in sorted(victim_nets):
        ok, extra, vt, vv = close_net(board, model, victim, guard,
                                      margins, deadline, depth, max_depth)
        created += extra
        tracks += vt
        vias += vv
        if not ok:
            undo(board, model, created)
            restore(board, model, recipes)
            return False, [], 0, 0, 0
    return True, created, tracks, vias, len(recipes)


def connect(board, model: Model, net: str, source, target, protected,
            margins=(4.0, 9.0, 18.0), deadline=45.0, allow_ripup=True,
            max_depth=1):
    found = route_once(board, model, net, source, target, protected,
                       margins, deadline)
    if found is not None:
        field, path, _region, widths = found
        t, v, _items = emit(board, model, field, net, path,
                            widths[0], widths[-1])
        return f"OK {t} tracks {v} vias 0 ripped"
    if not allow_ripup:
        return "FAILED"
    ok, _created, t, v, ripped = ripup_transaction(
        board, model, net, source, target, protected, margins, deadline,
        0, max_depth)
    if ok:
        return f"OK {t} tracks {v} vias {ripped} ripped"
    return "FAILED"


# --------------------------------------------------------------------------
def open_nets(board):
    out = []
    for _code, info in board.GetNetsByNetcode().items():
        name = info.GetNetname()
        if not name:
            continue
        groups = net_groups(board, name)
        if len(groups) > 1:
            out.append((name, groups))
    return out


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("board", type=Path)
    ap.add_argument("--out", type=Path)
    ap.add_argument("--net", action="append",
                    help="restrict to these nets (repeatable)")
    ap.add_argument("--skip", action="append", default=[],
                    help="leave these nets alone (repeatable)")
    ap.add_argument("--deadline", type=float, default=45.0,
                    help="seconds per connection attempt")
    ap.add_argument("--passes", type=int, default=6)
    ap.add_argument("--no-ripup", action="store_true",
                    help="never rip up autorouted copper")
    ap.add_argument("--ripup-depth", type=int, default=1,
                    help="how many levels of nested rip-up to allow")
    args = ap.parse_args()
    out = args.out or args.board

    board = pcbnew.LoadBoard(str(args.board))
    pcbnew.ZONE_FILLER(board).Fill(board.Zones())
    board.BuildConnectivity()
    model = Model(board)

    started = time.time()
    total_t = total_v = total_r = 0
    # Nets this run has already settled.  They are never ripped up again, so
    # two nets cannot take turns evicting each other.
    protected: set[str] = set()
    for attempt in range(args.passes):
        # Refill first.  These pours have island removal enabled, so a fill
        # can drop a fragment that was the only thing holding two pieces of a
        # rail together; groups must be read from the board as it will
        # actually be manufactured, not from a stale fill.
        pcbnew.ZONE_FILLER(board).Fill(board.Zones())
        board.BuildConnectivity()
        pending = [name for name, _groups in open_nets(board)]
        if args.net:
            pending = [n for n in pending if n in set(args.net)]
        pending = [n for n in pending if n not in set(args.skip)]
        if not pending:
            break
        print(f"\n--- pass {attempt + 1}: {len(pending)} nets with open groups ---")
        for net in pending:
            # Groups are re-derived before every single connection.  Ripping
            # up copper deletes board items, and any Python proxy captured
            # before that points at freed memory.
            while True:
                board.BuildConnectivity()
                groups = net_groups(board, net)
                if len(groups) < 2:
                    protected.add(net)
                    break
                groups.sort(key=len, reverse=True)
                extra, main_group = groups[-1], groups[0]
                centre = group_centre(extra)
                msg = connect(board, model, net, extra, main_group, protected,
                              deadline=args.deadline,
                              allow_ripup=not args.no_ripup,
                              max_depth=args.ripup_depth)
                print(f"  {net:16} {len(extra):3d} item(s) at "
                      f"({centre[0]:6.1f},{centre[1]:6.1f}) -> {msg}")
                if not msg.startswith("OK"):
                    break
                parts = msg.split()
                total_t += int(parts[1])
                total_v += int(parts[3])
                total_r += int(parts[5])

    # Rip-up leaves tails behind when a net comes back on a different path.
    # Fill first: a track that lands on a plane is only "connected" once the
    # pour around it has been recomputed, and the dangling test would
    # otherwise delete perfectly good copper.
    pcbnew.ZONE_FILLER(board).Fill(board.Zones())
    board.BuildConnectivity()
    import finalize_routing
    dead_vias, dead_tracks = finalize_routing.remove_dangling_copper(board)
    if dead_vias or dead_tracks:
        print(f"removed {dead_tracks} dangling stubs and {dead_vias} "
              f"one-layer vias")

    pcbnew.ZONE_FILLER(board).Fill(board.Zones())
    board.BuildConnectivity()
    remaining = board.GetConnectivity().GetUnconnectedCount(True)
    print(f"\nadded {total_t} tracks and {total_v} vias, ripped {total_r} "
          f"items, in {time.time() - started:.0f}s")
    print("unconnected endpoints:", remaining)

    pcbnew.SaveBoard(str(out), board)
    for suffix in (".kicad_pro", ".kicad_dru"):
        src = str(args.board).rsplit(".", 1)[0] + suffix
        dst = str(out).rsplit(".", 1)[0] + suffix
        if os.path.exists(src) and src != dst:
            shutil.copyfile(src, dst)
    print("saved", out)


if __name__ == "__main__":
    main()