Merge branch 'QuantumSavory:master' into i113

CHANGELOG.md

@@ -5,6 +5,12 @@
 
 # News
 
+## v0.9.11
+
+- `hcat` of Tableaux objects
+- `QuantumReedMuller` codes added to the ECC module
+- **(breaking)** change the convention for how to provide a representation function in the constructor of `LPCode` -- strictly speaking a breaking change, but this is not an API that is publicly used in practice
+
 ## v0.9.10 - 2024-09-26
 
 - The lifted product class of quantum LDPC codes is implemented in the ECC submodule.

Project.toml

@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <[email protected]> and QuantumSavory community members"]
-version = "0.9.10"
+version = "0.9.11"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -48,15 +48,15 @@ Combinatorics = "1.0"
 DataStructures = "0.18"
 DocStringExtensions = "0.9"
 Graphs = "1.9"
-Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33"
+Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
 HostCPUFeatures = "0.1.6"
 ILog2 = "0.2.3"
 InteractiveUtils = "1.9"
 LDPCDecoders = "0.3.1"
 LinearAlgebra = "1.9"
 MacroTools = "0.5.9"
 Makie = "0.20, 0.21"
-Nemo = "^0.42.1, 0.43, 0.44, 0.45, 0.46"
+Nemo = "0.42.1, 0.43, 0.44, 0.45, 0.46, 0.47"
 Plots = "1.38.0"
 PrecompileTools = "1.2"
 PyQDecoders = "0.2.1"

docs/src/references.bib

@@ -454,3 +454,36 @@ @article{raveendran2022finite
   issn = {2521-327X},
   doi = {10.22331/q-2022-07-20-767},
 }
+
+@article{steane1999quantum,
+  title={Quantum reed-muller codes},
+  author={Steane, Andrew M},
+  journal={IEEE Transactions on Information Theory},
+  volume={45},
+  number={5},
+  pages={1701--1703},
+  year={1999},
+  publisher={IEEE}
+}
+
+@article{campbell2012magic,
+  title={Magic-state distillation in all prime dimensions using quantum reed-muller codes},
+  author={Campbell, Earl T and Anwar, Hussain and Browne, Dan E},
+  journal={Physical Review X},
+  volume={2},
+  number={4},
+  pages={041021},
+  year={2012},
+  publisher={APS}
+}
+
+@article{anderson2014fault,
+  title={Fault-tolerant conversion between the steane and reed-muller quantum codes},
+  author={Anderson, Jonas T and Duclos-Cianci, Guillaume and Poulin, David},
+  journal={Physical review letters},
+  volume={113},
+  number={8},
+  pages={080501},
+  year={2014},
+  publisher={APS}
+}

docs/src/references.md

@@ -37,6 +37,9 @@ For quantum code construction routines:
 - [kitaev2003fault](@cite)
 - [fowler2012surface](@cite)
 - [knill1996concatenated](@cite)
+- [steane1999quantum](@cite)
+- [campbell2012magic](@cite)
+- [anderson2014fault](@cite)
 
 For classical code construction routines:
 - [muller1954application](@cite)

ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl

@@ -6,7 +6,8 @@ import QuantumClifford, LinearAlgebra
 import Hecke: Group, GroupElem, AdditiveGroup, AdditiveGroupElem,
     GroupAlgebra, GroupAlgebraElem, FqFieldElem, representation_matrix, dim, base_ring,
     multiplication_table, coefficients, abelian_group, group_algebra
-import Nemo: characteristic, matrix_repr, GF, ZZ
+import Nemo
+import Nemo: characteristic, matrix_repr, GF, ZZ, lift
 
 import QuantumClifford.ECC: AbstractECC, CSS, ClassicalCode,
     hgp, code_k, code_n, code_s, iscss, parity_checks, parity_checks_x, parity_checks_z, parity_checks_xz,

ext/QuantumCliffordHeckeExt/lifted.jl

