Merge branch 'master' into i208

CHANGELOG.md

@@ -5,6 +5,11 @@
 
 # News
 
+## dev
+
+- **(fix)** Bug fix to the `parity_checks(ReedMuller(r, m))` of classical Reed-Muller code (it was returning generator matrix).
+- `RecursiveReedMuller` code implementation as an alternative implementation of `ReedMuller`.
+
 ## v0.9.9 - 2024-08-05
 
 - `inv` is implemented for all Clifford operator types (symbolic, dense, sparse).

src/QuantumClifford.jl

@@ -767,10 +767,19 @@ function comm!(v, l::PauliOperator, r::Tableau)
     v
 end
 comm!(v, l::Tableau, r::PauliOperator) = comm!(v, r, l)
-@inline comm!(v, l::PauliOperator, r::Stabilizer, i::Int) = comm!(v, l, tab(r), i)
-@inline comm!(v, l::Stabilizer, r::PauliOperator, i::Int) = comm!(v, tab(l), r, i)
 @inline comm!(v, l::PauliOperator, r::Stabilizer) = comm!(v, l, tab(r))
 @inline comm!(v, l::Stabilizer, r::PauliOperator) = comm!(v, tab(l), r)
+function comm!(v, l::PauliOperator, r::Tableau, i)
+    v[i] = comm(l,r,i)
+    v
+end
+comm!(v, l::Tableau, r::PauliOperator, i) = comm!(v, r, l, i)
+@inline comm!(v, l::PauliOperator, r::Stabilizer, i::Int) = comm!(v, l, tab(r), i)
+@inline comm!(v, l::Stabilizer, r::PauliOperator, i::Int) = comm!(v, tab(l), r, i)
+function comm!(v, s::Tableau, l::Int, r::Int)
+    v[l] = comm(s, l, r)
+    v
+end
 @inline comm!(v, s::Stabilizer, l::Int, r::Int) = comm!(v, tab(s), l, r)
 
 

src/dense_cliffords.jl

@@ -82,7 +82,7 @@ function row_limit(str, limit=50)
 end
 
 digits_subchars = collect("₀₁₂₃₄₅₆₇₈₉")
-digits_substr(n,nwidth) = join(([digits_subchars[d+1] for d in reverse(digits(n, pad=nwidth))]))
+digits_substr(n::Int,nwidth::Int) = join(([digits_subchars[d+1] for d in reverse(digits(n, pad=nwidth))]))
 
 function Base.copy(c::CliffordOperator)
     CliffordOperator(copy(c.tab))

src/ecc/ECC.jl

@@ -1,14 +1,15 @@
 module ECC
 
 using LinearAlgebra
+using LinearAlgebra: I
 using QuantumClifford
 using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer
 import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
 using DocStringExtensions
 using Combinatorics: combinations
 using SparseArrays: sparse
 using Statistics: std
-using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible
+using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqPolyRingElem, FqFieldElem, is_zero, degree, defining_polynomial, is_irreducible, echelon_form
 
 abstract type AbstractECC end
 
