Computational supplement

def loss_domain_audit(): K,c = 4,3 sets = [()] for k in range(1,c+1): sets += [tuple(sorted(S,reverse=True)) for S in it.combinations(range(1,K+1),k) if max(S)-min(S)+1<=c] cfg = [(S,b) for S in sets for b in it.product(range(2),repeat=len(S))] M = 1+K+sum((K-d)*2**(d-1) for d in range(1,c)) D = 1+2*K+sum((K-d)*4*3**(d-1) for d in range(1,c)) need(len(sets)==M and len(cfg)==D,"loss-domain count") need((3,1) in sets,"missing nonconsecutive input") ix = z:i for i,z in enumerate(cfg) omega = 1:1.,2:2.,3:3.,4:5. gamma = r:.2+.12*r for r in range(1,K+1) A = np.zeros((D*D,D),complex) for col,(S,b) in enumerate(cfg): for mask in it.product(range(2),repeat=len(S)): ds = tuple(r for r,z in zip(S,mask) if z) fs = tuple(r for r,z in zip(S,mask) if not z) db = tuple(a for a,z in zip(b,mask) if z) fb = tuple(a for a,z in zip(b,mask) if not z) amp = np.prod([np.sqrt(gamma[r] if z else 1-gamma[r]) for r,z in zip(S,mask)]) A[ix[ds,db]*D+ix[fs,fb],col] += amp e = np.array([sum(omega[r] for r in S) for S,b in cfg]) eo = (e[:,None]+e[None,:]).reshape(-1) check("full_loss_domain_isometry",A.conj().T@A,np.eye(D)) check("nonuniform_loss_energy",eo[:,None]*A,A*e[None,:]) save("loss_domain",scattering=A,energy=e)

def put(out,key,v): # Exact-zero optimization only: no nonzero amplitude is truncated. if np.any(v!=0): out[key] = out.get(key,np.zeros(2,complex))+v

def localized(K,ell=0): J = np.zeros((2**K,2),complex) for m in range(2): J[m,m] = 1 U = np.eye(2**K,dtype=complex) for edge in range(K-2,K-ell-2,-1): U = embed(SW,(edge,edge+1),K)@U return U@J

def initial(J,K): return ((K,bits(i,K)),(),()):J[i] for i in range(2**K) if np.any(J[i]!=0)

def advance(entries,j,law,U,gamma): out = {} for ((n,b),dh,fh),v in entries.items(): G = U.get((j,n),np.eye(2**n,dtype=complex)) for a in range(2**n): z = G[a,integer(b)] if n else 1 if z==0: continue bb = bits(a,n) for k,p in enumerate(law[n]): if p==0: continue modes = tuple(range(n,n-k,-1)) for mask in it.product(range(2),repeat=k): dp = tuple((r,x) for r,x,s in zip(modes,bb[:k],mask) if s) fp = tuple((r,x) for r,x,s in zip(modes,bb[:k],mask) if not s) amp = z*np.sqrt(float(p))*np.prod([ np.sqrt(gamma(j,r) if s else 1-gamma(j,r)) for r,s in zip(modes,mask)]) put(out,((n-k,bb[k:]),dh+(dp,),fh+(fp,)),amp*v) return out

def dense(entries): dk = sorted(d for core,d,f in entries) ek = sorted((core,f) for core,d,f in entries) di,ei = d:i for i,d in enumerate(dk),e:i for i,e in enumerate(ek) A = np.zeros((len(dk),len(ek),2),complex) for (core,d,f),v in entries.items(): A[di[d],ei[core,f]] += v return A,dk,ek