hillclimb

You are improving the current best solution for the problem below.

# Problem

{{description}}

# Objective

Optimize **{{metric_name}}** ({{direction}}).

# Current best solution

`solution.py` in your working directory — validation {{metric_name}}: **{{best_score}}**.

Its execution stdout (tail):
```
{{stdout_tail}}
```

{{evaluation_report}}
# Already tried from this solution (do NOT repeat)

{{sibling_summaries}}

{{live_experience}}
{{prior_ablations}}
# Instruction

{{ablation_cue}}
Make **exactly ONE measurable change** to improve validation {{metric_name}}: one engineered feature, one model swap, one hyperparameter adjustment, one validation-scheme fix — a single attributable modification. State your hypothesis as the first line of `notes.md`. Keep everything else identical so the score change is attributable to your change.

{{contract}}

"""Heilbronn (convex-hull normalised), n=13.

Symmetry-reduced smoothing homotopy: the max-min objective is replaced by a
log-sum-exp soft-min of |triangle area| / hull area, maximised with L-BFGS
(analytic gradients) while the sharpness beta is annealed upward.  Starts are
drawn in symmetry-reduced coordinates (mirror, 3-fold rotation, or none), the
symmetry is then released, and promising results get an exact epigraph polish.
"""
import os
import time
from itertools import combinations
from multiprocessing import Pool

import numpy as np
from scipy.optimize import minimize
from scipy.spatial import ConvexHull

try:
    from hillclimb import spaces
    P = spaces.params()
except Exception:  # noqa: BLE001
    P = {}
TIME_BUDGET = float(P.get("time_budget", 45.0))
BETA0 = float(P.get("beta_start", 300.0))
BETA_GROWTH = float(P.get("beta_growth", 3.0))
SYM_FRAC = float(P.get("sym_frac", 0.6))
EXPLORE_FRAC = float(P.get("explore_frac", 0.5))  # restarts, then basin hopping

N = 13
TRI = np.array(list(combinations(range(N), 3)))
I, J, K = TRI[:, 0], TRI[:, 1], TRI[:, 2]
NT = len(TRI)


def signed_areas(X):
    x, y = X[:, 0], X[:, 1]
    return 0.5 * ((x[J] - x[I]) * (y[K] - y[I]) - (x[K] - x[I]) * (y[J] - y[I]))


def hull_info(X):
    h = ConvexHull(X).vertices  # counter-clockwise in 2-D
    x, y = X[h, 0], X[h, 1]
    H = 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
    return h, H


def score(X):
    try:
        _, H = hull_info(X)
    except Exception:  # noqa: BLE001
        return 0.0
    return float(np.min(np.abs(signed_areas(X))) / H)


def area_grad_terms(X, coef):
    """Gradient of sum_i coef_i * A_i (signed areas) w.r.t. X."""
    x, y = X[:, 0], X[:, 1]
    G = np.zeros((N, 2))
    np.add.at(G[:, 0], I, coef * 0.5 * (y[J] - y[K]))
    np.add.at(G[:, 0], J, coef * 0.5 * (y[K] - y[I]))
    np.add.at(G[:, 0], K, coef * 0.5 * (y[I] - y[J]))
    np.add.at(G[:, 1], I, coef * 0.5 * (x[K] - x[J]))
    np.add.at(G[:, 1], J, coef * 0.5 * (x[I] - x[K]))
    np.add.at(G[:, 1], K, coef * 0.5 * (x[J] - x[I]))
    return G


def hull_grad(X, h):
    G = np.zeros((N, 2))
    x, y = X[h, 0], X[h, 1]
    G[h, 0] = 0.5 * (np.roll(y, -1) - np.roll(y, 1))
    G[h, 1] = 0.5 * (np.roll(x, 1) - np.roll(x, -1))
    return G


