Merge branch 'master' into rs

.typos.toml

@@ -1,5 +1,6 @@
 [default.extend-words]
 ket = "ket"
+anc = "anc"
 
 [type.ipynb]
 # It detects false possitives in the base64 encoded images inside notebooks

CHANGELOG.md

@@ -5,10 +5,12 @@
 
 # News
 
-## v0.9.4 - 2024-04-14
-
+## v0.9.4 - dev
 - Added the classical Reed-Solomon code to the ECC module
 - Added the classical Bose–Chaudhuri–Hocquenghem (BCH) code to the ECC module
+- Gate errors are now conveniently supported by the various ECC benchmark setups in the `ECC` module.
+- Remove printing of spurious debug info from the PyBP decoder. 
+
 
 ## v0.9.3 - 2024-04-10
 

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org> and QuantumSavory community members"]
-version = "0.9.3"
+version = "0.9.4"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -41,7 +41,7 @@ QuantumCliffordQOpticsExt = "QuantumOpticsBase"
 QuantumCliffordQuantikzExt = "Quantikz"
 
 [compat]
-CUDA = "4.4.0"
+CUDA = "4.4.0, 5"
 Combinatorics = "1.0"
 DataStructures = "0.18"
 DocStringExtensions = "0.9"

docs/.gitignore

@@ -1,2 +1,3 @@
 build/
 site/
+.CondaPkg/

docs/src/graphs.md

