refactor: rewritten the reachability benchmark

examples/reachability.rs

@@ -8,11 +8,15 @@ fn main() {
     let representation = args[1].clone();
     let sbml_path = args[2].clone();
 
+    let now = std::time::Instant::now();
+
     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),
     }
+
+    println!("Time: {}s", now.elapsed().as_secs());
 }

examples/rewritten_reachability.rs

@@ -0,0 +1,27 @@
+use biodivine_lib_logical_models::{
+    benchmarks::rewritten_reachability::reachability_benchmark,
+    prelude::symbolic_domain::{
+        BinaryIntegerDomain, GrayCodeIntegerDomain, PetriNetIntegerDomain, UnaryIntegerDomain,
+    },
+};
+// use biodivine_lib_logical_models::prelude::old_symbolic_domain::{
+//     BinaryIntegerDomain, GrayCodeIntegerDomain, PetriNetIntegerDomain, UnaryIntegerDomain,
+// };
+
+fn main() {
+    let args = std::env::args().collect::<Vec<_>>();
+    let representation = args[1].clone();
+    let sbml_path = args[2].clone();
+
+    let now = std::time::Instant::now();
+
+    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),
+    }
+
+    println!("Time: {}s", now.elapsed().as_secs());
+}

src/benchmarks/mod.rs

@@ -1 +1,3 @@
 pub mod reachability;
+pub mod rewritten_reachability;
+mod utils;

src/benchmarks/rewritten_reachability.rs

@@ -0,0 +1,165 @@
+use biodivine_lib_bdd::Bdd;
+use std::fmt::Debug;
+
+use crate::{
+    benchmarks::utils::{count_states, log_percent, pick_state_bdd, unit_vertex_set},
+    prelude::find_start_of,
+    symbolic_domains::symbolic_domain::SymbolicDomainOrd,
+    update::update_fn::SmartSystemUpdateFn as RewrittenSmartSystemUpdateFn,
+};
+
+pub fn reachability_benchmark<DO: SymbolicDomainOrd<u8> + Debug>(sbml_path: &str) {
+    let smart_system_update_fn = {
+        let mut xml = xml::reader::EventReader::new(std::io::BufReader::new(
+            std::fs::File::open(sbml_path).expect("should be able to open file"),
+        ));
+
+        find_start_of(&mut xml, "listOfTransitions")
+            .expect("Cannot find transitions in the SBML file.");
+
+        RewrittenSmartSystemUpdateFn::<DO, u8>::try_from_xml(&mut xml)
+            .expect("Loading system fn update failed.")
+    };
+
+    let unit = unit_vertex_set(&smart_system_update_fn);
+    let system_var_count = smart_system_update_fn
+        .variables_transition_relation_and_domain
+        .len();
+    println!(
+        "Variables: {}, expected states {}",
+        system_var_count,
+        1 << system_var_count
+    );
+    println!(
+        "Computed state count: {}",
+        count_states(&smart_system_update_fn, &unit)
+    );
+    let mut universe = unit.clone();
+    while !universe.is_false() {
+        let mut weak_scc = pick_state_bdd(&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 (states={}, size={})",
+                    count_states(&smart_system_update_fn, &weak_scc),
+                    weak_scc.size()
+                );
+                weak_scc = fwd_bwd_reachable;
+            } else {
+                break;
+            }
+        }
+        println!(
+            " + Found weak SCC (states={}, size={})",
+            count_states(&smart_system_update_fn, &weak_scc),
+            weak_scc.size()
+        );
+        // Remove the SCC from the universe set and start over.
+        universe = universe.and_not(&weak_scc);
+        println!(
+            " + Remaining states: {}/{}",
+            count_states(&smart_system_update_fn, &universe),
+            count_states(&smart_system_update_fn, &unit),
+        );
+    }
+}
+
+/// 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: SymbolicDomainOrd<u8> + Debug>(
+    system: &RewrittenSmartSystemUpdateFn<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
+        .variables_transition_relation_and_domain
+        .iter()
+        .map(|(var_name, _)| var_name)
+        .collect::<Vec<_>>();
+    let mut result = initial.clone();
+    println!(
+        "Start forward reachability: (states={}, size={})",
+        count_states(system, &result),
+        result.size()
+    );
+    'fwd: loop {
+        for var in sorted_variables.iter().rev() {
+            let successors = system.successors_async(var.as_str(), &result);
+
+            // Should be equivalent to "successors \not\subseteq result".
+            if !successors.imp(&result).is_true() {
+                result = result.or(&successors);
+                println!(
+                    " >> (progress={:.2}%%, states={}, size={})",
+                    log_percent(&result, universe),
+                    count_states(system, &result),
+                    result.size()
+                );
+                continue 'fwd;
+            }
+        }
+
+        // No further successors were computed across all variables. We are done.
+        println!(
+            " >> Done. (states={}, size={})",
+            count_states(system, &result),
+            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: SymbolicDomainOrd<u8> + Debug>(
+    system: &RewrittenSmartSystemUpdateFn<D, u8>,
+    initial: &Bdd,
+    universe: &Bdd,
+) -> Bdd {
+    let sorted_variables = system
+        .variables_transition_relation_and_domain
+        .iter()
+        .map(|(var_name, _)| var_name)
+        .collect::<Vec<_>>();
+    let mut result = initial.clone();
+    println!(
+        "Start backward reachability: (states={}, size={})",
+        count_states(system, &result),
+        result.size()
+    );
+    'bwd: loop {
+        for var in sorted_variables.iter().rev() {
+            let predecessors = system.predecessors_async(var.as_str(), result.clone());
+
+            // Should be equivalent to "predecessors \not\subseteq result".
+            if !predecessors.imp(&result).is_true() {
+                result = result.or(&predecessors);
+                println!(
+                    " >> (progress={:.2}%%, states={}, size={})",
+                    log_percent(&result, universe),
+                    count_states(system, &result),
+                    result.size()
+                );
+                continue 'bwd;
+            }
+        }
+
+        // No further predecessors were computed across all variables. We are done.
+        println!(
+            " >> Done. (states={}, size={})",
+            count_states(system, &result),
+            result.size()
+        );
+        return result;
+    }
+}

