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_function_table.rs b/src/symbolic_async_graph/_impl_function_table.rs index 5b87a01..883b868 100644 --- a/src/symbolic_async_graph/_impl_function_table.rs +++ b/src/symbolic_async_graph/_impl_function_table.rs @@ -14,7 +14,8 @@ impl FunctionTable { bdd_builder.make_variable(bdd_var_name.as_str()) }) .collect(); - FunctionTable { arity, rows } + let name = name.to_string(); + FunctionTable { arity, rows, name } } /// List the symbolic variables that appear in this function table. @@ -25,6 +26,11 @@ impl FunctionTable { pub fn symbolic_variables(&self) -> &Vec { &self.rows } + + /// True if this [FunctionTable] contains the provided [BddVariable]. + pub fn contains(&self, var: BddVariable) -> bool { + self.rows.contains(&var) + } } /// Converts a `FunctionTable` into an iterator of `Vec` (function table row) and 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/_impl_symbolic_context.rs b/src/symbolic_async_graph/_impl_symbolic_context.rs index 7e668e6..2bb6e31 100644 --- a/src/symbolic_async_graph/_impl_symbolic_context.rs +++ b/src/symbolic_async_graph/_impl_symbolic_context.rs @@ -6,6 +6,7 @@ use biodivine_lib_bdd::{ }; use std::collections::HashMap; use std::convert::TryInto; +use std::fmt::{Debug, Formatter}; impl SymbolicContext { /// Create a new `SymbolicContext` that encodes the given `BooleanNetwork`, but otherwise has @@ -114,14 +115,10 @@ impl SymbolicContext { // Finally, collect all parameter BddVariables into one vector. let mut parameter_variables: Vec = Vec::new(); for table in &explicit_function_tables { - for p in &table.rows { - parameter_variables.push(*p); - } + parameter_variables.extend_from_slice(&table.rows); } for table in implicit_function_tables.iter().flatten() { - for p in &table.rows { - parameter_variables.push(*p); - } + parameter_variables.extend_from_slice(&table.rows); } Ok(SymbolicContext { @@ -175,6 +172,92 @@ impl SymbolicContext { result } + /// Create a copy of this [SymbolicContext] which only contains the relevant state and parameter variables. + /// + /// In other words, the new context is equivalent to calling [SymbolicContext::new] with the original network + /// that was also used to create this context. + /// + /// In general, if you are using extra variables in your algorithm, it is preferred to export the result in + /// such a way that the BDDs are compatible with this "canonical" context whenever possible. To do that, you + /// can use this method to create the canonical context, and then use [SymbolicContext::transfer_from] + /// to translate the output BDDs. + pub fn as_canonical_context(&self) -> SymbolicContext { + let mut builder = BddVariableSetBuilder::new(); + let mut state_variables = Vec::new(); + let mut explicit_function_tables = Vec::new(); + let mut implicit_function_tables: Vec> = + vec![None; self.num_state_variables()]; + // Now, a key problem is that we need to recreate all objects in the correct order, such that we are + // truly compatible with the original representation. This is a bit complicated, because we don't have that + // order saved anywhere explicitly. However, we can "reconstruct" it from the current data-structures. + let mut to_skip = 0usize; // Variables are skipped when they represent a single uninterpreted function. + let mut next_var_id = 0usize; // We use this to tell which variable are we processing now/have processed last. + 'var_loop: for var in self.bdd.variables() { + if to_skip > 0 { + to_skip -= 1; + continue 'var_loop; + } + if self.state_variables.contains(&var) { + // Recreate state variable. + let variable_name = self.bdd.name_of(var); + let state_variable = builder.make_variable(variable_name.as_str()); + state_variables.push(state_variable); + next_var_id += 1; + continue 'var_loop; + } + if self.all_extra_state_variables.contains(&var) { + // Extra state variables are skipped. + continue 'var_loop; + } + + // Find the variable in the uninterpreted functions. + for function in &self.explicit_function_tables { + if function.contains(var) { + // Recreate function + let new_table = + FunctionTable::new(function.name.as_str(), function.arity, &mut builder); + to_skip = new_table.rows.len() - 1; + explicit_function_tables.push(new_table); + continue 'var_loop; + } + } + for function in &self.implicit_function_tables { + if let Some(function) = function.as_ref() { + if function.contains(var) { + let new_table = FunctionTable::new( + function.name.as_str(), + function.arity, + &mut builder, + ); + to_skip = new_table.rows.len() - 1; + assert!(next_var_id > 0); + implicit_function_tables[next_var_id - 1] = Some(new_table); + continue 'var_loop; + } + } + } + unreachable!("There shouldn't be any other symbolic variables."); + } + + let mut parameter_variables: Vec = Vec::new(); + for table in &explicit_function_tables { + parameter_variables.extend_from_slice(&table.rows); + } + for table in implicit_function_tables.iter().flatten() { + parameter_variables.extend_from_slice(&table.rows); + } + + SymbolicContext { + bdd: builder.build(), + state_variables, + extra_state_variables: vec![vec![]; self.num_state_variables()], + all_extra_state_variables: Vec::new(), + parameter_variables, + explicit_function_tables, + implicit_function_tables, + } + } + /// The number of state variables (should be equal to the number of network variables). pub fn num_state_variables(&self) -> usize { self.state_variables.len() @@ -496,11 +579,58 @@ fn network_symbolic_size(network: &BooleanNetwork, extra: &HashMap) -> std::fmt::Result { + write!( + f, + "SymbolicContext(states={}, extras={}, params={})", + self.num_state_variables(), + self.all_extra_state_variables.len(), + self.parameter_variables.len() + ) + } +} + +impl PartialEq for SymbolicContext { + fn eq(&self, other: &Self) -> bool { + if self.bdd.variables() != other.bdd.variables() { + return false; + } + + for var in self.bdd.variables() { + if self.bdd.name_of(var) != other.bdd.name_of(var) { + return false; + } + } + + if self.state_variables != other.state_variables { + return false; + } + + if self.extra_state_variables != other.extra_state_variables { + return false; + } + + if self.explicit_function_tables != other.explicit_function_tables { + return false; + } + + if self.implicit_function_tables != other.implicit_function_tables { + return false; + } + + true + } +} + +impl Eq for SymbolicContext {} + #[cfg(test)] mod tests { use crate::biodivine_std::traits::Set; use crate::symbolic_async_graph::{SymbolicAsyncGraph, SymbolicContext}; - use crate::BooleanNetwork; + use crate::trap_spaces::SymbolicSpaceContext; + use crate::{BooleanNetwork, VariableId}; use biodivine_lib_bdd::BddPartialValuation; use std::collections::{HashMap, HashSet}; use std::convert::TryFrom; @@ -631,4 +761,64 @@ mod tests { assert!(!negative_bdd.is_false()); } } + + #[test] + fn test_canonical_context() { + let bn = BooleanNetwork::try_from( + r" + a -> b + a -> c + b -> c + b -> a + c -> a + c -> b + $a: b | c + $c: g(a, b) + ", + ) + .unwrap(); + + let extra = bn + .variables() + .enumerate() + .map(|(i, var)| (var, i as u16)) + .collect::>(); + let extra_context = SymbolicContext::with_extra_state_variables(&bn, &extra).unwrap(); + let canonical_context = SymbolicContext::new(&bn).unwrap(); + + assert_ne!(canonical_context, extra_context); + let normalized = extra_context.as_canonical_context(); + assert_eq!( + normalized.bdd.variables(), + canonical_context.bdd.variables() + ); + assert_eq!( + normalized.state_variables, + canonical_context.state_variables + ); + assert_eq!( + normalized.extra_state_variables, + canonical_context.extra_state_variables + ); + assert_eq!( + normalized.explicit_function_tables, + canonical_context.explicit_function_tables + ); + assert_eq!( + normalized.implicit_function_tables, + canonical_context.implicit_function_tables + ); + assert_eq!(canonical_context, extra_context.as_canonical_context()); + + let space_ctx = SymbolicSpaceContext::new(&bn); + let space_stg = SymbolicAsyncGraph::with_space_context(&bn, &space_ctx).unwrap(); + let unit_set = space_stg + .unit_colored_vertices() + .fix_network_variable(VariableId(0), true); + + let normal_stg = SymbolicAsyncGraph::new(&bn).unwrap(); + let translated = normal_stg.transfer_from(&unit_set, &space_stg).unwrap(); + + assert_eq!(translated.exact_cardinality(), unit_set.exact_cardinality()); + } } diff --git a/src/symbolic_async_graph/mod.rs b/src/symbolic_async_graph/mod.rs index a53130e..52b967d 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. /// @@ -167,8 +169,11 @@ pub struct SymbolicContext { /// The main functionality of a `FunctionTable` is that it provides an iterator over /// pairs of `Vec` (function input assignment) and `BddVariable` /// (corresponding symbolic variable). -#[derive(Debug, Clone)] +#[derive(Debug, Clone, PartialEq, Eq)] pub struct FunctionTable { + /// Original name of the function. This can be derived from the names of the symbolic variables, but it + /// is useful to have it around if we ever want to translate between representations. + name: String, pub arity: u16, rows: Vec, } 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; + } + } +}