Unexecuted reproducibility supplement

def psd_metrics(A): A = np.asarray(A) need(np.all(np.isfinite(A)),"nonfinite PSD matrix") vals = np.linalg.eigvalsh(herm(A)) return "hermiticity_op":float(np.linalg.norm(A-A.conj().T,2)), "minimum_eigenvalue":float(vals.min()), "negative_part_fro":float(np.linalg.norm(np.minimum(vals,0)))

def accept(name,psds=None,equalities=None,tol=STOL): psds,equalities = psds or {},equalities or {} record = {"PSD":k:psd_metrics(v) for k,v in psds.items(), "affine_fro":k:float(np.linalg.norm(np.asarray(v))) for k,v in equalities.items(), "tolerance":tol,"certified":False} DATA["residuals"][name] = record dump() for k,m in record["PSD"].items(): need(m["hermiticity_op"]<=tol and m["negative_part_fro"]<=tol, name+" PSD residual "+k) for k,v in record["affine_fro"].items(): need(np.isfinite(v) and v<=tol,name+" affine residual "+k)

def kron(*ops): A = np.array([[1]],complex) for B in ops: A = np.kron(A,B) return A

def bits(i,n): return tuple((i>>(n-1-j))&1 for j in range(n))

def integer(b): n = 0 for z in b: n = 2*n+int(z) return n

def embed(G,pos,n): A = np.zeros((2**n,2**n),complex) for col in range(2**n): b = list(bits(col,n)) gc = integer([b[p] for p in pos]) for row in range(2**len(pos)): a = b.copy() for p,z in zip(pos,bits(row,len(pos))): a[p] = z A[integer(a),col] += G[row,gc] return A

def units(d=2): for i,j in it.product(range(d),repeat=2): E = np.zeros((d,d),complex) E[i,j] = 1 yield i,j,E

def dep(A): return np.trace(A)*I/2

def ptout(A,o,b): return sum(A[a*b:(a+1)*b,a*b:(a+1)*b] for a in range(o))

def ptrace(A,keep,n): keep = tuple(keep) drop = tuple(i for i in range(n) if i not in keep) axes = keep+drop+tuple(n+i for i in keep)+tuple(n+i for i in drop) A = A.reshape((2,)*(2*n)).transpose(axes) A = A.reshape(2**len(keep),2**len(drop), 2**len(keep),2**len(drop)) return np.einsum("abcb->ac",A)

def choi(N,d=2): o = len(N[0,0]) return sum((np.kron(N[i,j],E) for i,j,E in units(d)), np.zeros((o*d,o*d),complex))

def channel_check(N,name,tol=TOL): d = int(round(np.sqrt(len(N)))) o = len(N[0,0]) J = choi(N,d) save(name+"_Choi",J=J) jwrite(name+"_Choi_basis","ordering":"output then input", "output_dimension":o,"input_dimension":d,"normalized":False) accept(name,"Choi":J,"TP":ptout(J,o,d)-np.eye(d),tol) return J

def entropy(A): accept("entropy_"+str(len(DATA["residuals"])), "state":A,"trace":np.array(np.trace(A)-1),TOL) vals = np.linalg.eigvalsh(herm(A)) vals = vals[vals>1e-13] return float(-np.sum(vals*np.log2(vals)))

def occupation_basis(J,N): return [n for n in it.product(range(N+1),repeat=J) if sum(n)<=N]

def permanent(A): n = len(A) if n==0: return 1.+0j return sum(np.prod([A[i,p[i]] for i in range(n)]) for p in it.permutations(range(n)))