Simulate circuits containing PauliMeasurements using PauliFrames

src/pauli_frames.jl

@@ -11,12 +11,14 @@ struct PauliFrame{T,S} <: AbstractQCState
     measurements::S # TODO check if when looping over this we are actually looping over the fast axis
 end
 
-nqubits(f::PauliFrame) = nqubits(f.frame)
+nqubits(f::PauliFrame) = nqubits(f.frame) - 1 # dont count the ancilla qubit
 Base.length(f::PauliFrame) = size(f.measurements, 1)
 Base.eachindex(f::PauliFrame) = 1:length(f)
 Base.copy(f::PauliFrame) = PauliFrame(copy(f.frame), copy(f.measurements))
 Base.view(frame::PauliFrame, r) = PauliFrame(view(frame.frame, r), view(frame.measurements, r, :))
 
+tab(f::PauliFrame) = Tableau(f.frame.tab.phases, nqubits(f), f.frame.tab.xzs)
+
 fastrow(s::PauliFrame) = PauliFrame(fastrow(s.frame), s.measurements)
 fastcolumn(s::PauliFrame) = PauliFrame(fastcolumn(s.frame), s.measurements)
 
@@ -26,7 +28,8 @@ $(TYPEDSIGNATURES)
 Prepare an empty set of Pauli frames with the given number of `frames` and `qubits`. Preallocates spaces for `measurement` number of measurements.
 """
 function PauliFrame(frames, qubits, measurements)
-    stab = fastcolumn(zero(Stabilizer, frames, qubits)) # TODO this should really be a Tableau
+    # one extra qubit for ancilla measurements
+    stab = fastcolumn(zero(Stabilizer, frames, qubits + 1)) # TODO this should really be a Tableau
     bits = zeros(Bool, frames, measurements)
     frame = PauliFrame(stab, bits)
     initZ!(frame)
@@ -120,6 +123,24 @@ function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRY, and faster sMRX
     return frame
 end
 
+function apply!(frame::PauliFrame, op::PauliMeasurement)
+    # this is inspired by ECC.naive_syndrome_circuit
+    n = nqubits(op.pauli)
+    for qubit in 1:n
+        if op.pauli[qubit] == (1, 0)
+            apply!(frame, sXCZ(qubit, n + 1))
+        elseif op.pauli[qubit] == (0, 1)
+            apply!(frame, sCNOT(qubit, n + 1))
+        elseif op.pauli[qubit] == (1, 1)
+            apply!(frame, sYCX(qubit, n + 1))
+        end
+    end
+    op.pauli.phase[] == 0 || apply!(frame, sX(n + 1))
+    apply!(frame, sMRZ(n + 1, op.bit))
+
+    return frame
+end
+
 function applynoise!(frame::PauliFrame,noise::UnbiasedUncorrelatedNoise,i::Int)
     p = noise.p
     T = eltype(frame.frame.tab.xzs)

src/sumtypes.jl

@@ -223,6 +223,12 @@ function concrete_typeparams(t::Type{NoiseOp})
     ]
 end
 
+function concrete_typeparams(::Type{PauliMeasurement})
+    return [
+        (Array{UInt8,0}, Vector{UInt64}),
+    ]
+end
+
 
 # XXX This has to happen after defining all the `concrete_typeparams` methods
 

test/test_ecc_syndromes.jl

@@ -32,8 +32,8 @@
             p = random_pauli(dataqubits, realphase=true)
             pₙ = embed(naive_qubits, 1:dataqubits, p)
             pₛ = embed(shor_qubits, 1:dataqubits, p)
-            mul_left!(naive_frames.frame, pₙ)
-            mul_left!(shor_frames.frame, pₛ)
+            mul_left!(tab(naive_frames), pₙ)
+            mul_left!(tab(shor_frames), pₛ)
             # run the syndrome circuits using the public API
             pftrajectories(naive_frames, naive_scirc)
             pftrajectories(shor_frames, shor_scirc)

test/test_pauliframe.jl

@@ -88,4 +88,29 @@
             @test all(0.25.*[1 0 0 0 1] .<= (sum(m, dims=1)[:,1:5])./n .<= 0.75.*[1 0 0 0 1])
         end
     end
+
+    @testset "PauliMeasurements" begin
+        n = 2000
+        state = Register(one(MixedDestabilizer, 3), 5)
+        frame = PauliFrame(n, 3, 5)
+
+        glassy_ghz_circuit = [
+            sHadamard(1), sHadamard(2), sHadamard(3),
+            PauliMeasurement(P"ZZ_", 1), PauliMeasurement(P"_ZZ", 2),
+            sMZ(1, 3), sMZ(2, 4), sMZ(3, 5)
+        ]
+        for m in [pfmeasurements(pftrajectories(copy(frame), glassy_ghz_circuit)),
+            pfmeasurements(pftrajectories(glassy_ghz_circuit; trajectories=n, threads=false)),
+            pfmeasurements(pftrajectories(glassy_ghz_circuit; trajectories=n, threads=true))]
+
+            # decode based on measurement outcomes
+            for r in eachrow(m)
+                r[4] ⊻= r[1]
+                r[5] ⊻= r[1] ⊻ r[2]
+            end
+
+            # check that the correct correlations are present
+            @test all(m[:, 3] .== m[:, 4] .== m[:, 5])
+        end
+    end
 end