embedded ldpc decoder 50% working

5 files changed, 542 insertions, 62 deletions

diff --git a/controller/fw/Makefile b/controller/fw/Makefile
index e2f0f4c..291f840 100644
--- a/controller/fw/Makefile
+++ b/controller/fw/Makefile
@@ -50,7 +50,7 @@ CFLAGS	+= -nostdlib -ffreestanding
 CFLAGS	+= -mthumb -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16
 CFLAGS	+= -fno-common -ffunction-sections -fdata-sections
 
-INT_CFLAGS	+= -Wall -Wextra -Wshadow -Wimplicit-function-declaration -Wundef
+INT_CFLAGS	+= -Wall -Wextra -Wpedantic -Wshadow -Wimplicit-function-declaration -Wundef
 INT_CFLAGS	+= -Wredundant-decls -Wmissing-prototypes -Wstrict-prototypes
 
 CXXFLAGS	+= -Os -g
@@ -97,6 +97,9 @@ $(BUILDDIR)/tinyaes/aes.o:
 	mkdir -p $(@D)
 	make -C $(@D) -f $(TINYAES_DIR_ABS)/Makefile VPATH=$(TINYAES_DIR_ABS) CFLAGS="$(CFLAGS) -c" CC=$(CC) LD=$(LD) AR=$(AR) aes.o
 
+ldpc_decoder_test.so: ldpc_decoder.c
+	gcc -fPIC -shared -Wall -Wextra -Wpedantic -std=gnu11 -O0 -g -o $@ $^
+
 clean:
 	-rm -r $(BUILDDIR)/src
 	-rm $(BUILDDIR)/$(BINARY)
