QuantumSavory · Krastanov · Jul 10, 2023 · Jul 10, 2023 · Jul 10, 2023 · Jul 10, 2023
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,9 +2,16 @@
 
 - `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`
+- some support for T gates and other non-clifford behavior
 
 # News
 
+## v0.8.11 - 2023-07-10
+
+- `petrajectories`, for (potentially symbolic) perturbative expansions of the result of a circuit, is moved out of `Experimental` into the public part of the interface. The underlying `petrajectory` is not made public yet due to the ad-hoc low-level return conventions for it.
+- `mctrajectory` and `petrajectory` can now optionally report the end state of each trajectory, not just the circuit status (i.e. "success", "detected failure", etc).
+- Internally we now use a trait system to distinguish deterministic from non-deterministic operations. Not part of the public API yet.
+
 ## v0.8.10 - 2023-07-05
 
 - Remove Polyester.jl multithreading, leading to simpler and better compiled single-thread code. Now single-thread performance is much higher. Multithreading was generally not useful at the inner loops where it was deployed.

diff --git a/Project.toml b/Project.toml
@@ -5,6 +5,7 @@ version = "0.8.10"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HostCPUFeatures = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
@@ -31,6 +32,7 @@ QuantumCliffordQuantikzExt = "Quantikz"
 
 [compat]
 Combinatorics = "1.0"
+DataStructures = "0.18.0"
 DocStringExtensions = "0.8, 0.9"
 Graphs = "1.4.1"
 HostCPUFeatures = "0.1.6"

diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
@@ -9,6 +9,10 @@ module QuantumClifford
 
 import LinearAlgebra
 using LinearAlgebra: inv, mul!, rank, Adjoint
+import DataStructures
+using DataStructures: DefaultDict, Accumulator
+using Combinatorics: combinations
+using Base.Cartesian
 using DocStringExtensions
 
 import QuantumInterface: tensor, ⊗, tensor_pow, apply!, nqubits, expect, project!, reset_qubits!, traceout!, ptrace, apply!, projectX!, projectY!, projectZ!, entanglement_entropy
@@ -72,6 +76,7 @@ export
     # mctrajectories
     CircuitStatus, continue_stat, true_success_stat, false_success_stat, failure_stat,
     mctrajectory!, mctrajectories, applywstatus!,
+    petrajectories, applybranches,
     # makie plotting -- defined only when extension is loaded
     stabilizerplot, stabilizerplot_axis,
     # sum types
@@ -1250,26 +1255,37 @@ function mixed_destab_looks_good(destabilizer)
     return true
 end
 
+# base tableaux handling
 include("mul_leftright.jl")
 include("canonicalization.jl")
-include("dense_cliffords.jl")
 include("project_trace_reset.jl")
-include("linalg.jl")
+include("fastmemlayout.jl")
+# dense clifford operator tableaux
+include("dense_cliffords.jl")
+# special one- and two- qubit operators
 include("symbolic_cliffords.jl")
+include("linalg.jl")
+# circuits
+include("operator_traits.jl")
 include("mctrajectory.jl")
+include("petrajectory.jl")
 include("misc_ops.jl")
 include("classical_register.jl")
-include("enumeration.jl")
-include("randoms.jl")
-include("useful_states.jl")
 include("noise.jl")
 include("affectedqubits.jl")
 include("pauli_frames.jl")
+# common states and operators
+include("enumeration.jl")
+include("randoms.jl")
+include("useful_states.jl")
+#
 include("experimental/Experimental.jl")
+#
 include("graphs.jl")
+#
 include("entanglement.jl")
+#
 include("tableau_show.jl")
-include("fastmemlayout.jl")
 include("sumtypes.jl")
 include("precompiles.jl")
 include("ecc/ECC.jl")

diff --git a/src/affectedqubits.jl b/src/affectedqubits.jl
@@ -10,6 +10,7 @@ affectedqubits(n::NoiseOp) = n.indices
 affectedqubits(g::PauliMeasurement) = 1:length(g.pauli)
 affectedqubits(p::PauliOperator) = 1:length(p)
 affectedqubits(m::Union{AbstractMeasurement,sMRX,sMRY,sMRZ}) = (m.qubit,)
+affectedqubits(v::VerifyOp) = v.indices
 
 affectedbits(o) = ()
 affectedbits(m::sMRZ) = (m.bit,)

