Simulate circuits containing PauliMeasurements using PauliFrames (#412)

Co-authored-by: Stefan Krastanov <[email protected]>
Co-authored-by: Stefan Krastanov <[email protected]>

Author: J-C-Q

CHANGELOG.md

@@ -5,6 +5,10 @@
 
 # News
 
+## v0.9.15 - 2024-12-22
+
+- `pftrajectories` now supports fast multiqubit measurements with `PauliMeasurement` in addition to the already supported single qubit measurements `sMX/Z/Y` and workarounds like `naive_syndrome_circuit`.
+
 ## v0.9.14 - 2024-11-03
 
 - **(fix)** `affectedqubits()` on `sMX`, `sMY`, and `sMR*`
@@ -75,7 +79,7 @@
 - Gate errors are now conveniently supported by the various ECC benchmark setups in the `ECC` module.
 - Significant improvements to the low-level circuit compiler (the sumtype compactifier), leading to faster Pauli frame simulation of noisy circuits.
 - Bump `QuantumOpticsBase.jl` package extension compat bound.
-- **(fix)** Remove printing of spurious debug info from the PyBP decoder. 
+- **(fix)** Remove printing of spurious debug info from the PyBP decoder.
 - **(fix)** Failed compactification of gates now only raises a warning instead of throwing an error. Defaults to slower non-compactified gates.
 
 ## v0.9.3 - 2024-04-10
@@ -93,7 +97,7 @@
 - Implemented `iscss` function to identify whether a given code is known to be a CSS (Calderbank-Shor-Steane) code.
 - Added the classical Reed-Muller code in the ECC module.
 - Added the surface code to the ECC module.
- 
+
 ## v0.9.0 - 2024-03-19
 
 - **(breaking)** The defaults in `random_pauli` are now `realphase=true` and `nophase=true`.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <[email protected]> and QuantumSavory community members"]
-version = "0.9.14"
+version = "0.9.15"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

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)
@@ -112,6 +115,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
     xzs = tab(frame.frame).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_bitpack.jl

@@ -1,6 +1,6 @@
 @testitem "Alternative bit packing" tags=[:bitpack] begin
     using Random
-    using QuantumClifford: Tableau
+    using QuantumClifford: Tableau, mul_left!
 
     @testset "alternative bit packing" begin
         for n in [1,3] # can not go higher than 4 (limitation from SIMD acting on transposed/strided arrays)
@@ -20,6 +20,8 @@
             after_cliff = stab_to_gf2(_after_clif);
             after_cliff_phases = phases(_after_clif);
 
+            after_mul = stab_to_gf2(mul_left!(copy(s64), p64))
+
             for int in [UInt8, UInt16, UInt32, UInt64]
                 p = PauliOperator(p64.phase, N, collect(reinterpret(int,p64.xz)));
                 xzs = collect(reinterpret(int, collect(xzs64)));
@@ -45,8 +47,38 @@
                         @test after_cliff == stab_to_gf2(after_clifford)
                         @test after_cliff_phases == phases(after_clifford)
                     end
+
+                    @test after_mul == stab_to_gf2(mul_left!(copy(s), p))
                 end
             end
         end
     end
+
+    @testset "fast column and fast row mul_left correctness" begin
+        reinterpret_stab(s) = Stabilizer(Tableau(copy(phases(s)), nqubits(s), collect(reinterpret(UInt8, collect(s.tab.xzs)))))
+        reinterpret_p(p) = PauliOperator(p.phase, nqubits(p), collect(reinterpret(UInt8, p.xz)))
+        for N in [7,8,9,33,61,62,63,64,65,66]
+
+            s0 = random_stabilizer(2,N)
+            p = random_pauli(N)
+            s = copy(s0)
+            sr = QuantumClifford.fastrow(copy(s))
+            sc = QuantumClifford.fastcolumn(copy(s))
+
+            mul_left!(s, p)
+            mul_left!(sr, p)
+            mul_left!(sc, p)
+
+            s8 = reinterpret_stab(copy(s0))
+            s8r = QuantumClifford.fastrow(copy(s8))
+            s8c = QuantumClifford.fastcolumn(copy(s8))
+            p8 = reinterpret_p(copy(p))
+            mul_left!(s8, p8)
+            mul_left!(s8r, p8)
+            mul_left!(s8c, p8)
+
+            @test stab_to_gf2(s) == stab_to_gf2(sr) == stab_to_gf2(sc) == stab_to_gf2(s8) == stab_to_gf2(s8r) == stab_to_gf2(s8c)
+        end
+    end
+    end
 end

