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}