Solar MPPT Power System

#!/usr/bin/env python3
"""Clearance-correct multi-layer autorouter (rewrite of route.py).

Why a rewrite: the first router reserved occupancy as a small fixed halo in
GRID units and then tested only the single centre cell of a candidate step.
That is geometrically wrong in two ways -- it ignores the routed trace's own
width, and it quantises the clearance to the grid -- and it produced 86
different-net shorts plus 150 clearance violations whose measured "actual"
clearances clustered at 0.10-0.19mm against a 0.20mm rule (i.e. systematic
near-misses, not isolated mistakes).

The model here is exact instead:

  occ[layer]  = int32 raster, 0 = bare laminate, else the net id of the
                copper physically occupying that cell (pads, tracks, vias).
  blocked     = binary_dilation(occ != 0 and occ != me, radius = my_half_width
                + CLEARANCE)

A cell that survives that dilation is a legal position for MY centreline:
everything within (their copper + clearance + my half width) has been
excluded. Dilation runs in C via scipy, once per (net, layer), so the A*
inner loop is a single array lookup rather than a per-step disk test.

Trace width is chosen from each net's real current (derived from the
schematic's topology and the electrothermal model's own case list), not from
keyword-matching the net name -- the bug the first simulation run exposed.
Nets neck down automatically if the target width cannot be routed.
"""
import re, os, math, heapq, json, uuid as U
from collections import defaultdict

import numpy as np
from scipy import ndimage

def uid(): return str(U.uuid4())

HERE = os.path.dirname(os.path.abspath(__file__))
PCB_PATH = "/mnt/c/Users/kirkt/Documents/mppt-power-system-mcu/mppt-power-system-mcu.kicad_pcb"

GRID = 0.1                 # mm per raster cell
CLEARANCE = 0.2            # mm, different-net copper to copper (board netclass)
HOLE_CLEARANCE = 0.25      # mm, drilled hole to foreign copper (board setup)
HOLE_TO_HOLE = 0.25        # mm, edge-to-edge between two drilled holes
ROUTE_LAYERS = ["F.Cu", "In2.Cu"]        # In1.Cu and B.Cu stay solid GND planes
ESCAPE_LAYER = "B.Cu"                     # last-resort only, per link
# In1.Cu stays a solid ground plane and is NOT a routing layer. Opening it up
# for the handful of links that will not route was tried, and is recorded here
# so it is not tried again: it placed none of them -- those links do not fail
# for want of somewhere to go, they fail because they cannot get a via off
# their own pad -- and putting a wholly empty layer into the search space
# slowed ordinary routing to a standstill, because A* then has an unobstructed
# plane to wander over on every link.
LAST_RESORT_LAYER = ESCAPE_LAYER
ALL_LAYERS = ["F.Cu", "In1.Cu", "In2.Cu", "B.Cu"]
VIA_PENALTY_MM = 4.0       # cost of a layer change, in mm-equivalent
VIA_SIZE, VIA_DRILL = 0.6, 0.3
# A 0.6mm via does not fit a 0.65mm-pitch fanout: the neighbouring escape
# trace's copper edge is 0.525mm from the pad axis, leaving 0.225mm to a
# 0.6mm via's ring -- legal by 0.025mm, i.e. no manufacturing margin at all.
# Signal nets therefore fan out on a smaller via, which every prototype fab
# supports (0.4mm pad / 0.2mm drill), and current-carrying nets keep the big
# one where there is room for it.
# 0.45/0.20 gives a 0.125mm annular ring, which is exactly the minimum via a
# standard prototype process guarantees. The 0.40mm pad tried first left only
# 0.100mm of ring -- inside the fab's drill-registration tolerance, so a
# fraction of those 96 vias would have broken out of their own pad on a real
# panel. It still fits the 0.65mm fanout: a 0.225mm ring plus 0.2mm clearance
# plus the raster margin is 0.475mm against the 0.525mm a neighbouring escape
# leaves, so there is 0.05mm in hand rather than the 0.025mm a 0.6mm via has.
SIG_VIA_SIZE, SIG_VIA_DRILL = 0.45, 0.2


def via_for(net):
    """(pad, drill) for this net's vias -- sized by current, not by name."""
    return ((VIA_SIZE, VIA_DRILL) if NET_CURRENT_A.get(net, 0.0) >= 1.0
            else (SIG_VIA_SIZE, SIG_VIA_DRILL))
GND_VIA_SIZE, GND_VIA_DRILL = 0.6, 0.3

from boardsize import BOARD_W, BOARD_H   # single source of truth

# --- current per net, from the schematic's real topology -------------------
# (the same case list the electrothermal model uses; anything unlisted is a
# control/sense signal and takes SIGNAL width)
NET_CURRENT_A = {
    "Net-(J1-Pin_1)": 5.0,     # panel in -> fuse
    "Net-(D1-K)":     5.0,     # fuse -> ideal-diode FET
    "/PANEL_PROT":    5.0,     # protected panel rail
    "Net-(D1-A)":     5.0,     # ideal-diode / TVS anode side
    "Net-(Q2-D)":     5.0,     # shunt -> buck-boost bridge input
    "Net-(Q2-S)":     5.0,     # buck-side switch node
    "Net-(Q4-S)":     5.0,     # boost-side switch node
    "Net-(Q4-D)":     2.5,     # raw VBATT (boost leg out, pre-shunt)
    "/VBATT":         2.5,     # battery bus
    "Net-(L3-Pad1)":  2.0,     # 5V POL switch node
    "Net-(U8-VOS)":   2.0,     # 5V rail
    "/V5V_OUT":       2.0,     # 5V output
    "Net-(U9-SW)":    1.0,     # 3V3 POL switch node
    "Net-(C32-Pad1)": 1.0,     # 3V3 rail
    "/V3V3_OUT":      1.0,     # 3V3 output
    "Net-(U3-VI)":    0.5,     # bias rail
}
SIGNAL_W = 0.25
MIN_W = 0.15               # neck-down floor (fab min is 0.09mm)

# A raster cell is "occupied" only when its CENTRE falls inside the shape, and
# the dilation disk likewise measures centre-to-centre. Both therefore
# under-state real copper by up to half a cell diagonal, and the two errors
# add: the model's effective clearance bottoms out at
# CLEARANCE - 2*(GRID*sqrt(2)/2) = 0.2 - 0.141 = 0.059mm.
# That is not a guess -- the previous pass's tightest measured different-net
# clearance was 0.060mm, and every one of its 96 clearance violations fell
# inside that 0.141mm band. Inflating every stamp and every dilation by this
# exact amount makes the raster conservative instead of optimistic.
# Half a cell. The full half-cell-DIAGONAL bound (0.0707) is only reachable
# when both shapes meet at 45 degrees; pads, tracks and the dilation disk are
# overwhelmingly axis-aligned, where the bound is half a cell. 0.0707 was
# tried first and cost 19 of 28 escapes on U4 -- the fanout corridor in a
# 0.65mm pin field is narrower than that margin. Residual raster error is
# caught exactly afterwards by the rip-up pass, which is what actually
# guarantees the board, so this number only has to be a good default.
RASTER_MARGIN = GRID / 2                   # 0.05mm

# Fanout comb geometry. Swept empirically against the number of fine-pitch
# escapes that actually route -- see the sweep in this module's history.
STUB_BASE = 0.5            # mm, shortest tooth
STUB_STEP = 0.55           # mm, extra length per step toward the row's middle
STUB_WIDTHS = (SIGNAL_W, 0.2, MIN_W)

def vs_drill_half(via, net):
    """Half the drill diameter of the via this net will use."""
    return (via if via else via_for(net))[1] / 2


def width_for(net):
    a = NET_CURRENT_A.get(net, 0.0)
    if a >= 5.0:  return 1.2
    if a >= 2.5:  return 0.8
    if a >= 2.0:  return 0.6
    if a >= 1.0:  return 0.4
    if a >= 0.5:  return 0.3
    return SIGNAL_W


# --------------------------------------------------------------------------
# board parsing
# --------------------------------------------------------------------------
def split_blocks(text, tag):
    out, i, n, needle = [], 0, len(text), '(' + tag
    while True:
        i = text.find(needle, i)
        if i == -1: break
        j = i + len(needle)
        if j < n and (text[j].isalnum() or text[j] == '_'):
            i = j; continue
        depth = 0; k = i; ins = False; esc = False
        while k < n:
            ch = text[k]
            if esc: esc = False
            elif ch == '\\': esc = True
            elif ch == '"': ins = not ins
            elif not ins:
                if ch == '(': depth += 1
                elif ch == ')':
                    depth -= 1
                    if depth == 0: k += 1; break
            k += 1
        out.append(text[i:k]); i = k
    return out


