A binary which tests the implementation by computing weakly connected…

… components using different encodings.

examples/test_reachability.rs

@@ -0,0 +1,15 @@
+use biodivine_lib_logical_models::{BinaryIntegerDomain, GrayCodeIntegerDomain, PetriNetIntegerDomain, reachability_benchmark, UnaryIntegerDomain};
+
+fn main() {
+    let args = std::env::args().collect::<Vec<_>>();
+    let representation = args[1].clone();
+    let sbml_path = args[2].clone();
+
+    match representation.as_str() {
+        "unary" => reachability_benchmark::<UnaryIntegerDomain>(sbml_path.as_str()),
+        "binary" => reachability_benchmark::<BinaryIntegerDomain<u8>>(sbml_path.as_str()),
+        "petri_net" => reachability_benchmark::<PetriNetIntegerDomain>(sbml_path.as_str()),
+        "gray" | "grey" => reachability_benchmark::<GrayCodeIntegerDomain<u8>>(sbml_path.as_str()),
+        _ => panic!("Unknown representation: {}.", representation),
+    }
+}

src/lib.rs

@@ -11,6 +11,9 @@ pub use symbolic_domain::{
     // GenericStateSpaceDomain,
     SymbolicDomain,
     UnaryIntegerDomain,
+    PetriNetIntegerDomain,
+    BinaryIntegerDomain,
+    GrayCodeIntegerDomain,
 };
 
 pub fn add(x: i32, y: i32) -> i32 {

src/prototype/mod.rs

@@ -21,3 +21,6 @@ pub use smart_system_update_fn::*;
 
 mod symbolic_transition_fn;
 pub use symbolic_transition_fn::*;
+
+mod reachability;
+pub use reachability::*;

src/prototype/reachability.rs

@@ -0,0 +1,114 @@
+use std::fmt::Debug;
+use biodivine_lib_bdd::{Bdd, BddPartialValuation};
+use crate::{SmartSystemUpdateFn, SymbolicDomain};
+
+pub fn reachability_benchmark<D: SymbolicDomain<u8> + Debug>(sbml_path: &str) {
+    let smart_system_update_fn = {
+        let file = std::fs::File::open(sbml_path.clone())
+            .expect("Cannot open SBML file.");
+        let reader = std::io::BufReader::new(file);
+        let mut xml = xml::reader::EventReader::new(reader);
+
+        crate::find_start_of(&mut xml, "listOfTransitions")
+            .expect("Cannot find transitions in the SBML file.");
+
+        let smart_system_update_fn = SmartSystemUpdateFn::<D, u8>::try_from_xml(&mut xml)
+                .expect("Loading system fn update failed.");
+
+        smart_system_update_fn
+    };
+
+    let mut universe = smart_system_update_fn.get_bdd_variable_set().mk_true();
+    while !universe.is_false() {
+        let mut weak_scc = pick_state(&smart_system_update_fn, &universe);
+        loop {
+            let bwd_reachable = reach_bwd(&smart_system_update_fn, &weak_scc, &universe);
+            let fwd_bwd_reachable = reach_fwd(&smart_system_update_fn, &bwd_reachable, &universe);
+
+            // FWD/BWD reachable set is not a subset of weak SCC, meaning the SCC can be expanded.
+            if !fwd_bwd_reachable.imp(&weak_scc).is_true() {
+                println!(" + SCC increased to ({}%%, size={})", log_percent(&weak_scc, &universe), weak_scc.size());
+                weak_scc = fwd_bwd_reachable;
+            } else {
+                break;
+            }
+        }
+        println!(" + Found weak SCC (card={}%%, size={})", log_percent(&weak_scc, &universe), weak_scc.size());
+
+        // Remove the SCC from the universe set and start over.
+        universe = universe.and_not(&weak_scc);
+    }
+}
+
+/// Compute the set of vertices that are forward-reachable from the `initial` set.
+///
+/// The result BDD contains a vertex `x` if and only if there is a (possibly zero-length) path
+/// from some vertex `x' \in initial` into `x`, i.e. `x' -> x`.
+pub fn reach_fwd<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D, u8>, initial: &Bdd, universe: &Bdd) -> Bdd {
+    // The list of system variables, sorted in descending order (i.e. opposite order compared
+    // to the ordering inside BDDs).
+    let sorted_variables = system.get_system_variables();
+    let mut result = initial.clone();
+    println!("Start forward reachability: (card={}%%, size={})", log_percent(&result, &universe), result.size());
+    'fwd: loop {
+        for var in sorted_variables.iter().rev() {
+            let successors = system.transition_under_variable(var.as_str(), &result);
+
+            // Should be equivalent to "successors \not\subseteq result".
+            if !successors.imp(&result).is_true() {
+                result = result.or(&successors);
+                println!(" >> (card={}%%, size={})", log_percent(&result, &universe), result.size());
+                continue 'fwd;
+            }
+        }
+
+        // No further successors were computed across all variables. We are done.
+        println!(" >> Done. (card={}%%, size={})", log_percent(&result, &universe), result.size());
+        return result;
+    }
+}
+
+/// Compute the set of vertices that are backward-reachable from the `initial` set.
+///
+/// The result BDD contains a vertex `x` if and only if there is a (possibly zero-length) path
+/// from `x` into some vertex `x' \in initial`, i.e. `x -> x'`.
+pub fn reach_bwd<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D, u8>, initial: &Bdd, universe: &Bdd) -> Bdd {
+    let sorted_variables = system.get_system_variables();
+    let mut result = initial.clone();
+    println!("Start backward reachability: (card={}%%, size={})", log_percent(&result, &universe), result.size());
+    'bwd: loop {
+        for var in sorted_variables.iter().rev() {
+            let predecessors = system.predecessors_under_variable(var.as_str(), &result);
+
+            // Should be equivalent to "predecessors \not\subseteq result".
+            if !predecessors.imp(&result).is_true() {
+                result = result.or(&predecessors);
+                println!(" >> (card={}%%, size={})", log_percent(&result, &universe), result.size());
+                continue 'bwd;
+            }
+        }
+
+        // No further predecessors were computed across all variables. We are done.
+        println!(" >> Done. (card={}%%, size={})", log_percent(&result, &universe), result.size());
+        return result;
+    }
+}
+
+/// Compute a [Bdd] which represents a single (un-primed) state within the given symbolic `set`.
+pub fn pick_state<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D, u8>, set: &Bdd) -> Bdd {
+    // Unfortunately, this is now a bit more complicated than it needs to be, because
+    // we have to ignore the primed variables, but it shouldn't bottleneck anything outside of
+    // truly extreme cases.
+    let standard_variables = system.standard_variables();
+    let valuation = set.sat_witness()
+        .expect("Cannot pick state from an empty set.");
+    let mut state_data = BddPartialValuation::empty();
+    for var in standard_variables {
+        state_data.set_value(var, valuation.value(var))
+    }
+    system.get_bdd_variable_set().mk_conjunctive_clause(&state_data)
+}
+
+pub fn log_percent(set: &Bdd, universe: &Bdd) -> f64 {
+    set.cardinality().log2() / universe.cardinality().log2() * 100.0
+}

