penalty relaxation + radius bisection, hex/random starts, hole-move basin hopping, stall-free final polish

-"""Baseline: a lattice of candidate points inside the container, spacing found by bisection so
-that at least n fit. Deliberately naive; scores well below the records. Beat it."""
+"""Equal circles in the unit disk: penalty relaxation (heavy-ball gradient descent on the overlap
+energy) with a bisection on the radius. Random and hexagonal seeds are optimised coarsely first
+and only promising ones are polished; the incumbent is perturbed and repolished (basin hopping)."""
 
-import itertools
 import math
+import random
+import time
 
-DIM = 2
 
+def _derived(x, y, n):
+    """Largest common radius the eval will derive from these centres (O(n^2), used sparingly)."""
+    r = min(1.0 - math.hypot(x[i], y[i]) for i in range(n))
+    d2min = 4.0 * r * r
+    for i in range(n):
+        xi, yi = x[i], y[i]
+        for j in range(i + 1, n):
+            dx = xi - x[j]
+            dy = yi - y[j]
+            d2 = dx * dx + dy * dy
+            if d2 < d2min:
+                d2min = d2
+    return min(r, 0.5 * math.sqrt(d2min))
 
-def _boundary(c):
-    return 1.0 - math.hypot(c[0], c[1])
 
+def _pairs(x, y, n, cut):
+    """Index pairs closer than cut, via a grid of cell size cut."""
+    if n <= 40:
+        return [(i, j) for i in range(n) for j in range(i + 1, n)]
+    grid = {}
+    inv = 1.0 / cut
+    for i in range(n):
+        grid.setdefault((int(math.floor(x[i] * inv)), int(math.floor(y[i] * inv))), []).append(i)
+    out = []
+    c2 = cut * cut
+    for (gx, gy), members in grid.items():
+        for ox in (-1, 0, 1):
+            for oy in (-1, 0, 1):
+                key = (gx + ox, gy + oy)
+                if key < (gx, gy) or key not in grid:
+                    continue
+                others = grid[key]
+                same = key == (gx, gy)
+                for i in members:
+                    xi, yi = x[i], y[i]
+                    for j in others:
+                        if same and j <= i:
+                            continue
+                        dx = xi - x[j]
+                        dy = yi - y[j]
+                        if dx * dx + dy * dy < c2:
+                            out.append((i, j))
+    return out
 
-def _weight(i):
-    return 1.0
 
+def _relax(x, y, n, r, iters, tol, deadline, stall=0.7):
+    """Heavy-ball gradient descent on the overlap penalty for a fixed target radius r.
+    Returns True once the worst overlap is below tol (relative to r)."""
+    two_r = 2.0 * r
+    tr2 = two_r * two_r
+    wall = 1.0 - r
+    eta, mu = 0.35, 0.6
+    vx = [0.0] * n
+    vy = [0.0] * n
+    pairs = None
+    tol *= r
+    prev = float("inf")
+    for it in range(iters):
+        if pairs is None or it % 12 == 0:
+            pairs = _pairs(x, y, n, two_r * 1.3)
+        fx = [0.0] * n
+        fy = [0.0] * n
+        worst = 0.0
+        for i, j in pairs:
+            dx = x[i] - x[j]
+            dy = y[i] - y[j]
+            d2 = dx * dx + dy * dy
+            if d2 < tr2:
+                d = math.sqrt(d2)
+                if d < 1e-12:
+                    dx, dy, d = 1e-6, 0.0, 1e-6
+                ov = two_r - d
+                if ov > worst:
+                    worst = ov
+                k = ov / d
+                fx[i] += k * dx
+                fy[i] += k * dy
+                fx[j] -= k * dx
+                fy[j] -= k * dy
+        for i in range(n):
+            xi, yi = x[i], y[i]
+            rho = math.hypot(xi, yi)
+            if rho > wall:
+                ov = rho - wall
+                if ov > worst:
+                    worst = ov
+                k = ov / rho
+                fx[i] -= k * xi
+                fy[i] -= k * yi
+            v = mu * vx[i] + eta * fx[i]
+            w = mu * vy[i] + eta * fy[i]
+            vx[i] = v
+            vy[i] = w
+            x[i] = xi + v
+            y[i] = yi + w
+        if worst < tol:
+            return True
+        if (it & 15) == 15:
+            # stalled: overlaps are not shrinking geometrically, treat the target as infeasible
+            if worst > stall * prev:
+                return False
+            prev = worst
+            if time.perf_counter() > deadline:
+                break
+    return False
 
