diff --git a/src/bin/bench_reach_symbolic.rs b/src/bin/bench_reach_symbolic.rs new file mode 100644 index 0000000..a15152d --- /dev/null +++ b/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.") +} diff --git a/src/bin/bench_reach_symbolic_basic.rs b/src/bin/bench_reach_symbolic_basic.rs new file mode 100644 index 0000000..6f01d96 --- /dev/null +++ b/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.") +} diff --git a/src/symbolic_async_graph/_impl_symbolic_async_graph_algorithm.rs b/src/symbolic_async_graph/_impl_symbolic_async_graph_algorithm.rs index 33825fe..8c010fb 100644 --- a/src/symbolic_async_graph/_impl_symbolic_async_graph_algorithm.rs +++ b/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] diff --git a/src/symbolic_async_graph/mod.rs b/src/symbolic_async_graph/mod.rs index a53130e..f00e621 100644 --- a/src/symbolic_async_graph/mod.rs +++ b/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. /// diff --git a/src/symbolic_async_graph/reachability.rs b/src/symbolic_async_graph/reachability.rs new file mode 100644 index 0000000..0be6a90 --- /dev/null +++ b/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::>(); + (regulator, targets) + }) + .collect::>(); + + 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; + } + } +}