Merge pull request #2 from sybila/reachability-test

Reachability test

Cargo.toml

@@ -10,11 +10,12 @@ name = "biodivine_lib_logical_models"
 path = "src/lib.rs"
 
 [dependencies]
-biodivine-lib-bdd = "0.5.1"
+biodivine-lib-bdd = "0.5.3"
 debug-ignore = "1.0.5"
 dyn-clonable = "0.9.0"
 rayon = "1.8.0"
 serde = { version = "1.0", features = ["derive"] }
 serde-xml-rs = "0.6.0"
 thiserror = "1.0.40"
 xml-rs = "0.8.14"
+num-bigint = "0.4.4"

README.md

@@ -21,3 +21,34 @@ A Rust library for working with logical models (Boolean/multi-valued networks) i
 - [ ] Can represent unknown/uncertain behaviour within the logical model.
 - [ ] Can represent and manipulate the state-transition graph of a logical model symbolically (maybe using multiple different encodings?).
 - [ ] Provides some basic utility algorithms for (a) exploring the structural properties of the model (feedback vertex sets, cycles, etc.) (b) exploring the model dynamics (reachability, fixed-points, trap spaces, etc.).
+
+
+### Running reachability integration test
+
+To test the correctness of the implementation, we now have a simple reachability consistency check.
+This check compares the results of all symbolic representations on a basic reachability exploration.
+To run the test, execute the following:
+
+```bash
+# First argument is a timeout, here 1 minute. This uses standard unix `timeout`.
+# Second argument is fwd/bwd to indicate which reachability direction you want to test.
+python3 reachability_integration_test.py 1m fwd
+```
+
+The process dumps all results into an appropriate `./data/results-*` directory.
+The script will not overwrite existing results, 
+so you have to delete the directory manually before running the test again.
+
+You should see an output similar to the following:
+
+```
+[PASS] No error discovered in `146_BUDDING-YEAST-FAURE-2009.sbml` in less than 1m.
+[PASS] No error discovered in `148_AGS-cell-fate-decision.sbml` in less than 1m.
+[PASS] No error discovered in `151_TCR-REDOX-METABOLISM.sbml` in less than 1m.
+[PASS] No error discovered in `155_CONTROL-OF-TH1-TH2-TH17-TREG-DIFFERENTATION.sbml` in less than 1m.
+[PASS] No error discovered in `157_CONTROL-OF-TH-DIFFERENTATION.sbml` in less than 1m.
+...
+```
+
+If you see a `[FAIL]` somewhere, it means an inconsistency has been detected.
+You can then check the results to examine the specific case.

data/.gitignore

@@ -0,0 +1 @@
+results-*

examples/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),
+    }
+}

examples/test_reachability_bwd.rs

@@ -0,0 +1,26 @@
+use biodivine_lib_logical_models::test_utils::ComputationStep;
+
+
+/// This binary is testing the implementation correctness by running reachability on the
+/// input model and validating that the set of reachable states has the same cardinality
+/// in every step.
+///
+/// For larger models, this is almost sure to take "forever", hence if you want to use it as
+/// an automated test, you should always run it with a timeout, and ideally with optimizations.
+/// This is also the reason why we don't use it as a normal integration test: because those
+/// run unoptimized by default, and timeout can be only used to fail tests.
+fn main() {
+    let args = std::env::args().collect::<Vec<_>>();
+    let sbml_path = args[1].clone();
+
+    let mut cmp = ComputationStep::new(sbml_path.as_str());
+    while !cmp.is_done() {
+        cmp.initialize();
+        while !cmp.can_initialize() {
+            cmp.perform_bwd_step();
+            cmp.check_consistency();
+        }
+        println!("Completed one wave of reachability. Reinitializing with {} states remaining.", cmp.remaining());
+    }
+    println!("Test completed successfully. Whole state space explored.");
+}

examples/test_reachability_fwd.rs

@@ -0,0 +1,26 @@
+use biodivine_lib_logical_models::test_utils::ComputationStep;
+
+
+/// This binary is testing the implementation correctness by running reachability on the
+/// input model and validating that the set of reachable states has the same cardinality
+/// in every step.
+///
+/// For larger models, this is almost sure to take "forever", hence if you want to use it as
+/// an automated test, you should always run it with a timeout, and ideally with optimizations.
+/// This is also the reason why we don't use it as a normal integration test: because those
+/// run unoptimized by default, and timeout can be only used to fail tests.
+fn main() {
+    let args = std::env::args().collect::<Vec<_>>();
+    let sbml_path = args[1].clone();
+
+    let mut cmp = ComputationStep::new(sbml_path.as_str());
+    while !cmp.is_done() {
+        cmp.initialize();
+        while !cmp.can_initialize() {
+            cmp.perform_fwd_step();
+            cmp.check_consistency();
+        }
+        println!("Completed one wave of reachability. Reinitializing with {} states remaining.", cmp.remaining());
+    }
+    println!("Test completed successfully. Whole state space explored.");
+}

reachability_integration_test.py

@@ -0,0 +1,37 @@
+import os
+import subprocess
+import sys
+
+TIMEOUT = sys.argv[1]
+DIR = sys.argv[2]
+
+if DIR != "fwd" and DIR != "bwd":
+    print("Invalid direction. Allowed `fwd` or `bwd`.")
+    exit(2)
+
+code = os.system(f"cargo build --release --example test_reachability_{DIR}")
+if code != 0:
+    print("Compilation error.")
+    exit(2)
+
+if os.path.isdir(f"./data/results-{DIR}-test"):
+    print("Test results already exist. Won't overwrite existing data.")
+    exit(2)
+
+os.system(f"mkdir -p ./data/results-{DIR}-test")
+
+files = list(os.listdir("./data/test-models"))
+files = list(sorted(files))
+
+for file in files:
+    if not file.endswith(".sbml"):
+        continue
+    name = file.replace(".sbml", "")
+    cmd_run = f"./target/release/examples/test_reachability_{DIR} ./data/test-models/{file} &> ./data/results-{DIR}-test/{name}.txt"
+    code = os.system(f"timeout {TIMEOUT} {cmd_run}")
+    if code == 31744 or code == 124:
+        print(f"[PASS] No error discovered in `{file}` in less than {TIMEOUT}.")
+    elif code != 0:
+        print(f"[FAIL] Error ({code}) when testing `{file}`. See `./data/results-{DIR}-test/{name}.bwd_test.txt` for details.")
+    else:
+        print(f"[PASS] Completed `{file}`.")

src/lib.rs

@@ -5,12 +5,20 @@
 ///     In the final library, we should re-export the relevant types from this module here.
 mod symbolic_domain;
 
+// TODO:
+//   Once this becomes a library, this needs to become private, but for now it is convenient
+//   to have it accessible from outside binaries.
+pub mod test_utils;
+
 pub use symbolic_domain::{
     // todo uncomment one those working
     // GenericIntegerDomain,
     // 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,113 @@
+use std::fmt::Debug;
+use biodivine_lib_bdd::Bdd;
+use crate::{count_states, log_percent, pick_state_bdd, 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 unit = smart_system_update_fn.unit_vertex_set();
+    let system_var_count = smart_system_update_fn.get_system_variables().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: 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: (states={}, size={})", count_states(system, &result), 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!(
+                    " >> (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: 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: (states={}, size={})", count_states(system, &result), 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!(
+                    " >> (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;
+    }
+}