From b7a311ef66a08e56dfb4a538a2b6cb47d8ea02e8 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Thu, 8 Jun 2023 23:58:05 -0400
Subject: [PATCH 1/9] initial implementation for nonclifford gates

---
 Project.toml             |  1 +
 src/QuantumClifford.jl   |  5 ++-
 src/nonclifford.jl       | 97 ++++++++++++++++++++++++++++++++++++++++
 test/test_nonclifford.jl | 37 +++++++++++++++
 4 files changed, 139 insertions(+), 1 deletion(-)
 create mode 100644 src/nonclifford.jl
 create mode 100644 test/test_nonclifford.jl

diff --git a/Project.toml b/Project.toml
index f084c5b4c..942f7ccb3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.8.2"
 
 [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"
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 524b10c40..ffbe8766e 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -690,7 +690,9 @@ function MixedDestabilizer(d::Destabilizer)
     end
 end
 
-function MixedDestabilizer(d::MixedStabilizer) MixedDestabilizer(stabilizerview(d)) end
+MixedDestabilizer(d::MixedStabilizer) = MixedDestabilizer(stabilizerview(d))
+
+MixedDestabilizer(s::MixedDestabilizer) = s
 
 Base.length(d::MixedDestabilizer) = length(d.tab)÷2
 
@@ -1377,6 +1379,7 @@ include("tableau_show.jl")
 include("sumtypes.jl")
 include("precompiles.jl")
 include("ecc/ECC.jl")
+include("nonclifford.jl")
 include("plotting_extensions.jl")
 
 end #module
diff --git a/src/nonclifford.jl b/src/nonclifford.jl
new file mode 100644
index 000000000..4cbeabbad
--- /dev/null
+++ b/src/nonclifford.jl
@@ -0,0 +1,97 @@
+using LinearAlgebra: dot
+using DataStructures: DefaultDict
+
+mutable struct StabMixture{T,F}
+    stab::T
+    destabweights::DefaultDict{Tuple{BitVector, BitVector}, F, F}
+end
+
+function StabMixture(state)
+    n = nqubits(state)
+    StabMixture(MixedDestabilizer(state), DefaultDict(0.0im, (falses(n),falses(n))=>1.0+0.0im)) # TODO maybe it should default to Destabilizer, not MixedDestabilizer
+end
+
+StabMixture(s::StabMixture) = s
+
+function Base.show(io::IO, s::StabMixture)
+    println(io, "A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is")
+    show(io,s.stab)
+    println(io)
+    println(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
+    for ((di,dj), χ) in s.destabweights
+        println(io, "  ", χ, " | ", string(_stabmixdestab(s.stab, di)), " | ", string(_stabmixdestab(s.stab, dj)))
+    end
+end
+
+function _stabmixdestab(mixeddestab, d)
+    destab = destabilizerview(mixeddestab)
+    p = zero(PauliOperator, nqubits(mixeddestab))
+    for i in eachindex(d)
+        if d[i]
+            mul_left!(p, destab, i) # TODO check whether you are screwing up the ordering
+        end
+    end
+    p
+end
+
+function apply!(state::StabMixture, gate::AbstractCliffordOperator) # TODO conjugate also the destabs
+    apply!(state.stab, gate)
+    state
+end
+
+abstract type AbstractPauliChannel <: AbstractOperation end
+
+struct TGate <: AbstractPauliChannel
+    qubit::Int
+end
+
+struct PauliChannel{T,S} <: AbstractPauliChannel
+    paulis::T
+    weights::S
+end
+
+function apply!(state::StabMixture, gate::PauliChannel)
+    dict = state.destabweights
+    stab = state.stab
+    newdict = typeof(dict)(zero(eltype(dict).parameters[2])) # TODO jeez, this is ugly
+    for ((dᵢ,dⱼ), χ) in dict
+        for ((Pₗ,Pᵣ), w) in zip(gate.paulis,gate.weights)
+            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ,stab) # TODO this should be cached?
+            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ,stab) # TODO this should be cached?
+            c = (dot(dₗˢᵗᵃᵇ,dᵢ) + dot(dᵣˢᵗᵃᵇ,dⱼ))*2
+            dᵢ′ = dₗ .⊻ dᵢ
+            dⱼ′ = dᵣ .⊻ dⱼ
+            χ′ = χ * w * (-1)^c * (1im)^(phaseₗ-phaseᵣ)
+            newdict[(dᵢ′,dⱼ′)] += χ′
+        end
+    end
+    for (k,v) in newdict # TODO is it safe to modify a dict while iterating over it?
+        if abs(v) < 1e-14 # TODO parameterize this pruning parameter
+            delete!(newdict, k)
+        end
+    end
+    state.destabweights = newdict
+    state
+end
+
+function rowdecompose(pauli,state::Union{MixedDestabilizer, Destabilizer})
+    n = nqubits(pauli)
+    stab = stabilizerview(state)
+    dest = destabilizerview(state)
+    b = falses(n)
+    c = falses(n)
+    Pₜ = zero(PauliOperator, n) # TODO is this the best API choice!?
+    for i in 1:n
+        if comm(pauli, stab, i) != 0
+            b[i] = true
+            mul_left!(Pₜ, dest, i)
+        end
+    end
+    for i in 1:n
+        if comm(pauli, dest, i) != 0
+            c[i] = true
+            mul_left!(Pₜ, stab, i) # ?
+        end
+    end
+    return Pₜ.phase[], b, c
+end
diff --git a/test/test_nonclifford.jl b/test/test_nonclifford.jl
new file mode 100644
index 000000000..b1cc930fc
--- /dev/null
+++ b/test/test_nonclifford.jl
@@ -0,0 +1,37 @@
+using QuantumClifford
+using QuantumClifford: StabMixture, rowdecompose, PauliChannel
+using Test
+
+##
+
+for n in [1,2,63,64,65,66]
+    n = 6
+    p = random_pauli(n; nophase=true)
+    s = random_destabilizer(n)
+    phase, b, c = rowdecompose(p,s)
+    p0 = zero(p)
+    for (i,f) in pairs(b)
+        f && mul_left!(p0, destabilizerview(s), i)
+    end
+    for (i,f) in pairs(c)
+        f && mul_left!(p0, stabilizerview(s), i)
+    end
+    @test (-im)^phase*p0 == p
+end
+
+##
+
+tgate = PauliChannel(
+[(I, I), (I, Z), (Z, I), (Z, Z)],
+[cos(π/8)^2, -im*sin(π/8)*cos(π/8),  im*sin(π/8)*cos(π/8), sin(π/8)^2]
+)
+
+state = StabMixture(S"X")
+
+apply!(state, tgate)
+apply!(state, tgate)
+apply!(state, tgate)
+apply!(state, tgate)
+
+@test state.destabweights |> values |> collect == [1]
+@test state.destabweights |> keys |> collect == [([1],[1])]