-def _lattice(n, r):
-    """Cubic lattice points at spacing 2r whose distance to the boundary is at least r."""
-    lo, hi = ((-1.0,) * DIM, (1.0,) * DIM)
-    step = 2.0 * r
-    axes = []
-    for d in range(DIM):
-        k = int((hi[d] - lo[d]) / step) + 1
-        axes.append([lo[d] + r + i * step for i in range(k)])
-    pts = [p for p in itertools.product(*axes) if _boundary(p) >= r]
-    return pts
 
+def _optimise(x, y, n, deadline, r_start, precision, tol, stall=0.7):
+    """Bisection on the radius with penalty relaxation; returns (best_r, best_x, best_y)."""
+    r_ok = _derived(x, y, n)
+    best = (r_ok, x[:], y[:])
+    r_bad = None
+    target = max(r_ok * 1.001, r_start)
+    while time.perf_counter() < deadline:
+        ok = _relax(x, y, n, target, 600 if stall < 1 else 6000, tol, deadline, stall)
+        d = _derived(x, y, n)
+        if d > best[0]:
+            best = (d, x[:], y[:])
+        r_ok = max(r_ok, d)
+        if ok:
+            target = target * 1.02 if r_bad is None else 0.5 * (target + r_bad)
+        else:
+            r_bad = target
+            target = 0.5 * (r_ok + r_bad)
+        if r_bad is not None and (r_bad - r_ok) < precision * r_ok:
+            break
+    return best
 
+
+def _random_start(n, rng):
+    x, y = [], []
+    while len(x) < n:
+        a, b = rng.uniform(-1, 1), rng.uniform(-1, 1)
+        if a * a + b * b < 0.9:
+            x.append(a)
+            y.append(b)
+    return x, y
+
+
+def _hex_start(n, rng, r):
+    """Hexagonal lattice at spacing 2r, randomly rotated and shifted, n points nearest the centre."""
+    while True:
+        s = 2.0 * r
+        th = rng.uniform(0, math.pi / 3)
+        ct, st = math.cos(th), math.sin(th)
+        ox, oy = rng.uniform(-s, s), rng.uniform(-s, s)
+        pts = []
+        m = int(2.2 / s) + 2
+        for i in range(-m, m + 1):
+            for j in range(-m, m + 1):
+                px = ox + s * (i + 0.5 * j)
+                py = oy + s * (j * math.sqrt(3) / 2)
+                qx, qy = ct * px - st * py, st * px + ct * py
+                rho = math.hypot(qx, qy)
+                if rho <= 1.0 - r:
+                    pts.append((rho, qx, qy))
+        if len(pts) >= n:
+            pts.sort()
+            return [p[1] for p in pts[:n]], [p[2] for p in pts[:n]]
+        r *= 0.97
+
+
+def _to_hole(x, y, n, rng, r):
+    """Move the circle with the most slack to the emptiest spot found by random sampling."""
+    # slack: distance to the nearest other centre (the loosest circle is the best candidate to move)
+    loose, slack = 0, -1.0
+    for i in range(n):
+        xi, yi = x[i], y[i]
+        m = min(math.hypot(xi - x[j], yi - y[j]) for j in range(n) if j != i)
+        m = min(m, 2.0 * (1.0 - math.hypot(xi, yi)))
+        if m > slack:
+            loose, slack = i, m
+    bx, by, bv = x[loose], y[loose], -1.0
+    for _ in range(60):
+        a = rng.uniform(0, 2 * math.pi)
+        rr = math.sqrt(rng.random()) * (1 - r)
+        px, py = rr * math.cos(a), rr * math.sin(a)
+        v = min(math.hypot(px - x[j], py - y[j]) for j in range(n) if j != loose)
+        v = min(v, 2.0 * (1.0 - math.hypot(px, py)))
+        if v > bv:
+            bx, by, bv = px, py, v
+    x[loose], y[loose] = bx, by
+
+
 def pack(n, time_budget, seed):
