add transformation to systematic code

This adds a CLI command that transforms a parity check matrix to
one that can be used for systematic encoding with the first variables
by permuting columns of the parity check matrix.

The crate::encoder::gauss module has been moved to crate::linalg
because it is now used outside crate::encoder too.

src/cli.rs

@@ -15,6 +15,7 @@ pub mod dvbs2;
 pub mod encode;
 pub mod mackay_neal;
 pub mod peg;
+pub mod systematic;
 
 /// Trait to run a CLI subcommand
 pub trait Run {
@@ -41,6 +42,8 @@ pub enum Args {
     MackayNeal(mackay_neal::Args),
     /// peg subcommand
     PEG(peg::Args),
+    /// systematic subcommand
+    Systematic(systematic::Args),
 }
 
 impl Run for Args {
@@ -53,6 +56,7 @@ impl Run for Args {
             Args::Encode(x) => x.run(),
             Args::MackayNeal(x) => x.run(),
             Args::PEG(x) => x.run(),
+            Args::Systematic(x) => x.run(),
         }
     }
 }

src/cli/systematic.rs

@@ -0,0 +1,27 @@
+//! Systematic CLI subcommand.
+//!
+//! This command can be used to convert an n x m parity check matrix into one
+//! that supports systematic encoding using the first m - n columns by permuting
+//! columns in such a way that the n x n submatrix formed by the last n columns
+//! is invertible.
+
+use crate::{cli::Run, sparse::SparseMatrix, systematic::parity_to_systematic};
+use clap::Parser;
+use std::error::Error;
+
+/// Systematic CLI arguments.
+#[derive(Debug, Parser)]
+#[command(about = "Converts a parity check matrix into systematic form")]
+pub struct Args {
+    /// alist file for the code
+    alist: String,
+}
+
+impl Run for Args {
+    fn run(&self) -> Result<(), Box<dyn Error>> {
+        let h = SparseMatrix::from_alist(&std::fs::read_to_string(&self.alist)?)?;
+        let h_sys = parity_to_systematic(&h)?;
+        println!("{}", h_sys.alist());
+        Ok(())
+    }
+}

src/encoder.rs

@@ -23,12 +23,11 @@
 //! right) to obtain the n-k parity check bits. In this case, the encoding
 //! complexity is O(n^2).
 
-use crate::{gf2::GF2, sparse::SparseMatrix};
+use crate::{gf2::GF2, linalg, sparse::SparseMatrix};
 use ndarray::{s, Array1, Array2, ArrayBase, Data, Ix1};
 use num_traits::One;
 use thiserror::Error;
 
-mod gauss;
 mod staircase;
 
 /// LDPC encoder error.
@@ -85,9 +84,9 @@ impl Encoder {
                 a[[j, t]] = GF2::one();
             }
 
-            match gauss::gauss_reduction(&mut a) {
+            match linalg::gauss_reduction(&mut a) {
                 Ok(()) => (),
-                Err(gauss::Error::NotInvertible) => return Err(Error::SubmatrixNotInvertible),
+                Err(linalg::Error::NotInvertible) => return Err(Error::SubmatrixNotInvertible),
             };
 
             let gen_matrix = a.slice(s![.., n..]).to_owned();

src/lib.rs

@@ -21,5 +21,7 @@ pub mod peg;
 pub mod rand;
 pub mod simulation;
 pub mod sparse;
+pub mod systematic;
 
+mod linalg;
 mod util;

src/encoder/gauss.rs → src/linalg.rs

@@ -64,6 +64,48 @@ pub fn gauss_reduction<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) -> Re
     Ok(())
 }
 
+pub fn row_echelon_form<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) {
+    let (n, m) = array.dim();
+
+    // Reduce to upper triangular with ones on diagonal
+    let mut k = 0;
+    for j in 0..n {
+        // Find non-zero element in current column, at or below row k
+        let Some(s) = array
+            .slice(s![k.., j])
+            .iter()
+            .enumerate()
+            .find_map(|(t, x)| if x.is_zero() { None } else { Some(k + t) })
+        else {
+            // All the elements at or below row k are zero. Done with this
+            // column.
+            continue;
+        };
+
+        if s != k {
+            // Swap rows s and k
+            for t in j..m {
+                array.swap([s, t], [k, t]);
+            }
+        }
+
+        let x = array[[k, j]];
+
+        // Subtract to rows below to make zeros below row k
+        for t in (k + 1)..n {
+            let y = array[[t, j]];
+            if !y.is_zero() {
+                // avoid calculations if we're subtracting zero
+                for u in j..m {
+                    array[[t, u]] = array[[t, u]] - y * array[[k, u]] / x;
+                }
+            }
+        }
+
+        k += 1;
+    }
+}
+
 #[cfg(test)]
 mod test {
     use super::*;
@@ -88,4 +130,22 @@ mod test {
         ]);
         assert_eq!(&a, &expected);
     }
