Bug fix to the parity_checks(ReedMuller(r, m)) along with RecursiveReedMuller code for cross-reference

src/ecc/ECC.jl

@@ -1,6 +1,7 @@
 module ECC
 
 using LinearAlgebra
+using LinearAlgebra: I
 using QuantumClifford
 using QuantumClifford: AbstractOperation, AbstractStabilizer, Stabilizer
 import QuantumClifford: Stabilizer, MixedDestabilizer, nqubits
@@ -69,6 +70,8 @@ In a [polynomial code](https://en.wikipedia.org/wiki/Polynomial_code), the gener
 """
 function generator_polynomial end
 
+function generator end
+
 parity_checks(s::Stabilizer) = s
 Stabilizer(c::AbstractECC) = parity_checks(c)
 MixedDestabilizer(c::AbstractECC; kwarg...) = MixedDestabilizer(Stabilizer(c); kwarg...)
@@ -361,4 +364,5 @@ include("codes/surface.jl")
 include("codes/concat.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,59 @@
+"""
+The Plotkin 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.
+"""
+abstract type ClassicalCode end
+
+struct RecursiveReedMuller <: ClassicalCode
+    r::Int
+    m::Int
+
+    function RecursiveReedMuller(r, m)
+        0 ≤ r ≤ m || throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m in order to obtain a valid code."))
+        new(r, m)
+    end
+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)
+
+"""
+This function generates the generator matrix, `G`, for RecursiveReedMuller`(RecursiveReedMuller(r, m))` error-correcting codes. 
+
+`generator(RecursiveReedMuller(r, m))`:
+- `m`: Positive integer representing the message length.
+- `r`: The order of the Reed-Muller code. Must satisfy `0 ≤ r ≤ m`.
+"""
+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

src/ecc/codes/classical/reedmuller.jl

@@ -1,8 +1,14 @@
 """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 exist:
+    1. `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. `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 ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/reed_muller).
 """
 
 abstract type ClassicalCode end
@@ -12,27 +18,41 @@ struct ReedMuller <: ClassicalCode
     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 || m < 1 || m >= 11
+            throw(ArgumentError("Invalid parameters: r must be non-negative and r ≤ m. Additionally, m must be positive and < 11 in order to obtain a valid code and to remain tractable"))
         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)
+"""
+This function generates the generator matrix, `G`, for Reed-Muller `(RM(r, m))` error-correcting codes. 
+
+`generator(ReedMuller(r, m))`:
+- `m`: Positive integer representing the message length.
+- `r`: Nonnegative integer less than or equal to `m`, specifying the code's order.
+"""
+function generator(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)]
+    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
+
+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)

test/Project.toml

@@ -20,6 +20,7 @@ Quantikz = "b0d11df0-eea3-4d79-b4a5-421488cbf74b"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"

test/test_ecc_reedmuller.jl

@@ -1,68 +1,74 @@
 using Test
 using Nemo
-using Combinatorics
 using LinearAlgebra
+using RowEchelon: rref
 using QuantumClifford
 using QuantumClifford.ECC
-using QuantumClifford.ECC: AbstractECC, ReedMuller
+using QuantumClifford.ECC: AbstractECC, ReedMuller, generator, RecursiveReedMuller 
 
-function binomial_coeff_sum(r, m)
-    total = 0
-    for i in 0:r
-        total += length(combinations(1:m, i))
+function designed_distance(matrix, m, r)
+    for row in eachrow(matrix)
+        count = sum(row)
+        if count >= 2 ^ (m - r)
+            return true
+        end
     end
-    return total
+    return false
 end
 
-@testset "Test RM(r, m) Matrix Rank" begin
-    for m in 2:5
-        for r in 0:m - 1
-            H = parity_checks(ReedMuller(r, m))
+@testset "Test RM(r, m) properties" begin
+    for m in 3:10
+        for r in 0:m
+            H = generator(ReedMuller(r, m))
             mat = Nemo.matrix(Nemo.GF(2), H)
             computed_rank = LinearAlgebra.rank(mat)
-            expected_rank = binomial_coeff_sum(r, m)
+            expected_rank = sum(binomial.(m, 0:r))
             @test computed_rank == expected_rank
+            @test designed_distance(H, m, r) == true
+            @test rate(ReedMuller(r, m)) == sum(binomial.(m, 0:r)) / 2 ^ m
+            @test code_n(ReedMuller(r, m)) == 2 ^ m
+            @test code_k(ReedMuller(r, m)) == sum(binomial.(m, 0:r))
+            @test distance(ReedMuller(r, m)) == 2 ^ (m - r) 
+            @test rref(generator(ReedMuller(r, m))) == rref(generator(RecursiveReedMuller(r, m))) 
         end
     end
-end
-
-@testset "Testing common examples of RM(r,m) codes [raaphorst2003reed](@cite), [djordjevic2021quantum](@cite), [abbe2020reed](@cite)" begin
-
-    #RM(0,3)  
-    @test parity_checks(ReedMuller(0,3)) == [1 1 1 1 1 1 1 1]
+
+    # Testing common examples of RM(r,m) codes [raaphorst2003reed](@cite), [djordjevic2021quantum](@cite), [abbe2020reed](@cite).
+    # RM(0,3)  
+    @test generator(ReedMuller(0,3)) == [1 1 1 1 1 1 1 1]
 
     #RM(1,3) 
-    @test parity_checks(ReedMuller(1,3)) == [1 1 1 1 1 1 1 1;
-                                             1 1 1 1 0 0 0 0;
-                                             1 1 0 0 1 1 0 0;
-                                             1 0 1 0 1 0 1 0]
+    @test generator(ReedMuller(1,3)) == [1 1 1 1 1 1 1 1;
+                                         1 1 1 1 0 0 0 0;
+                                         1 1 0 0 1 1 0 0;
+                                         1 0 1 0 1 0 1 0]
     #RM(2,3)
-    @test parity_checks(ReedMuller(2,3)) == [1 1 1 1 1 1 1 1;
-                                             1 1 1 1 0 0 0 0;
-                                             1 1 0 0 1 1 0 0;
-                                             1 0 1 0 1 0 1 0;
-                                             1 1 0 0 0 0 0 0;
-                                             1 0 1 0 0 0 0 0;
-                                             1 0 0 0 1 0 0 0]
+    @test generator(ReedMuller(2,3)) == [1 1 1 1 1 1 1 1;
+                                         1 1 1 1 0 0 0 0;
+                                         1 1 0 0 1 1 0 0;
+                                         1 0 1 0 1 0 1 0;
+                                         1 1 0 0 0 0 0 0;
+                                         1 0 1 0 0 0 0 0;
+                                         1 0 0 0 1 0 0 0]
     #RM(3,3)
-    @test parity_checks(ReedMuller(3,3)) == [1 1 1 1 1 1 1 1;
-                                             1 1 1 1 0 0 0 0;
-                                             1 1 0 0 1 1 0 0;
-                                             1 0 1 0 1 0 1 0;
-                                             1 1 0 0 0 0 0 0;
-                                             1 0 1 0 0 0 0 0;
-                                             1 0 0 0 1 0 0 0;
-                                             1 0 0 0 0 0 0 0]
+    @test generator(ReedMuller(3,3)) == [1 1 1 1 1 1 1 1;
+                                         1 1 1 1 0 0 0 0;
+                                         1 1 0 0 1 1 0 0;
+                                         1 0 1 0 1 0 1 0;
+                                         1 1 0 0 0 0 0 0;
+                                         1 0 1 0 0 0 0 0;
+                                         1 0 0 0 1 0 0 0;
+                                         1 0 0 0 0 0 0 0]
     #RM(2,4)
-    @test parity_checks(ReedMuller(2,4)) == [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
-                                             1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0;
-                                             1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0;
-                                             1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0;
-                                             1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0;
-                                             1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
-                                             1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
-                                             1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0;
-                                             1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0;
-                                             1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0;
-                                             1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0]
+    @test generator(ReedMuller(2,4)) == [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
+                                         1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0;
+                                         1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0;
+                                         1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0;
+                                         1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0;
+                                         1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
+                                         1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
+                                         1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0;
+                                         1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0;
+                                         1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0;
+                                         1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0]
 end

test/test_ecc_throws.jl

@@ -4,6 +4,7 @@ using QuantumClifford.ECC: ReedMuller, BCH
 
 @test_throws ArgumentError ReedMuller(-1, 3)
 @test_throws ArgumentError ReedMuller(1, 0)
+@test_throws ArgumentError ReedMuller(4, 2)
 
 @test_throws ArgumentError BCH(2, 2)
 @test_throws ArgumentError BCH(3, 4)