diff --git a/CHANGELOG.md b/CHANGELOG.md
index 68c18278b..2c3bebd2c 100644
--- 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
index 6a06a8f09..cfdc6a6f2 100644
--- 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
index 34d95f4d4..a58760216 100644
--- 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
index f27bc2469..97a7dd23f 100644
--- 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
index a021dff98..4182c0e6c 100644
--- 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
index 8dc61c64f..60cd766fe 100644
--- 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
index b9d5eb461..3a6f6e71b 100644
--- 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)
diff --git a/src/mctrajectory.jl b/src/mctrajectory.jl
index c514ebdd8..75343fc23 100644
--- a/src/mctrajectory.jl
+++ b/src/mctrajectory.jl
@@ -48,7 +48,7 @@ function mctrajectory!(state,circuit)
     return state, continue_stat
 end
 
-function countmap(samples::Vector{CircuitStatus}) # A simpler faster version of StatsBase.countmap
+function countmap(samples) # A simpler faster version of StatsBase.countmap that only works for CircuitStatus
     counts = zeros(length(registered_statuses))
     for s in samples
         counts[s.status+1] += 1
@@ -58,8 +58,16 @@ end
 
 """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
-    return counts
+See also: [`pftrajectories`](@ref), [`petrajectories`](@ref)"""
+mctrajectories(initialstate,circuit;trajectories=500,keepstates::Bool=false) = _mctrajectories(initialstate,circuit;trajectories,keepstates=Val(keepstates))
+
+function _mctrajectories(initialstate,circuit;trajectories=500,keepstates::Val{B}) where {B}
+    if B
+        counter = DefaultDict{Tuple{typeof(initialstate),CircuitStatus},Int}
+        counts = counter((mctrajectory!(copy(initialstate),circuit)[2] for i in 1:trajectories)) # TODO use threads or at least a generator
+        return counts
+    else
+        counts = countmap((mctrajectory!(copy(initialstate),circuit)[2] for i in 1:trajectories)) # TODO use threads or at least a generator
+        return counts
+    end
 end
diff --git a/src/misc_ops.jl b/src/misc_ops.jl
index fb59f678e..5361d21e7 100644
--- a/src/misc_ops.jl
+++ b/src/misc_ops.jl
@@ -17,6 +17,7 @@ function apply!(state::AbstractStabilizer, m::PauliMeasurement)
     state
 end
 
+
 """A Clifford gate, applying the given `cliff` operator to the qubits at the selected `indices`.
 
 `apply!(state, cliff, indices)` and `apply!(state, SparseGate(cliff, indices))` give the same result."""
@@ -68,6 +69,40 @@ function applywstatus!(s::AbstractQCState, m::BellMeasurement)
     end
 end
 
+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
+
+
 """A "probe" to verify that the state of the qubits corresponds to a desired `good_state`, e.g. at the end of the execution of a circuit."""
 struct VerifyOp <: AbstractOperation
     good_state::Stabilizer
@@ -89,3 +124,5 @@ function applywstatus!(s::AbstractQCState, v::VerifyOp) # XXX It assumes the oth
     end
     return s, true_success_stat
 end
+
+operatordeterminism(::Type{VerifyOp}) = DeterministicOperatorTrait()
diff --git a/src/noise.jl b/src/noise.jl
index 05b189d4c..9291bebe9 100644
--- a/src/noise.jl
+++ b/src/noise.jl
@@ -113,3 +113,58 @@ function apply!(r::Register, n::NoiseOp)
     apply!(quantumstate(r), n)
     return r
 end
+
+##
+# petrajectories
+##
+
+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
diff --git a/src/operator_traits.jl b/src/operator_traits.jl
new file mode 100644
index 000000000..38e5de20b
--- /dev/null
+++ b/src/operator_traits.jl
@@ -0,0 +1,7 @@
+abstract type OperatorDeterminismTrait end
+
+struct DeterministicOperatorTrait <: OperatorDeterminismTrait end
+struct NondeterministicOperatorTrait <: OperatorDeterminismTrait end
+
+operatordeterminism(::Type{<:AbstractCliffordOperator}) = DeterministicOperatorTrait()
+operatordeterminism(::Type{<:AbstractOperation}) = NondeterministicOperatorTrait()
diff --git a/src/pauli_frames.jl b/src/pauli_frames.jl
index e7bc17295..141dc4142 100644
--- a/src/pauli_frames.jl
+++ b/src/pauli_frames.jl
@@ -17,6 +17,9 @@ 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, :))
 
