Unexecuted reproducibility supplement

def annihilators(basis): lookup = n:i for i,n in enumerate(basis) out = [] d = len(basis)+1 for a in range(len(basis[0])): B = np.zeros((d,d),complex) for col,n in enumerate(basis): if n[a]: m = list(n) m[a] -= 1 B[lookup[tuple(m)],col] = np.sqrt(n[a]) out.append(B) return out

def gaussian_cutoff(K,N,name): J = len(K) accept(name+"_covariance","K":K,tol=TOL) vals,W = np.linalg.eigh(herm(K)) # No silent eigenvalue clipping in state construction. need(vals.min()>=0,"negative covariance eigenvalue") z = vals/(1+vals) basis = occupation_basis(J,N) T = [email protected](np.eye(J)+K) p0 = 1/float(np.linalg.det(np.eye(J)+K).real) rho = np.zeros((len(basis)+1,len(basis)+1),complex) for i,a in enumerate(basis): rows = [r for r,n in enumerate(a) for _ in range(n)] for j,b in enumerate(basis): if sum(a)!=sum(b): continue cols = [r for r,n in enumerate(b) for _ in range(n)] den = np.sqrt(np.prod([math.factorial(n) for n in a+b])) rho[i,j] = p0*permanent(T[np.ix_(rows,cols)])/den tail = 1-float(np.trace(rho).real) rho[-1,-1] = tail save(name,rho=rho,K=K,C=K.T,T=T,eigenvalues=vals) width = int(math.ceil(math.log2(len(rho)))) jwrite(name+"_basis", "matrix_axes":["ket occupation or flag","bra occupation or flag"], "occupation_order":basis,"flag_index":len(basis), "raw_width":width, "physical_to_binary":[format(i,"0 "unused_binary_indices":list(range(len(rho),2**width)), "unused_prepared_weight":0, "reverse_map":"pinch valid/invalid; invalid weight goes to flag", "annihilators_on_flag":0) accept(name+"_density","rho":rho, "trace":np.array(np.trace(rho)-1),TOL)

probs = np.zeros(N+1) for i,n in enumerate(basis): probs[sum(n)] += rho[i,i].real independent_probs = np.zeros(N+1) means = np.zeros(J) for n in basis: weight = float(np.prod((1-z)*z**np.array(n))) independent_probs[sum(n)] += weight means += np.array(n)*weight check(name+"_occupation_distribution",probs,independent_probs) B = annihilators(basis) Ctrunc = np.array([[np.trace([email protected]().T@Bop) for Bop in B] for A in B]) expectedC = ([email protected](means)@W.conj().T).T check(name+"_phase_sensitive_moments",Ctrunc,expectedC) save(name+"_moments",C_truncated=Ctrunc,C_expected=expectedC, occupation_probabilities=probs,normal_mode_retained_means=means) return rho,basis,tail,probs,Ctrunc