def neg_softmin(z, beta, A):
    """-softmin_beta(|A_i|/H) and gradient, with X = (A @ z) when A given."""
    X = (A @ z if A is not None else z).reshape(N, 2)
    try:
        h, H = hull_info(X)
    except Exception:  # noqa: BLE001
        return 0.0, np.zeros_like(z)
    a = signed_areas(X)
    s = np.sign(a)
    r = np.abs(a) / H
    m = r.min()
    e = np.exp(-beta * (r - m))
    S = e.sum()
    f = m - np.log(S / NT) / beta
    w = e / S
    G = area_grad_terms(X, w * s / H) - (np.dot(w, r) / H) * hull_grad(X, h)
    g = G.ravel()
    if A is not None:
        g = A.T @ g
    return -f, -g


def sym_matrix(kind, rng):
    """Linear map from reduced coordinates to the 26 point coordinates."""
    if kind == "mirror":
        # 6 mirrored pairs + 1 point on the axis (x = 0)
        A = np.zeros((2 * N, 13))
        for p in range(6):
            A[2 * (2 * p), 2 * p] = 1.0
            A[2 * (2 * p) + 1, 2 * p + 1] = 1.0
            A[2 * (2 * p + 1), 2 * p] = -1.0
            A[2 * (2 * p + 1) + 1, 2 * p + 1] = 1.0
        A[2 * 12 + 1, 12] = 1.0
        z = rng.uniform(-1, 1, 13)
        return A, z
    if kind == "rot3":
        # 4 orbits of the 3-fold rotation + the centre
        A = np.zeros((2 * N, 8))
        for o in range(4):
            for q in range(3):
                c, s = np.cos(2 * np.pi * q / 3), np.sin(2 * np.pi * q / 3)
                p = 3 * o + q
                A[2 * p, 2 * o] = c
                A[2 * p, 2 * o + 1] = -s
                A[2 * p + 1, 2 * o] = s
                A[2 * p + 1, 2 * o + 1] = c
        z = rng.uniform(-1, 1, 8)
        return A, z
    if kind == "point":
        # 6 antipodal pairs + the centre (kept slightly off-centre by release)
        A = np.zeros((2 * N, 12))
        for p in range(6):
            for d in range(2):
                A[2 * (2 * p) + d, 2 * p + d] = 1.0
                A[2 * (2 * p + 1) + d, 2 * p + d] = -1.0
        z = rng.uniform(-1, 1, 12)
        return A, z
    raise ValueError(kind)


def anneal(z, A, betas, maxiter):
    for beta in betas:
        res = minimize(neg_softmin, z, args=(beta, A), jac=True,
                       method="L-BFGS-B", options={"maxiter": maxiter})
        z = res.x
    return z


def polish(X, rounds=6):
    """Exact epigraph polish: max t s.t. s_i A_i >= t H, signs/hull frozen."""
    best, bs = X.copy(), score(X)
    for _ in range(rounds):
        X0 = best
        try:
            h, H0 = hull_info(X0)
        except Exception:  # noqa: BLE001
            break
        X0 = X0 / np.sqrt(H0)
        a0 = signed_areas(X0)
        s = np.sign(a0)
        t0 = np.abs(a0).min()
        act = np.abs(a0) < 2.0 * t0 + 1e-9  # only near-active triangles
        ia, ja, ka, sa = I[act], J[act], K[act], s[act]
        na = int(act.sum())
        rows = np.arange(na)

        def hull_area(X):
            x, y = X[h, 0], X[h, 1]
            return 0.5 * (np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))

        def cons(v):
            X = v[:-1].reshape(N, 2)
            x, y = X[:, 0], X[:, 1]
            a = 0.5 * ((x[ja] - x[ia]) * (y[ka] - y[ia])
                       - (x[ka] - x[ia]) * (y[ja] - y[ia]))
            return np.concatenate([sa * a - v[-1], [1.0 - hull_area(X)]])

        def cons_jac(v):
            X = v[:-1].reshape(N, 2)
            x, y = X[:, 0], X[:, 1]
            Jm = np.zeros((na + 1, 2 * N + 1))
            Jm[rows, 2 * ia] = sa * 0.5 * (y[ja] - y[ka])
            Jm[rows, 2 * ja] = sa * 0.5 * (y[ka] - y[ia])
            Jm[rows, 2 * ka] = sa * 0.5 * (y[ia] - y[ja])
            Jm[rows, 2 * ia + 1] = sa * 0.5 * (x[ka] - x[ja])
            Jm[rows, 2 * ja + 1] = sa * 0.5 * (x[ia] - x[ka])
            Jm[rows, 2 * ka + 1] = sa * 0.5 * (x[ja] - x[ia])
            Jm[:na, -1] = -1.0
            Jm[na, :-1] = -hull_grad(X, h).ravel()
            return Jm

        obj_g = np.zeros(2 * N + 1)
        obj_g[-1] = -1.0
        v0 = np.concatenate([X0.ravel(), [t0]])
        try:
            res = minimize(lambda v: -v[-1], v0, jac=lambda v: obj_g,
                           constraints=[{"type": "ineq", "fun": cons, "jac": cons_jac}],
                           method="SLSQP", options={"maxiter": 150, "ftol": 1e-14})
        except Exception:  # noqa: BLE001
            break
        Xn = res.x[:-1].reshape(N, 2)
        sn = score(Xn)
        if sn > bs + 1e-12:
            best, bs = Xn, sn
        else:
            break
    return best, bs


