Bicycle and Unicycle codes: Fixing doc and method errors

docs/src/references.bib

@@ -560,4 +560,4 @@ @article{haah2011local
   number={4},
   pages={042330},
   year={2011},
-}
+}

src/ecc/ECC.jl

@@ -27,6 +27,7 @@ export parity_checks, parity_checks_x, parity_checks_z, iscss,
     RepCode, LiftedCode,
     CSS,
     Shor9, Steane7, Cleve8, Perfect5, Bitflip3,
+    Unicycle, Bicycle,
     Toric, Gottesman, Surface, Concat, CircuitCode, QuantumReedMuller,
     LPCode, two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes,
     haah_cubic_codes,
@@ -392,5 +393,6 @@ include("codes/quantumreedmuller.jl")
 # qLDPC
 include("codes/classical/lifted.jl")
 include("codes/lifted_product.jl")
+include("codes/simple_sparse_codes.jl")
 
 end #module

src/ecc/codes/simple_sparse_codes.jl

@@ -0,0 +1,284 @@
+"""Struct for an LDPC bicycle code based on [mackay2004sparse](@cite). This code is produced by taking a circulant matrix based on a difference set, reducing the number of rows, and then using it as the base matrix for a CSS code.
+
+Parameters:
+$TYPEDFIELDS
+
+```jldoctest Bicycle
+julia> Bicycle(4, 2).n
+4
+
+julia> Bicycle(4, 2).m
+2
+
+julia> parity_checks(Bicycle(4, 2))
++ XXXX
++ ZZZZ
+
+julia> parity_checks(Bicycle(6, 4))
++ XX_X_X
++ X_X_XX
++ ZZ_Z_Z
++ Z_Z_ZZ
+```
+"""
+struct Bicycle
+    "Width of array, should be >= 2. Equal to number to number of physical qubits (N)"
+    n::Int
+    "Height of array, should be >= 2 and a multiple of 2. Equal to number of parity checks (M)"
+    m::Int
+end
+
+function parity_checks(b::Bicycle)
+    m = b.m
+    n = b.n
+    if m%2 == 1
+        throw(DomainError(m, " M should be a multiple for 2 for bicycle codes."))
+    end
+    if m < 2
+        throw(DomainError(m, " M is too small, make it greater than 1."))
+    end
+    if n%2 == 1
+        throw(DomainError(n, " N should be a multiple for 2 for bicycle codes."))
+    end
+    if n < 2
+        throw(DomainError(n, " N is too small, make it greater than 1."))
+    end
+    bs = bicycle_set_gen(Int(n/2))
+    bsc = circ_to_bicycle_h0(bs, Int(n/2))
+    while size(bsc)[1] > m/2
+        print(bsc)
+        bsc = reduce_bicycle(bsc)
+    end
+    return parity_checks(CSS(bsc, bsc))
+end
+
+"""Struct for an LDPC unicycle code based on [mackay2004sparse](@cite). The parity check matrix is produced by taking a perfect difference set, turning it into a circulant matrix, and removing the linearly dependent rows. Good performance when decoded by message passing decoders.
+
+Has the drawback of only being possible to make at sizes that correspond to existing perfect difference sets. Reference [mackay2004sparse](@cite) to find good sizes of perfect difference sets to use.
+
+Parameters:
+$TYPEDFIELDS
+
+```jldoctest Unicycle
+julia> typeof(Unicycle(21, [1, 3, 8, 9, 12]))
+Unicycle
+
+julia> parity_checks(Unicycle(21, [1, 3, 8, 9, 12]))
++ _X_________X_X____XX_X
++ __X_________X_X____XXX
++ X__X_________X_X____XX
++ XX__X_________X_X____X
+ ⋮
++ ____ZZ__Z_________Z_ZZ
++ Z____ZZ__Z_________Z_Z
++ _Z____ZZ__Z_________ZZ
+```
+"""
+struct Unicycle
+    "Size of parity check matrix, also max size of generated set array, should be >= 1. Equal to the number of physical qubits."
+    N::Int
+    "Perfect difference set modulo (N). Array of indices that are 'active' checks in the circulant code."
+    set::Vector{Int}
+end
+
+function parity_checks(u::Unicycle)
+    n = u.N
+    set = u.set
+    usc = circ_to_unicycle_h0(set, n)
+    rusc = reduce_unicycle(usc) # reduced unicycle code
+    return parity_checks(CSS(rusc, rusc))
+end
+
+"""Takes an untrimmed bicycle matrix and removes the row which keeps the spread of the column weights minimal.
+
+Required before the bicycle code can be used.
+
+Typical usage:
+```jldoctest reduce_bicycle
+
+julia> reduce_bicycle(circ_to_bicycle_h0([1, 2, 4], 7))
+6×14 Matrix{Bool}:
+ 0  0  1  1  0  1  0  1  0  0  0  1  0  1
+ 0  0  0  1  1  0  1  1  1  0  0  0  1  0
+ 1  0  0  0  1  1  0  0  1  1  0  0  0  1
+ 0  1  0  0  0  1  1  1  0  1  1  0  0  0
+ 1  0  1  0  0  0  1  0  1  0  1  1  0  0
+ 1  1  0  1  0  0  0  0  0  1  0  1  1  0
+
+```
+"""
+function reduce_bicycle(H0::Matrix{Bool})
+    m, n = size(H0)
+    r_i = 0
+    std_min = Inf
+    for i in 1:m
+        t_H0 = vcat(H0[1:i-1, :], H0[i+1:end, :])
+        std_temp = std(convert(Array, sum(t_H0, dims = 1)))
+        if std_temp < std_min
+            std_min = std_temp
+            r_i = i
+        end
+    end
+    return vcat(H0[1:r_i-1, :], H0[r_i+1:end, :])
+end
+
+"""Takes a list of indices and creates the base of the bicycle matrix.
+
+For example:
+
+```jldoctest circ_to_bicycle_h0
+
+julia> circ_to_bicycle_h0([1, 2, 4], 7)
+7×14 Matrix{Bool}:
+ 0  1  1  0  1  0  0  0  0  0  1  0  1  1
+ 0  0  1  1  0  1  0  1  0  0  0  1  0  1
+ 0  0  0  1  1  0  1  1  1  0  0  0  1  0
+ 1  0  0  0  1  1  0  0  1  1  0  0  0  1
+ 0  1  0  0  0  1  1  1  0  1  1  0  0  0
+ 1  0  1  0  0  0  1  0  1  0  1  1  0  0
+ 1  1  0  1  0  0  0  0  0  1  0  1  1  0
+```
+
+See [mackay2004sparse](@cite) for more details
+"""
+function circ_to_bicycle_h0(circ_indices, n::Int)
+    circ_arr = Vector{Bool}(undef, n)
+    circ_matrix = Matrix{Bool}(undef, n, n)
+    comp_matrix = Matrix{Bool}(undef, n, 2*n)
+    for i = 1:n
+        if Int(i-1) in circ_indices
+            circ_arr[i] = true
+        else
+            circ_arr[i] = false
+        end
+    end
+    for i = 1:n
+        circ_matrix[i,1:n] = circ_arr
+        li = circ_arr[end]
+        circ_arr[2:end] = circ_arr[1:end-1]
+        circ_arr[1] = li
+    end
+    comp_matrix[1:n,1:n] = circ_matrix
+    comp_matrix[1:n,n+1:2*n] = transpose(circ_matrix)
+    return comp_matrix
+end
+
+"""Takes an untrimmed unicycle matrix and removes linearly dependent rows.
+
+Required before the unicycle code can be used.
+
+Typical usage:
+
+```jldoctest reduce_unicycle
+julia> reduce_unicycle(circ_to_unicycle_h0([1, 2, 4], 7))
+4×8 Matrix{Bool}:
+ 0  0  0  1  1  0  1  1
+ 1  0  0  0  1  1  0  1
+ 0  1  0  0  0  1  1  1
+ 1  0  1  0  0  0  1  1
+
+```
+
+"""
+function reduce_unicycle(m::Matrix{Bool})
+    rrzz, = residue_ring(ZZ, 2)
+    nm = matrix(rrzz, m)
+    r = LinearAlgebra.rank(nm)
+    i = 1
+    while size(m)[1] > r
+        tm = vcat(m[1:i-1,:], m[i+1:end,:])
+        tr = LinearAlgebra.rank(matrix(rrzz, tm))
+        if(tr == r)
+            m = tm
+            if(size(m)[1] == r)
+                break
+            end
+        else
+            i += 1
+        end
+    end
+    return m
+end
+
+"""Takes a list of indices and creates the base of the unicycle matrix.
+
+For example:
+
+```jldoctest circ_to_unicycle_h0
+julia> circ_to_unicycle_h0([1, 2, 4], 7)
+7×8 Matrix{Bool}:
+ 1  1  0  1  0  0  0  1
+ 0  1  1  0  1  0  0  1
+ 0  0  1  1  0  1  0  1
+ 0  0  0  1  1  0  1  1
+ 1  0  0  0  1  1  0  1
+ 0  1  0  0  0  1  1  1
+ 1  0  1  0  0  0  1  1
+
+```
+See [mackay2004sparse](@cite) for more details"""
+function circ_to_unicycle_h0(circ_indices::Vector{Int}, n::Int)
+    circ_arr = fill(false, n)
+    one_col = transpose(fill(true, n))
+    circ_matrix = Matrix{Bool}(undef, n, n)
+    comp_matrix = Matrix{Bool}(undef, n, n+1)
+    for i = 1:n
+        if i in circ_indices
+            circ_arr[i] = true
+        else
+            circ_arr[i] = false
+        end
+    end
+    for i = 1:n
+        circ_matrix[i,1:n] = circ_arr
+        li = circ_arr[end]
+        circ_arr[2:end] = circ_arr[1:end-1]
+        circ_arr[1] = li
+    end
+    comp_matrix[1:n,1:n] = circ_matrix
+    comp_matrix[1:n,n+1] = one_col
+    return comp_matrix
+end
+
+"""Attempts to generate a list of indices to be used in a bicycle code using a search method given a matrix width (N)"""
+function bicycle_set_gen(N::Int)
+    circ_arr = Int[0]
+    diff_arr = Int[]
+    circ_arr[1] = 0
+    # test new elements
+    for add_i = (circ_arr[end] + 1):N - 1
+        valid = true
+        temp_circ_arr = copy(circ_arr)
+        temp_diff_arr = []
+        push!(temp_circ_arr, add_i)
+        for j = 1:size(temp_circ_arr)[1]
+            temp_arr = copy(temp_circ_arr)
+            # add lesser elements + N to temp_arr
+            for k = 1:size(temp_circ_arr)[1]
+                if k < j
+                    push!(temp_arr, temp_circ_arr[k] + N)
+                else
+                    break
+                end
+            end
+            # test if new index is valid
+            for k = 1:(size(temp_circ_arr)[1] - 2)
+                t_diff = (temp_arr[j + k] - temp_arr[j]) % N
+                if ((t_diff) in temp_diff_arr)
+                    valid = false
+                    break
+                else
+                    push!(temp_diff_arr, t_diff)
+                end
+            end
+            if !valid
+                break
+            end
+        end
+        if valid
+            circ_arr = copy(temp_circ_arr)
+            diff_arr = copy(temp_diff_arr)
+        end
+    end
+    return circ_arr
+end