src/prototype/smart_system_update_fn.rs

@@ -4,7 +4,7 @@
 
 use std::{collections::HashMap, fmt::Debug, io::BufRead};
 
-use biodivine_lib_bdd::{Bdd, BddPartialValuation, BddVariableSet, BddVariableSetBuilder};
+use biodivine_lib_bdd::{Bdd, BddPartialValuation, BddVariable, BddVariableSet, BddVariableSetBuilder};
 use debug_ignore::DebugIgnore;
 
 use crate::{
@@ -13,7 +13,7 @@ use crate::{
 };
 
 #[derive(Debug)]
-struct SmartSystemUpdateFn<D: SymbolicDomain<T>, T> {
+pub struct SmartSystemUpdateFn<D: SymbolicDomain<T>, T> {
     pub update_fns: HashMap<String, SymbolicTransitionFn<D, T>>,
     pub named_symbolic_domains: HashMap<String, D>,
     bdd_variable_set: DebugIgnore<BddVariableSet>,
@@ -328,21 +328,21 @@ impl<D: SymbolicDomain<u8> + Debug> SmartSystemUpdateFn<D, u8> {
             .transition_function
             .and(&current_state_primed);
 
-        println!(
+        /*println!(
             "is true {}",
             states_capable_of_transitioning_into_current.is_true()
-        );
+        );*/
 
         target_symbolic_domain_primed
             .symbolic_variables()
             .into_iter()
             .fold(
                 states_capable_of_transitioning_into_current,
                 |acc, primed_variable| {
-                    println!(
+                    /*println!(
                         "restricting variable with name {}",
                         self.bdd_variable_set.name_of(primed_variable)
-                    );
+                    );*/
                     acc.var_exists(primed_variable)
                 },
             )
@@ -426,6 +426,40 @@ impl<D: SymbolicDomain<u8> + Debug> SmartSystemUpdateFn<D, u8> {
         // constrain this specific sym variable to its specific value (& leave others unrestricted)
         vars_and_bits.fold(const_true, |acc, (var, bit)| acc.var_select(var, bit))
     }
+
+    /// The list of system variables, sorted in ascending order (i.e. the order in which they
+    /// also appear within the BDDs).
+    pub fn get_system_variables(&self) -> Vec<String> {
+        let mut variables = self.update_fns.keys()
+            .cloned()
+            .collect::<Vec<_>>();
+        variables.sort();
+        variables
+    }
+
+    /// Returns a list of [BddVariable]-s corresponding to the encoding of the "primed"
+    /// system variables.
+    pub fn primed_variables(&self) -> Vec<BddVariable> {
+        let mut result = Vec::new();
+        for (name, domain) in &self.named_symbolic_domains {
+            if name.contains("\'") {
+                result.append(&mut domain.symbolic_variables());
+            }
+        }
+        result
+    }
+
+    /// Returns a list of [BddVariable]-s corresponding to the encoding of the standard
+    /// (i.e. "un-primed") system variables.
+    pub fn standard_variables(&self) -> Vec<BddVariable> {
+        let mut result = Vec::new();
+        for (name, domain) in &self.named_symbolic_domains {
+            if !name.contains("\'") {
+                result.append(&mut domain.symbolic_variables());
+            }
+        }
+        result
+    }
 }
 
 /// expects the xml reader to be at the start of the <listOfTransitions> element