#!/usr/bin/env python3
# dn_router.py — DIMENSIONAL KURAMOTO LANE ROUTER (HDGL substrate)
#
# The realization: a PCB "crossing" is a 2D artifact. The HDGL Dn(r) operator
#   Dn(r) = sqrt(phi * F_n * 2^n * P_n * Omega) * r^((n+1)/8)
# lifts a coordinate into dimension n. At fixed r, the 8 dims form disjoint
# amplitude bands -> two nets that must cross in-plane DON'T cross if they live
# in different dim_n (= different copper layer / via-depth). Within one dim, the
# Kuramoto phase field (compute_derivatives_Dn) relaxes angular lanes to mutual
# separation. Conflicts that angular relaxation can't clear promote to dim_n+1.
#
# Two solvers, ONE substrate, checking each other:
#   (A) kuramoto_relax  — analog: settle phases, promote on residual conflict.
#   (B) exact_assign    — combinatorial: interval-graph / greedy-optimal lanes
#                          per dimension; ground-truth verifier for (A).
#
# This module is the LANE ORACLE used by phyllo_v3: given a set of (native_angle,
# radius_in, radius_out) demands in a cell, it returns for each a
# (dim_n, lane_angle, [bend]) assignment that is guaranteed conflict-free, using
# as few dimensions (layers) as the demand forces.

import math
PHI=1.6180339887498948
FIB=[1,1,2,3,5,8,13,21]; PRM=[2,3,5,7,11,13,17,19]
def delta(n,beta=0.0):
    # X+1=0 bridge: phase-advance coupled to a magnitude (layer) transition.
    # delta(i) = |cos(pi*beta*phi)| * ln(P_n) / phi^(n+beta)
    if n<1: n=1
    return abs(math.cos(math.pi*beta*PHI))*math.log(PRM[(n-1)%8])/(PHI**(n+beta))

def Dn(n,r,om=1.0):
    if n<1 or n>8: return 0.0
    i=n-1; k=(n+1)/8.0
    return math.sqrt(PHI*FIB[i]*(2.0**n)*PRM[i]*om)*(abs(r)**k)

# ---------------------------------------------------------------------------
# A track demand: a radial run from (r0) to (r1) that must terminate at angle
# ~theta (the pad direction). It owns an angular lane; two demands in the SAME
# dimension conflict iff their radial spans overlap AND their lanes are closer
# than DMIN. Different dimensions never conflict (orthogonal bands).
# ---------------------------------------------------------------------------
class Demand:
    __slots__=('idx','theta','r0','r1','pin_angle','net','dim','lane','bend')
    def __init__(self,idx,theta,r0,r1,net,pin_angle=None):
        self.idx=idx; self.theta=theta%(2*math.pi)
        self.r0=min(r0,r1); self.r1=max(r0,r1)
        self.net=net; self.pin_angle=pin_angle   # if set, lane MUST equal this
        self.dim=1; self.lane=self.theta; self.bend=None

def _sep(a,b):
    d=abs(a-b)%(2*math.pi); return min(d,2*math.pi-d)
def _radial_overlap(a,b,pad=0.30):
    if a.net==b.net: return False      # same copper never self-conflicts
    return not (a.r1+pad<b.r0 or b.r1+pad<a.r0)

# ===========================================================================
# (B) EXACT combinatorial assignment — ground truth.
# Per dimension, this is INTERVAL GRAPH COLORING on the angular circle with a
# hard min-separation. We assign lanes greedily in a globally-sorted order that
# is provably optimal for unit-demand interval separation (earliest-deadline).
# Demands that cannot fit the current dimension's angular budget promote to the
# next dimension. Returns True if all placed within max_dim dims.
# ===========================================================================
def exact_assign(demands, DMIN=0.075, max_dim=8):
    # group by radial-overlap conflict; within each conflict clique, lanes must
    # differ by >= DMIN. Model as: sort by native angle, sweep, pack into the
    # circle; if the circle is full for this dim, bump remainder to next dim.
    for d in demands: d.dim=1; d.lane=d.theta; d.bend=None
    remaining=sorted(demands,key=lambda d:d.theta)
    dim=1
    while remaining and dim<=max_dim:
        placed=[]; overflow=[]
        # taken lanes with their radial spans, per this dimension
        taken=[]  # (lane, r0, r1)
        for d in remaining:
            want = d.pin_angle if d.pin_angle is not None else d.theta
            want%=2*math.pi
            lane=_find_lane(want,d.r0,d.r1,taken,DMIN,d.net,pinned=(d.pin_angle is not None))
            if lane is None:
                overflow.append(d)
            else:
                d.dim=dim; d.lane=lane; taken.append((lane,d.r0,d.r1,d.net)); placed.append(d)
        remaining=overflow
        dim+=1
    ok = not remaining
    return ok, (dim-1)