@@ -70,6 +71,9 @@ In a [polynomial code](https://en.wikipedia.org/wiki/Polynomial_code), the gener
 """
 function generator_polynomial end
 
+"""The generator matrix of a code."""
+function generator end
+
 parity_checks(s::Stabilizer) = s
 Stabilizer(c::AbstractECC) = parity_checks(c)
 MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); kwarg...)
@@ -335,8 +339,9 @@ isdegenerate(c::AbstractECC, args...) = isdegenerate(parity_checks(c), args...)
 isdegenerate(c::AbstractStabilizer, args...) = isdegenerate(stabilizerview(c), args...)
 
 function isdegenerate(H::Stabilizer, errors) # Described in https://quantumcomputing.stackexchange.com/questions/27279
-    syndromes = comm.((H,), errors) # TODO This can be optimized by having something that always returns bitvectors
-    return length(Set(syndromes)) != length(errors)
+    syndromes = map(e -> comm(H,e), errors) # TODO This can be optimized by having something that always returns bitvectors
+    syndrome_set = Set(syndromes)
+    return length(syndrome_set) != length(errors)
 end
 
 function isdegenerate(H::Stabilizer, d::Int=1)
@@ -363,4 +368,5 @@ include("codes/concat.jl")
 include("codes/random_circuit.jl")
 include("codes/classical/reedmuller.jl")
 include("codes/classical/bch.jl")
+include("codes/classical/recursivereedmuller.jl")
 end #module

src/ecc/codes/classical/recursivereedmuller.jl

@@ -0,0 +1,64 @@
+"""
+A construction of the Reed-Muller class of codes using the recursive definition.
+
+The Plotkin `(u, u + v)` construction defines a recursive relation between generator matrices of Reed-Muller `(RM)` codes [abbe2020reed](@cite). To derive the generator matrix `G(m, r)` for `RM(r, m)`, the generator matrices of lower-order codes are utilized:
+- `G(r - 1, m - 1)`: Generator matrix of `RM(r - 1, m - 1)`
+- `G(r, m - 1)`: Generator matrix of `RM(r, m - 1)`
+
+The generator matrix `G(m, r)` of `RM(m, r)` is formulated as follows in matrix notation:
+
+```math
+G(m, r) = \begin{bmatrix}
+G(r, m - 1) & G(r, m - 1) \\
+0 & G(r - 1, m - 1)
+\end{bmatrix}
+```
+
+Here, the matrix 0 denotes an all-zero matrix with dimensions matching `G(r - 1, m - 1)`. This recursive approach facilitates the construction of higher-order Reed-Muller codes based on the generator matrices of lower-order codes.
+
+In addition, the dimension of `RM(m - r - 1, m)` equals the dimension of the dual of `RM(r, m)`. Thus, `RM(m - r - 1, m) = RM(r, m)^⊥` shows that the [dual code](https://en.wikipedia.org/wiki/Dual_code) of `RM(r, m)` is `RM(m − r − 1, m)`, indicating the parity check matrix of `RM(r, m)` is the generator matrix for `RM(m - r - 1, m)`.
+
+See also: `ReedMuller`
+"""
+struct RecursiveReedMuller <: ClassicalCode
+    r::Int
+    m::Int
+
+    function RecursiveReedMuller(r, m)
+        if r < 0 || r > m
+            throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to valid code."))
+        end
+        new(r, m)
+    end
+end
+
+function _recursiveReedMuller(r::Int, m::Int)
+    if r == 1 && m == 1
+        return Matrix{Int}([1 1; 0 1])
+    elseif r == m
+        return Matrix{Int}(I, 2^m, 2^m)
+    elseif r == 0
+        return Matrix{Int}(ones(1, 2^m))
+    else
+        Gᵣₘ₋₁ = _recursiveReedMuller(r, m - 1)
+        Gᵣ₋₁ₘ₋₁ = _recursiveReedMuller(r - 1, m - 1)
+        return vcat(hcat(Gᵣₘ₋₁, Gᵣₘ₋₁), hcat(zeros(Int, size(Gᵣ₋₁ₘ₋₁)...), Gᵣ₋₁ₘ₋₁))
+    end
+end
+
+function generator(c::RecursiveReedMuller)
+    return _recursiveReedMuller(c.r, c.m)
+end
+
+function parity_checks(c::RecursiveReedMuller)
+    H = generator(RecursiveReedMuller(c.m - c.r - 1, c.m))
+    return H
+end
+
+code_n(c::RecursiveReedMuller) = 2 ^ c.m
+
+code_k(c::RecursiveReedMuller) = sum(binomial.(c.m, 0:c.r))
+
+distance(c::RecursiveReedMuller) = 2 ^ (c.m - c.r)
+
+rate(c::RecursiveReedMuller) = code_k(c::RecursiveReedMuller) / code_n(c::RecursiveReedMuller)

src/ecc/codes/classical/reedmuller.jl

@@ -1,38 +1,60 @@
 """The family of Reed-Muller codes, as discovered by Muller in his 1954 paper [muller1954application](@cite) and Reed who proposed the first efficient decoding algorithm [reed1954class](@cite).
 
+Let `m` be a positive integer and `r` a nonnegative integer with `r ≤ m`. These linear codes, denoted as `RM(r, m)`, have order `r` (where `0 ≤ r ≤ m`) and codeword length `n` of `2ᵐ`.
+
+Two special cases of generator(RM(r, m)) exist:
+    1. `generator(RM(0, m))`: This is the `0ᵗʰ`-order `RM` code, similar to the binary repetition code with length `2ᵐ`. It's characterized by a single basis vector containing all ones.
+    2. `generator(RM(m, m))`: This is the `mᵗʰ`-order `RM` code. It encompasses the entire field `F(2ᵐ)`, representing all possible binary strings of length `2ᵐ`.
+
 You might be interested in consulting [raaphorst2003reed](@cite), [abbe2020reed](@cite), and [djordjevic2021quantum](@cite) as well.
 
-The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller)
-"""
+The dimension of `RM(m - r - 1, m)` equals the dimension of the dual of `RM(r, m)`. Thus, `RM(m - r - 1, m) = RM(r, m)^⊥` shows that the [dual code](https://en.wikipedia.org/wiki/Dual_code) of `RM(r, m)` is `RM(m − r − 1, m)`, indicating the parity check matrix of `RM(r, m)` is the generator matrix for `RM(m - r - 1, m)`.
 
-abstract type ClassicalCode end
+The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller).
 
+See also: `RecursiveReedMuller`
+"""
 struct ReedMuller <: ClassicalCode
     r::Int
     m::Int
 
     function ReedMuller(r, m)
-        if r < 0 || m < 1 || m >= 11
-            throw(ArgumentError("Invalid parameters: r must be non-negative and m must be positive and < 11 in order to obtain a valid code and to remain tractable"))
+        if r < 0 || r > m
+            throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to valid code."))
         end
         new(r, m)
     end
 end
 
-function variables_xi(m, i)
-    return repeat([fill(1, 2^(m - i - 1)); fill(0, 2^(m - i - 1))], outer = 2^i)
+function _variablesₓᵢ_rm(m, i)
+    return repeat([fill(1, 2 ^ (m - i - 1)); fill(0, 2 ^ (m - i - 1))], outer = 2 ^ i)
 end
 
-function vmult(vecs...)
+function _vmult_rm(vecs...)
     return [reduce(*, a, init=1) for a in zip(vecs...)]
 end
 
-function parity_checks(c::ReedMuller)
-    r=c.r
-    m=c.m
-    xi = [variables_xi(m, i) for i in 0:m - 1]
-    row_matrices = [reduce(vmult, [xi[i + 1] for i in S], init = ones(Int, 2^m)) for s in 0:r for S in combinations(0:m - 1, s)]
+function generator(c::ReedMuller)
+    r = c.r
+    m = c.m
+    xᵢ = [_variablesₓᵢ_rm(m, i) for i in 0:m - 1]
+    row_matrices = [reduce(_vmult_rm, [xᵢ[i + 1] for i in S], init = ones(Int, 2 ^ m)) for s in 0:r for S in combinations(0:m - 1, s)]
     rows = length(row_matrices)
     cols = length(row_matrices[1])
-    H = reshape(vcat(row_matrices...), cols, rows)'
-end
+    G = reshape(vcat(row_matrices...), cols, rows)'
+    G = Matrix{Bool}(G)
+    return G 
+end
+
+function parity_checks(c::ReedMuller)
+    H = generator(ReedMuller(c.m - c.r - 1, c.m))
+    return H
+end
+
+code_n(c::ReedMuller) = 2 ^ c.m
+
+code_k(c::ReedMuller) = sum(binomial.(c.m, 0:c.r))
+
+distance(c::ReedMuller) = 2 ^ (c.m - c.r)
+
+rate(c::ReedMuller) = code_k(c::ReedMuller) / code_n(c::ReedMuller)

src/ecc/decoder_pipeline.jl

@@ -132,7 +132,7 @@ function evaluate_decoder(d::AbstractSyndromeDecoder, setup::AbstractECCSetup, n
 
     physical_noisy_circ, syndrome_bits, n_anc = physical_ECC_circuit(H, setup)
     encoding_circ = naive_encoding_circuit(H)
-    preX = [sHadamard(i) for i in n-k+1:n]
+    preX = sHadamard[sHadamard(i) for i in n-k+1:n]
 
     mdH = MixedDestabilizer(H)
     logX_circ, _, logX_bits = naive_syndrome_circuit(logicalxview(mdH), n_anc+1, last(syndrome_bits)+1)

src/misc_ops.jl

@@ -29,11 +29,21 @@ end
 struct SparseGate{T<:Tableau} <: AbstractCliffordOperator # TODO simplify type parameters (remove nesting)
     cliff::CliffordOperator{T}
     indices::Vector{Int}
+    function SparseGate(cliff::CliffordOperator{T}, indices::Vector{Int}) where T<:Tableau
+        if length(indices) != nqubits(cliff)
+            throw(ArgumentError("The number of target qubits (indices) must match the qubit count in the CliffordOperator."))
+        end
+        new{T}(cliff, indices)
+    end
 end
 
 SparseGate(c,t::Tuple) = SparseGate(c,collect(t))
 
 function apply!(state::AbstractStabilizer, g::SparseGate; kwargs...)
+    m = maximum(g.indices)
+    if m > nqubits(state)
+        throw(ArgumentError(lazy"SparseGate was attempted on invalid qubit index $(m) when the state contains only $(nqubits(state)) qubits."))
+    end
     apply!(state, g.cliff, g.indices; kwargs...)
 end