ldpc_toolbox/
systematic.rs

1//! Systematic code constructions.
2//!
3//! This module contains a function [`parity_to_systematic`] that can be used to
4//! convert a full-rank parity check matrix into one that supports systematic
5//! encoding using the first variables (as done by the systematic encoder in the
6//! [`encoder`](crate::encoder) module) by permuting the columns of the parity
7//! check matrix.
8
9use crate::{gf2::GF2, linalg, sparse::SparseMatrix};
10use ndarray::Array2;
11use num_traits::{One, Zero};
12use thiserror::Error;
13
14/// Systematic construction error.
15#[derive(Debug, Copy, Clone, Eq, PartialEq, Hash, Error)]
16pub enum Error {
17    /// The parity check matrix has more rows than columns.
18    #[error("the parity check matrix has more rows than columns")]
19    ParityOverdetermined,
20    /// The parity check matrix does not have full rank.
21    #[error("the parity check matrix does not have full rank")]
22    NotFullRank,
23}
24
25/// Permutes the columns of the parity check matrix to obtain a parity check
26/// matrix that supports systematic encoding using the first variables.
27///
28/// This function returns a parity check matrix obtaining by permuting the
29/// columns of `h` in such a way that the square submatrix formed by the last
30/// columns is invertible.
31pub fn parity_to_systematic(h: &SparseMatrix) -> Result<SparseMatrix, Error> {
32    let n = h.num_rows();
33    let m = h.num_cols();
34    if n > m {
35        return Err(Error::ParityOverdetermined);
36    }
37    let mut a = Array2::zeros((n, m));
38    for (j, k) in h.iter_all() {
39        a[[j, k]] = GF2::one();
40    }
41    linalg::row_echelon_form(&mut a);
42    // Check that the matrix has full rank by checking that there is a non-zero
43    // element in the last row (we start looking by the end, since chances are
44    // higher to find a non-zero element there).
45    if !(0..m).rev().any(|j| a[[n - 1, j]] != GF2::zero()) {
46        return Err(Error::NotFullRank);
47    }
48    // write point for columns that do not "go down" in the row echelon form
49    let mut k = 0;
50    let mut h_new = SparseMatrix::new(n, m);
51    let mut j0 = 0;
52    for j in 0..n {
53        assert!(k < m - n);
54        let mut found = false;
55        for s in j0..m {
56            if a[[j, s]] == GF2::zero() {
57                // Column does not "go down" on row echelon form. Place it at the current write point.
58                for &u in h.iter_col(s) {
59                    h_new.insert(u, k);
60                }
61                k += 1;
62            } else {
63                // Column goes down on row echelon form. Move to its appropriate
64                // position in the last columns.
65                let col = m - n + j;
66                for &u in h.iter_col(s) {
67                    h_new.insert(u, col);
68                }
69                found = true;
70                j0 = s + 1;
71                break;
72            }
73        }
74        assert!(found);
75    }
76    // Insert remaining columns at the write point
77    for j in j0..m {
78        assert!(k < m - n);
79        for &u in h.iter_col(j) {
80            h_new.insert(u, k);
81        }
82        k += 1;
83    }
84    Ok(h_new)
85}
86
87#[cfg(test)]
88mod test {
89    use super::*;
90
91    #[test]
92    fn to_systematic() {
93        let mut h = SparseMatrix::new(3, 9);
94        h.insert_col(0, [0, 1, 2].into_iter());
95        h.insert_col(1, [0, 2].into_iter());
96        // h.insert_col(2, [].into_iter()); this does nothing and does not compile
97        h.insert_col(3, [1].into_iter());
98        h.insert_col(4, [0, 1].into_iter());
99        h.insert_col(5, [1, 2].into_iter());
100        h.insert_col(6, [0, 2].into_iter());
101        h.insert_col(7, [1].into_iter());
102        h.insert_col(8, [0, 2].into_iter());
103        let mut expected = SparseMatrix::new(3, 9);
104        expected.insert_col(6, [0, 1, 2].into_iter());
105        expected.insert_col(7, [0, 2].into_iter());
106        // expected.insert_col(0, [].into_iter()); this does nothing and does not compile
107        expected.insert_col(1, [1].into_iter());
108        expected.insert_col(8, [0, 1].into_iter());
109        expected.insert_col(2, [1, 2].into_iter());
110        expected.insert_col(3, [0, 2].into_iter());
111        expected.insert_col(4, [1].into_iter());
112        expected.insert_col(5, [0, 2].into_iter());
113        assert_eq!(parity_to_systematic(&h).unwrap(), expected);
114    }
115}