test/test_ecc.jl

@@ -68,5 +68,9 @@
         @test isdegenerate(Steane7()) == false
         @test isdegenerate(Steane7(), 2) == true
         @test isdegenerate(Bitflip3()) == true
+        @test isdegenerate(Bicycle(6, 4)) == true
+        @test isdegenerate(Bicycle(18, 16)) == true
+        @test isdegenerate(Bicycle(200, 120)) == true
+        @test isdegenerate(Unicycle(21, [1, 3, 8, 9, 12])) == true
     end
-end
+end

test/test_ecc_base.jl

@@ -152,6 +152,8 @@ const code_instance_args = Dict(
     :Surface => [(3,3), (4,4), (3,6), (4,3), (5,5)],
     :Gottesman => [3, 4, 5],
     :CSS => (c -> (parity_checks_x(c), parity_checks_z(c))).([Shor9(), Steane7(), Toric(4, 4)]),
+    :Bicycle => [(6, 4), (10, 6)],
+    :Unicycle => [(21, [1, 3, 8, 9, 12])],
     :Concat => [(Perfect5(), Perfect5()), (Perfect5(), Steane7()), (Steane7(), Cleve8()), (Toric(2, 2), Shor9())],
     :CircuitCode => random_circuit_code_args,
     :LPCode => (c -> (c.A, c.B)).(vcat(LP04, LP118, test_gb_codes, test_bb_codes, test_mbb_codes, test_coprimeBB_codes, test_hcubic_codes, other_lifted_product_codes)),