Merge branch 'optimize-inner-gpu-type' of github.com:Krastanov/Quantu…

…mClifford.jl into optimize-inner-gpu-type

CHANGELOG.md

@@ -2,9 +2,15 @@
 
 - `permute` will be a wrapper around to `QuantumInterface.permutesubsystems`. Documentation for `permute!` would be similarly updated
 - reworking the rest of `NoisyCircuits` and moving it out of `Experimental`
+- `random_pauli(...; realphase=true)` should be the default
 
 # News
 
+## v0.8.16 - 2023-08-31
+
+- **(breaking)** Changes to the circuit generation in the non-public ECC module. Adding a Shor syndrome measurement circuit.
+- Added a `ClassicalXOR` gate, for now supporting only `PauliFrame`s.
+
 ## v0.8.15 - 2023-08-16
 
 - Initial support for GPU accelerated circuit simulation (with CUDA).

Project.toml

@@ -47,7 +47,7 @@ Makie = "0.19.7"
 Nemo = "0.34, 0.35"
 Plots = "1.38.0"
 PrecompileTools = "1"
-Quantikz = "1.3"
+Quantikz = "1.3.1"
 QuantumInterface = "0.3.0"
 QuantumOpticsBase = "0.4"
 SIMD = "3.4.0"

ext/QuantumCliffordGPUExt/apply_noise.jl

@@ -18,16 +18,17 @@ end
 
 function applynoise_kernel(xzs::DeviceMatrix{Tme},
     p::Real,
-    ibig::Int, # todo why is this Int and others Tme?
-    ismallm::Tme, # why is rows Unsigned sometimes and sometimes Int?
+    ibig::Int,
+    ismallm::Tme,
     rows::Unsigned) where {Tme <: Unsigned} 
 
     f = (blockIdx().x - 1) * blockDim().x + threadIdx().x;
     if f > rows
         return nothing
     end
 
-    r = rand() # todo add type to it? is it uniform?
+
+    r = rand()
     if  r < p # X error
         xzs[ibig,f] ⊻= ismallm
     elseif r < 2p # Z error

ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl

@@ -72,6 +72,7 @@ function Quantikz.QuantikzOp(op::Reset) # TODO This is complicated because quant
 end
 Quantikz.QuantikzOp(op::NoiseOp) = Quantikz.Noise(op.indices)
 Quantikz.QuantikzOp(op::NoiseOpAll) = Quantikz.NoiseAll()
+Quantikz.QuantikzOp(op::ClassicalXOR) = Quantikz.ClassicalDecision(sort([op.store, op.bits...]))
 
 function lstring(pauli::PauliOperator)
     v = join(("\\mathtt{$(o)}" for o in replace(string(pauli)[3:end],"_"=>"I")),"\\\\")

src/QuantumClifford.jl

@@ -53,7 +53,7 @@ export
     # Misc Ops
     SparseGate,
     sMX, sMY, sMZ, PauliMeasurement, Reset, sMRX, sMRY, sMRZ,
-    BellMeasurement,
+    BellMeasurement, ClassicalXOR,
     VerifyOp,
     Register,
     # Enumeration and Randoms

src/ecc/ECC.jl

@@ -44,15 +44,26 @@ end
 """Generate the non-fault-tolerant stabilizer measurement cicuit for a given code instance or parity check tableau.
 
 Use the `ancillary_index` and `bit_index` arguments to offset where the corresponding part the circuit starts.
+
+Returns the circuit, the number of ancillary qubits that were added, and a list of bit indices that will store the measurement results.
+
+See also: [`shor_syndrome_circuit`](@ref)
 """
 function naive_syndrome_circuit end
 
 function naive_syndrome_circuit(code_type::AbstractECC, ancillary_index=1, bit_index=1)
     naive_syndrome_circuit(parity_checks(code_type), ancillary_index, bit_index)
 end
 
