Reachability method with less overhead (#55)

* Proposal for new, faster reachability procedures.

* Adapt to `0.5.0`.

* Nicer documentation.

* Remove old reachability methods.

* Better reachability test.

src/bin/bench_reach_symbolic.rs

@@ -0,0 +1,39 @@
+use biodivine_lib_param_bn::biodivine_std::traits::Set;
+use biodivine_lib_param_bn::symbolic_async_graph::reachability::Reachability;
+use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph;
+use biodivine_lib_param_bn::BooleanNetwork;
+
+fn main() {
+    let args = Vec::from_iter(std::env::args());
+    let path = &args[1];
+    let model = BooleanNetwork::try_from_file(path).unwrap();
+    let model = model.inline_inputs(true, true);
+
+    println!(
+        "Loaded model with {} variables and {} parameters.",
+        model.num_vars(),
+        model.num_parameters()
+    );
+
+    let stg = SymbolicAsyncGraph::new(&model).unwrap();
+
+    let mut count = 0;
+    let mut universe = stg.mk_unit_colored_vertices();
+    while !universe.is_empty() {
+        let reduced_stg = stg.restrict(&universe);
+        let mut component = universe.pick_vertex();
+        loop {
+            let update = Reachability::reach_bwd(&reduced_stg, &component);
+            let update = Reachability::reach_fwd(&reduced_stg, &update);
+            if update.is_subset(&component) {
+                break;
+            } else {
+                component = update;
+            }
+        }
+        println!("Found weak component: {}", component.approx_cardinality());
+        universe = universe.minus(&component);
+        count += 1;
+    }
+    println!("Total {count} components.")
+}

src/bin/bench_reach_symbolic_basic.rs

@@ -0,0 +1,39 @@
+use biodivine_lib_param_bn::biodivine_std::traits::Set;
+use biodivine_lib_param_bn::symbolic_async_graph::reachability::Reachability;
+use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph;
+use biodivine_lib_param_bn::BooleanNetwork;
+
+fn main() {
+    let args = Vec::from_iter(std::env::args());
+    let path = &args[1];
+    let model = BooleanNetwork::try_from_file(path).unwrap();
+    let model = model.inline_inputs(true, true);
+
+    println!(
+        "Loaded model with {} variables and {} parameters.",
+        model.num_vars(),
+        model.num_parameters()
+    );
+
+    let stg = SymbolicAsyncGraph::new(&model).unwrap();
+
+    let mut count = 0;
+    let mut universe = stg.mk_unit_colored_vertices();
+    while !universe.is_empty() {
+        let reduced_stg = stg.restrict(&universe);
+        let mut component = universe.pick_vertex();
+        loop {
+            let update = Reachability::reach_bwd_basic(&reduced_stg, &component);
+            let update = Reachability::reach_fwd_basic(&reduced_stg, &update);
+            if update.is_subset(&component) {
+                break;
+            } else {
+                component = update;
+            }
+        }
+        println!("Found weak component: {}", component.approx_cardinality());
+        universe = universe.minus(&component);
+        count += 1;
+    }
+    println!("Total {count} components.")
+}

src/symbolic_async_graph/_impl_symbolic_async_graph_algorithm.rs

@@ -1,4 +1,5 @@
 use crate::biodivine_std::traits::Set;
+use crate::symbolic_async_graph::reachability::Reachability;
 use crate::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph};
 use crate::{ExtendedBoolean, Space, VariableId};
 
@@ -13,36 +14,14 @@ impl SymbolicAsyncGraph {
     ///
     /// In other words, the result is the smallest forward-closed superset of `initial`.
     pub fn reach_forward(&self, initial: &GraphColoredVertices) -> GraphColoredVertices {
-        let mut result = initial.clone();
-        'fwd: loop {
-            for var in self.variables().rev() {
-                let step = self.var_post_out(var, &result);
-                if !step.is_empty() {
-                    result = result.union(&step);
-                    continue 'fwd;
-                }
-            }
-
-            return result;
-        }
+        Reachability::reach_fwd(self, initial)
     }
 
     /// Compute the set of backward-reachable vertices from the given `initial` set.
     ///
     /// In other words, the result is the smallest backward-closed superset of `initial`.
     pub fn reach_backward(&self, initial: &GraphColoredVertices) -> GraphColoredVertices {
-        let mut result = initial.clone();
-        'bwd: loop {
-            for var in self.variables().rev() {
-                let step = self.var_pre_out(var, &result);
-                if !step.is_empty() {
-                    result = result.union(&step);
-                    continue 'bwd;
-                }
-            }
-
-            return result;
-        }
+        Reachability::reach_bwd(self, initial)
     }
 
     /// Compute the subset of `initial` vertices that can only reach other vertices within