diff --git a/src/classical_register.jl b/src/classical_register.jl
@@ -114,3 +114,64 @@ function traceout!(r::Register, arg)
     traceout!(q,arg)
     q
 end
+
+##
+# petrajectories, applynoise_branches
+##
+
+function applynoise_branches(state::Register, noise, indices; max_order=1)
+    [(Register(newstate,copy(state.bits)), prob, order)
+     for (newstate, prob, order) in applynoise_branches(quantumstate(state), noise, indices; max_order=max_order)]
+end
+
+function applybranches(s::Register, op::PauliMeasurement; max_order=1) # TODO this is almost the same as `applybranches(s::Register, op::AbstractMeasurement; max_order=1)` defined below
+    stab = s.stab
+    stab, anticom, r = project!(stab, op.pauli)
+    new_branches = []
+    if isnothing(r)
+        s1 = s
+        phases(stabilizerview(s1.stab))[anticom] = 0x00
+        s1.bits[op.bit] = false
+        s2 = copy(s)
+        phases(stabilizerview(s2.stab))[anticom] = 0x02
+        s2.bits[op.bit] = true
+        push!(new_branches, (s1,continue_stat,1/2,0))
+        push!(new_branches, (s2,continue_stat,1/2,0))
+    else
+        s.bits[op.bit] = r==0x02
+        push!(new_branches, (s,continue_stat,1,0))
+    end
+    new_branches
+end
+
+function applybranches(s::Register, op::AbstractMeasurement; max_order=1)
+    stab = s.stab
+    stab, anticom, r = project!(stab, op)
+    new_branches = []
+    if isnothing(r)
+        s1 = s
+        phases(stabilizerview(s1.stab))[anticom] = 0x00
+        s1.bits[op.bit] = false
+        s2 = copy(s)
+        phases(stabilizerview(s2.stab))[anticom] = 0x02
+        s2.bits[op.bit] = true
+        push!(new_branches, (s1,continue_stat,1/2,0))
+        push!(new_branches, (s2,continue_stat,1/2,0))
+    else
+        s.bits[op.bit] = r==0x02
+        push!(new_branches, (s,continue_stat,1,0))
+    end
+    new_branches
+end
+
+#=
+function applybranches(state::Register, op::SparseMeasurement; max_order=1)
+    n = nqubits(state.stab) # TODO implement actual sparse measurements
+    p = zero(typeof(op.pauli), n)
+    for (ii,i) in enumerate(op.indices)
+        p[i] = op.pauli[ii]
+    end
+    dm = PauliMeasurement(p,op.bit)
+    applybranches(state,dm, max_order=max_order)
+end
+=#
diff --git a/src/experimental/NoisyCircuits.jl b/src/experimental/NoisyCircuits.jl
@@ -4,19 +4,13 @@ module NoisyCircuits
 
 using QuantumClifford
 using QuantumClifford: AbstractQCState, AbstractStabilizer, AbstractOperation, AbstractMeasurement, AbstractCliffordOperator, apply_single_x!, apply_single_y!, apply_single_z!, registered_statuses
-import QuantumClifford: applywstatus!, affectedqubits
+import QuantumClifford: applywstatus!, affectedqubits, applybranches, applynoise_branches
 
 using Combinatorics: combinations
 using Base.Cartesian
 
-export AbstractOperation,
-       UnbiasedUncorrelatedNoise,
-       NoisyBellMeasurement,
-       DecisionGate, ConditionalGate,
-       affectedqubits, applynoise!,# TODO rename applyop to apply_mc
-       applybranches, applynoise_branches,# TODO rename applybranches to apply_pe
-       petrajectory, petrajectories,
-       ConditionalGate, DecisionGate
+export NoisyBellMeasurement,
+       DecisionGate, ConditionalGate
 
 #TODO all these structs should use specified types
 #TODO all of the methods should better specified type signatures
@@ -53,68 +47,6 @@ function applywstatus!(s::AbstractQCState, m::NoisyBellMeasurement)
     end
 end
 
