"""Independent decoder tests. Python 3 and NumPy. No quantum binary encoder is assumed to run here.""" from fractions import Fraction as F import itertools import json from pathlib import Path import math import numpy as np def rational_unit(z): return (1-z*z)/(1+z*z), 2*z/(1+z*z) def exact_fixture(s=F(1,256)): # A mixture of three orthonormal bases; each outcome has effect (1/3)P_y. U=2*s; r=s/4 a,b = rational_unit(F(1,65536)) y=[] for z in [F(0),F(1,64),F(1,2)]: c,e=rational_unit(z) y.extend([(c,e,0,0),(-e,c,0,0),(0,0,1,0),(0,0,0,1)]) w=F(1,3) t0=[F(c[0])**2 for c in y] born=[w*(a*c[0]+b*c[1])**2 for c in y] assert sum(born)==1 near=[i for i,v in enumerate(t0) if v<=s] far=[i for i in range(len(y)) if i not in near] big=sum(w*t0[i] for i in far) u=s*a*a+b*b v=[1+(1/s-1)*h for h in t0] bt=[(a*c[0]+b*c[1])**2/(u*v[i]) for i,c in enumerate(y)] # All-three-vertex local barycenter is the true state as a matrix. weights=[1-b*b/r-a*b/r, b*b/r, a*b/r] assert all(f>0 for f in weights) and sum(weights)==1 assert 0 < u/U < 1 cuts=sorted({F(0),F(1),*(bt[i] for i in near)}) mean=[F(0)]*len(y); trace=[]; broken=0; neg_vertices=0 for k,vw in enumerate(weights): f=[] for c,e,*_ in y: val=F(c)**2 if k==1: val+=r*(e*e-c*c) if k==2: val+=2*r*c*e f.append(w*val) neg_vertices+=int(min(f)<0) assert sum(f)==1 for flag,pflag in [(0,1-u/U),(1,u/U)]: for l,h in zip(cuts[:-1],cuts[1:]): mid=(l+h)/2 q=[F(0)]*len(y) for i in near: q[i]=w*U*v[i]*flag*(mid=0 and correct_total==1 assert faulty_total-correct_total==near_mass broken+=int(faulty_total!=1) prob=vw*pflag*(h-l) assert prob>0 for i,val in enumerate(q): mean[i]+=prob*val trace.append(dict(vertex=k,gate=flag,interval=[str(l),str(h)], probability=str(prob),near_mass=str(near_mass), correct_total=str(correct_total),faulty_total=str(faulty_total))) assert sum(F(row["probability"]) for row in trace)==1 assert mean==born and len(trace)==18 and broken==6 # An interior formal vertex with negative eigenvalue can occur even if # this fixture's particular effect list doesn't reveal that eigenvalue. assert (1+r)*0-r*r<0 return dict(scale=str(s),outcomes=len(y),enumerated_blocks=len(trace),exact_residual="0", mutant_mass_failures=broken,vertices_with_negative_effect_scores=neg_vertices, mean=[str(value) for value in mean],born=[str(value) for value in born],blocks=trace) def unitary(rng,n): z=rng.normal(size=(n,n))+1j*rng.normal(size=(n,n)) q,r=np.linalg.qr(z) return q*(np.diag(r)/abs(np.diag(r)))[None,:] def frame(rng,d,h): # First row has prescribed weights. Remaining rows are Haar-scrambled. a=np.sqrt(h); a[0]-=1 H=np.eye(d)-2*np.outer(a,a)/(a@a) if a@a>0 else np.eye(d) H=H.astype(complex) H[1:]=unitary(rng,d-1)@H[1:] return H.T def stress(): rng=np.random.default_rng(625170) result=dict(frames=0,fixed_tape_vectors=0,largest_mass_error=0.0, largest_mean_error=0.0,minimum_component=1.0, min_far_ratio=1.0,near_cases=0,far_cases=0, negative_raw_score_cases=0,mutant_max_mass_error=0.0) for d in [4,5,7,8]: s=1/(32*d); U=2*s; r=s/4; delta=r/math.sqrt(2*(d*d-1)) p0=np.zeros((d,d),complex);p0[0,0]=1 for j in range(180): # Includes exact zeros, both sides of the near cutoff, repeated # outcomes, union of rotated bases, and general nonprojective POVMs. h=rng.dirichlet(np.ones(d)); h[1]=s*[0,0.5,1-1e-8,1,1+1e-8,2][j%6] others=np.arange(d)!=1; h[others]*=(1-h[1])/sum(h[others]) Y=frame(rng,d,h) if j%4==2: V=unitary(rng,2*d)[:,:d].conj() weights=np.sum(abs(V)**2,axis=1) Y=V/np.sqrt(weights[:,None]) else: Y=np.concatenate([Y,np.eye(d)],axis=0); weights=np.full(2*d,0.5) assert np.max(abs((Y.T*weights)@Y.conj()-np.eye(d)))<1e-12 t0=abs(Y[:,0])**2; near=t0<=s; far=~near z=unitary(rng,d-1)[:,0]; rad=delta*rng.random() if j%3 else 0. x=np.r_[math.sqrt(1-rad**2),rad*z] born=weights*abs(Y.conj()@x)**2 ax=x.copy(); ax[0]*=math.sqrt(s); u=np.vdot(ax,ax).real iy=Y.copy();iy[:,0]/=math.sqrt(s); v=np.sum(abs(iy)**2,axis=1) transformed=abs((iy/np.sqrt(v[:,None])).conj()@(ax/np.sqrt(u)))**2 R=weights[far]*t0[far];R/=R.sum() def scores(Q): return weights*np.einsum('ij,jk,ik->i',Y.conj(),Q,Y).real def complete(f,g): q=np.empty(len(Y));q[near]=g q[far]=f[far]+R*(sum(f[near])-sum(g)) return q # Expected formula from the explicitly transformed Born probabilities. P=np.outer(x,x.conj()) g=weights[near]*U*v[near]*(u/U)*transformed[near] err=np.max(abs(complete(scores(P),g)-born)) result['largest_mean_error']=max(result['largest_mean_error'],float(err)) # Vertices are extreme perturbations of the reference. Random tests # plus a targeted negative shift in a near direction. H=unitary(rng,d) for k in [0,1,2]: if k==0: Q=p0 elif k==1: Q=p0+r*(np.outer(H[:,0],H[:,0].conj())-np.outer(H[:,1],H[:,1].conj())) else: Q=p0+r*(p0-np.diag([0,1]+[0]*(d-2))) f=scores(Q) result['negative_raw_score_cases']+=int(np.min(f)<-1e-12) for response in (np.zeros(sum(near)),np.ones(sum(near)),rng.random(sum(near))): g=weights[near]*U*v[near]*response q=complete(f,g); masserr=abs(q.sum()-1) result['minimum_component']=min(result['minimum_component'],float(q.min())) result['largest_mass_error']=max(result['largest_mass_error'],float(masserr)) result['min_far_ratio']=min(result['min_far_ratio'],float(np.min(q[far]/(weights[far]*t0[far])))) faulty=q.copy() faulty[far]=f[far]+R*sum(f[near]) faulty_error=abs(float(faulty.sum())-1.0) result['mutant_max_mass_error']=max(result['mutant_max_mass_error'],faulty_error) result['fixed_tape_vectors']+=1 result['frames']+=1; result['near_cases']+=int(sum(near));result['far_cases']+=int(sum(far)) assert result['minimum_component']>=-1e-13 and result['largest_mass_error']<2e-12 assert result['largest_mean_error']<2e-12 and result['negative_raw_score_cases']>0 assert result['mutant_max_mass_error']>0.01 return result if __name__=='__main__': if not __debug__: raise SystemExit("These checks require assertions. Run Python without -O.") print(json.dumps(dict(exact=exact_fixture(),exact_theorem_scale=exact_fixture(F(1,128)),stress=stress()),indent=2))