From 2c16e2c2a38577f7fb88831db1517ad9ce920470 Mon Sep 17 00:00:00 2001
From: Andreas Tsouchlos <an.tsouchlos@gmail.com>
Date: Fri, 4 Nov 2022 22:05:16 +0100
Subject: [PATCH] Implemented projection; Added TODOs and fixed docstrings

---
 sw/decoders/proximal.py  | 21 ++++++++++++++++++++-
 sw/decoders/utility.py   |  6 ++++--
 sw/main.py               |  4 ++--
 sw/test/test_proximal.py |  2 +-
 sw/test/test_utility.py  |  2 +-
 5 files changed, 28 insertions(+), 7 deletions(-)

diff --git a/sw/decoders/proximal.py b/sw/decoders/proximal.py
index 06abbf1..216f817 100644
--- a/sw/decoders/proximal.py
+++ b/sw/decoders/proximal.py
@@ -6,24 +6,30 @@ class ProximalDecoder:
     """Class implementing the Proximal Decoding algorithm. See "Proximal Decoding for LDPC Codes" by Tadashi
     Wadayama, and Satoshi Takabe.
     """
-    def __init__(self, H: np.array, K: int = 10, step_size: float = 0.01, gamma: float = 0.05):
+    # TODO: How large should K be?
+    # TODO: How large should eta be?
+    # TODO: How large should step_size be?
+    def __init__(self, H: np.array, K: int = 100, step_size: float = 0.5, gamma: float = 0.05, eta: float = 1.1):
         """Construct a new ProximalDecoder Object.
 
         :param H: Parity Check Matrix
         :param K: Max number of iterations to perform when decoding
         :param step_size: Step size for the gradient descent process
         :param gamma: Positive constant. Arises in the approximation of the prior PDF
+        :param eta: Positive constant slightly larger than one. See 3.2,  p. 3
         """
         self._H = H
         self._K = K
         self._step_size = step_size
         self._gamma = gamma
+        self._eta = eta
 
     @staticmethod
     def _L_awgn(s: np.array, y: np.array) -> np.array:
         """Variation of the negative log-likelihood for the special case of AWGN noise. See 4.1, p. 4."""
         return s - y
 
+    # TODO: Is this correct?
     def _grad_h(self, x: np.array) -> np.array:
         """Gradient of the code-constraint polynomial. See 2.3, p. 2."""
         # Calculate first term
@@ -56,6 +62,16 @@ class ProximalDecoder:
 
         return np.array(result)
 
+    # TODO: Is this correct?
+    def _projection(self, x):
+        """Project a vector onto [-eta, eta]^n in order to avoid numerical instability.
+        Detailed in 3.2, p. 3 (Equation (15)).
+
+        :param x:
+        :return:
+        """
+        return np.clip(x, -self._eta, self._eta)
+
     def _check_parity(self, y_hat: np.array) -> bool:
         """Perform a parity check for a given codeword.
 
@@ -77,7 +93,10 @@ class ProximalDecoder:
         x_hat = 0
         for k in range(self._K):
             r = s - self._step_size * self._L_awgn(s, y)
+
             s = r - self._gamma * self._grad_h(r)
+            s = self._projection(s)  # Equation (15)
+
             x_hat = (np.sign(s) == 1) * 1
 
             if self._check_parity(x_hat):
diff --git a/sw/decoders/utility.py b/sw/decoders/utility.py
index 0fe573e..60986d3 100644
--- a/sw/decoders/utility.py
+++ b/sw/decoders/utility.py
@@ -7,7 +7,7 @@ from tqdm import tqdm
 
 
 def _get_noise_amp_from_SNR(SNR: float, signal_amp: float = 1) -> float:
-    """Calculate the amplitude of the noise from an SNR and the signal amplitude
+    """Calculate the amplitude of the noise from an SNR and the signal amplitude.
 
     :param SNR: Signal-to-Noise-Ratio in dB
     :param signal_amp: Signal Amplitude (linear)
@@ -49,7 +49,7 @@ def test_decoder(decoder: typing.Any,
         -> typing.Tuple[np.array, np.array]:
     """Calculate the Bit Error Rate (BER) for a given decoder for a number of SNRs.
 
-    This function prints it's progress to stdout
+    This function prints its progress to stdout.
 
     :param decoder: Instance of the decoder to be tested
     :param c: Codeword whose transmission is to be simulated
@@ -69,6 +69,8 @@ def test_decoder(decoder: typing.Any,
                       leave=False,
                       bar_format="{l_bar}{bar}| {n_fmt}/{total_fmt}"):
 
+            # TODO: Is this a valid simulation? Can we just add AWGN to the codeword, ignoring and modulation and (
+            #  e.g. matched) filtering?
             y = add_awgn(c, SNR)
             y_hat = decoder.decode(y)
 
diff --git a/sw/main.py b/sw/main.py
index 29f56e2..d4ef3fe 100644
--- a/sw/main.py
+++ b/sw/main.py
@@ -26,8 +26,8 @@ def main():
 
     print(f"Simulating with c = {c}")
 
-    decoder = proximal.ProximalDecoder(H, K=100)
-    SNRs, BERs = utility.test_decoder(decoder, c, N=100)
+    decoder = proximal.ProximalDecoder(H, K=100, gamma=0.01)
+    SNRs, BERs = utility.test_decoder(decoder, c, SNRs=[1, 3, 20], N=1000)
 
     data = pd.DataFrame({"SNR": SNRs, "BER": BERs})
 
diff --git a/sw/test/test_proximal.py b/sw/test/test_proximal.py
index ed3afb8..d5f6cd6 100644
--- a/sw/test/test_proximal.py
+++ b/sw/test/test_proximal.py
@@ -35,7 +35,7 @@ class CheckParityTestCase(unittest.TestCase):
 
 
 class GradientTestCase(unittest.TestCase):
-    """Test case for the calculation of the gradient of the code-constraint-polynomial"""
+    """Test case for the calculation of the gradient of the code-constraint-polynomial."""
     def test_grad_h(self):
         H = np.array([[1, 0, 0],
                       [0, 1, 0]])
diff --git a/sw/test/test_utility.py b/sw/test/test_utility.py
index 1a2dc35..2c67263 100644
--- a/sw/test/test_utility.py
+++ b/sw/test/test_utility.py
@@ -21,7 +21,7 @@ class CountBitErrorsTestCase(unittest.TestCase):
 
 
 class NoiseAmpFromSNRTestCase(unittest.TestCase):
-    """Test case for noise amplitude calculation"""
+    """Test case for noise amplitude calculation."""
     def test_get_noise_amp_from_SNR(self):
         SNR1 = 0
         SNR2 = 6