From 7a1b728b1835113c2de691e5f2a64587db445343 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Fri, 9 Jun 2023 00:23:04 -0400
Subject: [PATCH 2/9] missing compat

---
 Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Project.toml b/Project.toml
index 942f7ccb3..57a1d1e7d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,6 +33,7 @@ QuantumCliffordQuantikzExt = "Quantikz"
 
 [compat]
 Combinatorics = "1.0"
+DataStructures = "0.18"
 DocStringExtensions = "0.8, 0.9"
 Graphs = "1.4.1"
 HostCPUFeatures = "0.1.6"

From e7ecbfbdbe693cf3769daa432bbddc66ded66ebb Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Sun, 11 Jun 2023 15:17:43 -0400
Subject: [PATCH 3/9] MixedDestabilizer(::StabMixture)

---
 src/nonclifford.jl | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index 4cbeabbad..56f8122a7 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -13,6 +13,20 @@ end
 
 StabMixture(s::StabMixture) = s
 
+function MixedDestabilizer(s::StabMixture)
+    if length(s.destabweights) != 1
+        throw(DomainError("Trying to convert a non-Clifford state (instance of type StabMixture with more than one term in its sum representation) to a Clifford state (of type MixedDestabilizer). This is not possible. Consider whether you have performed some (unforeseen) non-Clifford operation on the state."))
+    else
+        ((dᵢ, dⱼ), χ) = first(s.destabweights) # TODO phases/order seem wrong here
+        dᵢ = _stabmixdestab(s.stab, dᵢ)
+        dⱼ = _stabmixdestab(s.stab, dⱼ)
+        stab = copy(s.stab)
+        apply!(stab, dᵢ)
+        apply!(stab, dⱼ)
+        return stab
+    end
+end
+
 function Base.show(io::IO, s::StabMixture)
     println(io, "A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is")
     show(io,s.stab)

From 2bd7325867f9c9d043de3d31495b9cb4b85fef79 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Mon, 12 Jun 2023 01:24:07 -0400
Subject: [PATCH 4/9] some more mul_left / mul_right tooling

---
 src/mul_leftright.jl       | 26 ++++++++++++++++++++++++++
 src/nonclifford.jl         | 16 ++++++++--------
 test/test_mul_leftright.jl |  2 ++
 3 files changed, 36 insertions(+), 8 deletions(-)

diff --git a/src/mul_leftright.jl b/src/mul_leftright.jl
index f21ab427f..5be86c0fd 100644
--- a/src/mul_leftright.jl
+++ b/src/mul_leftright.jl
@@ -144,3 +144,29 @@ end
     n = nqubits(t)
     mul_left!(t, i+n, j+n; phases=phases)
 end
+
+@inline function mul_left!(s::Tableau, p::PauliOperator; phases::Val{B}=Val(true)) where B # TODO multithread
+    @inbounds @simd for m in eachindex(s)
+        extra_phase = mul_left!((@view s.xzs[:,m]), p.xz; phases=phases)
+        B && (s.phases[m] = (extra_phase+s.phases[m]+p.phase[])&0x3)
+    end
+    s
+end
+
+@inline function mul_left!(s::AbstractStabilizer, p::PauliOperator; phases::Val{B}=Val(true)) where B
+    mul_left!(tab(s), p; phases=phases)
+    s
+end
+
+@inline function mul_right!(s::Tableau, p::PauliOperator; phases::Val{B}=Val(true)) where B # TODO multithread
+    @inbounds @simd for m in eachindex(s)
+        extra_phase = mul_right!((@view s.xzs[:,m]), p.xz; phases=phases)
+        B && (s.phases[m] = (extra_phase+s.phases[m]+p.phase[])&0x3)
+    end
+    s
+end
+
+@inline function mul_right!(s::AbstractStabilizer, p::PauliOperator; phases::Val{B}=Val(true)) where B
+    mul_right!(tab(s), p; phases=phases)
+    s
+end
diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index 56f8122a7..a77542420 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -17,12 +17,12 @@ function MixedDestabilizer(s::StabMixture)
     if length(s.destabweights) != 1
         throw(DomainError("Trying to convert a non-Clifford state (instance of type StabMixture with more than one term in its sum representation) to a Clifford state (of type MixedDestabilizer). This is not possible. Consider whether you have performed some (unforeseen) non-Clifford operation on the state."))
     else
-        ((dᵢ, dⱼ), χ) = first(s.destabweights) # TODO phases/order seem wrong here
-        dᵢ = _stabmixdestab(s.stab, dᵢ)
-        dⱼ = _stabmixdestab(s.stab, dⱼ)
-        stab = copy(s.stab)
-        apply!(stab, dᵢ)
-        apply!(stab, dⱼ)
+        ((dᵢ, dⱼ), χ) = first(s.destabweights)
+        dᵢ = _stabmixdestab(s.stab, dᵢ) # TODO
+        dⱼ = _stabmixdestab(s.stab, dⱼ) # make
+        stab = copy(s.stab)             # a faster
+        mul_left!(stab, dᵢ)             # in-place
+        mul_right!(stab, dⱼ)            # version
         return stab
     end
 end
@@ -70,8 +70,8 @@ function apply!(state::StabMixture, gate::PauliChannel)
     newdict = typeof(dict)(zero(eltype(dict).parameters[2])) # TODO jeez, this is ugly
     for ((dᵢ,dⱼ), χ) in dict
         for ((Pₗ,Pᵣ), w) in zip(gate.paulis,gate.weights)
-            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ,stab) # TODO this should be cached?
-            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ,stab) # TODO this should be cached?
+            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ,stab)
+            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ,stab)
             c = (dot(dₗˢᵗᵃᵇ,dᵢ) + dot(dᵣˢᵗᵃᵇ,dⱼ))*2
             dᵢ′ = dₗ .⊻ dᵢ
             dⱼ′ = dᵣ .⊻ dⱼ