src/benchmarks/utils.rs

@@ -0,0 +1,67 @@
+use std::fmt::Debug;
+
+use biodivine_lib_bdd::{Bdd, BddPartialValuation};
+
+use crate::{
+    symbolic_domains::symbolic_domain::SymbolicDomainOrd,
+    update::update_fn::SmartSystemUpdateFn as RewrittenSmartSystemUpdateFn,
+};
+
+pub fn states<D: SymbolicDomainOrd<u8>>(
+    system: &RewrittenSmartSystemUpdateFn<D, u8>,
+    set: &Bdd,
+) -> f64 {
+    let symbolic_var_count = system.bdd_variable_set.num_vars() as i32;
+    // TODO:
+    //   Here we assume that exactly half of the variables are primed, which may not be true
+    //   in the future, but should be good enough for now.
+    assert_eq!(symbolic_var_count % 2, 0);
+    let primed_vars = symbolic_var_count / 2;
+    set.cardinality() / 2.0f64.powi(primed_vars)
+}
+
+pub fn unit_vertex_set<D: SymbolicDomainOrd<u8>>(
+    system: &RewrittenSmartSystemUpdateFn<D, u8>,
+) -> Bdd {
+    system
+        .variables_transition_relation_and_domain
+        .iter()
+        .fold(system.bdd_variable_set.mk_true(), |acc, (_, var_info)| {
+            acc.and(&var_info.domain.unit_collection(&system.bdd_variable_set))
+        })
+}
+
+/// Compute an (approximate) count of state in the given `set` using the encoding of `system`.
+pub fn count_states<D: SymbolicDomainOrd<u8> + Debug>(
+    system: &RewrittenSmartSystemUpdateFn<D, u8>,
+    set: &Bdd,
+) -> f64 {
+    let symbolic_var_count = system.variables_transition_relation_and_domain.len() as i32;
+    set.cardinality() / 2.0f64.powi(symbolic_var_count)
+}
+
+/// Compute a [Bdd] which represents a single (un-primed) state within the given symbolic `set`.
+pub fn pick_state_bdd<D: SymbolicDomainOrd<u8> + Debug>(
+    system: &RewrittenSmartSystemUpdateFn<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
+        .variables_transition_relation_and_domain
+        .iter()
+        .flat_map(|transition| transition.1.domain.raw_bdd_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.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/utils.rs

@@ -11,6 +11,7 @@ use xml::{
 };
 
 use crate::prototype::symbolic_domain::SymbolicDomain;
+use crate::symbolic_domains::symbolic_domain::SymbolicDomain as RewrittenSymbolicDomain;
 
 use super::{SmartSystemUpdateFn, UpdateFn};
 

src/symbolic_domains/symbolic_domain.rs

@@ -267,6 +267,7 @@ impl SymbolicDomainOrd<u8> for UnaryIntegerDomain {
     }
 }
 
+#[derive(Debug)]
 pub struct PetriNetIntegerDomain {
     /// invariant: sorted
     variables: Vec<BddVariable>,
@@ -355,6 +356,7 @@ impl SymbolicDomainOrd<u8> for PetriNetIntegerDomain {
     }
 }
 
+#[derive(Debug)]
 pub struct BinaryIntegerDomain<T> {
     /// invariant: sorted
     variables: Vec<BddVariable>,
@@ -460,6 +462,7 @@ impl SymbolicDomainOrd<u8> for BinaryIntegerDomain<u8> {
     }
 }
 
+#[derive(Debug)]
 pub struct GrayCodeIntegerDomain<T> {
     /// invariant: sorted
     variables: Vec<BddVariable>,

src/update/update_fn.rs

@@ -326,12 +326,13 @@ where
     }
 }
 
-struct VarInfo<D, T>
+pub struct VarInfo<D, T>
+// todo do not keep pub; just for benchmarks
 where
     D: SymbolicDomain<T>,
 {
     primed_name: String,
-    domain: D,
+    pub domain: D, // todo do not keep pub; just for benchmarks
     primed_domain: D,
     transition_relation: Bdd,
     _marker: std::marker::PhantomData<T>,
@@ -342,8 +343,8 @@ where
     D: SymbolicDomain<T>,
 {
     /// ordered by variable name // todo add a method to get the update function by name (hash map or binary search)
-    variables_transition_relation_and_domain: Vec<(String, VarInfo<D, T>)>,
-    bdd_variable_set: BddVariableSet,
+    pub variables_transition_relation_and_domain: Vec<(String, VarInfo<D, T>)>, // todo do not keep pub; just here for benchmarking
+    pub bdd_variable_set: BddVariableSet, // todo do not keep pub; just here for testing
     _marker: std::marker::PhantomData<T>,
 }
 

tests/some_test.rs

@@ -745,7 +745,7 @@ fn predecessors_consistency_check() {
 
                     assert!(initial_state.are_same(&the_four), "initial states are same");
 
-                    let transitioned = the_four.successors_async(variable, initial_state);
+                    let transitioned = the_four.predecessors_async(variable, initial_state);
 
                     assert!(transitioned.are_same(&the_four), "all are same");