Finishing up freq meas

controller/fw/src/test_decoder.py

@@ -0,0 +1,168 @@
+"""Decoding module."""
+import numpy as np
+import warnings
+import test_pyldpc_utils as utils
+
+from numba import njit, int64, types, float64
+
+np.set_printoptions(linewidth=180, threshold=1000, edgeitems=20)
+
+def decode(H, y, snr, maxiter=100):
+    """Decode a Gaussian noise corrupted n bits message using BP algorithm.
+
+    Decoding is performed in parallel if multiple codewords are passed in y.
+
+    Parameters
+    ----------
+    H: array (n_equations, n_code). Decoding matrix H.
+    y: array (n_code, n_messages) or (n_code,). Received message(s) in the
+        codeword space.
+    maxiter: int. Maximum number of iterations of the BP algorithm.
+
+    Returns
+    -------
+    x: array (n_code,) or (n_code, n_messages) the solutions in the
+        codeword space.
+
+    """
+    m, n = H.shape
+
+    bits_hist, bits_values, nodes_hist, nodes_values = utils.bitsandnodes(H)
+
+    var = 10 ** (-snr / 10)
+
+    if y.ndim == 1:
+        y = y[:, None]
+    # step 0: initialization
+
+    Lc = 2 * y / var
+    _, n_messages = y.shape
+
+    Lq = np.zeros(shape=(m, n, n_messages))
+
+    Lr = np.zeros(shape=(m, n, n_messages))
+
+    for n_iter in range(maxiter):
+        #print(f'============================ iteration {n_iter} ============================')
+        Lq, Lr, L_posteriori = _logbp_numba(bits_hist, bits_values, nodes_hist,
+                                            nodes_values, Lc, Lq, Lr, n_iter)
+        #print("Lq=", Lq.flatten())
+        #print("Lr=", Lr.flatten())
+        #print("L_posteriori=", L_posteriori.flatten())
+        #print('L_posteriori=[')
+        #for row in L_posteriori.reshape([-1, 16]):
+        #    for val in row:
+        #        cc = '\033[91m' if val < 0 else ('\033[92m' if val > 0 else '\033[94m')
+        #        print(f"{cc}{val: 012.6g}\033[38;5;240m", end=', ')
+        #    print()
+        x = np.array(L_posteriori <= 0).astype(int)
+
+        product = utils.incode(H, x)
+
+        if product:
+            print(f'found, n_iter={n_iter}')
+            break
+
+    if n_iter == maxiter - 1:
+        warnings.warn("""Decoding stopped before convergence. You may want
+                       to increase maxiter""")
+    return x.squeeze()
+
+
+output_type_log2 = types.Tuple((float64[:, :, :], float64[:, :, :],
+                               float64[:, :]))
+
+
+#@njit(output_type_log2(int64[:], int64[:], int64[:], int64[:], float64[:, :],
+#                       float64[:, :, :],  float64[:, :, :], int64), cache=True)
+def _logbp_numba(bits_hist, bits_values, nodes_hist, nodes_values, Lc, Lq, Lr,
+                 n_iter):
+    """Perform inner ext LogBP solver."""
+    m, n, n_messages = Lr.shape
+    # step 1 : Horizontal
+
+    bits_counter = 0
+    nodes_counter = 0
+    for i in range(m):
+        #print(f'=== i={i}')
+        ff = bits_hist[i]
+        ni = bits_values[bits_counter: bits_counter + ff]
+        bits_counter += ff
+        for j_iter, j in enumerate(ni):
+            nij = ni[:]
+            #print(f'\033[38;5;240mj={j:04d}', end=' ')
+
+            X = np.ones(n_messages)
+            if n_iter == 0:
+                for kk in range(len(nij)):
+                    if nij[kk] != j:
+                        lcv = Lc[nij[kk],0]
+                        lcc = '\033[91m' if lcv < 0 else ('\033[92m' if lcv > 0 else '\033[94m')
+                        #print(f'nij={nij[kk]:04d} Lc={lcc}{lcv:> 8f}\033[38;5;240m', end=' ')
+                        X *= np.tanh(0.5 * Lc[nij[kk]])
+            else:
+                for kk in range(len(nij)):
+                    if nij[kk] != j:
+                        X *= np.tanh(0.5 * Lq[i, nij[kk]])
+            #print(f'\n==== {i:03d} {j_iter:01d} {X[0]:> 8f}')
+            num = 1 + X
+            denom = 1 - X
+            for ll in range(n_messages):
+                if num[ll] == 0:
+                    Lr[i, j, ll] = -1
+                elif denom[ll] == 0:
+                    Lr[i, j, ll] = 1
+                else:
+                    Lr[i, j, ll] = np.log(num[ll] / denom[ll])
+    # step 2 : Vertical
+
+    for j in range(n):
+        ff = nodes_hist[j]
+        mj = nodes_values[bits_counter: nodes_counter + ff]
+        nodes_counter += ff
+        for i in mj:
+            mji = mj[:]
+            Lq[i, j] = Lc[j]
+
+            for kk in range(len(mji)):
+                if mji[kk] != i:
+                    Lq[i, j] += Lr[mji[kk], j]
+
+    # LLR a posteriori:
+    L_posteriori = np.zeros((n, n_messages))
+    nodes_counter = 0
+    for j in range(n):
+        ff = nodes_hist[j]
+        mj = nodes_values[bits_counter: nodes_counter + ff]
+        nodes_counter += ff
+        L_posteriori[j] = Lc[j] + Lr[mj, j].sum(axis=0)
+
+    return Lq, Lr, L_posteriori
+
+
+def get_message(tG, x):
+    """Compute the original `n_bits` message from a `n_code` codeword `x`.
+
+    Parameters
+    ----------
+    tG: array (n_code, n_bits) coding matrix tG.
+    x: array (n_code,) decoded codeword of length `n_code`.
+
+    Returns
+    -------
+    message: array (n_bits,). Original binary message.
+
+    """
+    n, k = tG.shape
+
+    rtG, rx = utils.gausselimination(tG, x)
+
+    message = np.zeros(k).astype(int)
+
+    message[k - 1] = rx[k - 1]
+    for i in reversed(range(k - 1)):
+        message[i] = rx[i]
+        message[i] -= utils.binaryproduct(rtG[i, list(range(i+1, k))],
+                                          message[list(range(i+1, k))])
+
+    return abs(message)