diff --git a/test/test_mul_leftright.jl b/test/test_mul_leftright.jl
index 05f6ca512..ecf2b36eb 100644
--- a/test/test_mul_leftright.jl
+++ b/test/test_mul_leftright.jl
@@ -17,6 +17,8 @@ test_sizes = [1,2,10,63,64,65,127,128,129] # Including sizes that would test off
             @test mul_left!(copy(p2), p1) == (-1)^comm(p1,p2) * mul_right!(copy(p2), p1)
             @test mul_left!(copy(p2), s[i]) == mul_left!(copy(p2), s, i) == s[i]*p2
             @test mul_right!(copy(p2), s[i]) == mul_right!(copy(p2), s, i) == p2*s[i]
+            @test mul_left!(copy(s), p2)[i] == p2*s[i]
+            @test mul_right!(copy(s), p2)[i] == s[i]*p2
         end
     end
 end

From 464978f4e47bf239b83824b92358c02a30160b6f Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Mon, 12 Jun 2023 13:19:10 -0400
Subject: [PATCH 5/9] minor test cleanup

---
 test/test_nonclifford.jl | 52 ++++++++++++++++++++++------------------
 1 file changed, 29 insertions(+), 23 deletions(-)

diff --git a/test/test_nonclifford.jl b/test/test_nonclifford.jl
index b1cc930fc..baebf77bc 100644
--- a/test/test_nonclifford.jl
+++ b/test/test_nonclifford.jl
@@ -1,37 +1,43 @@
 using QuantumClifford
 using QuantumClifford: StabMixture, rowdecompose, PauliChannel
 using Test
+using InteractiveUtils
+using Random
 
 ##
 
-for n in [1,2,63,64,65,66]
-    n = 6
-    p = random_pauli(n; nophase=true)
-    s = random_destabilizer(n)
-    phase, b, c = rowdecompose(p,s)
-    p0 = zero(p)
-    for (i,f) in pairs(b)
-        f && mul_left!(p0, destabilizerview(s), i)
+@testset "Pauli decomposition into destabilizers" begin
+    for n in [1,2,63,64,65,66]
+        n = 6
+        p = random_pauli(n; nophase=true)
+        s = random_destabilizer(n)
+        phase, b, c = rowdecompose(p,s)
+        p0 = zero(p)
+        for (i,f) in pairs(b)
+            f && mul_left!(p0, destabilizerview(s), i)
+        end
+        for (i,f) in pairs(c)
+            f && mul_left!(p0, stabilizerview(s), i)
+        end
+        @test (-im)^phase*p0 == p
     end
-    for (i,f) in pairs(c)
-        f && mul_left!(p0, stabilizerview(s), i)
-    end
-    @test (-im)^phase*p0 == p
 end
 
 ##
 
-tgate = PauliChannel(
-[(I, I), (I, Z), (Z, I), (Z, Z)],
-[cos(π/8)^2, -im*sin(π/8)*cos(π/8),  im*sin(π/8)*cos(π/8), sin(π/8)^2]
-)
+@testset "PauliChannel T gate ^4 = Id" begin
+    tgate = PauliChannel(
+    [(I, I), (I, Z), (Z, I), (Z, Z)],
+    [cos(π/8)^2, -im*sin(π/8)*cos(π/8),  im*sin(π/8)*cos(π/8), sin(π/8)^2]
+    )
 
-state = StabMixture(S"X")
+    state = StabMixture(S"-Z")
 
-apply!(state, tgate)
-apply!(state, tgate)
-apply!(state, tgate)
-apply!(state, tgate)
+    apply!(state, tgate)
+    apply!(state, tgate)
+    apply!(state, tgate)
+    apply!(state, tgate)
 
-@test state.destabweights |> values |> collect == [1]
-@test state.destabweights |> keys |> collect == [([1],[1])]
+    @test state.destabweights |> values |> collect == [1]
+    @test state.destabweights |> keys |> collect == [([1],[1])]
+end

From 106a32eac2328f9319ce21830775e3930af71227 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Mon, 10 Jul 2023 12:15:04 -0400
Subject: [PATCH 6/9] cleanup

---
 src/QuantumClifford.jl   |  2 +-
 src/nonclifford.jl       | 58 +++++++++++++++++++++++++++++++++-------
 test/test_nonclifford.jl | 14 +++++-----
 3 files changed, 56 insertions(+), 18 deletions(-)

diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index c6b041d7c..3a3477272 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -8,7 +8,7 @@ module QuantumClifford
 # TODO Significant performance improvements: many operations do not need phase=true if the Pauli operations commute
 
 import LinearAlgebra
-using LinearAlgebra: inv, mul!, rank, Adjoint
+using LinearAlgebra: inv, mul!, rank, Adjoint, dot
 import DataStructures
 using DataStructures: DefaultDict, Accumulator
 using Combinatorics: combinations
diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index a77542420..4d6c7b831 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -1,11 +1,28 @@
-using LinearAlgebra: dot
-using DataStructures: DefaultDict
+#=
+1. adding tests for basic correctness
+2. single qubit gates / channels (and tests)
+3. embedding single qubit gates - Stefan
+4. pretty printing - Stefan
+5. good docstrings
+6. some superficial documentation
+7. picking names
+8. conversion into density matrices (QuantumOptics.jl) - Stefan
+9. special small gates
+10. make an overleaf for a paper
+=#
+
 
 mutable struct StabMixture{T,F}
     stab::T
     destabweights::DefaultDict{Tuple{BitVector, BitVector}, F, F}
 end
 
+# TODO: figure out a way to compute whether a StabMixture instance is pure or mixed
+# TODO which hane of these (and others)
+# GeneralizedStabilizer
+# StabilizerFrame
+# GenStabFrame
+
 function StabMixture(state)
     n = nqubits(state)
     StabMixture(MixedDestabilizer(state), DefaultDict(0.0im, (falses(n),falses(n))=>1.0+0.0im)) # TODO maybe it should default to Destabilizer, not MixedDestabilizer
@@ -42,7 +59,7 @@ function _stabmixdestab(mixeddestab, d)
     p = zero(PauliOperator, nqubits(mixeddestab))
     for i in eachindex(d)
         if d[i]
-            mul_left!(p, destab, i) # TODO check whether you are screwing up the ordering
+            mul_right!(p, destab, i) # TODO check whether you are screwing up the ordering
         end
     end
     p
@@ -67,15 +84,17 @@ end
 function apply!(state::StabMixture, gate::PauliChannel)
     dict = state.destabweights
     stab = state.stab