def _find_lane(want,r0,r1,taken,DMIN,net,pinned=False,step=0.012,maxk=260):
    def free(cm):
        for (l,tr0,tr1,tnet) in taken:
            if tnet==net: continue
            if not (r1+0.30<tr0 or tr1+0.30<r0):
                if _sep(cm,l)<DMIN: return False
        return True
    if pinned:
        return want if free(want) else None
    for k in range(maxk):
        for cand in ((want+k*step)%(2*math.pi),(want-k*step)%(2*math.pi)):
            if free(cand): return cand
    return None

# ===========================================================================
# (A) KURAMOTO relaxation — analog solver on the HDGL substrate.
# Each demand is an oscillator with phase = lane. Same-dimension demands whose
# radial spans overlap are neighbors; coupling K*sin(phi_j - phi_i) with a
# repulsive sign pushes overlapping lanes apart (anti-phase attractor at the
# separation set). We integrate the same dynamics as compute_derivatives_Dn
# (phase + dimensional amplitude), and when a pair stays under DMIN after
# settling, the higher-index one promotes dim -> dim+1 (the Dn ladder catches
# it in an orthogonal band). Returns (ok, dims_used, iters).
# ===========================================================================
def kuramoto_relax(demands, DMIN=0.075, K=0.9, dt=0.05, iters=600, max_dim=8, om=1.0):
    for d in demands: d.dim=1; d.lane=d.theta; d.bend=None
    phase={d.idx:d.theta for d in demands}
    velo ={d.idx:0.0 for d in demands}
    dimv ={d.idx:1 for d in demands}
    Dnv  ={d.idx:Dn(1,_r_of(d),om) for d in demands}
    home ={d.idx:d.theta for d in demands}   # weak pull to native angle
    pin  ={d.idx:d.pin_angle for d in demands}
    by=lambda: _neighbors(demands,dimv)
    for it in range(iters):
        nb=by()
        newp={}; newv={}
        for d in demands:
            i=d.idx; ph=phase[i]
            # repulsive phase coupling within same dimension + overlap
            force=0.0; nn=0
            for j in nb[i]:
                dphi=phase[j]-ph
                s=math.sin(dphi)
                # repulsive: push apart when too close; scale by Dn coupling
                gap=_sep(phase[j],ph)
                w = (1.0/(0.15+gap))               # stronger when closer
                Dn_factor = Dnv[j]/(1.0+abs(Dnv[i]))
                force += -K*w*s*(0.5+0.5*Dn_factor)
                nn+=1
            # home pull (keep near native so routes stay short) unless pinned
            if pin[i] is not None:
                force += -2.0*(((ph-pin[i]+math.pi)%(2*math.pi))-math.pi)
            else:
                force += -0.06*(((ph-home[i]+math.pi)%(2*math.pi))-math.pi)
            v=velo[i]+dt*force
            v*=math.exp(-0.05*dt)                   # LAMBDA damping
            newv[i]=v; newp[i]=(ph+dt*v)%(2*math.pi)
        phase=newp; velo=newv
        # periodic promotion: unresolved overlaps kick to next dim
        if it%40==39:
            promoted=_promote(demands,phase,dimv,Dnv,DMIN,max_dim,om)
            if not promoted and _all_clear(demands,phase,dimv,DMIN):
                break
    for d in demands:
        d.lane=phase[d.idx]; d.dim=dimv[d.idx]
    ok=_all_clear(demands,phase,dimv,DMIN)
    return ok, max(dimv.values()), it+1

