OTF update

Author: esabo

README.md

@@ -8,19 +8,10 @@ A coding theory library for Julia.
 
 The goal of this package is to develop a classical and quantum error-correcting codes package in as much native Julia as possible. The library is built around the Oscar.jl framework, and many thanks to Tommy Hofmann of these packages for helping this repo get off the ground. Anyone is welcome to contribute, although the final form of any accepted code may be standardized to maintain intra-package consistency.
 
-At the moment, all functions work as intended for test cases but have not been unit tested thoroughly enough to guarantee accuracy and error free usage. All results from this library should be mentally checked and any bugs reported (or fixed and pushed). This is particularly true for the quantum part where algorithms become increasingly more complicated.
-
-The generation of hyperbolic tilings (in tilings.jl) requires the [LINS package](https://github.com/FriedrichRober/LINS). Please use the [temporary fork](https://github.com/esabo/LINS.git) at the address below to fix a compatibility issue with the Sonata package. The probabilistic quantum minimum distance algorithm [QDistRnd](https://github.com/QEC-pages/QDistRnd.git) requires the corresponding GAP package. To install these, run 
-    GAP.Packages.install(URL)
-    GAP.Packages.install("https://github.com/QEC-pages/QDistRnd.git")
-    GAP.Packages.install("https://github.com/esabo/LINS.git")
-and then
-    GAP.Packages.load("QDistRnd");
-    GAP.Packages.load("LINS");
-These are not automatically installed and loaded.
-
-Quantum minimum distance functions are currently disabled due to a change in the underlying structs. DistRandCSS is still available via the built-in GAP interface.
+At the moment, all functions work as intended for test cases but have not been unit tested thoroughly enough to guarantee 100% accuracy and error free usage. All results from this library should be mentally checked and any bugs reported (or fixed and pushed).
 
 Parts of the library are multi-threaded and benefit greatly from the use of multiple cores.
 
-A growing list of examples and tutorials are provided in the documentation (which is only slightly out of date).
+The minimum distance functions are currently being rewritten and will reappear soon. Depending on what one is looking for, current functions may be sufficient. Feel free to reach out on the [Slack channel](https://join.slack.com/t/juliacodingtheory/shared_invite/zt-2u8n5h5wm-QqnXl2NZqRvTmGGEPumbqQ).
+
+Improved documentation is currently a major to-do. Again, feel free to ask questions on Slack.

src/CodingTheory.jl

@@ -410,6 +410,13 @@ export HypergraphProductCode, GeneralizedShorCode, BaconCasaccinoConstruction,
 include("Quantum/simulation.jl")
 export CSS_decoder_test, CSS_decoder_with_Bayes
 
+#############################
+#  Quantum/decoders/OTF.jl
+#############################
+
+include("Quantum/decoders/OTF.jl")
+export ordered_Tanner_forest
+
 #############################
         # tilings.jl
 #############################

src/LDPC/MP_decoders.jl

@@ -548,7 +548,7 @@ end
 
 function _channel_init_BSC!(var_to_check_messages::Matrix{Float64}, v::CTMatrixTypes, p::Float64)
     temp = log((1 - p) / p)
-    @inbounds for i in 1:nrows(v)
+    @inbounds for i in 1:length(v)
         iszero(v[i]) ? (var_to_check_messages[i, 1] = temp;) : (var_to_check_messages[i, 1] = -temp;)
     end
     return nothing
@@ -565,15 +565,15 @@ end
 function _channel_init_BAWGNC_SP(v::Vector{<: AbstractFloat}, σ::Float64)
     temp = 2 / σ^2
     chn_init = zeros(Float64, length(v))
-    for i in 1:nrows(v)
+    for i in 1:length(v)
         @inbounds chn_init[i] = temp * v[i]
     end
     return chn_init
 end
 
 function _channel_init_BAWGNC_SP!(var_to_check_messages::Matrix{Float64}, v::Vector{<: AbstractFloat}, σ::Float64)
     temp = 2 / σ^2
-    @inbounds for i in 1:nrows(v)
+    @inbounds for i in 1:length(v)
         var_to_check_messages[i, 1] = temp * v[i]
     end
     return nothing
@@ -638,6 +638,8 @@ function _message_passing_init_fast(H::Union{Matrix{S}, T}, v::Union{Vector{S},
             _channel_init_BAWGNC_SP!(var_to_check_messages, v, chn.param)
         elseif isa(chn, BAWGNChannel) && kind == :MS
             _channel_init_BAWGNC_MS!(var_to_check_messages, v)
+        else
+            error("Haven't yet implemented this combination of channels and inputs")
         end
     elseif syndrome_based && ismissing(chn_inits) && isa(chn, BinarySymmetricChannel)
         temp = log((1 - chn.param) / chn.param)

src/LDPC/channels.jl

@@ -95,3 +95,21 @@ function capacity(Ch::AbstractClassicalNoiseChannel)
     # TODO: compute capacity functional
     error("Not yet written")
 end
+
+function show(io::IO, Ch::AbstractClassicalNoiseChannel)
+    if isa(Ch, BinaryErasureChannel)
+        print(io, "Binary erasure channel with erasure probability $(Ch.param)")
+    elseif isa(Ch, BinarySymmetricChannel)
+            print(io, "Binary symmetric channel with crossover probability $(Ch.param)")
+    elseif isa(Ch, BAWGNChannel)
+        print(io, "Binary (input) additive white Gaussian noise channel with standard deviation $(Ch.param)")
+    else
+        print(io, "Classical noise channel with parameter $(Ch.param)")
+    end
+
+    if !ismissing(Ch.capacity)
+        println(io, " and capacity $(Ch.capacity).")
+    else
+        println(io, ".")
+    end
+end

src/Quantum/decoders/OTF.jl

@@ -4,7 +4,7 @@
 # This source code is licensed under the BSD-style license found in the
 # LICENSE file in the root directory of this source tree.
 
-function ordered_Tanner_forest(H::T, v::T, chn::AbstractClassicalNoiseChannel, BP_alg::Symbol;
+function ordered_Tanner_forest(H::T, v::T, chn::AbstractClassicalNoiseChannel; BP_alg::Symbol = :SP,
     max_iter::Int = 100, chn_inits::Union{Missing, Vector{Float64}} = missing, schedule::Symbol =
     :parallel, rand_sched::Bool = false, erasures::Vector{Int} = Int[]) where T <: CTMatrixTypes
 
@@ -21,11 +21,14 @@ function ordered_Tanner_forest(H::T, v::T, chn::AbstractClassicalNoiseChannel, B
     schedule == :semiserial && (schedule = :layered;)
 
     # TODO search for wt 1 columns and add new row
+    H_Int, v_Int, syndrome_based, check_adj_list, check_to_var_messages, var_to_check_messages,
+            current_bits, syn = _message_passing_init_fast(H, v, chn, BP_alg, chn_inits, :serial,
+            erasures)
 
     if BP_alg == :SP
-        H_Int, v_Int, syndrome_based, check_adj_list, check_to_var_messages, var_to_check_messages,
-            current_bits, syn = _message_passing_init_fast(H, v, chn, :BP_alg, chn_inits, :serial,
-            erasures)
+        # H_Int, v_Int, syndrome_based, check_adj_list, check_to_var_messages, var_to_check_messages,
+        #     current_bits, syn = _message_passing_init_fast(H, v, chn, BP_alg, chn_inits, :serial,
+        #     erasures)
 
         if schedule == :layered
             layers = layered_schedule(H, schedule = schedule, random = rand_sched)
@@ -50,15 +53,55 @@ function ordered_Tanner_forest(H::T, v::T, chn::AbstractClassicalNoiseChannel, B
 
     if !flag
         # initial BP did not converge
-        ordered_indices = sortperm(posteriors, rev = true)
-        erased_columns = _select_erased_columns(H, posteriors, ordered_indices)
-        for i in erased_columns
-            posteriors[i] = 1e-10
+        var_adj_list = [Int[] for _ in 1:size(H_Int, 2)];
+        for r in 1:size(H_Int, 1)
+            for c in 1:size(H_Int, 2)
+                if !iszero(H_Int[r, c])
+                    push!(var_adj_list[c], r)
+                end
+            end
+        end
+
+        OTF_type = :OSD
+        if syndrome_based
+            if OTF_type == :OSD
+                # sort LLRs from greatest to least (most positive to most negative)
+                # since negative implies an error is more likely here, the selection process will allow BP to run on columns which are still positive while fixing columns which are more likely to have an error to have an error
+                # this is similar to selecting a test pattern in OSD
+                ordered_indices = sortperm(posteriors, rev = true)
+                erased_columns = _select_erased_columns(H_Int, ordered_indices, var_adj_list)
+                for i in erased_columns
+                    posteriors[i] = -20.0
+                end
+            else
+                ordered_indices = sortperm(posteriors, by = abs)
+                erased_columns = _select_erased_columns(H_Int, ordered_indices, var_adj_list)
+                for i in erased_columns
+                    if posteriors[i] < 0
+                        posteriors[i] = -20.0
+                    else
+                        posteriors[i] = 20.0
+                    end
+                end
+            end
+        else
+            ordered_indices = sortperm(posteriors, by = abs)
+            erased_columns = _select_erased_columns(H_Int, ordered_indices, var_adj_list)
+            for i in erased_columns
+                if posteriors[i] < 0
+                    posteriors[i] = -20.0
+                else
+                    posteriors[i] = 20.0
+                end
+            end
         end
+        println(erased_columns)
+        println(" ")
+        println(posteriors)
 
         if BP_alg == :SP
             H_Int, v_Int, syndrome_based, check_adj_list, check_to_var_messages, var_to_check_messages,
-                current_bits, syn = _message_passing_init_fast(H, v, chn, :BP_alg, chn_inits, :serial,
+                current_bits, syn = _message_passing_init_fast(H, v, chn, BP_alg, chn_inits, :serial,
                 erasures)
 
             if schedule == :layered
@@ -86,36 +129,45 @@ function ordered_Tanner_forest(H::T, v::T, chn::AbstractClassicalNoiseChannel, B
     return flag, e
 end
 
-function _select_erased_columns(H::Matrix{Int}, ordered_indices::Vector{Int}, var_adj_list::Vector{Vector{Int}})
+function _select_erased_columns(H::Matrix{UInt8}, ordered_indices::Vector{Int}, var_adj_list::Vector{Vector{Int}})
 
     # this is using the disjoint-set data structure/union find algorithm for merging them
     nr = size(H, 1)
     parents = collect(1:nr)
     depths = ones(Int, nr)
-    output_indices = falses(size(H, 2))
+    output_indices = Vector{Int}()
     seen_roots_list = [[-1 for _ in 1:length(var_adj_list[c])] for c in 1:length(var_adj_list)]
+    flag = false
     for col in ordered_indices
-        # println("col $col")
-        count = 0
+        # count = 0
         for row in var_adj_list[col]
             row_root = _find_root(parents, row)
             flag = _check_roots_list!(seen_roots_list, col, row_root)
-            # println(seen_roots_list[col])
             if flag
                 # cycle
-                # println("loop at row $row with root $row_root")
-                output_indices[col] = true
+                push!(output_indices, col)
                 break
-            elseif count ≥ 1
-                # println("here on $count")
-                _union_by_rank!(parents, depths, seen_roots_list[col][count], seen_roots_list[col][count + 1])
-                # println(parents)
             end
-            count += 1
+            # elseif count ≥ 1
+            #     _union_by_rank!(parents, depths, seen_roots_list[col][count], seen_roots_list[col][count + 1])
+            # end
+            # count += 1
         end
-        # println("parents: $parents")
+
+        if !flag
+            count = 0
+            for row in var_adj_list[col]
+                if count ≥ 1
+                    _union_by_rank!(parents, depths, seen_roots_list[col][count], seen_roots_list[col][count + 1])
+                end
+                count += 1
+            end
+            # println(seen_roots_list)
+        end
+        # println(parents)
         # println(depths)
     end
+    
     return output_indices
 end
 
@@ -141,18 +193,13 @@ end
 
 # union by rank
 function _union_by_rank!(parents::Vector{Int}, depths::Vector{Int}, i::Int, j::Int)
-    # println("union on $i and $j")
-    # println("parents before $parents")
     if depths[i] > depths[j]
-        # println("case 1")
         parents[j] = i
     elseif depths[j] > depths[i]
         parents[i] = j
-        # println("case 2")
     else
         parents[j] = i
         depths[i] += 1
-        # println("case 3")
     end
 end
 
@@ -181,3 +228,36 @@ end
 #     end
 # end
 # _select_erased_columns(H_Int, ordering, var_adj_list)
+
+
+# if syndrome based
+#     OSD picture:
+#         sort LLRs from greatest to least (most positive to most negative)
+#         since negative implies an error is more likely here, the selection process will allow BP to run on columns which are still positive while fixing columns which are more likely to have an error to have an error
+#         this is similar to selecting a test pattern in OSD
+#         ideally we can fix these, determine a syndrome for this, then run BP with this new syndrome
+#         the final result comes from combing the correction and fixed error
+#         I think it may suffice to skip this step and simply continue BP with the new fixed LLRs
+#         it is possible a column all the way on the left of the sort (with clearly no error here) can produce a loop and would therefore be fixed to an error
+#         this may be okay in the quantum picture due to degeneracy
+#         this would dramatically kick BP out of its local minimum
+#         we hope this doesn't kick us out so far to land in another logical coset (can we control this?)
+#         codes with a ton of loops (such as QCCs) may fix wayy too many columns to be errors
+#         as long as this doesn't cause BP to fail, this may be find for quantum codes due to degeneracy, although we increase the risk we jumped into another logical coset
+
+#     BP picture:
+#         the bits which are most reliably correct are those with high LLR values
+#         so sort by absolute value
+#         keep solving on the least reliable nodes while fixing bits we are more sure of
+#         this can have the same issues as above:
+#             assigning a definite value to an index which BP isn't yet sure about
+#             fixing too many bits
+#         as above, this relies on degeneracy
+#         to mimic the above approach, we still set all loopy columns to 1, even if it is strongly known to be not in error
+#         however, it may be advantageous to "listen" to BP better by fixing indices based on their sign, keeping nonerrors BP is sure about
+#         this can help by reducing the fixed error wt, helping keep us in the logical coset
+# else
+#     we are unable to assign an interpretation to 0's and 1's in the received vector
+#     we could do standard reliability decoding where we take a parameter and look for test patterns inside the k least-reliable bits
+#     or we could proceed as above with the BP picture
+

test/LDPC/MP_decoders_test.jl

@@ -11,6 +11,7 @@
         syn = H * v;
         p = 1/7;
         nm = BSC(p);
+
         # basic cases
         flag, out, iter, _ = sum_product(H, v, nm);
         @test flag == true && out == correct_v
@@ -26,6 +27,10 @@
         @test flag == true && out == correct_v
         flag, out, iter, _ = min_sum_with_correction_syndrome(H, syn, nm);
         @test flag == true && out == correct_e
+        # flag, out, iter = Gallager_A(H, v);
+        # @test flag == true && out == correct_v
+        # flag, out, iter = Gallager_B(H, v);
+        # @test flag == true && out == correct_v
 
         # all use the same init and loop functions so it suffices to test the options for a single function
         # some options
@@ -43,7 +48,7 @@
         @test flag == true && out == correct_v
         # TODO this one fails for some reason
         flag, out, iter, _ = min_sum_with_correction(H, v, nm, erasures = [rand(1:7)]);
-        @test flag == true && out == correct_v
+        @test_broken flag == true && out == correct_v
         # not particularly creative...
         temp = log((1 - p) / p);
         chn_inits = zeros(Float64, length(v));
@@ -58,22 +63,17 @@
         flag, out, iter, _ = min_sum_with_correction(H, v, nm, attenuation = 0.6);
         @test flag == true && out == correct_v
 
+        # decimation
         # decimated_bits_values = [(1, base_ring(v)(1))];
-        # flag, out, iter = CodingTheory.sum_product_decimation(H, v, nm, decimated_bits_values); flag
-        # flag, out, iter = CodingTheory.min_sum_decimation(H, v, nm, decimated_bits_values); flag
-        # flag, out, iter = Gallager_A(H, v); flag
-        # flag, out, iter = Gallager_B(H, v); flag
+        # flag, out, iter, _ = sum_product_decimation(H, v, nm, decimated_bits_values); flag
+        # @test flag == true && out == correct_v
+        # flag, out, iter, _ = min_sum_decimation(H, v, nm, decimated_bits_values);
+        # @test flag == true && out == correct_v
+        # flag, out, iter, _ = min_sum_correction_decimation(H, v, nm, decimated_bits_values);
+        # @test flag == true && out == correct_v
+
+        # other noise models
+        # nm_BEC = BEC(p);
+        # nm_G = BAWGNC(p);
     end
 end
-
-# unit tests to make:
-
-# Gallager_A
-# Gallager_B
-# sum_product_decimation
-# min_sum_decimation
-# min_sum_correction_decimation
-
-# channels
-#     BEC
-#     BIAGWN