-"""Circuit that measures the corresponding PauliOperator by using conditional gates into an ancillary
-qubit at index `nqubits(p)+ancillary_index` and stores the measurement result into classical bit `bit_index`."""
+"""Naive Pauli measurement circuit using a single ancillary qubit.
+
+Not a fault-tolerant circuit, but useful for testing purposes.
+
+Measures the corresponding `PauliOperator` by using conditional gates into an ancillary
+qubit at index `nqubits(p)+ancillary_index` and stores the measurement result
+into classical bit `bit_index`.
+
+See also: [`naive_syndrome_circuit`](@ref), [`shor_ancillary_paulimeasurement`](@ref)"""
 function naive_ancillary_paulimeasurement(p::PauliOperator, ancillary_index=1, bit_index=1)
     circuit = AbstractOperation[]
     numQubits = nqubits(p)
@@ -73,14 +84,106 @@ end
 
 function naive_syndrome_circuit(parity_check_tableau, ancillary_index=1, bit_index=1)
     naive_sc = AbstractOperation[]
+    ancillaries = 0
+    bits = 0
+    for check in parity_check_tableau
+        append!(naive_sc,naive_ancillary_paulimeasurement(check, ancillary_index+ancillaries, bit_index+bits))
+        ancillaries +=1
+        bits +=1
+    end
+
+    return naive_sc, ancillaries, bit_index:bit_index+bits-1
+end
+
+"""Generate the Shor fault-tolerant stabilizer measurement cicuit for a given code instance or parity check tableau.
+
+Use the `ancillary_index` and `bit_index` arguments to offset where the corresponding part the circuit starts.
+Ancillary qubits
+
+Returns:
+  - The cat state preparation circuit.
+  - The Shor syndrome measurement circuit.
+  - The number of ancillary qubits that were added.
+  - The list of bit indices that store the final measurement results.
+
+See also: [`naive_syndrome_circuit`](@ref)
+"""
+function shor_syndrome_circuit end
+
+function shor_syndrome_circuit(code_type::AbstractECC, ancillary_index=1, bit_index=1)
+    shor_syndrome_circuit(parity_checks(code_type), ancillary_index, bit_index)
+end
+
+"""Shor's Pauli measurement circuit using a multiple entangled ancillary qubits.
+
+A fault-tolerant circuit.
+
+Measures the corresponding PauliOperator by using conditional gates into multiple ancillary
+entangled qubits starting at index `nqubits(p)+ancillary_index`
+and stores the measurement result into classical bits starting at `bit_index`.
+The final measurement result is the XOR of all the bits.
+
+Returns:
+  - The cat state preparation circuit.
+  - The Shor syndrome measurement circuit.
+  - One more than the index of the last added ancillary qubit.
+  - One more than the index of the last added classical bit.
+
+See also: [`naive_syndrome_circuit`](@ref), [`naive_ancillary_paulimeasurement`](@ref)"""
+function shor_ancillary_paulimeasurement(p::PauliOperator, ancillary_index=1, bit_index=1)
+    init_ancil_index = ancillary_index
+    circuit = AbstractOperation[]
+    measurements = AbstractOperation[]
+    numQubits = nqubits(p)
+    for qubit in 1:numQubits
+        if p[qubit] == (1,0)
+            push!(circuit, sXCZ(qubit, numQubits + ancillary_index))
+        elseif p[qubit] == (0,1)
+            push!(circuit, sZCZ(qubit, numQubits + ancillary_index))
+        elseif p[qubit] == (1,1)
+            push!(circuit, sYCZ(qubit, numQubits + ancillary_index))
+        end
+        if p[qubit] != (0,0)
+            push!(measurements, sMRX(numQubits + ancillary_index, bit_index))
+            ancillary_index +=1
+            bit_index +=1
+        end
+    end
+    circuit = vcat(circuit,measurements)
+
+    cat_state_circuit = AbstractOperation[]
+    push!(cat_state_circuit, sHadamard(numQubits + init_ancil_index))
+    for i in (init_ancil_index+1):(ancillary_index -1)
+        push!(cat_state_circuit, sCNOT(numQubits + (i - 1), numQubits + i))
+    end
+
+    return cat_state_circuit, circuit, ancillary_index, bit_index
+end
+
+function shor_syndrome_circuit(parity_check_tableau, ancillary_index=1, bit_index=1)
+    shor_sc = AbstractOperation[]
+    xor_indices = []
+    cat_circuit = AbstractOperation[]
+    initial_ancillary_index = ancillary_index
 
     for check in parity_check_tableau
