UX improvements to pftrajectories (#136)

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:
 
@@ -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)
 ```

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`"""
 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)