import json,time
import numpy as np
L,W=4,2

def make(seed,K):
 r=np.random.default_rng(seed); E=[]
 for t in range(L):
  src=range(t*W,(t+1)*W) if t else [-1]; dst=range((t+1)*W,(t+2)*W) if t<L-1 else [-2]
  E += [(a,b) for a in src for b in dst]
 F=[]
 for _ in E:
  lo=np.maximum(.03,r.normal(1.2,.25,K)); hi=lo+r.uniform(.05,.75,K); m=r.dirichlet(np.ones(K)); F.append((lo,hi,m))
 return E,F

def allroutes(E):
 ix={p:i for i,p in enumerate(E)}; R=[]
 for q in np.ndindex(*(W,)*(L-1)):
  ns=[-1]+[t*W+q[t-1] for t in range(1,L)]+[-2]; x=np.zeros(len(E),dtype=np.int8)
  for a,b in zip(ns[:-1],ns[1:]): x[ix[a,b]]=1
  R.append(x)
 return np.array(R)

def base(F): return np.array([np.dot(m,lo) for lo,hi,m in F])
def blockcost(lo,m,g):
 # singleton is exact; merged block is conservatively raised to max lower.
 mass=m[g].sum(); raw=np.dot(m[g],lo[g])
 return float(raw if len(g)==1 else max(raw,mass*max(lo[g])))

def compress(F,X,target,mode):
 G=[[[i] for i in range(len(F[e][0]))] for e in range(len(F))]
 c=base(F); x0=X[np.argmin(X@c)]; total=sum(map(len,G[0]))*len(G)
 while sum(len(g) for z in G for g in z)>target*len(G):
  best=(1e99,None,None)
  for e,gs in enumerate(G):
   lo,hi,m=F[e]
   if len(gs)<2: continue
   old={id(g):blockcost(lo,m,g) for g in gs}
   for i in range(len(gs)):
    for j in range(i+1,len(gs)):
     a,b=gs[i],gs[j]; merged=a+b; d=blockcost(lo,m,merged)-old[id(a)]-old[id(b)]
     if mode=='regret': s=d*(1 if x0[e] else 0)+.05*d
     elif mode=='random': s=float((e*997+i*31+j)%1000)
     elif mode=='mass': s=-(m[a].sum()+m[b].sum())
     else:
      span=max(max(hi[a]),max(hi[b]))-min(min(lo[a]),min(lo[b])); overlap=max(0,min(max(hi[a]),max(hi[b]))-max(min(lo[a]),min(lo[b])))
      s=-overlap/max(span,1e-9)
     if s<best[0]: best=(s,e,(i,j,d))
  if best[1] is None: break
  _,e,(i,j,_)=best; G[e]=G[e][:i]+[G[e][i]+G[e][j]]+G[e][i+1:j]+G[e][j+1:]
 return np.array([sum(blockcost(lo,m,g) for g in G[e]) for e,(lo,hi,m) in enumerate(F)])

def run():
 result={'theorem':{},'sweep':{},'baselines':{}}
 maxviol=0.; ratios=[]
 for seed in range(10):
  E,F=make(seed,8); X=allroutes(E); c=base(F); ch=compress(F,X,3,'regret'); delta=X@(ch-c); s=np.argmin(X@c); h=np.argmin(X@ch); reg=(X[h]-X[s])@c
  maxviol=max(maxviol,float(reg-delta[s])); ratios.append(reg/(delta[s]+1e-12))
 result['theorem']={'max_bound_violation':float(maxviol),'mean_regret_over_bound':float(np.mean(ratios))}
 for K in (4,8,12):
  for target in (1,2,3):
   a=[]
   for seed in range(3):
    E,F=make(100+seed,K); X=allroutes(E); c=base(F); s=np.argmin(X@c); ch=compress(F,X,target,'regret'); h=np.argmin(X@ch); d=X@(ch-c); a.append(((X[h]-X[s])@c,h!=s,d[s],sum(ch-c)))
   result['sweep'][f'K{K}_to{target}']={'regret':float(np.mean([z[0] for z in a])),'flip':float(np.mean([z[1] for z in a])),'bound':float(np.mean([z[2] for z in a])),'inflation':float(np.mean([z[3] for z in a]))}
 for mode in ('regret','jaccard','random'):
  a=[]
  for seed in range(5):
   E,F=make(900+seed,8); X=allroutes(E); c=base(F); s=np.argmin(X@c); ch=compress(F,X,2,mode); h=np.argmin(X@ch); a.append(((X[h]-X[s])@c,h!=s))
  result['baselines'][mode]={'regret':float(np.mean([z[0] for z in a])),'flip':float(np.mean([z[1] for z in a]))}
 print(json.dumps(result,indent=2))
if __name__=='__main__': run()