files / problems / fekete-s2 / compute / optimize.py
#!/usr/bin/env python3
"""Deterministic multi-start Riemannian descent for log-energy on S^2.
Seeds are fixed (Fibonacci / spiral / cube-plus-noise / RNG). Checkpoints
are written as JSON. A printed Table 3 value is beaten only if the
independent verifier agrees and the gap is larger than rounding.
"""
from __future__ import annotations
import argparse
import json
import math
from pathlib import Path
import numpy as np
from energy import log_energy, project_to_sphere
HERE = Path(__file__).resolve().parent
def fibonacci(n: int) -> np.ndarray:
# Deterministic spherical Fibonacci lattice.
i = np.arange(n, dtype=np.float64)
z = 1.0 - 2.0 * (i + 0.5) / n
r = np.sqrt(np.clip(1.0 - z * z, 0.0, 1.0))
phi = math.pi * (3.0 - math.sqrt(5.0)) * i
return np.column_stack([r * np.cos(phi), r * np.sin(phi), z])
def spiral(n: int) -> np.ndarray:
# Rakhmanov–Saff–Zhou generalized spiral.
i = np.arange(1, n + 1, dtype=np.float64)
z = -1.0 + 2.0 * (i - 1.0) / (n - 1.0) if n > 1 else np.array([1.0])
r = np.sqrt(np.clip(1.0 - z * z, 0.0, 1.0))
# cumulative longitude ~ 3.6 / sqrt(N) style
theta = np.zeros(n)
if n > 2:
mid = z[1:-1]
dtheta = 3.6 / np.sqrt(np.clip(n * (1.0 - mid * mid), 1e-15, None))
theta[1:-1] = np.cumsum(dtheta)
return np.column_stack([r * np.cos(theta), r * np.sin(theta), z])
def rng_points(n: int, seed: int) -> np.ndarray:
rng = np.random.default_rng(seed)
g = rng.standard_normal((n, 3))
return project_to_sphere(g)
def energy_and_riemannian_grad(pts: np.ndarray):
"""E and the tangent-space gradient of E at unit points."""
n = len(pts)
diff = pts[:, None, :] - pts[None, :, :]
sq = np.einsum("ijk,ijk->ij", diff, diff)
np.fill_diagonal(sq, np.inf)
inv = 1.0 / sq
grad = -(diff * inv[:, :, None]).sum(axis=1)
i, j = np.triu_indices(n, k=1)
e = float(-0.5 * np.log(sq[i, j]).sum())
grad = grad - (np.sum(grad * pts, axis=1))[:, None] * pts
return e, grad
def descend(pts: np.ndarray, max_iter: int = 4000, tol: float = 1e-14) -> np.ndarray:
pts = project_to_sphere(pts)
step = 0.25
e, g = energy_and_riemannian_grad(pts)
for _ in range(max_iter):
gnorm = float(np.linalg.norm(g))
if gnorm < 1e-12:
break
improved = False
for _try in range(20):
cand = project_to_sphere(pts - step * g)
e2, g2 = energy_and_riemannian_grad(cand)
if e2 < e - 1e-16:
pts, e, g = cand, e2, g2
step = min(step * 1.2, 2.0)
improved = True
break
step *= 0.5
if not improved:
break
if abs(e - e2) < tol and gnorm < 1e-10:
break
return pts
def square_antiprism(n: int = 8, h: float = 0.55) -> np.ndarray:
# N=8 putative: two squares, one rotated 45°, at z=±h.
assert n == 8
r = math.sqrt(max(1.0 - h * h, 1e-12))
pts = []
for sign, rot in ((1.0, 0.0), (-1.0, math.pi / 4.0)):
for k in range(4):
ang = rot + math.pi * k / 2.0
pts.append([r * math.cos(ang), r * math.sin(ang), sign * h])
return np.asarray(pts)
def seeds_for(n: int, n_random: int, base_seed: int):
yield "fibonacci", fibonacci(n)
yield "spiral", spiral(n)
if n == 8:
for k, h in enumerate(np.linspace(0.35, 0.70, 8)):
yield f"antiprism{k}", square_antiprism(8, float(h))
try:
from known import KNOWN
if n in KNOWN:
yield "known", KNOWN[n][1]()
except Exception:
pass
rathbun = Path("/tmp/fekete-data/rathbun") / f"log.3.{n}.80"
if rathbun.exists():
try:
from replay_published import parse_pari
pts, _ = parse_pari(rathbun)
yield "rathbun", pts
rng = np.random.default_rng(base_seed + 17 * n)
for k, eps in enumerate((1e-3, 3e-3, 1e-2, 3e-2)):
kick = project_to_sphere(pts + eps * rng.standard_normal(pts.shape))
yield f"kick{k}", kick
except Exception:
pass
for k in range(n_random):
yield f"rng{k}", rng_points(n, base_seed + 10007 * k + n)
def search(n: int, n_random: int = 24, base_seed: int = 4102, max_iter: int = 4000):
best_e = float("inf")
best_pts = None
best_src = None
log = []
for src, pts0 in seeds_for(n, n_random, base_seed):
pts = descend(pts0, max_iter=max_iter)
e = log_energy(pts)
log.append({"seed": src, "E": e})
if e < best_e:
best_e, best_pts, best_src = e, pts, src
return {
"N": n,
"E": best_e,
"best_seed": best_src,
"points": best_pts.tolist(),
"starts": log,
}
def main() -> None:
ap = argparse.ArgumentParser()
ap.add_argument("N", type=int, nargs="+")
ap.add_argument("--random", type=int, default=24)
ap.add_argument("--seed", type=int, default=4102)
ap.add_argument("--iters", type=int, default=4000)
ap.add_argument("--out", type=Path, default=HERE / "checkpoints")
args = ap.parse_args()
args.out.mkdir(parents=True, exist_ok=True)
table = json.loads((HERE / "ridgway2018.json").read_text())["globals"]
for n in args.N:
rec = search(n, n_random=args.random, base_seed=args.seed, max_iter=args.iters)
# Re-verify from stored points, independently of the descent energy.
e = log_energy(rec["points"])
rec["E"] = e
rec["E_verifier"] = e
published = table.get(str(n))
rec["published_ridgway2018"] = published
rec["delta_vs_published"] = None if published is None else e - published
path = args.out / f"N{n:03d}.json"
path.write_text(json.dumps(rec, indent=2) + "\n")
delta = rec["delta_vs_published"]
flag = ""
if delta is not None and delta < -1e-7:
flag = " CANDIDATE BEAT"
print(
f"N={n:3d} E={e:.12f} pub={published} "
f"delta={delta} seed={rec['best_seed']}{flag}"
)
if __name__ == "__main__":
main()