ldpc_toolbox/
sparse.rs

1//! # Sparse binary matrix representation and functions
2//!
3//! This module implements a representation for sparse binary matrices based on
4//! the alist format used to handle LDPC parity check matrices.
5
6use std::borrow::Borrow;
7use std::slice::Iter;
8
9mod bfs;
10mod girth;
11
12pub use bfs::BFSResults;
13
14/// A [`String`] with an description of the error.
15pub type Error = String;
16/// A [`Result`] type containing an error [`String`].
17pub type Result<T> = std::result::Result<T, Error>;
18
19/// A sparse binary matrix
20///
21/// The internal representation for this matrix is based on the alist format.
22#[derive(PartialEq, Eq, Debug, Clone)]
23pub struct SparseMatrix {
24    rows: Vec<Vec<usize>>,
25    cols: Vec<Vec<usize>>,
26}
27
28impl SparseMatrix {
29    /// Create a new sparse matrix of a given size
30    ///
31    /// The matrix is inizialized to the zero matrix.
32    ///
33    /// # Examples
34    /// ```
35    /// # use ldpc_toolbox::sparse::SparseMatrix;
36    /// let h = SparseMatrix::new(10, 30);
37    /// assert_eq!(h.num_rows(), 10);
38    /// assert_eq!(h.num_cols(), 30);
39    /// ```
40    pub fn new(nrows: usize, ncols: usize) -> SparseMatrix {
41        use std::iter::repeat_with;
42        let rows = repeat_with(Vec::new).take(nrows).collect();
43        let cols = repeat_with(Vec::new).take(ncols).collect();
44        SparseMatrix { rows, cols }
45    }
46
47    /// Returns the number of rows of the matrix
48    pub fn num_rows(&self) -> usize {
49        self.rows.len()
50    }
51
52    /// Returns the number of columns of the matrix
53    pub fn num_cols(&self) -> usize {
54        self.cols.len()
55    }
56
57    /// Returns the row weight of `row`
58    ///
59    /// The row weight is defined as the number of entries equal to
60    /// one in a particular row. Rows are indexed starting by zero.
61    pub fn row_weight(&self, row: usize) -> usize {
62        self.rows[row].len()
63    }
64
65    /// Returns the column weight of `column`
66    ///
67    /// The column weight is defined as the number of entries equal to
68    /// one in a particular column. Columns are indexed starting by zero.
69    pub fn col_weight(&self, col: usize) -> usize {
70        self.cols[col].len()
71    }
72
73    /// Returns `true` if the entry corresponding to a particular
74    /// row and column is a one
75    pub fn contains(&self, row: usize, col: usize) -> bool {
76        // typically columns are shorter, so we search in the column
77        self.cols[col].contains(&row)
78    }
79
80    /// Inserts a one in a particular row and column.
81    ///
82    /// If there is already a one in this row and column, this function does
83    /// nothing.
84    ///
85    /// # Examples
86    /// ```
87    /// # use ldpc_toolbox::sparse::SparseMatrix;
88    /// let mut h = SparseMatrix::new(10, 30);
89    /// assert!(!h.contains(3, 7));
90    /// h.insert(3, 7);
91    /// assert!(h.contains(3, 7));
92    /// ```
93    pub fn insert(&mut self, row: usize, col: usize) {
94        if !self.contains(row, col) {
95            self.rows[row].push(col);
96            self.cols[col].push(row);
97        }
98    }
99
100    /// Removes a one in a particular row and column.
101    ///
102    /// If there is no one in this row and column, this function does nothing.
103    ///
104    /// # Examples
105    /// ```
106    /// # use ldpc_toolbox::sparse::SparseMatrix;
107    /// let mut h = SparseMatrix::new(10, 30);
108    /// h.insert(3, 7);
109    /// assert!(h.contains(3, 7));
110    /// h.remove(3, 7);
111    /// assert!(!h.contains(3, 7));
112    /// ```
113    pub fn remove(&mut self, row: usize, col: usize) {
114        self.rows[row].retain(|&c| c != col);
115        self.cols[col].retain(|&r| r != row);
116    }
117
118    /// Toggles the 0/1 in a particular row and column.
119    ///
120    /// If the row and column contains a zero, this function sets a one, and
121    /// vice versa. This is useful to implement addition modulo 2.
122    pub fn toggle(&mut self, row: usize, col: usize) {
123        match self.contains(row, col) {
124            true => self.remove(row, col),
125            false => self.insert(row, col),
126        }
127    }
128
129    /// Inserts ones in particular columns of a row
130    ///
131    /// This effect is as calling `insert()` on each of the elements
132    /// of the iterator `cols`.
133    ///
134    /// # Examples
135    /// ```
136    /// # use ldpc_toolbox::sparse::SparseMatrix;
137    /// let mut h1 = SparseMatrix::new(10, 30);
138    /// let mut h2 = SparseMatrix::new(10, 30);
139    /// let c = vec![3, 7, 9];
140    /// h1.insert_row(0, c.iter());
141    /// for a in &c {
142    ///     h2.insert(0, *a);
143    /// }
144    /// assert_eq!(h1, h2);
145    /// ```
146    pub fn insert_row<T, S>(&mut self, row: usize, cols: T)
147    where
148        T: Iterator<Item = S>,
149        S: Borrow<usize>,
150    {
151        for col in cols {
152            self.insert(row, *col.borrow());
153        }
154    }
155
156    /// Inserts ones in a particular rows of a column
157    ///
158    /// This works like `insert_row()`.
159    pub fn insert_col<T, S>(&mut self, col: usize, rows: T)
160    where
161        T: Iterator<Item = S>,
162        S: Borrow<usize>,
163    {
164        for row in rows {
165            self.insert(*row.borrow(), col);
166        }
167    }
168
169    /// Remove all the ones in a particular row
170    pub fn clear_row(&mut self, row: usize) {
171        for &col in &self.rows[row] {
172            self.cols[col].retain(|r| *r != row);
173        }
174        self.rows[row].clear();
175    }
176
177    /// Remove all the ones in a particular column
178    pub fn clear_col(&mut self, col: usize) {
179        for &row in &self.cols[col] {
180            self.rows[row].retain(|c| *c != col);
181        }
182        self.cols[col].clear();
183    }
184
185    /// Set the elements that are equal to one in a row
186    ///
187    /// The effect of this is like calling `clear_row()` followed
188    /// by `insert_row()`.
189    pub fn set_row<T, S>(&mut self, row: usize, cols: T)
190    where
191        T: Iterator<Item = S>,
192        S: Borrow<usize>,
193    {
194        self.clear_row(row);
195        self.insert_row(row, cols);
196    }
197
198    /// Set the elements that are equal to one in a column
199    pub fn set_col<T, S>(&mut self, col: usize, rows: T)
200    where
201        T: Iterator<Item = S>,
202        S: Borrow<usize>,
203    {
204        self.clear_col(col);
205        self.insert_col(col, rows);
206    }
207
208    /// Returns an [Iterator] over the indices entries equal to one in all the
209    /// matrix.
210    pub fn iter_all(&self) -> impl Iterator<Item = (usize, usize)> + '_ {
211        self.rows
212            .iter()
213            .enumerate()
214            .flat_map(|(j, r)| r.iter().map(move |&k| (j, k)))
215    }
216
217    /// Returns an [Iterator] over the entries equal to one
218    /// in a particular row
219    pub fn iter_row(&self, row: usize) -> Iter<'_, usize> {
220        self.rows[row].iter()
221    }
222
223    /// Returns an [Iterator] over the entries equal to one
224    /// in a particular column
225    pub fn iter_col(&self, col: usize) -> Iter<'_, usize> {
226        self.cols[col].iter()
227    }
228
229    fn write_alist_maybe_padding<W: std::fmt::Write>(
230        &self,
231        w: &mut W,
232        use_padding: bool,
233    ) -> std::fmt::Result {
234        writeln!(w, "{} {}", self.num_cols(), self.num_rows())?;
235        let directions = [&self.cols, &self.rows];
236        let mut direction_lengths = [0, 0];
237        for (dir, len) in directions.iter().zip(direction_lengths.iter_mut()) {
238            *len = dir.iter().map(|el| el.len()).max().unwrap_or(0);
239        }
240        writeln!(w, "{} {}", direction_lengths[0], direction_lengths[1])?;
241        for dir in directions.iter() {
242            let mut lengths = dir.iter().map(|el| el.len());
243            if let Some(len) = lengths.next() {
244                write!(w, "{}", len)?;
245            }
246            for len in lengths {
247                write!(w, " {}", len)?;
248            }
249            writeln!(w)?;
250        }
251        for (dir, &dirlen) in directions.iter().zip(direction_lengths.iter()) {
252            for el in *dir {
253                let mut v = el.clone();
254                v.sort_unstable();
255                let vlen = v.len();
256                let mut v = v.iter().map(|x| x + 1);
257                if let Some(x) = v.next() {
258                    write!(w, "{}", x)?;
259                }
260                for x in v {
261                    write!(w, " {}", x)?;
262                }
263                if use_padding {
264                    if vlen == 0 {
265                        write!(w, "0")?;
266                    }
267                    // .max(1) because we've added one padding element if vlen
268                    // was zero
269                    let num_padding = dirlen - vlen.max(1);
270                    for _ in 0..num_padding {
271                        write!(w, " 0")?;
272                    }
273                }
274                writeln!(w)?;
275            }
276        }
277        Ok(())
278    }
279
280    /// Writes the matrix in alist format to a writer.
281    ///
282    /// This function includes zeros as padding for irregular codes, as
283    /// originally defined by MacKay.
284    ///
285    /// # Errors
286    /// If a call to `write!()` returns an error, this function returns
287    /// such an error.
288    pub fn write_alist<W: std::fmt::Write>(&self, w: &mut W) -> std::fmt::Result {
289        self.write_alist_maybe_padding(w, true)
290    }
291
292    /// Writes the matrix in alist format to a writer.
293    ///
294    /// This function does not include zeros as padding for irregular codes.
295    ///
296    /// # Errors
297    /// If a call to `write!()` returns an error, this function returns
298    /// such an error.
299    pub fn write_alist_no_padding<W: std::fmt::Write>(&self, w: &mut W) -> std::fmt::Result {
300        self.write_alist_maybe_padding(w, false)
301    }
302
303    /// Returns a [`String`] with the alist representation of the matrix.
304    ///
305    /// This function includes zeros as padding for irregular codes, as
306    /// originally defined by MacKay.
307    pub fn alist(&self) -> String {
308        let mut s = String::new();
309        self.write_alist(&mut s).unwrap();
310        s
311    }
312
313    /// Returns a [`String`] with the alist representation of the matrix.
314    ///
315    /// This function does not include zeros as padding for irregular codes.
316    pub fn alist_no_padding(&self) -> String {
317        let mut s = String::new();
318        self.write_alist_no_padding(&mut s).unwrap();
319        s
320    }
321
322    /// Constructs and returns a sparse matrix from its alist representation.
323    ///
324    /// This function is able to read alists that use zeros for padding in the
325    /// case of an irregular code (as was defined originally by MacKay), as well
326    /// as alists that omit these zeros.
327    ///
328    /// # Errors
329    /// `alist` should hold a valid alist representation. If an error is found
330    /// while parsing `alist`, a `String` describing the error will be returned.
331    pub fn from_alist(alist: &str) -> Result<SparseMatrix> {
332        let mut alist = alist.split('\n');
333        let sizes = alist
334            .next()
335            .ok_or_else(|| String::from("alist first line not found"))?;
336        let mut sizes = sizes.split_whitespace();
337        let ncols = sizes
338            .next()
339            .ok_or_else(|| String::from("alist first line does not contain enough elements"))?
340            .parse()
341            .map_err(|_| String::from("ncols is not a number"))?;
342        let nrows = sizes
343            .next()
344            .ok_or_else(|| String::from("alist first line does not contain enough elements"))?
345            .parse()
346            .map_err(|_| String::from("nrows is not a number"))?;
347        let mut h = SparseMatrix::new(nrows, ncols);
348        alist.next(); // skip max weights
349        alist.next();
350        alist.next(); // skip weights
351        for col in 0..ncols {
352            let col_data = alist
353                .next()
354                .ok_or_else(|| String::from("alist does not contain expected number of lines"))?;
355            let col_data = col_data.split_whitespace();
356            for row in col_data {
357                let row: usize = row
358                    .parse()
359                    .map_err(|_| String::from("row value is not a number"))?;
360                // row == 0 is used for padding in irregular codes
361                if row != 0 {
362                    h.insert(row - 1, col);
363                }
364            }
365        }
366        // we do not need to process the rows of the alist
367        Ok(h)
368    }
369
370    /// Returns the girth of the bipartite graph defined by the matrix
371    ///
372    /// The girth is the length of the shortest cycle. If there are no
373    /// cycles, `None` is returned.
374    ///
375    /// # Examples
376    /// The following shows that a 2 x 2 matrix whose entries are all
377    /// equal to one has a girth of 4, which is the smallest girth that
378    /// a bipartite graph can have.
379    /// ```
380    /// # use ldpc_toolbox::sparse::SparseMatrix;
381    /// let mut h = SparseMatrix::new(2, 2);
382    /// for j in 0..2 {
383    ///     for k in 0..2 {
384    ///         h.insert(j, k);
385    ///     }
386    /// }
387    /// assert_eq!(h.girth(), Some(4));
388    /// ```
389    pub fn girth(&self) -> Option<usize> {
390        self.girth_with_max(usize::MAX)
391    }
392
393    /// Returns the girth of the bipartite graph defined by the matrix
394    /// as long as it is smaller than a maximum
395    ///
396    /// By imposing a maximum value in the girth search algorithm,
397    /// the execution time is reduced, since paths longer than the
398    /// maximum do not need to be explored.
399    ///
400    /// Often it is only necessary to check that a graph has at least
401    /// some minimum girth, so it is possible to use `girth_with_max()`.
402    ///
403    /// If there are no cycles with length smaller or equal to `max`, then
404    /// `None` is returned.
405    pub fn girth_with_max(&self, max: usize) -> Option<usize> {
406        (0..self.num_cols())
407            .filter_map(|c| self.girth_at_node_with_max(Node::Col(c), max))
408            .min()
409    }
410
411    /// Returns the local girth at a particular node
412    ///
413    /// The local girth at a node of a graph is defined as the minimum
414    /// length of the cycles containing that node.
415    ///
416    /// This function returns the local girth of at the node correponding
417    /// to a column or row of the matrix, or `None` if there are no cycles containing
418    /// that node.
419    pub fn girth_at_node(&self, node: Node) -> Option<usize> {
420        self.girth_at_node_with_max(node, usize::MAX)
421    }
422
423    /// Returns the girth at a particular node with a maximum
424    ///
425    /// This function works like `girth_at_node()` but imposes a maximum in the
426    /// length of the cycles considered. `None` is returned if there are no
427    /// cycles containing the node with length smaller or equal than `max`.
428    pub fn girth_at_node_with_max(&self, node: Node, max: usize) -> Option<usize> {
429        bfs::BFSContext::new(self, node).local_girth(max)
430    }
431
432    /// Run the BFS algorithm
433    ///
434    /// This uses a node of the graph associated to the matrix as the root
435    /// for the BFS algorithm and finds the distances from each of the nodes
436    /// of the graph to that root using breadth-first search.
437    /// # Examples
438    /// Run BFS on a matrix that has two connected components.
439    /// ```
440    /// # use ldpc_toolbox::sparse::SparseMatrix;
441    /// # use ldpc_toolbox::sparse::Node;
442    /// let mut h = SparseMatrix::new(4, 4);
443    /// for j in 0..4 {
444    ///     for k in 0..4 {
445    ///         if (j % 2) == (k % 2) {
446    ///             h.insert(j, k);
447    ///         }
448    ///     }
449    /// }
450    /// println!("{:?}", h.bfs(Node::Col(0)));
451    /// ```
452    pub fn bfs(&self, node: Node) -> BFSResults {
453        bfs::BFSContext::new(self, node).bfs()
454    }
455}
456
457/// A node in the graph associated to a sparse matrix
458///
459/// A node can represent a row or a column of the graph.
460#[derive(Debug, Copy, Clone, Eq, PartialEq)]
461pub enum Node {
462    /// Node representing row number `n`
463    Row(usize),
464    /// Node representing column number `n`
465    Col(usize),
466}
467
468impl Node {
469    fn iter(self, h: &SparseMatrix) -> impl Iterator<Item = Node> + '_ {
470        match self {
471            Node::Row(n) => h.iter_row(n),
472            Node::Col(n) => h.iter_col(n),
473        }
474        .map(move |&x| match self {
475            Node::Row(_) => Node::Col(x),
476            Node::Col(_) => Node::Row(x),
477        })
478    }
479}
480
481#[cfg(test)]
482mod tests {
483    use super::*;
484
485    #[test]
486    fn test_insert() {
487        let mut h = SparseMatrix::new(100, 300);
488        assert!(!h.contains(27, 154));
489        h.insert(27, 154);
490        assert!(h.contains(27, 154));
491        assert!(!h.contains(28, 154));
492    }
493
494    #[test]
495    fn test_insert_twice() {
496        let mut h = SparseMatrix::new(100, 300);
497        h.insert(27, 154);
498        h.insert(43, 28);
499        h.insert(53, 135);
500        let h2 = h.clone();
501        h.insert(43, 28);
502        assert_eq!(h, h2);
503    }
504
505    #[test]
506    fn iter_all() {
507        use std::collections::HashSet;
508
509        let mut h = SparseMatrix::new(10, 20);
510        let entries = [
511            (7, 8),
512            (5, 14),
513            (6, 6),
514            (6, 7),
515            (8, 10),
516            (0, 4),
517            (0, 0),
518            (0, 15),
519        ];
520        for entry in &entries {
521            h.insert(entry.0, entry.1);
522        }
523        let result = h.iter_all().collect::<HashSet<_>>();
524        assert_eq!(result, HashSet::from(entries));
525    }
526
527    #[test]
528    fn test_alist() {
529        let mut h = SparseMatrix::new(4, 12);
530        for j in 0..4 {
531            h.insert(j, j);
532            h.insert(j, j + 4);
533            h.insert(j, j + 8);
534        }
535        let expected = "12 4
5361 3
5371 1 1 1 1 1 1 1 1 1 1 1
5383 3 3 3
5391
5402
5413
5424
5431
5442
5453
5464
5471
5482
5493
5504
5511 5 9
5522 6 10
5533 7 11
5544 8 12
555";
556        assert_eq!(h.alist(), expected);
557
558        let h2 = SparseMatrix::from_alist(expected).unwrap();
559        assert_eq!(h2.alist(), expected);
560    }
561
562    #[test]
563    fn test_alist_irregular() {
564        let mut h = SparseMatrix::new(4, 12);
565        for j in 0..4 {
566            h.insert(j, j);
567            h.insert(j, j + 4);
568            if j < 2 {
569                h.insert(j, j + 8);
570            }
571        }
572
573        // with zero padding
574
575        let expected = "12 4
5761 3
5771 1 1 1 1 1 1 1 1 1 0 0
5783 3 2 2
5791
5802
5813
5824
5831
5842
5853
5864
5871
5882
5890
5900
5911 5 9
5922 6 10
5933 7 0
5944 8 0
595";
596        let expected_no_padding = "12 4
5971 3
5981 1 1 1 1 1 1 1 1 1 0 0
5993 3 2 2
6001
6012
6023
6034
6041
6052
6063
6074
6081
6092
610
611
6121 5 9
6132 6 10
6143 7
6154 8
616";
617
618        assert_eq!(h.alist(), expected);
619        assert_eq!(h.alist_no_padding(), expected_no_padding);
620        let h2 = SparseMatrix::from_alist(expected).unwrap();
621        assert_eq!(h2.alist(), expected);
622        assert_eq!(h2.alist_no_padding(), expected_no_padding);
623        let h3 = SparseMatrix::from_alist(expected_no_padding).unwrap();
624        assert_eq!(h3.alist(), expected);
625        assert_eq!(h3.alist_no_padding(), expected_no_padding);
626    }
627}