LDPC codes with belief propagation decoding
irmin: clean up dead files, document structure
+446
+17
test/interop/ccsds131/dune
+17
test/interop/ccsds131/dune
···
1
+
(test
2
+
(name test)
3
+
(libraries ldpc csvt alcotest)
4
+
(deps
5
+
(source_tree traces)
6
+
(source_tree scripts)))
7
+
8
+
; Regenerate traces against numpy: dune build @regen-traces
9
+
10
+
(rule
11
+
(alias regen-traces)
12
+
(deps
13
+
(source_tree scripts))
14
+
(action
15
+
(chdir
16
+
scripts
17
+
(run bash generate.sh))))
+314
test/interop/ccsds131/scripts/generate.py
+314
test/interop/ccsds131/scripts/generate.py
···
1
+
"""Generate CCSDS AR4JA LDPC interop traces using numpy.
2
+
3
+
The parity check matrix H is constructed from CCSDS 131.0-B-4 Section 7
4
+
(AR4JA rate 1/2, k=1024, n_full=2560, punctured=512, n=2048).
5
+
6
+
The oracle is the syndrome check: for a valid codeword c, H*c mod 2 = 0.
7
+
Reference codewords are derived via GF(2) Gaussian elimination on H,
8
+
which is standard linear algebra -- not a reimplementation of any specific
9
+
encoding algorithm.
10
+
11
+
Traces generated by: numpy (GF(2) linear algebra on CCSDS 131.0-B-4 H matrix)
12
+
Regenerate: dune build @regen-traces
13
+
"""
14
+
15
+
import csv
16
+
import os
17
+
import sys
18
+
19
+
import numpy as np
20
+
21
+
# ---------- CCSDS 131.0-B-4 AR4JA rate 1/2 H matrix construction ----------
22
+
23
+
# Block types
24
+
HZ = 0 # zero block
25
+
HI = 1 # identity block (possibly shifted)
26
+
HP = 2 # permutation block using theta/phi tables
27
+
28
+
# Prototype matrix layers (3 rows x 5 used columns, 3 layers for
29
+
# superimposed edges). Each entry is (block_type, block_value).
30
+
31
+
PROTO_0 = [
32
+
# row 0
33
+
[(HZ, 0), (HZ, 0), (HI, 0), (HZ, 0), (HI, 0)],
34
+
# row 1
35
+
[(HI, 0), (HI, 0), (HZ, 0), (HI, 0), (HP, 1)],
36
+
# row 2
37
+
[(HI, 0), (HP, 4), (HZ, 0), (HP, 6), (HI, 0)],
38
+
]
39
+
40
+
PROTO_1 = [
41
+
# row 0
42
+
[(0, 0), (0, 0), (0, 0), (0, 0), (HP, 0)],
43
+
# row 1
44
+
[(0, 0), (0, 0), (0, 0), (0, 0), (HP, 2)],
45
+
# row 2
46
+
[(0, 0), (HP, 5), (0, 0), (HP, 7), (0, 0)],
47
+
]
48
+
49
+
PROTO_2 = [
50
+
# row 0
51
+
[(0, 0), (0, 0), (0, 0), (0, 0), (0, 0)],
52
+
# row 1
53
+
[(0, 0), (0, 0), (0, 0), (0, 0), (HP, 3)],
54
+
# row 2
55
+
[(0, 0), (0, 0), (0, 0), (0, 0), (0, 0)],
56
+
]
57
+
58
+
# THETA_K: row permutation index (CCSDS 131.0-B-4 Table 7-2)
59
+
THETA_K = [
60
+
3, 0, 1, 2, 2, 3, 0, 1, 0, 1, 2, 0, 2, 3, 0, 1,
61
+
2, 0, 1, 2, 0, 1, 2, 1, 2, 3,
62
+
]
63
+
64
+
# PHI_J_K for circulant_size=128 (CCSDS 131.0-B-4 Table 7-5)
65
+
# phi[j][k] for j in 0..3, k in 0..25
66
+
PHI_J_K = [
67
+
[1, 22, 0, 26, 0, 10, 5, 18, 3, 22, 3, 8, 25, 25, 2, 27, 7, 7, 15, 10, 4, 19, 7, 9, 26, 17],
68
+
[0, 27, 30, 28, 7, 1, 8, 20, 26, 24, 4, 12, 23, 15, 15, 22, 31, 3, 29, 21, 2, 5, 11, 26, 9, 17],
69
+
[0, 12, 30, 18, 10, 16, 13, 9, 7, 15, 16, 18, 4, 23, 5, 3, 29, 11, 4, 8, 2, 11, 11, 3, 15, 13],
70
+
[0, 13, 19, 14, 15, 20, 17, 4, 4, 11, 17, 20, 8, 22, 19, 15, 5, 21, 17, 9, 20, 18, 31, 13, 2, 18],
71
+
]
72
+
73
+
M_SUB = 512 # submatrix size
74
+
CIRC = M_SUB // 4 # circulant size = 128
75
+
LOG_CIRC = 7 # log2(128)
76
+
NUM_ROWS = 3
77
+
NUM_COLS = 5
78
+
H_M = NUM_ROWS * M_SUB # 1536
79
+
H_N = NUM_COLS * M_SUB # 2560 (n_full)
80
+
K = 1024
81
+
PUNCTURED = 512
82
+
N = H_N - PUNCTURED # 2048
83
+
84
+
85
+
def build_h_sparse():
86
+
"""Build the H matrix as a sparse set of (row, col) edges."""
87
+
edges = set()
88
+
89
+
def process_layer(proto):
90
+
for row_idx in range(NUM_ROWS):
91
+
for col_idx in range(NUM_COLS):
92
+
btype, bval = proto[row_idx][col_idx]
93
+
if btype == HI:
94
+
for c in range(M_SUB):
95
+
chk = row_idx * M_SUB + c
96
+
var = col_idx * M_SUB + ((c + bval) & (M_SUB - 1))
97
+
edges.add((chk, var))
98
+
elif btype == HP:
99
+
for c in range(M_SUB):
100
+
j = c >> LOG_CIRC
101
+
pi = (
102
+
(((THETA_K[bval] + j) % 4) << LOG_CIRC)
103
+
+ ((PHI_J_K[j][bval] + c) & (CIRC - 1))
104
+
)
105
+
chk = row_idx * M_SUB + c
106
+
var = col_idx * M_SUB + pi
107
+
edges.add((chk, var))
108
+
109
+
process_layer(PROTO_0)
110
+
process_layer(PROTO_1)
111
+
process_layer(PROTO_2)
112
+
return edges
113
+
114
+
115
+
def build_h_dense(edges):
116
+
"""Build the H matrix as a dense numpy uint8 array."""
117
+
H = np.zeros((H_M, H_N), dtype=np.uint8)
118
+
for r, c in edges:
119
+
H[r, c] = 1
120
+
return H
121
+
122
+
123
+
def gf2_systematic(H, m, n_full, k):
124
+
"""Put H into systematic form [A | I_m] via GF(2) Gaussian elimination.
125
+
126
+
Returns (perm, gen_p) where:
127
+
- perm[i] = original column index at position i after permutation
128
+
- gen_p[j] = list of info-bit indices that XOR to produce parity bit j
129
+
"""
130
+
mat = H.copy()
131
+
perm = list(range(n_full))
132
+
133
+
for pivot in range(m):
134
+
col = k + pivot
135
+
# Find a row >= pivot with a 1 in column col
136
+
found_row = -1
137
+
for r in range(pivot, m):
138
+
if mat[r, col] == 1:
139
+
found_row = r
140
+
break
141
+
142
+
if found_row < 0:
143
+
# Find a column in 0..k-1 that has a 1 in some row >= pivot
144
+
for c in range(k):
145
+
for r in range(pivot, m):
146
+
if mat[r, c] == 1:
147
+
# Swap columns c and col
148
+
mat[:, [c, col]] = mat[:, [col, c]]
149
+
perm[c], perm[col] = perm[col], perm[c]
150
+
found_row = r
151
+
break
152
+
if found_row >= 0:
153
+
break
154
+
155
+
if found_row >= 0:
156
+
# Swap rows
157
+
if found_row != pivot:
158
+
mat[[pivot, found_row]] = mat[[found_row, pivot]]
159
+
# Eliminate
160
+
for r2 in range(m):
161
+
if r2 != pivot and mat[r2, col] == 1:
162
+
mat[r2] ^= mat[pivot]
163
+
164
+
# Extract generator parity sub-matrix
165
+
gen_p = []
166
+
for j in range(m):
167
+
cols = [i for i in range(k) if mat[j, i] == 1]
168
+
gen_p.append(cols)
169
+
170
+
return perm, gen_p
171
+
172
+
173
+
def systematic_encode(info_bits, k, m, n_full, punctured, gen_p):
174
+
"""Systematic LDPC encoding.
175
+
176
+
info_bits: array of k bits (0/1)
177
+
Returns: array of n bits (n_full - punctured), MSB-first packed bytes.
178
+
"""
179
+
n = n_full - punctured
180
+
# Compute parity bits
181
+
parity = np.zeros(m, dtype=np.uint8)
182
+
for j in range(m):
183
+
s = 0
184
+
for i in gen_p[j]:
185
+
s ^= info_bits[i]
186
+
parity[j] = s
187
+
188
+
# Full codeword: [info | parity]
189
+
full = np.zeros(n_full, dtype=np.uint8)
190
+
full[:k] = info_bits
191
+
full[k:k+m] = parity
192
+
193
+
# Puncture: take first n bits
194
+
transmitted = full[:n]
195
+
196
+
# Pack MSB-first into bytes
197
+
n_bytes = (n + 7) // 8
198
+
result = bytearray(n_bytes)
199
+
for i in range(n):
200
+
if transmitted[i]:
201
+
byte_idx = i // 8
202
+
bit_pos = 7 - (i % 8)
203
+
result[byte_idx] |= (1 << bit_pos)
204
+
205
+
return bytes(result)
206
+
207
+
208
+
def verify_syndrome(H, codeword_bits_full):
209
+
"""Verify H * c mod 2 = 0 for the full (unpunctured) codeword."""
210
+
syndrome = H @ codeword_bits_full % 2
211
+
return np.all(syndrome == 0)
212
+
213
+
214
+
# ---------- Test vector generation ----------
215
+
216
+
def make_info_bytes(seed, length=128):
217
+
"""Generate deterministic test data from a seed."""
218
+
return bytes(((i * 37 + seed) * 53 + 7) & 0xFF for i in range(length))
219
+
220
+
221
+
def bytes_to_bits(data, n_bits):
222
+
"""Unpack bytes to bit array (MSB-first)."""
223
+
bits = np.zeros(n_bits, dtype=np.uint8)
224
+
for i in range(min(n_bits, len(data) * 8)):
225
+
byte_idx = i // 8
226
+
bit_pos = 7 - (i % 8)
227
+
bits[i] = (data[byte_idx] >> bit_pos) & 1
228
+
return bits
229
+
230
+
231
+
def main():
232
+
if len(sys.argv) < 2:
233
+
print("Usage: generate.py <trace_dir>", file=sys.stderr)
234
+
sys.exit(1)
235
+
236
+
trace_dir = sys.argv[1]
237
+
os.makedirs(trace_dir, exist_ok=True)
238
+
239
+
print("Building CCSDS AR4JA rate 1/2 H matrix...", file=sys.stderr)
240
+
edges = build_h_sparse()
241
+
H = build_h_dense(edges)
242
+
print(f" H shape: {H.shape}, nnz: {np.sum(H)}", file=sys.stderr)
243
+
244
+
print("Computing systematic form...", file=sys.stderr)
245
+
perm, gen_p = gf2_systematic(H, H_M, H_N, K)
246
+
247
+
# Rebuild H with permuted columns for syndrome checking
248
+
inv_perm = [0] * H_N
249
+
for new_pos, orig in enumerate(perm):
250
+
inv_perm[orig] = new_pos
251
+
H_perm = np.zeros_like(H)
252
+
for r, c in edges:
253
+
H_perm[r, inv_perm[c]] = 1
254
+
255
+
# Test vectors
256
+
vectors = [
257
+
("all_zeros", bytes(128)),
258
+
("all_ones", bytes([0xFF] * 128)),
259
+
("alternating_01", bytes([0x55] * 128)),
260
+
("alternating_10", bytes([0xAA] * 128)),
261
+
("single_bit_0", bytes([0x80] + [0x00] * 127)),
262
+
("single_bit_last", bytes([0x00] * 127 + [0x01])),
263
+
("single_bit_middle", bytes([0x00] * 64 + [0x80] + [0x00] * 63)),
264
+
("seed_13", make_info_bytes(13)),
265
+
("seed_42", make_info_bytes(42)),
266
+
("seed_99", make_info_bytes(99)),
267
+
("seed_128", make_info_bytes(128)),
268
+
("seed_200", make_info_bytes(200)),
269
+
("seed_255", make_info_bytes(255)),
270
+
("counting", bytes(i & 0xFF for i in range(128))),
271
+
("prng_lcg", bytes(((i * 1103515245 + 12345) >> 16) & 0xFF for i in range(128))),
272
+
]
273
+
274
+
codewords_path = os.path.join(trace_dir, "codewords.csv")
275
+
print(f"Generating {len(vectors)} test vectors...", file=sys.stderr)
276
+
277
+
with open(codewords_path, "w", newline="") as f:
278
+
writer = csv.writer(f)
279
+
writer.writerow(["name", "info_hex", "codeword_hex"])
280
+
281
+
for name, info_bytes in vectors:
282
+
info_bits = bytes_to_bits(info_bytes, K)
283
+
codeword_bytes = systematic_encode(info_bits, K, H_M, H_N, PUNCTURED, gen_p)
284
+
285
+
# Verify syndrome: unpack codeword to full n_full bits
286
+
# (punctured bits are zero since they come from encoding)
287
+
cw_bits = bytes_to_bits(codeword_bytes, N)
288
+
full_bits = np.zeros(H_N, dtype=np.uint8)
289
+
full_bits[:N] = cw_bits
290
+
# Recompute punctured bits from encoding
291
+
full_bits[:K] = info_bits
292
+
parity = np.zeros(H_M, dtype=np.uint8)
293
+
for j in range(H_M):
294
+
s = 0
295
+
for i in gen_p[j]:
296
+
s ^= info_bits[i]
297
+
parity[j] = s
298
+
full_bits[K:K+H_M] = parity
299
+
300
+
assert verify_syndrome(H_perm, full_bits), \
301
+
f"Syndrome check failed for {name}"
302
+
303
+
writer.writerow([
304
+
name,
305
+
info_bytes.hex(),
306
+
codeword_bytes.hex(),
307
+
])
308
+
print(f" {name}: ok (syndrome=0)", file=sys.stderr)
309
+
310
+
print(f"Wrote {codewords_path}", file=sys.stderr)
311
+
312
+
313
+
if __name__ == "__main__":
314
+
main()
+13
test/interop/ccsds131/scripts/generate.sh
+13
test/interop/ccsds131/scripts/generate.sh
···
1
+
#!/bin/bash
2
+
set -euo pipefail
3
+
SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)"
4
+
TRACE_DIR="$SCRIPT_DIR/../traces"
5
+
mkdir -p "$TRACE_DIR"
6
+
TRACE_DIR="$(cd "$TRACE_DIR" && pwd)"
7
+
8
+
cd "$SCRIPT_DIR"
9
+
if [ ! -d .venv ]; then
10
+
python3 -m venv .venv
11
+
.venv/bin/pip install -r requirements.txt
12
+
fi
13
+
.venv/bin/python3 generate.py "$TRACE_DIR"
+1
test/interop/ccsds131/scripts/requirements.txt
+1
test/interop/ccsds131/scripts/requirements.txt
···
1
+
numpy==2.2.4
+101
test/interop/ccsds131/test.ml
+101
test/interop/ccsds131/test.ml
···
1
+
(** CCSDS 131.0-B-4 interop tests for ocaml-ldpc.
2
+
3
+
The oracle is numpy performing GF(2) linear algebra on the CCSDS AR4JA
4
+
rate 1/2 parity check matrix H (constructed from the spec tables in
5
+
CCSDS 131.0-B-4 Section 7). Reference codewords satisfy H*c = 0 by
6
+
construction.
7
+
8
+
Traces generated by: numpy 2.2.4 (GF(2) Gaussian elimination on H)
9
+
Regenerate: dune build @regen-traces *)
10
+
11
+
let trace path = Filename.concat "traces" path
12
+
13
+
(* {1 Helpers} *)
14
+
15
+
let hex_to_bytes hex =
16
+
let len = String.length hex / 2 in
17
+
Bytes.init len (fun i ->
18
+
let digit c =
19
+
if c >= '0' && c <= '9' then Char.code c - Char.code '0'
20
+
else if c >= 'a' && c <= 'f' then Char.code c - Char.code 'a' + 10
21
+
else Char.code c - Char.code 'A' + 10
22
+
in
23
+
Char.chr ((digit hex.[i * 2] lsl 4) lor digit hex.[(i * 2) + 1]))
24
+
25
+
let bytes_to_hex b =
26
+
let buf = Buffer.create (Bytes.length b * 2) in
27
+
Bytes.iter
28
+
(fun c -> Buffer.add_string buf (Printf.sprintf "%02x" (Char.code c)))
29
+
b;
30
+
Buffer.contents buf
31
+
32
+
(* {1 Trace parsing} *)
33
+
34
+
type codeword_row = { name : string; info_hex : string; codeword_hex : string }
35
+
36
+
let codeword_codec =
37
+
Csvt.(
38
+
Row.(
39
+
obj (fun name info_hex codeword_hex -> { name; info_hex; codeword_hex })
40
+
|> col "name" string ~enc:(fun r -> r.name)
41
+
|> col "info_hex" string ~enc:(fun r -> r.info_hex)
42
+
|> col "codeword_hex" string ~enc:(fun r -> r.codeword_hex)
43
+
|> finish))
44
+
45
+
let parse_codewords () =
46
+
match Csvt.decode_file codeword_codec (trace "codewords.csv") with
47
+
| Ok rows -> rows
48
+
| Error e ->
49
+
Alcotest.failf "failed to parse codewords.csv: %a" Csvt.pp_error e
50
+
51
+
(* {1 Tests} *)
52
+
53
+
let ldpc = Ldpc.ccsds_rate_1_2
54
+
55
+
let test_encode row () =
56
+
let info = hex_to_bytes row.info_hex in
57
+
let got = Ldpc.encode ldpc info in
58
+
let expected = hex_to_bytes row.codeword_hex in
59
+
if got <> expected then
60
+
Alcotest.failf "%s: encode mismatch\n expected: %s\n got: %s"
61
+
row.name (bytes_to_hex expected) (bytes_to_hex got)
62
+
63
+
let test_decode row () =
64
+
let codeword = hex_to_bytes row.codeword_hex in
65
+
let expected_info = hex_to_bytes row.info_hex in
66
+
match Ldpc.decode ldpc codeword with
67
+
| Ok recovered ->
68
+
if recovered <> expected_info then
69
+
Alcotest.failf "%s: decode mismatch\n expected: %s\n got: %s"
70
+
row.name
71
+
(bytes_to_hex expected_info)
72
+
(bytes_to_hex recovered)
73
+
| Error e -> Alcotest.failf "%s: decode failed: %s" row.name e
74
+
75
+
let test_roundtrip row () =
76
+
let info = hex_to_bytes row.info_hex in
77
+
let codeword = Ldpc.encode ldpc info in
78
+
match Ldpc.decode ldpc codeword with
79
+
| Ok recovered ->
80
+
if recovered <> info then
81
+
Alcotest.failf "%s: roundtrip mismatch\n expected: %s\n got: %s"
82
+
row.name (bytes_to_hex info) (bytes_to_hex recovered)
83
+
| Error e -> Alcotest.failf "%s: roundtrip decode failed: %s" row.name e
84
+
85
+
let () =
86
+
let rows = parse_codewords () in
87
+
Alcotest.run "ldpc-interop"
88
+
[
89
+
( "encode",
90
+
List.map
91
+
(fun r -> Alcotest.test_case r.name `Quick (test_encode r))
92
+
rows );
93
+
( "decode",
94
+
List.map
95
+
(fun r -> Alcotest.test_case r.name `Quick (test_decode r))
96
+
rows );
97
+
( "roundtrip",
98
+
List.map
99
+
(fun r -> Alcotest.test_case r.name `Quick (test_roundtrip r))
100
+
rows );
101
+
]