-        naive_sc = append!(naive_sc,naive_ancillary_paulimeasurement(check, ancillary_index, bit_index))
-        ancillary_index +=1
-        bit_index +=1
+        cat_circ, circ, new_ancillary_index, new_bit_index = shor_ancillary_paulimeasurement(check, ancillary_index, bit_index)
+        push!(xor_indices, collect(bit_index:new_bit_index-1))
+
+        append!(shor_sc, circ)
+        append!(cat_circuit, cat_circ)
+
+        ancillary_index = new_ancillary_index
+        bit_index = new_bit_index
+    end
+
+    final_bits = 0
+    for indices in xor_indices
+        push!(shor_sc, QuantumClifford.ClassicalXOR(indices, bit_index+final_bits))
+        final_bits += 1
     end
 
-    return naive_sc
+    return cat_circuit, shor_sc, ancillary_index-initial_ancillary_index, bit_index:bit_index+final_bits-1
 end
 
 """The distance of a code."""

src/misc_ops.jl

@@ -126,3 +126,12 @@ function applywstatus!(s::AbstractQCState, v::VerifyOp) # XXX It assumes the oth
 end
 
 operatordeterminism(::Type{VerifyOp}) = DeterministicOperatorTrait()
+
+"""Applies an XOR gate to classical bits. Currently only implemented for funcitonality with pauli frames."""
+struct ClassicalXOR{N} <: AbstractOperation
+    "The indices of the classical bits to be xor-ed"
+    bits::NTuple{N,Int}
+    "The index of the classical bit that will store the results"
+    store::Int
+end
+ClassicalXOR(bits::Vector, store::Int) = ClassicalXOR(tuple(bits...),store)

src/pauli_frames.jl

@@ -57,7 +57,28 @@ function apply!(f::PauliFrame, op::AbstractCliffordOperator)
     return f
 end
 
-function apply!(frame::PauliFrame, op::sMZ) # TODO sMX, sMY
+function apply!(frame::PauliFrame, xor::ClassicalXOR)
+    for f in eachindex(frame)
+        value = frame.measurements[f,xor.bits[1]]
+        for i in xor.bits[2:end]
+            value ⊻= frame.measurements[f,i]
+        end
+        frame.measurements[f, xor.store] = value
+    end
+end
+
+function apply!(frame::PauliFrame, op::sMX) # TODO implement a faster direct version
+    op.bit == 0 && return frame
+    apply!(frame, sHadamard(op.qubit))
+    apply!(frame, sMZ(op.qubit, op.bit))
+end
+
+function apply!(frame::PauliFrame, op::sMRX) # TODO implement a faster direct version
+    apply!(frame, sHadamard(op.qubit))
+    apply!(frame, sMRZ(op.qubit, op.bit))
+end
+
+function apply!(frame::PauliFrame, op::sMZ) # TODO sMY, and faster sMX
     op.bit == 0 && return frame
     i = op.qubit
     xzs = frame.frame.tab.xzs
@@ -75,7 +96,7 @@ function apply!(frame::PauliFrame, op::sMZ) # TODO sMX, sMY
     return frame
 end
 
-function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRX, sMRY
+function apply!(frame::PauliFrame, op::sMRZ) # TODO sMRY, and faster sMRX
     i = op.qubit
     xzs = frame.frame.tab.xzs
     T = eltype(xzs)