diff --git a/controller/fw/ldpc_decoder.c b/controller/fw/ldpc_decoder.c
index 5200e8f..10b5a08 100644
--- a/controller/fw/ldpc_decoder.c
+++ b/controller/fw/ldpc_decoder.c
@@ -1,86 +1,160 @@
-"""Decoding module."""
-
-/* H: decoding matrix as shape (m, n) array
- *      m: number of equations represented by H
- *      n: number of code words
- * y: received message in codeword space, length n
- *
- * out: output array of length n where the solutions in codeword space will be written
- *
- * maxiter: Maximum iteration number
- */
-int decode(uint8_t H[], const int m, n, uint8_t y[], uint8_t out[], int maxiter) {
-
-    /* fixme: preprocess the following.
-def bitsandnodes(H):
-    """Return bits and nodes of a parity-check matrix H."""
-    bits_indices, bits = scipy.sparse.find(H)[:2]
-    nodes_indices, nodes = scipy.sparse.find(H.T)[:2]
-    bits_histogram = np.bincount(bits_indices)
-    nodes_histogram = np.bincount(nodes_indices)
-    return bits_histogram, bits, nodes_histogram, nodes
-    bits_hist, bits_values, nodes_hist, nodes_values = utils.bitsandnodes(H)
-    */
 
+#include <stdint.h>
+#include <unistd.h>
+#include <stdbool.h>
+#include <math.h>
+#include <stdio.h>
+
+void inner_logbp(
+        size_t m, size_t n,
+        size_t bits_count, size_t nodes_count, const uint32_t bits_values[], const uint32_t nodes_values[],
+        int8_t Lc[],
+        float Lq[], float Lr[],
+        unsigned int n_iter,
+        float L_posteriori_out[]);
+
+//decode(384, 6, 8, ...)
+int decode(size_t n, size_t nodes_count, size_t bits_count, uint32_t bits[], int8_t y[], int8_t out[], unsigned int maxiter) {
+    const size_t m = n * nodes_count / bits_count;
     float Lq[m*n];
     float Lr[m*n];
     float L_posteriori[n];
 
-    for (int n_iter=0; n_iter<maxiter; n_iter++) {
-        Lq, Lr, L_posteriori = inner_logbp(bits_hist, bits_values, nodes_hist,
-                                            nodes_values, y, Lq, Lr, n_iter, L_posteriori)
+    /* Calculate column bit positions from row bit positions */
+    int32_t bits_transposed[nodes_count * n];
+    for (size_t i=0; i<nodes_count * n; i++)
+        bits_transposed[i] = -1;
+
+    for (size_t i=0; i<m; i++) {
+        for (size_t j=0; j<bits_count; j++) {
+            int32_t *base = bits_transposed + bits[i*bits_count + j] * nodes_count;
+            for (; *base != -1; base++)
+                ;
+            *base = i;
+        }
+    }
+
+    /*
+    printf("Row positions: [");
+    for (size_t i=0; i<m*bits_count; i++) {
+        if (i)
+            printf(", ");
+        if (i%32 == 0)
+            printf("\n    ");
+        printf("%4d", bits[i]);
+    }
+    printf("\n]\n");
+
+    printf("Column positions: [");
+    for (size_t i=0; i<n*nodes_count; i++) {
+        if (i)
+            printf(", ");
+        if (i%32 == 0)
+            printf("\n    ");
+        printf("%4d", bits_transposed[i]);
+    }
+    printf("\n]\n");
+    */
+
+    /* Run iterative optimization algorithm */
+    for (unsigned int n_iter=0; n_iter<maxiter; n_iter++) {
+        inner_logbp(m, n, bits_count, nodes_count, bits, (uint32_t*)bits_transposed, y, Lq, Lr, n_iter, L_posteriori);
+
+        float *arrs[3] = {Lq, Lr, L_posteriori};
+        const char *names[3] = {"Lq", "Lr", "L_posteriori"};
+        size_t lens[3] = {m*n, m*n, n};
+        const size_t head_tail = 10;
+        for (int j=0; j<3; j++) {
+            printf("%s=[", names[j]);
+            bool ellipsis = false;
+            const int w = 16;
+            for (size_t i=0; i<lens[j]; i++) {
+                if (lens[j] > 1000 && i/w > head_tail && i/w < m*n/w-head_tail) {
+                    if (!ellipsis) {
+                        ellipsis = true;
+                        printf("\n    ...");
+                    }
+                    continue;
+                }
+                if (i)
+                    printf(", ");
+                if (i%w == 0)
+                    printf("\n    ");
+                float outf = arrs[j][i];
+                char *s = outf < 0 ? "\033[91m" : (outf > 0 ? "\033[92m" : "\033[94m");
+                printf("%s% 012.6g\033[38;5;240m", s, outf);
+            }
+            printf("\n]\n");
+        }
+
         for (size_t i=0; i<n; i++)
             out[i] = L_posteriori[i] <= 0.0f;
 
         for (size_t i=0; i<m; i++) {
             bool sum = 0;
-
-            for (size_t j=0; j<n; j++)
-                sum ^= H[i*n + j] * out[j];
-
+            for (size_t j=0; j<bits_count; j++)
+                sum ^= out[bits[i*bits_count + j]];
             if (sum)
                 continue;
         }
 
-        return 1;
+        fflush(stdout);
+        return n_iter;
     }
-    
+
+    fflush(stdout);
     return -1;
 }
 
 /* Perform inner ext LogBP solver */
 void inner_logbp(
-        size_t m, n,
-        int bits_hist[], bits_values[], nodes_hist[], nodes_values[],
-        uint8_t Lc[],
-        float Lq[], Lr[],
-        int n_iter,
+        size_t m, size_t n,
+        size_t bits_count, size_t nodes_count, uint32_t const bits_values[], const uint32_t nodes_values[],
+        int8_t Lc[],
+        float Lq[], float Lr[],
+        unsigned int n_iter,
         float L_posteriori_out[]) {
 
+    /*
+    printf("Input data: [");
+    for (size_t i=0; i<n; i++) {
+        if (i)
+            printf(", ");
+        if (i%32 == 0)
+            printf("\n    ");
+        printf("%4d", Lc[i]);
+    }
+    printf("\n]\n");
+    */
+
     /* step 1 : Horizontal */
-    int bits_counter = 0;
-    int nodes_counter = 0;
+    unsigned int bits_counter = 0;
     for (size_t i=0; i<m; i++) {
-        int ff = bits_hist[i];
-        bits_counter += ff;
-        
-        for (size_t p=bits_counter; p<bits_counter+ff; p++) {
-            int j = bits_values[p];
+        printf("=== i=%zu\n", i);
+        for (size_t p=bits_counter; p<bits_counter+bits_count; p++) {
+            size_t j = bits_values[p];
+            printf("\033[38;5;240mj=%04zd ", j);
 
             float x = 1;
             if (n_iter == 0) {
-                for (size_t q=bits_counter; q<bits_counter+ff; q++) {
-                    if (bits_values[q] != j)
+                for (size_t q=bits_counter; q<bits_counter+bits_count; q++) {
+                    if (bits_values[q] != j) {
+                        int lcv = Lc[bits_values[q]];
+                        char *s = lcv < 0 ? "\033[91m" : (lcv > 0 ? "\033[92m" : "\033[94m");
+                        printf("nij=%04u Lc=%s%3d\033[38;5;240m ", bits_values[q], s, lcv);
                         x *= tanhf(0.5f * Lc[bits_values[q]]);
+                    }
                 }
 
             } else {
-                for (size_t q=bits_counter; q<bits_counter+ff; q++) {
+                for (size_t q=bits_counter; q<bits_counter+bits_count; q++) {
                     if (bits_values[q] != j)
-                        x *= tanhf(0.5f * Lq[i, bits_values[q]]);
+                        x *= tanhf(0.5f * Lq[i*n + bits_values[q]]);
                 }
             }
 
+            printf("\n==== i=%03zd p=%01zd x=%08f\n", i, p-bits_counter, x);
+
             float num = 1 + x;
             float denom = 1 - x;
             if (num == 0)
@@ -90,40 +164,40 @@ void inner_logbp(
             else
                 Lr[i*n + j] = logf(num/denom);
         }
+
+        bits_counter += bits_count;
     }
 
     /* step 2 : Vertical */
+    unsigned int nodes_counter = 0;
     for (size_t j=0; j<n; j++) {
-        ff = nodes_hist[j];
-
-        for (size_t p=bits_counter; p<nodes_counter+ff; p++) {
-            int i = nodes_values[p];
+        for (size_t p=bits_counter; p<nodes_counter+nodes_count; p++) {
+            size_t i = nodes_values[p];
 
-            Lq[i, j] = Lc[j]
+            Lq[i*n + j] = Lc[j];
 
-            for (size_t q=bits_counter; q<nodes_counter+ff; q++) {
+            for (size_t q=bits_counter; q<nodes_counter+nodes_count; q++) {
                 if (nodes_values[q] != i)
-                    Lq[i, j] += Lr[nodes_values[q], j];
+                    Lq[i*n + j] += Lr[nodes_values[q]*n + j];
             }
         }
 
-        nodes_counter += ff;
+        nodes_counter += nodes_count;
     }
 
     /* LLR a posteriori */
-    size_t nodes_counter = 0;
+    nodes_counter = 0;
     for (size_t j=0; j<n; j++) {
-        size_t ff = nodes_hist[j];
         float sum = 0;
-        for (size_t k=bits_counter; k<nodes_counter+ff; k++)
+        for (size_t k=bits_counter; k<nodes_counter+nodes_count; k++)
             sum += Lr[nodes_values[k]*n + j];
-        nodes_counter += ff;
+        nodes_counter += nodes_count;
 
         L_posteriori_out[j] = Lc[j] + sum;
     }
 }
 
-
+/* TODO
 def get_message(tG, x):
     """Compute the original `n_bits` message from a `n_code` codeword `x`.
 
@@ -150,3 +224,4 @@ def get_message(tG, x):
                                           message[list(range(i+1, k))])
 
     return abs(message)
+*/
diff --git a/controller/fw/ldpc_decoder_test.py b/controller/fw/ldpc_decoder_test.py
new file mode 100644
index 0000000..d5cbb24
--- /dev/null
+++ b/controller/fw/ldpc_decoder_test.py
@@ -0,0 +1,53 @@
+
+import pyldpc
+import scipy.sparse
+import numpy as np
+import test_decoder
+import os, sys
+import ctypes as C
+import argparse
+
+if __name__ != '__main__':
+    raise RuntimeError("Please don't import this module, this is a command-line program.")
+
+parser = argparse.ArgumentParser()
+parser.add_argument('-r', '--reference', action='store_true', default=False, help='Run reference decoder instead of C implemention')
+args = parser.parse_args()
+
+lib = C.CDLL('./ldpc_decoder_test.so')
+
+n = 384
+nodes, bits = 6, 8
+H, G = pyldpc.make_ldpc(n, nodes, bits, systematic=False, seed=0)
+_1, bits_pos, _2 = scipy.sparse.find(H)
+_, k = G.shape
+
+st = np.random.RandomState(seed=0)
+test_data = st.randint(0, 2, k)
+d = np.dot(G, test_data) % 2
+x = (-1) ** d
+
+bits_pos = bits_pos.astype(np.uint32)
+x = x.astype(np.int8)
+
+lib.decode.argtypes = [C.c_size_t, C.c_size_t, C.c_size_t, C.POINTER(C.c_size_t), C.POINTER(C.c_int8), C.POINTER(C.c_int8), C.c_uint]
+
+if args.reference:
+    ref_out = test_decoder.decode(H, x, 3)
+    print('decoder output:', ref_out, flush=True)
+    print('msg reconstruction:', pyldpc.get_message(G, ref_out))
+    print('reference decoder: ', np.all(np.equal(pyldpc.get_message(G, ref_out), test_data)), flush=True)
+    print(ref_out.dtype)
+
+else:
+    out = np.zeros(n, dtype=np.uint8)
+    # print('python data:', x, flush=True)
+    print('decoder iterations:', lib.decode(n, nodes, bits,
+        bits_pos.ctypes.data_as(C.POINTER(C.c_ulong)),
+        x.ctypes.data_as(C.POINTER(C.c_int8)),
+        out.ctypes.data_as(C.POINTER(C.c_int8)),
+        25), flush=True)
+    print('decoder output:', out)
+    print('msg reconstruction:', pyldpc.get_message(G, out))
+    print('decoder under test:', np.all(np.equal(pyldpc.get_message(G, out.astype(np.int64)), test_data)))
+    print(out.dtype)
diff --git a/controller/fw/test_decoder.py b/controller/fw/test_decoder.py
new file mode 100644
index 0000000..806125b
--- /dev/null
+++ b/controller/fw/test_decoder.py
@@ -0,0 +1,167 @@
+"""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:
+            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)
diff --git a/controller/fw/test_pyldpc_utils.py b/controller/fw/test_pyldpc_utils.py
new file mode 100644
index 0000000..6b14532
--- /dev/null
+++ b/controller/fw/test_pyldpc_utils.py
@@ -0,0 +1,182 @@
+"""Conversion tools."""
+import math
+import numbers
+import numpy as np
+import scipy
+from scipy.stats import norm
+pi = math.pi
+
+
+def int2bitarray(n, k):
+    """Change an array's base from int (base 10) to binary (base 2)."""
+    binary_string = bin(n)
+    length = len(binary_string)
+    bitarray = np.zeros(k, 'int')
+    for i in range(length - 2):
+        bitarray[k - i - 1] = int(binary_string[length - i - 1])
+
+    return bitarray
+
+
+def bitarray2int(bitarray):
+    """Change array's base from binary (base 2) to int (base 10)."""
+    bitstring = "".join([str(i) for i in bitarray])
+
+    return int(bitstring, 2)
+
+
+def binaryproduct(X, Y):
+    """Compute a matrix-matrix / vector product in Z/2Z."""
+    A = X.dot(Y)
+    try:
+        A = A.toarray()
+    except AttributeError:
+        pass
+    return A % 2
+
+
+def gaussjordan(X, change=0):
+    """Compute the binary row reduced echelon form of X.
+
+    Parameters
+    ----------
+    X: array (m, n)
+    change : boolean (default, False). If True returns the inverse transform
+
+    Returns
+    -------
+    if `change` == 'True':
+        A: array (m, n). row reduced form of X.
+        P: tranformations applied to the identity
+    else:
+        A: array (m, n). row reduced form of X.
+
+    """
+    A = np.copy(X)
+    m, n = A.shape
+
+    if change:
+        P = np.identity(m).astype(int)
+
+    pivot_old = -1
+    for j in range(n):
+        filtre_down = A[pivot_old+1:m, j]
+        pivot = np.argmax(filtre_down)+pivot_old+1
+
+        if A[pivot, j]:
+            pivot_old += 1
+            if pivot_old != pivot:
+                aux = np.copy(A[pivot, :])
+                A[pivot, :] = A[pivot_old, :]
+                A[pivot_old, :] = aux
+                if change:
+                    aux = np.copy(P[pivot, :])
+                    P[pivot, :] = P[pivot_old, :]
+                    P[pivot_old, :] = aux
+
+            for i in range(m):
+                if i != pivot_old and A[i, j]:
+                    if change:
+                        P[i, :] = abs(P[i, :]-P[pivot_old, :])
+                    A[i, :] = abs(A[i, :]-A[pivot_old, :])
+
+        if pivot_old == m-1:
+            break
+
+    if change:
+        return A, P
+    return A
+
+
+def binaryrank(X):
+    """Compute rank of a binary Matrix using Gauss-Jordan algorithm."""
+    A = np.copy(X)
+    m, n = A.shape
+
+    A = gaussjordan(A)
+
+    return sum([a.any() for a in A])
+
+
+def f1(y, sigma):
+    """Compute normal density N(1,sigma)."""
+    f = norm.pdf(y, loc=1, scale=sigma)
+    return f
+
+
+def fm1(y, sigma):
+    """Compute normal density N(-1,sigma)."""
+
+    f = norm.pdf(y, loc=-1, scale=sigma)
+    return f
+
+
+def bitsandnodes(H):
+    """Return bits and nodes of a parity-check matrix H."""
+    if type(H) != scipy.sparse.csr_matrix:
+        bits_indices, bits = np.where(H)
+        nodes_indices, nodes = np.where(H.T)
+    else:
+        bits_indices, bits = scipy.sparse.find(H)[:2]
+        nodes_indices, nodes = scipy.sparse.find(H.T)[:2]
+    bits_histogram = np.bincount(bits_indices)
+    nodes_histogram = np.bincount(nodes_indices)
+
+    return bits_histogram, bits, nodes_histogram, nodes
+
+
+def incode(H, x):
+    """Compute Binary Product of H and x."""
+    return (binaryproduct(H, x) == 0).all()
+
+
+def gausselimination(A, b):
+    """Solve linear system in Z/2Z via Gauss Gauss elimination."""
+    if type(A) == scipy.sparse.csr_matrix:
+        A = A.toarray().copy()
+    else:
+        A = A.copy()
+    b = b.copy()
+    n, k = A.shape
+
+    for j in range(min(k, n)):
+        listedepivots = [i for i in range(j, n) if A[i, j]]
+        if len(listedepivots):
+            pivot = np.min(listedepivots)
+        else:
+            continue
+        if pivot != j:
+            aux = (A[j, :]).copy()
+            A[j, :] = A[pivot, :]
+            A[pivot, :] = aux
+
+            aux = b[j].copy()
+            b[j] = b[pivot]
+            b[pivot] = aux
+
+        for i in range(j+1, n):
+            if A[i, j]:
+                A[i, :] = abs(A[i, :]-A[j, :])
+                b[i] = abs(b[i]-b[j])
+
+    return A, b
+
+
+def check_random_state(seed):
+    """Turn seed into a np.random.RandomState instance
+    Parameters
+    ----------
+    seed : None | int | instance of RandomState
+        If seed is None, return the RandomState singleton used by np.random.
+        If seed is an int, return a new RandomState instance seeded with seed.
+        If seed is already a RandomState instance, return it.
+        Otherwise raise ValueError.
+    """
+    if seed is None or seed is np.random:
+        return np.random.mtrand._rand
+    if isinstance(seed, numbers.Integral):
+        return np.random.RandomState(seed)
+    if isinstance(seed, np.random.RandomState):
+        return seed
+    raise ValueError('%r cannot be used to seed a numpy.random.RandomState'
+                     ' instance' % seed)