diff --git a/docs/src/noisycircuits_perturb.md b/docs/src/noisycircuits_perturb.md
index 54866ae2b..04f040d81 100644
--- a/docs/src/noisycircuits_perturb.md
+++ b/docs/src/noisycircuits_perturb.md
@@ -42,7 +42,52 @@ circuit = [n,g1,g2,m,v]
petrajectories(initial_state, circuit)
```
-
+```@raw html
+
+
+flowchart TD
+ state[XX ZZ] -->|Simulation 2|gate1a([CNOT 1,3])
+ state[XX ZZ] -->|Simulation 3|gate1b([CNOT 1,3])
+ state[XX ZZ] -->|Simulation 1|gate1c([CNOT 1,3])
+ gate1a -->|Success| gate2a([CNOT 2,4])
+ gate1b -->|Success| gate2b([CNOT 2,4])
+ gate1c -.->|Error| gate2c([CNOT 2,4])
+ gate2a -->|Success|result2[Result]
+ gate2b -.->|Error|result3[Result]
+ gate2c -->|Success|result1[Result]
+ result1-->sum([Sum results and symbolic probabilities])
+ result2-->sum
+ result3-->sum
+ sum-->values[Probabilities
+ True Success = 0.903546
+ Failure = 0.0365069
+ False Success = 0.0547604]
+ style state fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style gate1a fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style gate1b fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style gate1c fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style gate2a fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style gate2b fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style gate2c fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style result1 fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style result2 fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style result3 fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style sum fill:#ffffff, font-family: 'Serif',stroke:#000000
+ style values fill:#ffffff, font-family: 'Serif',stroke:#000000
+ linkStyle 0 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 1 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 2 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 3 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 4 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 5 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 6 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 7 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 8 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 9 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 10 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 11 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+ linkStyle 12 stroke:#000000,stroke-width:1px,color:black,font-family: 'Serif'
+
```
For more examples, see the [notebook comparing the Monte Carlo and Perturbative method](https://nbviewer.jupyter.org/github/QuantumSavory/QuantumClifford.jl/blob/master/docs/src/notebooks/Perturbative_Expansions_vs_Monte_Carlo_Simulations.ipynb) or this tutorial on [entanglement purification](https://github.com/QuantumSavory/QuantumClifford.jl/blob/master/docs/src/notebooks/Noisy_Circuits_Tutorial_with_Purification_Circuits.ipynb).
## Symbolic expansions
diff --git a/src/petrajectory.jl b/src/petrajectory.jl
index 125ecc30e..8649d0843 100644
--- a/src/petrajectory.jl
+++ b/src/petrajectory.jl
@@ -12,6 +12,29 @@ function applybranches(::DeterministicOperatorTrait, state, op; max_order=1)
[(applywstatus!(copy(state),op)...,1,0)]
end
+function applybranches(::NondeterministicOperatorTrait, state, op::T; max_order=1) where {T<:Union{sMZ,sMX,sMY}}
+ s = copy(state)
+
+ if T===sMZ
+ d, anticom, res = projectZ!(s, op.qubit)
+ elseif T===sMX
+ d, anticom, res = projectX!(s, op.qubit)
+ elseif T===sMY
+ d, anticom, res = projectY!(s, op.qubit)
+ else
+ error("not reachable")
+ end
+
+ if isnothing(res)
+ tab(stabilizerview(s)).phases[anticom] = 0x0
+ s1 = copy(s)
+ tab(stabilizerview(s)).phases[anticom] = 0x2
+ s2 = s
+ return [(s1, continue_stat, .5, 0), (s2, continue_stat, .5, 0)]
+ else
+ return [(s, continue_stat, 1, 0)] end
+end
+
function applybranches(::NondeterministicOperatorTrait, state, op; max_order=1)
throw(ArgumentError(lazy"""
You are trying to apply a non-deterministic operator $(typeof(op)) in a perturbative expansion, but this particular operator does not have a `applybranches` method defined for it.
@@ -45,9 +68,12 @@ end
function petrajectory_keep(state, circuit; branch_weight=1.0, current_order=0, max_order=1) # TODO a lot of repetition with petrajectory - dry out
A = Accumulator{Tuple{typeof(state),CircuitStatus},typeof(branch_weight)}
+ dict = A()
if size(circuit)[1] == 0
- return A()
+ DataStructures.inc!(dict, (state,continue_stat), branch_weight)
+ return dict
end
+
next_op = circuit[1]
rest_of_circuit = circuit[2:end]
@@ -67,7 +93,8 @@ function petrajectory_keep(state, circuit; branch_weight=1.0, current_order=0, m
end
"""Run a perturbative expansion to a given order. This is the main public function for the perturbative expansion approach.
-
+Currently only has defined behavior for Clifford operators(AbstractCliffordOperator) and single-qubit X, Y and Z measurements(sMX, sMY, sMZ)
+If using the measurements, set keepstates to true. When keepstates is false, will return 0 for true success probability.
See also: [`pftrajectories`](@ref), [`mctrajectories`](@ref)"""
function petrajectories(initialstate, circuit; branch_weight=1.0, max_order=1, keepstates::Bool=false)
if keepstates
diff --git a/test/test_noisycircuits.jl b/test/test_noisycircuits.jl
index 4f5e92c3e..c203f6553 100644
--- a/test/test_noisycircuits.jl
+++ b/test/test_noisycircuits.jl
@@ -233,6 +233,34 @@
@test stabilizerview(state.stab) == S"ZZ
ZI"
end
+
+ @testset "Symbolic Measurements" begin
+ state = MixedDestabilizer(S"Z")
+ had = sHadamard(1)
+ Zmeas = sMZ(1)
+ Xmeas = sMX(1)
+ Ymeas = sMY(1)
+ circuit = [Zmeas]
+ pe = petrajectories(state, circuit, keepstates=true)
+ for i in pe
+ @test i[2] == 1.0
+ end
+ circuit = [had, Zmeas]
+ pe = petrajectories(state, circuit, keepstates=true)
+ for i in pe
+ @test i[2] == 0.5
+ end
+ circuit = [had, Xmeas]
+ pe = petrajectories(state, circuit, keepstates=true)
+ for i in pe
+ @test i[2] == 1.0
+ end
+ circuit = [had, Ymeas]
+ pe = petrajectories(state, circuit, keepstates=true)
+ for i in pe
+ @test i[2] == 0.5
+ end
+ end
end
@testset "Classical Bits" begin