-    # treat every object as the largest one when choosing the lattice spacing
-    wmax = max(_weight(i + 1) for i in range(n))
-    # find a feasible spacing by halving, then bisect between it and the last infeasible one
-    a = 1.0
-    while len(_lattice(n, a * wmax)) < n and a > 1e-9:
-        a /= 2
-    b = 2 * a
-    for _ in range(40):
-        m = (a + b) / 2
-        if len(_lattice(n, m * wmax)) >= n:
-            a = m
+    rng = random.Random(seed)
+    final = time.perf_counter() + time_budget * 0.93
+    end = final - time_budget * 0.06   # keep the tail for a patient polish of the incumbent
+    # rough radius estimate from hexagonal density, used to seed targets and the hex lattice
+    r_est = math.sqrt(0.9069 / n) * (1 - 0.9 / math.sqrt(n))
+    big = n >= 100
+    best = (0.0, None, None)
+    k = 0
+    while time.perf_counter() < end:
+        now = time.perf_counter()
+        slot = min(end, now + max((end - now) / 4, 0.05))
+        hop = best[1] is not None and k >= 2 and rng.random() < (0.6 if big else 0.5)
+        if hop:
+            x, y = best[1][:], best[2][:]
+            u = rng.random()
+            if big and u < 0.33:
+                # shake: small gaussian kick to every centre
+                sig = best[0] * rng.uniform(0.05, 0.3)
+                for i in range(n):
+                    x[i] += rng.gauss(0, sig)
+                    y[i] += rng.gauss(0, sig)
+            elif u < 0.66:
+                _to_hole(x, y, n, rng, best[0])
+            else:
+                # relocate one to three circles at random
+                for _ in range(1 + int(rng.random() * min(3, n))):
+                    i = rng.randrange(n)
+                    a = rng.uniform(0, 2 * math.pi)
+                    rr = rng.uniform(0, 1 - best[0])
+                    x[i], y[i] = rr * math.cos(a), rr * math.sin(a)
+            r_start = best[0] * 0.99
+        elif k % 2 == 0 or big:
+            x, y = _hex_start(n, rng, r_est)
+            r_start = r_est * 0.9
         else:
-            b = m
-    pts = _lattice(n, a * wmax)
-    pts.sort(key=lambda p: -_boundary(p))   # keep the most interior points
-    return [tuple(p) for p in pts[:n]]
+            x, y = _random_start(n, rng)
+            r_start = r_est * 0.85
+        # coarse pass; polish only if it is within reach of the incumbent
+        res = _optimise(x, y, n, slot, r_start, 2e-3, 1e-5)
+        if res[0] > best[0] * 0.995:
+            x, y = res[1][:], res[2][:]
+            res = _optimise(x, y, n, min(end, slot + (end - now) / 4), res[0] * 0.999, 1e-7, 1e-9)
+        if res[0] > best[0]:
+            best = res
+        k += 1
+        if n <= 2:
+            break
+    if best[1] is None:
+        x, y = _hex_start(n, rng, r_est)
+        return list(zip(x, y))
+    # stall detection speeds the search but can leave the last digits; polish without it
+    res = _optimise(best[1][:], best[2][:], n, final, best[0] * 0.999, 1e-9, 1e-11, stall=1.0)
+    if res[0] > best[0]:
+        best = res
+    return list(zip(best[1], best[2]))