"""n33 — does the SAME edge bath produce dissipator + coherent tau->eps transfer?

FULL microscopic model (all Gaussian, exact):
  clocks: [tau_i, eps_j] = i delta_ij on ring of N nodes
  H_clocks = sum_i (1 + lam*h_i) eps_i        (h_i: static matter energy profile)
  per edge e: damped auxiliary a_e, H_aux = Delta a_e†a_e, damping kappa D[a_e]
  Q_e = (B tau)_e + i eta (B eps)_e
  H_sb = g sum_e (a_e† Q_e + Q_e† a_e)  =  sqrt(2) g sum_e (x_ae x_e + eta p_ae p_e)

CANDIDATE effective model (to derive & kill or confirm):
  gamma_eff = g^2 kappa / (Delta^2 + kappa^2/4),  Omega_eff = -g^2 Delta / (Delta^2 + kappa^2/4)
  H_LS = Omega_eff sum_e Q_e†Q_e = Omega_eff (tau^T L_G tau + eta^2 eps^T L_G eps) + const
  first moments:
    dtau/dt = -(gamma_eff*eta) L tau + 2 Omega_eff eta^2 L eps + lam*(h - hbar)
    deps/dt = -(gamma_eff*eta) L eps - 2 Omega_eff L tau
  NESS (non-zero modes):  eps_ss = -(2 lam Omega)/(Gamma^2 + 4 Omega^2 eta^2) L^+ (h-hbar)
GATES: full-vs-eff transient <1% in regime; eps_ss vs closed form <1e-5 (eff) and
  vs FULL NESS; force from same H_int (HF vs grad U) <1%; N3 <1%;
  controls: Delta=0 -> F~0 with Gamma>0 ; Delta -> -Delta flips sign ; scaling <5%.
NESS of full model = direct linear solve of A R + b = 0 (no integration)."""
import numpy as np

N=12; ETA=0.5; LAM=0.1
B=np.zeros((N,N))
for e in range(N): B[e,e]=1.0; B[e,(e+1)%N]=-1.0
LG=B.T@B; E=N  # edges

def full_drift(g,Delta,kappa):
    """R = (tau[N], eps[N], xa[E], pa[E]). Returns A, and source map for b."""
    n=2*N+2*E
    A=np.zeros((n,n))
    # Hamiltonian part: dR/dt = S dH/dR with S the symplectic form.
    # H = sum eps_i (source handled in b) + Delta/2 sum (xa^2+pa^2)  [a†a up to const]
    #   + sqrt2 g sum_e ( xa_e * (B tau)_e + eta * pa_e * (B eps)_e )
    # canonical pairs: (tau_i, eps_i), (xa_e, pa_e)
    # dtau/dt = dH/deps ; deps/dt = -dH/dtau ; dxa/dt = dH/dpa ; dpa/dt = -dH/dxa
    sq2g=np.sqrt(2)*g
    # dtau_i/dt += dH/deps_i = sq2*g*eta * (B^T pa)_i
    A[0:N,2*N+E:2*N+2*E]+= sq2g*ETA*B.T
    # deps_i/dt -= dH/dtau_i = sq2*g*(B^T xa)_i
    A[N:2*N,2*N:2*N+E]-= sq2g*B.T
    # dxa_e/dt = dH/dpa_e = Delta*pa_e + sq2*g*eta*(B eps)_e ; damping -kappa/2 xa
    A[2*N:2*N+E,2*N+E:2*N+2*E]+=Delta*np.eye(E)
    A[2*N:2*N+E,N:2*N]+= sq2g*ETA*B
    A[2*N:2*N+E,2*N:2*N+E]+= -kappa/2*np.eye(E)
    # dpa_e/dt = -dH/dxa_e = -Delta*xa_e - sq2*g*(B tau)_e ; damping -kappa/2 pa
    A[2*N+E:2*N+2*E,2*N:2*N+E]+= -Delta*np.eye(E)
    A[2*N+E:2*N+2*E,0:N]+= -sq2g*B
    A[2*N+E:2*N+2*E,2*N+E:2*N+2*E]+= -kappa/2*np.eye(E)
    return A

def eff_drift(g,Delta,kappa):
    gam=g*g*kappa/(Delta**2+kappa**2/4)
    Om=-g*g*Delta/(Delta**2+kappa**2/4)
    Gam=gam*ETA
    A=np.zeros((2*N,2*N))
    A[0:N,0:N]=-Gam*LG; A[0:N,N:2*N]=2*Om*ETA**2*LG
    A[N:2*N,N:2*N]=-Gam*LG; A[N:2*N,0:N]=-2*Om*LG
    return A,gam,Om,Gam

def ness(A,b):
    return np.linalg.lstsq(A,-b,rcond=None)[0]

h=np.zeros(N); h[3]=1.0; h[4]=1.0; hs=LAM*(h-h.mean())

