LDPC codes with belief propagation decoding
irmin: add commit/checkout and proof encode/decode
+375
+17
test/interop/commpy/dune
+17
test/interop/commpy/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 commpy: 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))))
+258
test/interop/commpy/scripts/generate.py
+258
test/interop/commpy/scripts/generate.py
···
1
+
#!/usr/bin/env python3
2
+
"""Generate CCSDS LDPC interop traces for ocaml-ldpc.
3
+
4
+
Oracle: standalone CCSDS AR4JA encoder
5
+
Regenerate: dune build @regen-traces
6
+
7
+
Implements the CCSDS 131.0-B-4 Section 7 AR4JA rate 1/2 LDPC encoder:
8
+
- k=1024 information bits, n=2048 transmitted bits
9
+
- n_full=2560 (512 punctured bits)
10
+
- Circulant-based parity-check matrix with two-step expansion
11
+
- Systematic encoding via Gaussian elimination over GF(2)
12
+
13
+
commpy provides a generic LDPC BP decoder/encoder but does NOT include
14
+
CCSDS AR4JA code constructions, so we implement the full CCSDS H-matrix
15
+
construction and systematic encoder here, matching the OCaml implementation.
16
+
"""
17
+
import os
18
+
import sys
19
+
20
+
import numpy as np
21
+
22
+
23
+
# --- Constants ---
24
+
25
+
M_SUB = 512 # submatrix size
26
+
CIRC = M_SUB // 4 # circulant size = 128
27
+
NUM_ROWS = 3 # prototype check rows
28
+
NUM_COLS = 5 # prototype variable columns (only first 5 non-zero)
29
+
H_M = NUM_ROWS * M_SUB # 1536 check equations
30
+
H_N = NUM_COLS * M_SUB # 2560 full codeword length
31
+
K = 1024 # information bits
32
+
N = 2048 # transmitted bits
33
+
PUNCTURED = 512 # punctured bits
34
+
35
+
HZ = 0 # zero block
36
+
HI = 1 # identity (shifted)
37
+
HP = 2 # permutation block
38
+
39
+
# Prototype layers (3 rows x 11 columns, but only first 5 cols matter)
40
+
# Layer 0
41
+
PROTO_0 = [
42
+
[(HZ,0), (HZ,0), (HI,0), (HZ,0), (HI,0)],
43
+
[(HI,0), (HI,0), (HZ,0), (HI,0), (HP,1)],
44
+
[(HI,0), (HP,4), (HZ,0), (HP,6), (HI,0)],
45
+
]
46
+
47
+
# Layer 1
48
+
PROTO_1 = [
49
+
[(0,0), (0,0), (0,0), (0,0), (HP,0)],
50
+
[(0,0), (0,0), (0,0), (0,0), (HP,2)],
51
+
[(0,0), (HP,5), (0,0), (HP,7), (0,0)],
52
+
]
53
+
54
+
# Layer 2
55
+
PROTO_2 = [
56
+
[(0,0), (0,0), (0,0), (0,0), (0,0)],
57
+
[(0,0), (0,0), (0,0), (0,0), (HP,3)],
58
+
[(0,0), (0,0), (0,0), (0,0), (0,0)],
59
+
]
60
+
61
+
# THETA_K: row permutation index
62
+
THETA_K = [
63
+
3, 0, 1, 2, 2, 3, 0, 1, 0, 1, 2, 0, 2, 3, 0, 1,
64
+
2, 0, 1, 2, 0, 1, 2, 1, 2, 3,
65
+
]
66
+
67
+
# PHI_J_K for M=128 (circulant_size = 32, values mod 32)
68
+
PHI_J_K = [
69
+
[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],
70
+
[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],
71
+
[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],
72
+
[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],
73
+
]
74
+
75
+
76
+
def build_ar4ja_h():
77
+
"""Construct the CCSDS AR4JA parity-check matrix H for rate 1/2, k=1024.
78
+
79
+
Returns a sparse representation: list of (row, col) pairs.
80
+
"""
81
+
edges = [[] for _ in range(H_M)]
82
+
mod_m = M_SUB - 1
83
+
mod_circ = CIRC - 1
84
+
log_circ = 7 # log2(128)
85
+
86
+
def process_layer(proto):
87
+
for row_idx in range(NUM_ROWS):
88
+
for col_idx in range(NUM_COLS):
89
+
btype, bval = proto[row_idx][col_idx]
90
+
if btype == HI:
91
+
for c in range(M_SUB):
92
+
chk = row_idx * M_SUB + c
93
+
var = col_idx * M_SUB + ((c + bval) & mod_m)
94
+
edges[chk].append(var)
95
+
elif btype == HP:
96
+
for c in range(M_SUB):
97
+
j = c >> log_circ # j = c // 128
98
+
pi = (((THETA_K[bval] + j) % 4) << log_circ) + \
99
+
((PHI_J_K[j][bval] + c) & mod_circ)
100
+
chk = row_idx * M_SUB + c
101
+
var = col_idx * M_SUB + pi
102
+
edges[chk].append(var)
103
+
104
+
process_layer(PROTO_0)
105
+
process_layer(PROTO_1)
106
+
process_layer(PROTO_2)
107
+
108
+
# Sort and deduplicate
109
+
h_rows = [sorted(set(row)) for row in edges]
110
+
return h_rows
111
+
112
+
113
+
def systematic_form(h_rows, m, n):
114
+
"""Gaussian elimination over GF(2) to get H into [A | I_m] form.
115
+
116
+
Returns (perm, gen_p) where:
117
+
- perm[i] = original column at position i
118
+
- gen_p[j] = list of data column indices for parity bit j
119
+
"""
120
+
k = n - m
121
+
122
+
# Build dense matrix
123
+
mat = np.zeros((m, n), dtype=np.int8)
124
+
for i, cols in enumerate(h_rows):
125
+
for j in cols:
126
+
mat[i, j] = 1
127
+
128
+
perm = list(range(n))
129
+
130
+
def swap_cols(c1, c2):
131
+
if c1 != c2:
132
+
perm[c1], perm[c2] = perm[c2], perm[c1]
133
+
mat[:, [c1, c2]] = mat[:, [c2, c1]]
134
+
135
+
# Put identity in the last m columns
136
+
for pivot in range(m):
137
+
col = k + pivot
138
+
# Find a row >= pivot with a 1 in column col
139
+
found_row = -1
140
+
for r in range(pivot, m):
141
+
if mat[r, col] == 1:
142
+
found_row = r
143
+
break
144
+
145
+
if found_row < 0:
146
+
# Find a column in 0..k-1 with a 1 in some row >= pivot
147
+
for c in range(k):
148
+
for r in range(pivot, m):
149
+
if mat[r, c] == 1:
150
+
swap_cols(c, col)
151
+
found_row = r
152
+
break
153
+
if found_row >= 0:
154
+
break
155
+
156
+
if found_row >= 0:
157
+
if found_row != pivot:
158
+
mat[[pivot, found_row]] = mat[[found_row, pivot]]
159
+
160
+
# Eliminate
161
+
for r2 in range(m):
162
+
if r2 != pivot and mat[r2, col] == 1:
163
+
mat[r2] ^= mat[pivot]
164
+
165
+
# Extract generator parity sub-matrix
166
+
gen_p = []
167
+
for j in range(m):
168
+
cols = [i for i in range(k) if mat[j, i] == 1]
169
+
gen_p.append(cols)
170
+
171
+
return perm, gen_p
172
+
173
+
174
+
def ldpc_encode(data_bytes, perm, gen_p):
175
+
"""Systematic LDPC encode.
176
+
177
+
Args:
178
+
data_bytes: 128 bytes (k=1024 bits)
179
+
perm: column permutation from systematic_form
180
+
gen_p: generator parity matrix from systematic_form
181
+
182
+
Returns:
183
+
256 bytes (n=2048 transmitted bits)
184
+
"""
185
+
assert len(data_bytes) == K // 8, \
186
+
f"Expected {K // 8} bytes, got {len(data_bytes)}"
187
+
188
+
# Extract k information bits (MSB-first)
189
+
info = np.zeros(K, dtype=int)
190
+
for i in range(K):
191
+
byte_idx = i // 8
192
+
bit_pos = 7 - (i % 8)
193
+
info[i] = (data_bytes[byte_idx] >> bit_pos) & 1
194
+
195
+
# Compute parity bits
196
+
parity = np.zeros(H_M, dtype=int)
197
+
for j in range(H_M):
198
+
s = 0
199
+
for i in gen_p[j]:
200
+
s ^= info[i]
201
+
parity[j] = s
202
+
203
+
# Build full codeword [info | parity]
204
+
full_bits = np.concatenate([info, parity])
205
+
206
+
# Puncture: keep only first N bits
207
+
tx_bits = full_bits[:N]
208
+
209
+
# Pack into bytes (MSB-first)
210
+
result = []
211
+
for i in range(0, N, 8):
212
+
byte = 0
213
+
for j in range(8):
214
+
if i + j < N:
215
+
byte |= int(tx_bits[i + j]) << (7 - j)
216
+
result.append(byte)
217
+
return bytes(result)
218
+
219
+
220
+
# --- Generate test vectors ---
221
+
222
+
TRACE_DIR = (sys.argv[1] if len(sys.argv) > 1
223
+
else os.path.join(os.path.dirname(__file__), "..", "traces"))
224
+
225
+
# Test input patterns
226
+
DATA_BYTES_LEN = K // 8 # 128 bytes
227
+
228
+
inputs = [
229
+
("zeros", bytes(DATA_BYTES_LEN)),
230
+
("ones", bytes([0xFF] * DATA_BYTES_LEN)),
231
+
("alternating_01", bytes([0x55] * DATA_BYTES_LEN)),
232
+
("alternating_10", bytes([0xAA] * DATA_BYTES_LEN)),
233
+
("ascending", bytes([i & 0xFF for i in range(DATA_BYTES_LEN)])),
234
+
]
235
+
236
+
if __name__ == '__main__':
237
+
print("Building CCSDS AR4JA H-matrix...", file=sys.stderr)
238
+
h_rows = build_ar4ja_h()
239
+
print("Computing systematic form...", file=sys.stderr)
240
+
perm, gen_p = systematic_form(h_rows, H_M, H_N)
241
+
print("Generating test vectors...", file=sys.stderr)
242
+
243
+
os.makedirs(TRACE_DIR, exist_ok=True)
244
+
with open(os.path.join(TRACE_DIR, "vectors.csv"), "w") as f:
245
+
f.write("name,input_hex,encoded_hex\n")
246
+
for name, data in inputs:
247
+
coded = ldpc_encode(data, perm, gen_p)
248
+
f.write(f"{name},{data.hex()},{coded.hex()}\n")
249
+
250
+
print(f"Wrote {TRACE_DIR}/vectors.csv ({len(inputs)} vectors)")
251
+
252
+
# Verify systematic property: first 128 bytes of codeword = data
253
+
for name, data in inputs:
254
+
coded = ldpc_encode(data, perm, gen_p)
255
+
if coded[:128] != data:
256
+
print(f"WARNING: {name}: systematic property violated!", file=sys.stderr)
257
+
else:
258
+
print(f" {name}: systematic property OK", file=sys.stderr)
+13
test/interop/commpy/scripts/generate.sh
+13
test/interop/commpy/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"
+87
test/interop/commpy/test.ml
+87
test/interop/commpy/test.ml
···
1
+
(** CCSDS LDPC interop tests for ocaml-ldpc.
2
+
3
+
Traces generated by: Python CCSDS AR4JA encoder (same H-matrix
4
+
construction as the OCaml implementation)
5
+
Regenerate: dune build @regen-traces *)
6
+
7
+
let trace path = Filename.concat "traces" path
8
+
9
+
let bytes_of_hex hex =
10
+
let len = String.length hex / 2 in
11
+
Bytes.init len (fun i ->
12
+
Char.chr (int_of_string ("0x" ^ String.sub hex (i * 2) 2)))
13
+
14
+
let hex_of_bytes b =
15
+
let buf = Buffer.create (Bytes.length b * 2) in
16
+
Bytes.iter
17
+
(fun c -> Buffer.add_string buf (Printf.sprintf "%02x" (Char.code c)))
18
+
b;
19
+
Buffer.contents buf
20
+
21
+
type vector = { name : string; input : bytes; expected_coded : bytes }
22
+
23
+
let vector_codec =
24
+
Csvt.(
25
+
Row.(
26
+
obj (fun name input_hex encoded_hex ->
27
+
{
28
+
name;
29
+
input = bytes_of_hex input_hex;
30
+
expected_coded = bytes_of_hex encoded_hex;
31
+
})
32
+
|> col "name" string ~enc:(fun v -> v.name)
33
+
|> col "input_hex" string ~enc:(fun v -> hex_of_bytes v.input)
34
+
|> col "encoded_hex" string ~enc:(fun v -> hex_of_bytes v.expected_coded)
35
+
|> finish))
36
+
37
+
let parse_vectors () =
38
+
match Csvt.decode_file vector_codec (trace "vectors.csv") with
39
+
| Ok rows -> rows
40
+
| Error e -> Alcotest.failf "CSV: %a" Csvt.pp_error e
41
+
42
+
let ldpc = Ldpc.ccsds_rate_1_2
43
+
44
+
let test_encode vec () =
45
+
let coded = Ldpc.encode ldpc vec.input in
46
+
if coded <> vec.expected_coded then
47
+
Alcotest.failf "%s: encode mismatch\n expected: %s\n got: %s"
48
+
vec.name
49
+
(hex_of_bytes vec.expected_coded)
50
+
(hex_of_bytes coded)
51
+
52
+
let test_roundtrip vec () =
53
+
let coded = Ldpc.encode ldpc vec.input in
54
+
match Ldpc.decode ldpc coded with
55
+
| Ok recovered ->
56
+
if recovered <> vec.input then
57
+
Alcotest.failf "%s: roundtrip mismatch\n expected: %s\n got: %s"
58
+
vec.name (hex_of_bytes vec.input) (hex_of_bytes recovered)
59
+
| Error e -> Alcotest.failf "%s: decode failed: %s" vec.name e
60
+
61
+
let test_systematic vec () =
62
+
let coded = Ldpc.encode ldpc vec.input in
63
+
(* For the rate 1/2 code: first k/8 = 128 bytes should be the data *)
64
+
let k_bytes = 128 in
65
+
let prefix = Bytes.sub coded 0 k_bytes in
66
+
if prefix <> vec.input then
67
+
Alcotest.failf
68
+
"%s: systematic property violated\n data: %s\n prefix: %s" vec.name
69
+
(hex_of_bytes vec.input) (hex_of_bytes prefix)
70
+
71
+
let () =
72
+
let vectors = parse_vectors () in
73
+
Alcotest.run "ldpc-interop"
74
+
[
75
+
( "encode",
76
+
List.map
77
+
(fun v -> Alcotest.test_case v.name `Quick (test_encode v))
78
+
vectors );
79
+
( "systematic",
80
+
List.map
81
+
(fun v -> Alcotest.test_case v.name `Quick (test_systematic v))
82
+
vectors );
83
+
( "roundtrip",
84
+
List.map
85
+
(fun v -> Alcotest.test_case v.name `Quick (test_roundtrip v))
86
+
vectors );
87
+
]