@@ -139,6 +118,7 @@ impl SymbolicAsyncGraph {
 #[cfg(test)]
 mod tests {
     use crate::biodivine_std::traits::Set;
+    use crate::symbolic_async_graph::reachability::Reachability;
     use crate::symbolic_async_graph::SymbolicAsyncGraph;
     use crate::ExtendedBoolean::Zero;
     use crate::{BooleanNetwork, ExtendedBoolean, Space};
@@ -165,6 +145,8 @@ mod tests {
         let pivot = stg.unit_colored_vertices().pick_vertex();
         let fwd = stg.reach_forward(&pivot);
         let bwd = stg.reach_backward(&pivot);
+        let fwd_basic = Reachability::reach_fwd_basic(&stg, &pivot);
+        let bwd_basic = Reachability::reach_bwd_basic(&stg, &pivot);
         let scc = fwd.intersect(&bwd);
 
         // Should contain all cases where pivot is in an attractor.
@@ -179,6 +161,13 @@ mod tests {
         // For some reason, there is only a very small number of parameter valuations
         // where this is not empty -- less than in the pivot in fact.
         assert!(!repeller_basin.is_empty());
+
+        assert!(pivot.is_subset(&fwd));
+        assert!(pivot.is_subset(&bwd));
+        assert!(pivot.is_subset(&scc));
+
+        assert_eq!(fwd, fwd_basic);
+        assert_eq!(bwd, bwd_basic);
     }
 
     #[test]

src/symbolic_async_graph/mod.rs

@@ -58,6 +58,8 @@ mod _impl_symbolic_context;
 /// that are used to iterate through various projections of symbolic sets.
 pub mod projected_iteration;
 
+pub mod reachability;
+
 /// A module with a trait that describes common methods shared by all set representations
 /// based on BDDs.
 ///

src/symbolic_async_graph/reachability.rs

@@ -0,0 +1,184 @@
+use crate::biodivine_std::traits::Set;
+use crate::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph};
+use crate::VariableId;
+use std::cmp::max;
+use std::collections::{HashMap, HashSet};
+
+pub struct Reachability {
+    _dummy: (),
+}
+
+impl Reachability {
+    /// A FWD reachability procedure that uses "structural saturation".
+    pub fn reach_fwd(
+        graph: &SymbolicAsyncGraph,
+        initial: &GraphColoredVertices,
+    ) -> GraphColoredVertices {
+        Self::reach(graph, initial, |g, s, v| g.var_post_out(v, s))
+    }
+
+    /// A BWD reachability procedure that uses "structural saturation".
+    pub fn reach_bwd(
+        graph: &SymbolicAsyncGraph,
+        initial: &GraphColoredVertices,
+    ) -> GraphColoredVertices {
+        Self::reach(graph, initial, |g, s, v| g.var_pre_out(v, s))
+    }
+
+    /// "Basic" saturation FWD reachability procedure.
+    pub fn reach_fwd_basic(
+        graph: &SymbolicAsyncGraph,
+        initial: &GraphColoredVertices,
+    ) -> GraphColoredVertices {
+        Self::reach_basic_saturation(graph, initial, |g, s, v| g.var_post_out(v, s))
+    }
+
+    /// "Basic" saturation BWD reachability procedure.
+    pub fn reach_bwd_basic(
+        graph: &SymbolicAsyncGraph,
+        initial: &GraphColoredVertices,
+    ) -> GraphColoredVertices {
+        Self::reach_basic_saturation(graph, initial, |g, s, v| g.var_pre_out(v, s))
+    }
+
+    /// A basic saturation procedure which performs reachability over all transition functions of a graph.
+    pub fn reach_basic_saturation<
+        F: Fn(&SymbolicAsyncGraph, &GraphColoredVertices, VariableId) -> GraphColoredVertices,
+    >(
+        graph: &SymbolicAsyncGraph,
+        initial: &GraphColoredVertices,
+        step: F,
+    ) -> GraphColoredVertices {
+        if cfg!(feature = "print-progress") {
+            println!(
+                "Start basic symbolic reachability with {}[nodes:{}] candidates.",
+                initial.approx_cardinality(),
+                initial.symbolic_size()
+            );
+        }
+
+        let mut result = initial.clone();
+
+        let mut max_size = 0;
+
+        'reach: loop {
+            for var in graph.variables().rev() {
+                let step = step(graph, &result, var);
+                if !step.is_empty() {
+                    result = result.union(&step);
+                    max_size = max(max_size, result.symbolic_size());
+
+                    if cfg!(feature = "print-progress") {
+                        let current = result.approx_cardinality();
+                        let max = graph.unit_colored_vertices().approx_cardinality();
+                        println!(
+                            " >> Reach progress: {}[nodes:{}] candidates ({:.2} log-%).",
+                            result.approx_cardinality(),
+                            result.symbolic_size(),
+                            (current.log2() / max.log2()) * 100.0
+                        );
+                    }
+
+                    continue 'reach;
+                }
+            }
+
+            if cfg!(feature = "print-progress") {
+                println!(
+                    "Basic reachability done: {}[nodes:{}] candidates. Max intermediate size: {}.",
+                    result.approx_cardinality(),
+                    result.symbolic_size(),
+                    max_size
+                );
+            }
+
+            return result;
+        }
+    }
+
+    /// A "structural" reachability procedure which uses the dependencies between the update functions to reduce
+    /// the number of transitions that need to be tested.
+    ///
+    /// This sometimes increases the size of the BDDs a little bit, but is overall beneficial in the vast majority
+    /// of cases.
+    pub fn reach<
+        F: Fn(&SymbolicAsyncGraph, &GraphColoredVertices, VariableId) -> GraphColoredVertices,
+    >(
+        graph: &SymbolicAsyncGraph,
+        initial: &GraphColoredVertices,
+        step: F,
+    ) -> GraphColoredVertices {
+        if cfg!(feature = "print-progress") {
+            println!(
+                "Start symbolic reachability with {}[nodes:{}] candidates.",
+                initial.approx_cardinality(),
+                initial.symbolic_size()
+            );
+        }
+
+        // Construct a symbolic regulators relation (this is more accurate than the normal regulatory graph,
+        // and does not require an underlying Boolean network).
+        let targets_of_var = graph
+            .variables()
+            .map(|regulator| {
+                let targets = graph
+                    .variables()
+                    .filter(|var| {
+                        let update = graph.get_symbolic_fn_update(*var);
+                        let state_variable = graph.symbolic_context().get_state_variable(regulator);
+                        update.support_set().contains(&state_variable)
+                    })
+                    .collect::<Vec<_>>();
+                (regulator, targets)
+            })
+            .collect::<HashMap<_, _>>();
+
+        let mut result = initial.clone();
+        let mut max_size = 0;
+
+        // A set of saturated variables. We can only remove a variable from here if it's regulator has been changed.
+        let mut saturated = HashSet::new();
+        'reach: loop {
+            for var in graph.variables().rev() {
+                let next_step = step(graph, &result, var);
+                if next_step.is_empty() {
+                    saturated.insert(var);
+                } else {
+                    result = result.union(&next_step);
+                    max_size = max(max_size, result.symbolic_size());
+
+                    if cfg!(feature = "print-progress") {
+                        let current = result.approx_cardinality();
+                        let max = graph.unit_colored_vertices().approx_cardinality();
+                        println!(
+                            " >> [{}/{} saturated] Reach progress: {}[nodes:{}] candidates ({:.2} log-%).",
+                            saturated.len(),
+                            graph.num_vars(),
+                            result.approx_cardinality(),
+                            result.symbolic_size(),
+                            (current.log2() / max.log2()) * 100.0
+                        );
+                    }
+
+                    if !saturated.contains(&var) {
+                        for t in &targets_of_var[&var] {
+                            saturated.remove(t);
+                        }
+                        continue 'reach;
+                    }
+                }
+            }
+
+            if cfg!(feature = "print-progress") {
+                println!(
+                    "Structural reachability done: {}[nodes:{}] candidates. Max intermediate size: {}.",
+                    result.approx_cardinality(),
+                    result.symbolic_size(),
+                    max_size
+                );
+            }
+
+            return result;
+        }
+    }
+}