+
+    #[test]
+    fn row_echelon() {
+        let i = GF2::one();
+        let o = GF2::zero();
+        let mut a = arr2(&[
+            [i, i, o, o, i, o, i, o, i],
+            [i, o, o, i, i, i, o, i, o],
+            [i, i, o, o, o, i, i, o, i],
+        ]);
+        row_echelon_form(&mut a);
+        let expected = arr2(&[
+            [i, i, o, o, i, o, i, o, i],
+            [o, i, o, i, o, i, i, i, i],
+            [o, o, o, o, i, i, o, o, o],
+        ]);
+        assert_eq!(&a, &expected);
+    }
 }

src/systematic.rs

@@ -0,0 +1,112 @@
+//! Systematic code constructions.
+//!
+//! This module contains a function [`parity_to_systematic`] that can be used to
+//! convert a full-rank parity check matrix into one that supports systematic
+//! encoding using the first variables (as done by the systematic encoder in the
+//! [`encoder`](crate::encoder) module) by permuting the columns of the parity
+//! check matrix.
+
+use crate::{gf2::GF2, linalg, sparse::SparseMatrix};
+use ndarray::Array2;
+use num_traits::{One, Zero};
+use thiserror::Error;
+
+/// Systematic construction error.
+#[derive(Debug, Copy, Clone, Eq, PartialEq, Hash, Error)]
+pub enum Error {
+    /// The parity check matrix has more rows than columns.
+    #[error("the parity check matrix has more rows than columns")]
+    ParityOverdetermined,
+    /// The parity check matrix does not have full rank.
+    #[error("the parity check matrix does not have full rank")]
+    NotFullRank,
+}
+
+/// Permutes the columns of the parity check matrix to obtain a parity check
+/// matrix that supports systematic encoding using the first variables.
+///
+/// This function returns a parity check matrix obtaining by permuting the
+/// columns of `h` in such a way that the square submatrix formed by the last
+/// columns is invertible.
+pub fn parity_to_systematic(h: &SparseMatrix) -> Result<SparseMatrix, Error> {
+    let n = h.num_rows();
+    let m = h.num_cols();
+    if n > m {
+        return Err(Error::ParityOverdetermined);
+    }
+    let mut a = Array2::zeros((n, m));
+    for (j, k) in h.iter_all() {
+        a[[j, k]] = GF2::one();
+    }
+    linalg::row_echelon_form(&mut a);
+    // write point for columns that do not "go down" in the row echelon form
+    let mut k = 0;
+    let mut h_new = SparseMatrix::new(n, m);
+    let mut j0 = 0;
+    for j in 0..n {
+        assert!(k < m - n);
+        let mut found = false;
+        for s in j0..m {
+            if a[[j, s]] == GF2::zero() {
+                // Column does not "go down" on row echelon form. Place it at the current write point.
+                for &u in h.iter_col(s) {
+                    h_new.insert(u, k);
+                }
+                k += 1;
+            } else {
+                // Column goes down on row echelon form. Move to its appropriate
+                // position in the last columns.
+                let col = m - n + j;
+                for &u in h.iter_col(s) {
+                    h_new.insert(u, col);
+                }
+                found = true;
+                j0 = s + 1;
+                break;
+            }
+        }
+        if !found {
+            // No column "going down" found. The matrix doesn't have full rank.
+            return Err(Error::NotFullRank);
+        }
+    }
+    // Insert remaining columns at the write point
+    for j in j0..m {
+        assert!(k < m - n);
+        for &u in h.iter_col(j) {
+            h_new.insert(u, k);
+        }
+        k += 1;
+    }
+    Ok(h_new)
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    #[test]
+    fn to_systematic() {
+        let mut h = SparseMatrix::new(3, 9);
+        h.insert_col(0, [0, 1, 2].into_iter());
+        h.insert_col(1, [0, 2].into_iter());
+        // h.insert_col(2, [].into_iter()); this does nothing and does not compile
+        h.insert_col(3, [1].into_iter());
+        h.insert_col(4, [0, 1].into_iter());
+        h.insert_col(5, [1, 2].into_iter());
+        h.insert_col(6, [0, 2].into_iter());
+        h.insert_col(7, [1].into_iter());
+        h.insert_col(8, [0, 2].into_iter());
+        let mut expected = SparseMatrix::new(3, 9);
+        expected.insert_col(6, [0, 1, 2].into_iter());
+        expected.insert_col(7, [0, 2].into_iter());
+        // expected.insert_col(0, [].into_iter()); this does nothing and does not compile
+        expected.insert_col(1, [1].into_iter());
+        expected.insert_col(8, [0, 1].into_iter());
+        expected.insert_col(2, [1, 2].into_iter());
+        expected.insert_col(3, [0, 2].into_iter());
+        expected.insert_col(4, [1].into_iter());
+        expected.insert_col(5, [0, 2].into_iter());
+        assert_eq!(parity_to_systematic(&h).unwrap(), expected);
+    }
+}