From 983f9a1c9582caa67babe9d893603784c05ff618 Mon Sep 17 00:00:00 2001 From: wnhack Date: Sun, 11 Jul 2021 23:07:39 -0400 Subject: [PATCH] added hamming code and gallager ldpc eccs and commandline tools for benchmarking --- benchmarking/benchmarking.go | 364 ++++++++++++++++++ benchmarking/benchmarking_test.go | 170 ++++++++ benchmarking/randomdata.go | 83 ++++ cmd/create.go | 72 ++++ cmd/internal/create/gallager/gallager.go | 66 ++++ cmd/internal/create/hamming/hamming.go | 63 +++ cmd/internal/tools/bec/bec.go | 74 ++++ cmd/internal/tools/bec/simple/simple.go | 145 +++++++ cmd/internal/tools/bsc/bsc.go | 73 ++++ cmd/internal/tools/bsc/dwbf/dwbf.go | 152 ++++++++ cmd/internal/tools/bsc/gallager/gallager.go | 147 +++++++ cmd/internal/tools/csv/csv.go | 90 +++++ cmd/internal/tools/tools.go | 114 ++++++ cmd/root.go | 18 + cmd/tools.go | 128 ++++++ data/linearblock/hamming/hamming_15_11.json | 1 + .../hamming/hamming_15_11_bec_results.json | 1 + .../hamming_15_11_bsc_dwbf_results.json | 1 + .../hamming_15_11_bsc_gallager_results.json | 1 + data/linearblock/hamming/hamming_7_4.json | 1 + .../hamming/hamming_7_4_bec_results.json | 1 + .../hamming/hamming_7_4_bsc_dwbf_results.json | 1 + .../hamming_7_4_bsc_gallager_results.json | 1 + go.mod | 23 ++ linearblock/cycle.go | 243 ++++++++++++ linearblock/cycle_test.go | 67 ++++ linearblock/girth.go | 160 ++++++++ linearblock/girth_test.go | 85 ++++ linearblock/hamming/hamming.go | 46 +++ linearblock/hamming/hamming_test.go | 27 ++ linearblock/internal/gaussianjordan.go | 226 +++++++++++ linearblock/internal/gaussianjordan_test.go | 40 ++ linearblock/internal/utils.go | 107 +++++ linearblock/internal/utils_test.go | 40 ++ linearblock/ldpc/gallager/gallager.go | 132 +++++++ linearblock/linearblock.go | 181 +++++++++ linearblock/linearblock_test.go | 31 ++ linearblock/messagepassing/bec/bec.go | 25 ++ .../messagepassing/bec/iterative/simple.go | 123 ++++++ .../bec/iterative/simple_test.go | 33 ++ .../bitflipping/harddecision/bitflipping.go | 22 ++ .../bitflipping/harddecision/dwbf.go | 156 ++++++++ .../bitflipping/harddecision/dwdf_test.go | 64 +++ .../bitflipping/harddecision/gallager.go | 91 +++++ .../bitflipping/harddecision/gallager_test.go | 61 +++ main.go | 7 + 46 files changed, 3757 insertions(+) create mode 100644 benchmarking/benchmarking.go create mode 100644 benchmarking/benchmarking_test.go create mode 100644 benchmarking/randomdata.go create mode 100644 cmd/create.go create mode 100644 cmd/internal/create/gallager/gallager.go create mode 100644 cmd/internal/create/hamming/hamming.go create mode 100644 cmd/internal/tools/bec/bec.go create mode 100644 cmd/internal/tools/bec/simple/simple.go create mode 100644 cmd/internal/tools/bsc/bsc.go create mode 100644 cmd/internal/tools/bsc/dwbf/dwbf.go create mode 100644 cmd/internal/tools/bsc/gallager/gallager.go create mode 100644 cmd/internal/tools/csv/csv.go create mode 100644 cmd/internal/tools/tools.go create mode 100644 cmd/root.go create mode 100644 cmd/tools.go create mode 100644 data/linearblock/hamming/hamming_15_11.json create mode 100644 data/linearblock/hamming/hamming_15_11_bec_results.json create mode 100644 data/linearblock/hamming/hamming_15_11_bsc_dwbf_results.json create mode 100644 data/linearblock/hamming/hamming_15_11_bsc_gallager_results.json create mode 100644 data/linearblock/hamming/hamming_7_4.json create mode 100644 data/linearblock/hamming/hamming_7_4_bec_results.json create mode 100644 data/linearblock/hamming/hamming_7_4_bsc_dwbf_results.json create mode 100644 data/linearblock/hamming/hamming_7_4_bsc_gallager_results.json create mode 100644 go.mod create mode 100644 linearblock/cycle.go create mode 100644 linearblock/cycle_test.go create mode 100644 linearblock/girth.go create mode 100644 linearblock/girth_test.go create mode 100644 linearblock/hamming/hamming.go create mode 100644 linearblock/hamming/hamming_test.go create mode 100644 linearblock/internal/gaussianjordan.go create mode 100644 linearblock/internal/gaussianjordan_test.go create mode 100644 linearblock/internal/utils.go create mode 100644 linearblock/internal/utils_test.go create mode 100644 linearblock/ldpc/gallager/gallager.go create mode 100644 linearblock/linearblock.go create mode 100644 linearblock/linearblock_test.go create mode 100644 linearblock/messagepassing/bec/bec.go create mode 100644 linearblock/messagepassing/bec/iterative/simple.go create mode 100644 linearblock/messagepassing/bec/iterative/simple_test.go create mode 100644 linearblock/messagepassing/bitflipping/harddecision/bitflipping.go create mode 100644 linearblock/messagepassing/bitflipping/harddecision/dwbf.go create mode 100644 linearblock/messagepassing/bitflipping/harddecision/dwdf_test.go create mode 100644 linearblock/messagepassing/bitflipping/harddecision/gallager.go create mode 100644 linearblock/messagepassing/bitflipping/harddecision/gallager_test.go create mode 100644 main.go diff --git a/benchmarking/benchmarking.go b/benchmarking/benchmarking.go new file mode 100644 index 0000000..daacf2f --- /dev/null +++ b/benchmarking/benchmarking.go @@ -0,0 +1,364 @@ +package benchmarking + +import ( + "context" + "fmt" + "github.com/cheggaaa/pb/v3" + "github.com/nathanhack/avgstd" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + mat "github.com/nathanhack/sparsemat" + "github.com/nathanhack/threadpool" + mat2 "gonum.org/v1/gonum/mat" + "math" + "sync" +) + +type Stats struct { + ChannelCodewordError avgstd.AvgStd // probability of a bit error after channel errors are fixed + ChannelMessageError avgstd.AvgStd // probability of a bit error after channel errors are fixed + ChannelParityError avgstd.AvgStd // probability of a bit error after channel errors are fixed +} + +func (s Stats) String() string { + return fmt.Sprintf("{Codeword:%0.02f(+/-%0.02f), Message:%0.02f(+/-%0.02f), Parity:%0.02f(+/-%0.02f)}", + s.ChannelCodewordError.Mean, math.Sqrt(s.ChannelCodewordError.SampledVariance()), + s.ChannelMessageError.Mean, math.Sqrt(s.ChannelMessageError.SampledVariance()), + s.ChannelParityError.Mean, math.Sqrt(s.ChannelParityError.SampledVariance()), + ) +} + +type Checkpoints func(updatedStats Stats) + +type BinaryMessageConstructor func(trial int) (message mat.SparseVector) + +//specfic to BSC +type BinarySymmetricChannelEncoder func(message mat.SparseVector) (codeword mat.SparseVector) +type BinarySymmetricChannel func(codeword mat.SparseVector) (channelInducedCodeword mat.SparseVector) +type BinarySymmetricChannelCorrection func(originalCodeword, channelInducedCodeword mat.SparseVector) (fixedChannelInducedCodeword mat.SparseVector) +type BinarySymmetricChannelMetrics func(originalMessage, originalCodeword, fixedChannelInducedCodeword mat.SparseVector) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) + +//specific to BEC +type BinaryErasureChannelEncoder func(message mat.SparseVector) (codeword []bec.ErasureBit) +type BinaryErasureChannel func(codeword []bec.ErasureBit) (channelInducedCodeword []bec.ErasureBit) +type BinaryErasureChannelCorrection func(originalCodeword, channelInducedCodeword []bec.ErasureBit) (fixedChannelInducedCodeword []bec.ErasureBit) +type BinaryErasureChannelMetrics func(originalMessage mat.SparseVector, originalCodeword, fixedChannelInducedCodeword []bec.ErasureBit) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) + +//specific to BPSK +type BPSKChannelEncoder func(message mat.SparseVector) (codeword mat2.Vector) +type BPSKChannel func(codeword mat2.Vector) (channelInducedCodeword mat2.Vector) +type BPSKChannelCorrection func(originalCodeword, channelInducedCodeword mat2.Vector) (fixedChannelInducedCodeword mat2.Vector) +type BPSKChannelMetrics func(originalMessage mat.SparseVector, originalCodeword, fixedChannelInducedCodeword mat2.Vector) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) + +func BenchmarkBSC(ctx context.Context, + trials int, threads int, + createMessage BinaryMessageConstructor, + encode BinarySymmetricChannelEncoder, + channel BinarySymmetricChannel, + codewordRepair BinarySymmetricChannelCorrection, + metrics BinarySymmetricChannelMetrics, + checkpoints Checkpoints, + showProgress bool) Stats { + return BenchmarkBSCContinueStats(ctx, trials, threads, createMessage, encode, channel, codewordRepair, metrics, checkpoints, Stats{}, showProgress) +} + +func BenchmarkBSCContinueStats(ctx context.Context, + trials int, threads int, + createMessage BinaryMessageConstructor, + encode BinarySymmetricChannelEncoder, + channel BinarySymmetricChannel, + codewordRepair BinarySymmetricChannelCorrection, + metrics BinarySymmetricChannelMetrics, + checkpoints Checkpoints, + previousStats Stats, + showProgress bool) Stats { + trialsToRun := trials - previousStats.ChannelCodewordError.Count + if trialsToRun <= 0 { + return previousStats + } + + var bar *pb.ProgressBar + if showProgress { + bar = pb.StartNew(trialsToRun) + } + + pool := threadpool.New(ctx, threads, trialsToRun) + statsMux := sync.Mutex{} + + trial := func(i int) { + if showProgress { + bar.Increment() + } + //we create a random message + message := createMessage(i) + + // encode to get our codeword + codeword := encode(message) + + // send through the channel to get channel induced errors + channelInducedCodeword := channel(codeword) + + // repair the codeword (if possible) + repaired := codewordRepair(codeword, channelInducedCodeword) + + // get metrics + percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors := metrics(message, codeword, repaired) + + statsMux.Lock() + previousStats.ChannelCodewordError.Update(percentFixedCodewordErrors) + previousStats.ChannelMessageError.Update(percentFixedMessageErrors) + previousStats.ChannelParityError.Update(percentFixedParityErrors) + if checkpoints != nil { + checkpoints(previousStats) //give them the updated checkpoint + } + statsMux.Unlock() + } + + for i := previousStats.ChannelCodewordError.Count; i < trials; i++ { + tmp := i + pool.Add(func() { trial(tmp) }) + } + pool.Wait() + if showProgress { + bar.Finish() + } + return previousStats +} + +func BenchmarkBEC(ctx context.Context, + trials, threads int, + createMessage BinaryMessageConstructor, + encode BinaryErasureChannelEncoder, + channel BinaryErasureChannel, + codewordRepair BinaryErasureChannelCorrection, + metrics BinaryErasureChannelMetrics, + checkpoints Checkpoints, showBar bool) Stats { + return BenchmarkBECContinueStats(ctx, trials, threads, createMessage, encode, channel, codewordRepair, metrics, checkpoints, Stats{}, showBar) +} + +func BenchmarkBECContinueStats( + ctx context.Context, + trials, threads int, + createMessage BinaryMessageConstructor, + encode BinaryErasureChannelEncoder, + channel BinaryErasureChannel, + codewordRepair BinaryErasureChannelCorrection, + metrics BinaryErasureChannelMetrics, + checkpoints Checkpoints, + previousStats Stats, + showProgressBar bool) Stats { + trialsToRun := trials - previousStats.ChannelCodewordError.Count + if trialsToRun <= 0 { + return previousStats + } + + var bar *pb.ProgressBar + if showProgressBar { + bar = pb.StartNew(trialsToRun) + } + + pool := threadpool.New(ctx, threads, trialsToRun) + statsMux := sync.Mutex{} + + trial := func(i int) { + if showProgressBar { + bar.Increment() + } + //we create a random message + message := createMessage(i) + + // encode to get our codeword + codeword := encode(message) + + // send through the channel to get channel induced errors + channelInducedCodeword := channel(codeword) + + // repair the codeword (if possible) and return metrics + repaired := codewordRepair(codeword, channelInducedCodeword) + + // get metrics + percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors := metrics(message, codeword, repaired) + + statsMux.Lock() + previousStats.ChannelCodewordError.Update(percentFixedCodewordErrors) + previousStats.ChannelMessageError.Update(percentFixedMessageErrors) + previousStats.ChannelParityError.Update(percentFixedParityErrors) + + if checkpoints != nil { + checkpoints(previousStats) //give them the updated checkpoint + } + statsMux.Unlock() + } + + for i := previousStats.ChannelCodewordError.Count; i < trials; i++ { + pool.Add(func() { trial(i) }) + } + + pool.Wait() + if showProgressBar { + bar.Finish() + } + return previousStats +} + +func BenchmarkBPSK(ctx context.Context, + trials int, threads int, + createMessage BinaryMessageConstructor, + encode BPSKChannelEncoder, + channel BPSKChannel, + codewordRepair BPSKChannelCorrection, + metrics BPSKChannelMetrics, + checkpoints Checkpoints, showProgress bool) Stats { + return BenchmarkBPSKContinueStats(ctx, trials, threads, createMessage, encode, channel, codewordRepair, metrics, checkpoints, Stats{}, showProgress) +} + +func BenchmarkBPSKContinueStats(ctx context.Context, + trials int, threads int, + createMessage BinaryMessageConstructor, + encode BPSKChannelEncoder, + channel BPSKChannel, + codewordRepair BPSKChannelCorrection, + metrics BPSKChannelMetrics, + checkpoints Checkpoints, + previousStats Stats, + showProgress bool) Stats { + trialsToRun := trials - previousStats.ChannelCodewordError.Count + if trialsToRun <= 0 { + return previousStats + } + + var bar *pb.ProgressBar + if showProgress { + bar = pb.StartNew(trialsToRun) + } + pool := threadpool.New(ctx, threads, trialsToRun) + statsMux := sync.Mutex{} + + trial := func(i int) { + if showProgress { + bar.Increment() + } + //we create a random message + message := createMessage(i) + + // encode to get our codeword + codeword := encode(message) + + // send through the channel to get channel induced errors + channelInducedCodeword := channel(codeword) + + // repair the codeword (if possible) + repaired := codewordRepair(codeword, channelInducedCodeword) + + // get metrics + percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors := metrics(message, codeword, repaired) + + statsMux.Lock() + previousStats.ChannelCodewordError.Update(percentFixedCodewordErrors) + previousStats.ChannelMessageError.Update(percentFixedMessageErrors) + previousStats.ChannelParityError.Update(percentFixedParityErrors) + + if checkpoints != nil { + checkpoints(previousStats) //give them the updated checkpoint + } + statsMux.Unlock() + } + + for i := previousStats.ChannelCodewordError.Count; i < trials; i++ { + pool.Add(func() { trial(i) }) + } + pool.Wait() + if showProgress { + bar.Finish() + } + return previousStats +} + +//HammingDistanceErasuresToBits calculates number of bits different. +// If a and b are different sizes it assumes they are +// both aligned with the zero index (the difference is at the end) +func HammingDistanceErasuresToBits(a []bec.ErasureBit, b mat.SparseVector) int { + min := len(a) + max := b.Len() + if min > max { + min = b.Len() + max = len(a) + } + + count := 0 + for i := 0; i < min; i++ { + if int(a[i]) != b.At(i) { + count++ + } + } + return max - min + count +} + +//BitsToBPSK converts a [0,1] matrix to a [-1,1] matrix +func BitsToBPSK(a mat.SparseVector) mat2.Vector { + output := mat2.NewVecDense(a.Len(), nil) + + for i := 0; i < a.Len(); i++ { + if a.At(i) > 0 { + output.SetVec(i, 1) + } else { + output.SetVec(i, -1) + } + } + + return output +} + +//BPSKToBits conversts a BPSK vector [-1,1] to sparse vector [0,1]. +// Values >= boundary will be considered a 1, otherwise a 0. +func BPSKToBits(a mat2.Vector, boundary float64) mat.SparseVector { + result := mat.CSRVec(a.Len()) + + for i := 0; i < a.Len(); i++ { + if a.AtVec(i) >= boundary { + result.Set(i, 1) + } + } + return result +} + +//HammingDistanceBPSK calculates number of bits different. +// Assumes >=0 is 1 and <0 is 0 +// If a and b are different sizes it assumes they are +// both aligned with the zero index (the difference is at the end) +func HammingDistanceBPSK(a, b mat2.Vector) int { + min := a.Len() + max := b.Len() + if min > max { + min = b.Len() + max = a.Len() + } + + count := 0 + for i := 0; i < min; i++ { + aOne := a.AtVec(i) >= 0 + bOne := b.AtVec(i) >= 0 + if aOne != bOne { + count++ + } + } + return max - min + count +} + +//BitsToErased creates a slice of ErasureBits from the codeword passed in +func BitsToErased(codeword mat.SparseVector) []bec.ErasureBit { + output := make([]bec.ErasureBit, codeword.Len()) + for i := 0; i < codeword.Len(); i++ { + output[i] = bec.ErasureBit(codeword.At(i)) + } + return output +} + +//ErasedCount returns the number of erased bits +func ErasedCount(base []bec.ErasureBit) (count int) { + for _, e := range base { + if e == bec.Erased { + count++ + } + } + return +} diff --git a/benchmarking/benchmarking_test.go b/benchmarking/benchmarking_test.go new file mode 100644 index 0000000..90cdab8 --- /dev/null +++ b/benchmarking/benchmarking_test.go @@ -0,0 +1,170 @@ +package benchmarking + +import ( + "context" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/linearblock/hamming" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec/iterative" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision" + mat "github.com/nathanhack/sparsemat" + mat2 "gonum.org/v1/gonum/mat" + "runtime" +) + +func ExampleBenchmarkBSC() { + linearBlock, _ := hamming.New(context.Background(), 3, 0) + + createMessage := func(trial int) mat.SparseVector { + t := trial % 64 + message := mat.CSRVec(4) + for i := 0; i < 4; i++ { + message.Set(i, (t&(1<>i) + } + return message + } + + encode := func(message mat.SparseVector) (codeword mat.SparseVector) { + return linearBlock.Encode(message) + } + + channel := func(originalCodeword mat.SparseVector) (erroredCodeword mat.SparseVector) { + //since hamming can fix only one bit wrong we'll just flip one bit per codeword + return RandomFlipBitCount(originalCodeword, 1) + } + repair := func(originalCodeword, channelInducedCodeword mat.SparseVector) (fixed mat.SparseVector) { + alg := &harddecision.Gallager{ + H: linearBlock.H, + } + + return harddecision.BitFlipping(alg, linearBlock.H, channelInducedCodeword, 50) + } + + metrics := func(originalMessage, originalCodeword, fixedChannelInducedCodeword mat.SparseVector) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) { + codewordErrors := originalCodeword.HammingDistance(fixedChannelInducedCodeword) + message := linearBlock.Decode(fixedChannelInducedCodeword) + messageErrors := message.HammingDistance(originalMessage) + parityErrors := codewordErrors - messageErrors + + percentFixedCodewordErrors = float64(codewordErrors) / float64(linearBlock.CodewordLength()) + percentFixedMessageErrors = float64(messageErrors) / float64(linearBlock.MessageLength()) + percentFixedParityErrors = float64(parityErrors) / float64(linearBlock.ParitySymbols()) + return + } + + checkpoint := func(updatedStats Stats) {} + + stats := BenchmarkBSC(context.Background(), 1, 1, createMessage, encode, channel, repair, metrics, checkpoint, false) + + fmt.Println("Bit Error Probability :", stats) + //Output: + // Bit Error Probability : {Codeword:0.00(+/-0.00), Message:0.00(+/-0.00), Parity:0.00(+/-0.00)} +} + +func ExampleBenchmarkBEC() { + linearBlock, _ := hamming.New(context.Background(), 3, 0) + + createMessage := func(trial int) mat.SparseVector { + t := trial % 64 + message := mat.CSRVec(4) + for i := 0; i < 4; i++ { + message.Set(i, (t&(1<>i) + } + return message + } + encode := func(message mat.SparseVector) (codeword []bec.ErasureBit) { + return BitsToErased(linearBlock.Encode(message)) + } + + channel := func(originalCodeword []bec.ErasureBit) (erroredCodeword []bec.ErasureBit) { + //since hamming can fix only one bit wrong under BSC but for BEC this code can fix 2 errors!! + return RandomEraseCount(originalCodeword, 2) + } + + repair := func(originalCodeword, channelInducedCodeword []bec.ErasureBit) (fixed []bec.ErasureBit) { + alg := &iterative.Simple{ + H: linearBlock.H, + } + return bec.Flipping(alg, channelInducedCodeword) + } + + metrics := func(originalMessage mat.SparseVector, originalCodeword, fixedChannelInducedCodeword []bec.ErasureBit) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) { + codewordErrors := ErasedCount(fixedChannelInducedCodeword) + message := linearBlock.DecodeBE(fixedChannelInducedCodeword) + messageErrors := ErasedCount(message) + parityErrors := codewordErrors - messageErrors + + percentFixedCodewordErrors = float64(codewordErrors) / float64(linearBlock.CodewordLength()) + percentFixedMessageErrors = float64(messageErrors) / float64(linearBlock.MessageLength()) + percentFixedParityErrors = float64(parityErrors) / float64(linearBlock.ParitySymbols()) + return + } + + checkpoint := func(updatedStats Stats) { + } + + stats := BenchmarkBEC(context.Background(), 10000, 1, createMessage, encode, channel, repair, metrics, checkpoint, false) + + fmt.Println("Bit Error Probability :", stats) + //Output: + // Bit Error Probability : {Codeword:0.00(+/-0.00), Message:0.00(+/-0.00), Parity:0.00(+/-0.00)} +} + +func ExampleBenchmarkBPSK() { + threads := runtime.NumCPU() + linearBlock, _ := hamming.New(context.Background(), 3, threads) + + createMessage := func(trial int) mat.SparseVector { + t := trial % 64 + message := mat.CSRVec(4) + for i := 0; i < 4; i++ { + message.Set(i, (t&(1<>i) + } + return message + } + + encode := func(message mat.SparseVector) (codeword mat2.Vector) { + return BitsToBPSK(linearBlock.Encode(message)) + } + + channel := func(codeword mat2.Vector) (channelInducedCodeword mat2.Vector) { + //since hamming can fix only one bit wrong, should have zero errors around 2Eb but will sometimes fail do to rounding + return RandomNoiseBPSK(codeword, 2.0) + } + + repair := func(originalCodeword, channelInducedCodeword mat2.Vector) (codeword mat2.Vector) { + //we're going to simulate a hard decision of >=0 is 1 + // and <0 will be 0 on the output codeword + + tmp := mat.CSRVec(channelInducedCodeword.Len()) + for i := 0; i < channelInducedCodeword.Len(); i++ { + if channelInducedCodeword.AtVec(i) >= 0 { + tmp.Set(i, 1) + } + } + alg := &harddecision.Gallager{H: linearBlock.H} + //next we'll do the simple Gallager hard decision bit flipping with a max of 20 iterations + tmp = harddecision.BitFlipping(alg, linearBlock.H, tmp, 20) + return BitsToBPSK(tmp) + } + + metrics := func(message mat.SparseVector, originalCodeword, fixedChannelInducedCodeword mat2.Vector) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) { + codewordErrors := HammingDistanceBPSK(originalCodeword, fixedChannelInducedCodeword) + decoded := linearBlock.Decode(BPSKToBits(fixedChannelInducedCodeword, 0)) + messageErrors := decoded.HammingDistance(message) + parityErrors := codewordErrors - messageErrors + + percentFixedCodewordErrors = float64(codewordErrors) / float64(linearBlock.CodewordLength()) + percentFixedMessageErrors = float64(messageErrors) / float64(linearBlock.MessageLength()) + percentFixedParityErrors = float64(parityErrors) / float64(linearBlock.ParitySymbols()) + return + } + + checkpoint := func(updatedStats Stats) {} + + stats := BenchmarkBPSK(context.Background(), 100_000, threads, createMessage, encode, channel, repair, metrics, checkpoint, false) + + fmt.Println("Bit Error Probability :", stats) + //Output: + // Bit Error Probability : {Codeword:0.00(+/-0.04), Message:0.00(+/-0.05), Parity:0.00(+/-0.05)} +} diff --git a/benchmarking/randomdata.go b/benchmarking/randomdata.go new file mode 100644 index 0000000..bdf605b --- /dev/null +++ b/benchmarking/randomdata.go @@ -0,0 +1,83 @@ +package benchmarking + +import ( + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + mat "github.com/nathanhack/sparsemat" + mat2 "gonum.org/v1/gonum/mat" + "math" + "math/rand" +) + +//RandomMessage creates a random message of length len. +func RandomMessage(len int) mat.SparseVector { + message := mat.CSRVec(len) + for i := 0; i < len; i++ { + message.Set(i, rand.Intn(2)) + } + return message +} + +//RandomMessage creates a random message o lenght len with a hamming weight equal to onesCount +func RandomMessageOnesCount(len int, onesCount int) mat.SparseVector { + message := mat.CSRVec(len) + for message.HammingWeight() < onesCount { + message.Set(rand.Intn(len), 1) + } + return message +} + +//RandomFlipBitCount randomly flips min(numberOfBitsToFlip,len(input)) number of bits. +func RandomFlipBitCount(input mat.SparseVector, numberOfBitsToFlip int) mat.SparseVector { + output := mat.CSRVecCopy(input) + + flip := make(map[int]bool) + for len(flip) < numberOfBitsToFlip && len(flip) < input.Len() { + flip[rand.Intn(input.Len())] = true + } + + for i := range flip { + output.Set(i, output.At(i)+1) + } + return output +} + +//RandomErase creates a new slice of ErasureBits with some of them set to Erased given the probabilityOfErasure +func RandomErase(codeword []bec.ErasureBit, probabilityOfErasure float64) []bec.ErasureBit { + return RandomEraseCount(codeword, int(math.Round(probabilityOfErasure*float64(len(codeword))))) +} + +//RandomErase creates a copy of the codeword and randomly sets numberOfBitsToFlip of them to Erased +func RandomEraseCount(codeword []bec.ErasureBit, numberOfBitsToFlip int) []bec.ErasureBit { + output := make([]bec.ErasureBit, len(codeword)) + + //randomly pick indices to erase + flip := make(map[int]bool) + for len(flip) < numberOfBitsToFlip { + flip[rand.Intn(len(codeword))] = true + } + + //copy the old data + for i := range codeword { + output[i] = codeword[i] + } + + //set the erased + for i := range flip { + output[i] = bec.Erased + } + + return output +} + +//RandomNoiseBPSK creates a randomizes version of the bpsk vector using the E_b/N_0 passed in +func RandomNoiseBPSK(bpsk mat2.Vector, E_bPerN_0 float64) mat2.Vector { + //using σ^2 = N_0/2 and E_b=1 + // we get σ = sqrt(1/(2*E_bPerN_0)) + σ := math.Sqrt(1 / (2 * E_bPerN_0)) + result := mat2.NewVecDense(bpsk.Len(), nil) + for i := 0; i < bpsk.Len(); i++ { + result.SetVec(i, rand.NormFloat64()*σ) + } + result.AddVec(result, bpsk) + return result +} diff --git a/cmd/create.go b/cmd/create.go new file mode 100644 index 0000000..3c87a5a --- /dev/null +++ b/cmd/create.go @@ -0,0 +1,72 @@ +package cmd + +import ( + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/create/gallager" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/create/hamming" + "github.com/spf13/cobra" +) + +// createCmd represents the create command +var createCmd = &cobra.Command{ + Use: "create", + Aliases: []string{"c"}, + Short: "used to create a new ECC", + Long: `create provides the ability to make a new ECC from the list of built-in ECCs and save them so they can be used later by the tools.`, +} + +// createlinearblockCmd represents the linearblock command +var createlinearblockCmd = &cobra.Command{ + Use: "linearblock", + Aliases: []string{"lb", "l"}, + Short: "creates linearblock ECCs", + Long: `Creates linearblock ECCs.`, +} + +// createldpcCmd represents the ldpc command +var createldpcCmd = &cobra.Command{ + Use: "ldpc", + Aliases: []string{"l"}, + Short: "creates LDPC", + Long: `Creates linearblock ECCs known as Low Density Parity Check (LDPC)`, +} + +// createGallagerCmd represents the gallager command +var createGallagerCmd = &cobra.Command{ + Use: "gallager OUTPUT_LDPC_JSON", + Aliases: []string{"g"}, + Short: "Creates a new Gallager based ECC", + Long: `Creates a new Gallager based ECC. Note a small cycle has a negative effect on the effectiveness of the LDPC.`, + Args: cobra.ExactArgs(1), + Run: gallager.GallagerRun, +} + +//createHammingCmd represents the Hamming command +var createHammingCmd = &cobra.Command{ + Use: "hamming OUTPUT_HAMMING_JSON", + Aliases: []string{"h", "ham"}, + Short: "Creates a new Hamming code based ECC", + Long: `Creates a new Hamming code based ECC.`, + Args: cobra.ExactArgs(1), + Run: hamming.HammingRun, +} + +func init() { + rootCmd.AddCommand(createCmd) + createCmd.AddCommand(createlinearblockCmd) + createlinearblockCmd.AddCommand(createldpcCmd) + createldpcCmd.AddCommand(createGallagerCmd) + + createGallagerCmd.Flags().UintVarP(&gallager.Message, "message", "m", 1000, "the number of bits in the message") + createGallagerCmd.Flags().UintVarP(&gallager.Wc, "column", "c", 3, "the column weight (number of ones in the H matrix column) (>=3)") + createGallagerCmd.Flags().UintVarP(&gallager.Wr, "row", "r", 4, "the row weight (number of ones in the H matrix row) (column < row)") + createGallagerCmd.Flags().UintVarP(&gallager.Smallest, "smallest", "s", 4, "the smallest allowed cycle: 4, 6, 8...") + createGallagerCmd.Flags().UintVarP(&gallager.Iter, "iter", "i", 10000, "the number of iterations to try before terminating the search") + createGallagerCmd.Flags().UintVarP(&gallager.Threads, "threads", "t", 0, "the number of threads to use; note 0 means use the number of cpus") + createGallagerCmd.Flags().BoolVarP(&gallager.Verbose, "verbose", "v", false, "enable verbose info") + + createlinearblockCmd.AddCommand(createHammingCmd) + + createHammingCmd.Flags().UintVarP(&hamming.ParityBits, "parity", "p", 4, "the parity >=2, sets codeword size (cs) == 2^parity-1 and message size == cs-parity") + createHammingCmd.Flags().UintVarP(&hamming.Threads, "threads", "t", 0, "the number of threads to use; note 0 means use the number of cpus") + +} diff --git a/cmd/internal/create/gallager/gallager.go b/cmd/internal/create/gallager/gallager.go new file mode 100644 index 0000000..18ce4af --- /dev/null +++ b/cmd/internal/create/gallager/gallager.go @@ -0,0 +1,66 @@ +package gallager + +import ( + "context" + "encoding/json" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/linearblock/ldpc/gallager" + "github.com/sirupsen/logrus" + "github.com/spf13/cobra" + "io/ioutil" + "math/rand" + "os" + "os/signal" + "syscall" + "time" +) + +var Message uint +var Wc uint +var Wr uint +var Smallest uint +var Iter uint +var Threads uint +var Verbose bool + +var GallagerRun = func(cmd *cobra.Command, args []string) { + //we seed the randomizer so we get something different every time + rand.Seed(time.Now().Unix()) + + if Verbose { + logrus.SetLevel(logrus.DebugLevel) + } else { + logrus.SetLevel(logrus.InfoLevel) + } + sigs := make(chan os.Signal, 1) + signal.Notify(sigs, syscall.SIGINT, syscall.SIGTERM) + ctx, cancel := context.WithCancel(context.Background()) + go func() { + sig := <-sigs + fmt.Println() + fmt.Println(sig) + cancel() + }() + + g, err := gallager.Search(ctx, int(Message), int(Wc), int(Wr), int(Smallest), int(Iter), int(Threads)) + if err != nil { + fmt.Println("Unable to create gallager LDPC: ", err) + return + } + + if g == nil { + fmt.Println("Unable to create gallager LDPC try again") + return + } + + bs, err := json.Marshal(g) + if err != nil { + fmt.Println("Unable to serialize the gallager LDPC: ", err) + return + } + + err = ioutil.WriteFile(args[0], bs, 0644) + if err != nil { + fmt.Println("unable to write file: ", err) + } +} diff --git a/cmd/internal/create/hamming/hamming.go b/cmd/internal/create/hamming/hamming.go new file mode 100644 index 0000000..5cf1072 --- /dev/null +++ b/cmd/internal/create/hamming/hamming.go @@ -0,0 +1,63 @@ +package hamming + +import ( + "context" + "encoding/json" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/linearblock/hamming" + "github.com/sirupsen/logrus" + "github.com/spf13/cobra" + "io/ioutil" + "math/rand" + "os" + "os/signal" + "syscall" + "time" +) + +var ( + ParityBits uint + Threads uint + Verbose bool +) +var HammingRun = func(cmd *cobra.Command, args []string) { + //we seed the randomizer so we get something different every time + rand.Seed(time.Now().Unix()) + + if Verbose { + logrus.SetLevel(logrus.DebugLevel) + } else { + logrus.SetLevel(logrus.InfoLevel) + } + sigs := make(chan os.Signal, 1) + signal.Notify(sigs, syscall.SIGINT, syscall.SIGTERM) + ctx, cancel := context.WithCancel(context.Background()) + go func() { + sig := <-sigs + fmt.Println() + fmt.Println(sig) + cancel() + }() + + g, err := hamming.New(ctx, int(ParityBits), int(Threads)) + if err != nil { + fmt.Println("Unable to create gallager LDPC: ", err) + return + } + + if g == nil { + fmt.Println("Unable to create gallager LDPC try again") + return + } + + bs, err := json.Marshal(g) + if err != nil { + fmt.Println("Unable to serialize the gallager LDPC: ", err) + return + } + + err = ioutil.WriteFile(args[0], bs, 0644) + if err != nil { + fmt.Println("unable to write file: ", err) + } +} diff --git a/cmd/internal/tools/bec/bec.go b/cmd/internal/tools/bec/bec.go new file mode 100644 index 0000000..fa699b6 --- /dev/null +++ b/cmd/internal/tools/bec/bec.go @@ -0,0 +1,74 @@ +package bec + +import ( + "context" + "github.com/nathanhack/errorcorrectingcodes/benchmarking" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + mat "github.com/nathanhack/sparsemat" + "math" + "math/rand" + "sync" +) + +func RunBEC(ctx context.Context, + l *linearblock.LinearBlock, + percentage float64, trials, threads int, + correctionAlg benchmarking.BinaryErasureChannelCorrection, + previousStats benchmarking.Stats, + checkpoints benchmarking.Checkpoints, + showProgressBar bool) benchmarking.Stats { + messageHistory := make(map[string]bool) + messageHistoryMux := sync.RWMutex{} + messageHistoryMax := math.Pow(2, float64(l.MessageLength())) + + createMessage := func(trial int) mat.SparseVector { + message := mat.CSRVec(l.MessageLength()) + messageHistoryMux.RLock() + _, has := messageHistory[message.String()] + messageHistoryMux.RUnlock() + for has { + reset := false + for i := 0; i < l.MessageLength(); i++ { + message.Set(i, rand.Intn(2)) + } + messageHistoryMux.RLock() + _, has = messageHistory[message.String()] + reset = float64(len(messageHistory)) >= messageHistoryMax + messageHistoryMux.RUnlock() + + if reset { + messageHistoryMux.Lock() + messageHistory = make(map[string]bool) + messageHistoryMux.Unlock() + } + } + messageHistoryMux.Lock() + messageHistory[message.String()] = true + messageHistoryMux.Unlock() + return message + } + + encode := func(message mat.SparseVector) (codeword []bec.ErasureBit) { + return l.EncodeBE(message) + } + + channel := func(originalCodeword []bec.ErasureBit) (erroredCodeword []bec.ErasureBit) { + count := int(percentage * float64(len(originalCodeword))) + return benchmarking.RandomEraseCount(originalCodeword, count) + } + + metrics := func(originalMessage mat.SparseVector, originalCodeword, fixedChannelInducedCodeword []bec.ErasureBit) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) { + codewordErrors := benchmarking.ErasedCount(fixedChannelInducedCodeword) + message := l.DecodeBE(fixedChannelInducedCodeword) + messageErrors := benchmarking.ErasedCount(message) + parityErrors := codewordErrors - messageErrors + + percentFixedCodewordErrors = float64(codewordErrors) / float64(l.CodewordLength()) + percentFixedMessageErrors = float64(messageErrors) / float64(l.MessageLength()) + percentFixedParityErrors = float64(parityErrors) / float64(l.ParitySymbols()) + return + } + + return benchmarking.BenchmarkBECContinueStats(ctx, trials, threads, createMessage, encode, channel, correctionAlg, metrics, checkpoints, previousStats, showProgressBar) +} diff --git a/cmd/internal/tools/bec/simple/simple.go b/cmd/internal/tools/bec/simple/simple.go new file mode 100644 index 0000000..7d29417 --- /dev/null +++ b/cmd/internal/tools/bec/simple/simple.go @@ -0,0 +1,145 @@ +package simple + +import ( + "context" + "fmt" + "github.com/cheggaaa/pb/v3" + "github.com/nathanhack/errorcorrectingcodes/benchmarking" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/bec" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + bec2 "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec/iterative" + "github.com/spf13/cobra" + "os" + "os/signal" + "reflect" + "runtime" + "sync" + "syscall" +) + +var ( + Trials uint + ErrorProbability []float64 + Threads uint +) +var BecRun = func(cmd *cobra.Command, args []string) { + if len(args) != 2 { + fmt.Println("requires both ECC_JSON_FILE RESULT_JSON") + return + } + + //first get the ECC to use + ecc, err := tools.LoadLinearBlockECC(args[0]) + if err != nil { + fmt.Println(err) + return + } + + //next we see if the RESULT_JSON exists if so we load it and validate we're running it against the right thing + data, err := tools.LoadResults(args[1]) + if err != nil { + fmt.Println(err) + return + } + + //if data is nil then we create it + if data == nil { + data = &tools.SimulationStats{ + TypeInfo: typeInfo(), + ECCInfo: tools.Md5Sum(ecc.H), + Stats: make(map[float64]benchmarking.Stats), + } + } + + //in either case lets validate it + if data.TypeInfo != typeInfo() { + fmt.Printf("csv loaded does not match the same type expected %v but found %v\n", typeInfo(), data.TypeInfo) + return + } + if data.ECCInfo != tools.Md5Sum(ecc.H) { + fmt.Printf("csv laoded does not match the ECC") + return + } + + // handle ctrl-C's to kill in a nice way + sigs := make(chan os.Signal, 1) + signal.Notify(sigs, syscall.SIGINT, syscall.SIGTERM) + ctx, cancel := context.WithCancel(context.Background()) + go func() { + sig := <-sigs + fmt.Println() + fmt.Println(sig) + cancel() + }() + + runSimulation(ctx, data, ecc, args[1]) + + err = tools.SaveResults(args[1], data) + if err != nil { + fmt.Println(err) + } +} + +func typeInfo() string { + t := reflect.TypeOf(iterative.Simple{}) + return fmt.Sprintf("BEC:%v/%v", t.PkgPath(), t.Name()) + +} + +func min(a, b int) int { + if a < b { + return a + } + return b +} + +func runSimulation(ctx context.Context, data *tools.SimulationStats, ecc *linearblock.LinearBlock, outputFilename string) { + checkpointMux := sync.Mutex{} + checkpointCount := 0 + + alg := &iterative.Simple{ + H: ecc.H, + } + correctionAlg := func(originalCodeword, channelInducedCodeword []bec2.ErasureBit) (fixedChannelInducedCodeword []bec2.ErasureBit) { + return bec2.Flipping(alg, channelInducedCodeword) + } + + numberOfThread := int(Threads) + if numberOfThread == 0 { + numberOfThread = runtime.NumCPU() + } + + trialsPerIter := numberOfThread * 10 + bar := pb.StartNew(int(Trials) * len(ErrorProbability)) +trialLoops: + for t := 0; t <= int(Trials); t += trialsPerIter { + select { + case <-ctx.Done(): + break trialLoops + default: + } + + for _, p := range ErrorProbability { + checkpoint := func(stats benchmarking.Stats) { + //we want to save the checkpoint + checkpointMux.Lock() + defer checkpointMux.Unlock() + + data.Stats[p] = stats + + if checkpointCount%trialsPerIter == 0 { + err := tools.SaveResults(outputFilename, data) + if err != nil { + fmt.Println(err) + } + } + checkpointCount++ + } + data.Stats[p] = bec.RunBEC(ctx, ecc, p, min(t, int(Trials)), numberOfThread, correctionAlg, data.Stats[p], checkpoint, false) + bar.Add(trialsPerIter) + } + } + bar.Finish() +} diff --git a/cmd/internal/tools/bsc/bsc.go b/cmd/internal/tools/bsc/bsc.go new file mode 100644 index 0000000..64a9538 --- /dev/null +++ b/cmd/internal/tools/bsc/bsc.go @@ -0,0 +1,73 @@ +package bsc + +import ( + "context" + "github.com/nathanhack/errorcorrectingcodes/benchmarking" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + mat "github.com/nathanhack/sparsemat" + "math" + "math/rand" + "sync" +) + +func RunBSC(ctx context.Context, + l *linearblock.LinearBlock, + crossoverProbability float64, trials, threads int, + correctionAlg benchmarking.BinarySymmetricChannelCorrection, + previousStats benchmarking.Stats, + checkpoints benchmarking.Checkpoints, + showProgress bool) benchmarking.Stats { + messageHistory := make(map[string]bool) + messageHistoryMux := sync.RWMutex{} + messageHistoryMax := math.Pow(2, float64(l.MessageLength())) + + createMessage := func(trial int) mat.SparseVector { + message := mat.CSRVec(l.MessageLength()) + messageHistoryMux.RLock() + _, has := messageHistory[message.String()] + messageHistoryMux.RUnlock() + for has { + reset := false + for i := 0; i < l.MessageLength(); i++ { + message.Set(i, rand.Intn(2)) + } + messageHistoryMux.RLock() + _, has = messageHistory[message.String()] + reset = float64(len(messageHistory)) >= messageHistoryMax + messageHistoryMux.RUnlock() + + if reset { + messageHistoryMux.Lock() + messageHistory = make(map[string]bool) + messageHistoryMux.Unlock() + } + } + messageHistoryMux.Lock() + messageHistory[message.String()] = true + messageHistoryMux.Unlock() + return message + } + + encode := func(message mat.SparseVector) (codeword mat.SparseVector) { + return l.Encode(message) + } + + channel := func(originalCodeword mat.SparseVector) (erroredCodeword mat.SparseVector) { + count := int(crossoverProbability * float64(originalCodeword.Len())) + return benchmarking.RandomFlipBitCount(originalCodeword, count) + } + + metrics := func(originalMessage, originalCodeword, fixedChannelInducedCodeword mat.SparseVector) (percentFixedCodewordErrors, percentFixedMessageErrors, percentFixedParityErrors float64) { + codewordErrors := originalCodeword.HammingDistance(fixedChannelInducedCodeword) + message := l.Decode(fixedChannelInducedCodeword) + messageErrors := message.HammingDistance(originalMessage) + parityErrors := codewordErrors - messageErrors + + percentFixedCodewordErrors = float64(codewordErrors) / float64(l.CodewordLength()) + percentFixedMessageErrors = float64(messageErrors) / float64(l.MessageLength()) + percentFixedParityErrors = float64(parityErrors) / float64(l.ParitySymbols()) + return + } + + return benchmarking.BenchmarkBSCContinueStats(ctx, trials, threads, createMessage, encode, channel, correctionAlg, metrics, checkpoints, previousStats, showProgress) +} diff --git a/cmd/internal/tools/bsc/dwbf/dwbf.go b/cmd/internal/tools/bsc/dwbf/dwbf.go new file mode 100644 index 0000000..5382df7 --- /dev/null +++ b/cmd/internal/tools/bsc/dwbf/dwbf.go @@ -0,0 +1,152 @@ +package dwbf + +import ( + "context" + "fmt" + "github.com/cheggaaa/pb/v3" + "github.com/nathanhack/errorcorrectingcodes/benchmarking" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/bsc" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision" + mat "github.com/nathanhack/sparsemat" + "github.com/spf13/cobra" + "os" + "os/signal" + "reflect" + "runtime" + "sync" + "syscall" +) + +var ( + Trials uint + ErrorProbability []float64 + Threads uint + Verbose bool + MaxIter uint + Alpha float64 + EtaThreshold float64 +) + +var DwbfRun = func(cmd *cobra.Command, args []string) { + if len(args) != 2 { + fmt.Println("requires both ECC_JSON_FILE RESULT_JSON") + return + } + + //first get the ECC to use + ecc, err := tools.LoadLinearBlockECC(args[0]) + if err != nil { + fmt.Println(err) + return + } + + //next we see if the RESULT_JSON exists if so we load it and validate we're running it against the right thing + var data *tools.SimulationStats + data, err = tools.LoadResults(args[1]) + if err != nil { + fmt.Println(err) + return + } + + //if data is nil then we create it + if data == nil { + data = &tools.SimulationStats{ + TypeInfo: typeInfo(), + ECCInfo: tools.Md5Sum(ecc.H), + Stats: make(map[float64]benchmarking.Stats), + } + } + + //in either case lets validate it + if data.TypeInfo != typeInfo() { + fmt.Printf("csv loaded does not match the same type expected %v but found %v\n", typeInfo(), data.TypeInfo) + return + } + if data.ECCInfo != tools.Md5Sum(ecc.H) { + fmt.Printf("csv laoded does not match the ECC") + return + } + + sigs := make(chan os.Signal, 1) + signal.Notify(sigs, syscall.SIGINT, syscall.SIGTERM) + ctx, cancel := context.WithCancel(context.Background()) + go func() { + sig := <-sigs + fmt.Println() + fmt.Println(sig) + cancel() + }() + + runSimulation(ctx, data, ecc, args[1]) + + err = tools.SaveResults(args[1], data) + if err != nil { + fmt.Println(err) + } +} + +func typeInfo() string { + t := reflect.TypeOf(harddecision.DWBF_F{}) + return fmt.Sprintf("BSC:%v/%v", t.PkgPath(), t.Name()) +} + +func min(a, b int) int { + if a < b { + return a + } + return b +} + +func runSimulation(ctx context.Context, data *tools.SimulationStats, ecc *linearblock.LinearBlock, outputFilename string) { + checkpointMux := sync.Mutex{} + checkpointCount := 0 + + correctionAlg := func(originalCodeword, channelInducedCodeword mat.SparseVector) (fixedChannelInducedCodeword mat.SparseVector) { + //since this is parallel there is no way to isolate data from one codeword from the next + // this alg has internal state + alg := &harddecision.DWBF_F{ + AlphaFactor: Alpha, + H: ecc.H, + } + return harddecision.BitFlipping(alg, ecc.H, channelInducedCodeword, int(MaxIter)) + } + + numberOfThread := int(Threads) + if numberOfThread == 0 { + numberOfThread = runtime.NumCPU() + } + + trialsPerIter := numberOfThread * 10 + bar := pb.StartNew(int(Trials) * len(ErrorProbability)) +trialLoops: + for t := 0; t <= int(Trials); t += trialsPerIter { + select { + case <-ctx.Done(): + break trialLoops + default: + } + + for _, p := range ErrorProbability { + checkpoint := func(stats benchmarking.Stats) { + //we want to save the checkpoint + checkpointMux.Lock() + defer checkpointMux.Unlock() + + data.Stats[p] = stats + + if checkpointCount%trialsPerIter == 0 { + err := tools.SaveResults(outputFilename, data) + if err != nil { + fmt.Println(err) + } + } + checkpointCount++ + } + data.Stats[p] = bsc.RunBSC(ctx, ecc, p, min(t, int(Trials)), numberOfThread, correctionAlg, data.Stats[p], checkpoint, false) + bar.Add(trialsPerIter) + } + } + bar.Finish() +} diff --git a/cmd/internal/tools/bsc/gallager/gallager.go b/cmd/internal/tools/bsc/gallager/gallager.go new file mode 100644 index 0000000..a1784e8 --- /dev/null +++ b/cmd/internal/tools/bsc/gallager/gallager.go @@ -0,0 +1,147 @@ +package gallager + +import ( + "context" + "fmt" + "github.com/cheggaaa/pb/v3" + "github.com/nathanhack/errorcorrectingcodes/benchmarking" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/bsc" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision" + mat "github.com/nathanhack/sparsemat" + "github.com/spf13/cobra" + "os" + "os/signal" + "reflect" + "runtime" + "sync" + "syscall" +) + +var ( + Trials uint + ErrorProbability []float64 + Threads uint + MaxIter uint +) + +var GallagerRun = func(cmd *cobra.Command, args []string) { + if len(args) != 2 { + fmt.Println("requires both ECC_JSON_FILE RESULT_JSON") + return + } + + //first get the ECC to use + ecc, err := tools.LoadLinearBlockECC(args[0]) + if err != nil { + fmt.Println(err) + return + } + + //next we see if the RESULT_JSON exists if so we load it and validate we're running it against the right thing + data, err := tools.LoadResults(args[1]) + if err != nil { + fmt.Println(err) + return + } + + //if data is nil then we create it + if data == nil { + data = &tools.SimulationStats{ + TypeInfo: typeInfo(), + ECCInfo: tools.Md5Sum(ecc.H), + Stats: make(map[float64]benchmarking.Stats), + } + } + + //in either case lets validate it + if data.TypeInfo != typeInfo() { + fmt.Printf("csv loaded does not match the same type expected %v but found %v\n", typeInfo(), data.TypeInfo) + return + } + if data.ECCInfo != tools.Md5Sum(ecc.H) { + fmt.Printf("csv laoded does not match the ECC") + return + } + + sigs := make(chan os.Signal, 1) + signal.Notify(sigs, syscall.SIGINT, syscall.SIGTERM) + ctx, cancel := context.WithCancel(context.Background()) + go func() { + sig := <-sigs + fmt.Println() + fmt.Println(sig) + cancel() + }() + + runSimulation(ctx, data, ecc, args[1]) + + err = tools.SaveResults(args[1], data) + if err != nil { + fmt.Println(err) + } +} + +func typeInfo() string { + t := reflect.TypeOf(harddecision.Gallager{}) + return fmt.Sprintf("BSC:%v/%v", t.PkgPath(), t.Name()) +} + +func min(a, b int) int { + if a < b { + return a + } + return b +} + +func runSimulation(ctx context.Context, data *tools.SimulationStats, ecc *linearblock.LinearBlock, outputFilename string) { + checkpointMux := sync.Mutex{} + checkpointCount := 0 + + alg := &harddecision.Gallager{ + H: ecc.H, + } + + correctionAlg := func(originalCodeword, channelInducedCodeword mat.SparseVector) (fixedChannelInducedCodeword mat.SparseVector) { + alg.Reset() //must be called before correcting the next codeword + return harddecision.BitFlipping(alg, ecc.H, channelInducedCodeword, int(MaxIter)) + } + + numberOfThread := int(Threads) + if numberOfThread == 0 { + numberOfThread = runtime.NumCPU() + } + + trialsPerIter := numberOfThread * 10 + bar := pb.StartNew(int(Trials) * len(ErrorProbability)) +trialLoops: + for t := 0; t <= int(Trials); t += trialsPerIter { + select { + case <-ctx.Done(): + break trialLoops + default: + } + + for _, p := range ErrorProbability { + checkpoint := func(stats benchmarking.Stats) { + //we want to save the checkpoint + checkpointMux.Lock() + defer checkpointMux.Unlock() + + data.Stats[p] = stats + + if checkpointCount%trialsPerIter == 0 { + err := tools.SaveResults(outputFilename, data) + if err != nil { + fmt.Println(err) + } + } + checkpointCount++ + } + data.Stats[p] = bsc.RunBSC(ctx, ecc, p, min(t, int(Trials)), numberOfThread, correctionAlg, data.Stats[p], checkpoint, false) + bar.Add(trialsPerIter) + } + } + bar.Finish() +} diff --git a/cmd/internal/tools/csv/csv.go b/cmd/internal/tools/csv/csv.go new file mode 100644 index 0000000..adbe3ab --- /dev/null +++ b/cmd/internal/tools/csv/csv.go @@ -0,0 +1,90 @@ +package csv + +import ( + "encoding/csv" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools" + "github.com/spf13/cobra" + "os" + "path/filepath" + "sort" + "strings" +) + +var OutputFile string +var MessageError bool +var ParityError bool + +var CSVRun = func(cmd *cobra.Command, args []string) { + if len(args) < 1 { + fmt.Println("requires at least one RESULTS_JSON") + return + } + + stats := make([]*tools.SimulationStats, len(args)) + var err error + percentagesFloats := make(map[float64]bool) + for i, resultFile := range args { + stats[i], err = tools.LoadResults(resultFile) + if err != nil { + fmt.Println(err) + return + } + for p := range stats[i].Stats { + percentagesFloats[p] = true + } + } + + f, err := os.Create(OutputFile) + if err != nil { + fmt.Println(err) + return + } + defer f.Close() + w := csv.NewWriter(f) + defer w.Flush() + + //first write headers + percentagesList := make([]float64, 0, len(percentagesFloats)) + for p := range percentagesFloats { + percentagesList = append(percentagesList, p) + } + sort.Float64s(percentagesList) + + header := []string{"Results File"} + + for _, p := range percentagesList { + header = append(header, fmt.Sprintf("%v", p)) + } + + err = w.Write(header) + if err != nil { + fmt.Println(err) + return + } + + for i, s := range stats { + record := make([]string, len(header)) + record[0] = strings.TrimSuffix(args[i], filepath.Ext(args[i])) + + for i, p := range percentagesList { + v, has := s.Stats[p] + if has { + switch { + case MessageError: + record[i+1] = fmt.Sprintf("%v", v.ChannelMessageError.Mean) + case ParityError: + record[i+1] = fmt.Sprintf("%v", v.ChannelParityError.Mean) + default: + record[i+1] = fmt.Sprintf("%v", v.ChannelCodewordError.Mean) + } + } + } + + err = w.Write(record) + if err != nil { + fmt.Println(err) + return + } + } +} diff --git a/cmd/internal/tools/tools.go b/cmd/internal/tools/tools.go new file mode 100644 index 0000000..0ea2103 --- /dev/null +++ b/cmd/internal/tools/tools.go @@ -0,0 +1,114 @@ +package tools + +import ( + "crypto/md5" + "encoding/json" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/benchmarking" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + mat "github.com/nathanhack/sparsemat" + "io/ioutil" + "os" + "strconv" +) + +type SimulationStats struct { + TypeInfo string + ECCInfo string + Stats map[float64]benchmarking.Stats +} +type simulationStats struct { + TypeInfo string + ECCInfo string + Stats map[string]benchmarking.Stats +} + +func (s *SimulationStats) MarshalJSON() ([]byte, error) { + ss := simulationStats{ + TypeInfo: s.TypeInfo, + ECCInfo: s.ECCInfo, + Stats: map[string]benchmarking.Stats{}, + } + + for f, stat := range s.Stats { + ss.Stats[fmt.Sprintf("%v", f)] = stat + } + + return json.Marshal(ss) +} + +func (s *SimulationStats) UnmarshalJSON(bytes []byte) error { + var ss simulationStats + + err := json.Unmarshal(bytes, &ss) + if err != nil { + return err + } + + s.TypeInfo = ss.TypeInfo + s.ECCInfo = ss.ECCInfo + s.Stats = map[float64]benchmarking.Stats{} + + for fs, stat := range ss.Stats { + f, err := strconv.ParseFloat(fs, 64) + if err != nil { + return err + } + s.Stats[f] = stat + } + return nil +} + +func Md5Sum(H mat.SparseMat) string { + return fmt.Sprintf("%x", md5.Sum([]byte(H.String()))) +} + +func LoadLinearBlockECC(filepath string) (*linearblock.LinearBlock, error) { + if _, err := os.Stat(filepath); os.IsNotExist(err) { + return nil, fmt.Errorf("the ECC_JSON_FILE must exist") + } + + bs, err := ioutil.ReadFile(filepath) + if err != nil { + return nil, fmt.Errorf("error while reading file %v: %v\n", filepath, err) + } + + var ecc linearblock.LinearBlock + err = json.Unmarshal(bs, &ecc) + if err != nil { + return nil, fmt.Errorf("error while reading file %v: %v\n", filepath, err) + } + + return &ecc, nil +} + +func LoadResults(filepath string) (*SimulationStats, error) { + if _, err := os.Stat(filepath); os.IsNotExist(err) { + return nil, nil + } + + bs, err := ioutil.ReadFile(filepath) + if err != nil { + return nil, fmt.Errorf("error while reading file %v: %v\n", filepath, err) + } + + var stat SimulationStats + err = json.Unmarshal(bs, &stat) + if err != nil { + return nil, fmt.Errorf("error while unmarshalling file %v: %v\n", filepath, err) + } + return &stat, nil +} + +func SaveResults(filepath string, data *SimulationStats) error { + bs, err := json.Marshal(data) + if err != nil { + return fmt.Errorf("error serializing csv: %v\n", err) + } + + err = ioutil.WriteFile(filepath, bs, 0644) + if err != nil { + fmt.Errorf("error while saving csv to %v: %v\n", filepath, err) + } + return nil +} diff --git a/cmd/root.go b/cmd/root.go new file mode 100644 index 0000000..fb28b5e --- /dev/null +++ b/cmd/root.go @@ -0,0 +1,18 @@ +package cmd + +import ( + "github.com/spf13/cobra" +) + +// rootCmd represents the base command when called without any subcommands +var rootCmd = &cobra.Command{ + Use: "errorcorrectingcodes", + Short: "A suite of tools for ECCs", + Long: `A suite of tools for ECCs`, +} + +// Execute adds all child commands to the root command and sets flags appropriately. +// This is called by main.main(). It only needs to happen once to the rootCmd. +func Execute() { + cobra.CheckErr(rootCmd.Execute()) +} diff --git a/cmd/tools.go b/cmd/tools.go new file mode 100644 index 0000000..956cedc --- /dev/null +++ b/cmd/tools.go @@ -0,0 +1,128 @@ +package cmd + +import ( + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/bec/simple" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/bsc/dwbf" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/bsc/gallager" + "github.com/nathanhack/errorcorrectingcodes/cmd/internal/tools/csv" + + "github.com/spf13/cobra" +) + +// toolsCmd represents the tools command +var toolsCmd = &cobra.Command{ + Use: "tools", + Aliases: []string{"t"}, + Short: "Tools for ECCs", + Long: `Tools for ECCs`, +} + +// toolsChansimCmd represents the chansim command +var toolsChansimCmd = &cobra.Command{ + Use: "chansim", + Aliases: []string{"cs", "c"}, + Short: "Channel simulators", + Long: `Channel simulators for linearblock ECCs`, +} + +// toolsLinearblockCmd represents the linearblock command +var toolsLinearblockCmd = &cobra.Command{ + Use: "linearblock", + Aliases: []string{"lb", "l"}, + Short: "Linearblock channel simulators", + Long: `Channel simulators for linearblock ECCs`, +} + +// toolsHarddecisionCmd represents the harddecision command +var toolsHarddecisionCmd = &cobra.Command{ + Use: "harddecision", + Aliases: []string{"hard", "h"}, + Short: "Using hard decisions", + Long: `Channel simulators for linearblock ECCs using hard decisions`, +} + +// toolsBecCmd represents the bec command +var toolsBecCmd = &cobra.Command{ + Use: "bec ECC_JSON_FILE RESULT_JSON", + Short: "An erasure channel simulator", + Long: `A simple erasure channel simulator for linearblock ECCs`, + Run: simple.BecRun, +} + +// toolsBscCmd represents the bsc command +var toolsBscCmd = &cobra.Command{ + Use: "bsc", + Short: "A binary symmetric channel simulator", + Long: `A binary symmetric channel simulator for linearblock ECCs`, +} + +// toolsDwbfCmd represents the dwbf command +var toolsDwbfCmd = &cobra.Command{ + Use: "dwbf ECC_JSON_FILE RESULT_JSON", + Aliases: []string{"d"}, + Short: "A linearblock BSC simulator with dwbf based bit flipping algorithm", + Long: `A linearblock BSC simulator with dwbf based bit flipping algorithm`, + Run: dwbf.DwbfRun, +} + +// toolsGallagerCmd represents the gallager command +var toolsGallagerCmd = &cobra.Command{ + Use: "gallager ECC_JSON_FILE RESULT_JSON", + Aliases: []string{"g"}, + Short: "A linearblock BSC simulator with gallager based bit flipping algorithm", + Long: `A linearblock BSC simulator with gallager based bit flipping algorithm`, + Run: gallager.GallagerRun, +} + +// toolsResultsCmd represents the csv command +var toolsResultsCmd = &cobra.Command{ + Use: "results", + Aliases: []string{"r"}, + Short: "A tool to organize results for graphing and comparison", + Long: `A tool to organize results for graphing and comparison`, +} + +//toolsCSVCmd represents the csv command +var toolsCSVCmd = &cobra.Command{ + Use: "csv RESULTS_JSON [RESULTS_JSON] ...", + Aliases: []string{"c"}, + Short: "Export to a CSV file", + Long: `Export to a CSV file`, + Run: csv.CSVRun, +} + +func init() { + rootCmd.AddCommand(toolsCmd) + toolsCmd.AddCommand(toolsChansimCmd) + toolsCmd.AddCommand(toolsResultsCmd) + + toolsChansimCmd.AddCommand(toolsLinearblockCmd) + toolsLinearblockCmd.AddCommand(toolsHarddecisionCmd) + + toolsHarddecisionCmd.AddCommand(toolsBecCmd) + toolsBecCmd.Flags().UintVarP(&simple.Trials, "trials", "t", 1_000_000, "the number of trials per step") + toolsBecCmd.Flags().Float64SliceVarP(&simple.ErrorProbability, "probability", "p", []float64{0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99}, "probability of erasure [0, 1)") + toolsBecCmd.Flags().UintVar(&simple.Threads, "threads", 0, "number of threads to use (0 means to use the # of threads equal to the # of CPUs)") + + toolsHarddecisionCmd.AddCommand(toolsBscCmd) + + toolsBscCmd.AddCommand(toolsDwbfCmd) + toolsDwbfCmd.Flags().UintVarP(&dwbf.Trials, "trials", "t", 1_000_000, "the number of trials per step") + toolsDwbfCmd.Flags().Float64SliceVarP(&dwbf.ErrorProbability, "probability", "p", []float64{0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50}, "probability of crossover errors to test [0, 0.5]") + toolsDwbfCmd.Flags().UintVar(&dwbf.Threads, "threads", 0, "number of threads to use (0 means to use the # of threads equal to the # of CPUs)") + toolsDwbfCmd.Flags().UintVarP(&dwbf.MaxIter, "iters", "i", 20, "max number of iterations the bitflip algorithm is allowed") + toolsDwbfCmd.Flags().Float64VarP(&dwbf.Alpha, "alpha", "a", .5, "hyperparameter 0<α<1") + toolsDwbfCmd.Flags().Float64VarP(&dwbf.EtaThreshold, "eta", "e", 0.0, "hyperparameter η threshold: no requirement but frequently 0.0 is a good value") + + toolsBscCmd.AddCommand(toolsGallagerCmd) + + toolsGallagerCmd.Flags().UintVarP(&gallager.Trials, "trials", "t", 1_000_000, "the number of trials per step") + toolsGallagerCmd.Flags().Float64SliceVarP(&gallager.ErrorProbability, "probability", "p", []float64{0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50}, "probability of crossover errors to test [0, 0.5]") + toolsGallagerCmd.Flags().UintVar(&gallager.Threads, "threads", 0, "number of threads to use (0 means to use the # of threads equal to the # of CPUs)") + toolsGallagerCmd.Flags().UintVarP(&gallager.MaxIter, "iters", "i", 20, "max number of iterations the bitflip algorithm is allowed") + + toolsResultsCmd.AddCommand(toolsCSVCmd) + toolsCSVCmd.Flags().StringVarP(&csv.OutputFile, "output", "o", "results.csv", "filename of the combined csv") + toolsCSVCmd.Flags().BoolVarP(&csv.MessageError, "message", "m", false, "outputs the MessageError instead of CodewordError or ParityError") + toolsCSVCmd.Flags().BoolVarP(&csv.ParityError, "parity", "p", false, "outputs the ParityError instead of CodewordError or MessageError") +} diff --git a/data/linearblock/hamming/hamming_15_11.json b/data/linearblock/hamming/hamming_15_11.json new file mode 100644 index 0000000..6924bd3 --- /dev/null +++ b/data/linearblock/hamming/hamming_15_11.json @@ -0,0 +1 @@ +{"H":{"Rows":4,"Cols":15,"Data":[[0,2,4,6,8,10,12,14],[1,2,5,6,9,10,13,14],[3,4,5,6,11,12,13,14],[7,8,9,10,11,12,13,14]]},"Processing":{"HColumnOrder":[4,5,6,7,8,9,10,11,12,13,2,0,1,14,3],"G":{"Rows":11,"Cols":15,"Data":[[0,11,14],[1,12,14],[2,11,12,14],[3,11,12,13,14],[4,12,13,14],[5,11,13,14],[6,13,14],[7,11,12,13],[8,12,13],[9,11,13],[10,11,12]]}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_15_11_bec_results.json b/data/linearblock/hamming/hamming_15_11_bec_results.json new file mode 100644 index 0000000..6fc8cb0 --- /dev/null +++ b/data/linearblock/hamming/hamming_15_11_bec_results.json @@ -0,0 +1 @@ +{"TypeInfo":"BEC:github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec/iterative/Simple","ECCInfo":"111543035e239f539647c220bbabfd5a","Stats":{"0.01":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.1":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.2":{"ChannelCodewordError":{"Count":1000000,"Mean":0.030311399999999943,"S":5143.499030040308},"ChannelMessageError":{"Count":1000000,"Mean":0.032171545454545376,"S":6290.669349016466},"ChannelParityError":{"Count":1000000,"Mean":0.025196000000000156,"S":7308.786584000052}},"0.3":{"ChannelCodewordError":{"Count":1000000,"Mean":0.14615093333333073,"S":14726.02024135099},"ChannelMessageError":{"Count":1000000,"Mean":0.1531665454545421,"S":18645.546543603046},"ChannelParityError":{"Count":1000000,"Mean":0.12685799999999428,"S":32125.422836000726}},"0.4":{"ChannelCodewordError":{"Count":1000000,"Mean":0.39107000000001163,"S":515.5884333333495},"ChannelMessageError":{"Count":1000000,"Mean":0.3905850909090831,"S":7053.8652719335},"ChannelParityError":{"Count":1000000,"Mean":0.3924034999999984,"S":45495.86818774953}},"0.5":{"ChannelCodewordError":{"Count":1000000,"Mean":0.4643276666666691,"S":150.4624123333463},"ChannelMessageError":{"Count":1000000,"Mean":0.464196636363629,"S":6724.879482900737},"ChannelParityError":{"Count":1000000,"Mean":0.46468799999999816,"S":48451.06265599885}},"0.6":{"ChannelCodewordError":{"Count":1000000,"Mean":0.6,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.6001200909090976,"S":6235.851693876081},"ChannelParityError":{"Count":1000000,"Mean":0.5996697500000328,"S":47158.62843494023}},"0.7":{"ChannelCodewordError":{"Count":1000000,"Mean":0.6666666666666666,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.6665624545454347,"S":5785.9239424049765},"ChannelParityError":{"Count":1000000,"Mean":0.6669532499999958,"S":43756.04981443932}},"0.8":{"ChannelCodewordError":{"Count":1000000,"Mean":0.8,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.8000206363636302,"S":4159.832631991517},"ChannelParityError":{"Count":1000000,"Mean":0.7999432500000084,"S":31458.73427943751}},"0.9":{"ChannelCodewordError":{"Count":1000000,"Mean":0.8666666666666667,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.8667447272727224,"S":3003.478571371803},"ChannelParityError":{"Count":1000000,"Mean":0.8664520000000066,"S":22713.80669600095}},"0.99":{"ChannelCodewordError":{"Count":1000000,"Mean":0.9333333333333333,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.9334538181818238,"S":1621.2585779834733},"ChannelParityError":{"Count":1000000,"Mean":0.9330019999999811,"S":12260.767995999611}}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_15_11_bsc_dwbf_results.json b/data/linearblock/hamming/hamming_15_11_bsc_dwbf_results.json new file mode 100644 index 0000000..15ffc57 --- /dev/null +++ b/data/linearblock/hamming/hamming_15_11_bsc_dwbf_results.json @@ -0,0 +1 @@ +{"TypeInfo":"BSC:github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision/DWBF_F","ECCInfo":"111543035e239f539647c220bbabfd5a","Stats":{"0.01":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.05":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.1":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.15":{"ChannelCodewordError":{"Count":1000000,"Mean":0.2,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.19990563636364345,"S":4153.571260793176},"ChannelParityError":{"Count":1000000,"Mean":0.20025949999999848,"S":31411.382659749437}},"0.2":{"ChannelCodewordError":{"Count":1000000,"Mean":0.2615021333333232,"S":317.6298176711228},"ChannelMessageError":{"Count":1000000,"Mean":0.2615455454545368,"S":5322.4813718925425},"ChannelParityError":{"Count":1000000,"Mean":0.2613827500000011,"S":38168.87050243619}},"0.25":{"ChannelCodewordError":{"Count":1000000,"Mean":0.26154519999999465,"S":315.2016902933231},"ChannelMessageError":{"Count":1000000,"Mean":0.2615265454545322,"S":5323.66767550416},"ChannelParityError":{"Count":1000000,"Mean":0.2615965000000101,"S":38193.771187749226}},"0.3":{"ChannelCodewordError":{"Count":1000000,"Mean":0.2871956666666505,"S":3680.2357145557507},"ChannelMessageError":{"Count":1000000,"Mean":0.2871546363636297,"S":8913.892500818225},"ChannelParityError":{"Count":1000000,"Mean":0.2873085000000046,"S":43202.07582775158}},"0.35":{"ChannelCodewordError":{"Count":1000000,"Mean":0.3450296666666789,"S":4060.746897666593},"ChannelMessageError":{"Count":1000000,"Mean":0.3451319999999924,"S":9814.960427239894},"ChannelParityError":{"Count":1000000,"Mean":0.3447482499999966,"S":47557.45662193539}},"0.4":{"ChannelCodewordError":{"Count":1000000,"Mean":0.4181674666666718,"S":3864.840932693221},"ChannelMessageError":{"Count":1000000,"Mean":0.4182231818181955,"S":10070.386718801594},"ChannelParityError":{"Count":1000000,"Mean":0.4180142500000069,"S":50885.52429693832}},"0.45":{"ChannelCodewordError":{"Count":1000000,"Mean":0.4181525999999994,"S":3866.5408910176848},"ChannelMessageError":{"Count":1000000,"Mean":0.41812509090909733,"S":10094.631492759847},"ChannelParityError":{"Count":1000000,"Mean":0.4182282499999999,"S":50917.56840193852}},"0.5":{"ChannelCodewordError":{"Count":1000000,"Mean":0.4766141333333307,"S":4045.607906915413},"ChannelMessageError":{"Count":1000000,"Mean":0.47664245454546167,"S":10428.623417446597},"ChannelParityError":{"Count":1000000,"Mean":0.4765362499999966,"S":52240.51493593664}}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_15_11_bsc_gallager_results.json b/data/linearblock/hamming/hamming_15_11_bsc_gallager_results.json new file mode 100644 index 0000000..1621aaa --- /dev/null +++ b/data/linearblock/hamming/hamming_15_11_bsc_gallager_results.json @@ -0,0 +1 @@ +{"TypeInfo":"BSC:github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision/Gallager","ECCInfo":"111543035e239f539647c220bbabfd5a","Stats":{"0.01":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.05":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.1":{"ChannelCodewordError":{"Count":1000000,"Mean":0.0002745333333333369,"S":55.40018700444445},"ChannelMessageError":{"Count":1000000,"Mean":0.00026709090909090807,"S":58.969984760331585},"ChannelParityError":{"Count":1000000,"Mean":0.000294999999999997,"S":113.28797499999867}},"0.15":{"ChannelCodewordError":{"Count":1000000,"Mean":0.20004820000000603,"S":13.468787871112115},"ChannelMessageError":{"Count":1000000,"Mean":0.20006818181817465,"S":4167.433367768621},"ChannelParityError":{"Count":1000000,"Mean":0.19999324999999574,"S":31431.137454437685}},"0.2":{"ChannelCodewordError":{"Count":1000000,"Mean":0.2615398666666597,"S":321.7870328711009},"ChannelMessageError":{"Count":1000000,"Mean":0.26159027272728475,"S":5327.314338438109},"ChannelParityError":{"Count":1000000,"Mean":0.26140125000000036,"S":38198.94899843835}},"0.25":{"ChannelCodewordError":{"Count":1000000,"Mean":0.261529799999991,"S":322.150378626659},"ChannelMessageError":{"Count":1000000,"Mean":0.261611363636356,"S":5338.775408057654},"ChannelParityError":{"Count":1000000,"Mean":0.26130550000000374,"S":38260.68566974736}},"0.3":{"ChannelCodewordError":{"Count":1000000,"Mean":0.287255666666648,"S":3683.6797456665017},"ChannelMessageError":{"Count":1000000,"Mean":0.287185363636373,"S":8906.442946107396},"ChannelParityError":{"Count":1000000,"Mean":0.2874490000000169,"S":43203.07239899912}},"0.35":{"ChannelCodewordError":{"Count":1000000,"Mean":0.3450120666666644,"S":4069.567187728747},"ChannelMessageError":{"Count":1000000,"Mean":0.3451319090909185,"S":9839.312434701873},"ChannelParityError":{"Count":1000000,"Mean":0.3446825000000052,"S":47625.47419375111}},"0.4":{"ChannelCodewordError":{"Count":1000000,"Mean":0.41812846666665515,"S":3872.72314076008},"ChannelMessageError":{"Count":1000000,"Mean":0.41816845454545376,"S":10073.8461024048},"ChannelParityError":{"Count":1000000,"Mean":0.41801850000002616,"S":50925.03365774801}},"0.45":{"ChannelCodewordError":{"Count":1000000,"Mean":0.41828980000001104,"S":3866.972104848858},"ChannelMessageError":{"Count":1000000,"Mean":0.41818836363637,"S":10073.302436495378},"ChannelParityError":{"Count":1000000,"Mean":0.4185687499999955,"S":50879.76402343719}},"0.5":{"ChannelCodewordError":{"Count":1000000,"Mean":0.476727333333339,"S":4046.4007639997876},"ChannelMessageError":{"Count":1000000,"Mean":0.47673299999999486,"S":10425.25828124782},"ChannelParityError":{"Count":1000000,"Mean":0.4767117499999942,"S":52300.969911936685}}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_7_4.json b/data/linearblock/hamming/hamming_7_4.json new file mode 100644 index 0000000..01b9a12 --- /dev/null +++ b/data/linearblock/hamming/hamming_7_4.json @@ -0,0 +1 @@ +{"H":{"Rows":3,"Cols":7,"Data":[[0,2,4,6],[1,2,5,6],[3,4,5,6]]},"Processing":{"HColumnOrder":[3,4,5,2,0,1,6],"G":{"Rows":4,"Cols":7,"Data":[[0,4,5,6],[1,5,6],[2,4,6],[3,4,5]]}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_7_4_bec_results.json b/data/linearblock/hamming/hamming_7_4_bec_results.json new file mode 100644 index 0000000..6601479 --- /dev/null +++ b/data/linearblock/hamming/hamming_7_4_bec_results.json @@ -0,0 +1 @@ +{"TypeInfo":"BEC:github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec/iterative/Simple","ECCInfo":"b538a6813c0fd136b9f13efc194c03b5","Stats":{"0.01":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.1":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.2":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.3":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.4":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.5":{"ChannelCodewordError":{"Count":1000000,"Mean":0.12243300000000157,"S":37481.446225286294},"ChannelMessageError":{"Count":1000000,"Mean":0.12866099999999972,"S":47823.097078999825},"ChannelParityError":{"Count":1000000,"Mean":0.11412899999999784,"S":44003.12691455631}},"0.6":{"ChannelCodewordError":{"Count":1000000,"Mean":0.5225524285713955,"S":4593.428781857105},"ChannelMessageError":{"Count":1000000,"Mean":0.5220064999999997,"S":33529.46395775057},"ChannelParityError":{"Count":1000000,"Mean":0.5232803333333182,"S":52664.02607988657}},"0.7":{"ChannelCodewordError":{"Count":1000000,"Mean":0.5224564285714399,"S":4597.7496321427525},"ChannelMessageError":{"Count":1000000,"Mean":0.5214964999999961,"S":33468.900487749845},"ChannelParityError":{"Count":1000000,"Mean":0.5237363333333185,"S":52593.69759100301}},"0.8":{"ChannelCodewordError":{"Count":1000000,"Mean":0.7142857142857143,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.7144205000000197,"S":25483.474179748602},"ChannelParityError":{"Count":1000000,"Mean":0.7141059999999881,"S":45303.95409733253}},"0.9":{"ChannelCodewordError":{"Count":1000000,"Mean":0.8571428571428571,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.8568172499999847,"S":15294.38760243743},"ChannelParityError":{"Count":1000000,"Mean":0.857576999999981,"S":27190.022404334664}},"0.99":{"ChannelCodewordError":{"Count":1000000,"Mean":0.8571428571428571,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.8570747499999994,"S":15303.685412437639},"ChannelParityError":{"Count":1000000,"Mean":0.8572336666666631,"S":27206.55184433464}}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_7_4_bsc_dwbf_results.json b/data/linearblock/hamming/hamming_7_4_bsc_dwbf_results.json new file mode 100644 index 0000000..6d81d7e --- /dev/null +++ b/data/linearblock/hamming/hamming_7_4_bsc_dwbf_results.json @@ -0,0 +1 @@ +{"TypeInfo":"BSC:github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision/DWBF_F","ECCInfo":"b538a6813c0fd136b9f13efc194c03b5","Stats":{"0.01":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.05":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.1":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.15":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.2":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.25":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.3":{"ChannelCodewordError":{"Count":1000000,"Mean":0.42857142857142855,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.42859624999999263,"S":30645.691985936486},"ChannelParityError":{"Count":1000000,"Mean":0.4285383333333383,"S":54481.230197226374}},"0.35":{"ChannelCodewordError":{"Count":1000000,"Mean":0.42857142857142855,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.4287372499999994,"S":30645.432962437317},"ChannelParityError":{"Count":1000000,"Mean":0.4283503333333387,"S":54480.76971100154}},"0.4":{"ChannelCodewordError":{"Count":1000000,"Mean":0.42857142857142855,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0.42864675000001196,"S":30621.40121443599},"ChannelParityError":{"Count":1000000,"Mean":0.4284710000000026,"S":54438.0466034422}},"0.45":{"ChannelCodewordError":{"Count":1000000,"Mean":0.5428602857142862,"S":3265.036724816346},"ChannelMessageError":{"Count":1000000,"Mean":0.5428599999999956,"S":33877.645399998575},"ChannelParityError":{"Count":1000000,"Mean":0.542860666666677,"S":57669.62991955738}},"0.5":{"ChannelCodewordError":{"Count":1000000,"Mean":0.5428824285714199,"S":3263.138136142847},"ChannelMessageError":{"Count":1000000,"Mean":0.54290050000002,"S":33840.047099750685},"ChannelParityError":{"Count":1000000,"Mean":0.5428583333333209,"S":57677.38548611025}}}} \ No newline at end of file diff --git a/data/linearblock/hamming/hamming_7_4_bsc_gallager_results.json b/data/linearblock/hamming/hamming_7_4_bsc_gallager_results.json new file mode 100644 index 0000000..44b6968 --- /dev/null +++ b/data/linearblock/hamming/hamming_7_4_bsc_gallager_results.json @@ -0,0 +1 @@ +{"TypeInfo":"BSC:github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bitflipping/harddecision/Gallager","ECCInfo":"b538a6813c0fd136b9f13efc194c03b5","Stats":{"0.01":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.05":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.1":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.15":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.2":{"ChannelCodewordError":{"Count":1000000,"Mean":0,"S":0},"ChannelMessageError":{"Count":1000000,"Mean":0,"S":0},"ChannelParityError":{"Count":1000000,"Mean":0,"S":0}},"0.25":{"ChannelCodewordError":{"Count":1000000,"Mean":0.00021671428571428368,"S":93.4836471632645},"ChannelMessageError":{"Count":1000000,"Mean":0.00023000000000000196,"S":121.57210000000079},"ChannelParityError":{"Count":1000000,"Mean":0.0001989999999999872,"S":107.62706566666941}},"0.3":{"ChannelCodewordError":{"Count":1000000,"Mean":0.42854328571428923,"S":38.10124879591816},"ChannelMessageError":{"Count":1000000,"Mean":0.42863249999999287,"S":30673.679943750434},"ChannelParityError":{"Count":1000000,"Mean":0.42842433333332375,"S":54501.590607886246}},"0.35":{"ChannelCodewordError":{"Count":1000000,"Mean":0.42854714285714346,"S":37.42798163265286},"ChannelMessageError":{"Count":1000000,"Mean":0.42868900000000626,"S":30612.86627899866},"ChannelParityError":{"Count":1000000,"Mean":0.428358000000001,"S":54397.20161377632}},"0.4":{"ChannelCodewordError":{"Count":1000000,"Mean":0.4285580000000029,"S":34.57124824489785},"ChannelMessageError":{"Count":1000000,"Mean":0.4286065000000088,"S":30637.84315774983},"ChannelParityError":{"Count":1000000,"Mean":0.4284933333333348,"S":54439.90773333281}},"0.45":{"ChannelCodewordError":{"Count":1000000,"Mean":0.5427977142857061,"S":3270.3964682449014},"ChannelMessageError":{"Count":1000000,"Mean":0.5426662499999572,"S":33935.278610936926},"ChannelParityError":{"Count":1000000,"Mean":0.5429730000000028,"S":57771.54349321943}},"0.5":{"ChannelCodewordError":{"Count":1000000,"Mean":0.5428719999999888,"S":3264.2773302856212},"ChannelMessageError":{"Count":1000000,"Mean":0.5430697499999954,"S":33908.55913493762},"ChannelParityError":{"Count":1000000,"Mean":0.5426083333333102,"S":57676.307708334294}}}} \ No newline at end of file diff --git a/go.mod b/go.mod new file mode 100644 index 0000000..053223c --- /dev/null +++ b/go.mod @@ -0,0 +1,23 @@ +module github.com/nathanhack/errorcorrectingcodes + +go 1.14 + +require ( + github.com/cheggaaa/pb/v3 v3.0.5 + github.com/cloudfoundry/gosigar v1.1.0 + github.com/google/go-cmp v0.5.1 // indirect + github.com/mattn/go-runewidth v0.0.12 // indirect + github.com/nathanhack/avgstd v0.0.0-20210531154130-989d51af96a2 + github.com/nathanhack/intmat v0.0.0-20210228003141-935751eb5e16 + github.com/nathanhack/sparsemat v0.0.0-20210707154153-866acabf5d96 + github.com/nathanhack/threadpool v0.0.0-20210629001431-9d0d768b0f6d + github.com/onsi/ginkgo v1.14.1 // indirect + github.com/onsi/gomega v1.10.2 // indirect + github.com/pkg/errors v0.9.1 // indirect + github.com/rivo/uniseg v0.2.0 // indirect + github.com/sirupsen/logrus v1.7.0 + github.com/spf13/cobra v1.1.3 + github.com/stretchr/testify v1.6.1 // indirect + gonum.org/v1/gonum v0.8.2 + +) diff --git a/linearblock/cycle.go b/linearblock/cycle.go new file mode 100644 index 0000000..0f880cd --- /dev/null +++ b/linearblock/cycle.go @@ -0,0 +1,243 @@ +package linearblock + +import ( + "fmt" + mat "github.com/nathanhack/sparsemat" + "math" + "runtime" + "strings" + "sync" +) + +//Node is a structure containing the Index and state of whether it is a check node or a variable node +type Node struct { + Index int + Check bool +} + +//Cycle is a slice of Nodes +type Cycle []Node + +//String returns a standard rep of the cycle +func (c Cycle) String() string { + sb := strings.Builder{} + sb.WriteString("[") + for i, n := range c { + if n.Check { + sb.WriteString(fmt.Sprintf("c:%v", n.Index)) + } else { + sb.WriteString(fmt.Sprintf("v:%v", n.Index)) + } + if i < len(c)-1 { + sb.WriteString(" ") + } + } + sb.WriteString("]") + return sb.String() +} + +//Equal compares two cycles to see if they are equal. To be equal +// the first node must be equal but they could be in opposite orderVector. +func (c Cycle) Equal(c2 Cycle) bool { + if len(c) != len(c2) { + return false + } + + if len(c) == 0 { + return true + } + + if c[0] != c2[0] { + return false + } + + t1 := c[1:] + t2 := c2[1:] + + //we try forward direction first + failed := false + for i, n := range t1 { + if t2[i] != n { + failed = true + break + } + } + + if !failed { + return true + } + + //else we check the opposite direction + l := len(t1) + for i, n := range t1 { + if t2[l-1-i] != n { + return false + } + } + return true +} + +//SmallestCycle non-deterministically returns a cycle from the set of smallest cycles. +func SmallestCycle(m mat.SparseMat, parallel bool) Cycle { + rows, _ := m.Dims() + wg := sync.WaitGroup{} + threadSize := 1 + if parallel { + threadSize = runtime.NumCPU() + } + limit := make(chan bool, threadSize) + mux := sync.RWMutex{} + var smallest []Node + girth := -1 + for i := 0; i < rows; i++ { + wg.Add(1) + limit <- true + go func(index int) { + c := smallestCycle(m, index, girth-2) + + mux.Lock() + if len(smallest) == 0 || (len(c) != 0 && len(smallest) > len(c)) { + smallest = c + girth = len(smallest) + } + mux.Unlock() + wg.Done() + <-limit + }(i) + } + wg.Wait() + return smallest +} + +func smallestCycle(m mat.SparseMat, checkIndex int, minGirth int) []Node { + if minGirth < 0 { + minGirth = math.MaxInt64 + } + //we make a history that will alternate between variable nodes and check nodes + // as we extend to each new hop away from the checkIndex + history := make([]map[int]girthNode, 0) + rows, _ := m.Dims() + + //we prime the history + hop := make(map[int]girthNode) + for i := range m.Row(checkIndex).NonzeroMap() { + hop[i] = girthNode{parentIndex: checkIndex} + } + //if there was only one variable node then there is no way + // this will have a loop + if len(hop) == 1 { + return nil + } + history = append(history, hop) + + for level := 1; level < 2*rows && level < minGirth/2+1; level++ { + if level%2 == 0 { + //this round we look at adding check nodes + prevHop := history[level-1] + hop := make(map[int]girthNode) + for v, gn := range prevHop { + for i := range m.Row(v).NonzeroMap() { + if i == gn.parentIndex { + continue + } + _, has := hop[i] + if has { + //so we make two node lists both start from this ith index + // but one will be what's already in history/hop + // and the other will be from this new girthNode + // at the end we'll concatenate the lists (while removing dup the starting node) + + //first we get the lists (it contains all nodes except the current one) + a := path(history, hop[i].parentIndex, true) + b := path(history, v, true) + + //we add the connection to the checknode in questions + a = append([]Node{{ + Index: checkIndex, + Check: true, + }}, a...) + + // we add the node to one of them + a = append(a, Node{Index: i, Check: false}) + + //next we reverse the orderVector of b + reverse(b) + + //now concatenate the two together + return append(a, b...) + } + hop[i] = girthNode{parentIndex: v} + } + } + history = append(history, hop) + } else { + //this round we look at adding variable nodes + prevHop := history[level-1] + hop := make(map[int]girthNode) + for v, gn := range prevHop { + for i := range m.Column(v).NonzeroMap() { + if i == gn.parentIndex { + continue + } + _, has := hop[i] + if has || i == checkIndex { + //so we make two node lists both start from this ith index + // but one will be what's already in history/hop + // and the other will be from this new girthNode + // at the end we'll concatenate the lists (while removing dup the starting node) + + //first we get the lists (it contains all nodes except the current one) + a := path(history, hop[i].parentIndex, false) + b := path(history, v, false) + + //we add the connection to the checknode in questions + a = append([]Node{{ + Index: checkIndex, + Check: true, + }}, a...) + + // we add the node to one of them + if i != checkIndex { + a = append(a, Node{Index: i, Check: true}) + } + + //next we reverse the orderVector of b + reverse(b) + + //now concatenate the two together + return append(a, b...) + } + hop[i] = girthNode{parentIndex: v} + } + } + history = append(history, hop) + } + } + return nil +} + +func path(history []map[int]girthNode, index int, check bool) []Node { + //we start from the last item in the history + histLen := len(history) + + path := make([]Node, 0, histLen+1) + + for i := histLen - 1; i >= 0; i-- { + n := Node{ + Index: index, + Check: check, + } + + path = append(path, n) + index = history[i][index].parentIndex + check = !check + } + reverse(path) // now we correct the orderVector + return path +} + +func reverse(nodes []Node) { + for i, j := 0, len(nodes)-1; i < j; i, j = i+1, j-1 { + nodes[i], nodes[j] = nodes[j], nodes[i] + } +} diff --git a/linearblock/cycle_test.go b/linearblock/cycle_test.go new file mode 100644 index 0000000..459195e --- /dev/null +++ b/linearblock/cycle_test.go @@ -0,0 +1,67 @@ +package linearblock + +import ( + mat "github.com/nathanhack/sparsemat" + "strconv" + "testing" +) + +func TestSmallestCycle(t *testing.T) { + tests := []struct { + h mat.SparseMat + expected Cycle + }{ + {mat.CSRMat(2, 2, 1, 1, 1, 1), Cycle{ + {Index: 0, Check: true}, + {Index: 0, Check: false}, + {Index: 1, Check: true}, + {Index: 1, Check: false}, + }}, + {mat.CSRMat(3, 3, 1, 1, 0, 0, 1, 1, 1, 0, 1), Cycle{ + {Index: 0, Check: true}, + {Index: 0, Check: false}, + {Index: 2, Check: true}, + {Index: 2, Check: false}, + {Index: 1, Check: true}, + {Index: 1, Check: false}, + }}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + actual := SmallestCycle(test.h, false) + if !actual.Equal(test.expected) { + t.Fatalf("expected %v but found %v", test.expected, actual) + } + }) + } +} + +func TestCycle_Equal(t *testing.T) { + + c1 := Cycle{ + {Index: 0, Check: true}, + {Index: 0, Check: false}, + {Index: 1, Check: true}, + {Index: 1, Check: false}, + } + c2 := Cycle{ + {Index: 0, Check: true}, + {Index: 0, Check: false}, + {Index: 1, Check: true}, + {Index: 1, Check: false}, + } + c3 := Cycle{ + {Index: 0, Check: true}, + {Index: 1, Check: false}, + {Index: 1, Check: true}, + {Index: 0, Check: false}, + } + + if !c1.Equal(c2) { + t.Fatalf("expected equal") + } + + if !c1.Equal(c3) { + t.Fatalf("expected equal") + } +} diff --git a/linearblock/girth.go b/linearblock/girth.go new file mode 100644 index 0000000..45ce408 --- /dev/null +++ b/linearblock/girth.go @@ -0,0 +1,160 @@ +package linearblock + +import ( + "context" + mat "github.com/nathanhack/sparsemat" + "github.com/nathanhack/threadpool" + "math" + "sync" +) + +//CalculateGirthLowerBoundByEdges returns a bool if true it is possible for the matrix to be +// free of all cycles C_i where i \in 3<=i<=minGirth. If false, it is impossible to +// be free of all cycles C_i. +func CalculateGirthLowerBoundByEdges(m mat.SparseMat, minGirth int) bool { + // from the paper Fast Distributed Algorithms for Girth, Cycles and Small Subgraphs by K. Censor-Hillel, et al + // that if m is free of all C_i cycles 3<=i<=2k then m contains at most n^(1+1/k)+n edges. + + rows, cols := m.Dims() + edges := 0 + for i := 0; i < rows; i++ { + edges += len(m.Row(i).NonzeroMap()) + } + + n := float64(rows + cols) + n = math.Pow(n, 1+1/(float64(minGirth)/2)) + n + if int(n) < edges { + return false + } + return true +} + +type girthNode struct { + parentIndex int +} + +//CalculateGirth calculates the girth of the tanner csv induced by m. +// threads specifies the number of threads to use if <=0 will use runtime.NumCPU() +func CalculateGirth(m mat.SparseMat, threads int) int { + return CalculateGirthLowerBound(m, -1, threads) +} + +//CalculateGirthLowerBound returns the length of the smallest cycle. +// It searches for cycles with a length <= smallestGirth. If no cycles are found +// that are smaller or equal to smallestGirth then it returns -1. +// threads specifies the number of threads to use if <=0 will use runtime.NumCPU() +func CalculateGirthLowerBound(m mat.SparseMat, smallestGirth, threads int) int { + if smallestGirth != -1 && (smallestGirth < 4 || smallestGirth%2 != 0) { + panic("smallestGirth == -1 or smallestGirth must be a even number >=4") + } + + rows, _ := m.Dims() + + pool := threadpool.New(context.Background(), threads, rows) + calulated := -1 + mux := sync.RWMutex{} + for i := 0; i < rows; i++ { + index := i // + pool.Add(func() { + + mux.RLock() + g := CalculateCycleLowerBound(m, index, smallestGirth) + mux.RUnlock() + + mux.Lock() + if g > 0 && (g <= smallestGirth || smallestGirth == -1) { + smallestGirth = g + calulated = g + } + mux.Unlock() + }) + } + pool.Wait() + return calulated +} + +//HasGirthSmallerThan will search for cycle smaller than the given cycleLen. +// Return true if it found a cycle smaller than cycleLen, else returns false. +// threads specifies the number of threads to use if <=0 will use runtime.NumCPU() +func HasGirthSmallerThan(m mat.SparseMat, cycleLen, threads int) bool { + if cycleLen != -1 && cycleLen < 4 { + panic("cycleLen == -1 or cycleLen >=4 required") + } + + rows, _ := m.Dims() + pool := threadpool.New(context.Background(), threads, rows) + smaller := false + mux := sync.RWMutex{} + for i := 0; i < rows; i++ { + index := i + pool.Add(func() { + mux.RLock() + if smaller { + // nothing to do here + defer mux.RUnlock() + return + } + g := CalculateCycleLowerBound(m, index, cycleLen) + + mux.RUnlock() + mux.Lock() + if g > 0 && g < cycleLen { + smaller = true + } + mux.Unlock() + }) + } + pool.Wait() + return smaller +} + +//CalculateCycleLowerBound runs a BFS starting at the checkIndex check node, for maxGirth/2 steps +// if maxGirth ==-1 it will search until it finds a cycle +// in either case it returns the length of the cycle (up to maxGirth) or -1 if no cycle was found +func CalculateCycleLowerBound(m mat.SparseMat, checkIndex, maxGirth int) int { + if maxGirth == -1 { + maxGirth = math.MaxInt64 + } + //we make a history that will alternate between variable nodes and check nodes + // as we extend to each new hop away from the checkIndex + history := make([]map[int]girthNode, 0) + rows, _ := m.Dims() + + //we prime the history + hop := make(map[int]girthNode) + for _, i := range m.Row(checkIndex).NonzeroArray() { + hop[i] = girthNode{parentIndex: checkIndex} + } + //if there was only one variable node then there is no way + // this will have a loop + if len(hop) == 1 { + return -1 + } + history = append(history, hop) + + for level := 1; level < 2*rows && level < maxGirth/2+1; level++ { + prevHop := history[level-1] + hop := make(map[int]girthNode) + for v, gn := range prevHop { + levelHop := level % 2 + var indices []int + if levelHop == 0 { + indices = m.Row(v).NonzeroArray() + } else { + indices = m.Column(v).NonzeroArray() + } + for _, i := range indices { + if i == gn.parentIndex { + continue + } + _, has := hop[i] + if has || (levelHop == 1 && i == checkIndex) { + return (level + 1) * 2 + } + hop[i] = girthNode{parentIndex: v} + } + } + history = append(history, hop) + } + return -1 +} diff --git a/linearblock/girth_test.go b/linearblock/girth_test.go new file mode 100644 index 0000000..4195290 --- /dev/null +++ b/linearblock/girth_test.go @@ -0,0 +1,85 @@ +package linearblock + +import ( + mat "github.com/nathanhack/sparsemat" + "strconv" + "testing" +) + +func TestCalculateGirthLowerBound(t *testing.T) { + tests := []struct { + h mat.SparseMat + minGirth int + expected int + }{ + {mat.CSRIdentity(500), -1, -1}, + {mat.CSRIdentity(500), 4, -1}, + {mat.CSRMat(2, 2, 1, 1, 1, 1), -1, 4}, + {mat.CSRMat(2, 2, 1, 1, 1, 1), 6, 4}, + {mat.CSRMat(2, 2, 1, 0, 0, 1), -1, -1}, + {mat.CSRMat(2, 2, 1, 0, 0, 1), 4, -1}, + {mat.CSRMat(4, 8, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1), -1, 8}, + {mat.CSRMat(3, 6, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1), -1, 6}, + {mat.CSRMat(3, 6, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1), 6, 6}, + {mat.CSRMat(3, 6, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1), 4, -1}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + actual := CalculateGirthLowerBound(test.h, test.minGirth, -1) + if actual != test.expected { + t.Fatalf("expected %v but found %v", test.expected, actual) + } + }) + } +} +func TestCalculateCycleLowerBound(t *testing.T) { + tests := []struct { + h mat.SparseMat + checkIndex int + minGirth int + expected int + }{ + {mat.CSRMat(2, 2, 1, 1, 1, 1), 0, -1, 4}, + {mat.CSRMat(2, 2, 1, 1, 1, 1), 0, 6, 4}, + {mat.CSRMat(2, 2, 1, 0, 0, 1), 0, -1, -1}, + {mat.CSRMat(2, 2, 1, 0, 0, 1), 0, 4, -1}, + {mat.CSRMat(2, 2, 1, 0, 1, 0), 0, -1, -1}, + {mat.CSRMat(2, 2, 1, 0, 1, 0), 0, 4, -1}, + {mat.CSRIdentity(500), 0, -1, -1}, + {mat.CSRIdentity(500), 0, 4, -1}, + {mat.CSRMat(4, 8, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1), 0, -1, 8}, + {mat.CSRMat(3, 6, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1), 0, -1, 6}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + actual := CalculateCycleLowerBound(test.h, test.checkIndex, test.minGirth) + if actual != test.expected { + t.Fatalf("expected %v but found %v", test.expected, actual) + } + }) + } +} +func BenchmarkCalculateGirthLowerBound(b *testing.B) { + h := mat.CSRMat(2, 2, 1, 1, 1, 1) + for i := 0; i < b.N; i++ { + CalculateGirthLowerBound(h, -1, 1) + } +} + +func BenchmarkCalculateGirthLowerBound2(b *testing.B) { + h := mat.CSRIdentity(10000) + b.ResetTimer() + for i := 0; i < b.N; i++ { + CalculateGirthLowerBound(h, -1, 0) + } +} +func BenchmarkCalculateGirthLowerBound3(b *testing.B) { + h := mat.CSRIdentity(10000) + for i := 0; i < 10000-2; i++ { + h.SetMatrix(mat.CSRMat(2, 2, 1, 1, 1, 1), i, i) + } + b.ResetTimer() + for i := 0; i < b.N; i++ { + CalculateGirthLowerBound(h, -1, 0) + } +} diff --git a/linearblock/hamming/hamming.go b/linearblock/hamming/hamming.go new file mode 100644 index 0000000..4f1c2af --- /dev/null +++ b/linearblock/hamming/hamming.go @@ -0,0 +1,46 @@ +package hamming + +import ( + "context" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + "github.com/nathanhack/errorcorrectingcodes/linearblock/internal" + mat "github.com/nathanhack/sparsemat" +) + +//New creates the systematic hamming code with paritySymbols number of parity symbols. +// Hamming codes can detect up to two-bit errors or correct one-bit errors without +// detection of uncorrected errors. +func New(ctx context.Context, paritySymbols int, threads int) (*linearblock.LinearBlock, error) { + if paritySymbols < 3 { + panic("hamming codes require >=3 parity symbols") + } + n := 1< [1,n] (note they're nonzero) + for i := 1; i <= n; i++ { + vec := mat.CSRVec(paritySymbols) + for j := 0; j < paritySymbols; j++ { + if i&(1< 0 { + vec.Set(j, 1) + } + } + H.SetColumn(i-1, vec) + } + + order, g := internal.NewFromH(ctx, H, threads) + if order == nil { + return nil, fmt.Errorf("unable to create generator for H matrix") + } + + return &linearblock.LinearBlock{ + H: H, + Processing: &linearblock.Systemic{ + HColumnOrder: order, + G: g, + }, + }, nil +} diff --git a/linearblock/hamming/hamming_test.go b/linearblock/hamming/hamming_test.go new file mode 100644 index 0000000..e26276e --- /dev/null +++ b/linearblock/hamming/hamming_test.go @@ -0,0 +1,27 @@ +package hamming + +import ( + "context" + "strconv" + "testing" +) + +func TestNew(t *testing.T) { + tests := []struct { + paritySymbols int + }{ + {3}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + actual, err := New(context.Background(), test.paritySymbols, 0) + if err != nil { + t.Fatalf("expected no error found :%v", err) + } + + if !actual.Validate() { + t.Fatalf("expected valid linearblock code") + } + }) + } +} diff --git a/linearblock/internal/gaussianjordan.go b/linearblock/internal/gaussianjordan.go new file mode 100644 index 0000000..4c22f7b --- /dev/null +++ b/linearblock/internal/gaussianjordan.go @@ -0,0 +1,226 @@ +package internal + +import ( + "context" + "github.com/cheggaaa/pb/v3" + mat "github.com/nathanhack/sparsemat" + "github.com/nathanhack/threadpool" + "github.com/sirupsen/logrus" + "os" + "sync" +) + +func swapColOrder(i, j int, colIndices []int) { + x := len(colIndices) + if 0 <= i && i < x && 0 <= j && j < x { + idx := colIndices[i] + colIndices[i] = colIndices[j] + colIndices[j] = idx + } +} + +func findPivotColGF2(H mat.SparseMat, forRow int) int { + rows, _ := H.Dims() + + for r := forRow; r < rows; r++ { + row := H.Row(r).NonzeroArray() + if len(row) == 0 { + continue + } + + col := row[len(row)-1] + if col > forRow { + return col + } + } + return -1 +} + +func GaussianJordanEliminationGF2(ctx context.Context, H mat.SparseMat, threads int) (mat.SparseMat, []int) { + rows, cols := H.Dims() + result := mat.CSRMatCopy(H) + columnSwapHistory := make([]int, cols) + + //initialize the columnIndices + for c := 0; c < cols; c++ { + columnSwapHistory[c] = c + } + + if cols < rows { + //null space must equal the rank + return nil, nil + } + + //to fail fast we first do the lower triangle then do the upper + if lowerTrianglar(ctx, rows, result, columnSwapHistory, threads, logrus.GetLevel() == logrus.DebugLevel) != rows { + logrus.Debugf("All rows not linearly independant") + return nil, nil + } + + // at this point we know we had all linearly independent rows! + // and the lower triangle is done so we'll take care of the top + if !upperTrianglar(ctx, rows, result, threads, logrus.GetLevel() == logrus.DebugLevel) { + return nil, nil + } + + logrus.Debugf("Gaussian-Jordan Elimination complete") + return result, columnSwapHistory +} + +func upperTrianglar(ctx context.Context, rows int, H mat.SparseMat, threads int, showProgressBar bool) bool { + bar := pb.Full.New(rows) + logrus.Debugf("Reduced row echelon") + bar.Set("prefix", "Processing Row ") + bar.SetWriter(os.Stdout) + if showProgressBar { + bar.Start() + } + + for r := 0; r < rows; r++ { + bar.Increment() + select { + case <-ctx.Done(): + return false + default: + } + eliminateOtherRows(ctx, r, H, threads) + } + bar.SetTemplateString(`{{string . "prefix"}}{{counters . }}{{string . "suffix"}}`) + bar.Set("suffix", " Done") + bar.Finish() + return true +} + +func pivotsSwapReturn(H mat.SparseMat, rowIndex int, columnsSwapHistory []int) []int { + pivots := H.Column(rowIndex).NonzeroArray() + if len(pivots) == 0 || pivots[len(pivots)-1] < rowIndex { + // there aren't any where we need them + // so we'll do a columns swap to get what we need + colPivot := findPivotColGF2(H, rowIndex) + if colPivot == -1 { + //we get here when there aren't any more non zero rows + //so this matrix null space doesn't span the rank + return nil + } + + H.SwapColumns(rowIndex, colPivot) //this functions for CSR has a problem + swapColOrder(rowIndex, colPivot, columnsSwapHistory) + + //now redo pivots + pivots = H.Column(rowIndex).NonzeroArray() + } + return pivots +} + +func eliminateOtherRows(ctx context.Context, rowIndex int, result mat.SparseMat, threads int) { + //create a pool with 1 less pivot + pivots := result.Column(rowIndex).NonzeroArray() + pool := threadpool.New(ctx, threads, len(pivots)-1) + rrow := result.Row(rowIndex) + mut := sync.RWMutex{} + + //for all pivots except the one equal to r subtract it (in GF2 subtract is add) + for _, index := range pivots { + pIndex := index + if index != rowIndex { + pool.Add(func() { + func(p int) { + mut.RLock() + prow := result.Row(p) + mut.RUnlock() + prow.Add(prow, rrow) + mut.Lock() + result.SetRow(p, prow) + mut.Unlock() + }(pIndex) + }) + } + } + pool.Wait() +} + +func eliminateLowerRows(ctx context.Context, rowIndex int, result mat.SparseMat, threads int) { + //create a pool with 1 less pivot + pivots := result.Column(rowIndex).NonzeroArray() + pool := threadpool.New(ctx, threads, len(pivots)) + rrow := result.Row(rowIndex) + mut := sync.RWMutex{} + + //for all pivots except the one equal to r subtract it (in GF2 subtract is add) + for _, index := range pivots { + pIndex := index + pool.Add(func() { + if pIndex > rowIndex { + func(p int) { + mut.RLock() + prow := result.Row(p) + mut.RUnlock() + prow.Add(prow, rrow) + mut.Lock() + result.SetRow(p, prow) + mut.Unlock() + }(pIndex) + } + }) + } + pool.Wait() +} + +func CalculateRank(ctx context.Context, H mat.SparseMat, threads int, showProgressBar bool) int { + if H == nil { + return -1 + } + + tmp := mat.CSRMatCopy(H) + + rows, cols := H.Dims() + + min := rows + if cols < rows { + min = cols + } + columnSwapHistory := make([]int, cols) + + return lowerTrianglar(ctx, min, tmp, columnSwapHistory, threads, showProgressBar) +} + +func lowerTrianglar(ctx context.Context, rows int, H mat.SparseMat, columnSwapHistory []int, threads int, showProgressBar bool) int { + bar := pb.Full.New(rows) + logrus.Debugf("Row echelon") + bar.Set("prefix", "Processing Row ") + bar.SetWriter(os.Stdout) + if showProgressBar { + bar.Start() + } + + for r := 0; r < rows; r++ { + select { + case <-ctx.Done(): + return -1 + default: + } + bar.Increment() + //we process both upper and lower triangle at the same time + //first find pivots for the column equal to r + pivots := pivotsSwapReturn(H, r, columnSwapHistory) + if pivots == nil { + return r + } + + //when here we know pivots have values and the last one + //should be a row to switch with + pivot := pivots[len(pivots)-1] + H.SwapRows(r, pivot) + + // now the rth row has a pivot in the rth row and rth column + // so we now subtract it from all other rows with a 1 + // in the rth column + eliminateLowerRows(ctx, r, H, threads) + } + + bar.SetTemplateString(`{{string . "prefix"}}{{counters . }}{{string . "suffix"}}`) + bar.Set("suffix", " Done") + bar.Finish() + + return rows +} diff --git a/linearblock/internal/gaussianjordan_test.go b/linearblock/internal/gaussianjordan_test.go new file mode 100644 index 0000000..646090f --- /dev/null +++ b/linearblock/internal/gaussianjordan_test.go @@ -0,0 +1,40 @@ +package internal + +import ( + "context" + mat "github.com/nathanhack/sparsemat" + "strconv" + "testing" +) + +func TestGaussianJordanEliminationGF2(t *testing.T) { + tests := []struct { + input mat.SparseMat + expected mat.SparseMat + }{ + { //Hamming 7 + mat.CSRMat(3, 7, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1), + mat.CSRMat(3, 7, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1), + }, + { //Random - one linearly dependent row + mat.CSRMat(4, 5, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1), + nil, + }, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + + gen, _ := GaussianJordanEliminationGF2(context.Background(), test.input, 0) + + if test.expected != nil { + if !test.expected.Equals(gen) { + t.Fatalf("expected \n%v\n but found \n%v\n", test.expected, gen) + } + } else { + if gen != nil { + t.Fatalf("expecte nil but found \n%v\n", gen) + } + } + }) + } +} diff --git a/linearblock/internal/utils.go b/linearblock/internal/utils.go new file mode 100644 index 0000000..cd2fe98 --- /dev/null +++ b/linearblock/internal/utils.go @@ -0,0 +1,107 @@ +package internal + +import ( + "context" + "fmt" + mat "github.com/nathanhack/sparsemat" + "github.com/sirupsen/logrus" +) + +func ExtractAFromH(ctx context.Context, H mat.SparseMat, threads int) (A mat.SparseMat, columnOrdering []int) { + m, N := H.Dims() + + gje, ordering := GaussianJordanEliminationGF2(ctx, H, threads) + + if gje == nil { + return nil, nil + } + + //let's check if we got a [ I, * ] format + actual := gje.Slice(0, 0, m, m) + ident := mat.CSRIdentity(m) + if !actual.Equals(ident) { + logrus.Errorf("failed to create transform H matrix into [I,*]") + for i := 0; i < m; i++ { + x := actual.Row(i).NonzeroArray() + if len(x) > 0 { + fmt.Println(i, ":", x) + } + } + return nil, nil + } + + //we need to convert gje from [ I, A] to [ A, I] (while keeping track) + // and then extract A + + // first the keeping track part + columnOrdering = make([]int, len(ordering)) + copy(columnOrdering[0:N-m], ordering[m:N]) + copy(columnOrdering[N-m:N], ordering[0:m]) + + //finally extract the A + A = gje.Slice(0, m, m, N-m) + return +} + +//NewFromH creates a derived H an G matrix containing the systematic form of parity matrix H and generator G. +// Note: the LinearBlock parity matrix's columns may be swapped. +func NewFromH(ctx context.Context, H mat.SparseMat, threads int) (HColumnOrder []int, G mat.SparseMat) { + hrows, hcols := H.Dims() + if hrows >= hcols { + panic("H matrix shape == (rows, cols) where rows < cols required") + } + // So we now take the current H matrix + // convert H=[*] -> H=[A,I] + // then extract out the A and keep track of columnSwaps during it + logrus.Debugf("Creating generator matrix from H matrix") + A, columnSwaps := ExtractAFromH(ctx, H, threads) + if A == nil { + logrus.Debugf("Unable to create generator matrix from H") + return nil, nil + } + + AT := A.T() // transpose of A + atRows, atCols := AT.Dims() + + //Next using A make G=[I, A^T] where A^T is the transpose of A + G = mat.DOKMat(atRows, atRows+atCols) + G.SetMatrix(mat.CSRIdentity(atRows), 0, 0) + G.SetMatrix(AT, 0, atRows) + + logrus.Debugf("Generator Matrix complete") + return columnSwaps, G +} + +func ColumnSwapped(H mat.SparseMat, order []int) mat.SparseMat { + rows, cols := H.Dims() + result := mat.CSRMat(rows, cols) + + for c, c1 := range order { + result.SetColumn(c, H.Column(c1)) + } + return result +} + +//ValidateHGMatrices tests if G*H.T ==0 where H.T is the transpose of H +func ValidateHGMatrices(G, H mat.SparseMat) bool { + rows, _ := G.Dims() + cols, _ := H.Dims() + + //we cache the H.T hopefully this is in CSR so this should be way + // faster than taking the actual H.T() then doing this + cache := make([]mat.SparseVector, cols) + for i := 0; i < cols; i++ { + cache[i] = H.Row(i) + } + for i := 0; i < rows; i++ { + row := G.Row(i) + for j := 0; j < cols; j++ { + //equiv to G*H.T + if row.Dot(cache[j]) > 0 { + return false + } + } + } + + return true +} diff --git a/linearblock/internal/utils_test.go b/linearblock/internal/utils_test.go new file mode 100644 index 0000000..a5c74d9 --- /dev/null +++ b/linearblock/internal/utils_test.go @@ -0,0 +1,40 @@ +package internal + +import ( + "context" + mat "github.com/nathanhack/sparsemat" + "strconv" + "testing" +) + +func TestExtractAFromH(t *testing.T) { + tests := []struct { + input mat.SparseMat + expected mat.SparseMat + }{ + { //Hamming 7 + mat.CSRMat(3, 7, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1), + mat.CSRMat(3, 4, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1), + }, + { //Random - one linearly dependent row + mat.CSRMat(4, 5, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1), + nil, + }, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + + A, _ := ExtractAFromH(context.Background(), test.input, 0) + + if test.expected != nil { + if !test.expected.Equals(A) { + t.Fatalf("expected \n%v\n but found \n%v\n", test.expected, A) + } + } else { + if A != nil { + t.Fatalf("expecte nil but found \n%v\n", A) + } + } + }) + } +} diff --git a/linearblock/ldpc/gallager/gallager.go b/linearblock/ldpc/gallager/gallager.go new file mode 100644 index 0000000..23fd237 --- /dev/null +++ b/linearblock/ldpc/gallager/gallager.go @@ -0,0 +1,132 @@ +package gallager + +import ( + "context" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/linearblock" + "github.com/nathanhack/errorcorrectingcodes/linearblock/internal" + mat "github.com/nathanhack/sparsemat" + "github.com/sirupsen/logrus" + "math/rand" +) + +//GallagerRateInput takes in the message input size in bits (m), the column weight (wc), and row weight (wr) +func Search(ctx context.Context, m, wc, wr, smallestCycleAllowed, maxIter, threads int) (lb *linearblock.LinearBlock, err error) { + if 3 > wc { + return nil, fmt.Errorf("wc must be greater than or equal to 3") + } + if wc >= wr { + return nil, fmt.Errorf("wc (%v) must be less than wr (%v)", wc, wr) + } + if m%wc != 0 { + return nil, fmt.Errorf("wc (%v) must divide m (%v)", wc, m) + } + if smallestCycleAllowed%2 != 0 { + return nil, fmt.Errorf("smallestCycle must be an even number") + } + if smallestCycleAllowed < 4 { + return nil, fmt.Errorf("smallestCycle must at least 4") + } + + N := m / wc * wr + K := m / wc + // Prepare a HPrime used to create all H'subs + HPrime := mat.DOKMat(K, N) + for i := 0; i < K; i++ { + offset := i * wr + for col := 0; col < wr; col++ { + HPrime.Set(i, col+offset, 1) + } + } + iter := maxIter + + for iter > 0 { + iter, lb, err = search(ctx, N, K, m, wc, iter, smallestCycleAllowed, threads, HPrime) + } + return +} + +func search(ctx context.Context, N, K, m, wc, iter, smallestCycleAllowed, threads int, HPrime mat.SparseMat) (int, *linearblock.LinearBlock, error) { + //make the real parity matrix, we'll fill it with the + // correct data next + H := mat.DOKMat(m, N) + + // H is made of three subs + // the first sub == HPrime + // the others are permuted(HPrime) + // we're guaranteed that we can remove + // 4 cycles so if not allow we just + // loop again eventually it'll provide + // permutation that works + + s := 0 + for s < wc && iter > 0 { + iter-- + logrus.Debugf("Iterations remaining %v", iter) + sub := HPrime + if s > 0 { + sub = permuteColumns(HPrime) + } + setSubH(H, sub, s) + + calGirth := linearblock.CalculateGirthLowerBound(H, smallestCycleAllowed, threads) + if -1 < calGirth && calGirth < smallestCycleAllowed { + continue + } + + rank := internal.CalculateRank(ctx, H, threads, false) + if rank != (s+1)*K { + continue + } + s++ + } + if s != wc { + return iter, nil, fmt.Errorf("failed to find a solution") + } + logrus.Debugf("Gallager H Matrix found") + + order, g := internal.NewFromH(ctx, H, threads) + if order == nil { + return iter, nil, fmt.Errorf("unable to create generator for H matrix") + } + + return 0, &linearblock.LinearBlock{ + H: H, + Processing: &linearblock.Systemic{ + HColumnOrder: order, + G: g, + }, + }, nil +} + +func permuteColumns(H mat.SparseMat) mat.SparseMat { + rows, cols := H.Dims() + result := mat.DOKMat(rows, cols) + + //make indices to do permutation + idx := make([]int, cols) + for i := 0; i < cols; i++ { + idx[i] = i + } + + //shuffle them + rand.Shuffle(len(idx), func(i, j int) { + tmp := idx[i] + idx[i] = idx[j] + idx[j] = tmp + }) + + //now set the columns + for i, col := range idx { + tmp := H.Column(col) + result.SetColumn(i, tmp) + } + + return result +} + +func setSubH(H, Hsub mat.SparseMat, index int) { + K, _ := Hsub.Dims() + offset := index * K + H.SetMatrix(Hsub, offset, 0) +} diff --git a/linearblock/linearblock.go b/linearblock/linearblock.go new file mode 100644 index 0000000..72d5f2c --- /dev/null +++ b/linearblock/linearblock.go @@ -0,0 +1,181 @@ +package linearblock + +import ( + "encoding/json" + "fmt" + "github.com/nathanhack/errorcorrectingcodes/linearblock/internal" + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + mat "github.com/nathanhack/sparsemat" + "strings" +) + +type Systemic struct { + HColumnOrder []int + G mat.SparseMat +} + +//LinearBlock contains matrices for the original H matrix and the systemic G generator. +type LinearBlock struct { + H mat.SparseMat //the original H(parity) matrix + Processing *Systemic // contains systemic generator matrix +} + +//// For JSON unmarshalling +type systemic struct { + HColumnOrder []int + G mat.CSRMatrix +} +type linearblock struct { + H mat.CSRMatrix + Processing *systemic +} + +//UnmarshalJSON is needed because LinearBlock has a mat.SparseMat and requires special handling +func (l *LinearBlock) UnmarshalJSON(bytes []byte) error { + var lb linearblock + err := json.Unmarshal(bytes, &lb) + if err != nil { + return err + } + + l.H = &lb.H + if lb.Processing == nil { + return nil + } + + l.Processing = &Systemic{ + HColumnOrder: lb.Processing.HColumnOrder, + G: &lb.Processing.G, + } + + return nil +} + +//Encode take in a message and encodes it using the linear block, returning a codeword +func (l *LinearBlock) Encode(message mat.SparseVector) (codeword mat.SparseVector) { + G := l.Processing.G + rows, cols := G.Dims() + if message.Len() != rows { + panic(fmt.Sprintf("message length == %v is required but found %v", rows, message.Len())) + } + + codeword = mat.DOKVec(cols) + codeword.MulMat(message, G) + + return unorderVector(codeword, l.Processing.HColumnOrder) +} + +//EncodeBE is encode for Binary Erasure channels +func (l *LinearBlock) EncodeBE(message mat.SparseVector) (codeword []bec.ErasureBit) { + tmp := l.Encode(message) + + codeword = make([]bec.ErasureBit, tmp.Len()) + for i := 0; i < tmp.Len(); i++ { + codeword[i] = bec.ErasureBit(tmp.At(i)) + } + return codeword +} + +func unorderVector(codeword mat.SparseVector, ordering []int) mat.SparseVector { + if len(ordering) > 0 && codeword.Len() != len(ordering) { + panic("vector length must equal ordering length") + } + result := mat.DOKVec(codeword.Len()) + + for c, c1 := range ordering { + result.Set(c1, codeword.At(c)) + } + + return result +} + +func orderVector(codeword mat.SparseVector, ordering []int) mat.SparseVector { + if len(ordering) > 0 && codeword.Len() != len(ordering) { + panic("vector length must equal ordering length") + } + result := mat.DOKVec(codeword.Len()) + + for c, c1 := range ordering { + result.Set(c, codeword.At(c1)) + } + + return result +} + +func orderVectorBE(codeword []bec.ErasureBit, ordering []int) []bec.ErasureBit { + if len(ordering) == 0 { + panic("ordering length must be >0") + } + if len(ordering) > 0 && len(codeword) != len(ordering) { + panic("vector length must equal ordering length") + } + result := make([]bec.ErasureBit, len(codeword)) + + for c, c1 := range ordering { + result[c] = codeword[c1] + } + + return result +} + +//Decode takes in a codeword and returns the message contained in it +func (l *LinearBlock) Decode(codeword mat.SparseVector) (message mat.SparseVector) { + if codeword.Len() != l.CodewordLength() { + panic(fmt.Sprintf("codeword length == %v required but found %v", l.CodewordLength(), codeword.Len())) + } + + ml := l.MessageLength() + + codeword = orderVector(codeword, l.Processing.HColumnOrder) + return codeword.Slice(0, ml) +} + +func (l *LinearBlock) DecodeBE(codeword []bec.ErasureBit) (message []bec.ErasureBit) { + if len(codeword) != l.CodewordLength() { + panic(fmt.Sprintf("codeword length == %v required but found %v", l.CodewordLength(), len(codeword))) + } + + ml := l.MessageLength() + + codeword = orderVectorBE(codeword, l.Processing.HColumnOrder) + return codeword[0:ml] +} + +func (l *LinearBlock) Syndrome(codeword mat.SparseVector) (syndrome mat.SparseVector) { + syndrome = mat.CSRVec(l.ParitySymbols()) + syndrome.MatMul(l.H, codeword) + return +} + +func (l *LinearBlock) MessageLength() int { + k, _ := l.Processing.G.Dims() + return k +} +func (l *LinearBlock) ParitySymbols() int { + m, _ := l.H.Dims() + return m +} +func (l *LinearBlock) CodewordLength() int { + _, n := l.H.Dims() + return n +} +func (l *LinearBlock) CodeRate() float64 { + return float64(l.MessageLength()) / float64(l.CodewordLength()) +} + +//Validate will test if this linearblock satisfies G*H.T=0, where G is the generator matrix and H.T is the transpose of H +func (l *LinearBlock) Validate() bool { + //now we validate it + return internal.ValidateHGMatrices(l.Processing.G, internal.ColumnSwapped(l.H, l.Processing.HColumnOrder)) +} + +func (l *LinearBlock) String() string { + buf := strings.Builder{} + buf.WriteString("{\nH:\n") + buf.WriteString(l.H.String()) + buf.WriteString(fmt.Sprintf("Order: %v", l.Processing.HColumnOrder)) + buf.WriteString("\nG:\n") + buf.WriteString(l.Processing.G.String()) + buf.WriteString("\n}\n") + return buf.String() +} diff --git a/linearblock/linearblock_test.go b/linearblock/linearblock_test.go new file mode 100644 index 0000000..1b58add --- /dev/null +++ b/linearblock/linearblock_test.go @@ -0,0 +1,31 @@ +package linearblock + +import ( + mat "github.com/nathanhack/sparsemat" + "math/rand" + "reflect" + "testing" +) + +func TestOrderUnorderVector(t *testing.T) { + vec := mat.DOKVec(100) + columns := make([]int, vec.Len()) + for i := 0; i < vec.Len(); i++ { + vec.Set(i, rand.Intn(2)) + columns[i] = i + } + + rand.Shuffle(len(columns), func(i, j int) { + tmp := columns[i] + columns[i] = columns[j] + columns[j] = tmp + }) + + swapped := orderVector(vec, columns) + + actual := unorderVector(swapped, columns) + + if !reflect.DeepEqual(vec, actual) { + t.Fatalf("expected %v but found %v", vec, actual) + } +} diff --git a/linearblock/messagepassing/bec/bec.go b/linearblock/messagepassing/bec/bec.go new file mode 100644 index 0000000..df4d4a9 --- /dev/null +++ b/linearblock/messagepassing/bec/bec.go @@ -0,0 +1,25 @@ +package bec + +type ErasureBit int + +const ( + Zero ErasureBit = iota + One + Erased +) + +type BECFlippingAlg interface { + Flip(currentCodeword []ErasureBit) (nextCodeword []ErasureBit, done bool) +} + +func Flipping(alg BECFlippingAlg, codeword []ErasureBit) (result []ErasureBit) { + done := false + result = make([]ErasureBit, len(codeword)) + for i := 0; i < len(codeword); i++ { + result[i] = codeword[i] + } + for !done { + result, done = alg.Flip(result) + } + return result +} diff --git a/linearblock/messagepassing/bec/iterative/simple.go b/linearblock/messagepassing/bec/iterative/simple.go new file mode 100644 index 0000000..4d8e0c6 --- /dev/null +++ b/linearblock/messagepassing/bec/iterative/simple.go @@ -0,0 +1,123 @@ +package iterative + +import ( + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + mat "github.com/nathanhack/sparsemat" +) + +type Simple struct { + H mat.SparseMat + cache [][]int +} + +func (s *Simple) Flip(currentCodeword []bec.ErasureBit) (nextCodeword []bec.ErasureBit, done bool) { + if s.H == nil { + panic("Simple BEC flipping algorithm must have the H parity matrix set before using") + } + if s.cache == nil { + s.init() + } + + nextCodeword = copyErasureBits(currentCodeword) + + if !hasErasedBits(nextCodeword) { + return nextCodeword, true + } + + progress := false + for _, B := range s.cache { + if progressM(nextCodeword, B) { + progress = true + } + } + + done = !(progress && hasErasedBits(currentCodeword)) + return nextCodeword, done +} + +func (s *Simple) init() { + if s.cache != nil { + return + } + rows, _ := s.H.Dims() + s.cache = make([][]int, rows) + for i := range s.cache { + s.cache[i] = s.H.Row(i).NonzeroArray() + } +} + +func copyErasureBits(m []bec.ErasureBit) []bec.ErasureBit { + M := make([]bec.ErasureBit, len(m)) + for i, r := range m { + M[i] = r + } + return M +} + +////BEC using the linear block it will try and recover any erased (Erased) bits. All received values of +//// One(1) or Zero(0) will be assumed correct. Every element in received is expected to contain one +//// of the following: One(1), Zero(0), or Erased(2). Any other values will have undetermined csv. +//func BEC(lb *linearblock.LinearBlock, received []bec.ErasureBit) []bec.ErasureBit { +// +// //M will be used to hold the message,and will be ultimately returned +// M := make([]bec.ErasureBit, len(received)) +// for i, r := range received { +// M[i] = r +// } +// +// rows, _ := lb.H.Dims() +// cache := make([][]int, rows) +// for i := range cache { +// cache[i] = lb.H.Row(i).NonzeroArray() +// } +// +// progress := true +// for progress && hasErasedBits(M) { +// progress = false +// for i := 0; i < rows; i++ { +// B := cache[i] +// if progressM(M, B) { +// progress = true +// } +// } +// } +// +// //we make no claim on what is in M and return as is +// // hopefully some or all of the erased bit were recovered +// return M +//} + +func progressM(M []bec.ErasureBit, B []int) bool { + count := 0 + missing := -1 + value := 0 + for _, b := range B { + if M[b] != bec.Erased { + value += int(M[b]) + continue + } + + count++ + missing = b + //we can only fix a check node with only 1 missing value + if count > 1 { + return false + } + } + + if count != 1 { + return false + } + M[missing] = bec.ErasureBit(value % 2) + + return true +} + +func hasErasedBits(M []bec.ErasureBit) bool { + for i := 0; i < len(M); i++ { + if M[i] == bec.Erased { + return true + } + } + return false +} diff --git a/linearblock/messagepassing/bec/iterative/simple_test.go b/linearblock/messagepassing/bec/iterative/simple_test.go new file mode 100644 index 0000000..ea3b1cc --- /dev/null +++ b/linearblock/messagepassing/bec/iterative/simple_test.go @@ -0,0 +1,33 @@ +package iterative + +import ( + "github.com/nathanhack/errorcorrectingcodes/linearblock/messagepassing/bec" + mat "github.com/nathanhack/sparsemat" + "reflect" + "strconv" + "testing" +) + +func TestLinearBlock_BEC(t *testing.T) { + tests := []struct { + H mat.SparseMat + codeword []bec.ErasureBit + expected []bec.ErasureBit + }{ + {mat.DOKMat(4, 6, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1), []bec.ErasureBit{0, 0, 1, bec.Erased, bec.Erased, bec.Erased}, []bec.ErasureBit{0, 0, 1, 0, 1, 1}}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + alg := &Simple{ + H: test.H, + } + + actual := bec.Flipping(alg, test.codeword) + if !reflect.DeepEqual(actual, test.expected) { + t.Log("H") + t.Log(test.H) + t.Fatalf("expected %v but found %v", test.expected, actual) + } + }) + } +} diff --git a/linearblock/messagepassing/bitflipping/harddecision/bitflipping.go b/linearblock/messagepassing/bitflipping/harddecision/bitflipping.go new file mode 100644 index 0000000..6580144 --- /dev/null +++ b/linearblock/messagepassing/bitflipping/harddecision/bitflipping.go @@ -0,0 +1,22 @@ +package harddecision + +import ( + mat "github.com/nathanhack/sparsemat" +) + +type BitFlippingAlg interface { + Flip(currentSyndromes mat.SparseVector, currentCodeword mat.SparseVector) (nextCodeword mat.SparseVector, done bool) + Reset() //resets internal state for next codeword +} + +func BitFlipping(bitFlippingAlg BitFlippingAlg, H mat.SparseMat, codeword mat.SparseVector, maxIter int) (result mat.SparseVector) { + done := false + rows, _ := H.Dims() + result = mat.CSRVecCopy(codeword) + syndrome := mat.CSRVec(rows) + for i := 0; i < maxIter && !done; i++ { + syndrome.MatMul(H, result) + result, done = bitFlippingAlg.Flip(syndrome, result) + } + return result +} diff --git a/linearblock/messagepassing/bitflipping/harddecision/dwbf.go b/linearblock/messagepassing/bitflipping/harddecision/dwbf.go new file mode 100644 index 0000000..eb786bc --- /dev/null +++ b/linearblock/messagepassing/bitflipping/harddecision/dwbf.go @@ -0,0 +1,156 @@ +package harddecision + +import ( + "fmt" + mat "github.com/nathanhack/sparsemat" + mat2 "gonum.org/v1/gonum/mat" +) + +//DWBF_F is a single bit flipping hard decision alg based on the paper +//"Dynamic Weighted Bit-Flipping Decoding Algorithms for LDPC Codes" +// by Tofar C.-Y. Chang and Yu T. Su +type DWBF_F struct { + AlphaFactor float64 //α: 0 < α < 1 + EtaThreshold float64 //η: no requirement but frequently 0.0 is a good value + H mat.SparseMat + z mat.SparseVector //original codeword + r *mat2.Dense + e_n []float64 + columnCache [][]int + rowCache [][]int +} + +func argMaxFloat(values []float64) int { + result := 0 + max := values[0] + for i := 1; i < len(values); i++ { + v := values[i] + if max < v { + result = i + max = v + } + } + return result +} + +func (D *DWBF_F) init(codeword mat.SparseVector) { + if D.AlphaFactor <= 0 || 1 <= D.AlphaFactor { + panic(fmt.Sprintf("0<α<1 is required but found %v ", D.AlphaFactor)) + } + D.z = mat.CSRVecCopy(codeword) + + rows, cols := D.H.Dims() + D.r = mat2.NewDense(rows, cols, nil) + + for i := 0; i < rows; i++ { + for j := 0; j < cols; j++ { + D.r.Set(i, j, 1) + } + } + D.e_n = make([]float64, cols) + + if D.rowCache == nil { + D.rowCache = make([][]int, cols) + for i := 0; i < rows; i++ { + D.rowCache[i] = D.H.Row(i).NonzeroArray() + } + } + + if D.columnCache == nil { + D.columnCache = make([][]int, cols) + for j := 0; j < cols; j++ { + D.columnCache[j] = D.H.Column(j).NonzeroArray() + } + } +} +func (D *DWBF_F) Reset() { + D.z = nil + D.r = nil +} +func (D *DWBF_F) Flip(currentSyndromes mat.SparseVector, currentCodeword mat.SparseVector) (nextCodeword mat.SparseVector, done bool) { + //first we check if the syndromes was zero + if currentSyndromes.IsZero() { + return currentCodeword, true + } + + //internally we keep state + // and if the first time we init() + if D.z == nil || D.r == nil { + D.init(currentCodeword) + } + + // if there are errors then we calculate the + // flipping function E_n vector and r matrix + D.nextE_n(currentSyndromes, currentCodeword) + + // with the updated E_n we now determine the bit set + // B = {n|n arg max_i E_i} + n := argMaxFloat(D.e_n) + + // then let E_n = -E_n + for i, e := range D.e_n { + D.e_n[i] = -e + } + + // and we flip that bit + nextCodeword = mat.CSRVecCopy(currentCodeword) + nextCodeword.Set(n, nextCodeword.At(n)+1) + + D.nextR() + return nextCodeword, false +} + +func (D *DWBF_F) nextE_n(syndromes mat.SparseVector, codeword mat.SparseVector) { + // where l is the iteration + // E^(l)_n = -(1-2*z_n)*(1-2*u_n)-α * sum(r^(l-1)_{mn}*(1-2*s_m), m ∈ M(n)) + + synds := syndromes.NonzeroArray() + syndsLen := len(synds) + for n := 0; n < D.z.Len(); n++ { + sum := 0.0 + cacheRow := D.columnCache[n] + cacheRowLen := len(cacheRow) + for i, j := 0, 0; i < cacheRowLen && j < syndsLen; { + if cacheRow[i] == synds[j] { + sum += -D.r.At(cacheRow[i], n) + i++ + j++ + } else if cacheRow[i] < synds[j] { + sum += D.r.At(cacheRow[i], n) + i++ + } else { + j++ + } + } + D.e_n[n] = -float64((1-2*D.z.At(n))*(1-2*codeword.At(n))) - D.AlphaFactor*sum + } +} + +func (D *DWBF_F) nextR() { + // r^(l)_{mn} = min( thresh(-E^(l)_n'), n' ∈ N(m)\n) + rows, cols := D.r.Dims() + for m := 0; m < rows; m++ { + for n := 0; n < cols; n++ { + min := 0.0 + minIndex := -1 + for _, n1 := range D.rowCache[m] { + if n == n1 { + continue + } + v := threshold(-D.e_n[n1], D.EtaThreshold) + if minIndex == -1 || min > v { + min = v + minIndex = n1 + } + } + D.r.Set(m, n, min) + } + } +} + +func threshold(value, thresh float64) float64 { + if value >= thresh { + return value - thresh + } + return 0 +} diff --git a/linearblock/messagepassing/bitflipping/harddecision/dwdf_test.go b/linearblock/messagepassing/bitflipping/harddecision/dwdf_test.go new file mode 100644 index 0000000..d942508 --- /dev/null +++ b/linearblock/messagepassing/bitflipping/harddecision/dwdf_test.go @@ -0,0 +1,64 @@ +package harddecision + +import ( + "context" + "github.com/nathanhack/errorcorrectingcodes/linearblock/hamming" + mat "github.com/nathanhack/sparsemat" + "strconv" + "testing" +) + +func TestDWDF_BitFlippingHammingCodes(t *testing.T) { + block, err := hamming.New(context.Background(), 3, 0) + if err != nil { + t.Fatalf("expected no error but found: %v", err) + } + tests := []struct { + message mat.SparseVector + flipCodewordBits []int + maxIter int + }{ + {mat.DOKVec(4, 1, 0, 1, 1), []int{0}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{1}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{2}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{3}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{4}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{5}, 20}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + alg := &DWBF_F{ + AlphaFactor: .5, + EtaThreshold: 0, + H: block.H, + } + + //create codeword + codeword := block.Encode(test.message) + expected := mat.CSRVecCopy(codeword) + + for _, index := range test.flipCodewordBits { + codeword.Set(index, codeword.At(index)+1) + } + actual := BitFlipping(alg, block.H, codeword, test.maxIter) + + if !actual.Equals(expected) { + t.Fatalf("expected %v but found %v", expected, actual) + } + }) + } +} + +func BenchmarkDWDF_BitFlipping(b *testing.B) { + h := mat.CSRMat(4, 6, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1) + d := &DWBF_F{ + AlphaFactor: .5, + EtaThreshold: 0, + H: h, + } + input := mat.CSRVec(6, 1, 0, 1, 0, 1, 1) + b.ResetTimer() + for i := 0; i < b.N; i++ { + BitFlipping(d, h, input, 1) + } +} diff --git a/linearblock/messagepassing/bitflipping/harddecision/gallager.go b/linearblock/messagepassing/bitflipping/harddecision/gallager.go new file mode 100644 index 0000000..054a254 --- /dev/null +++ b/linearblock/messagepassing/bitflipping/harddecision/gallager.go @@ -0,0 +1,91 @@ +package harddecision + +import ( + mat "github.com/nathanhack/sparsemat" +) + +func argMaxInt(values []int) int { + result := 0 + max := values[0] + for i := 1; i < len(values); i++ { + v := values[i] + if max < v { + result = i + max = v + } + } + return result +} + +type Gallager struct { + H mat.SparseMat + e_n []int + rowCache [][]int +} + +func (g *Gallager) Reset() { + //nothing to do here +} + +func (g *Gallager) Flip(currentSyndromes mat.SparseVector, currentCodeword mat.SparseVector) (nextCodeword mat.SparseVector, done bool) { + if g.H == nil { + panic("Gallager H matrix must be set before calling Algorithm") + } + //first we check if the syndromes was zero + if currentSyndromes.IsZero() { + return currentCodeword, true + } + + if g.e_n == nil { + g.init(currentCodeword.Len()) + } + + //since there was errors(syndrome is not zero) we'll + // calculate the flipping function E_n vector + g.nextE_n(currentSyndromes) + + n := argMaxInt(g.e_n) + + // and we flip that bit + nextCodeword = mat.CSRVecCopy(currentCodeword) + nextCodeword.Set(n, nextCodeword.At(n)+1) + + return nextCodeword, false +} + +func (g *Gallager) nextE_n(syndromes mat.SparseVector) { + // E_n = -sum((1-2*s_m), m ∈ M(n)) + + synIndices := syndromes.NonzeroArray() + synIndicesLen := len(synIndices) + for n := 0; n < len(g.e_n); n++ { + sum := 0 + indices := g.rowCache[n] + indicesLen := len(indices) + for i, j := 0, 0; i < indicesLen && j < synIndicesLen; { + if indices[i] == synIndices[j] { + sum++ + i++ + j++ + } else if indices[i] < synIndices[j] { + i++ + } else { + j++ + } + } + + g.e_n[n] = -len(indices) + 2*sum + } +} + +func (g *Gallager) init(codewordLen int) { + if g.e_n != nil { + return + } + + g.e_n = make([]int, codewordLen) + g.rowCache = make([][]int, codewordLen) + for n := 0; n < len(g.e_n); n++ { + g.rowCache[n] = g.H.Column(n).NonzeroArray() + } +} diff --git a/linearblock/messagepassing/bitflipping/harddecision/gallager_test.go b/linearblock/messagepassing/bitflipping/harddecision/gallager_test.go new file mode 100644 index 0000000..ccc0b6f --- /dev/null +++ b/linearblock/messagepassing/bitflipping/harddecision/gallager_test.go @@ -0,0 +1,61 @@ +package harddecision + +import ( + "context" + "github.com/nathanhack/errorcorrectingcodes/linearblock/hamming" + mat "github.com/nathanhack/sparsemat" + "strconv" + "testing" +) + +func TestGallager_BitFlippingHammingCodes(t *testing.T) { + block, err := hamming.New(context.Background(), 3, 0) + if err != nil { + t.Fatalf("expected no error but found: %v", err) + } + tests := []struct { + message mat.SparseVector + flipCodewordBits []int + maxIter int + }{ + {mat.DOKVec(4, 1, 0, 1, 1), []int{0}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{1}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{2}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{3}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{4}, 20}, + {mat.DOKVec(4, 1, 0, 1, 1), []int{5}, 20}, + } + for i, test := range tests { + t.Run(strconv.Itoa(i), func(t *testing.T) { + alg := &Gallager{ + H: block.H, + } + + //create codeword + codeword := block.Encode(test.message) + expected := mat.CSRVecCopy(codeword) + + for _, index := range test.flipCodewordBits { + codeword.Set(index, codeword.At(index)+1) + } + + actual := BitFlipping(alg, block.H, codeword, test.maxIter) + + if !actual.Equals(expected) { + t.Fatalf("expected %v but found %v", expected, actual) + } + }) + } +} + +func BenchmarkGallager_BitFlipping(b *testing.B) { + h := mat.CSRMat(4, 6, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1) + g := &Gallager{ + H: h, + } + input := mat.CSRVec(6, 1, 0, 1, 0, 1, 1) + b.ResetTimer() + for i := 0; i < b.N; i++ { + BitFlipping(g, h, input, 1) + } +} diff --git a/main.go b/main.go new file mode 100644 index 0000000..1845072 --- /dev/null +++ b/main.go @@ -0,0 +1,7 @@ +package main + +import "github.com/nathanhack/errorcorrectingcodes/cmd" + +func main() { + cmd.Execute() +}