Add balanced quadtree grid engine (adaptive cells, phase 1)

fill_resistance/quadtree.py decomposes a layer's fine copper mask into
2:1-balanced power-of-two leaves: boundary and keep-fine cells stay at
the fine size, interiors coarsen with their Chebyshev distance to the
nearest feature (guard factor, default 4), and an explicit enforcement
pass splits any leaf more than twice an edge-adjacent neighbor. Face
conductances use the series-half-cell rule, which reduces to the
production harmonic mean for equal sizes and EXACTLY to sigma in the
uniform limit - verified edge-for-edge against solver.build_edges and
to rel 1e-12 in R against run_solve, so the exact-value test suite
stays authoritative for this engine.

Measured (feature-dense 120x120 plate, h=50um, production AMG solver):
uniform 5.58M unknowns ~35s; adaptive guard=4 823k / ~8s at -1.1%;
guard=8 1.74M / ~12s at -0.47%. tools/adaptive_proto.py now benchmarks
the engine itself.

Not yet wired into the pipeline: phase 2 ports electrodes, barrels,
1D chains, buildup and field output onto leaves.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

fill_resistance/quadtree.py

@@ -0,0 +1,171 @@
+"""Adaptive quadtree grid engine (phase 1 of the adaptive-cell work).
+
+Decomposes a layer's fine copper mask into a 2:1-BALANCED set of square
+leaves (power-of-two sizes in fine-cell units): cells at copper
+boundaries and in keep_fine regions stay at the fine size, interiors
+coarsen with their Chebyshev distance to the nearest feature (guard
+factor), and an explicit enforcement pass splits any leaf more than
+twice the size of an edge-adjacent neighbor.
+
+The uniform limit (max_block=1) reproduces the fine grid EXACTLY: one
+leaf per copper cell and face conductance sigma - the production
+solver's grid is the special case, which keeps the exact-value test
+suite authoritative for this engine.
+
+Face conductance between edge-adjacent leaves a, b sharing w fine-cell
+faces is the series-half-cell expression
+
+    g = w / (size_a / (2 sigma_a) + size_b / (2 sigma_b))
+
+which reduces to the harmonic mean 2 sigma_a sigma_b / (sigma_a +
+sigma_b) for equal sizes (the production buildup/mouth face rule) and
+to sigma in the uniform limit.
+
+Phase 2 (not here) wires this into the pipeline: electrodes, barrels,
+1D chains, buildup and field output still run on the uniform grid.
+"""
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+from scipy import ndimage
+
+
+@dataclass
+class LeafGrid:
+    """Square leaves over one layer, in fine-cell units."""
+    y0: np.ndarray                # int32, aligned: y0 % size == 0
+    x0: np.ndarray
+    size: np.ndarray              # int32, power of two
+    id_grid: np.ndarray           # (ny, nx) int32; -1 = not copper
+
+    @property
+    def n(self) -> int:
+        return len(self.size)
+
+
+def build_leaves(mask: np.ndarray, keep_fine: np.ndarray | None = None,
+                 max_block: int = 32, guard: int = 4) -> LeafGrid:
+    """Balanced leaf decomposition of a boolean copper mask. keep_fine
+    marks cells that must stay at the fine size (electrodes, and later
+    via mouths / buildup edges). guard scales how much clearance a block
+    of size s needs (Chebyshev distance >= guard * s)."""
+    ny, nx = mask.shape
+    if keep_fine is None:
+        keep_fine = np.zeros_like(mask)
+    coarsenable = mask & ~keep_fine
+
+    # Chebyshev distance to the nearest non-coarsenable cell (array
+    # border padded as background so edges never look like interior)
+    pad = np.pad(coarsenable, 1)
+    D = ndimage.distance_transform_cdt(pad, metric="chessboard")[1:-1, 1:-1]
+
+    S = np.zeros((ny, nx), dtype=np.int32)
+    S[mask] = 1
+    s = 2
+    while s <= max_block:
+        S[D >= guard * s] = s
+        s *= 2
+
+    grid = _emit(S, mask, max_block)
+    for _ in range(32):
+        if _split_unbalanced(grid, S):
+            grid = _emit(S, mask, max_block)
+        else:
+            return grid
+    raise RuntimeError("quadtree balance did not converge")
+
+
+def _emit(S: np.ndarray, mask: np.ndarray, max_block: int) -> LeafGrid:
+    """Greedy top-down emission: an aligned s-block becomes a leaf where
+    every cell allows size s and nothing larger claimed it."""
+    ny, nx = mask.shape
+    py, px = (-ny) % max_block, (-nx) % max_block
+    Sp = np.pad(S, ((0, py), (0, px)))
+    NY, NX = Sp.shape
+    id_grid = np.full((NY, NX), -1, dtype=np.int32)
+    covered = np.zeros((NY, NX), dtype=bool)
+    y0s, x0s, sizes = [], [], []
+    nid = 0
+    s = max_block
+    while s >= 2:
+        min_S = Sp.reshape(NY // s, s, NX // s, s).min(axis=(1, 3))
+        free = ~covered.reshape(NY // s, s, NX // s, s).any(axis=(1, 3))
+        cand = (min_S >= s) & free
+        k = int(cand.sum())
+        if k:
+            lvl = np.full(cand.shape, -1, dtype=np.int32)
+            lvl[cand] = nid + np.arange(k, dtype=np.int32)
+            up = np.repeat(np.repeat(lvl, s, axis=0), s, axis=1)
+            sel = up >= 0
+            id_grid[sel] = up[sel]
+            covered |= sel
+            ii, jj = np.nonzero(cand)
+            y0s.append(ii.astype(np.int32) * s)
+            x0s.append(jj.astype(np.int32) * s)
+            sizes.append(np.full(k, s, dtype=np.int32))
+            nid += k
+        s //= 2
+    fi, fj = np.nonzero(np.pad(mask, ((0, py), (0, px))) & ~covered)
+    id_grid[fi, fj] = nid + np.arange(len(fi), dtype=np.int32)
+    y0s.append(fi.astype(np.int32))
+    x0s.append(fj.astype(np.int32))
+    sizes.append(np.ones(len(fi), dtype=np.int32))
+    return LeafGrid(
+        y0=np.concatenate(y0s) if y0s else np.zeros(0, np.int32),
+        x0=np.concatenate(x0s) if x0s else np.zeros(0, np.int32),
+        size=np.concatenate(sizes) if sizes else np.zeros(0, np.int32),
+        id_grid=id_grid[:ny, :nx],
+    )
+
+
+def _split_unbalanced(grid: LeafGrid, S: np.ndarray) -> bool:
+    """Cap S over any leaf more than 2x an edge-adjacent neighbor, so the
+    next emission splits it. Returns True if anything was capped."""
+    ia, ib, _ = leaf_faces(grid)
+    sa, sb = grid.size[ia], grid.size[ib]
+    big = np.unique(np.concatenate([ia[sa > 2 * sb], ib[sb > 2 * sa]]))
+    for lid in big:
+        y, x, s = int(grid.y0[lid]), int(grid.x0[lid]), int(grid.size[lid])
+        np.minimum(S[y:y + s, x:x + s], s // 2, out=S[y:y + s, x:x + s])
+    return len(big) > 0
+
+
+def leaf_faces(grid: LeafGrid) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """(ia, ib, w) for every edge-adjacent leaf pair, each pair once;
+    w = shared face length in fine-cell units."""
+    n = grid.n
+    if n == 0:
+        z = np.zeros(0, dtype=np.int64)
+        return z, z, z
+    out = []
+    for sl_a, sl_b in ((np.s_[:, :-1], np.s_[:, 1:]),
+                       (np.s_[:-1, :], np.s_[1:, :])):
+        a = grid.id_grid[sl_a].ravel().astype(np.int64)
+        b = grid.id_grid[sl_b].ravel().astype(np.int64)
+        ok = (a >= 0) & (b >= 0) & (a != b)
+        key, counts = np.unique(a[ok] * n + b[ok], return_counts=True)
+        out.append((key // n, key % n, counts))
+    ia = np.concatenate([o[0] for o in out])
+    ib = np.concatenate([o[1] for o in out])
+    w = np.concatenate([o[2] for o in out])
+    return ia, ib, w
+
+
+def leaf_edges(grid: LeafGrid, sigma) -> tuple[np.ndarray, np.ndarray,
+                                               np.ndarray]:
+    """Face conductances [S]: sigma is a scalar or per-leaf array of
+    sheet conductance. Uniform limit -> exactly sigma per face."""
+    ia, ib, w = leaf_faces(grid)
+    sig = np.broadcast_to(np.asarray(sigma, dtype=float), (grid.n,))
+    g = w / (grid.size[ia] / (2.0 * sig[ia])
+             + grid.size[ib] / (2.0 * sig[ib]))
+    return ia, ib, g
+
+
+def balanced(grid: LeafGrid) -> bool:
+    """2:1 balance invariant: edge-adjacent leaves differ <= 2x in size."""
+    ia, ib, _ = leaf_faces(grid)
+    sa, sb = grid.size[ia], grid.size[ib]
+    return bool((np.maximum(sa, sb) <= 2 * np.minimum(sa, sb)).all())

tests/test_quadtree.py

@@ -0,0 +1,141 @@
+"""Quadtree grid engine tests (phase 1): exact uniform limit against the
+production graph and solver, partition/alignment/balance invariants, and
+adaptive-vs-fine R agreement."""
+import numpy as np
+import pytest
+from scipy import ndimage
+
+from fill_resistance import quadtree, raster, solver
+from tests.util import NM, make_problem, strip_problem
+
+
+def _plate_with_holes(n=5, size_mm=40.0):
+    holes = []
+    pitch = size_mm / n
+    for i in range(n):
+        for j in range(n):
+            x, y = pitch * (i + 0.4), pitch * (j + 0.4)
+            holes.append([(x, y), (x + 1, y), (x + 1, y + 1), (x, y + 1)])
+    outline = [(0, 0), (size_mm, 0), (size_mm, size_mm), (0, size_mm)]
+    return make_problem([(outline, holes)],
+                        rect1_mm=(0, 15, 2, 25),
+                        rect2_mm=(size_mm - 2, 15, size_mm, 25))
+
+
+def _solve_on_leaves(problem, stack, e1, e2, grid):
+    """Equipotential mini-solve on the leaf graph, using the production
+    assembly and linear solver."""
+    ia, ib, g = quadtree.leaf_edges(grid, problem.sigma_s(0))
+    state = np.ones(grid.n, dtype=np.uint8)
+    for e, code in ((e1[0], 2), (e2[0], 3)):
+        ids = grid.id_grid[e]
+        state[ids[ids >= 0]] = code
+    edges = solver.Edges(a=ia, b=ib, w=g,
+                         via_index=np.full(len(ia), -1, dtype=np.int32))
+    A, rhs, _ = solver._assemble(state, edges, None)
+    x, _ = solver.solve_system(A, rhs)
+    V = np.zeros(grid.n)
+    V[state == 2] = 1.0
+    V[state == 1] = x
+    Ie = g * (V[ia] - V[ib])
+    sa, sb = state[ia], state[ib]
+    I1 = float(Ie[sa == 2].sum() - Ie[sb == 2].sum())
+    I2 = float(Ie[sb == 3].sum() - Ie[sa == 3].sum())
+    return 1.0 / (0.5 * (I1 + I2))
+
+
+def test_uniform_limit_graph_identical():
+    """max_block=1: one leaf per cell, and the edge list matches the
+    production in-plane graph exactly (same pairs, conductance sigma)."""
+    p = _plate_with_holes()
+    stack = raster.rasterize_stack(p, 0.5 * NM)
+    grid = quadtree.build_leaves(stack.masks[0], max_block=1)
+    assert int(stack.masks.sum()) == grid.n
+    assert (grid.size == 1).all()
+
+    ny, nx = stack.shape2d
+    flat_of_leaf = grid.y0.astype(np.int64) * nx + grid.x0
+    ia, ib, g = quadtree.leaf_edges(grid, p.sigma_s(0))
+    ours = np.sort(np.stack([
+        np.minimum(flat_of_leaf[ia], flat_of_leaf[ib]),
+        np.maximum(flat_of_leaf[ia], flat_of_leaf[ib])], axis=1), axis=0)
+
+    edges = solver.build_edges(stack, p, [p.sigma_s(0)])
+    ref = np.sort(np.stack([np.minimum(edges.a, edges.b),
+                            np.maximum(edges.a, edges.b)], axis=1), axis=0)
+    assert ours.shape == ref.shape
+    assert np.array_equal(np.sort(ours.view("i8,i8"), order=["f0", "f1"],
+                                  axis=0),
+                          np.sort(ref.view("i8,i8"), order=["f0", "f1"],
+                                  axis=0))
+    assert np.allclose(g, p.sigma_s(0), rtol=0, atol=0)
+
+
+def test_uniform_limit_R_matches_production():
+    p = strip_problem(length=50, width=10, e_len=5)
+    stack = raster.rasterize_stack(p, 0.25 * NM)
+    e1, e2 = raster.electrode_masks(stack, p)
+    ref = solver.run_solve(p, stack, e1, e2, 1.0,
+                           contact_model="equipotential")
+    stack2 = raster.rasterize_stack(p, 0.25 * NM)
+    e1b, e2b = raster.electrode_masks(stack2, p)
+    grid = quadtree.build_leaves(stack2.masks[0], max_block=1)
+    R = _solve_on_leaves(p, stack2, e1b, e2b, grid)
+    assert R == pytest.approx(ref.R_ohm, rel=1e-12)
+
+
+def test_partition_alignment_and_balance():
+    p = _plate_with_holes()
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    mask = stack.masks[0]
+    grid = quadtree.build_leaves(mask, max_block=32)
+
+    # exact partition of the copper
+    assert int((grid.size.astype(np.int64) ** 2).sum()) == int(mask.sum())
+    assert (grid.id_grid >= 0).sum() == int(mask.sum())
+    assert not (grid.id_grid[~mask] >= 0).any()
+    counts = np.bincount(grid.id_grid[grid.id_grid >= 0], minlength=grid.n)
+    assert np.array_equal(counts, grid.size.astype(np.int64) ** 2)
+
+    # power-of-two sizes, aligned to their own size
+    assert np.array_equal(grid.size & (grid.size - 1),
+                          np.zeros_like(grid.size))
+    assert (grid.y0 % grid.size == 0).all()
+    assert (grid.x0 % grid.size == 0).all()
+
+    # 2:1 balance and real coarsening (max size is geometry-limited by
+    # the guard distance, not by max_block, on this feature-dense plate)
+    assert quadtree.balanced(grid)
+    assert grid.n < 0.5 * int(mask.sum())
+    assert grid.size.max() >= 4
+
+
+def test_boundary_and_keep_fine_stay_fine():
+    p = _plate_with_holes()
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    mask = stack.masks[0]
+    keep = np.zeros_like(mask)
+    keep[50:60, 50:60] = True
+    grid = quadtree.build_leaves(mask, keep_fine=keep, max_block=32)
+
+    boundary = mask & ndimage.binary_dilation(~mask)
+    assert (grid.size[grid.id_grid[boundary]] == 1).all()
+    assert (grid.size[grid.id_grid[keep & mask]] == 1).all()
+
+
+def test_adaptive_R_close_to_fine():
+    """Adaptive leaves reproduce the fine-uniform R within 1% on the
+    holey plate (features everywhere - the adversarial case)."""
+    p = _plate_with_holes()
+    stack = raster.rasterize_stack(p, 0.1 * NM)
+    e1, e2 = raster.electrode_masks(stack, p)
+    ref = solver.run_solve(p, stack, e1, e2, 1.0,
+                           contact_model="equipotential")
+
+    stack2 = raster.rasterize_stack(p, 0.1 * NM)
+    e1b, e2b = raster.electrode_masks(stack2, p)
+    grid = quadtree.build_leaves(stack2.masks[0],
+                                 keep_fine=(e1b[0] | e2b[0]))
+    R = _solve_on_leaves(p, stack2, e1b, e2b, grid)
+    assert grid.n < 0.4 * ref.solve_info.n_unknowns
+    assert R == pytest.approx(ref.R_ohm, rel=0.01)

tools/adaptive_proto.py

@@ -1,27 +1,25 @@
-"""Prototype: quadtree-adaptive grid for the fill-resistance solver.
+"""Benchmark for the adaptive quadtree grid engine
+(fill_resistance/quadtree.py, phase 1 of the adaptive-cell work).
 
-Coarsens the existing fine raster bottom-up (power-of-two blocks that are
+    .venv/Scripts/python tools/adaptive_proto.py
-fully copper away from electrodes, erosion-graded per level), builds the
-leaf graph fully vectorized via an id-grid (face conductance
-g = sigma * overlap / mean-size, which reduces EXACTLY to the production
-sigma in the uniform limit), and reuses the production assembly + AMG
-solver. Run from the repo root: .venv/Scripts/python tools/adaptive_proto.py
 
-Measured 2026-07-15 (120x120 mm plate, 400 holes, h = 50 um; adversarial:
+Solves a feature-dense 120x120 mm plate (400 holes) at h = 50 um on the
-features everywhere, so geometric refinement has no smooth interior):
+uniform grid and on the balanced quadtree (engine defaults: guard=4,
+max_block=32), using the production assembly + solver for both.
 
-    uniform      R = 0.508504 mOhm   5.58M unknowns   27 s   (reference)
+Measured 2026-07-15 on this machine:
-    max_block=4  R = 0.506368 mOhm   474k  unknowns    3 s   -0.42%  12x
+
-    max_block=8  R = 0.503386 mOhm   253k  unknowns    1 s   -1.0%   22x
+    uniform          R = 0.508504 mOhm   5.58M unknowns   ~30-40 s
-    max_block=16 R = 0.497873 mOhm   213k  unknowns    1 s   -2.1%   26x
+    adaptive guard=4 R = 0.502793 mOhm    823k unknowns    ~8 s   -1.1%
+    adaptive guard=8 R = 0.506102 mOhm   1.74M unknowns   ~12 s   -0.47%
 
     uniform-limit check (max_block=1): rel diff 0.00e+00 vs production.
 
-Coarsening biases R low (coarse cells overestimate conductance where the
+This is the ADVERSARIAL case (features at 6 mm pitch everywhere, no
-field curves); a production version needs true 2:1 balancing + a guard
+smooth interior); big-pour boards coarsen far more aggressively. The
-band, and optionally one residual-driven refine pass, to push the
+residual bias is the first-order coarse-fine interface flux - phase 4
-max_block=4 accuracy to larger blocks. On big-pour boards (smooth
+(gradient-corrected fluxes / solution-adaptive refinement) is the
-interiors) the unknown ratios are far higher than on this geometry.
+lever if tighter accuracy per leaf is needed.
 """
 import sys
 import time
@@ -31,127 +29,49 @@ import numpy as np
 
 sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
 
-from fill_resistance import raster, solver
+from fill_resistance import quadtree, raster, solver
 from tests.util import NM, make_problem, strip_problem
 
-MAX_BLOCK = 32          # coarsest leaf = 32 x 32 fine cells
 
-
+def solve_adaptive(problem, h_mm, **kw):
-def build_leaves(mask, e1, e2, max_block=MAX_BLOCK):
-    """Greedy top-down coarsening. Returns (y0, x0, size) per leaf plus
-    an id_grid at fine resolution (-1 = empty)."""
-    ny, nx = mask.shape
-    pad_y = (-ny) % max_block
-    pad_x = (-nx) % max_block
-    m = np.pad(mask, ((0, pad_y), (0, pad_x)))
-    coars = m & ~np.pad(e1 | e2, ((0, pad_y), (0, pad_x)))
-    NY, NX = m.shape
-
-    from scipy import ndimage
-    id_grid = np.full((NY, NX), -1, dtype=np.int64)
-    covered = np.zeros((NY, NX), dtype=bool)
-    y0s, x0s, sizes = [], [], []
-    nid = 0
-    levels = []
-    red = coars.copy()
-    s = 1
-    while s < max_block:
-        red = red.reshape(red.shape[0] // 2, 2, red.shape[1] // 2, 2
-                          ).all(axis=(1, 3))
-        s *= 2
-        # grading: a size-s block must sit in an all-copper 3x3 block
-        # neighborhood at its own level, so leaf sizes step down smoothly
-        # toward boundaries (guard band + approximate 2:1 balance)
-        graded = ndimage.binary_erosion(red, np.ones((3, 3), dtype=bool))
-        levels.append((s, graded))
-    for s, allc in reversed(levels):
-        cov_k = covered.reshape(NY // s, s, NX // s, s).any(axis=(1, 3))
-        cand = allc & ~cov_k
-        ii, jj = np.nonzero(cand)
-        for i_, j_ in zip(ii, jj):
-            id_grid[i_ * s:(i_ + 1) * s, j_ * s:(j_ + 1) * s] = nid
-            y0s.append(i_ * s)
-            x0s.append(j_ * s)
-            sizes.append(s)
-            nid += 1
-        covered |= np.repeat(np.repeat(cand, s, axis=0), s, axis=1)
-    fi, fj = np.nonzero(m & ~covered)
-    n_fine = len(fi)
-    id_grid[fi, fj] = nid + np.arange(n_fine)
-    y0s.extend(fi.tolist())
-    x0s.extend(fj.tolist())
-    sizes.extend([1] * n_fine)
-    return (np.array(y0s), np.array(x0s), np.array(sizes),
-            id_grid[:ny, :nx])
-
-
-def leaf_edges(id_grid, sizes, sigma):
-    """All leaf-leaf face conductances, vectorized: count shared fine
-    faces per leaf pair (= overlap length w), g = sigma * w / mean(sa, sb).
-    Uniform limit: w = 1, sizes 1 -> g = sigma (identical to production)."""
-    aa, bb, ww = [], [], []
-    n = len(sizes)
-    for sl_a, sl_b in ((np.s_[:, :-1], np.s_[:, 1:]),
-                       (np.s_[:-1, :], np.s_[1:, :])):
-        a = id_grid[sl_a].ravel()
-        b = id_grid[sl_b].ravel()
-        ok = (a >= 0) & (b >= 0) & (a != b)
-        key = a[ok] * n + b[ok]
-        uniq, counts = np.unique(key, return_counts=True)
-        ia = uniq // n
-        ib = uniq % n
-        g = sigma * counts / (0.5 * (sizes[ia] + sizes[ib]))
-        aa.append(ia)
-        bb.append(ib)
-        ww.append(g)
-    return (np.concatenate(aa), np.concatenate(bb), np.concatenate(ww))
-
-
-def solve_adaptive(problem, h_mm, max_block=MAX_BLOCK):
     stack = raster.rasterize_stack(problem, h_mm * NM)
     e1, e2 = raster.electrode_masks(stack, problem)
     t0 = time.perf_counter()
-    y0, x0, sizes, id_grid = build_leaves(stack.masks[0], e1[0], e2[0],
+    grid = quadtree.build_leaves(stack.masks[0],
-                                          max_block)
+                                 keep_fine=(e1[0] | e2[0]), **kw)
-    a, b, w = leaf_edges(id_grid, sizes, problem.sigma_s(0))
+    ia, ib, g = quadtree.leaf_edges(grid, problem.sigma_s(0))
     t_build = time.perf_counter() - t0
 
-    n = len(sizes)
+    state = np.ones(grid.n, dtype=np.uint8)
-    state = np.ones(n, dtype=np.uint8)
+    for e, code in ((e1[0], 2), (e2[0], 3)):
-    fine_ids = id_grid[e1[0]]
+        ids = grid.id_grid[e]
-    state[fine_ids[fine_ids >= 0]] = 2
+        state[ids[ids >= 0]] = code
-    fine_ids = id_grid[e2[0]]
+    edges = solver.Edges(a=ia, b=ib, w=g,
-    state[fine_ids[fine_ids >= 0]] = 3
+                         via_index=np.full(len(ia), -1, dtype=np.int32))
-
-    edges = solver.Edges(a=a, b=b, w=w,
-                         via_index=np.full(len(a), -1, dtype=np.int32))
     t0 = time.perf_counter()
-    A, rhs, idx = solver._assemble(state, edges, None)
+    A, rhs, _ = solver._assemble(state, edges, None)
     x, info = solver.solve_system(A, rhs)
     t_solve = time.perf_counter() - t0
 
-    V = np.zeros(n)
+    V = np.zeros(grid.n)
     V[state == 2] = 1.0
     V[state == 1] = x
-    Ie = w * (V[a] - V[b])
+    Ie = g * (V[ia] - V[ib])
-    sa, sb = state[a], state[b]
+    sa, sb = state[ia], state[ib]
     I1 = float(Ie[sa == 2].sum() - Ie[sb == 2].sum())
     I2 = float(Ie[sb == 3].sum() - Ie[sa == 3].sum())
-    R = 1.0 / (0.5 * (I1 + I2))
+    return 1.0 / (0.5 * (I1 + I2)), grid, info, t_build, t_solve
-    mismatch = abs(I1 - I2) / max(abs(I1), abs(I2))
-    return R, n, info, t_build, t_solve, mismatch
 
 
-# --- 1) uniform-limit correctness: max_block=1 must equal production ---
+# --- 1) uniform-limit correctness ---
 p = strip_problem(length=50, width=10, e_len=5)
 stack = raster.rasterize_stack(p, 0.25 * NM)
 e1, e2 = raster.electrode_masks(stack, p)
 res = solver.run_solve(p, stack, e1, e2, 1.0, contact_model="equipotential")
-R_u, n_u, *_ = solve_adaptive(p, 0.25, max_block=1)
+R_u, *_ = solve_adaptive(p, 0.25, max_block=1)
-print(f"uniform-limit check: production R={res.R_ohm:.12g}, "
+print(f"uniform-limit check: rel diff {abs(R_u / res.R_ohm - 1):.2e}")
-      f"prototype R={R_u:.12g}, rel diff {abs(R_u / res.R_ohm - 1):.2e}")
 
-# --- 2) the payoff case: 120x120 plate, 400 holes, h = 50 um ---
+# --- 2) feature-dense plate at h = 50 um ---
 holes = []
 for i in range(20):
     for j in range(20):
@@ -167,13 +87,14 @@ e1, e2 = raster.electrode_masks(stack, big)
 ref = solver.run_solve(big, stack, e1, e2, 1.0,
                        contact_model="equipotential")
 t_ref = time.perf_counter() - t0
-print(f"\nuniform 50um : R = {ref.R_ohm * 1e3:.6g} mOhm, "
+print(f"uniform : R = {ref.R_ohm * 1e3:.6g} mOhm, "
-      f"{ref.solve_info.n_unknowns} unknowns, {t_ref:.1f} s total "
+      f"{ref.solve_info.n_unknowns} unknowns, {t_ref:.1f} s "
       f"({ref.solve_info.method})")
 
-for mb in (4, 8, 16, 32):
+R_a, grid, info, t_b, t_s = solve_adaptive(big, 0.05)
-    R_a, n_a, info, t_b, t_s, mm = solve_adaptive(big, 0.05, max_block=mb)
+print(f"adaptive: R = {R_a * 1e3:.6g} mOhm, {info.n_unknowns} unknowns, "
-    print(f"adaptive mb={mb:2d}: R = {R_a * 1e3:.6g} mOhm, "
+      f"build {t_b:.1f} s + solve {t_s:.1f} s ({info.method})")
-          f"{info.n_unknowns:8d} unknowns, build {t_b:.1f} s + "
+print(f"rel diff {abs(R_a / ref.R_ohm - 1):.2e}, "
-          f"solve {t_s:.1f} s, rel diff {abs(R_a / ref.R_ohm - 1):.2e}, "
+      f"{ref.solve_info.n_unknowns / info.n_unknowns:.1f}x fewer unknowns, "
-          f"ratio {ref.solve_info.n_unknowns / info.n_unknowns:.0f}x")
+      f"max leaf {int(grid.size.max())} cells, balanced="
+      f"{quadtree.balanced(grid)}")