test/test_ecc_base.jl

@@ -160,11 +160,14 @@ const code_instance_args = Dict(
 
 function all_testablable_code_instances(;maxn=nothing)
     codeinstances = []
+    i = 1
     for t in subtypes(QuantumClifford.ECC.AbstractECC)
         for c in get(code_instance_args, t.name.name, [])
             codeinstance = t(c...)
             !isnothing(maxn) && nqubits(codeinstance) > maxn && continue
             push!(codeinstances, codeinstance)
+            #@show i, t, code_n(codeinstance), code_k(codeinstance), code_s(codeinstance), code_n(codeinstance)-code_k(codeinstance)
+            i += 1
         end
     end
     return codeinstances

test/test_ecc_syndromes.jl

@@ -5,6 +5,11 @@
 
     include("test_ecc_base.jl")
 
+    using QuantumClifford: Tableau
+    reinterpret_frame(frame) = PauliFrame(reinterpret_stab(frame.frame), copy(frame.measurements))
+    reinterpret_stab(s) = Stabilizer(Tableau(copy(phases(s)), nqubits(s), collect(reinterpret(UInt8, collect(s.tab.xzs)))[[1:1+(nqubits(s)-1)÷8;end÷2+1:end÷2+1+(nqubits(s)-1)÷8],:]))
+    reinterpret_p(p) = PauliOperator(p.phase, nqubits(p), collect(reinterpret(UInt8, p.xz))[[1:1+(nqubits(p)-1)÷8;end÷2+1:end÷2+1+(nqubits(p)-1)÷8]])
+
     function pframe_naive_vs_shor_syndrome(code)
         ecirc = naive_encoding_circuit(code)
         naive_scirc, naive_ancillaries = naive_syndrome_circuit(code)
@@ -30,19 +35,35 @@
             pftrajectories(shor_frames, vcat(ecirc, shor_cat_scirc))
             # manually injecting the same type of noise in the frames -- not really a user accessible API
             p = random_pauli(dataqubits, realphase=true)
-            pₙ = embed(naive_qubits, 1:dataqubits, p)
-            pₛ = embed(shor_qubits, 1:dataqubits, p)
+            pₙ = embed(naive_qubits+1, 1:dataqubits, p) # +1 to account for the buffer qubit hidden in pauli frames
+            pₛ = embed(shor_qubits+1, 1:dataqubits, p)  # +1 to account for the buffer qubit hidden in pauli frames
             mul_left!(naive_frames.frame, pₙ)
             mul_left!(shor_frames.frame, pₛ)
             # run the syndrome circuits using the public API
             pftrajectories(naive_frames, naive_scirc)
             pftrajectories(shor_frames, shor_scirc)
             @test pfmeasurements(naive_frames) == pfmeasurements(shor_frames)[:,shor_bits]
+
+            # just for completeness, let's also try bitpacking in UInt8 instead of the default UInt
+            _naive_frames = PauliFrame(nframes, naive_qubits, syndromebits)
+            _shor_frames = PauliFrame(nframes, shor_qubits, last(shor_bits))
+            naive_uint8 = reinterpret_frame(_naive_frames)
+            shor_uint8 = reinterpret_frame(_shor_frames)
+            pftrajectories(naive_uint8, ecirc)
+            pftrajectories(shor_uint8, vcat(ecirc, shor_cat_scirc))
+            p_uint8 = reinterpret_p(p)
+            pₙ_uint8 = embed(naive_qubits+1, 1:dataqubits, p_uint8)
+            pₛ_uint8 = embed(shor_qubits+1, 1:dataqubits, p_uint8)
+            mul_left!(naive_uint8.frame, pₙ_uint8)
+            mul_left!(shor_uint8.frame, pₛ_uint8)
+            pftrajectories(naive_uint8, naive_scirc)
+            pftrajectories(shor_uint8, shor_scirc)
+            @test pfmeasurements(shor_uint8)[:,shor_bits] == pfmeasurements(shor_frames)[:,shor_bits] == pfmeasurements(naive_frames) == pfmeasurements(naive_uint8)
         end
     end
 
     @testset "naive and shor measurement circuits" begin
-        for c in all_testablable_code_instances()
+        for (i,c) in enumerate(all_testablable_code_instances())
             pframe_naive_vs_shor_syndrome(c)
         end
     end

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