-"""Compute all possible new states after the application of the given operator. Reports the probability of each one of them. Deterministic, part of the Perturbative Expansion interface."""
-function applybranches end
-
-#TODO is the use of copy here necessary?
-function applybranches(state, op; max_order=1)
-    [(applywstatus!(copy(state),op)...,1,0)]
-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; max_order=1)
-    n = nqubits(s)
-    l = length(indices)
-    infid = noise.errprobthird
-    if l==0
-        return [s,one(infid)]
-    end
-    error1 = 3*infid
-    no_error1 = 1-error1
-    no_error = no_error1^l
-    results = [(copy(s),no_error,0)] # state, prob, order
-    for order in 1:min(max_order,l)
-        error_prob = no_error1^(l-order)*infid^order
-        for error_indices in combinations(indices, order)
-            _applynoise_branches_unbiased_uncorrelated(Val(order),s,error_indices,results,error_prob)
-        end
-    end
-    results
-end
-
-@generated function _applynoise_branches_unbiased_uncorrelated(::Val{order},s,error_indices,results,error_prob) where {order}
-    error_calls = Expr(:block)
-    for i in 1:order
-        call = quote (apply_single_x!,apply_single_y!,apply_single_z!)[$(Symbol(:i_,i))](new_state,error_indices[$i]) end
-        push!(error_calls.args, call)
-    end
-    # n nested loops, one for each affected qubit, each loop dedicated to the 3 possible errors (X, Y, or Z)
-    quote
-        @nloops $order i d->1:3 begin
-            new_state = copy(s)
-            $error_calls
-            push!(results,(new_state, error_prob, order))
-        end
-        results
-    end
-end
-
-function applybranches(s::AbstractQCState, nop::NoiseOpAll; max_order=1)
-    n = nqubits(s)
-    return [(state, continue_stat, prob, order) for (state, prob, order) in applynoise_branches(s, nop.noise, 1:n, max_order=max_order)]
-end
-
-function applybranches(s::AbstractQCState, nop::NoiseOp; max_order=1)
-    return [(state, continue_stat, prob, order) for (state, prob, order) in applynoise_branches(s, nop.noise, affectedqubits(nop), max_order=max_order)]
-end
-
-function applybranches(s::AbstractQCState, g::NoisyGate; max_order=1)
-    news, _,_,_ = applybranches(s,g.gate,max_order=max_order)[1] # TODO this assumes only one always successful branch for the gate
-    return [(state, continue_stat, prob, order) for (state, prob, order) in applynoise_branches(news, g.noise, affectedqubits(g), max_order=max_order)]
-end
-
 function applybranches(s::AbstractQCState, m::NoisyBellMeasurement; max_order=1)
     measurement_branches = applybranches(s, m.meas, max_order=max_order)
     if max_order==0
@@ -133,74 +65,6 @@ function applybranches(s::AbstractQCState, m::NoisyBellMeasurement; max_order=1)
     end
 end
 
-# TODO a lot of repetition with applywstatus!
-function applybranches(s::AbstractQCState, m::BellMeasurement; max_order=1)
-    n = nqubits(s)
-    [(ns,iseven(r>>1 + m.parity) ? continue_stat : failure_stat, p,0)
-     for (ns,r,p) in _applybranches_measurement([(s,0x0,1.0)],m.measurements,n)]
-end
-
-# TODO XXX THIS IS PARTICULARLY INEFFICIENT recurrent implementation
-function _applybranches_measurement(branches, measurements, n)
-    if length(measurements) == 0
-        return branches
-    end
-
-    new_branches = []
-    pauli = measurements[1]
-    otherpaulis = measurements[2:end]
-
-    for (s,r0,p) in branches
-        _,anticom,r = project!(quantumstate(s),pauli)
-        if isnothing(r) # TODO anticom could be zero if there was a rank change
-            s1 = s
-            s2 = copy(s)
-            r1 = phases(stabilizerview(s1))[anticom] = 0x00
-            r2 = phases(stabilizerview(s2))[anticom] = 0x02
-            push!(new_branches, (s1,r0+r1,p/2))
-            push!(new_branches, (s2,r0+r2,p/2))
-        else
-            push!(new_branches, (s,r0+r,p))
-        end
-    end
-
-    return _applybranches_measurement(new_branches, otherpaulis, n)
-end
-
-"""Run a perturbative expansion to a given order. Uses applybranches under the hood."""
-function petrajectory(state, circuit; branch_weight=1.0, current_order=0, max_order=1)
-    if size(circuit)[1] == 0
-        return fill(zero(branch_weight), length(registered_statuses)-1)
-    end
-    next_op = circuit[1]
-    rest_of_circuit = circuit[2:end]
-
-    status_probs = fill(zero(branch_weight), length(registered_statuses)-1)
-
-    for (i,(newstate, status, prob, order)) in enumerate(applybranches(state, next_op, max_order=max_order-current_order))
-        if status==continue_stat # TODO is the copy below necessary?
-            out_probs = petrajectory(copy(newstate), rest_of_circuit,
-                branch_weight=branch_weight*prob, current_order=current_order+order, max_order=max_order)
-            status_probs .+= out_probs
-        else
-            status_probs[status.status] += prob*branch_weight
-        end
-    end
-
-    return status_probs
-end
-
-"""Run a perturbative expansion to a given order. This is the main public fuction for the perturbative expansion approach."""
-function petrajectories(initialstate, circuit; branch_weight=1.0, max_order=1)
-    status_probs = petrajectory(initialstate, circuit; branch_weight=branch_weight, current_order=0, max_order=max_order)
-    Dict([CircuitStatus(i)=>status_probs[i] for i in eachindex(status_probs)])
-end
-
-function applynoise_branches(state::Register, noise, indices; max_order=1)
-    [(Register(newstate,copy(state.bits)), prob, order)
-     for (newstate, prob, order) in applynoise_branches(quantumstate(state), noise, indices; max_order=max_order)]
-end
-
 # TODO tests for this
 function QuantumClifford.apply!(state::Register, op::ConditionalGate)
     if state.bits[op.controlbit]
