sim_lib_discrete_algebra/
closure.rs

1//! Kleene closure `A* = I + A + A^2 + ...` over a closed semiring.
2//!
3//! This is the single derivation point for several graph algorithms:
4//! boolean closure yields transitive reachability, min-plus closure yields
5//! all-pairs shortest paths, and so on. The graph crate wraps this; it never
6//! re-implements Floyd-Warshall or Warshall directly.
7
8use crate::error::AlgebraError;
9use crate::matrix::{AlgebraLimits, Matrix};
10use crate::semiring::Semiring;
11
12impl<S: Semiring> Matrix<S> {
13    /// Compute the Kleene closure `A* = I + A + A^2 + ...`.
14    ///
15    /// Uses the generalized Floyd-Warshall (Lehmann) asteration: for each
16    /// pivot `k`, paths may pass through `k` and loop there `star(A[k][k])`
17    /// times. The identity is added at the end so `A*` includes the empty path
18    /// on the diagonal (distance-0 for min-plus, reflexive for boolean).
19    ///
20    /// Returns [`AlgebraError::NoStar`] when a pivot's diagonal star does not
21    /// converge (e.g. a negative cycle in min-plus, any directed cycle in
22    /// counting, or a semiring with no `star` such as `RealF64`).
23    ///
24    /// # Examples
25    ///
26    /// Boolean closure of a directed chain `0 -> 1 -> 2` yields reflexive
27    /// reachability: node `i` reaches node `j` exactly when `j >= i`.
28    ///
29    /// ```
30    /// use sim_lib_discrete_algebra::{AlgebraLimits, BoolRing, Matrix};
31    ///
32    /// let mut a = Matrix::new(3, 3);
33    /// a.set(0, 1, BoolRing(true)).unwrap();
34    /// a.set(1, 2, BoolRing(true)).unwrap();
35    ///
36    /// let reach = a.closure(AlgebraLimits::default()).unwrap();
37    /// assert_eq!(reach.get(0, 2).unwrap(), &BoolRing(true)); // 0 reaches 2
38    /// assert_eq!(reach.get(1, 1).unwrap(), &BoolRing(true)); // reflexive
39    /// assert_eq!(reach.get(2, 0).unwrap(), &BoolRing(false)); // 2 cannot reach 0
40    /// ```
41    pub fn closure(&self, limits: AlgebraLimits) -> Result<Self, AlgebraError> {
42        self.validate()?;
43        if !self.is_square() {
44            return Err(AlgebraError::ShapeMismatch(format!(
45                "closure requires a square matrix, got {}x{}",
46                self.rows, self.cols
47            )));
48        }
49        let n = self.rows;
50        if n > limits.max_dim {
51            return Err(AlgebraError::LimitExceeded(format!(
52                "closure: dimension {n} exceeds max_dim {}",
53                limits.max_dim
54            )));
55        }
56        let mut c = self.clone();
57        for k in 0..n {
58            let s = c.data[k * n + k].star().ok_or(AlgebraError::NoStar)?;
59            let mut next = c.clone();
60            for i in 0..n {
61                let cik = c.data[i * n + k].mul(&s);
62                for j in 0..n {
63                    let via = cik.mul(&c.data[k * n + j]);
64                    next.data[i * n + j] = c.data[i * n + j].add(&via);
65                }
66            }
67            c = next;
68        }
69        // Add the identity so the closure contains the empty path A^0 = I.
70        for i in 0..n {
71            let updated = c.data[i * n + i].add(&S::one());
72            c.data[i * n + i] = updated;
73        }
74        Ok(c)
75    }
76}
77
78#[cfg(test)]
79mod tests {
80    use super::*;
81    use crate::boolean::BoolRing;
82    use crate::tropical_min::MinPlus;
83
84    #[test]
85    fn boolean_closure_of_directed_chain_is_upper_triangular_reachability() {
86        // Directed chain 0->1->2->3. Reachability is the reflexive upper
87        // triangle: i reaches j iff j >= i.
88        let mut a = Matrix::new(4, 4);
89        a.set(0, 1, BoolRing(true)).unwrap();
90        a.set(1, 2, BoolRing(true)).unwrap();
91        a.set(2, 3, BoolRing(true)).unwrap();
92        let star = a.closure(AlgebraLimits::default()).unwrap();
93        for i in 0..4 {
94            for j in 0..4 {
95                let expected = BoolRing(j >= i);
96                assert_eq!(star.get(i, j).unwrap(), &expected, "reachability ({i},{j})");
97            }
98        }
99    }
100
101    #[test]
102    fn min_plus_closure_matches_floyd_warshall() {
103        // Directed weighted graph: 0->1 (1), 1->2 (2), 0->2 (5), 2->3 (1).
104        let mut a = Matrix::filled(4, 4, MinPlus::Inf);
105        a.set(0, 1, MinPlus::Fin(1)).unwrap();
106        a.set(1, 2, MinPlus::Fin(2)).unwrap();
107        a.set(0, 2, MinPlus::Fin(5)).unwrap();
108        a.set(2, 3, MinPlus::Fin(1)).unwrap();
109        let star = a.closure(AlgebraLimits::default()).unwrap();
110        // shortest 0->2 is via 1: 1+2=3, beating the direct edge 5.
111        assert_eq!(star.get(0, 2).unwrap(), &MinPlus::Fin(3));
112        // 0->3 = 1+2+1 = 4.
113        assert_eq!(star.get(0, 3).unwrap(), &MinPlus::Fin(4));
114        // diagonal is 0 (empty path); unreachable stays Inf.
115        assert_eq!(star.get(0, 0).unwrap(), &MinPlus::Fin(0));
116        assert_eq!(star.get(3, 0).unwrap(), &MinPlus::Inf);
117    }
118
119    #[test]
120    fn closure_without_star_reports_no_star() {
121        use crate::real::RealF64;
122        let a: Matrix<RealF64> = Matrix::filled(2, 2, RealF64(0.5));
123        assert_eq!(
124            a.closure(AlgebraLimits::default()).unwrap_err(),
125            AlgebraError::NoStar
126        );
127    }
128}
sim_lib_discrete_algebra/closure.rs

sim_lib_discrete_algebra/
closure.rs
