proper testing of the logical error rate and move to column-major ord…

…er (#13)

Project.toml

@@ -1,7 +1,7 @@
 name = "LDPCDecoders"
 uuid = "3c486d74-64b9-4c60-8b1a-13a564e77efb"
 authors = ["Krishna Praneet Gudipaty", "Stefan Krastanov", "QuantumSavory contributors"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/LDPCDecoders.jl

@@ -10,7 +10,7 @@ using RowEchelon
 
 
 export
-    decode!, batchdecode!, 
+    decode!, batchdecode!,
     BeliefPropagationDecoder
 
 include("generator.jl")
@@ -20,8 +20,13 @@ include("bp_decoder.jl")
 include("bp_simulator.jl")
 include("it_decoder.jl")
 include("parity_generator.jl")
-include("syndrome_decoder.jl")
+
+include("decoders/abstract_decoder.jl")
+include("decoders/belief_propogation.jl")
+include("syndrome_bp_decoder.jl")
 include("syndrome_simulator.jl")
 include("syndrome_it_decoder.jl")
 include("syndrome_it_simulate.jl")
+
+
 end

src/decoders/belief_propogation.jl

@@ -1,15 +1,12 @@
 using SparseArrays
 
-include("abstract_decoder.jl")
-
-# TODO: Better nomenclature
 struct BeliefPropagationScratchSpace
   "Log probabilities"
   log_probabs::Vector{Float64}
-  
+
   "Channel probabilities"
   channel_probs::Vector{Float64}
-  
+
   "Bit to check message matrix"
   bit_2_check::Matrix{Float64}
 
@@ -18,33 +15,41 @@ struct BeliefPropagationScratchSpace
 
   "predicted error for each bit"
   err::Vector{Float64}
-end 
+end
+
+function BeliefPropagationScratchSpace(n, s, per)
+  return BeliefPropagationScratchSpace(zeros(n), fill(per, n), zeros(s, n), zeros(s, n), zeros(n))
+end
 
 struct BeliefPropagationDecoder <: AbstractDecoder
   "Physical error rate"
   per::Float64
 
   "Number of max iterations of Belief propagation decoder"
   max_iters::Int
-  
+
   "Num of parity checks i.e. number of rows of parity check matrix"
   s::Int
-  
+
   "Num of bits in the code i.e number of columns of parity check matrix"
   n::Int
-  
+
   "Sparse form of the parity check matrix"
   sparse_H::SparseArrays.SparseMatrixCSC{Bool,Int}
-  
+
   "Sparse form of the transform of the parity check matrix"
   sparse_HT::SparseArrays.SparseMatrixCSC{Bool,Int}
+
+  "Scratch space for the decoder"
+  scratch::BeliefPropagationScratchSpace
 end
 
 function BeliefPropagationDecoder(H, per::Float64, max_iters::Int)
   s, n = size(H)
   sparse_H = sparse(H)
   sparse_HT = sparse(H')
-  return BeliefPropagationDecoder(per, max_iters, s, n, sparse_H, sparse_HT)
+  scratch = BeliefPropagationScratchSpace(n, s, per)
+  return BeliefPropagationDecoder(per, max_iters, s, n, sparse_H, sparse_HT, scratch)
 end
 
 """
@@ -58,21 +63,23 @@ Function to reset the scratch space for Belief propagation decoding
 ```jldoctest
 julia> decoder = BeliefPropagationDecoder(LDPCDecoders.parity_check_matrix(1000, 10, 9), 0.01, 100);
 
-julia> scratch = BeliefPropagationScratchSpace(zeros(decoder.n), fill(decoder.per, decoder.n), 
+julia> scratch = BeliefPropagationScratchSpace(zeros(decoder.n), fill(decoder.per, decoder.n),
                 zeros(decoder.s, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.n));
 
-julia> reset!(scratch, decoder); 
+julia> reset!(scratch, decoder);
 ````
 """
-function reset!(bp_setup::BeliefPropagationScratchSpace, bp_decoder::BeliefPropagationDecoder)
+function reset!(bp_decoder::BeliefPropagationDecoder)
+  bp_setup = bp_decoder.scratch
   bp_setup.log_probabs .= 0.0
   bp_setup.channel_probs .= bp_decoder.per
   bp_setup.bit_2_check .= 0.0
   bp_setup.check_2_bit .= 0.0
   bp_setup.err .= 0.0
+  bp_decoder
 end
 
-# TODO: Belief decoder type, syndrome, error -> should be modified 
+# TODO: Belief decoder type, syndrome, error -> should be modified
 
 """
 Function to decode given the parity check matrix, syndrome and error
@@ -84,31 +91,26 @@ Function to decode given the parity check matrix, syndrome and error
 
 # Examples
 ```jldoctest
-julia> decoder = BeliefPropagationDecoder(LDPCDecoders.parity_check_matrix(1000, 10, 9), 0.01, 100);
+julia> H = LDPCDecoders.parity_check_matrix(1000, 10, 9);
+
+julia> decoder = BeliefPropagationDecoder(H, 0.01, 100);
 
 julia> error = rand(1000) .< 0.01;
 
-julia> syndrome::BitArray{1} = (H * error) .% 2;
+julia> syndrome = (H * error) .% 2;
 
-julia> decoded_error::BitArray{1} = zeros(1000);
+julia> guess, success = decode!(decoder, syndrome, zeros(1000));
 
-julia> decode!(decoder, syndrome, decoded_error)
+julia> error == guess
 true
 ```
 """
-
-function decode!(decoder::BeliefPropagationDecoder, syndrome::BitArray{1}, error::BitArray{1})
-
-  setup = BeliefPropagationScratchSpace(zeros(decoder.n), fill(decoder.per, decoder.n), 
-            zeros(decoder.s, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.n))
-
-  reset!(setup, decoder)
-  success = syndrome_decode!(decoder, setup, syndrome)
-  copyto!(error, setup.err)
-  return success
+function decode!(decoder::BeliefPropagationDecoder, syndrome::AbstractVector, error::AbstractVector) # TODO check if casting to bitarrays helps with performance -- if it does, set up warnings to the user for cases where they have not done the casting
+  reset!(decoder)
+  success = syndrome_decode!(decoder, decoder.scratch, syndrome) # TODO: Delete syndrome_decode! and just move its content here. There is no point in this indirection.
+  return decoder.scratch.err, success
 end
 
-
 """
 Function to batch decode given the parity check matrix, syndrome and error
 
@@ -123,32 +125,28 @@ Scratch space allocations are done once and re-used for better performance
 ```jldoctest
 julia> decoder = BeliefPropagationDecoder(LDPCDecoders.parity_check_matrix(1000, 10, 9), 0.01, 100);
 
-julia> errors::BitArray{2} = rand(Base.Float64, (10,1000)) .< 0.01;
+julia> samples = 100
 
-julia> syndromes::BitArray{2} = zeros(10, 900);
+julia> errors = rand(1000,samples) .< 0.01;
 
-julia> for i in axes(errors, 1)
-             syn = (H * errors[i, :]) .% 2;
-             syndromes[i, :] = syn;
-           end  
+julia> syndromes = zeros(900, samples);
 
-julia> batchdecode!(decoder, syndromes, errors) < 0.9
-true
-"""
-function batchdecode!(decoder::BeliefPropagationDecoder, syndromes::BitArray{2}, errors::BitArray{2})
-
-  @assert size(syndromes, 1) == size(errors, 1)
-  num_trials::Int = size(syndromes, 1)
-
-  setup = BeliefPropagationScratchSpace(zeros(decoder.n), fill(decoder.per, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.n))
+julia> syndromes = H*errors .% 2
 
-  success::BitArray{1} = zeros(num_trials)
+julia> guesses, successes = batchdecode!(decoder, syndromes, zero(errors));
 
-  for i in axes(syndromes, 1)
-    reset!(setup, decoder)
-    success[i] = syndrome_decode!(decoder, setup, syndromes[i, :])
-    errors[i, :] = setup.err
+julia> sum((guesses[:,i] == errors[:,i] for i in 1:samples)) > 0.995*samples
+"""
+function batchdecode!(decoder::BeliefPropagationDecoder, syndromes, errors)
+  @assert size(syndromes, 2) == size(errors, 2)
+  num_trials::Int = size(syndromes, 2)
+  success::AbstractVector{Bool} = zeros(num_trials)
+
+  for i in axes(syndromes, 2)
+    reset!(decoder)
+    success[i] = syndrome_decode!(decoder, decoder.scratch, syndromes[:, i])
+    errors[:, i] = decoder.scratch.err
   end
 
-  return success
-end
+  return errors, success
+end

src/syndrome_decoder.jl → src/syndrome_bp_decoder.jl

@@ -1,6 +1,4 @@
-include("decoders/belief_propogation.jl")
-
-# TODO: Depricated, to be removed. Only new implementation is required. 
+# TODO: Depricated, to be removed. Only new implementation is required.
 function syndrome_decode(pcm, pcmT, syndrome, max_iters, channel_probs, b2c, c2b, log_probabs, error)
 
   # Get size of Parity check matrix
@@ -103,8 +101,8 @@ end
 
 # In-place variant
 # TODO: Currently very slow. Need to improve this
-function syndrome_decode!(decoder::BeliefPropagationDecoder, setup::BeliefPropagationScratchSpace, syndrome::BitArray{1})
-  
+function syndrome_decode!(decoder::BeliefPropagationDecoder, setup::BeliefPropagationScratchSpace, syndrome::AbstractVector) # TODO check if casting to bitarrays helps with performance, and raise warnings for the user if there is performance left on the table
+
   # Get size of Parity check matrix
   m, n = size(decoder.sparse_H)
   rows = rowvals(decoder.sparse_H)
@@ -141,7 +139,7 @@ function syndrome_decode!(decoder::BeliefPropagationDecoder, setup::BeliefPropag
       end
     end
 
-    
+
     # Bit to check messages
     for j in 1:n
       temp = setup.channel_probs[j] / (1 - setup.channel_probs[j])
@@ -158,7 +156,7 @@ function syndrome_decode!(decoder::BeliefPropagationDecoder, setup::BeliefPropag
       setup.log_probabs[j] = log(1 / temp)
       if temp >= 1
         setup.err[j] = 1
-      else 
+      else
         setup.err[j] = 0
       end
 
@@ -181,7 +179,7 @@ function syndrome_decode!(decoder::BeliefPropagationDecoder, setup::BeliefPropag
       return true
     end
   end
-  
+
   return false
 
 end

test/test_bp_decoder.jl

@@ -5,60 +5,52 @@ using LDPCDecoders
 
   """Test for belief propagation decoder"""
   function test_bp_decoder()
-    H::BitArray{2} = LDPCDecoders.parity_check_matrix(1000, 10, 9)
+    H = LDPCDecoders.parity_check_matrix(1000, 10, 9)
     per = 0.01
-    err::BitArray{1} = rand(1000) .< per
-    syn::BitArray{1} = (H * err) .% 2
+    err = rand(1000) .< per
+    syn = (H * err) .% 2
 
     bpd = BeliefPropagationDecoder(H, per, 100)
-    decoded_error::BitArray{1} = zeros(1000)
+    decoded_error = zeros(1000)
 
-    success = decode!(bpd, syn, decoded_error)
+    guess, success = decode!(bpd, syn, decoded_error)
 
-    return success
+    return guess == err
   end
 
   """Test for batch belief propagation decoder"""
   function test_bp_decoder_batch()
-    H::BitArray{2} = LDPCDecoders.parity_check_matrix(1000, 10, 9)
-    per = 0.1
+    H = LDPCDecoders.parity_check_matrix(1000, 10, 9)
+    per = 0.01
     num_trials = 100
-    errors::BitArray{2} = rand(Base.Float64, (num_trials,1000)) .< per
-    syndromes::BitArray{2} = zeros(num_trials, 900)
-
-    for i in axes(errors, 1)
-      syn = (H * errors[i, :]) .% 2
-      syndromes[i, :] = syn
-    end
-
+    errors = rand(1000,num_trials) .< per
+    syndromes = H*errors .% 2
     bpd = BeliefPropagationDecoder(H, per, 100)
-    success::BitArray{1} = batchdecode!(bpd, syndromes, errors)
-    ler = (num_trials - sum(success)) / num_trials
+    guesses, successes = batchdecode!(bpd, syndromes, zero(errors))
+    actual_successes = [guesses[:,i]==errors[:,i] for i in 1:num_trials]
+    ler = (num_trials - sum(actual_successes)) / num_trials
     return ler
   end
 
   """Test for batch belief propagation decoder"""
-  function test_depricated_syndrome_decoder()
-    H::BitArray{2} = LDPCDecoders.parity_check_matrix(1000, 10, 9)
-    per = 0.1
+  function test_deprecated_syndrome_decoder()
+    H = LDPCDecoders.parity_check_matrix(1000, 10, 9)
+    per = 0.01
     num_trials = 100
-    errors::BitArray{2} = rand(Base.Float64, (num_trials,1000)) .< per
-    syndromes::BitArray{2} = zeros(num_trials, 900)
-
-    for i in axes(errors, 1)
-      syn = (H * errors[i, :]) .% 2
-      syndromes[i, :] = syn
-    end
+    errors = rand(Base.Float64, (1000,num_trials)) .< per
+    syndromes = H*errors .% 2
 
     decoder = BeliefPropagationDecoder(H, per, 100)
-    success::BitArray{1} = batchdecode!(decoder, syndromes, errors)
+    guess_errors, success = batchdecode!(decoder, syndromes, zero(errors))
+    actual_successes = 0
 
-    for i in axes(syndromes, 1)
-      args = decoder.sparse_H, decoder.sparse_HT, syndromes[i, :], 100, fill(decoder.per, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.n), zeros(decoder.n)
+    for i in axes(syndromes, 2)
+      args = decoder.sparse_H, decoder.sparse_HT, syndromes[:, i], 100, fill(decoder.per, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.s, decoder.n), zeros(decoder.n), zeros(decoder.n)
       decoded_error, converged = LDPCDecoders.syndrome_decode(args...)
+      actual_successes += decoded_error == errors[:, i]
 
       # Test whether the new decoder works similar to the old one
-      @assert decoded_error == errors[i, :]
+      @assert decoded_error == guess_errors[:, i]
     end
 
     ler = (num_trials - sum(success)) / num_trials
@@ -68,6 +60,6 @@ using LDPCDecoders
   @test test_bp_decoder()
 
   # There is a low possibility of these tests failing
-  @test test_bp_decoder_batch() < 0.9
-  @test test_depricated_syndrome_decoder() < 0.9
+  @test test_bp_decoder_batch() < 0.005
+  @test test_deprecated_syndrome_decoder() < 0.005
 end