@@ -22,10 +22,11 @@ The default `GA` is the group algebra of `A[1, 1]`, the default representation `
 
 ## The representation function `repr`
 
-In this struct, we use the default representation function `default_repr` to convert a `GF(2)`-group algebra element to a binary matrix.
+We use the default representation function `Hecke.representation_matrix` to convert a `GF(2)`-group algebra element to a binary matrix.
 The default representation, provided by `Hecke`, is the permutation representation.
 
-We also accept a custom representation function.
+We also accept a custom representation function (the `repr` field of the constructor).
+Whatever the representation, the matrix elements need to be convertible to Integers (e.g. permit `lift(ZZ, ...)`).
 Such a customization would be useful to reduce the number of bits required by the code construction.
 
 For example, if we use a D4 group for lifting, our default representation will be `8×8` permutation matrices,
@@ -54,14 +55,12 @@ struct LiftedCode <: ClassicalCode
     end
 end
 
-default_repr(y::GroupAlgebraElem{FqFieldElem, <: GroupAlgebra}) = Matrix((x -> Bool(Int(lift(ZZ, x)))).(representation_matrix(y)))
-
 """
 `LiftedCode` constructor using the default `GF(2)` representation (coefficients converted to a permutation matrix by `representation_matrix` provided by Hecke).
 """ # TODO doctest example
 function LiftedCode(A::Matrix{GroupAlgebraElem{FqFieldElem, <: GroupAlgebra}}; GA::GroupAlgebra=parent(A[1,1]))
     !(characteristic(base_ring(A[1, 1])) == 2) && error("The default permutation representation applies only to GF(2) group algebra; otherwise, a custom representation function should be provided")
-    LiftedCode(A; GA=GA, repr=default_repr)
+    LiftedCode(A; GA=GA, repr=representation_matrix)
 end
 
 # TODO document and doctest example
@@ -71,7 +70,7 @@ function LiftedCode(group_elem_array::Matrix{<: GroupOrAdditiveGroupElem}; GA::G
         A[i, j] = GA[A[i, j]]
     end
     if repr === nothing
-        return LiftedCode(A; GA=GA, repr=default_repr)
+        return LiftedCode(A; GA=GA, repr=representation_matrix)
     else
         return LiftedCode(A; GA=GA, repr=repr)
     end
@@ -85,11 +84,16 @@ function LiftedCode(shift_array::Matrix{Int}, l::Int; GA::GroupAlgebra=group_alg
             A[i, j] = GA[shift_array[i, j]%l+1]
         end
     end
-    return LiftedCode(A; GA=GA, repr=default_repr)
+    return LiftedCode(A; GA=GA, repr=representation_matrix)
 end
 
-function lift(repr::Function, mat::GroupAlgebraElemMatrix)
-    vcat([hcat([repr(mat[i, j]) for j in axes(mat, 2)]...) for i in axes(mat, 1)]...)
+lift_to_bool(x) = Bool(Int(lift(ZZ,x)))
+
+function concat_lift_repr(repr, mat)
+    x = repr.(mat)
+    y = hvcat(size(x,2), transpose(x)...)
+    z = Matrix(lift_to_bool.(y))
+    return z
 end
 
 function parity_checks(c::LiftedCode)

ext/QuantumCliffordHeckeExt/lifted_product.jl

@@ -72,7 +72,7 @@ julia> code_n(c2), code_k(c2)
 
 ## The representation function
 
-In this struct, we use the default representation function `default_repr` to convert a `GF(2)`-group algebra element to a binary matrix.
+We use the default representation function `Hecke.representation_matrix` to convert a `GF(2)`-group algebra element to a binary matrix.
 The default representation, provided by `Hecke`, is the permutation representation.
 
 We also accept a custom representation function as detailed in [`LiftedCode`](@ref).
@@ -107,24 +107,24 @@ end
 
 # TODO document and doctest example
 function LPCode(A::FqFieldGroupAlgebraElemMatrix, B::FqFieldGroupAlgebraElemMatrix; GA::GroupAlgebra=parent(A[1,1]))
-    LPCode(LiftedCode(A; GA=GA, repr=default_repr), LiftedCode(B; GA=GA, repr=default_repr); GA=GA, repr=default_repr)
+    LPCode(LiftedCode(A; GA=GA, repr=representation_matrix), LiftedCode(B; GA=GA, repr=representation_matrix); GA=GA, repr=representation_matrix)
 end
 
 # TODO document and doctest example
 function LPCode(group_elem_array1::Matrix{<: GroupOrAdditiveGroupElem}, group_elem_array2::Matrix{<: GroupOrAdditiveGroupElem}; GA::GroupAlgebra=group_algebra(GF(2), parent(group_elem_array1[1,1])))
-    LPCode(LiftedCode(group_elem_array1; GA=GA), LiftedCode(group_elem_array2; GA=GA); GA=GA, repr=default_repr)
+    LPCode(LiftedCode(group_elem_array1; GA=GA), LiftedCode(group_elem_array2; GA=GA); GA=GA, repr=representation_matrix)
 end
 
 # TODO document and doctest example
 function LPCode(shift_array1::Matrix{Int}, shift_array2::Matrix{Int}, l::Int; GA::GroupAlgebra=group_algebra(GF(2), abelian_group(l)))
-    LPCode(LiftedCode(shift_array1, l; GA=GA), LiftedCode(shift_array2, l; GA=GA); GA=GA, repr=default_repr)
+    LPCode(LiftedCode(shift_array1, l; GA=GA), LiftedCode(shift_array2, l; GA=GA); GA=GA, repr=representation_matrix)
 end
 
 iscss(::Type{LPCode}) = true
 
 function parity_checks_xz(c::LPCode)
     hx, hz = hgp(c.A, c.B')
-    hx, hz = lift(c.repr, hx), lift(c.repr, hz)
+    hx, hz = concat_lift_repr(c.repr,hx), concat_lift_repr(c.repr,hz)
     return hx, hz
 end
 

ext/QuantumCliffordHeckeExt/types.jl

@@ -15,7 +15,7 @@ Compute the adjoint of a group algebra element.
 The adjoint is defined as the conjugate of the element in the group algebra,
 i.e. the inverse of the element in the associated group.
 """
-function _adjoint(a::GroupAlgebraElem{T}) where T # TODO Is this used? Should it be deleted?
+function Base.adjoint(a::GroupAlgebraElem{T}) where T # TODO we would like to use Base.adjoint, but that would be type piracy. Upstream this to Nemo or Hecke or AbstractAlgebra
     A = parent(a)
     d = dim(A)
     v = Vector{T}(undef, d)

ext/QuantumCliffordLDPCDecodersExt/QuantumCliffordLDPCDecodersExt.jl

@@ -74,16 +74,16 @@ parity_checks(d::BeliefPropDecoder) = d.H
 parity_checks(d::BitFlipDecoder) = d.H
 
 function decode(d::BeliefPropDecoder, syndrome_sample)
-    row_x = syndrome_sample[1:d.cx]
-    row_z = syndrome_sample[d.cx+1:d.cx+d.cz]
+    row_x = @view syndrome_sample[1:d.cx]
+    row_z = @view syndrome_sample[d.cx+1:d.cx+d.cz]
     guess_z, success = LDPCDecoders.decode!(d.bpdecoderx, row_x)
     guess_x, success = LDPCDecoders.decode!(d.bpdecoderz, row_z)
     return vcat(guess_x, guess_z)
 end
 
 function decode(d::BitFlipDecoder, syndrome_sample)
-    row_x = syndrome_sample[1:d.cx]
-    row_z = syndrome_sample[d.cx+1:d.cx+d.cz]
+    row_x = @view syndrome_sample[1:d.cx]
+    row_z = @view syndrome_sample[d.cx+1:d.cx+d.cz]
     guess_z, success = LDPCDecoders.decode!(d.bfdecoderx, row_x)
     guess_x, success = LDPCDecoders.decode!(d.bfdecoderz, row_z)
     return vcat(guess_x, guess_z)

ext/QuantumCliffordPyQDecodersExt/QuantumCliffordPyQDecodersExt.jl

@@ -69,8 +69,8 @@ end
 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]
+    row_x = @view syndrome_sample[1:d.nx]
+    row_z = @view syndrome_sample[d.nx+1:end]
     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)
@@ -106,18 +106,19 @@ end
 parity_checks(d::PyMatchingDecoder) = d.H
 
 function decode(d::PyMatchingDecoder, syndrome_sample)
-    row_x = syndrome_sample[1:d.nx] # TODO This copy is costly!
-    row_z = syndrome_sample[d.nx+1:end]
+    row_x = @view syndrome_sample[1:d.nx]
+    row_z = @view syndrome_sample[d.nx+1:end]
     guess_z_errors = PythonCall.PyArray(d.pyx.decode(row_x))
     guess_x_errors = PythonCall.PyArray(d.pyz.decode(row_z))
     vcat(guess_x_errors, guess_z_errors)
 end
 
 function batchdecode(d::PyMatchingDecoder, syndrome_samples)
-    row_x = syndrome_samples[:,1:d.nx] # TODO This copy is costly!
-    row_z = syndrome_samples[:,d.nx+1:end]
+    row_x = @view syndrome_samples[:,1:d.nx]
+    row_z = @view syndrome_samples[:,d.nx+1:end]
     guess_z_errors = PythonCall.PyArray(d.pyx.decode_batch(row_x))
     guess_x_errors = PythonCall.PyArray(d.pyz.decode_batch(row_z))
+    n_cols_x = size(guess_x_errors, 2)
     hcat(guess_x_errors, guess_z_errors)
 end
 

src/QuantumClifford.jl

@@ -963,14 +963,61 @@ function check_allrowscommute(stabilizer::Tableau)
 end
 check_allrowscommute(stabilizer::Stabilizer)=check_allrowscommute(tab(stabilizer))
 
+"""
+Vertically concatenates tableaux.
+
+```jldoctest
+julia> vcat(ghz(2), ghz(2))
++ XX
++ ZZ
++ XX
++ ZZ
+```
+
+See also: [`hcat`](@ref)
+"""
 function Base.vcat(tabs::Tableau...)
     Tableau(
         vcat((s.phases for s in tabs)...),
         tabs[1].nqubits,
         hcat((s.xzs for s in tabs)...))
 end
 
-Base.vcat(stabs::Stabilizer...) = Stabilizer(vcat((tab(s) for s in stabs)...))
+Base.vcat(stabs::Stabilizer{T}...) where {T} = Stabilizer(vcat((tab(s) for s in stabs)...))
+
+"""
+Horizontally concatenates tableaux.
+
+```jldoctest
+julia> hcat(ghz(2), ghz(2))
++ XXXX
++ ZZZZ
+```
+
+See also: [`vcat`](@ref)
+"""
+function Base.hcat(tabs::Tableau...) # TODO this implementation is slow as it unpacks each bitvector into bits and repacks them -- reuse the various tableau inset functionality we have to speed this up
+    rows = size(tabs[1], 1)
+    cols = sum(map(nqubits, tabs))
+    newtab = zero(Tableau, rows, cols)
+    cols_idx = 1
+    for tab in tabs
+        rows_tab, cols_tab = size(tab)
+        if rows_tab != rows
+            throw(ArgumentError("All input Tableaux/Stabilizers must have the same number of rows."))
+        end
+        for i in 1:rows
+            for j in 1:cols_tab
+                newtab[i, cols_idx+j-1]::Tuple{Bool,Bool} = tab[i, j]::Tuple{Bool,Bool}
+            end
+            newtab.phases[i] = (newtab.phases[i]+tab.phases[i])%4
+        end
+        cols_idx += cols_tab
+    end
+    return newtab
+end
+
+Base.hcat(stabs::Stabilizer{T}...) where {T} = Stabilizer(hcat((tab(s) for s in stabs)...))
 
 ##############################
 # Unitary Clifford Operations
@@ -1017,6 +1064,16 @@ function _apply!(stab::AbstractStabilizer, p::PauliOperator, indices; phases::Va
     stab
 end
 
+##############################
+# Conversion and promotion
+##############################
+
+Base.promote_rule(::Type{<:Destabilizer{T}}   , ::Type{<:MixedDestabilizer{T}}) where {T<:Tableau} = MixedDestabilizer{T}
+Base.promote_rule(::Type{<:MixedStabilizer{T}}, ::Type{<:MixedDestabilizer{T}}) where {T<:Tableau} = MixedDestabilizer{T}
+Base.promote_rule(::Type{<:Stabilizer{T}}     , ::Type{<:S}                   ) where {T<:Tableau, S<:Union{MixedStabilizer{T}, Destabilizer{T}, MixedDestabilizer{T}}} = S
+
+Base.convert(::Type{<:MixedDestabilizer{T}}, x::Union{Destabilizer{T}, MixedStabilizer{T}, Stabilizer{T}}) where {T <: Tableau} = MixedDestabilizer(x)
+
 ##############################
 # Helpers for binary codes
 ##############################

src/ecc/ECC.jl

@@ -20,7 +20,7 @@ export parity_checks, parity_checks_x, parity_checks_z, iscss,
     RepCode, LiftedCode,
     CSS,
     Shor9, Steane7, Cleve8, Perfect5, Bitflip3,
-    Toric, Gottesman, Surface, Concat, CircuitCode,
+    Toric, Gottesman, Surface, Concat, CircuitCode, QuantumReedMuller,
     LPCode, two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes,
     random_brickwork_circuit_code, random_all_to_all_circuit_code,
     evaluate_decoder,
@@ -376,10 +376,10 @@ include("codes/gottesman.jl")
 include("codes/surface.jl")
 include("codes/concat.jl")
 include("codes/random_circuit.jl")
-
 include("codes/classical/reedmuller.jl")
 include("codes/classical/recursivereedmuller.jl")
 include("codes/classical/bch.jl")
+include("codes/quantumreedmuller.jl")
 
 # qLDPC
 include("codes/classical/lifted.jl")