def parse_pads(text):
    """Every pad, in world coordinates, with its real size and layer set."""
    pads = []
    for blk in split_blocks(text, 'footprint'):
        m_ref = re.search(r'\(property "Reference" "([^"]*)"', blk)
        ref = m_ref.group(1) if m_ref else "?"
        m_at = re.search(
            r'\(layer "[^"]*"\)\s*(?:\([^()]*\)\s*)*?\(uuid "[^"]*"\)'
            r'\s*(?:\([^()]*\)\s*)*?\(at ([-\d.]+) ([-\d.]+)(?: ([-\d.]+))?\)', blk)
        if not m_at: continue
        fx, fy = float(m_at.group(1)), float(m_at.group(2))
        th = math.radians(float(m_at.group(3) or 0))
        ct, st = math.cos(th), math.sin(th)
        for pblk in split_blocks(blk, 'pad'):
            mp = re.match(r'\(pad "(\w*)"\s+(\w+)\s+(\w+)', pblk)
            if not mp or not mp.group(1): continue
            num = mp.group(1)
            m_lat = re.search(r'\(at ([-\d.]+) ([-\d.]+)(?: ([-\d.]+))?\)', pblk)
            if not m_lat: continue
            lx, ly = float(m_lat.group(1)), float(m_lat.group(2))
            pad_rot = float(m_lat.group(3) or 0)
            m_sz = re.search(r'\(size ([-\d.]+) ([-\d.]+)\)', pblk)
            pw, ph = (float(m_sz.group(1)), float(m_sz.group(2))) if m_sz else (0.5, 0.5)
            narrow = min(pw, ph)
            m_ly = re.search(r'\(layers ([^)]*)\)', pblk)
            lys = m_ly.group(1) if m_ly else ''
            any_layer = '*.Cu' in lys
            # the layers this pad's copper actually occupies. Recording only
            # "is it *.Cu" made every B.Cu-only pad invisible to the router:
            # U4's HTSSOP exposed pad has a bottom-side thermal pad, and the
            # B.Cu escape routes ran straight through it -- 28 real overlaps
            # that no raster check could see, because the copper was never
            # stamped in the first place.
            pad_layers = (tuple(ALL_LAYERS) if any_layer else
                          tuple(L for L in ALL_LAYERS if f'"{L}"' in lys))
            m_dr = re.search(r'\(drill ([\d.]+)\)', pblk)
            drill = float(m_dr.group(1)) if m_dr else 0.0
            m_net = re.search(r'\(net (?:\d+ )?"([^"]*)"\)', pblk)
            net = m_net.group(1) if m_net else ""
            pads.append({
                "ref": ref, "num": num, "net": net,
                "x": fx + lx * ct + ly * st,
                "y": fy + ly * ct - lx * st,
                "any_layer": any_layer,
                "layers": pad_layers or ("F.Cu",),
                "drill": drill,
                # TRUE pad geometry. An earlier version approximated every pad
                # by a circle of radius hypot(w,h)/2; for U4's 1.55x0.4mm
                # HTSSOP pads that is a 1.60mm envelope against a 0.65mm pitch,
                # so neighbouring pads' envelopes overlapped by ~0.95mm and
                # buried the entire pin field. Every escape from U4 was then
                # geometrically "impossible" and 19 of 27 nets failed to route.
                "w": pw, "h": ph,
                # ABSOLUTE, not footprint + pad: KiCad orients a pad's shape by the
                # pad's own angle alone (see geomcheck for the test that pins this)
                "angle": pad_rot,
                "radius": math.hypot(pw, ph) / 2,   # only for coarse distance tests
                # narrow dimension: a trace wider than this cannot land on the
                # pad without spilling over its neighbours on a fine-pitch part
                "narrow": narrow,
            })
    return pads


# --------------------------------------------------------------------------
# raster
# --------------------------------------------------------------------------
NX = int(BOARD_W / GRID) + 2
NY = int(BOARD_H / GRID) + 2

def to_cell(x, y):
    return int(round(x / GRID)), int(round(y / GRID))

def disk(radius_mm):
    r = max(1, int(math.ceil(radius_mm / GRID)))
    yy, xx = np.ogrid[-r:r + 1, -r:r + 1]
    return (xx * xx + yy * yy) <= r * r

def stamp(arr, cx, cy, radius_mm, value):
    """Paint a filled disk of net id `value` into a raster (conservatively)."""
    radius_mm += RASTER_MARGIN
    r = max(1, int(math.ceil(radius_mm / GRID)))
    x0, x1 = max(0, cx - r), min(NX, cx + r + 1)
    y0, y1 = max(0, cy - r), min(NY, cy + r + 1)
    if x1 <= x0 or y1 <= y0: return
    d = disk(radius_mm)
    dy0 = y0 - (cy - r); dx0 = x0 - (cx - r)
    sub = d[dy0:dy0 + (y1 - y0), dx0:dx0 + (x1 - x0)]
    region = arr[y0:y1, x0:x1]
    region[sub] = value

def stamp_rect(arr, x, y, w, h, angle_deg, value):
    """Paint an oriented rectangle (a real SMD pad) into a raster, exactly.

    Pads are rectangles, and on fine-pitch parts the difference between the
    rectangle and its circumscribed circle is the difference between a
    routable and an unroutable pin field."""
    hw, hh = w / 2 + RASTER_MARGIN, h / 2 + RASTER_MARGIN
    reach = math.hypot(hw, hh)
    r = int(math.ceil(reach / GRID)) + 1
    cx, cy = to_cell(x, y)
    x0, x1 = max(0, cx - r), min(NX, cx + r + 1)
    y0, y1 = max(0, cy - r), min(NY, cy + r + 1)
    if x1 <= x0 or y1 <= y0: return
    ys, xs = np.mgrid[y0:y1, x0:x1]
    dx = xs * GRID - x
    dy = ys * GRID - y
    th = math.radians(angle_deg)
    ct, st = math.cos(th), math.sin(th)
    # rotate the sample points into the pad's own frame
    rx = dx * ct + dy * st
    ry = -dx * st + dy * ct
    inside = (np.abs(rx) <= hw) & (np.abs(ry) <= hh)
    region = arr[y0:y1, x0:x1]
    region[inside] = value


def stamp_segment(arr, c0, c1, radius_mm, value):
    x0, y0 = c0; x1, y1 = c1
    n = max(abs(x1 - x0), abs(y1 - y0), 1)
    for i in range(n + 1):
        stamp(arr, int(round(x0 + (x1 - x0) * i / n)),
                   int(round(y0 + (y1 - y0) * i / n)), radius_mm, value)


def along_of(p, ux, uy):
    """Coordinate of a pad along its own pin row."""
    return -p["x"] * uy + p["y"] * ux


def turn_sign(p, ux, uy, rowmid):
    """+1/-1: which way along the row is the nearer way out."""
    mid = rowmid.get((p["ref"], round(ux), round(uy)), 0.0)
    return 1.0 if along_of(p, ux, uy) >= mid else -1.0



def hard_link_pass(R, hard, rounds=4, radius=4.0):
    """Route the links that had nowhere to go, by moving what is in the way.

    Every link here failed because its endpoints are walled in -- almost always
    a pin on a fine-pitch package whose escape corridor was taken by a net that
    happened to be routed earlier. No amount of re-trying the same link against
    the same board changes that; something has to move.

    So: take the copper belonging to OTHER links that sits near the blocked
    pads, rip it up, route the hard link into the space that frees, then route
    the ripped links again -- they are ordinary links with the whole rest of
    the board still available to them, so they nearly always come back. If the
    hard link still cannot be placed, the ripped links are restored and nothing
    is lost. This is the rip-up-and-reroute every commercial router does, scoped
    to the handful of links that actually need it.
    """
    fixed, still = 0, []
    for rnd in range(1, rounds + 1):
        if not hard:
            break
        # A fixed eviction radius, deliberately. Widening it each round was
        # tried: it evicts more, rolls back more, and measured worse -- 4 links
        # left unrouted against 3. The links that survive this pass are not
        # short of nearby space; they cannot get a via off their own pad at
        # all, so evicting copper further away buys nothing.
        radius_now = radius
        print(f"\nhard-link pass {rnd}: {len(hard)} link(s) with no route "
              f"(eviction radius {radius_now:.0f}mm)", flush=True)
        progressed = False
        remaining = []
        for (name, a, b, eff_w) in hard:
            me = R.nid(name)
            # which committed links own copper near either blocked pad?
            blockers = set()
            for li, (link, s0, s1, v0, v1) in enumerate(R.link_span):
                if link[0] == name:
                    continue
                for i in range(s0, s1):
                    (x0, y0), (x1, y1), _w, _L, _n = R.segments[i]
                    for (px, py) in ((a["x"], a["y"]), (b["x"], b["y"])):
                        if (min(math.hypot(x0 - px, y0 - py),
                                math.hypot(x1 - px, y1 - py)) < radius_now):
                            blockers.add(li); break
                    if li in blockers: break
                if li in blockers: continue
                for i in range(v0, v1):
                    x, y, _n, _s, _d = R.vias[i]
                    if any(math.hypot(x - px, y - py) < radius_now
                           for (px, py) in ((a["x"], a["y"]), (b["x"], b["y"]))):
                        blockers.add(li); break
            # Before disturbing anything, try the whole ladder including the
            # plane layer. An eviction that is not needed is pure risk: the
            # evicted links may not come back, and then a perfectly good
            # attempt gets rolled back for nothing.
            ALL3 = ROUTE_LAYERS + [ESCAPE_LAYER]
            free_try = w0 = None
            for w0, lys, mm in [(SIGNAL_W, ALL3, 40.0),
                                (MIN_W, ALL3, 90.0)]:
                if w0 > eff_w: continue
                free_try = R.route(name, a, b, w0, lys, margin_mm=mm)
                if free_try: break
            if free_try:
                lys0 = (ROUTE_LAYERS + [ESCAPE_LAYER]
                        if max(q[0] for q in free_try) >= len(ROUTE_LAYERS)
                        else ROUTE_LAYERS)
                s0, v0 = len(R.segments), len(R.vias)
                R.commit(name, free_try, w0, lys0, a, b)
                R.link_span.append(((name, a, b, eff_w), s0, len(R.segments),
                                    v0, len(R.vias)))
                fixed += 1; progressed = True
                print(f"   {name}: placed on the plane layer, nothing evicted",
                      flush=True)
                continue

            if not blockers:
                remaining.append((name, a, b, eff_w)); continue
            blockers = sorted(blockers)
            snap = R.snapshot_links(blockers)
            R.rip(blockers)

            got = w = None
            for w, lys, mm in [(eff_w, ALL3, 25.0),
                               (SIGNAL_W, ALL3, 25.0),
                               (MIN_W, ALL3, 60.0)]:
                if w > eff_w: continue
                got = R.route(name, a, b, w, lys, margin_mm=mm)
                if got: break

            if not got:
                R.restore_links(snap)          # exact undo, nothing lost
                remaining.append((name, a, b, eff_w))
                continue

            lys2 = (ROUTE_LAYERS + [ESCAPE_LAYER]
                    if max(q[0] for q in got) >= len(ROUTE_LAYERS)
                    else ROUTE_LAYERS)
            s0, v0 = len(R.segments), len(R.vias)
            R.commit(name, got, w, lys2, a, b)
            hard_idx = len(R.link_span)
            R.link_span.append(((name, a, b, eff_w), s0, len(R.segments),
                                v0, len(R.vias)))

            # the hard link is placed; now see whether the evicted ones fit
            # around it
            back, lost = 0, []
            for (en, ea, eb, ew) in [sn[0] for sn in snap]:
                p2 = w2 = None
                for w2, lys, mm in [(ew, ROUTE_LAYERS + [ESCAPE_LAYER], 25.0),
                                    (SIGNAL_W, ROUTE_LAYERS + [ESCAPE_LAYER], 40.0),
                                    (MIN_W, ROUTE_LAYERS + [ESCAPE_LAYER], 90.0)]:
                    if w2 > ew: continue
                    p2 = R.route(en, ea, eb, w2, lys, margin_mm=mm)
                    if p2: break
                if p2:
                    lys3 = (ROUTE_LAYERS + [ESCAPE_LAYER]
                            if max(q[0] for q in p2) >= len(ROUTE_LAYERS)
                            else ROUTE_LAYERS)
                    t0, u0 = len(R.segments), len(R.vias)
                    R.commit(en, p2, w2, lys3, ea, eb)
                    R.link_span.append(((en, ea, eb, ew), t0, len(R.segments),
                                        u0, len(R.vias)))
                    back += 1
                else:
                    lost.append((en, ea, eb, ew))

            if lost:
                # Placing one link by losing one or more is not progress. Undo
                # the whole attempt exactly and move on.
                keep = [i for i, (lk, *_r) in enumerate(R.link_span)
                        if i >= hard_idx]
                R.rip(keep)
                R.restore_links(snap)
                remaining.append((name, a, b, eff_w))
                print(f"   {name}: would have cost {len(lost)} link(s) -- "
                      f"attempt rolled back", flush=True)
            else:
                fixed += 1
                progressed = True

        hard = remaining
        print(f"   routed {fixed} so far; {len(hard)} still unrouted",
              flush=True)
        if not progressed:
            break
    return fixed, hard