@@ -20,7 +20,7 @@ graphstate(ghz(4))[1]
 
 ```@eval
 using Random; Random.seed!(1); using QuantumClifford, GraphMakie, CairoMakie;
-f = Figure(resolution=(200,200))
+f = Figure(size=(200,200))
 a = Axis(f[1,1])
 graphplot!(a,graphstate(ghz(4))[1])
 hidedecorations!(a); hidespines!(a)

docs/src/references.bib

@@ -272,6 +272,8 @@ @article{kitaev2003fault
 	abstract = {A two-dimensional quantum system with anyonic excitations can be considered as a quantum computer. Unitary transformations can be performed by moving the excitations around each other. Measurements can be performed by joining excitations in pairs and observing the result of fusion. Such computation is fault-tolerant by its physical nature.},
 	pages = {2--30},
 	number = {1},
+  year = {2003},
+  journal={Annals of Physics},
 	journaltitle = {Annals of Physics},
 	shortjournal = {Annals of Physics},
 	author = {Kitaev, A.Yu.}
@@ -286,6 +288,8 @@ @article{fowler2012surface
 	shorttitle = {Surface codes},
 	pages = {032324},
 	number = {3},
+  year = {2012},
+  journal={Physical Review A},
 	journaltitle = {Physical Review A},
 	shortjournal = {Phys. Rev. A},
 	author = {Fowler, Austin G. and Mariantoni, Matteo and Martinis, John M. and Cleland, Andrew N.},

ext/QuantumCliffordGPUExt/apply.jl

@@ -4,19 +4,19 @@ using QuantumClifford: getxbit, getzbit, setxbit, setzbit
 # SingleQubitOperator Kernel
 ##############################
 
-function single_qubit_gpu_kernel(xzs::DeviceMatrix{Tme},
-                                 phases::DeviceVector{Tmz},
+function single_qubit_gpu_kernel(xzs::DeviceMatrix{Tₘₑ},
+                                 phases::DeviceVector{Tₚₑ},
                                  op::SingleQubitOperator,
                                  rows::Int,
-                                 compute_phases::Bool) where {Tme <: Unsigned, Tmz <: Unsigned}
+                                 compute_phases::Bool) where {Tₘₑ <: Unsigned, Tₚₑ <: Unsigned}
     idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     if idx > rows
         return nothing
     end
     c = op.q
     r = idx
-    sh = QuantumClifford.getshift(Tme, c)
-    xx,zx,xz,zz = Tme.((op.xx,op.zx,op.xz,op.zz)) .<< sh # maybe putting this in parent function call will improve speed?
+    sh = QuantumClifford.getshift(Tₘₑ, c)
+    xx,zx,xz,zz = Tₘₑ.((op.xx,op.zx,op.xz,op.zz)) .<< sh # maybe putting this in parent function call will improve speed?
     anticom = ~iszero((~zz & xz & ~xx & zx) | ( zz & ~xz & xx & zx) | (zz &  xz & xx & ~zx))
 
     x = getxbit(xzs, r, c)
@@ -45,20 +45,20 @@ end
 # AbstractSingleQubitOperator Kernel
 ##############################
 
-function abstract_single_qubit_gpu_kernel(xzs::DeviceMatrix{Tme},
-                                 phases::DeviceVector{Tmz},
+function abstract_single_qubit_gpu_kernel(xzs::DeviceMatrix{Tₘₑ},
+                                 phases::DeviceVector{Tₚₑ},
                                  gate::QuantumClifford.AbstractSingleQubitOperator,
                                  rows::Int,
-                                 compute_phases::Bool) where {Tme <: Unsigned, Tmz <: Unsigned}
+                                 compute_phases::Bool) where {Tₘₑ <: Unsigned, Tₚₑ <: Unsigned}
     idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     if idx > rows
         return nothing
     end
     c = gate.q
     r = idx
 
-    x::Tme = getxbit(xzs, r, c)
-    z::Tme = getzbit(xzs, r, c)
+    x::Tₘₑ = getxbit(xzs, r, c)
+    z::Tₘₑ = getzbit(xzs, r, c)
     x,z,phase::Bool = QuantumClifford.qubit_kernel(gate,x,z)
     setxbit(xzs, r, c, x)
     setzbit(xzs, r, c, z)
@@ -70,26 +70,26 @@ end
 # AbstractTwoQubitOperator Kernel
 ##############################
 
-function two_qubit_gpu_kernel(xzs::DeviceMatrix{Tme},
-                              phases::DeviceVector{Tze},
+function two_qubit_gpu_kernel(xzs::DeviceMatrix{Tₘₑ},
+                              phases::DeviceVector{Tₚ},
                               gate::QuantumClifford.AbstractTwoQubitOperator,
                               rows::Int,
-                              compute_phases::Bool=true) where {Tme <: Unsigned, Tze <: Unsigned}
+                              compute_phases::Bool=true) where {Tₘₑ <: Unsigned, Tₚ <: Unsigned}
     idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
     if idx > rows
         return nothing
     end
 
     q1 = gate.q1
     q2 = gate.q2
-    shift = QuantumClifford.getshift(Tme, q1) - QuantumClifford.getshift(Tme, q2)
+    shift = QuantumClifford.getshift(Tₘₑ, q1) - QuantumClifford.getshift(Tₘₑ, q2)
     r = idx;
 
-    _x1::Tme = getxbit(xzs, r, q1)
-    _z1::Tme = getzbit(xzs, r, q1)
-    _x2::Tme = getxbit(xzs, r, q2)<<shift
-    _z2::Tme = getzbit(xzs, r, q2)<<shift
-    x1::Tme,z1::Tme,x2::Tme,z2::Tme,phase::Bool = QuantumClifford.qubit_kernel(gate,_x1,_z1,_x2,_z2) # Most `qubit_kernel` functions are defined by a `qubitop2` macro
+    _x1::Tₘₑ = getxbit(xzs, r, q1)
+    _z1::Tₘₑ = getzbit(xzs, r, q1)
+    _x2::Tₘₑ = getxbit(xzs, r, q2)<<shift
+    _z2::Tₘₑ = getzbit(xzs, r, q2)<<shift
+    x1::Tₘₑ,z1::Tₘₑ,x2::Tₘₑ,z2::Tₘₑ,phase::Bool = QuantumClifford.qubit_kernel(gate,_x1,_z1,_x2,_z2) # Most `qubit_kernel` functions are defined by a `qubitop2` macro
     setxbit(xzs, r, q1, x1, 0)
     setzbit(xzs, r, q1, z1, 0)
     setxbit(xzs, r, q2, x2, -shift)
@@ -125,8 +125,8 @@ function _apply!(stab::StabilizerGPU{T},
     return stab
 end
 
-function _apply!(stab::StabilizerGPU{T}, 
-                 gate::G; 
+function _apply!(stab::StabilizerGPU{T},
+                 gate::G;
                  phases::Val{B}=Val(true)) where {B, G<:QuantumClifford.AbstractTwoQubitOperator, T <: Unsigned}
 
     rows = size(stab, 1)

ext/QuantumCliffordGPUExt/apply_noise.jl

@@ -18,11 +18,11 @@ function applynoise!(frame::PauliFrameGPU{T},noise::UnbiasedUncorrelatedNoise,i:
 end
 
 
-function applynoise_kernel(xzs::DeviceMatrix{Tme},
+function applynoise_kernel(xzs::DeviceMatrix{Tₘₑ},
     p::Real,
     ibig::Int,
-    ismallm::Tme,
-    rows::Int) where {Tme <: Unsigned}
+    ismallm::Tₘₑ,
+    rows::Int) where {Tₘₑ <: Unsigned}
 
     f = (blockIdx().x - 1) * blockDim().x + threadIdx().x;
     if f > rows

ext/QuantumCliffordGPUExt/pauli_frames.jl

@@ -2,16 +2,16 @@
 # sMZ
 ##############################
 
-function apply_sMZ_kernel!(xzs::DeviceMatrix{Tme},
+function apply_sMZ_kernel!(xzs::DeviceMatrix{Tₘₑ},
                           measurements::DeviceMatrix{Bool},
                           op::sMZ,
                           ibig::Int,
-                          ismallm::Tme,
-                          rows::Int) where {Tme <: Unsigned} 
+                          ismallm::Tₘₑ,
+                          rows::Int) where {Tₘₑ <: Unsigned}
     f = (blockIdx().x - 1) * blockDim().x + threadIdx().x;
     if f > rows
         return nothing
-    end    
+    end
     should_flip = !iszero(xzs[ibig,f] & ismallm)
     measurements[f,op.bit] = should_flip
     return nothing
@@ -33,12 +33,12 @@ end
 # sMRZ
 ##############################
 
-function apply_sMRZ_kernel!(xzs::DeviceMatrix{Tme},
+function apply_sMRZ_kernel!(xzs::DeviceMatrix{Tₘₑ},
                           measurements::DeviceMatrix{Bool},
                           op::QuantumClifford.sMRZ,
                           ibig::Int, # todo change to Int
-                          ismallm::Tme,
-                          rows::Int) where {Tme <: Unsigned} 
+                          ismallm::Tₘₑ,
+                          rows::Int) where {Tₘₑ <: Unsigned}
     f = (blockIdx().x - 1) * blockDim().x + threadIdx().x;
     if f > rows
         return nothing

ext/QuantumCliffordGPUExt/types.jl

@@ -20,13 +20,13 @@ DeviceMatrix{T} = Union{CuDeviceArray{T, 2, 1}, Adjoint{T, CuDeviceArray{T, 2, 1
 
 
 function getTableauGPU(GPUValue, GPUVector, GPUMatrix)
-    TableauGPU{T} = QuantumClifford.Tableau{Tzv, Tm} where {T <: Unsigned, Tzv <: GPUVector{PhaseInnerType}, Tm <: GPUMatrix{T}}
+    TableauGPU{T} = QuantumClifford.Tableau{Tₚᵥ, Tₘ} where {T <: Unsigned, Tₚᵥ <: GPUVector{PhaseInnerType}, Tₘ <: GPUMatrix{T}}
 end
 function getStabilizerGPU(GPUValue, GPUVector, GPUMatrix, GPUTableau)
     StabilizerGPU{T} = QuantumClifford.Stabilizer{<:GPUTableau{T}} where {T <: Unsigned}
 end
 function getPauliOperatorGPU(GPUValue, GPUVector, GPUMatrix, GPUTableau)
-    PauliOperatorGPU{T} = QuantumClifford.PauliOperator{Tz, Tv} where {T <: Unsigned, Tz<:GPUValue{PhaseInnerType}, Tv<:GPUVector{T}}
+    PauliOperatorGPU{T} = QuantumClifford.PauliOperator{Tₚ, Tᵥ} where {T <: Unsigned, Tₚ<:GPUValue{PhaseInnerType}, Tᵥ<:GPUVector{T}}
 end
 
 # todo. type definition here is stronger than the code in pauliframes.jl  this will cause serious problems

ext/QuantumCliffordMakieExt/QuantumCliffordMakieExt.jl

@@ -57,8 +57,8 @@ function stabilizerplot_axis(subfig, s; colorbar=true, args...)
     Makie.hidedecorations!(ax)
     Makie.hidespines!(ax)
     ax.aspect = Makie.DataAspect()
-    colorbar && Makie.Colorbar(subfig[2, 1], p, ticks = (0:3, ["I", "X", "Z", "Y"]), vertical = false, flipaxis = false)
-    Makie.colsize!(subfig.layout, 1, Makie.Aspect(1, min(1,size(s,2)/size(s,1))))
+    colorbar && Makie.Colorbar(subfig[1, 2], p, ticks = ((0.5:3.51)*3/4, ["I", "X", "Z", "Y"]), vertical = true, flipaxis = true)
+    #Makie.colsize!(subfig.layout, 1, Makie.Aspect(1, min(1,size(s,2)/size(s,1))))
     subfig,ax,p
 end
 

ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl

@@ -71,7 +71,6 @@ parity_checks(d::PyBP) = d.H
 function decode(d::PyBP, syndrome_sample)
     row_x = syndrome_sample[1:d.nx] # TODO These copies and indirections might be costly!
     row_z = syndrome_sample[d.nx+1:end]
-    @show (size(row_x), size(row_z))
     guess_z_errors = PythonCall.PyArray(d.pyx.decode(np.array(row_x)))
     guess_x_errors = PythonCall.PyArray(d.pyz.decode(np.array(row_z)))
     vcat(guess_x_errors, guess_z_errors)

ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl

@@ -71,7 +71,7 @@ function Quantikz.QuantikzOp(op::Reset) # TODO This is complicated because quant
         Quantikz.Initialize("$str",affectedqubits(op)) # TODO make Quantikz work with tuples and remove the collect
     end
 end
-Quantikz.QuantikzOp(op::NoiseOp) = Quantikz.Noise(op.indices)
+Quantikz.QuantikzOp(op::NoiseOp) = Quantikz.Noise(collect(op.indices))
 Quantikz.QuantikzOp(op::NoiseOpAll) = Quantikz.NoiseAll()
 Quantikz.QuantikzOp(op::ClassicalXORConcreteWorkaround) = Quantikz.ClassicalDecision(sort([op.store, op.bits...]))
 

src/QuantumClifford.jl

@@ -130,18 +130,18 @@ include("pauli_operator.jl")
 """Internal Tableau type for storing a list of Pauli operators in a compact form.
 No special semantic meaning is attached to this type, it is just a convenient way to store a list of Pauli operators.
 E.g. it is not used to represent a stabilizer state, or a stabilizer group, or a Clifford circuit."""
-struct Tableau{Tzv<:AbstractVector{UInt8}, Tm<:AbstractMatrix{<:Unsigned}}
-    phases::Tzv
+struct Tableau{Tₚᵥ<:AbstractVector{UInt8}, Tₘ<:AbstractMatrix{<:Unsigned}}
+    phases::Tₚᵥ
     nqubits::Int
-    xzs::Tm
+    xzs::Tₘ
 end
 
-function Tableau(paulis::AbstractVector{PauliOperator{Tz,Tv}}) where {Tz<:AbstractArray{UInt8,0},Tve<:Unsigned,Tv<:AbstractVector{Tve}}
+function Tableau(paulis::AbstractVector{PauliOperator{Tₚ,Tᵥ}}) where {Tₚ<:AbstractArray{UInt8,0},Tᵥₑ<:Unsigned,Tᵥ<:AbstractVector{Tᵥₑ}}
     r = length(paulis)
     n = nqubits(paulis[1])
-    tab = zero(Tableau{Vector{UInt8},Matrix{Tve}},r,n)::Tableau{Vector{UInt8},Matrix{Tve}} # typeassert for JET
+    tab = zero(Tableau{Vector{UInt8},Matrix{Tᵥₑ}},r,n)::Tableau{Vector{UInt8},Matrix{Tᵥₑ}} # typeassert for JET
     for i in eachindex(paulis)
-        tab[i] = paulis[i]::PauliOperator{Tz,Tv} # typeassert for JET
+        tab[i] = paulis[i]::PauliOperator{Tₚ,Tᵥ} # typeassert for JET
     end
     tab
 end
@@ -172,9 +172,9 @@ end
 Base.getindex(tab::Tableau, i::Int) = PauliOperator(tab.phases[i], nqubits(tab), tab.xzs[:,i])
 @inline function Base.getindex(tab::Tableau, r::Int, c::Int)
      # TODO this has code repetition with the Pauli getindex
-     Tme = eltype(tab.xzs)
-     x = (tab.xzs[_div(Tme,c-1)+1,r] & Tme(0x1)<<_mod(Tme,c-1))!=0x0
-     z = (tab.xzs[end÷2+_div(Tme,c-1)+1,r] & Tme(0x1)<<_mod(Tme,c-1))!=0x0
+     Tₘₑ = eltype(tab.xzs)
+     x = (tab.xzs[_div(Tₘₑ,c-1)+1,r] & Tₘₑ(0x1)<<_mod(Tₘₑ,c-1))!=0x0
+     z = (tab.xzs[end÷2+_div(Tₘₑ,c-1)+1,r] & Tₘₑ(0x1)<<_mod(Tₘₑ,c-1))!=0x0
      return (x,z)
 end
 Base.getindex(tab::Tableau, r) = Tableau(tab.phases[r], nqubits(tab), tab.xzs[:,r])
@@ -200,16 +200,16 @@ function Base.setindex!(tab::Tableau, t::Tableau, i)
     tab
 end
 
-function Base.setindex!(tab::Tableau{Tzv,Tm}, (x,z)::Tuple{Bool,Bool}, i, j) where {Tzv<:AbstractVector{UInt8}, Tme<:Unsigned, Tm<:AbstractMatrix{Tme}} # TODO this has code repetitions with the Pauli setindex
+function Base.setindex!(tab::Tableau{Tₚᵥ,Tₘ}, (x,z)::Tuple{Bool,Bool}, i, j) where {Tₚᵥ<:AbstractVector{UInt8}, Tₘₑ<:Unsigned, Tₘ<:AbstractMatrix{Tₘₑ}} # TODO this has code repetitions with the Pauli setindex
     if x
-        tab.xzs[_div(Tme,j-1)+1,        i] |= Tme(0x1)<<_mod(Tme,j-1)
+        tab.xzs[_div(Tₘₑ,j-1)+1,        i] |= Tₘₑ(0x1)<<_mod(Tₘₑ,j-1)
     else
-        tab.xzs[_div(Tme,j-1)+1,        i] &= ~(Tme(0x1)<<_mod(Tme,j-1))
+        tab.xzs[_div(Tₘₑ,j-1)+1,        i] &= ~(Tₘₑ(0x1)<<_mod(Tₘₑ,j-1))
     end
     if z
-        tab.xzs[end÷2+_div(Tme,j-1)+1, i] |= Tme(0x1)<<_mod(Tme,j-1)
+        tab.xzs[end÷2+_div(Tₘₑ,j-1)+1, i] |= Tₘₑ(0x1)<<_mod(Tₘₑ,j-1)
     else
-        tab.xzs[end÷2+_div(Tme,j-1)+1, i] &= ~(Tme(0x1)<<_mod(Tme,j-1))
+        tab.xzs[end÷2+_div(Tₘₑ,j-1)+1, i] &= ~(Tₘₑ(0x1)<<_mod(Tₘₑ,j-1))
     end
     tab
 end
@@ -235,7 +235,7 @@ Base.hash(t::Tableau, h::UInt) = hash(t.nqubits, hash(t.phases, hash(t.xzs, h)))
 
 Base.copy(t::Tableau) = Tableau(copy(t.phases), t.nqubits, copy(t.xzs))
 
-function Base.zero(::Type{Tableau{Tzv, Tm}}, r, q) where {Tzv,T<:Unsigned,Tm<:AbstractMatrix{T}}
+function Base.zero(::Type{Tableau{Tₚᵥ, Tₘ}}, r, q) where {Tₚᵥ,T<:Unsigned,Tₘ<:AbstractMatrix{T}}
     phases = zeros(UInt8,r)
     xzs = zeros(UInt, _nchunks(q,T), r)
     Tableau(phases, q, xzs)::Tableau{Vector{UInt8},Matrix{UInt}}
@@ -354,8 +354,8 @@ struct Stabilizer{T<:Tableau} <: AbstractStabilizer
     tab::T
 end
 
-Stabilizer(phases::Tzv, nqubits::Int, xzs::Tm) where {Tzv<:AbstractVector{UInt8}, Tm<:AbstractMatrix{<:Unsigned}} = Stabilizer(Tableau(phases, nqubits, xzs))
-Stabilizer(paulis::AbstractVector{PauliOperator{Tz,Tv}}) where {Tz,Tv} = Stabilizer(Tableau(paulis))
+Stabilizer(phases::Tₚᵥ, nqubits::Int, xzs::Tₘ) where {Tₚᵥ<:AbstractVector{UInt8}, Tₘ<:AbstractMatrix{<:Unsigned}} = Stabilizer(Tableau(phases, nqubits, xzs))
+Stabilizer(paulis::AbstractVector{PauliOperator{Tₚ,Tᵥ}}) where {Tₚ,Tᵥ} = Stabilizer(Tableau(paulis))
 Stabilizer(phases::AbstractVector{UInt8}, xs::AbstractMatrix{Bool}, zs::AbstractMatrix{Bool}) = Stabilizer(Tableau(phases, xs, zs))
 Stabilizer(phases::AbstractVector{UInt8}, xzs::AbstractMatrix{Bool}) = Stabilizer(Tableau(phases, xzs))
 Stabilizer(xs::AbstractMatrix{Bool}, zs::AbstractMatrix{Bool}) = Stabilizer(Tableau(xs,zs))
@@ -832,13 +832,13 @@ end
 # TODO is this used anywhere?
 """Swap two columns of a stabilizer in place."""
 @inline function colswap!(s::Tableau, i, j)
-    Tme = eltype(s.xzs)
-    lowbit = Tme(1)
-    ibig = _div(Tme,i-1)+1
-    ismall = _mod(Tme,i-1)
+    Tₘₑ = eltype(s.xzs)
+    lowbit = Tₘₑ(1)
+    ibig = _div(Tₘₑ,i-1)+1
+    ismall = _mod(Tₘₑ,i-1)
     ismallm = lowbit<<(ismall)
-    jbig = _div(Tme,j-1)+1
-    jsmall = _mod(Tme,j-1)
+    jbig = _div(Tₘₑ,j-1)+1
+    jsmall = _mod(Tₘₑ,j-1)
     jsmallm = lowbit<<(jsmall)
     for off in [0,size(s.xzs,2)÷2]
         ibig += off

src/dense_cliffords.jl

@@ -136,7 +136,6 @@ function _apply!(stab::AbstractStabilizer, c::CliffordOperator; phases::Val{B}=V
 end
 
 # TODO Added a lot of type assertions to help Julia infer types, but they are much too strict for cases where bitpacking varies (check tests)
-#@inline function apply_row_kernel!(new_stabrow::PauliOperator{Array{UInt8,0},Vector{Tme}}, row::Int, s_tab::Stabilizer{Tv,Tm}, c_tab::Stabilizer{Tv,Tm}; phases=true) where {Tme,Tv<:AbstractVector{UInt8},Tm<:AbstractMatrix{Tme}}
 @inline function apply_row_kernel!(new_stabrow, row, s_tab, c_tab; phases::Val{B}=Val(true)) where B
     B && (new_stabrow.phase[] = s_tab.phases[row])
     n = nqubits(c_tab)