test/runtests.jl

@@ -60,6 +60,7 @@ end
 @doset "enumerate"
 @doset "quantumoptics"
 @doset "ecc"
+@doset "ecc_syndromes"
 @doset "precompile"
 @doset "pauliframe"
 @doset "allocations"

test/test_ecc.jl

@@ -67,7 +67,7 @@ function test_naive_syndrome(c::AbstractECC)
     s1 = copy(logicalqubits)
     syndrome1 = [project!(s1, check)[3] for check in parity_checks(c)]
     # measure using `naive_syndrome_circuit`
-    naive_circuit = naive_syndrome_circuit(c)
+    naive_circuit, _ = naive_syndrome_circuit(c)
     ancillaryqubits = one(Stabilizer,code_s(c))
     s2 = copy(logicalqubits)
     syndrome2 = Register(s2⊗ancillaryqubits, falses(code_s(c)))
@@ -89,7 +89,7 @@ end
 
 function test_with_pframes(code)
     ecirc = encoding_circuit(code)
-    scirc = naive_syndrome_circuit(code)
+    scirc, _ = naive_syndrome_circuit(code)
     nframes = 10
     dataqubits = code_n(code)
     ancqubits = code_s(code)

test/test_ecc_syndromes.jl

@@ -0,0 +1,54 @@
+using Test
+using QuantumClifford
+using QuantumClifford: mul_left!
+using QuantumClifford.ECC: AbstractECC, Steane7, Shor9, Bitflip3, naive_syndrome_circuit, encoding_circuit, shor_syndrome_circuit, code_n, code_s
+
+codes = [
+    Bitflip3(),
+    Steane7(),
+    Shor9(),
+]
+
+##
+
+function pframe_naive_vs_shor_syndrome(code)
+    ecirc = encoding_circuit(code)
+    naive_scirc, naive_ancillaries = naive_syndrome_circuit(code)
+    shor_cat_scirc, shor_scirc, shor_ancillaries, shor_bits = shor_syndrome_circuit(code)
+    nframes = 10
+    dataqubits = code_n(code)
+    syndromebits = code_s(code)
+    naive_qubits = dataqubits + syndromebits
+    shor_qubits = dataqubits + shor_ancillaries
+    # no noise
+    naive_frames = PauliFrame(nframes, naive_qubits, syndromebits)
+    shor_frames = PauliFrame(nframes, shor_qubits, last(shor_bits))
+    naive_circuit = [ecirc..., naive_scirc...]
+    shor_circuit = [ecirc..., shor_cat_scirc..., shor_scirc...]
+    pftrajectories(naive_frames, naive_circuit)
+    pftrajectories(shor_frames, shor_circuit)
+    @test pfmeasurements(naive_frames) == pfmeasurements(shor_frames)[:,shor_bits]
+    # with errors
+    for _ in 1:10
+        naive_frames = PauliFrame(nframes, naive_qubits, syndromebits)
+        shor_frames = PauliFrame(nframes, shor_qubits, last(shor_bits))
+        pftrajectories(naive_frames, ecirc)
+        pftrajectories(shor_frames, [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)
+        pn = embed(naive_qubits, 1:dataqubits, p)
+        ps = embed(shor_qubits, 1:dataqubits, p)
+        mul_left!(naive_frames.frame, pn)
+        mul_left!(shor_frames.frame, ps)
+        # 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]
+    end
+end
+
+@testset "naive and shor measurement circuits" begin
+    for c in codes
+        pframe_naive_vs_shor_syndrome(c)
+    end
+end

test/test_jet.jl

@@ -36,5 +36,5 @@ end
     )
     @show rep
     @test_broken length(JET.get_reports(rep)) == 0
-    @test length(JET.get_reports(rep)) <= 7
+    @test length(JET.get_reports(rep)) <= 9
 end