#!/usr/bin/env python3
"""
hn_hoffman_bound.py — Hadwiger-Nelson Hoffman bound via spectral method.

Computes χ(G) ≥ 1 - λ_max / λ_min for a unit-distance graph G.

Built-in graphs:
  - Moser spindle (7 vertices, χ=4)
  - de Grey 1581-vertex graph (χ=5) using coordinates from arXiv:1804.02385v3 §5.1

Usage:
  python3 hn_hoffman_bound.py [--graph moser|degrey]
"""
import sys, json, math, hashlib, time
from pathlib import Path

HERE = Path(__file__).resolve().parent
OUT_DIR = HERE.parent / ".openresearch" / "artifacts"
OUT_DIR.mkdir(parents=True, exist_ok=True)

s3 = math.sqrt(3)
s11 = math.sqrt(11)
s33 = math.sqrt(33)
s12 = 2 * s3

# ── Coordinate set S from de Grey §5.1 ────────────────────────────

S_COORDS = [
    (0, 0),
    (1/3, 0),
    (1, 0),
    (2, 0),
    ((s33 - 3)/6, 0),
    (1/2, 1/s12),
    (1, 1/s3),
    (3/2, s3/2),
    (7/6, s11/6),
    (1/6, (s12 - s11)/6),
    (5/6, (s12 - s11)/6),
    (2/3, (s11 - s3)/6),
    (2/3, (3*s3 - s11)/6),
    (s33/6, 1/s12),
    ((s33 + 3)/6, 1/s3),
    ((s33 + 1)/6, (3*s3 - s11)/6),
    ((s33 - 1)/6, (3*s3 - s11)/6),
    ((s33 + 1)/6, (s11 - s3)/6),
    ((s33 - 1)/6, (s11 - s3)/6),
    ((s33 - 2)/6, (2*s3 - s11)/6),
    ((s33 - 4)/6, (2*s3 - s11)/6),
    ((s33 + 13)/12, (s11 - s3)/12),
    ((s33 + 11)/12, (s3 + s11)/12),
    ((s33 + 9)/12, (s11 - s3)/4),
    ((s33 + 9)/12, (3*s3 + s11)/12),
    ((s33 + 7)/12, (s3 + s11)/12),
    ((s33 + 7)/12, (3*s3 - s11)/12),
    ((s33 + 5)/12, (5*s3 - s11)/12),
    ((s33 + 5)/12, (s11 - s3)/12),
    ((s33 + 3)/12, (3*s11 - 5*s3)/12),
    ((s33 + 3)/12, (s3 + s11)/12),
    ((s33 + 3)/12, (3*s3 - s11)/12),
    ((s33 + 1)/12, (s11 - s3)/12),
    ((s33 - 1)/12, (3*s3 - s11)/12),
    ((s33 - 3)/12, (s11 - s3)/12),
    ((15 - s33)/12, (s11 - s3)/4),
    ((15 - s33)/12, (7*s3 - 3*s11)/12),
    ((13 - s33)/12, (3*s3 - s11)/12),
    ((11 - s33)/12, (s11 - s3)/12),
]

def rotate(pt, angle):
    c, s = math.cos(angle), math.sin(angle)
    return (pt[0]*c - pt[1]*s, pt[0]*s + pt[1]*c)

def translate(pt, dx, dy):
    return (pt[0]+dx, pt[1]+dy)

def negate_y(pt):
    return (pt[0], -pt[1])

def is_unit(p, q, eps=1e-9):
    d2 = (p[0]-q[0])**2 + (p[1]-q[1])**2
    return abs(d2 - 1.0) < eps

def build_adjacency(points, eps=1e-9):
    n = len(points)
    adj = [[] for _ in range(n)]
    for i in range(n):
        pi = points[i]
        for j in range(i+1, n):
            if is_unit(pi, points[j], eps):
                adj[i].append(j)
                adj[j].append(i)
    return adj

def remove_duplicates(points, eps=1e-9):
    uniq = []
    for p in points:
        if not any(abs(p[0]-q[0])<eps and abs(p[1]-q[1])<eps for q in uniq):
            uniq.append(p)
    return uniq