print("== D: full vs effective (transient + NESS), sweep g ==")
kappa,Delta=2.0,1.0
for g in (0.05,0.1,0.2,0.4):
    Af=full_drift(g,Delta,kappa); Ae,gam,Om,Gam=eff_drift(g,Delta,kappa)
    bf=np.zeros(2*N+2*E); bf[0:N]=hs
    be=np.zeros(2*N); be[0:N]=hs
    Rf=ness(Af,bf); Re=ness(Ae,be)
    ef=Rf[N:2*N]-Rf[N:2*N].mean(); ee=Re[N:2*N]-Re[N:2*N].mean()
    pred=-(2*LAM*Om)/(Gam**2+4*Om**2*ETA**2)*(np.linalg.pinv(LG)@(h-h.mean()))
    dev_fe=np.abs(ef-ee).max()/(np.abs(ee).max()+1e-30)
    dev_ep=np.abs(ee-pred).max()/(np.abs(pred).max()+1e-30)
    print(f"  g={g}: eps_ss full-vs-eff dev={dev_fe:.2%}  eff-vs-closedform dev={dev_ep:.2e}  (Gam={Gam:.4f} Om={Om:+.4f})")

print()
print("== F/G: closure + force from the SAME H_int (g=0.1) ==")
g=0.1
Af=full_drift(g,Delta,kappa)
def eps_field(source_profile):
    b=np.zeros(2*N+2*E); b[0:N]=LAM*(source_profile-source_profile.mean())
    R=ness(Af,b); return R[N:2*N]-R[N:2*N].mean()
def gauss_prof(c,wd=1.0):
    x=np.arange(N); d=np.minimum(np.abs(x-c),N-np.abs(x-c))
    p=np.exp(-0.5*(d/wd)**2); return p/p.sum()   # unit total energy
# object B fixed at site 3; object A at site xA: U(xA) = lam * h_A(xA) . eps_B
hB=gauss_prof(3); epsB=eps_field(hB)
_,gam,Om,Gam=eff_drift(g,Delta,kappa)
C=(2*LAM*Om)/(Gam**2+4*Om**2*ETA**2)
U_meas=[]; U_pred=[]
for xA in range(N):
    hA=gauss_prof(xA)
    U_meas.append(LAM*float(hA@epsB))
    U_pred.append(-LAM*C*float(hA@np.linalg.pinv(LG)@(hB-hB.mean())))
U_meas=np.array(U_meas); U_pred=np.array(U_pred)
dev=np.abs((U_meas-U_meas.mean())-(U_pred-U_pred.mean())).max()/np.abs(U_pred-U_pred.mean()).max()
print(f"  U_AB(x) measured (full model) vs predicted -lam^2 C h_A L^+ h_B: dev = {dev:.2e}")
F_meas=-(U_meas[6]-U_meas[4])/2; F_pred=-(U_pred[6]-U_pred[4])/2
print(f"  force at d=2 sites from B: measured {F_meas:+.4e}  predicted {F_pred:+.4e}")
sign="ATTRACTIVE" if U_meas[4]<U_meas[8] else "REPULSIVE"
print(f"  sign at Delta=+1: potential near B is {'lower' if U_meas[4]<U_meas[8] else 'higher'} -> {sign}")

print()
print("== N3: F(A<-B) vs F(B<-A), full model ==")
hA5=gauss_prof(7)
epsA=eps_field(hA5)
FaB=-(LAM*float(gauss_prof(8)@epsB)-LAM*float(gauss_prof(6)@epsB))/2
FbA=-(LAM*float(gauss_prof(4)@epsA)-LAM*float(gauss_prof(2)@epsA))/2
print(f"  F_A<-B = {FaB:+.5e}   F_B<-A = {FbA:+.5e}   Delta_N3 = {abs(FaB+FbA)/(0.5*(abs(FaB)+abs(FbA))):.2e}")

print()
print("== H: controls ==")
for Dl in (0.0,-1.0):
    Af2=full_drift(g,Dl,kappa)
    b=np.zeros(2*N+2*E); b[0:N]=LAM*(hB-hB.mean())
    R2=ness(Af2,b); e2=R2[N:2*N]-R2[N:2*N].mean()
    _,gam2,Om2,Gam2=eff_drift(g,Dl,kappa)
    print(f"  Delta={Dl:+.1f}: max|eps_ss| = {np.abs(e2).max():.3e}   (Om={Om2:+.4f}, Gam={Gam2:.4f})")
print("  scaling sweep (g,Delta,kappa) vs formula:")
devs=[]
rng=[(0.05,0.5,1.5),(0.1,2.0,3.0),(0.2,-0.8,2.5),(0.08,1.5,1.0),(0.15,-2.0,4.0)]
for (gg,DD,kk) in rng:
    Afx=full_drift(gg,DD,kk); b=np.zeros(2*N+2*E); b[0:N]=LAM*(hB-hB.mean())
    ex=ness(Afx,b)[N:2*N]; ex-=ex.mean()
    _,gamx,Omx,Gamx=eff_drift(gg,DD,kk)
    px=-(2*LAM*Omx)/(Gamx**2+4*Omx**2*ETA**2)*(np.linalg.pinv(LG)@(hB-hB.mean()))
    devs.append(np.abs(ex-px).max()/np.abs(px).max())
    print(f"    g={gg} Delta={DD} kappa={kk}: dev vs formula = {devs[-1]:.2%}")
print(f"  max scaling dev = {max(devs):.2%}  (gate < 5%)")
print("done.")