def pad_key(p):
    return (p["ref"], p["num"], round(p["x"], 3), round(p["y"], 3))


def union_find_pads(R):
    """Group every pad into the piece of copper it is actually joined to.

    Two pads are joined if a committed link runs between them, or if their
    copper physically touches (multi-pad pins, tab pads).
    """
    parent = {}

    def find(k):
        parent.setdefault(k, k)
        while parent[k] != k:
            parent[k] = parent[parent[k]]; k = parent[k]
        return k

    def union(x, y):
        rx, ry = find(x), find(y)
        if rx != ry: parent[rx] = ry

    for p in R.pads:
        find(pad_key(p))
    for (link, _s0, _s1, _v0, _v1) in R.link_span:
        _n, a, b, _w = link
        union(pad_key(a), pad_key(b))
    by_net = defaultdict(list)
    for p in R.pads:
        if p["net"]: by_net[p["net"]].append(p)
    for ps in by_net.values():
        for i, a in enumerate(ps):
            for b in ps[i + 1:]:
                if (abs(a["x"] - b["x"]) < 6 and abs(a["y"] - b["y"]) < 6
                        and pads_abut(a, b)):
                    union(pad_key(a), pad_key(b))
    return {pad_key(p): find(pad_key(p)) for p in R.pads}


def pads_abut(a, b):
    """Do these two pads' real copper shapes actually meet?"""
    import geomcheck as G
    ra = G.ORect(a["x"], a["y"], a["w"], a["h"], a["angle"])
    oa = G.Obj("pad", a["net"], ("F.Cu",), ax=a["x"], ay=a["y"],
               bx=a["x"], by=a["y"], rect=ra)
    ob = G.Obj("pad", b["net"], ("F.Cu",), ax=b["x"], ay=b["y"], bx=b["x"],
               by=b["y"], rect=G.ORect(b["x"], b["y"], b["w"], b["h"], b["angle"]))
    return oa.copper_gap(ob) <= 0.0