-    newdict = typeof(dict)(zero(eltype(dict).parameters[2])) # TODO jeez, this is ugly
-    for ((dᵢ,dⱼ), χ) in dict
-        for ((Pₗ,Pᵣ), w) in zip(gate.paulis,gate.weights)
+    tzero = zero(eltype(dict).parameters[2])
+    tone = one(eltype(dict).parameters[2])
+    newdict = typeof(dict)(tzero) # TODO jeez, this is ugly
+    for ((dᵢ,dⱼ), χ) in dict # the state
+        for ((Pₗ,Pᵣ), w) in zip(gate.paulis,gate.weights) # the channel
             phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ,stab)
             phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ,stab)
             c = (dot(dₗˢᵗᵃᵇ,dᵢ) + dot(dᵣˢᵗᵃᵇ,dⱼ))*2
             dᵢ′ = dₗ .⊻ dᵢ
             dⱼ′ = dᵣ .⊻ dⱼ
-            χ′ = χ * w * (-1)^c * (1im)^(phaseₗ-phaseᵣ)
+            χ′ = χ * w * (-tone)^c * (tone*im)^(-phaseₗ+phaseᵣ+4)
             newdict[(dᵢ′,dⱼ′)] += χ′
         end
     end
@@ -88,6 +107,24 @@ function apply!(state::StabMixture, gate::PauliChannel)
     state
 end
 
+"""Decompose a Pauli ``P`` in terms of stabilizer and destabilizer rows from a given tableaux.
+
+For given tableaux of rows destabilizer rows ``\\{d_i\\}`` and stabilizer rows ``\\{s_i\\}``,
+there are boolean vectors ``b`` and ``c`` such that
+``P = i^p \\prod_i d_i^{b_i} \\prod_i s_i^{c_i}``.
+
+This function returns `p`, `b`, `c`.
+
+```
+julia> s = MixedDestabilizer(ghz(2));
+
+julia> rowdecompose(P"XY", s)
+(3, Bool[1, 0], Bool[1, 1])
+
+julia> -im * P"Z_" * P"XX" * P"ZZ"
++ XY
+```
+"""
 function rowdecompose(pauli,state::Union{MixedDestabilizer, Destabilizer})
     n = nqubits(pauli)
     stab = stabilizerview(state)
@@ -98,14 +135,15 @@ function rowdecompose(pauli,state::Union{MixedDestabilizer, Destabilizer})
     for i in 1:n
         if comm(pauli, stab, i) != 0
             b[i] = true
-            mul_left!(Pₜ, dest, i)
+            mul_right!(Pₜ, dest, i)
         end
     end
     for i in 1:n
         if comm(pauli, dest, i) != 0
             c[i] = true
-            mul_left!(Pₜ, stab, i) # ?
+            mul_right!(Pₜ, stab, i)
         end
     end
-    return Pₜ.phase[], b, c
+    p = mod(-Pₜ.phase[],4) # complex conjugate
+    return p, b, c
 end
diff --git a/test/test_nonclifford.jl b/test/test_nonclifford.jl
index baebf77bc..ba802b98c 100644
--- a/test/test_nonclifford.jl
+++ b/test/test_nonclifford.jl
@@ -1,25 +1,25 @@
 using QuantumClifford
-using QuantumClifford: StabMixture, rowdecompose, PauliChannel
+using QuantumClifford: StabMixture, rowdecompose, PauliChannel, mul_left!, mul_right!
 using Test
 using InteractiveUtils
 using Random
 
 ##
 
+
 @testset "Pauli decomposition into destabilizers" begin
-    for n in [1,2,63,64,65,66]
-        n = 6
+    for n in [1,2,63,64,65,66,300]
         p = random_pauli(n; nophase=true)
         s = random_destabilizer(n)
         phase, b, c = rowdecompose(p,s)
         p0 = zero(p)
         for (i,f) in pairs(b)
-            f && mul_left!(p0, destabilizerview(s), i)
+            f && mul_right!(p0, destabilizerview(s), i)
         end
         for (i,f) in pairs(c)
-            f && mul_left!(p0, stabilizerview(s), i)
+            f && mul_right!(p0, stabilizerview(s), i)
         end
-        @test (-im)^phase*p0 == p
+        @test (im)^phase*p0 == p
     end
 end
 
@@ -31,7 +31,7 @@ end
     [cos(π/8)^2, -im*sin(π/8)*cos(π/8),  im*sin(π/8)*cos(π/8), sin(π/8)^2]
     )
 
-    state = StabMixture(S"-Z")
+    state = StabMixture(S"X")
 
     apply!(state, tgate)
     apply!(state, tgate)

From fc6ef91aa81ee2893638e8cb58619acf694af451 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Tue, 11 Jul 2023 15:44:44 -0400
Subject: [PATCH 7/9] fixup

---
 CHANGELOG.md             |   4 +
 src/QuantumClifford.jl   | 137 ++----------------------------
 src/nonclifford.jl       |  98 +++++++++++++++++++---
 src/pauli_operator.jl    | 176 +++++++++++++++++++++++++++++++++++++++
 test/runtests.jl         |   1 +
 test/test_embed.jl       |  10 +++
 test/test_nonclifford.jl |   8 +-
 7 files changed, 286 insertions(+), 148 deletions(-)
 create mode 100644 src/pauli_operator.jl
 create mode 100644 test/test_embed.jl

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2c3bebd2c..4c861ae80 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,10 @@
 
 # News
 
+## v0.8.12 - dev
+
+- Initial implementation of non-Clifford simulation (mainly for circuits that are slightly non-Clifford, e.g. containing T gates). See `StabMixture`, `PauliChannel`, and `tT`.
+
 ## 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.
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 3a3477272..6b201ffc3 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -15,7 +15,7 @@ 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
+import QuantumInterface: tensor, ⊗, tensor_pow, apply!, nqubits, expect, project!, reset_qubits!, traceout!, ptrace, apply!, projectX!, projectY!, projectZ!, entanglement_entropy, embed
 
 export
     @P_str, PauliOperator, ⊗, I, X, Y, Z,
@@ -41,7 +41,7 @@ export
     generate!, project!, reset_qubits!, traceout!,
     projectX!, projectY!, projectZ!,
     projectrand!, projectXrand!, projectYrand!, projectZrand!,
-    puttableau!,
+    puttableau!, embed,
     # Clifford Ops
     CliffordOperator, @C_str, permute,
     tCNOT, tCPHASE, tSWAP, tHadamard, tPhase, tId1,
@@ -76,7 +76,10 @@ export
     # mctrajectories
     CircuitStatus, continue_stat, true_success_stat, false_success_stat, failure_stat,
     mctrajectory!, mctrajectories, applywstatus!,
+    # petrajectories
     petrajectories, applybranches,
