UX improvements to pftrajectories

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`

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org>"]
-version = "0.8.6"
+version = "0.8.7"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

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:
 
@@ -64,7 +55,6 @@ 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)
 ```

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)

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

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,)

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

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

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`](@ref)`]"""
 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

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)

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

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)