From 136ae94ddb04279888073050cf433fceac1cd52d Mon Sep 17 00:00:00 2001 From: Stefan Krastanov Date: Thu, 22 Jun 2023 01:32:20 -0400 Subject: [PATCH] UX improvements to pftrajectories (#136) --- CHANGELOG.md | 6 ++ Project.toml | 2 +- docs/src/ecc_example_sim.md | 26 +++---- .../QuantumCliffordQuantikzExt.jl | 30 ++++---- src/QuantumClifford.jl | 8 ++- src/affectedqubits.jl | 12 ++-- src/dense_cliffords.jl | 2 +- src/experimental/NoisyCircuits.jl | 2 +- src/mctrajectory.jl | 4 +- src/pauli_frames.jl | 51 ++++++++++++- src/sumtypes.jl | 2 + src/symbolic_cliffords.jl | 12 +++- test/test_pauliframe.jl | 71 +++++++++++++------ 13 files changed, 157 insertions(+), 71 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 630597930..f63730503 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,12 @@ # News +## v0.8.7 - 2023-06-22 + +- Better UX and threading support for `pftrajectories`. +- `affectedqubits` now more consistently returns tuples instead of vectors. +- Private `affectedbits` is now implemented. +- Many operation constructors now throw an error if they are given negative qubit indices. ## v0.8.6 - 2023-06-20 - Fixes to Quantikz circuit plotting of empty circuits and `PauliOperator` diff --git a/Project.toml b/Project.toml index e49adab05..b6aa42b17 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumClifford" uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1" authors = ["Stefan Krastanov "] -version = "0.8.6" +version = "0.8.7" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" diff --git a/docs/src/ecc_example_sim.md b/docs/src/ecc_example_sim.md index 0db8ca546..5feef8edf 100644 --- a/docs/src/ecc_example_sim.md +++ b/docs/src/ecc_example_sim.md @@ -32,24 +32,15 @@ scirc = naive_syndrome_circuit(code) The most straightforward way to start sampling syndromes is to set up a table of Pauli frames. ```@example 1 -nframes = 4 -dataqubits = code_n(code) -ancqubits = code_s(code) -regbits = ancqubits -frames = PauliFrame(nframes, dataqubits+ancqubits, regbits) circuit = [ecirc..., scirc...] -pftrajectories(frames, circuit) # run the sims -pfmeasurements(frames) # extract the measurements +nframes = 4 +frames = pftrajectories(circuit; trajectories=nframes) # run the sims +pfmeasurements(frames) # extract the measurements ``` -To avoid runtime dispatch, it is a good practice to compactify the circuit: - -```@example 1 -frames = PauliFrame(nframes, dataqubits+ancqubits, regbits) -ccircuit = compactify_circuit(circuit) -pftrajectories(frames, ccircuit) -pfmeasurements(frames) -``` +The [`pftrajectories`](@ref) function is multithreaded. +If you want more low-level control over these Pauli frame simulations, check out the [`PauliFrame`](@ref) structure, +the other methods of [`pftrajectories`](@ref), and the circuit compactifaction function [`compactify_circuit`](@ref). If you want to model Pauli errors, use: @@ -58,13 +49,12 @@ If you want to model Pauli errors, use: ```@example 1 errprob = 0.1 -errors = [PauliError(i,errprob) for i in 1:dataqubits] +errors = [PauliError(i,errprob) for i in 1:code_n(code)] fullcircuit = [ecirc..., errors..., scirc...] ``` And running this noisy simulation: ```@example 1 -frames = PauliFrame(nframes, dataqubits+ancqubits, regbits) -pftrajectories(frames, fullcircuit) +frames = pftrajectories(fullcircuit; trajectories=nframes) pfmeasurements(frames) ``` \ No newline at end of file diff --git a/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl b/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl index db9ea6d04..d65eb6e47 100644 --- a/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl +++ b/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl @@ -19,8 +19,8 @@ function Quantikz.QuantikzOp(op::SparseGate) return Quantikz.MultiControlU([],[],op.indices) # TODO Permit skipping the string end end -Quantikz.QuantikzOp(op::AbstractOperation) = Quantikz.MultiControlU(affectedqubits(op)) -Quantikz.QuantikzOp(op::PauliOperator) = Quantikz.MultiControlU("\\begin{array}{c}$(lstring(op))\\end{array}", affectedqubits(op)) +Quantikz.QuantikzOp(op::AbstractOperation) = Quantikz.MultiControlU(collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect +Quantikz.QuantikzOp(op::PauliOperator) = Quantikz.MultiControlU("\\begin{array}{c}$(lstring(op))\\end{array}", collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect Quantikz.QuantikzOp(op::sId1) = Quantikz.Id(affectedqubits(op)...) function Quantikz.QuantikzOp(op::AbstractSingleQubitOperator) T = typeof(op) @@ -40,37 +40,37 @@ Quantikz.QuantikzOp(op::sXCZ) = Quantikz.CNOT(reverse(affectedqubits(op))...) Quantikz.QuantikzOp(op::sCPHASE) = Quantikz.CPHASE(affectedqubits(op)...) Quantikz.QuantikzOp(op::sSWAP) = Quantikz.SWAP(affectedqubits(op)...) Quantikz.QuantikzOp(op::sZCZ) = Quantikz.CPHASE(affectedqubits(op)...) -Quantikz.QuantikzOp(op::sXCX) = Quantikz.MultiControl([],[],affectedqubits(op),[]) +Quantikz.QuantikzOp(op::sXCX) = Quantikz.MultiControl([],[],collect(affectedqubits(op)),[]) # TODO make Quantikz work with tuples and remove the collect Quantikz.QuantikzOp(op::sZCY) = Quantikz.MultiControlU("Y",affectedqubits(op)...) Quantikz.QuantikzOp(op::sYCZ) = Quantikz.MultiControlU("Y",reverse(affectedqubits(op))...) function Quantikz.QuantikzOp(op::AbstractTwoQubitOperator) T = typeof(op) if T in (sXCY,sYCX,sYCY,sYCZ,sZCY) - return Quantikz.U(string(T),affectedqubits(op)) + return Quantikz.U(string(T),collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect else - return Quantikz.U("{U}",affectedqubits(op)...) + return Quantikz.U("{U}",collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect end end -Quantikz.QuantikzOp(op::BellMeasurement) = Quantikz.ParityMeasurement(["\\mathtt{$(string(typeof(o))[3])}" for o in op.measurements], affectedqubits(op)) +Quantikz.QuantikzOp(op::BellMeasurement) = Quantikz.ParityMeasurement(["\\mathtt{$(string(typeof(o))[3])}" for o in op.measurements], collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect Quantikz.QuantikzOp(op::NoisyBellMeasurement) = Quantikz.QuantikzOp(op.meas) -Quantikz.QuantikzOp(op::ConditionalGate) = Quantikz.ClassicalDecision(affectedqubits(op),op.controlbit) -Quantikz.QuantikzOp(op::DecisionGate) = Quantikz.ClassicalDecision(affectedqubits(op),Quantikz.ibegin:Quantikz.iend) +Quantikz.QuantikzOp(op::ConditionalGate) = Quantikz.ClassicalDecision(affectedqubits(op),op.controlbit) # TODO make Quantikz work with tuples and remove the collect +Quantikz.QuantikzOp(op::DecisionGate) = Quantikz.ClassicalDecision(affectedqubits(op),Quantikz.ibegin:Quantikz.iend) # TODO make Quantikz work with tuples and remove the collect #Quantikz.QuantikzOp(op::DenseGate) = Quantikz.MultiControlU(affectedqubits(op)) -Quantikz.QuantikzOp(op::PauliMeasurement) = Quantikz.Measurement("\\begin{array}{c}$(lstring(op.pauli))\\end{array}",affectedqubits(op),op.bit) -Quantikz.QuantikzOp(op::sMX) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{X}\\end{array}",affectedqubits(op), iszero(op.bit) ? nothing : op.bit) -Quantikz.QuantikzOp(op::sMY) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Y}\\end{array}",affectedqubits(op), iszero(op.bit) ? nothing : op.bit) -Quantikz.QuantikzOp(op::sMZ) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Z}\\end{array}",affectedqubits(op), iszero(op.bit) ? nothing : op.bit) +Quantikz.QuantikzOp(op::PauliMeasurement) = Quantikz.Measurement("\\begin{array}{c}$(lstring(op.pauli))\\end{array}",collect(affectedqubits(op)),op.bit) # TODO make Quantikz work with tuples and remove the collect +Quantikz.QuantikzOp(op::Union{sMX,sMRX}) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{X}\\end{array}",collect(affectedqubits(op)), iszero(op.bit) ? nothing : op.bit) # TODO make Quantikz work with tuples and remove the collect +Quantikz.QuantikzOp(op::Union{sMY,sMRY}) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Y}\\end{array}",collect(affectedqubits(op)), iszero(op.bit) ? nothing : op.bit) # TODO make Quantikz work with tuples and remove the collect +Quantikz.QuantikzOp(op::Union{sMZ,sMRZ}) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Z}\\end{array}",collect(affectedqubits(op)), iszero(op.bit) ? nothing : op.bit) # TODO make Quantikz work with tuples and remove the collect #Quantikz.QuantikzOp(op::SparseMeasurement) = Quantikz.Measurement("\\begin{array}{c}$(lstring(op.pauli))\\end{array}",affectedqubits(op),op.bit) Quantikz.QuantikzOp(op::NoisyGate) = Quantikz.QuantikzOp(op.gate) -Quantikz.QuantikzOp(op::VerifyOp) = Quantikz.MultiControlU("\\begin{array}{c}\\mathrm{Verify:}\\\\$(lstring(op.good_state))\\end{array}",affectedqubits(op)) +Quantikz.QuantikzOp(op::VerifyOp) = Quantikz.MultiControlU("\\begin{array}{c}\\mathrm{Verify:}\\\\$(lstring(op.good_state))\\end{array}",collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect function Quantikz.QuantikzOp(op::Reset) # TODO This is complicated because quantikz requires $$ for some operators but not all of them... Fix in Quantikz.jl m,M = extrema(op.indices) indices = sort(op.indices) str = "\\begin{array}{c}\\\\$(lstring(op.resetto))\\end{array}" if collect(m:M)==indices - Quantikz.Initialize("\$$str\$",affectedqubits(op)) + Quantikz.Initialize("\$$str\$",affectedqubits(op)) # TODO make Quantikz work with tuples and remove the collect else - Quantikz.Initialize("$str",affectedqubits(op)) + Quantikz.Initialize("$str",affectedqubits(op)) # TODO make Quantikz work with tuples and remove the collect end end Quantikz.QuantikzOp(op::NoiseOp) = Quantikz.Noise(op.indices) diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl index 32831bb24..efbc87a14 100644 --- a/src/QuantumClifford.jl +++ b/src/QuantumClifford.jl @@ -27,7 +27,7 @@ export stabilizerview, destabilizerview, logicalxview, logicalzview, phases, bitview, quantumstate, tab, BadDataStructure, - affectedqubits, + affectedqubits, #TODO move to QuantumInterface? # GF2 stab_to_gf2, gf2_gausselim!, gf2_isinvertible, gf2_invert, gf2_H_to_G, # Canonicalization @@ -91,6 +91,12 @@ function __init__() BIG_INT_FOUR[] = BigInt(4) end +const MINBATCH1Q = 100 +const MINBATCH2Q = 100 +const MINBATCHDENSE = 25 + +const NoZeroQubit = ArgumentError("Qubit indices have to be larger than zero, but you attempting are creating a gate acting on a qubit with a non-positive index. Ensure indexing always starts from 1.") + # Predefined constants representing the permitted phases encoded # in the low bits of UInt8. const _p = 0x00 diff --git a/src/affectedqubits.jl b/src/affectedqubits.jl index 712a7f712..f27bc2469 100644 --- a/src/affectedqubits.jl +++ b/src/affectedqubits.jl @@ -1,12 +1,16 @@ """A method giving the qubits acted upon by a given operation. Part of the Noise interface.""" function affectedqubits end -affectedqubits(g::AbstractSingleQubitOperator) = [g.q,] -affectedqubits(g::AbstractTwoQubitOperator) = [g.q1, g.q2] +affectedqubits(g::AbstractSingleQubitOperator) = (g.q) +affectedqubits(g::AbstractTwoQubitOperator) = (g.q1, g.q2) affectedqubits(g::NoisyGate) = affectedqubits(g.gate) affectedqubits(g::SparseGate) = g.indices -affectedqubits(b::BellMeasurement) = [m.qubit for m in b.measurements] +affectedqubits(b::BellMeasurement) = map(m->m.qubit, b.measurements) affectedqubits(r::Reset) = r.indices affectedqubits(n::NoiseOp) = n.indices affectedqubits(g::PauliMeasurement) = 1:length(g.pauli) affectedqubits(p::PauliOperator) = 1:length(p) -affectedqubits(m::AbstractMeasurement) = [m.qubit] +affectedqubits(m::Union{AbstractMeasurement,sMRX,sMRY,sMRZ}) = (m.qubit,) + +affectedbits(o) = () +affectedbits(m::sMRZ) = (m.bit,) +affectedbits(m::sMZ) = (m.bit,) diff --git a/src/dense_cliffords.jl b/src/dense_cliffords.jl index 7c810fcb5..a67bd2013 100644 --- a/src/dense_cliffords.jl +++ b/src/dense_cliffords.jl @@ -127,7 +127,7 @@ function _apply!(stab::AbstractStabilizer, c::CliffordOperator; phases::Val{B}=V nqubits(stab)==nqubits(c) || throw(DimensionMismatch("The tableau and the Clifford operator need to act on the same number of qubits. Consider specifying an array of indices as a third argument to the `apply!` function to avoid this error.")) s_tab = tab(stab) c_tab = tab(c) - @batch minbatch=25 threadlocal=zero(c_tab[1]) for row_stab in eachindex(s_tab) + @batch minbatch=MINBATCHDENSE threadlocal=zero(c_tab[1]) for row_stab in eachindex(s_tab) zero!(threadlocal) # a new stabrow for temporary storage apply_row_kernel!(threadlocal, row_stab, s_tab, c_tab, phases=phases) end diff --git a/src/experimental/NoisyCircuits.jl b/src/experimental/NoisyCircuits.jl index 3cee5d18d..8dc61c64f 100644 --- a/src/experimental/NoisyCircuits.jl +++ b/src/experimental/NoisyCircuits.jl @@ -64,7 +64,7 @@ end """Compute all possible new states after the application of the given noise model. Reports the probability of each one of them. Deterministic, part of the Noise interface.""" function applynoise_branches end -function applynoise_branches(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,indices::AbstractVector{Int}; max_order=1) +function applynoise_branches(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,indices; max_order=1) n = nqubits(s) l = length(indices) infid = noise.errprobthird diff --git a/src/mctrajectory.jl b/src/mctrajectory.jl index bc3ff6157..c6bf74404 100644 --- a/src/mctrajectory.jl +++ b/src/mctrajectory.jl @@ -56,7 +56,9 @@ function countmap(samples::Vector{CircuitStatus}) # A simpler faster version of Dict(CircuitStatus(i-1)=>counts[i] for i in eachindex(counts)) end -"""Run multiple Monte Carlo trajectories and report the aggregate final statuses of each.""" +"""Run multiple Monte Carlo trajectories and report the aggregate final statuses of each. + +See also: [`pftrajectories`](@ref), `petrajectories`""" function mctrajectories(initialstate,circuit;trajectories=500) counts = countmap([mctrajectory!(copy(initialstate),circuit)[2] for i in 1:trajectories]) # TODO use threads or at least a generator, but without breaking Polyester return counts diff --git a/src/pauli_frames.jl b/src/pauli_frames.jl index b41696a66..fa6c96f7d 100644 --- a/src/pauli_frames.jl +++ b/src/pauli_frames.jl @@ -6,15 +6,16 @@ although it does conjugate the same under Clifford operations. Each row in the tableau refers to a single frame. The row represents the Pauli operation by which the frame and the reference differ. """ -struct PauliFrame{T} <: AbstractQCState - frame::T # TODO this should really be a Tableau - measurements::Matrix{Bool} +struct PauliFrame{T,S} <: AbstractQCState + frame::T # TODO this should really be a Tableau, but now most of the code seems to be assuming a Stabilizer + measurements::S # TODO check if when looping over this we are actually looping over the fast axis end nqubits(f::PauliFrame) = nqubits(f.frame) 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, :)) """ $(TYPEDSIGNATURES) @@ -124,6 +125,50 @@ function pftrajectories end """ $(TYPEDSIGNATURES) +The main method for running Pauli frame simulations of circuits. +See the other methods for lower level access. + +Multithreading is enabled by default, but can be disabled by setting `threads=false`. +Do not forget to launch Julia with multiple threads enabled, e.g. `julia -t4`, if you want +to use multithreading. + +Note for advanced users: Much of the underlying QuantumClifford.jl functionaly is capable of +using Polyester.jl threads, but they are fully dissabled here as this is an embarassingly +parallel problem. If you want to use Polyester.jl threads, use the lower level methods. +The `threads` keyword argument controls whether standard Julia threads are used. +""" +function pftrajectories(circuit;trajectories=5000,threads=true) + Polyester.disable_polyester_threads() do + _pftrajectories(circuit;trajectories,threads) + end +end + +function _pftrajectories(circuit;trajectories=5000,threads=true) + ccircuit = if eltype(circuit) <: CompactifiedGate + circuit + else + compactify_circuit(circuit) + end + qmax=maximum((maximum(affectedqubits(g)) for g in ccircuit)) + bmax=maximum((maximum(affectedbits(g),init=1) for g in ccircuit)) + frames = PauliFrame(trajectories, qmax, bmax) + nthr = min(Threads.nthreads(),trajectories÷(MINBATCH1Q)) + if threads && nthr>1 + batchsize = trajectories÷nthr + Threads.@threads for i in 1:nthr + b = (i-1)*batchsize+1 + e = i==nthr ? trajectories : i*batchsize + pftrajectories((@view frames[b:e]), ccircuit) + end + else + pftrajectories(frames, ccircuit) + end + return frames +end + +""" +$(TYPEDSIGNATURES) + Evolve each frame stored in [`PauliFrame`](@ref) by the given circuit. """ function pftrajectories(state::PauliFrame, circuit) diff --git a/src/sumtypes.jl b/src/sumtypes.jl index 2f6167196..a3824c9d6 100644 --- a/src/sumtypes.jl +++ b/src/sumtypes.jl @@ -114,6 +114,8 @@ function make_all_sumtype_infrastructure() (:applywstatus!, (:(s::Register),), ()), (:apply!, (:(s::PauliFrame),), ()), (:applywstatus!, (:(s::PauliFrame),), ()), + (:affectedqubits, (), ()), + (:affectedbits, (), ()), ] ) |> eval end diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl index 8972d7e0a..8f84fd5a5 100644 --- a/src/symbolic_cliffords.jl +++ b/src/symbolic_cliffords.jl @@ -9,9 +9,6 @@ abstract type AbstractTwoQubitOperator <: AbstractSymbolicOperator end """Supertype of all symbolic single-qubit measurements.""" abstract type AbstractMeasurement <: AbstractOperation end -const MINBATCH1Q = 100 -const MINBATCH2Q = 100 - # Stim has a good list of specialized single and two qubit operations at https://github.com/quantumlib/Stim/blob/e51ea66d213b25920e72c08e53266ec56fd14db4/src/stim/stabilizers/tableau_specialized_prepend.cc # Note that their specialized operations are for prepends (right multiplications), while we implement append (left multiplication) operations. @@ -67,6 +64,7 @@ macro qubitop1(name, kernel) quote struct $(esc(prefixname)) <: AbstractSingleQubitOperator q::Int + $(esc(prefixname))(q) = if q<=0 throw(NoZeroQubit) else new(q) end end @doc $docstring $prefixname @inline $(esc(:qubit_kernel))(::$prefixname, x, z) = $kernel @@ -132,6 +130,7 @@ struct SingleQubitOperator <: AbstractSingleQubitOperator zz::Bool px::Bool pz::Bool + SingleQubitOperator(q,args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end function _apply!(stab::AbstractStabilizer, op::SingleQubitOperator; phases::Val{B}=Val(true)) where B # TODO Generated functions that simplify the whole `if phases` branch might be a good optimization, but a quick benchmakr comparing sHadamard to SingleQubitOperator(sHadamard) did not show a worthwhile difference. s = tab(stab) @@ -249,6 +248,7 @@ macro qubitop2(name, kernel) struct $(esc(prefixname)) <: AbstractTwoQubitOperator q1::Int q2::Int + $(esc(prefixname))(q1,q2) = if q1<=0 || q2<=0 throw(NoZeroQubit) else new(q1,q2) end end @doc $docstring $prefixname @inline $(esc(:qubit_kernel))(::$prefixname, x1, z1, x2, z2) = $kernel @@ -363,18 +363,21 @@ end struct sMX <: AbstractMeasurement qubit::Int bit::Int + sMX(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end """Symbolic single qubit Y measurement. See also [`Register`](@ref), [`projectYrand!`](@ref), [`sMX`](@ref), [`sMZ`](@ref)""" struct sMY <: AbstractMeasurement qubit::Int bit::Int + sMY(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end """Symbolic single qubit Z measurement. See also [`Register`](@ref), [`projectZrand!`](@ref), [`sMX`](@ref), [`sMY`](@ref)""" struct sMZ <: AbstractMeasurement qubit::Int bit::Int + sMZ(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end sMX(i) = sMX(i,0) @@ -471,6 +474,7 @@ See also: [`Reset`](@ref), [`sMZ`](@ref)""" struct sMRZ <: AbstractOperation qubit::Int bit::Int + sMRZ(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end """Measure a qubit in the X basis and reset to the |+⟩ state. @@ -479,6 +483,7 @@ See also: [`sMRZ`](@ref), [`Reset`](@ref), [`sMZ`](@ref)""" struct sMRX <: AbstractOperation qubit::Int bit::Int + sMRX(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end """Measure a qubit in the Y basis and reset to the |i₊⟩ state. @@ -487,6 +492,7 @@ See also: [`sMRZ`](@ref), [`Reset`](@ref), [`sMZ`](@ref)""" struct sMRY <: AbstractOperation qubit::Int bit::Int + sMRY(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end end sMRX(i) = sMRX(i,0) diff --git a/test/test_pauliframe.jl b/test/test_pauliframe.jl index 14d0ccd8d..0669f2ac2 100644 --- a/test/test_pauliframe.jl +++ b/test/test_pauliframe.jl @@ -10,11 +10,15 @@ using Test ] frame = PauliFrame(100, 5, 2) frame = pftrajectories(frame, circuit) - m = pfmeasurements(frame) + frame1 = pftrajectories(circuit; trajectories=100, threads=false) + frame2 = pftrajectories(circuit; trajectories=100, threads=true) # If the x component is set on the second qubit, then the fourth and fifth qubits should also have it set - for (i, row) in enumerate(frame.frame) - if row[2] == (true, false) - @test row[4][1] && row[5][1] && m[i,1] && m[i,2] + for f in [frame, frame1, frame2] + m = pfmeasurements(f) + for (i, row) in enumerate(f.frame) + if row[2] == (true, false) + @test row[4][1] && row[5][1] && m[i,1] && m[i,2] + end end end end @@ -24,44 +28,65 @@ end sHadamard(1), sCNOT(1,2), sCNOT(1,3), # prepare a GHZ state sMZ(1,1), sMZ(2,2), sMZ(3,3) # measure each qubit ] - frame = PauliFrame(10^6, 3, 3) + n = 10^6 + frame = PauliFrame(n, 3, 3) f = pftrajectories(frame, ghz_circuit) m = pfmeasurements(f) - rowtotal_1s = sum(m, dims=2)[:,1] - @test all(rowtotal_1s .% 3 .== 0) - fractotal_1s = sum(rowtotal_1s)/3 / 10^6 - @test (fractotal_1s > 0.49) && (fractotal_1s < 0.51) + frame1 = pftrajectories(ghz_circuit; trajectories=n, threads=false) + m1 = pfmeasurements(frame1) + frame2 = pftrajectories(ghz_circuit; trajectories=n, threads=true) + m2 = pfmeasurements(frame2) + for _m in [m, m1, m2] + rowtotal_1s = sum(m, dims=2)[:,1] + @test all(rowtotal_1s .% 3 .== 0) + fractotal_1s = sum(rowtotal_1s)/3 / n + @test (fractotal_1s > 0.49) && (fractotal_1s < 0.51) + end end -@testset "sMZ vs sMZR - mctrajectories vs pftrajectories" begin +@testset "sMZ vs sMRZ - mctrajectories vs pftrajectories" begin + n = 2000 state = Register(one(MixedDestabilizer, 3), 6) - frame = PauliFrame(100, 3, 6) + frame = PauliFrame(n, 3, 6) ghz_circuit1 = [ sHadamard(1), sCNOT(1,2), sCNOT(1,3), # prepare a GHZ state sMZ(1,1), sMZ(2,2), sMZ(3,3), # measure each qubit sMZ(1,4), sMZ(2,5), sMZ(3,6) # measure each qubit again ] - ms1 = stack([bitview(mctrajectory!(copy(state), ghz_circuit1)[1]) for i in 1:100], dims=1) - mf1 = pfmeasurements(pftrajectories(copy(frame), ghz_circuit1)) - @test all(0.25 .< sum(ms1, dims=1)./100 .< 0.75) - @test all(0.25 .< sum(mf1, dims=1)./100 .< 0.75) + for m in [ + stack([bitview(mctrajectory!(copy(state), ghz_circuit1)[1]) for i in 1:n], dims=1), + pfmeasurements(pftrajectories(copy(frame), ghz_circuit1)), + pfmeasurements(pftrajectories(ghz_circuit1;trajectories=n,threads=false)), + pfmeasurements(pftrajectories(ghz_circuit1;trajectories=n,threads=true)), + ] + @test all(0.25 .< sum(m, dims=1)./n .< 0.75) + end ghz_circuit2 = [ sHadamard(1), sCNOT(1,2), sCNOT(1,3), # prepare a GHZ state sMRZ(1,1), sMZ(2,2), sMZ(3,3), # measure and reset each qubit sMZ(1,4), sMZ(2,5), sMZ(3,6) # measure each qubit again ] - ms2 = stack([bitview(mctrajectory!(copy(state), ghz_circuit2)[1]) for i in 1:100], dims=1) - mf2 = pfmeasurements(pftrajectories(copy(frame), ghz_circuit2)) - @test all(0.25.*[1 1 1 0 1 1] .<= sum(ms2, dims=1)./100 .<= 0.75.*[1 1 1 0 1 1]) - @test all(0.25.*[1 1 1 0 1 1] .<= sum(mf2, dims=1)./100 .<= 0.75.*[1 1 1 0 1 1]) + for m in [ + stack([bitview(mctrajectory!(copy(state), ghz_circuit2)[1]) for i in 1:n], dims=1), + pfmeasurements(pftrajectories(copy(frame), ghz_circuit2)), + pfmeasurements(pftrajectories(ghz_circuit2;trajectories=n,threads=false)), + pfmeasurements(pftrajectories(ghz_circuit2;trajectories=n,threads=true)), + ] + @test all(0.25.*[1 1 1 0 1 1] .<= sum(m, dims=1)./n .<= 0.75.*[1 1 1 0 1 1]) + end noncom_circuit = [ sHadamard(1), sMRZ(1,1), sX(1), sMZ(1,2), sMRZ(1,3), sMRZ(1,4), sHadamard(1), sMZ(1,5) ] - ms3 = stack([bitview(mctrajectory!(copy(state), noncom_circuit)[1]) for i in 1:100], dims=1) - mf3 = pfmeasurements(pftrajectories(copy(frame), noncom_circuit)) - @test all(0.25.*[1 4 4 0 1 0] .<= sum(ms3, dims=1)./100 .<= 0.75.*[1 2 2 0 1 0]) - @test all(0.25.*[1 0 0 0 1 0] .<= sum(mf3, dims=1)./100 .<= 0.75.*[1 0 0 0 1 0]) + ms3 = stack([bitview(mctrajectory!(copy(state), noncom_circuit)[1]) for i in 1:n], dims=1) + @test all(0.25.*[1 4 4 0 1 0] .<= sum(ms3, dims=1)./n .<= 0.75.*[1 2 2 0 1 0]) + for m in [ + pfmeasurements(pftrajectories(copy(frame), noncom_circuit)), + pfmeasurements(pftrajectories(noncom_circuit;trajectories=n,threads=false)), + pfmeasurements(pftrajectories(noncom_circuit;trajectories=n,threads=true)), + ] + @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