+    # nonclifford
+    StabMixture, PauliChannel, tT,
     # makie plotting -- defined only when extension is loaded
     stabilizerplot, stabilizerplot_axis,
     # sum types
@@ -114,135 +117,7 @@ include("macrotools.jl")
 abstract type AbstractOperation end
 abstract type AbstractCliffordOperator <: AbstractOperation end
 
-"""
-A multi-qubit Pauli operator (``±\\{1,i\\}\\{I,Z,X,Y\\}^{\\otimes n}``).
-
-A Pauli can be constructed with the `P` custom string macro or by building
-up one through products and tensor products of smaller operators.
-
-```jldoctest
-julia> pauli3 = P"-iXYZ"
--iXYZ
-
-julia> pauli4 = 1im * pauli3 ⊗ X
-+ XYZX
-
-julia> Z*X
-+iY
-```
-
-We use a typical F(2,2) encoding internally. The X and Z bits are stored
-in a single concatenated padded array of UInt chunks of a bit array.
-
-```jldoctest
-julia> p = P"-IZXY";
-
-
-julia> p.xz
-2-element Vector{UInt64}:
- 0x000000000000000c
- 0x000000000000000a
-```
-
-You can access the X and Z bits through getters and setters or through the
-`xview`, `zview`, `xbit`, and `zbit` functions.
-
-```jldoctest
-julia> p = P"XYZ"; p[1]
-(true, false)
-
-julia> p[1] = (true, true); p
-+ YYZ
-```
-"""
-struct PauliOperator{Tz<:AbstractArray{UInt8,0}, Tv<:AbstractVector{<:Unsigned}} <: AbstractCliffordOperator
-    phase::Tz
-    nqubits::Int
-    xz::Tv
-end
-
-PauliOperator(phase::UInt8, nqubits::Int, xz::Tv) where Tv<:AbstractVector{<:Unsigned} = PauliOperator(fill(UInt8(phase),()), nqubits, xz)
-PauliOperator(phase::UInt8, x::AbstractVector{Bool}, z::AbstractVector{Bool}) = PauliOperator(fill(UInt8(phase),()), length(x), vcat(reinterpret(UInt,BitVector(x).chunks),reinterpret(UInt,BitVector(z).chunks)))
-PauliOperator(x::AbstractVector{Bool}, z::AbstractVector{Bool}) = PauliOperator(0x0, x, z)
-PauliOperator(xz::AbstractVector{Bool}) = PauliOperator(0x0, (@view xz[1:end÷2]), (@view xz[end÷2+1:end]))
-
-"""Get a view of the X part of the `UInt` array of packed qubits of a given Pauli operator."""
-function xview(p::PauliOperator)
-    @view p.xz[1:end÷2]
-end
-"""Get a view of the Y part of the `UInt` array of packed qubits of a given Pauli operator."""
-function zview(p::PauliOperator)
-    @view p.xz[end÷2+1:end]
-end
-"""Extract as a new bit array the X part of the `UInt` array of packed qubits of a given Pauli operator."""
-function xbit(p::PauliOperator)
-    one = eltype(p.xz)(1)
-    size = sizeof(eltype(p.xz))*8
-    [(word>>s)&one==one for word in xview(p) for s in 0:size-1][begin:p.nqubits]
-end
-"""Extract as a new bit array the Z part of the `UInt` array of packed qubits of a given Pauli operator."""
-function zbit(p::PauliOperator)
-    one = eltype(p.xz)(1)
-    size = sizeof(eltype(p.xz))*8
-    [(word>>s)&one==one for word in zview(p) for s in 0:size-1][begin:p.nqubits]
-end
-
-function _P_str(a)
-    letters = filter(x->occursin(x,"_IZXY"),a)
-    phase = phasedict[strip(filter(x->!occursin(x,"_IZXY"),a))]
-    PauliOperator(phase, [l=='X'||l=='Y' for l in letters], [l=='Z'||l=='Y' for l in letters])
-end
-
-macro P_str(a)
-    quote _P_str($a) end
-end
-
-Base.getindex(p::PauliOperator{Tz,Tv}, i::Int) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}} = (p.xz[_div(Tve, i-1)+1] & Tve(0x1)<<_mod(Tve,i-1))!=0x0, (p.xz[end÷2+_div(Tve,i-1)+1] & Tve(0x1)<<_mod(Tve,i-1))!=0x0
-Base.getindex(p::PauliOperator{Tz,Tv}, r) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}} = PauliOperator(p.phase[], xbit(p)[r], zbit(p)[r])
-
-function Base.setindex!(p::PauliOperator{Tz,Tv}, (x,z)::Tuple{Bool,Bool}, i) where {Tz, Tve, Tv<:AbstractVector{Tve}}
-    if x
-        p.xz[_div(Tve,i-1)+1] |= Tve(0x1)<<_mod(Tve,i-1)
-    else
-        p.xz[_div(Tve,i-1)+1] &= ~(Tve(0x1)<<_mod(Tve,i-1))
-    end
-    if z
-        p.xz[end÷2+_div(Tve,i-1)+1] |= Tve(0x1)<<_mod(Tve,i-1)
-    else
-        p.xz[end÷2+_div(Tve,i-1)+1] &= ~(Tve(0x1)<<_mod(Tve,i-1))
-    end
-    p
-end
-
-Base.firstindex(p::PauliOperator) = 1
-
-Base.lastindex(p::PauliOperator) = p.nqubits
-
-Base.eachindex(p::PauliOperator) = 1:p.nqubits
-
-Base.size(pauli::PauliOperator) = (pauli.nqubits,)
-
-Base.length(pauli::PauliOperator) = pauli.nqubits
-
-nqubits(pauli::PauliOperator) = pauli.nqubits
-
-Base.:(==)(l::PauliOperator, r::PauliOperator) = r.phase==l.phase && r.nqubits==l.nqubits && r.xz==l.xz
-
-Base.hash(p::PauliOperator, h::UInt) = hash(p.phase,hash(p.nqubits,hash(p.xz, h)))
-
-Base.copy(p::PauliOperator) = PauliOperator(copy(p.phase),p.nqubits,copy(p.xz))
-
-_nchunks(i::Int,T::Type{<:Unsigned}) = 2*( (i-1) ÷ (8*sizeof(T)) + 1 )
-Base.zero(::Type{PauliOperator{Tz, Tv}}, q) where {Tz,T<:Unsigned,Tv<:AbstractVector{T}} = PauliOperator(zeros(UInt8), q, zeros(T, _nchunks(q,T)))
-Base.zero(::Type{PauliOperator}, q) = zero(PauliOperator{Array{UInt8, 0}, Vector{UInt}}, q)
-Base.zero(p::P) where {P<:PauliOperator} = zero(P, nqubits(p))
-
-"""Zero-out the phases and single-qubit operators in a [`PauliOperator`](@ref)"""
-@inline function zero!(p::PauliOperator{Tz,Tv}) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}}
-    fill!(p.xz, zero(Tve))
-    p.phase[] = 0x0
-    p
-end
+include("pauli_operator.jl")
 
 ##############################
 # Generic Tableaux
diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index 4d6c7b831..bade86ed2 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -11,18 +11,49 @@
 10. make an overleaf for a paper
 =#
 
+"""
+$(TYPEDEF)
+
+Represents mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is a pure stabilizer state.
+
+```jldoctest
+julia> StabMixture(S"-X")
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z
+𝒮𝓉𝒶𝒷
+- X
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 1.0+0.0im | + _ | + _
+
+julia> tT
+Pauli channel ρ ↦ ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† with the following branches:
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 0.853553+0.0im | + _ | + _
+ 0.0-0.353553im | + _ | + Z
+ 0.0+0.353553im | + Z | + _
+ 0.146447+0.0im | + Z | + Z
+
+julia> apply!(StabMixture(S"-X"), tT)
+A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z
+𝒮𝓉𝒶𝒷
+- X
+with ϕᵢⱼ | Pᵢ | Pⱼ:
+ 0.0-0.353553im | + _ | + Z
+ 0.0+0.353553im | + Z | + _
+ 0.853553+0.0im | + _ | + _
+ 0.146447+0.0im | + Z | + Z
+```
 
+See also: [`PauliChannel`](@ref)
+"""
 mutable struct StabMixture{T,F}
     stab::T
     destabweights::DefaultDict{Tuple{BitVector, BitVector}, F, F}
 end
 
-# TODO: figure out a way to compute whether a StabMixture instance is pure or mixed
-# TODO which hane of these (and others)
-# GeneralizedStabilizer
-# StabilizerFrame
-# GenStabFrame
-
 function StabMixture(state)
     n = nqubits(state)
     StabMixture(MixedDestabilizer(state), DefaultDict(0.0im, (falses(n),falses(n))=>1.0+0.0im)) # TODO maybe it should default to Destabilizer, not MixedDestabilizer