+fastrow(s::PauliFrame) = PauliFrame(fastrow(s.frame), s.measurements)
+fastcolumn(s::PauliFrame) = PauliFrame(fastcolumn(s.frame), s.measurements)
+
 """
 $(TYPEDSIGNATURES)
 
@@ -132,6 +135,8 @@ 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.
+
+See also: [`mctrajectories`](@ref), [`petrajectories`](@ref)
 """
 function pftrajectories(circuit;trajectories=5000,threads=true)
     _pftrajectories(circuit;trajectories,threads)
diff --git a/src/petrajectory.jl b/src/petrajectory.jl
new file mode 100644
index 000000000..deb8fb993
--- /dev/null
+++ b/src/petrajectory.jl
@@ -0,0 +1,79 @@
+"""Compute all possible new states after the application of the given operator. Reports the probability of each one of them. Deterministic (as it reports all branches of potentially random processes), part of the Perturbative Expansion interface."""
+function applybranches end
+
+"""Compute all possible new states after the application of the given noise model. Reports the probability of each one of them. Deterministic (as it reports all branches of potentially random processes), part of the Noise interface."""
+function applynoise_branches end
+
+#TODO is the use of copy here necessary?
+#TODO `copy` does not preserve the `fastcolumn` vs `fastrow` effect
+applybranches(state, op; max_order=1) = applybranches(operatordeterminism(typeof(op)), state, op; max_order=max_order)
+
+function applybranches(::DeterministicOperatorTrait, state, op; max_order=1)
+    [(applywstatus!(copy(state),op)...,1,0)]
+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.
+        If this is an operator type that you have defined yourself, please either implement an `applybranches` method for it.
+        If you know for a fact that this method is deterministic, then you can opt into the default by giving it a `DeterministicOperatorTrait` trait, e.g. `operatordeterminism(::Type{YourType}) = DeterministicOperatorTrait()`.
+        If you were not working with custom operators, please report this as a bug in QuantumClifford.jl -- we probably have not gotten around to expanding this functionality to that particular operator.
+    """))
+end
+
+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
+
+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)}
+    if size(circuit)[1] == 0
+        return A()
+    end
+    next_op = circuit[1]
+    rest_of_circuit = circuit[2:end]
+
+    dict = A()
+
+    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_dict = petrajectory_keep(copy(newstate), rest_of_circuit,
+                branch_weight=branch_weight*prob, current_order=current_order+order, max_order=max_order)
+            DataStructures.merge!(dict, out_dict)
+        else
+            DataStructures.inc!(dict, (newstate,status), prob*branch_weight)
+        end
+    end
+
+    return dict
+end
+
+"""Run a perturbative expansion to a given order. This is the main public fuction for the perturbative expansion approach.
+
+See also: [`pftrajectories`](@ref), [`mctrajectories`](@ref)"""
+function petrajectories(initialstate, circuit; branch_weight=1.0, max_order=1, keepstates::Bool=false)
+    if keepstates
+        return petrajectory_keep(initialstate, circuit; branch_weight=branch_weight, current_order=0, max_order=max_order)
+    else
+        status_probs = petrajectory(initialstate, circuit; branch_weight=branch_weight, current_order=0, max_order=max_order)
+        return Dict([CircuitStatus(i)=>status_probs[i] for i in eachindex(status_probs)])
+    end
+end
diff --git a/test/test_syndromemeas.jl b/test/test_syndromemeas.jl
index 279cd2b24..23526250d 100644
--- a/test/test_syndromemeas.jl
+++ b/test/test_syndromemeas.jl
@@ -1,17 +1,13 @@
 using QuantumClifford
+using QuantumClifford: AbstractOperation
 
-test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off-by-one errors in the bit encoding.
-
-using QuantumClifford.Experimental.NoisyCircuits
-import AbstractAlgebra
-
-@testset "Syndrome Measurements" begin
+@testset "Syndrome Measurements with mctrajectory!" begin # TODO this is a rather old test that is now done in a few other places, e.g. the ECC module -- double check and consider deleting
     codeˢᵗᵉᵃⁿᵉ = S"Z__Z_ZZ
-                    _Z_ZZ_Z
-                    __Z_ZZZ
-                    X__X_XX
-                    _X_XX_X
-                    __X_XXX";
+                   _Z_ZZ_Z
+                   __Z_ZZZ
+                   X__X_XX
+                   _X_XX_X
+                   __X_XXX";
     function get_check(code,index) # TODO a fairly ugly piece of code that should just be part of the library
         s,n = size(code)
         @assert index<=s