QuantumSavory · Krastanov · Jul 11, 2023 · Jun 9, 2023 · Jun 9, 2023 · Jun 11, 2023
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,6 +6,11 @@
 
 # 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`.
+- `embed` implemented for `PauliOperator` and `PauliChannel`
+
 ## 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/Project.toml b/Project.toml
@@ -32,7 +32,7 @@ QuantumCliffordQuantikzExt = "Quantikz"
 
 [compat]
 Combinatorics = "1.0"
-DataStructures = "0.18.0"
+DataStructures = "0.18"
 DocStringExtensions = "0.8, 0.9"
 Graphs = "1.4.1"
 HostCPUFeatures = "0.1.6"

diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
@@ -8,14 +8,14 @@ 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
 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
@@ -1289,6 +1164,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/mul_leftright.jl b/src/mul_leftright.jl
@@ -185,3 +185,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