@@ -221,45 +85,11 @@ function QuantumClifford.apply!(state::Register, op::DecisionGate)
     state
 end
 
-
 applybranches(s::Register, op::ConditionalGate; max_order=1) = [(applywstatus!(copy(s),op)...,1,0)]
 applybranches(s::Register, op::DecisionGate; max_order=1) = [(applywstatus!(copy(s),op)...,1,0)]
 
-function applybranches(s::Register, op::PauliMeasurement; max_order=1)
-    stab = s.stab
-    stab, anticom, r = project!(stab, op.pauli)
-    new_branches = []
-    if isnothing(r)
-        s1 = s
-        phases(stabilizerview(s1.stab))[anticom] = 0x00
-        s1.bits[op.bit] = false
-        s2 = copy(s)
-        phases(stabilizerview(s2.stab))[anticom] = 0x02
-        s2.bits[op.bit] = true
-        push!(new_branches, (s1,continue_stat,1/2,0))
-        push!(new_branches, (s2,continue_stat,1/2,0))
-    else
-        s.bits[op.bit] = r==0x02
-        push!(new_branches, (s,continue_stat,1,0))
-    end
-    new_branches
-end
-
-#= TODO switch to sMX etc
-function applybranches(state::Register, op::SparseMeasurement; max_order=1)
-    n = nqubits(state.stab) # TODO implement actual sparse measurements
-    p = zero(typeof(op.pauli), n)
-    for (ii,i) in enumerate(op.indices)
-        p[i] = op.pauli[ii]
-    end
-    dm = PauliMeasurement(p,op.bit)
-    applybranches(state,dm, max_order=max_order)
-end
-=#
-
 affectedqubits(m::NoisyBellMeasurement) = affectedqubits(m.meas)
 affectedqubits(d::ConditionalGate) = union(affectedqubits(d.truegate), affectedqubits(d.falsegate))
 affectedqubits(d::DecisionGate) = [(union(affectedqubits.(d.gates))...)...]
-affectedqubits(v::VerifyOp) = v.indices
 
 end
diff --git a/src/fastmemlayout.jl b/src/fastmemlayout.jl
@@ -31,6 +31,3 @@ fastcolumn(s::MixedStabilizer) = MixedStabilizer(fastcolumn(s.tab), s.rank)
 
 fastrow(s::MixedDestabilizer) = MixedDestabilizer(fastrow(s.tab), s.rank)
 fastcolumn(s::MixedDestabilizer) = MixedDestabilizer(fastcolumn(s.tab), s.rank)
-
-fastrow(s::PauliFrame) = PauliFrame(fastrow(s.frame), s.measurements)
-fastcolumn(s::PauliFrame) = PauliFrame(fastcolumn(s.frame), s.measurements)