class Router:
    def __init__(self, pads):
        self.pads = pads
        self.net_id = {"": 0}
        self.id_net = {0: ""}
        for p in pads:
            if p["net"] and p["net"] not in self.net_id:
                i = len(self.net_id)
                self.net_id[p["net"]] = i
                self.id_net[i] = p["net"]
        # occupancy raster per routable layer (+B.Cu escape)
        self.layers = ROUTE_LAYERS + [ESCAPE_LAYER]
        self.occ = {L: np.zeros((NY, NX), dtype=np.int32) for L in self.layers}
        # every drilled hole, regardless of net: the hole-to-hole rule does not
        # care whether two vias are on the same net, and 13 findings said so
        self.holes = np.zeros((NY, NX), dtype=bool)
        self.segments = []
        self.vias = []
        self.via_audit = []
        # (link, first_seg, last_seg, first_via, last_via) so a single link's
        # copper can be removed again without disturbing anything else
        self.link_span = []
        self.stub_segs = 0
        self.used_anchors = set()
        # which fanout tip each fanout segment belongs to: an L has two legs
        # and they must live or die together
        self.stub_owner = []

    def nid(self, net):
        if net not in self.net_id:
            i = len(self.net_id)
            self.net_id[net] = i; self.id_net[i] = net
        return self.net_id[net]

    def stamp_pads(self):
        for p in self.pads:
            v = self.nid(p["net"]) if p["net"] else -1   # -1 = unnetted copper
            for L in self.layers:
                if L in p["layers"]:
                    stamp_rect(self.occ[L], p["x"], p["y"],
                               p["w"], p["h"], p["angle"], v)
            if p.get("drill"):      # a through-hole pad drills a hole too
                stamp(self.holes, *to_cell(p["x"], p["y"]), p["drill"] / 2, True)


    def escape_stubs(self):
        """Give every crowded pad an exact-geometry escape trace.

        The raster carries a deliberate 0.141mm safety margin (RASTER_MARGIN),
        which is the right call in open board area and fatal at a fine-pitch
        pin field. On U4 (HTSSOP, 0.65mm pitch, 0.4mm pads) the legal escape
        corridor down the middle of a pad is about 0.1mm wider than the
        centreline needs -- real, but finer than the 0.1mm grid can express, so
        rounding the pad centre to a cell closes it. That is why 19 of U4's
        28 pins could not be escaped and the routed-net count fell from 50 to
        19 when the margin was introduced.

        So the pad escape is not done on the raster at all. Each crowded pad
        gets a short trace straight out along its own long axis, sized and
        verified with exact geometry (geomcheck, no grid anywhere). A* then
        starts from the TIP of that stub, out in open copper where the raster's
        conservatism costs nothing. This is the fanout pass a commercial
        autorouter runs before it routes, and for the same reason.
        """
        import geomcheck as G
        # exact geometry of every pad, for the clearance test
        gp = [G.Obj("pad", p["net"], ("F.Cu",), tag=f'{p["ref"]}.{p["num"]}',
                    ax=p["x"], ay=p["y"], bx=p["x"], by=p["y"],
                    rect=G.ORect(p["x"], p["y"], p["w"], p["h"], p["angle"]))
              for p in self.pads]
        idx = G.Index(gp, bucket=2.0)

        centre = {}
        for p in self.pads:
            centre.setdefault(p["ref"], []).append((p["x"], p["y"]))
        centre = {r: (sum(a for a, _ in v) / len(v), sum(b for _, b in v) / len(v))
                  for r, v in centre.items()}

        # position of each pad within its own pin row, counted from whichever
        # end of the row is nearer -- the comb's tooth length
        rows = {}
        for p in self.pads:
            if (not p["net"] or p["any_layer"] or "F.Cu" not in p["layers"]
                    or p["net"].startswith("unconnected-")):
                continue
            th = math.radians(p["angle"])
            if p["w"] >= p["h"]: ux, uy = math.cos(th), math.sin(th)
            else:                ux, uy = -math.sin(th), math.cos(th)
            cx, cy = centre.get(p["ref"], (p["x"], p["y"]))
            if (p["x"] - cx) * ux + (p["y"] - cy) * uy < 0:
                ux, uy = -ux, -uy
            side = (p["ref"], round(ux), round(uy))
            # distance along the row (perpendicular to the escape direction)
            along = -p["x"] * uy + p["y"] * ux
            rows.setdefault(side, []).append((along, p))
        rowpos, rowspan, rowmid, rowpitch = {}, {}, {}, {}
        for side, items in rows.items():
            items.sort(key=lambda t: t[0])       # sort BEFORE measuring the row
            rowspan[side] = items[-1][0] - items[0][0] if len(items) > 1 else 0.0
            rowpitch[side] = (rowspan[side] / (len(items) - 1)
                              if len(items) > 1 else 9.9)
            rowmid[side] = (items[0][0] + items[-1][0]) / 2 if items else 0.0
        for side, items in rows.items():
            n = len(items)
            for i, (_a, p) in enumerate(items):
                # Distance from the NEARER end of the row, so the teeth get
                # longer toward the middle where an escape has furthest to go
                # before it reaches open board. A plain 3-cycle was tried
                # instead, on the theory that it guarantees no two adjacent
                # teeth end level; it measured worse on the whole board
                # (44 of 60 nets routed against 46, 20 failed links against
                # 16), because keeping the teeth short matters less than
                # getting the middle of a long row out past its neighbours.
                # Distance from the nearer end of the row. Adding a half-step
                # for the far half was tried, to stop the two pins either side
                # of a row's centre turning at the same x and sealing whatever
                # is between them -- it fixes that, and costs five other links
                # their route, because it perturbs every fanout on the board.
                # The one place it mattered is handled locally instead.
                rowpos[(p["ref"], p["num"],
                        round(p["x"], 3), round(p["y"], 3))] = min(i, n - 1 - i)

        made = 0
        for p, g in zip(self.pads, gp):
            p["ax"], p["ay"] = p["x"], p["y"]        # default: the pad centre
            p["crowded"] = False
            if not p["net"] or p["any_layer"] or p["net"] == "GND":
                continue    # GND reaches the pour on its own layer
            if "F.Cu" not in p["layers"]:
                continue    # a bottom-side-only pad has no F.Cu fanout
            if p["net"].startswith("unconnected-"):
                # A no-connect pin has nothing to route, so its fanout is pure
                # obstacle: U9's NC pin was handed a 1.9mm L-shaped stub that
                # sat directly across the escape of U9.4, the 3.3V buck's
                # switch node, and helped make that link unroutable.
                continue
            # crowded? -- a different-net pad close enough to squeeze the escape
            if idx.min_gap(g, reach=1.6) > 0.9:
                continue
            p["crowded"] = True
            long_dim, short_dim = max(p["w"], p["h"]), min(p["w"], p["h"])
            th = math.radians(p["angle"])
            # the pad's long axis, in world coordinates
            if p["w"] >= p["h"]:
                ux, uy = math.cos(th), math.sin(th)
            else:
                ux, uy = -math.sin(th), math.cos(th)
            cx, cy = centre.get(p["ref"], (p["x"], p["y"]))
            if (p["x"] - cx) * ux + (p["y"] - cy) * uy < 0:
                ux, uy = -ux, -uy                     # point away from the part
            # Staggered comb, not a flat wall. Equal-length stubs put every tip
            # of a pin row on one line, and at a 0.5mm pitch (U2) that line has
            # no gap in it at all: neighbouring 0.25mm escapes touch, and no
            # via fits between them, so every pin on that side was trapped.
            # Making the stub longer toward the middle of the row means each
            # trace reaches open board where its neighbour has already
            # stopped, and can turn or drop a via there. This is the standard
            # fanout comb.
            j = rowpos.get((p["ref"], p["num"], round(p["x"], 3), round(p["y"], 3)), 0)
            base = STUB_BASE + STUB_STEP * j
            # ---- L-shaped fanout -------------------------------------
            # A radial stub gets the trace off its pad and no further: at the
            # tip it is still inside the pin field, flanked by its neighbours'
            # stubs, and on a 0.5mm-pitch part there is room neither to turn
            # nor to drop a via. Every one of U2's six failing pins failed
            # exactly there, with the A* start cell walled in on F.Cu and both
            # inner layers wide open underneath it.
            #
            # So the fanout is an L, which is what it is on a real board: out
            # along the pad's own axis far enough to clear the package, then a
            # turn along the row and out past the end of it. Because the
            # out-leg gets longer toward the middle of the row (base, above),
            # the turns happen at staggered radii and the along-legs stack in
            # parallel lanes STUB_STEP apart instead of colliding. A* then
            # starts from a point outside the package footprint entirely.
            side_key = (p["ref"], round(ux), round(uy))
            row_len = rowspan.get(side_key, 0.0)
            # Turn only where a straight escape cannot finish. A via needs
            # VIA/2 + CLEARANCE + margin = 0.45mm of room beside it, and at a
            # pin pitch P a neighbouring 0.25mm escape leaves P - 0.125mm. So
            # below ~0.575mm pitch the trace can never via at its own tip and
            # must be carried out of the pin field; above it, it can, and the
            # L only wastes board. Measured both ways: forcing the L on U4
            # (0.65mm pitch, 14 pins a side) stacks 14 parallel along-legs
            # through the very area its routes need, and took its escapes from
            # 4 of 5 down to 1 of 5. U2/U8/U9 are all 0.5mm and need it.
            pitch = rowpitch.get(side_key, 9.9)
            if pitch >= 0.575:
                row_len = 0.0
            done = False
            for extra in [e for e in ([base + 0.4]
                                      + [base - STUB_STEP * k for k in range(int(j) + 1)]
                                      + [0.5]) if e >= 0.45]:
                cx1 = p["x"] + ux * (long_dim / 2 + extra)
                cy1 = p["y"] + uy * (long_dim / 2 + extra)
                for w in STUB_WIDTHS:
                    if w > short_dim:
                        continue
                    legs = [((p["x"], p["y"]), (cx1, cy1))]
                    # the turn: along the row, away from its centre, far
                    # enough to leave the package behind
                    if row_len > 0.6:
                        run = row_len / 2 - abs(along_of(p, ux, uy)
                                                - rowmid.get((p["ref"], round(ux),
                                                              round(uy)), 0.0)) + 1.4
                        if run > 0.6:
                            vx, vy = -uy * turn_sign(p, ux, uy, rowmid), \
                                     ux * turn_sign(p, ux, uy, rowmid)
                            legs.append(((cx1, cy1),
                                         (cx1 + vx * run, cy1 + vy * run)))
                    ok = True
                    for (ax_, ay_), (bx_, by_) in legs:
                        cap = G.Obj("track", p["net"], ("F.Cu",), ax=ax_, ay=ay_,
                                    bx=bx_, by=by_, r=w / 2)
                        if idx.min_gap(cap, reach=2.0) < CLEARANCE + RASTER_MARGIN:
                            ok = False; break
                    if not ok:
                        # the turn is blocked -- keep the straight leg alone
                        legs = legs[:1]
                        cap = G.Obj("track", p["net"], ("F.Cu",), ax=p["x"],
                                    ay=p["y"], bx=cx1, by=cy1, r=w / 2)
                        if idx.min_gap(cap, reach=2.0) < CLEARANCE + RASTER_MARGIN:
                            continue
                    # a leg can come out zero-length when the turn lands on
                    # the corner it started from; KiCad calls those dangling
                    legs = [(q, r) for q, r in legs
                            if math.hypot(r[0] - q[0], r[1] - q[1]) > 1e-6]
                    if not legs:
                        continue
                    tip = (round(legs[-1][1][0], 3), round(legs[-1][1][1], 3))
                    for (ax_, ay_), (bx_, by_) in legs:
                        self.segments.append(((ax_, ay_), (bx_, by_), w,
                                              "F.Cu", p["net"]))
                        self.stub_owner.append(tip)
                        stamp_segment(self.occ["F.Cu"], to_cell(ax_, ay_),
                                      to_cell(bx_, by_), w / 2, self.nid(p["net"]))
                    p["ax"], p["ay"] = legs[-1][1]
                    made += 1
                    done = True
                    break
                if done:
                    break
        return made


    def as_geometry(self):
        """This router's copper, as exact shapes for geomcheck."""
        import geomcheck as G
        objs = []
        for p in self.pads:
            lys = p["layers"]
            objs.append(G.Obj("pad", p["net"], lys,
                              tag=f'{p["ref"]}.{p["num"]}',
                              ax=p["x"], ay=p["y"], bx=p["x"], by=p["y"],
                              rect=G.ORect(p["x"], p["y"], p["w"], p["h"],
                                           p["angle"])))
        seg_obj, via_obj = {}, {}
        for i, ((x0, y0), (x1, y1), w, L, net) in enumerate(self.segments):
            o = G.Obj("track", net, (L,), ax=x0, ay=y0, bx=x1, by=y1, r=w / 2)
            seg_obj[id(o)] = i; objs.append(o)
        for i, (x, y, net, size, drill) in enumerate(self.vias):
            o = G.Obj("via", net, G.ALL_CU, ax=x, ay=y, bx=x, by=y,
                      r=size / 2, hole=drill / 2)
            via_obj[id(o)] = i; objs.append(o)
        return objs, seg_obj, via_obj

    def offending_links(self):
        """Links whose copper breaks a real clearance rule, worst first.

        The raster is an approximation and always will be. This is the check
        that decides whether the board is actually legal: exact distances
        between exact shapes, no grid involved.
        """
        import geomcheck as G
        objs, seg_obj, via_obj = self.as_geometry()
        owner = {}
        for li, (link, s0, s1, v0, v1) in enumerate(self.link_span):
            for i in range(s0, s1): owner[("s", i)] = li
            for i in range(v0, v1): owner[("v", i)] = li
        bad = {}; unowned = 0
        for gap, rule, why, a, b in G.violations(objs):
            for o in (a, b):
                key = (("s", seg_obj[id(o)]) if id(o) in seg_obj else
                       ("v", via_obj[id(o)]) if id(o) in via_obj else None)
                if key is None: continue
                li = owner.get(key)
                if li is None:
                    unowned += 1              # a fanout stub or a GND via
                    continue
                bad[li] = min(bad.get(li, 1e9), gap - rule)
        if unowned:
            # These are real violations the rip-up loop cannot act on, because
            # the copper belongs to no rippable link. Saying so is the point:
            # silently skipping them is how 28 overlaps with U4's thermal pad
            # survived a "clean" exact-geometry check.
            print(f"   note: {unowned} violation(s) involve copper owned by no "
                  f"link (fanout stub or stitching via) -- not rippable",
                  flush=True)
        return sorted(bad.items(), key=lambda kv: kv[1])

    def rip(self, link_indices):
        """Remove those links' copper and rebuild the occupancy raster."""
        drop = set(link_indices)
        keep_s, keep_v, new_span = [], [], []
        # escape stubs live at the head of self.segments and are never ripped
        stub = self.segments[:self.stub_segs]
        for li, (link, s0, s1, v0, v1) in enumerate(self.link_span):
            if li in drop: continue
            ns0 = len(stub) + len(keep_s); nv0 = len(keep_v)
            keep_s.extend(self.segments[s0:s1])
            keep_v.extend(self.vias[v0:v1])
            new_span.append((link, ns0, len(stub) + len(keep_s),
                             nv0, len(keep_v)))
        self.segments = stub + keep_s
        self.vias = keep_v
        self.link_span = new_span
        # A ripped-up link's fanout is no longer attached to anything, so the
        # anchor set has to be rebuilt from the links that actually survive --
        # otherwise the pruner keeps its fanout and leaves it dangling.
        self.used_anchors = set()
        for (link, _s0, _s1, _v0, _v1) in self.link_span:
            _n, a, b, _w = link
            for e in (a, b):
                self.used_anchors.add((round(e["ax"], 3), round(e["ay"], 3)))
        for L in self.layers:
            self.occ[L][:] = 0
        self.holes[:] = False       # BEFORE stamp_pads -- it stamps PTH drills
        self.stamp_pads()
        for (x0, y0), (x1, y1), w, L, net in self.segments:
            stamp_segment(self.occ[L], to_cell(x0, y0), to_cell(x1, y1),
                          w / 2, self.nid(net))
        for x, y, net, size, drill in self.vias:
            cx, cy = to_cell(x, y)
            for L in self.layers:
                stamp(self.occ[L], cx, cy, size / 2, self.nid(net))
            stamp(self.holes, cx, cy, drill / 2, True)
        return len(drop)

    def snapshot_links(self, indices):
        """Copy those links' copper verbatim, so it can be put back exactly.

        Re-ROUTING an evicted link is not a restore: by the time it is evicted
        the board is far fuller than when it was first placed, so the router
        often cannot find the path it already had, and an "undo" that loses
        five links is worse than never having tried. Keeping the geometry means
        undo is exact and free.
        """
        snap = []
        for li in indices:
            link, s0, s1, v0, v1 = self.link_span[li]
            snap.append((link, list(self.segments[s0:s1]), list(self.vias[v0:v1])))
        return snap

    def restore_links(self, snap):
        """Put snapshotted copper back, exactly where it was."""
        for link, segs, vias in snap:
            s0, v0 = len(self.segments), len(self.vias)
            for sg in segs:
                (x0, y0), (x1, y1), w, L, net = sg
                self.segments.append(sg)
                stamp_segment(self.occ[L], to_cell(x0, y0), to_cell(x1, y1),
                              w / 2, self.nid(net))
            for vv in vias:
                x, y, net, size, drill = vv
                self.vias.append(vv)
                cx, cy = to_cell(x, y)
                for L in self.layers:
                    stamp(self.occ[L], cx, cy, size / 2, self.nid(net))
                stamp(self.holes, cx, cy, drill / 2, True)
            self.link_span.append((link, s0, len(self.segments),
                                   v0, len(self.vias)))
            # Re-arm the anchors. rip() rebuilds used_anchors from the links
            # that survive it, so a restored link's fanout tips are no longer
            # in that set -- and prune_unused_stubs then deletes the fanout
            # this link is standing on, silently disconnecting the pad. That
            # is how a route reporting 3 unroutable links produced a board
            # with 10 split nets.
            _n, a, b, _w = link
            for e in (a, b):
                self.used_anchors.add((round(e["ax"], 3), round(e["ay"], 3)))
        return len(snap)

    def prune_unused_stubs(self):
        """Throw away escape stubs that no route ever attached to.

        A stub whose link failed is a track with a free end: real copper going
        nowhere. KiCad calls it track_dangling and it is right to -- it is
        an antenna on a switching board and it confuses the connectivity
        check. The stubs are speculative by design, so the unused ones are
        simply removed once routing is finished.
        """
        stubs = self.segments[:self.stub_segs]
        rest = self.segments[self.stub_segs:]
        # Keyed on the fanout's TIP, not on each segment's own end point. An
        # L-shaped fanout has an out-leg whose end is the corner, which is
        # never an anchor -- matching per segment deleted the out-leg while
        # keeping the along-leg, leaving 17 floating tracks and taking the
        # unconnected count from 17 to 34.
        keep_i = [i for i, sg in enumerate(stubs)
                  if self.stub_owner[i] in self.used_anchors]
        keep = [stubs[i] for i in keep_i]
        dropped = len(stubs) - len(keep)
        if dropped:
            self.segments = keep + rest
            self.stub_owner = [self.stub_owner[i] for i in keep_i]
            self.stub_segs = len(keep)
            for L in self.layers:
                self.occ[L][:] = 0
            self.holes[:] = False   # BEFORE stamp_pads -- it stamps PTH drills
            self.stamp_pads()
            for (x0, y0), (x1, y1), w, L, net in self.segments:
                stamp_segment(self.occ[L], to_cell(x0, y0), to_cell(x1, y1),
                              w / 2, self.nid(net))
            for x, y, net, size, drill in self.vias:
                cx, cy = to_cell(x, y)
                for LL in self.layers:
                    stamp(self.occ[LL], cx, cy, size / 2, self.nid(net))
                stamp(self.holes, cx, cy, drill / 2, True)
        return dropped

    def via_fits(self, x, y, net, size=VIA_SIZE, drill=VIA_DRILL):
        """Can a via sit here without violating copper OR hole clearance?

        Two separate rules apply and the hole rule is the stricter one on this
        board: the annular ring must keep CLEARANCE from foreign copper, and
        the drilled hole must keep HOLE_CLEARANCE from it. Ignoring the hole
        rule produced 58 hole_clearance defects on the previous pass."""
        me = self.nid(net)
        r = max(size / 2 + CLEARANCE,
                drill / 2 + HOLE_CLEARANCE) + RASTER_MARGIN
        rc = max(1, int(math.ceil(r / GRID)))
        cx, cy = to_cell(x, y)
        if not (rc <= cx < NX - rc and rc <= cy < NY - rc):
            return False
        # hole-to-hole, against every existing hole regardless of net
        rh = drill / 2 + HOLE_TO_HOLE + RASTER_MARGIN
        rhc = max(1, int(math.ceil(rh / GRID)))
        if rhc <= cx < NX - rhc and rhc <= cy < NY - rhc:
            dh = disk(rh)
            hs = self.holes[cy - rhc:cy + rhc + 1, cx - rhc:cx + rhc + 1]
            if hs.shape == dh.shape and np.any(hs[dh]):
                return False
        d = disk(r)
        for L in self.layers:
            sub = self.occ[L][cy - rc:cy + rc + 1, cx - rc:cx + rc + 1]
            if sub.shape != d.shape:
                return False
            hit = sub[d]
            if np.any((hit != 0) & (hit != me)):
                return False
        return True

    def place_via(self, x, y, net, size=VIA_SIZE, drill=VIA_DRILL):
        me = self.nid(net)
        cx, cy = to_cell(x, y)
        for L in self.layers:
            stamp(self.occ[L], cx, cy, size / 2, me)
        stamp(self.holes, cx, cy, drill / 2, True)
        self.vias.append((x, y, net, size, drill))

    def thermal_pad_vias(self):
        """Put ground vias inside every exposed thermal pad that lacks them.

        On a QFN/WSON/HTSSOP the exposed pad is both the part's ground
        connection and its only real heat path, and a pad that only touches
        copper on its own layer has neither. U2's LM5165 pad had zero vias:
        electrically it was relying on the F.Cu pour alone, and thermally on
        nothing at all.

        These are vias in pad, which the assembler needs to know about -- they
        should be tented or plugged so the pad does not wick solder down the
        barrel. That is a normal note on a board with exposed-pad parts, and it
        is in the fabrication notes.
        """
        big = {}
        for p in self.pads:
            if p["net"] != "GND" or p["any_layer"]:
                continue
            if min(p["w"], p["h"]) < 1.6:      # not an exposed pad
                continue
            big.setdefault(p["ref"], []).append(p)

        added = 0
        for ref, pads in sorted(big.items()):
            for p in pads:
                # a grid inside the pad, inset so the ring stays on the pad
                inset = SIG_VIA_SIZE / 2 + 0.15
                hw, hh = p["w"] / 2 - inset, p["h"] / 2 - inset
                if hw <= 0 or hh <= 0:
                    continue
                nx = max(1, int(2 * hw / 1.0) + 1)
                ny = max(1, int(2 * hh / 1.0) + 1)
                for i in range(nx):
                    for j in range(ny):
                        x = p["x"] - hw + (2 * hw * i / (nx - 1) if nx > 1 else hw)
                        y = p["y"] - hh + (2 * hh * j / (ny - 1) if ny > 1 else hh)
                        # A 0.6mm via does not fit a 1.65mm exposed pad ringed
                        # by 0.5mm-pitch pins, which is exactly the LM5165's
                        # pad -- it got no vias at all. Fall back to the small
                        # via rather than leave the pad unstitched.
                        for vs, vd in ((GND_VIA_SIZE, GND_VIA_DRILL),
                                       (SIG_VIA_SIZE, SIG_VIA_DRILL)):
                            if self.via_fits(x, y, "GND", vs, vd):
                                self.place_via(x, y, "GND", vs, vd)
                                added += 1
                                break
        return added

    def gnd_stitch_vias(self):
        """Stitch the F.Cu GND pour to the inner planes.

        These are NOT punched through GND pads any more. With a GND pour on
        F.Cu, an SMD GND pad already connects to copper on its own layer, so
        the stitching vias only have to tie the pours together -- which means
        they can be placed wherever they legally fit.

        Placed AFTER signal routing, deliberately. The first router had to
        place them first, because nothing stopped a via from later landing on
        top of an already-routed track. via_fits() now checks against every
        piece of copper on the board, so the ordering can be inverted -- and
        it must be: seeding 561 stitching vias first walled off the escape
        corridors out of the fine-pitch ICs and stalled 19 of 27 nets."""
        placed = []
        # a coarse lattice over the board, keeping clear of the edge
        step = 10.0
        y = 8.0
        while y < BOARD_H - 8.0:
            x = 8.0
            while x < BOARD_W - 8.0:
                if self.via_fits(x, y, "GND", GND_VIA_SIZE, GND_VIA_DRILL):
                    self.place_via(x, y, "GND", GND_VIA_SIZE, GND_VIA_DRILL)
                    placed.append((x, y))
                x += step
            y += step
        # plus a via close to each SMD GND pad where one will fit, so the pad
        # has a short path to the inner planes rather than relying on the pour
        for p in self.pads:
            if p["net"] != "GND" or p["any_layer"]:
                continue
            if any(math.hypot(p["x"] - ox, p["y"] - oy) < 4.0 for ox, oy in placed):
                continue
            for rad in (0.9, 1.3, 1.8, 2.5):
                found = False
                for k in range(12):
                    ang = 2 * math.pi * k / 12
                    vx, vy = p["x"] + rad * math.cos(ang), p["y"] + rad * math.sin(ang)
                    if self.via_fits(vx, vy, "GND", GND_VIA_SIZE, GND_VIA_DRILL):
                        self.place_via(vx, vy, "GND", GND_VIA_SIZE, GND_VIA_DRILL)
                        placed.append((vx, vy)); found = True; break
                if found: break
        return len(placed)

    def stitch_pour_regions(self):
        """Make sure every region of the F.Cu ground pour reaches the planes.

        The lattice + per-pad stitching above covers the board evenly, but
        routing chops the top pour into regions, and a region that happens to
        contain GND pads is not an island -- KiCad keeps it, and it is then
        real ground copper with no path to In1/B.Cu. DRC reported six such
        regions. On a switching board that is not a cosmetic complaint: it is
        return current with nowhere to go.

        The pour is not in the raster (KiCad fills it), so it is reconstructed
        the way KiCad would fill it -- everything the pour can reach is
        everything that is not within clearance of another net's copper --
        then labelled, and any region holding a GND pad but no via gets one.
        """
        gid = self.nid("GND")
        o = self.occ["F.Cu"]
        foreign = (o != 0) & (o != gid)
        se = disk(CLEARANCE + 0.125 + RASTER_MARGIN)
        pour = ~ndimage.binary_dilation(foreign, structure=se)
        pour[:int(3 / GRID), :] = False; pour[-int(3 / GRID):, :] = False
        pour[:, :int(3 / GRID)] = False; pour[:, -int(3 / GRID):] = False
        lbl, n = ndimage.label(pour)
        if not n:
            return 0
        has_via = set()
        for x, y, net, size, drill in self.vias:
            if net != "GND": continue
            cx, cy = to_cell(x, y)
            if 0 <= cy < NY and 0 <= cx < NX and lbl[cy, cx]:
                has_via.add(int(lbl[cy, cx]))
        # which regions actually carry GND pads?
        want = {}
        for p in self.pads:
            if p["net"] != "GND": continue
            cx, cy = to_cell(p["x"], p["y"])
            for dx in range(-12, 13, 2):
                for dy in range(-12, 13, 2):
                    yy, xx = cy + dy, cx + dx
                    if 0 <= yy < NY and 0 <= xx < NX and lbl[yy, xx]:
                        want.setdefault(int(lbl[yy, xx]), (p["x"], p["y"]))
        added = 0
        for region, (px, py) in want.items():
            if region in has_via:
                continue
            ys, xs = np.where(lbl == region)
            if len(ys) < 60:            # too small to be worth a via
                continue
            d = (xs - to_cell(px, py)[0]) ** 2 + (ys - to_cell(px, py)[1]) ** 2
            for i in np.argsort(d)[:600]:
                x, y = xs[i] * GRID, ys[i] * GRID
                if self.via_fits(x, y, "GND", GND_VIA_SIZE, GND_VIA_DRILL):
                    self.place_via(x, y, "GND", GND_VIA_SIZE, GND_VIA_DRILL)
                    added += 1
                    break
        return added

    def blocked_masks(self, net, half_w, layers, win=None, via=None, extra=0.0):
        """Cells where MY centreline may not go: other-net copper dilated by
        (my half width + clearance).

        Only the search window is dilated. Dilating all 1952x1352 cells of
        every layer, twice, on every one of the ~3000 routing attempts was
        where essentially all of the run time went; a typical link needs a
        window a few percent of that area. The slice is grown by the
        structuring element's radius before dilating and cropped back after,
        so copper just outside the window still casts its shadow inside it and
        the result is identical to dilating the whole board.
        """
        me = self.nid(net)
        se = disk(half_w + CLEARANCE + RASTER_MARGIN + extra)
        # a via is wider than the trace and drills a hole, so layer changes get
        # their own, stricter mask -- checking a via against the trace's mask
        # is what let vias land too close to foreign copper previously
        vs, vd = via if via else via_for(net)
        se_via = disk(max(vs / 2 + CLEARANCE,
                          vd / 2 + HOLE_CLEARANCE) + RASTER_MARGIN + extra)
        pad = max(se.shape[0], se_via.shape[0]) // 2 + 1
        if win is None:
            wx0, wx1, wy0, wy1 = 0, NX - 1, 0, NY - 1
        else:
            wx0, wx1, wy0, wy1 = win
        sx0 = max(0, wx0 - pad); sx1 = min(NX, wx1 + pad + 1)
        sy0 = max(0, wy0 - pad); sy1 = min(NY, wy1 + pad + 1)
        cx0, cx1 = wx0 - sx0, wx1 - sx0 + 1
        cy0, cy1 = wy0 - sy0, wy1 - sy0 + 1
        # hole-to-hole applies between ALL holes, same net or not
        se_h = disk(vs_drill_half(via, net) + HOLE_TO_HOLE + RASTER_MARGIN)
        hsub = self.holes[sy0:sy1, sx0:sx1]
        hole_bad = np.zeros((NY, NX), dtype=bool)
        if hsub.any():
            hole_bad[wy0:wy1 + 1, wx0:wx1 + 1] = ndimage.binary_dilation(
                hsub, structure=se_h)[cy0:cy1, cx0:cx1]

        out, out_via = {}, {}
        # A via is a plated barrel through the whole stack: it must clear
        # copper on every layer, not just the two the route happens to be
        # using. Routing on F.Cu+In2.Cu alone left B.Cu unchecked, and the
        # via audit caught 12 vias placed straight through B.Cu copper.
        for L in set(layers) | set(self.layers):
            sub = self.occ[L][sy0:sy1, sx0:sx1]
            others = (sub != 0) & (sub != me)
            full = np.zeros((NY, NX), dtype=bool)
            full_v = np.zeros((NY, NX), dtype=bool)
            if others.any():
                full[wy0:wy1 + 1, wx0:wx1 + 1] = ndimage.binary_dilation(
                    others, structure=se)[cy0:cy1, cx0:cx1]
                full_v[wy0:wy1 + 1, wx0:wx1 + 1] = ndimage.binary_dilation(
                    others, structure=se_via)[cy0:cy1, cx0:cx1]
            out[L] = full; out_via[L] = full_v | hole_bad
        return out, out_via

    def route(self, net, a, b, width, layers, margin_mm=18.0, extra_clear=0.0):
        """A* between two pads. Endpoints are pinned to the layers the pad
        actually exists on (an SMD pad is F.Cu only -- letting a route 'start'
        on another layer with no via was the single biggest connectivity bug
        in the first router)."""
        half = width / 2
        sx, sy = to_cell(a["ax"], a["ay"]); tx, ty = to_cell(b["ax"], b["ay"])
        # Detour allowance around the straight line. 60mm was generous enough
        # to make every window most of the board; 18mm still permits a healthy
        # dogleg and the caller widens it if the link genuinely fails.
        margin = int(margin_mm / GRID)
        x0 = max(0, min(sx, tx) - margin); x1 = min(NX - 1, max(sx, tx) + margin)
        y0 = max(0, min(sy, ty) - margin); y1 = min(NY - 1, max(sy, ty) + margin)
        blocked, via_blocked = self.blocked_masks(net, half, layers,
                                                  (x0, x1, y0, y1),
                                                  via=via_for(net),
                                                  extra=extra_clear)
        src_layers = [i for i, L in enumerate(layers) if L in a["layers"]]
        dst_layers = [i for i, L in enumerate(layers) if L in b["layers"]]
        if not src_layers or not dst_layers: return None

        # the endpoints themselves sit on their own pads, which are blocked
        # for every other net but must be enterable for this one; the mask is
        # already net-aware, so nothing special is needed unless the pad is
        # boxed in -- in which case we punch a small hole at the exact pad.
        for li in src_layers: blocked[layers[li]][sy, sx] = False
        for li in dst_layers: blocked[layers[li]][ty, tx] = False

        goals = {(li, tx, ty) for li in dst_layers}
        step_diag = math.sqrt(2) * GRID
        openq = []
        for li in src_layers:
            openq.append((0.0, 0.0, (li, sx, sy), None))
        heapq.heapify(openq)
        came, best, seen = {}, {}, set()
        for li in src_layers: best[(li, sx, sy)] = 0.0
        MOVES = [(1,0,GRID),(-1,0,GRID),(0,1,GRID),(0,-1,GRID),
                 (1,1,step_diag),(1,-1,step_diag),(-1,1,step_diag),(-1,-1,step_diag)]
        while openq:
            f, g, node, parent = heapq.heappop(openq)
            if node in seen: continue
            seen.add(node); came[node] = parent
            if node in goals:
                path = []; cur = node
                while cur is not None:
                    path.append(cur); cur = came[cur]
                return path[::-1]
            li, cx, cy = node
            L = layers[li]
            for dx, dy, cost in MOVES:
                nx_, ny_ = cx + dx, cy + dy
                if not (x0 <= nx_ <= x1 and y0 <= ny_ <= y1): continue
                if blocked[L][ny_, nx_]: continue
                nn = (li, nx_, ny_); ng = g + cost
                if ng < best.get(nn, 1e18):
                    best[nn] = ng
                    h = math.hypot(nx_ - tx, ny_ - ty) * GRID
                    heapq.heappush(openq, (ng + h, ng, nn, node))
            for oli in range(len(layers)):
                if oli == li: continue
                # a layer change places a real via: it must clear foreign
                # copper on BOTH layers by the via rule, not the trace rule
                if any(via_blocked[LL][cy, cx] for LL in via_blocked):
                    continue
                nn = (oli, cx, cy); ng = g + VIA_PENALTY_MM
                if ng < best.get(nn, 1e18):
                    best[nn] = ng
                    h = math.hypot(cx - tx, cy - ty) * GRID
                    heapq.heappush(openq, (ng + h, ng, nn, node))
        return None

    def commit(self, net, path, width, layers, a, b):
        """Turn a cell path into real segments/vias and stamp it as copper."""
        me = self.nid(net); half = width / 2
        for e in (a, b):
            self.used_anchors.add((round(e["ax"], 3), round(e["ay"], 3)))
        # collapse collinear runs per layer
        i = 0; n = len(path)
        pts_out = []
        while i < n - 1:
            li = path[i][0]; j = i
            while j + 1 < n and path[j + 1][0] == li: j += 1
            run = path[i:j + 1]
            # simplify: keep direction changes only
            keep = [run[0]]
            for k in range(1, len(run) - 1):
                px, py = run[k - 1][1] - run[k][1], run[k - 1][2] - run[k][2]
                nx_, ny_ = run[k][1] - run[k + 1][1], run[k][2] - run[k + 1][2]
                if (px, py) != (nx_, ny_): keep.append(run[k])
            if len(run) > 1: keep.append(run[-1])
            xy = [(c[1] * GRID, c[2] * GRID) for c in keep]
            # Reach the exact pad centre with a short LANDING STUB rather than
            # by moving the first/last vertex onto it. Moving it (the previous
            # behaviour) shifted a whole simplified segment by up to half a
            # cell away from the centreline A* actually validated -- which is
            # why 71 pad-to-track clearances landed at a median of 0.163mm
            # against the 0.2mm rule, worst right at the pads. The stub is at
            # most one cell diagonal (~0.07mm) and stays inside the pad.
            if i == 0:
                xy.insert(0, (a["ax"], a["ay"]))
            if j == n - 1:
                xy.append((b["ax"], b["ay"]))
            for k in range(len(xy) - 1):
                # compare by distance, not tuple equality: the landing stub
                # reaches an anchor stored as a float, so 15.9 and
                # 15.900000000000002 are "different" points 0.000000000002mm
                # apart and KiCad reports the result as a dangling track
                if math.hypot(xy[k + 1][0] - xy[k][0],
                              xy[k + 1][1] - xy[k][1]) < 1e-6:
                    continue
                self.segments.append((xy[k], xy[k + 1], width, layers[li], net))
                stamp_segment(self.occ[layers[li]],
                              to_cell(*xy[k]), to_cell(*xy[k + 1]), half, me)
            if j + 1 < n:  # layer change -> via
                vx, vy = path[j][1] * GRID, path[j][2] * GRID
                vs, vd = via_for(net)
                # AUDIT: A* said this layer change was legal. Ask the
                # independent point-check whether it agrees. Any disagreement
                # is a bug in the mask, not a close call.
                if not self.via_fits(vx, vy, net, vs, vd):
                    self.via_audit.append((net, vx, vy, layers[li],
                                           layers[path[j + 1][0]]))
                self.vias.append((vx, vy, net, vs, vd))
                cx, cy = to_cell(vx, vy)
                for L in self.layers:
                    stamp(self.occ[L], cx, cy, vs / 2, me)
                stamp(self.holes, cx, cy, vd / 2, True)
            i = j + 1

    def emit(self):
        out = []
        for (x0, y0), (x1, y1), w, layer, net in self.segments:
            out.append(
                f'\t(segment\n\t\t(start {x0:.4f} {y0:.4f})\n\t\t(end {x1:.4f} {y1:.4f})\n'
                f'\t\t(width {w})\n\t\t(layer "{layer}")\n\t\t(net "{net}")\n'
                f'\t\t(uuid "{uid()}")\n\t)')
        for x, y, net, size, drill in self.vias:
            out.append(
                f'\t(via\n\t\t(at {x:.4f} {y:.4f})\n\t\t(size {size})\n\t\t(drill {drill})\n'
                f'\t\t(layers "F.Cu" "B.Cu")\n\t\t(net "{net}")\n\t\t(uuid "{uid()}")\n\t)')
        return out