def _r_of(d): return 0.5*(d.r0+d.r1)/ max(1.0,d.r1)  # normalized radial in [0,1)-ish
def _neighbors(demands,dimv):
    nb={d.idx:[] for d in demands}
    for a in range(len(demands)):
        for b in range(a+1,len(demands)):
            A,B=demands[a],demands[b]
            if dimv[A.idx]==dimv[B.idx] and _radial_overlap(A,B):
                nb[A.idx].append(B.idx); nb[B.idx].append(A.idx)
    return nb
def _all_clear(demands,phase,dimv,DMIN):
    for a in range(len(demands)):
        for b in range(a+1,len(demands)):
            A,B=demands[a],demands[b]
            if dimv[A.idx]==dimv[B.idx] and _radial_overlap(A,B):
                if _sep(phase[A.idx],phase[B.idx])<DMIN: return False
    return True
def _promote(demands,phase,dimv,Dnv,DMIN,max_dim,om):
    moved=False
    for a in range(len(demands)):
        for b in range(a+1,len(demands)):
            A,B=demands[a],demands[b]
            if dimv[A.idx]==dimv[B.idx] and _radial_overlap(A,B) \
               and _sep(phase[A.idx],phase[B.idx])<DMIN:
                # promote the one with the smaller angular room / higher idx
                victim = B if B.pin_angle is None else (A if A.pin_angle is None else B)
                if dimv[victim.idx]<max_dim:
                    dimv[victim.idx]+=1
                    Dnv[victim.idx]=Dn(dimv[victim.idx],_r_of(victim),om)
                    moved=True
    return moved

# ===========================================================================
# Cross-check harness
# ===========================================================================
def route(demands, DMIN=0.075, prefer='kuramoto', verify=True):
    import copy
    dem_k=[_clone(d) for d in demands]
    dem_e=[_clone(d) for d in demands]
    okE,dimsE=exact_assign(dem_e,DMIN)
    okK,dimsK,itK=kuramoto_relax(dem_k,DMIN)
    result = dem_k if (prefer=='kuramoto' and okK) else dem_e
    used = dimsK if result is dem_k else dimsE
    ok = okK if result is dem_k else okE
    info=dict(exact_ok=okE,exact_dims=dimsE,kur_ok=okK,kur_dims=dimsK,kur_iters=itK,
              chosen='kuramoto' if result is dem_k else 'exact',dims=used,ok=ok)
    # copy assignment back onto the originals
    for src,dst in zip(result,demands):
        dst.dim=src.dim; dst.lane=src.lane; dst.bend=src.bend
    return ok,info
def _clone(d):
    e=Demand(d.idx,d.theta,d.r0,d.r1,d.net,d.pin_angle); return e

# ---------------------------------------------------------------------------
if __name__=='__main__':
    import random
    random.seed(1)
    # stress: 24 demands, many radially-overlapping, some pinned to same pole
    dem=[]
    for i in range(24):
        th=random.uniform(0,2*math.pi); r0=random.uniform(3,8); r1=random.uniform(12,21)
        pin = 0.74 if i in (0,1,2,3) else None   # 4 taps pinned near one pole
        dem.append(Demand(i,th,r0,r1,'N%d'%i,pin))
    ok,info=route(dem)
    print("Dn ladder @r=1:",["%.1f"%Dn(n,1.0) for n in range(1,9)])
    print("exact:  ok=%s dims=%d"%(info['exact_ok'],info['exact_dims']))
    print("kuramoto: ok=%s dims=%d iters=%d"%(info['kur_ok'],info['kur_dims'],info['kur_iters']))
    print("chosen=%s final ok=%s dims=%d"%(info['chosen'],info['ok'],info['dims']))
    from collections import Counter
    print("dim usage:",dict(Counter(d.dim for d in dem)))
    # verify no residual conflict in chosen assignment
    bad=0
    for a in range(len(dem)):
        for b in range(a+1,len(dem)):
            A,B=dem[a],dem[b]
            if A.dim==B.dim and _radial_overlap(A,B) and _sep(A.lane,B.lane)<0.075: bad+=1
    print("residual conflicts in chosen:",bad)