@@ -50,7 +81,9 @@ function Base.show(io::IO, s::StabMixture)
     println(io)
     println(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
     for ((di,dj), χ) in s.destabweights
-        println(io, "  ", χ, " | ", string(_stabmixdestab(s.stab, di)), " | ", string(_stabmixdestab(s.stab, dj)))
+        print(io, " ")
+        print(IOContext(io, :compact => true), χ)
+        println(" | ", string(_stabmixdestab(s.stab, di)), " | ", string(_stabmixdestab(s.stab, dj)))
     end
 end
 
@@ -79,8 +112,36 @@ end
 struct PauliChannel{T,S} <: AbstractPauliChannel
     paulis::T
     weights::S
+    function PauliChannel(paulis, weights)
+        n = nqubits(paulis[1][1])
+        for p in paulis
+            n == nqubits(p[1]) == nqubits(p[2]) || throw(ArgumentError(lazy"""
+            You are attempting to construct a `PauliChannel` but have provided Pauli operators
+            $(p[1]) and $(p[2])
+            that are not all of the same size (same number of qubits).
+            Please ensure that all of the Pauli operators being provided of of the same size.
+            """))
+        end
+        new{typeof(paulis),typeof(weights)}(paulis,weights)
+    end
+end
+
+function Base.show(io::IO, pc::PauliChannel)
+    println(io, "Pauli channel ρ ↦ ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† with the following branches:")
+    println(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
+    for ((di,dj), χ) in zip(pc.paulis, pc.weights)
+        print(io, " ")
+        print(IOContext(io, :compact => true), χ)
+        println(io, " | ", di, " | ", dj)
+    end
+end
+
+function embed(n,idx,pc::PauliChannel)
+    PauliChannel(map(p->(embed(n,idx,p[1]),embed(n,idx,p[2])),pc.paulis), pc.weights)
 end
 
+nqubits(pc::PauliChannel) = nqubits(pc.paulis[1][1])
+
 function apply!(state::StabMixture, gate::PauliChannel)
     dict = state.destabweights
     stab = state.stab
@@ -116,12 +177,18 @@ there are boolean vectors ``b`` and ``c`` such that
 This function returns `p`, `b`, `c`.
 
 ```
-julia> s = MixedDestabilizer(ghz(2));
-
-julia> rowdecompose(P"XY", s)
+julia> s = MixedDestabilizer(ghz(2))
+𝒟ℯ𝓈𝓉𝒶𝒷
++ Z_
++ _X
+𝒮𝓉𝒶𝒷
++ XX
++ ZZ
+
+julia> phase, destab_rows, stab_rows = QuantumClifford.rowdecompose(P"XY", s)
 (3, Bool[1, 0], Bool[1, 1])
 
-julia> -im * P"Z_" * P"XX" * P"ZZ"
+julia> im^3 * P"Z_" * P"XX" * P"ZZ"
 + XY
 ```
 """
@@ -147,3 +214,12 @@ function rowdecompose(pauli,state::Union{MixedDestabilizer, Destabilizer})
     p = mod(-Pₜ.phase[],4) # complex conjugate
     return p, b, c
 end
+
+##
+# Predefined Pauli Channels
+##
+
+const tT = PauliChannel(
+    [(I, I), (I, Z), (Z, I), (Z, Z)],
+    [cos(π/8)^2, -im*sin(π/8)*cos(π/8),  im*sin(π/8)*cos(π/8), sin(π/8)^2]
+    )
diff --git a/src/pauli_operator.jl b/src/pauli_operator.jl
new file mode 100644
index 000000000..65f86a790
--- /dev/null
+++ b/src/pauli_operator.jl
@@ -0,0 +1,176 @@
+"""
+A multi-qubit Pauli operator (``±\\{1,i\\}\\{I,Z,X,Y\\}^{\\otimes n}``).
+
+A Pauli can be constructed with the `P` custom string macro or by building
+up one through products and tensor products of smaller operators.
+
+```jldoctest
+julia> pauli3 = P"-iXYZ"
+-iXYZ
+
+julia> pauli4 = 1im * pauli3 ⊗ X
++ XYZX
+
+julia> Z*X
++iY
+```
+
+We use a typical F(2,2) encoding internally. The X and Z bits are stored
+in a single concatenated padded array of UInt chunks of a bit array.
+
+```jldoctest
+julia> p = P"-IZXY";
+
+
+julia> p.xz
+2-element Vector{UInt64}:
+ 0x000000000000000c
+ 0x000000000000000a
+```
+
+You can access the X and Z bits through getters and setters or through the
+`xview`, `zview`, `xbit`, and `zbit` functions.
+
+```jldoctest
+julia> p = P"XYZ"; p[1]
+(true, false)
+
+julia> p[1] = (true, true); p
++ YYZ
+```
+"""
+struct PauliOperator{Tz<:AbstractArray{UInt8,0}, Tv<:AbstractVector{<:Unsigned}} <: AbstractCliffordOperator
+    phase::Tz
+    nqubits::Int
+    xz::Tv
+end
+
+PauliOperator(phase::UInt8, nqubits::Int, xz::Tv) where Tv<:AbstractVector{<:Unsigned} = PauliOperator(fill(UInt8(phase),()), nqubits, xz)
+PauliOperator(phase::UInt8, x::AbstractVector{Bool}, z::AbstractVector{Bool}) = PauliOperator(fill(UInt8(phase),()), length(x), vcat(reinterpret(UInt,BitVector(x).chunks),reinterpret(UInt,BitVector(z).chunks)))
+PauliOperator(x::AbstractVector{Bool}, z::AbstractVector{Bool}) = PauliOperator(0x0, x, z)
+PauliOperator(xz::AbstractVector{Bool}) = PauliOperator(0x0, (@view xz[1:end÷2]), (@view xz[end÷2+1:end]))
+
+"""Get a view of the X part of the `UInt` array of packed qubits of a given Pauli operator."""
+function xview(p::PauliOperator)
+    @view p.xz[1:end÷2]
+end
+"""Get a view of the Y part of the `UInt` array of packed qubits of a given Pauli operator."""
+function zview(p::PauliOperator)
+    @view p.xz[end÷2+1:end]
+end
+"""Extract as a new bit array the X part of the `UInt` array of packed qubits of a given Pauli operator."""
+function xbit(p::PauliOperator)
+    one = eltype(p.xz)(1)
+    size = sizeof(eltype(p.xz))*8
+    [(word>>s)&one==one for word in xview(p) for s in 0:size-1][begin:p.nqubits]
+end
+"""Extract as a new bit array the Z part of the `UInt` array of packed qubits of a given Pauli operator."""
+function zbit(p::PauliOperator)
+    one = eltype(p.xz)(1)
+    size = sizeof(eltype(p.xz))*8
+    [(word>>s)&one==one for word in zview(p) for s in 0:size-1][begin:p.nqubits]
+end
+
+function _P_str(a)
+    letters = filter(x->occursin(x,"_IZXY"),a)
+    phase = phasedict[strip(filter(x->!occursin(x,"_IZXY"),a))]
+    PauliOperator(phase, [l=='X'||l=='Y' for l in letters], [l=='Z'||l=='Y' for l in letters])
+end
+
+macro P_str(a)
+    quote _P_str($a) end
+end
+
+Base.getindex(p::PauliOperator{Tz,Tv}, i::Int) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}} = (p.xz[_div(Tve, i-1)+1] & Tve(0x1)<<_mod(Tve,i-1))!=0x0, (p.xz[end÷2+_div(Tve,i-1)+1] & Tve(0x1)<<_mod(Tve,i-1))!=0x0
+Base.getindex(p::PauliOperator{Tz,Tv}, r) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}} = PauliOperator(p.phase[], xbit(p)[r], zbit(p)[r])
+
+function Base.setindex!(p::PauliOperator{Tz,Tv}, (x,z)::Tuple{Bool,Bool}, i) where {Tz, Tve, Tv<:AbstractVector{Tve}}
+    if x
+        p.xz[_div(Tve,i-1)+1] |= Tve(0x1)<<_mod(Tve,i-1)
+    else
+        p.xz[_div(Tve,i-1)+1] &= ~(Tve(0x1)<<_mod(Tve,i-1))
+    end
+    if z
+        p.xz[end÷2+_div(Tve,i-1)+1] |= Tve(0x1)<<_mod(Tve,i-1)
+    else
+        p.xz[end÷2+_div(Tve,i-1)+1] &= ~(Tve(0x1)<<_mod(Tve,i-1))
+    end
+    p
+end
+
+Base.firstindex(p::PauliOperator) = 1
+
+Base.lastindex(p::PauliOperator) = p.nqubits
+
+Base.eachindex(p::PauliOperator) = 1:p.nqubits
+
+Base.size(pauli::PauliOperator) = (pauli.nqubits,)
+
+Base.length(pauli::PauliOperator) = pauli.nqubits
+
+nqubits(pauli::PauliOperator) = pauli.nqubits
+
+Base.:(==)(l::PauliOperator, r::PauliOperator) = r.phase==l.phase && r.nqubits==l.nqubits && r.xz==l.xz
+
+Base.hash(p::PauliOperator, h::UInt) = hash(p.phase,hash(p.nqubits,hash(p.xz, h)))
+
+Base.copy(p::PauliOperator) = PauliOperator(copy(p.phase),p.nqubits,copy(p.xz))
+
+_nchunks(i::Int,T::Type{<:Unsigned}) = 2*( (i-1) ÷ (8*sizeof(T)) + 1 )
+Base.zero(::Type{PauliOperator{Tz, Tv}}, q) where {Tz,T<:Unsigned,Tv<:AbstractVector{T}} = PauliOperator(zeros(UInt8), q, zeros(T, _nchunks(q,T)))
+Base.zero(::Type{PauliOperator}, q) = zero(PauliOperator{Array{UInt8, 0}, Vector{UInt}}, q)
+Base.zero(p::P) where {P<:PauliOperator} = zero(P, nqubits(p))
+
+"""Zero-out the phases and single-qubit operators in a [`PauliOperator`](@ref)"""
+@inline function zero!(p::PauliOperator{Tz,Tv}) where {Tz, Tve<:Unsigned, Tv<:AbstractVector{Tve}}
+    fill!(p.xz, zero(Tve))
+    p.phase[] = 0x0
+    p
+end
+
+"""
+Embed a Pauli operator in a larger Pauli operator.
+
+```jldoctest
+julia> embed(5, 3, P"-Y")
+- __Y__
+
+julia> embed(5, (3,5), P"-YX")
+- __Y_X
+```
+"""
+function embed(n::Int, i::Int, p::PauliOperator)
+    if nqubits(p) == 1
+        pout = zero(typeof(p), n)
+        pout[i] = p[1]
+        pout.phase[] = p.phase[]
+        return pout
+    else
+        throw(ArgumentError("""
+        You are trying to embed a small Pauli operator into a larger Pauli operator.
+        However, you have not given all the positions at which the operator needs to be embedded.
+        If you are directly calling `embed`, use the form `embed(nqubits, indices::Tuple, p::PauliOperator)`.
+        If you are not using `embed` directly, then `embed` must have been incorrectly called
+        by one of the functions you have called.
+        """))
+    end
+end
+
+function embed(n::Int, indices, p::PauliOperator)
+    if nqubits(p) == length(indices)
+        pout = zero(typeof(p), n)
+        @inbounds @simd for i in eachindex(indices)
+            pout[indices[i]] = p[i]
+        end
+        pout.phase[] = p.phase[]
+        return pout
+    else
+        throw(ArgumentError(lazy"""
+        You are trying to embed a small Pauli operator into a larger Pauli operator.
+        However, you have not given all the positions at which the operator needs to be embedded.
+        The operator you are embedding is of length $(length(p)), but you have specified $(length(indices)) indices.
+        If you are not using `embed` directly, then `embed` must have been incorrectly called
+        by one of the functions you have called.
+        """))
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index a98f02cb4..16708b619 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,6 +30,7 @@ println("Starting tests with $(Threads.nthreads()) threads out of `Sys.CPU_THREA
 @doset "stabcanon"
 @doset "mul_leftright"
 @doset "inner"
+@doset "embed"
 @doset "gf2"
 @doset "projections"
 @doset "expect"
diff --git a/test/test_embed.jl b/test/test_embed.jl
new file mode 100644
index 000000000..a4c6e39e6
--- /dev/null
+++ b/test/test_embed.jl
@@ -0,0 +1,10 @@
+using Test
+using QuantumClifford
+
+@testset "embed PauliOperator" begin
+    @test embed(5,3,P"-Y") == P"-__Y__"
+    @test embed(5,(3,5),P"-YZ") == P"-__Y_Z"
+    @test embed(5,[3,5],P"-YZ") == P"-__Y_Z"
+    @test_throws ArgumentError embed(5,[3,5],P"-YZX")
+    @test_throws ArgumentError embed(5,3,P"-ZX")
+end
diff --git a/test/test_nonclifford.jl b/test/test_nonclifford.jl
index ba802b98c..eea945d5f 100644
--- a/test/test_nonclifford.jl
+++ b/test/test_nonclifford.jl
@@ -1,5 +1,5 @@
 using QuantumClifford
-using QuantumClifford: StabMixture, rowdecompose, PauliChannel, mul_left!, mul_right!
+using QuantumClifford: StabMixture, rowdecompose, PauliChannel, mul_left!, mul_right!, tT
 using Test
 using InteractiveUtils
 using Random
@@ -26,11 +26,7 @@ end
 ##
 
 @testset "PauliChannel T gate ^4 = Id" begin
-    tgate = PauliChannel(
-    [(I, I), (I, Z), (Z, I), (Z, Z)],
-    [cos(π/8)^2, -im*sin(π/8)*cos(π/8),  im*sin(π/8)*cos(π/8), sin(π/8)^2]
-    )
-
+    tgate = tT
     state = StabMixture(S"X")
 
     apply!(state, tgate)

From d6c090617eaf1628b881a954ba841b6dc7953865 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Tue, 11 Jul 2023 18:07:32 -0400
Subject: [PATCH 8/9] fixup

---
 CHANGELOG.md       |  1 +
 src/nonclifford.jl | 11 +++++++----
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4c861ae80..024914fb6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@
 ## v0.8.12 - dev
 
 - Initial implementation of non-Clifford simulation (mainly for circuits that are slightly non-Clifford, e.g. containing T gates). See `StabMixture`, `PauliChannel`, and `tT`.
+- `embed` implemented for `PauliOperator` and `PauliChannel`
 
 ## v0.8.11 - 2023-07-10
 
diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index bade86ed2..505f147d6 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -79,11 +79,12 @@ function Base.show(io::IO, s::StabMixture)
     println(io, "A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is")
     show(io,s.stab)
     println(io)
-    println(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
+    print(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
     for ((di,dj), χ) in s.destabweights
+        println(io)
         print(io, " ")
         print(IOContext(io, :compact => true), χ)
-        println(" | ", string(_stabmixdestab(s.stab, di)), " | ", string(_stabmixdestab(s.stab, dj)))
+        print(io, " | ", string(_stabmixdestab(s.stab, di)), " | ", string(_stabmixdestab(s.stab, dj)))
     end
 end
 
@@ -109,6 +110,7 @@ struct TGate <: AbstractPauliChannel
     qubit::Int
 end
 
+"""A Pauli channel datastructure, mainly for use with [`StabMixture`](@ref)"""
 struct PauliChannel{T,S} <: AbstractPauliChannel
     paulis::T
     weights::S
@@ -128,11 +130,12 @@ end
 
 function Base.show(io::IO, pc::PauliChannel)
     println(io, "Pauli channel ρ ↦ ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† with the following branches:")
-    println(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
+    print(io, "with ϕᵢⱼ | Pᵢ | Pⱼ:")
     for ((di,dj), χ) in zip(pc.paulis, pc.weights)
+        println(io)
         print(io, " ")
         print(IOContext(io, :compact => true), χ)
-        println(io, " | ", di, " | ", dj)
+        print(io, " | ", di, " | ", dj)
     end
 end
 

From bef98600961ea82c171f2a8613bfefdc10fb96d4 Mon Sep 17 00:00:00 2001
From: Stefan Krastanov <stefan@krastanov.org>
Date: Tue, 11 Jul 2023 18:24:21 -0400
Subject: [PATCH 9/9] fix ambiguity

---
 src/nonclifford.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/nonclifford.jl b/src/nonclifford.jl
index 505f147d6..47cad9996 100644
--- a/src/nonclifford.jl
+++ b/src/nonclifford.jl
@@ -139,7 +139,7 @@ function Base.show(io::IO, pc::PauliChannel)
     end
 end
 
-function embed(n,idx,pc::PauliChannel)
+function embed(n::Int,idx,pc::PauliChannel)
     PauliChannel(map(p->(embed(n,idx,p[1]),embed(n,idx,p[2])),pc.paulis), pc.weights)
 end