diff options
author | jaseg <git-bigdata-wsl-arch@jaseg.de> | 2020-03-02 19:42:36 +0100 |
---|---|---|
committer | jaseg <git-bigdata-wsl-arch@jaseg.de> | 2020-03-02 19:42:36 +0100 |
commit | ca01d52a8658ed01301f14a4d294097d4db533cf (patch) | |
tree | 28b32499b6a849beda66c6ebc70239fbd9f21530 /controller/fw/src | |
parent | 5effadcbaf66f476c8fffefa2c349676f41c3f52 (diff) | |
download | master-thesis-ca01d52a8658ed01301f14a4d294097d4db533cf.tar.gz master-thesis-ca01d52a8658ed01301f14a4d294097d4db533cf.tar.bz2 master-thesis-ca01d52a8658ed01301f14a4d294097d4db533cf.zip |
Finishing up freq meas
Diffstat (limited to 'controller/fw/src')
-rw-r--r-- | controller/fw/src/freq_meas.c | 94 | ||||
-rw-r--r-- | controller/fw/src/freq_meas.h | 7 | ||||
-rw-r--r-- | controller/fw/src/gold_code.h | 4 | ||||
-rw-r--r-- | controller/fw/src/ldpc_decoder.c | 267 | ||||
-rw-r--r-- | controller/fw/src/ldpc_decoder_test.py | 72 | ||||
-rw-r--r-- | controller/fw/src/test_decoder.py | 168 | ||||
-rw-r--r-- | controller/fw/src/test_pyldpc_utils.py | 182 |
7 files changed, 794 insertions, 0 deletions
diff --git a/controller/fw/src/freq_meas.c b/controller/fw/src/freq_meas.c new file mode 100644 index 0000000..e6b3976 --- /dev/null +++ b/controller/fw/src/freq_meas.c @@ -0,0 +1,94 @@ + +#include <unistd.h> +#include <math.h> + +#include <arm_math.h> +#include <levmarq.h> + +#include "freq_meas.h" +#include "sr_global.h" + + +/* FTT window lookup table defined in generated/fmeas_fft_window.c */ +extern const float * const fmeas_fft_window_table; + +/* jury-rig some definitions for these functions since the ARM headers only export an over-generalized variable bin size + * variant. */ +extern arm_status arm_rfft_32_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_64_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_128_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_256_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_512_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_1024_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_2048_fast_init_f32(arm_rfft_fast_instance_f32 * S); +extern arm_status arm_rfft_4096_fast_init_f32(arm_rfft_fast_instance_f32 * S); + +#define CONCAT(A, B, C) A ## B ## C +#define arm_rfft_init_name(nbits) CONCAT(arm_rfft_, nbits, _fast_init_f32) + +float func_gauss_grad(float *out, float *params, int x, void *userdata); +float func_gauss(float *params, int x, void *userdata); + +int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out) { + int rc; + float in_buf[FMEAS_FFT_LEN]; + float out_buf[FMEAS_FFT_LEN]; + for (size_t i=0; i<FMEAS_FFT_LEN; i++) + in_buf[i] = (float)adc_buf[i] / (float)FMEAS_ADC_MAX * fmeas_fft_window_table[i]; + + arm_rfft_fast_instance_f32 fft_inst; + if ((rc = arm_rfft_init_name(FMEAS_FFT_LEN)(&fft_inst)) != ARM_MATH_SUCCESS) + return rc; + + arm_rfft_fast_f32(&fft_inst, in_buf, out_buf, 0); + +#define FMEAS_FFT_WINDOW_MIN_F 30.0f +#define FMEAS_FFT_WINDOW_MAX_F 70.0f + const float binsize = (float)FMEAS_ADC_SAMPLING_RATE / FMEAS_FFT_LEN; + const int first_bin = (int)(FMEAS_FFT_WINDOW_MIN_F / binsize); + const int last_bin = (int)(FMEAS_FFT_WINDOW_MAX_F / binsize + 0.5f); + const int nbins = last_bin - first_bin + 1; + + /* Copy real values of target data to front of output buffer */ + for (size_t i=0; i<nbins; i++) + out_buf[i] = out_buf[2 * (first_bin + i)]; + + LMstat lmstat; + levmarq_init(&lmstat); + + float a_max = 0.0f; + int i_max = 0; + for (size_t i=0; i<nbins; i++) { + if (out_buf[i] > a_max) { + a_max = out_buf[i]; + i_max = i; + } + } + + float par[3] = { + a_max, i_max, 1.0f + }; + + if (levmarq(3, ¶ms, nbins, out_buf, NULL, func_gauss, func_gauss_grad, NULL, &lmstat)) + return -1; + + *out = (params[1] + first_bin) * binsize; + + return 0; +} + +float func_gauss(float *params, int x, void *userdata) { + UNUSED(userdata); + float a = params[0]; + float mu = params[1]; + float sigma = params[2]; + return a*expf(-arm_power_f32((x-mu), 2.0f/(2.0f*(sigma*sigma)))); +} + +float func_gauss_grad(float *out, float *params, int x, void *userdata) { + UNUSED(userdata); + float a = params[0]; + float mu = params[1]; + float sigma = params[2]; + return -(x-mu) / ( sigma*sigma*sigma * 2.5066282746310002f) * a*expf(-arm_power_f32((x-mu), 2.0f/(2.0f*(sigma*sigma)))); +} diff --git a/controller/fw/src/freq_meas.h b/controller/fw/src/freq_meas.h new file mode 100644 index 0000000..1c083f8 --- /dev/null +++ b/controller/fw/src/freq_meas.h @@ -0,0 +1,7 @@ + +#ifndef __FREQ_MEAS_H__ +#define __FREQ_MEAS_H__ + +int adc_buf_measure_freq(uint16_t adc_buf[FMEAS_FFT_LEN], float *out); + +#endif /* __FREQ_MEAS_H__ */ diff --git a/controller/fw/src/gold_code.h b/controller/fw/src/gold_code.h new file mode 100644 index 0000000..739b477 --- /dev/null +++ b/controller/fw/src/gold_code.h @@ -0,0 +1,4 @@ + +/* header file for generated gold code tables */ + +extern const uint8_t * const gold_code_table; diff --git a/controller/fw/src/ldpc_decoder.c b/controller/fw/src/ldpc_decoder.c new file mode 100644 index 0000000..fe59d77 --- /dev/null +++ b/controller/fw/src/ldpc_decoder.c @@ -0,0 +1,267 @@ + +#include <stdint.h> +#include <unistd.h> +#include <stdbool.h> +#include <math.h> +#include <stdio.h> + + +void gausselimination(size_t n, size_t k, int8_t *A, int8_t *b); + +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]; + + /* 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<bits_count; j++) + sum ^= out[bits[i*bits_count + j]]; + if (sum) + continue; + } + + fflush(stdout); + return n_iter; + } + + fflush(stdout); + return -1; +} + +/* Perform inner ext LogBP solver */ +void inner_logbp( + 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 */ + unsigned int bits_counter = 0; + for (size_t i=0; i<m; i++) { + //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+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+bits_count; q++) { + if (bits_values[q] != j) + 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) + Lr[i*n + j] = -1.0f; + else if (denom == 0) + Lr[i*n + j] = 1.0f; + 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++) { + for (size_t p=bits_counter; p<nodes_counter+nodes_count; p++) { + size_t i = nodes_values[p]; + + Lq[i*n + j] = Lc[j]; + + for (size_t q=bits_counter; q<nodes_counter+nodes_count; q++) { + if (nodes_values[q] != i) + Lq[i*n + j] += Lr[nodes_values[q]*n + j]; + } + } + + nodes_counter += nodes_count; + } + + /* LLR a posteriori */ + nodes_counter = 0; + for (size_t j=0; j<n; j++) { + float sum = 0; + for (size_t k=bits_counter; k<nodes_counter+nodes_count; k++) + sum += Lr[nodes_values[k]*n + j]; + nodes_counter += nodes_count; + + L_posteriori_out[j] = Lc[j] + sum; + } +} + +/* Compute the original (k) bit message from a (n) bit codeword x. + * + * tG: (n, k)-matrix + * x: (n)-vector + * out: (k)-vector + */ +void get_message(size_t n, size_t k, int8_t *tG, int8_t *x, int8_t *out) { + + gausselimination(n, k, tG, x); + + out[k - 1] = x[k - 1]; + for (ssize_t i=k-2; i>=0; i--) { + out[i] = x[i]; + + uint8_t sum = 0; + for (size_t j=i+1; j<k; j++) + sum ^= tG[i*k + j] * out[j]; + + out[i] = !!(out[i] - sum); + } +} + +/* Solve linear system in Z/2Z via Gauss Gauss elimination. + * + * A: (n, k)-matrix + * b: (n)-vector + */ +void gausselimination(size_t n, size_t k, int8_t *A, int8_t *b) { + ssize_t d = k<n ? k : n; + for (ssize_t j=0; j<d; j++) { + + ssize_t pivot = -1; + for (size_t i=j; i<n; i++) { + if (A[i*k + j]) { + pivot = i; + break; + } + } + if (pivot == -1) + continue; + + if (pivot != j) { + for (size_t i=0; i<k; i++) { + int8_t tmp = A[j*k + i]; + A[j*k + i] = A[pivot*k + i]; + A[pivot*k + i] = tmp; + } + + int8_t tmp = b[j]; + b[j] = b[pivot]; + b[pivot] = tmp; + } + + for (size_t i=j+1; i<n; i++) { + if (A[i*k + j]) { + for (size_t p=0; p<k; p++) + A[i*k + p] = !!(A[i*k + p] - A[j*k + p]); + b[i] = !!(b[i] - b[j]); + } + } + } +} + diff --git a/controller/fw/src/ldpc_decoder_test.py b/controller/fw/src/ldpc_decoder_test.py new file mode 100644 index 0000000..3b91bba --- /dev/null +++ b/controller/fw/src/ldpc_decoder_test.py @@ -0,0 +1,72 @@ + +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 = 5*19 +nodes, bits = 17, 19 +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 +x[29:] = 0 + +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] +lib.get_message.argtypes = [C.c_size_t, C.c_size_t, C.POINTER(C.c_int8), C.POINTER(C.c_int8), C.POINTER(C.c_int8)] + +if args.reference: + ref_out = test_decoder.decode(H, x, 3) + print('decoder output:', ref_out, flush=True) + print('msg reconstruction:', test_decoder.get_message(G, ref_out)) + print('reference decoder: ', np.all(np.equal(test_decoder.get_message(G, ref_out), test_data)), flush=True) + np.set_printoptions(linewidth=220) + print(test_data) + print(test_decoder.get_message(G, ref_out)) + print(test_decoder.get_message(G, ref_out) ^ test_data) + +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:', test_decoder.get_message(G, out.astype(np.int64))) + print('decoder under test:', np.all(np.equal(test_decoder.get_message(G, out.astype(np.int64)), test_data))) + np.set_printoptions(linewidth=220) + print(test_data) + print(test_decoder.get_message(G, out.astype(np.int64))) + G = G.astype(np.int8) + msg = np.zeros(k, dtype=np.int8) + lib.get_message( + n, k, + G.ctypes.data_as(C.POINTER(C.c_int8)), + out.astype(np.int8).ctypes.data_as(C.POINTER(C.c_int8)), + msg.ctypes.data_as(C.POINTER(C.c_int8))) + print(msg) + print(msg ^ test_data) + +print('codeword length:', len(x)) +print('data length:', len(test_data)) diff --git a/controller/fw/src/test_decoder.py b/controller/fw/src/test_decoder.py new file mode 100644 index 0000000..8be5b02 --- /dev/null +++ b/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) diff --git a/controller/fw/src/test_pyldpc_utils.py b/controller/fw/src/test_pyldpc_utils.py new file mode 100644 index 0000000..6b14532 --- /dev/null +++ b/controller/fw/src/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) |