def build_degrey_1581():
    """Construct de Grey's 1581-vertex graph per §5.1."""
    # Step 1: S_a — rotate S by 60°, reflect y
    sa = set()
    angle60 = math.pi / 3
    angle2 = 2 * math.asin(1/4)
    angle_y = math.pi/2 + math.asin(1/8)
    
    for x, y in S_COORDS:
        for k in range(6):
            pt = rotate((x, y), k * angle60)
            sa.add((round(pt[0], 12), round(pt[1], 12)))
            pt_r = negate_y(pt)
            sa.add((round(pt_r[0], 12), round(pt_r[1], 12)))
    
    sa_list = list(sa)
    
    # Step 3: S_b = S_a rotated by 2*arcsin(1/4)
    sb_list = [(round(rotate(p, angle2)[0], 12),
                round(rotate(p, angle2)[1], 12)) for p in sa_list]
    
    # Step 4: Y = union minus (1/3, 0) and (-1/3, 0)
    y_set = set(sa_list) | set(sb_list)
    y_set.discard((round(1/3, 12), 0.0))
    y_set.discard((round(-1/3, 12), 0.0))
    y_list = list(y_set)
    
    # Step 5-6: Rotate Y about (-2, 0)
    ya = [translate(rotate(translate(p, 2, 0), angle_y), -2, 0) for p in y_list]
    yb = [translate(rotate(translate(p, 2, 0), math.pi - angle_y), -2, 0) for p in y_list]
    ya = [(round(x, 12), round(y, 12)) for x, y in ya]
    yb = [(round(x, 12), round(y, 12)) for x, y in yb]
    
    # Step 7: G = Y_a ∪ Y_b
    all_pts = remove_duplicates(ya + yb)
    
    print(f"  Total vertices before dedup: {len(ya)+len(yb)}")
    print(f"  After dedup: {len(all_pts)}")
    
    # Build adjacency
    adj = build_adjacency(all_pts)
    n_edges = sum(len(nbrs) for nbrs in adj) // 2
    
    return adj, all_pts, f"De Grey G ({len(all_pts)} vertices, {n_edges} edges)"

# ── Moser spindle ──────────────────────────────────────────────────

def moser_spindle():
    s3 = math.sqrt(3)
    pts = [(0,0), (1,0), (0.5,s3/2), (-0.5,s3/2), (-1,0), (-0.5,-s3/2), (0.5,-s3/2)]
    adj = build_adjacency(pts)
    n_edges = sum(len(nbrs) for nbrs in adj) // 2
    return adj, pts, f"Moser spindle (7 vertices, {n_edges} edges)"

# ── Hoffman bound ──────────────────────────────────────────────────

def hoffman_bound(adj):
    try:
        import numpy as np
    except ImportError:
        print("ERROR: numpy required")
        sys.exit(1)
    
    n = len(adj)
    A = np.zeros((n, n), dtype=np.float64)
    for i in range(n):
        for j in adj[i]:
            A[i, j] = 1.0
    
    eigenvals = np.linalg.eigvalsh(A)
    lambda_max = eigenvals[-1]
    lambda_min = eigenvals[0]
    hb = 1.0 - lambda_max / lambda_min if lambda_min < 0 else float('inf')
    
    return float(lambda_max), float(lambda_min), float(hb)

# ── Main ───────────────────────────────────────────────────────────

def main():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("--graph", choices=["moser", "degrey"], default="degrey")
    args = parser.parse_args()
    
    t0 = time.time()
    
    if args.graph == "moser":
        adj, pts, desc = moser_spindle()
    else:
        adj, pts, desc = build_degrey_1581()
    
    n_vertices = len(adj)
    n_edges = sum(len(nbrs) for nbrs in adj) // 2
    
    print(f"Graph: {desc}")
    if n_vertices < 3:
        print("  Too few vertices.")
        return
    
    lambda_max, lambda_min, hb = hoffman_bound(adj)
    chi_lower = math.ceil(hb) if math.isfinite(hb) else None
    
    print(f"\nSpectral analysis:")
    print(f"  λ_max = {lambda_max:.6f}")
    print(f"  λ_min = {lambda_min:.6f}")
    print(f"  λ_max/|λ_min| = {lambda_max / abs(lambda_min):.6f}")
    if chi_lower:
        print(f"\n  Hoffman bound: χ ≥ {hb:.6f}")
        print(f"  Rounded up:   χ ≥ {chi_lower}")
    
    elapsed = time.time() - t0
    result = {
        "experiment": "hn_hoffman_bound",
        "graph": args.graph,
        "description": desc,
        "n_vertices": n_vertices,
        "n_edges": n_edges,
        "lambda_max": round(lambda_max, 12),
        "lambda_min": round(lambda_min, 12),
        "hoffman_bound": round(hb, 12) if math.isfinite(hb) else None,
        "chi_lower_bound": chi_lower,
        "elapsed_s": round(elapsed, 2),
        "timestamp": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
    }
    content = json.dumps(result, indent=2)
    result["sha256"] = hashlib.sha256(content.encode()).hexdigest()
    
    out_path = OUT_DIR / "hn_hoffman_bound.json"
    with open(out_path, "w") as f:
        json.dump(result, f, indent=2)
    print(f"\nResults → {out_path}")
    print(f"SHA-256: {result['sha256']}")
    print(f"Elapsed: {elapsed:.2f}s")

if __name__ == "__main__":
    main()