def worker(args):
    seed, budget, start = args
    rng = np.random.default_rng(seed)
    t_end = time.time() + budget
    betas = [BETA0 * BETA_GROWTH ** k for k in range(6)]
    best, bs = (None, -1.0) if start is None else (start.copy(), score(start))
    while time.time() < t_end:
        u = rng.random()
        if start is not None:
            # exploitation: basin-hop around the shared (pool-wide) incumbent
            _, H = hull_info(best)
            X = best / np.sqrt(H)
            X = X + 10 ** rng.uniform(-2.5, -1.0) * rng.standard_normal(X.shape)
            z = anneal(X.ravel(), None, betas[3:], 200)
        elif u < SYM_FRAC:
            kind = rng.choice(["mirror", "rot3", "point"], p=[0.5, 0.35, 0.15])
            A, z = sym_matrix(kind, rng)
            z = anneal(z, A, betas[:4], 150)
            X = (A @ z).reshape(N, 2)
            X = X + 1e-3 * rng.standard_normal(X.shape)  # release symmetry
            z = anneal(X.ravel(), None, betas[2:], 200)
        else:
            ang = rng.uniform(0, 2 * np.pi, N)
            rad = np.sqrt(rng.uniform(0, 1, N))
            X = np.c_[rad * np.cos(ang), rad * np.sin(ang)]
            z = anneal(X.ravel(), None, betas, 200)
        X = z.reshape(N, 2)
        sc = score(X)
        if sc > 0.9 * bs:
            X, sc = polish(X)
        if sc > bs:
            best, bs = X.copy(), sc
    return bs, best


def main():
    seed = int(os.environ.get("HILLCLIMB_REPLICATE_SEED", "0"))
    nproc = min(8, os.cpu_count() or 1)
    t_explore = EXPLORE_FRAC * TIME_BUDGET
    with Pool(nproc) as pool:
        explored = pool.map(worker, [(seed * 1000 + 17 * w + 1, t_explore, None) for w in range(nproc)])
        # exploitation: every worker basin-hops around one of the pool's top-3
        # incumbents (half of them around the very best) instead of its own
        ranked = sorted([r for r in explored if r[1] is not None], key=lambda r: -r[0])
        starts = [ranked[0 if w % 2 == 0 else min(1 + (w // 2) % 2, len(ranked) - 1)][1]
                  for w in range(nproc)]
        results = pool.map(worker, [(seed * 1000 + 17 * w + 7919, TIME_BUDGET - t_explore, starts[w])
                                    for w in range(nproc)])
    print("explore per-worker:", sorted(round(r[0], 6) for r in explored))
    bs, best = max(results, key=lambda r: r[0])
    best, bs = polish(best, rounds=20)
    h, H = hull_info(best)
    best = (best - best.mean(axis=0)) / np.sqrt(H)
    with open("submission.csv", "w") as f:
        f.write("id,x,y\n")
        for i, (x, y) in enumerate(best):
            f.write(f"{i},{float(x)!r},{float(y)!r}\n")
    print(f"best normalized min triangle area: {score(best):.12f}")
    print("per-worker:", sorted(round(r[0], 6) for r in results))


if __name__ == "__main__":
    main()