def chain_order(pads):
    """Nearest-neighbour ordering over a net's pads."""
    rem = pads[1:]; order = [pads[0]]
    while rem:
        last = order[-1]
        rem.sort(key=lambda p: (p["x"] - last["x"])**2 + (p["y"] - last["y"])**2)
        order.append(rem.pop(0))
    return order


def main():
    text = open(PCB_PATH, errors='replace').read()
    pads = parse_pads(text)
    # de-dupe physically identical pad entries
    seen = set(); uniq = []
    for p in pads:
        # The layer set is part of the identity. U4's exposed pad has an
        # F.Cu pad and a B.Cu pad at the SAME coordinates with the same pad
        # number, and keying on (ref, num, x, y) alone threw the bottom one
        # away as a duplicate -- so the router never saw it, and the B.Cu
        # escape routes ran straight through it.
        k = (p["ref"], p["num"], round(p["x"], 3), round(p["y"], 3),
             p["layers"])
        if k in seen: continue
        seen.add(k); uniq.append(p)
    pads = uniq
    print(f"parsed {len(pads)} pads")

    R = Router(pads)
    R.stamp_pads()
    nstub = R.escape_stubs()
    R.stub_segs = len(R.segments)
    print(f"escape stubs (exact geometry, fine-pitch pads): {nstub}")

    by_net = defaultdict(list)
    for p in pads:
        if p["net"]: by_net[p["net"]].append(p)

    pins_per_ref = defaultdict(int)
    for p in pads: pins_per_ref[p["ref"]] += 1

    nets = [(n, ps) for n, ps in by_net.items() if n != "GND" and len(ps) >= 2]
    # route the highest-current and the densest-IC nets first: they have the
    # least freedom, and whatever routes first gets the good corridors
    def prio(item):
        """Widest current first, then densest package, then fewest pads.

        Putting crowded-pin-field nets first instead was tried, on the theory
        that fanout room is scarcer than trace width: it measured worse (44 of
        60 nets against 48, 19 failed links against 14). Routing signals early
        lets them take direct paths straight across the board, and the
        passive-to-passive links that come later -- the U4 feedback divider
        chain, the MCU's port resistors -- are then the ones with nowhere to
        go. Congestion is not concentrated at the pin fields; it is global.
        """
        n, ps = item
        cur = NET_CURRENT_A.get(n, 0.0)
        dense = max(pins_per_ref[p["ref"]] for p in ps)
        return (-cur, -dense, len(ps))
    nets.sort(key=prio)

    done = 0; failed = []; hard = []
    for name, ps in nets:
        target_w = width_for(name)
        order = chain_order(ps)
        ok_links = 0; need = 0
        for a, b in zip(order, order[1:]):
            # Two pads that physically touch need no trace between them.
            # Tested on real rectangles: the circular-envelope version of this
            # test called U8's three separate /VBATT pads "touching" because a
            # 0.25x0.88mm pad's envelope radius is 0.457mm and they sit 0.5mm
            # apart, so the links were skipped and the pads left unconnected.
            if pads_abut(a, b):
                continue
            need += 1
            # a trace cannot be wider than the narrower of the pads it lands
            # on, or it spills onto their neighbours on a fine-pitch part
            cap = max(MIN_W, min(a["narrow"], b["narrow"]))
            eff_w = min(target_w, cap)
            path = None; used_w = None
            # Escalation ladder, cheapest first. Each rung relaxes exactly one
            # thing -- width, then the B.Cu escape layer, then how far the
            # route may wander -- so a link that succeeds early never pays for
            # the expensive rungs, and a link that fails has demonstrably been
            # tried with everything.
            LADDER = [
                (eff_w,    ROUTE_LAYERS,                    18.0),
                (eff_w,    ROUTE_LAYERS + [ESCAPE_LAYER],   18.0),
                (SIGNAL_W, ROUTE_LAYERS + [ESCAPE_LAYER],   18.0),
                (MIN_W,    ROUTE_LAYERS + [ESCAPE_LAYER],   35.0),
                (MIN_W,    ROUTE_LAYERS + [ESCAPE_LAYER],   80.0),
            ]
            for w, lys, mm in LADDER:
                if w > eff_w: continue
                path = R.route(name, a, b, w, lys, margin_mm=mm)
                if path:
                    used_w = w; break
            if path:
                # commit needs the same layer list the path was found on
                layers = (ROUTE_LAYERS + [ESCAPE_LAYER]
                          if max(p[0] for p in path) >= len(ROUTE_LAYERS)
                          else ROUTE_LAYERS)
                s0, v0 = len(R.segments), len(R.vias)
                R.commit(name, path, used_w, layers, a, b)
                R.link_span.append(((name, a, b, eff_w), s0, len(R.segments),
                                    v0, len(R.vias)))
                ok_links += 1
            else:
                failed.append((name, a["ref"] + "." + a["num"], b["ref"] + "." + b["num"]))
                hard.append((name, a, b, eff_w))
        if ok_links == need: done += 1
        print(f"  {name:24s} w={target_w:.2f} links {ok_links}/{need}", flush=True)

    if hard:
        nfix, hard = hard_link_pass(R, hard)
        # `hard` is now the authoritative list of what is still unrouted
        failed = [(n, a2["ref"] + "." + a2["num"], b2["ref"] + "." + b2["num"])
                  for (n, a2, b2, _w) in hard]
        print(f"\nhard-link pass closed {nfix} link(s); "
              f"{len(failed)} still unrouted", flush=True)

    # ----------------------------------------------------------------------
    # Exact-geometry rip-up. Up to here the board is only as good as a 0.1mm
    # raster; from here it is checked against real geometry and made to comply.
    # A link that violates a clearance is removed, the raster is rebuilt
    # without it, and it is routed again into the space that is now free --
    # with a wider safety margin, since the first attempt demonstrably sat too
    # close to something. A link that still cannot be placed legally is left
    # out and reported, because an unrouted net is a visible, countable gap
    # whereas a short is a dead board.
    # ----------------------------------------------------------------------
    for round_no in range(1, 5):
        bad = R.offending_links()
        if not bad:
            print(f"\nexact-geometry check: clean after {round_no - 1} rip-up round(s)")
            break
        print(f"\nrip-up round {round_no}: {len(bad)} links violate exact clearance "
              f"(worst {bad[0][1]:+.3f}mm)")
        links = [R.link_span[li][0] for li, _ in bad]
        R.rip([li for li, _ in bad])
        redone = 0
        for (name, a, b, eff_w) in links:
            got = None
            for w, lys, mm in [(eff_w, ROUTE_LAYERS + [ESCAPE_LAYER], 25.0),
                               (SIGNAL_W, ROUTE_LAYERS + [ESCAPE_LAYER], 25.0),
                               (MIN_W, ROUTE_LAYERS + [ESCAPE_LAYER], 45.0)]:
                if w > eff_w: continue
                got = R.route(name, a, b, w, lys, margin_mm=mm,
                              extra_clear=RASTER_MARGIN)
                if got: break
            if got:
                lys2 = (ROUTE_LAYERS + [ESCAPE_LAYER]
                        if max(p[0] for p in got) >= len(ROUTE_LAYERS)
                        else ROUTE_LAYERS)
                s0, v0 = len(R.segments), len(R.vias)
                R.commit(name, got, w, lys2, a, b)
                R.link_span.append(((name, a, b, eff_w), s0, len(R.segments),
                                    v0, len(R.vias)))
                redone += 1
            else:
                failed.append((name, a["ref"] + "." + a["num"],
                               b["ref"] + "." + b["num"]))
        print(f"   re-routed {redone}/{len(links)} clear of the obstruction")
    else:
        left = R.offending_links()
        if left:
            print(f"\nrip-up gave up with {len(left)} violating links -- removing them")
            for (name, a, b, _w) in [R.link_span[li][0] for li, _ in left]:
                failed.append((name, a["ref"] + "." + a["num"],
                               b["ref"] + "." + b["num"]))
            R.rip([li for li, _ in left])

    if R.via_audit:
        print(f"\n!! VIA AUDIT: {len(R.via_audit)} vias A* accepted but the "
              f"point-check rejects")
        for a in R.via_audit[:12]: print("   ", a)
    else:
        print("\nvia audit: A* and the point-check agree on every via")

    # ----------------------------------------------------------------------
    # Connectivity repair. Routing works link by link along a chain, and the
    # rip-up pass removes links -- so a net can finish in two or more pieces
    # even though every individual link "succeeded". That is how Net-(Q4-D)
    # ended up split: its chain routed 8 of 8, then rip-up removed
    # C19.1->R13.1 to clear a 0.097mm violation and took the rest of the net
    # with it, leaving the charge shunt RS2 stranded. The electrothermal model
    # found it before DRC did, by refusing to solve a case whose terminal had
    # no copper.
    #
    # So connectivity is now checked as connectivity: group each net's pads
    # into the pieces that are actually joined, then route the shortest link
    # between two pieces, and repeat while that keeps working.
    # ----------------------------------------------------------------------
    for attempt in range(4):
        joined = union_find_pads(R)
        holes_left = []
        by_net = defaultdict(list)
        for pd in R.pads:
            if pd["net"] and pd["net"] != "GND":
                by_net[pd["net"]].append(pd)
        for name, ps in by_net.items():
            groups = defaultdict(list)
            for pd in ps:
                groups[joined[pad_key(pd)]].append(pd)
            if len(groups) < 2:
                continue
            gs = sorted(groups.values(), key=len, reverse=True)
            best = None
            for other in gs[1:]:
                for a in gs[0]:
                    for b in other:
                        d = math.hypot(a["x"] - b["x"], a["y"] - b["y"])
                        if best is None or d < best[0]:
                            best = (d, a, b)
            if best:
                holes_left.append((name, best[1], best[2]))
        if not holes_left:
            print(f"\nconnectivity repair: every net in one piece "
                  f"after {attempt} pass(es)")
            break
        print(f"\nconnectivity repair pass {attempt + 1}: "
              f"{len(holes_left)} nets split into pieces")
        fixed = 0
        for name, a, b in holes_left:
            got = w = None
            for w, lys, mm in [(width_for(name), ROUTE_LAYERS + [ESCAPE_LAYER], 30.0),
                               (SIGNAL_W, ROUTE_LAYERS + [ESCAPE_LAYER], 45.0),
                               (MIN_W, ROUTE_LAYERS + [ESCAPE_LAYER], 90.0)]:
                got = R.route(name, a, b, w, lys, margin_mm=mm)
                if got: break
            if got:
                lys2 = (ROUTE_LAYERS + [ESCAPE_LAYER]
                        if max(q[0] for q in got) >= len(ROUTE_LAYERS)
                        else ROUTE_LAYERS)
                s0, v0 = len(R.segments), len(R.vias)
                R.commit(name, got, w, lys2, a, b)
                R.link_span.append(((name, a, b, w), s0, len(R.segments),
                                    v0, len(R.vias)))
                fixed += 1
        print(f"   reconnected {fixed}/{len(holes_left)}")
        if not fixed:
            break

    dropped = R.prune_unused_stubs()
    print(f"unused escape stubs removed: {dropped}")

    nth = R.thermal_pad_vias()
    print(f"thermal vias placed in exposed ground pads: {nth}")

    ngnd = R.gnd_stitch_vias()
    npour = R.stitch_pour_regions()
    print(f"\nGND stitching vias (placed last, into leftover space): {ngnd}"
          f"  (+{npour} tying isolated pour regions to the planes)")
    print(f"nets fully routed: {done}/{len(nets)}   failed links: {len(failed)}")
    for f in failed: print("   FAIL", f)

    json.dump([[a, b, c] for a, b, c in failed],
              open(os.path.join(HERE, "route2_failed.json"), "w"))

    body = R.emit()
    stripped = text.rstrip()
    idx = stripped.rfind('\n')
    new = stripped[:idx + 1] + '\n'.join(body) + '\n' + stripped[idx + 1:] + '\n'
    open(PCB_PATH, "w").write(new)
    print(f"wrote {len(R.segments)} segments + {len(R.vias)} vias")


if __name__ == "__main__":
    main()