Add wider test for reachability match.

data/test-results/146_BUDDING-YEAST-FAURE-2009.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/148_AGS-cell-fate-decision.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/151_TCR-REDOX-METABOLISM.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/155_CONTROL-OF-TH1-TH2-TH17-TREG-DIFFERENTATION.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(3.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/157_CONTROL-OF-TH-DIFFERENTATION.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(3.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/159_BUDDING-YEAST-CORE.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/160_IL17-DIFFERENTIAL-EXPRESSION.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/161_DIFFERENTIATION-OF-MONOCYTES.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/167_DROSOPHILA-MESODERM.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/175_SEA-URCHIN.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(3.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/176-myelofibrotic-microenvironment.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/178-mast-cell-activation.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/179-microenvironment-control.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/183-alterations-in-bladder.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/190-BRAF-treatment-response.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/192-segment-polarity-6-cell.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/194-vulvar-precursor-cells.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/195-CTLA4-PD1-checkpoint-inhibitors.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/196-T-lymphocyte-specification.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(4.0) <> Some(2.0) <> Some(2.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/197-anterior-posterior-boundary.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/InVitro.free-inputs.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/InVivo.free-inputs.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/Leukaemia.free-inputs.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/Metabolism_demo.free-inputs.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(32.0) <> Some(8.0) <> Some(8.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/SkinModel.free-inputs.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(28.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

data/test-results/VPC.free-inputs.bwd_test.txt

@@ -0,0 +1,2 @@
+thread 'main' panicked at 'Error at step 1. Some(8.0) <> Some(4.0) <> Some(4.0) <> "??"', examples/test_reachability_bwd.rs:116:13
+note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

examples/test_reachability_bwd.rs

@@ -0,0 +1,136 @@
+use std::fmt::Debug;
+use biodivine_lib_bdd::Bdd;
+use biodivine_lib_logical_models::{BinaryIntegerDomain, count_states, find_start_of, GrayCodeIntegerDomain, PetriNetIntegerDomain, pick_state, SmartSystemUpdateFn, SymbolicDomain, SymbolicTransitionFn, UnaryIntegerDomain};
+
+/// 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());
+    loop {
+        cmp.perform_step();
+        cmp.check_consistency();
+        if cmp.result_unary.is_none() {
+            // The results must be equal because consistency check passed.
+            println!("Test completed successfully.");
+            return;
+        }
+    }
+}
+
+struct ComputationStep {
+    steps: usize,
+    result_unary: Option<Bdd>,
+    result_binary: Option<Bdd>,
+    result_gray: Option<Bdd>,
+    //result_petri_net: Option<Bdd>,
+    system_unary: SmartSystemUpdateFn<UnaryIntegerDomain, u8>,
+    system_binary: SmartSystemUpdateFn<BinaryIntegerDomain<u8>, u8>,
+    system_gray: SmartSystemUpdateFn<GrayCodeIntegerDomain<u8>, u8>,
+    //system_petri_net: SmartSystemUpdateFn<PetriNetIntegerDomain, u8>,
+}
+
+/// Perform one step of backward reachability procedure. Returns either a new [Bdd] value, or
+/// `None` if no new predecessors can be included.
+fn bwd_step<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D, u8>, set: &Bdd) -> Option<Bdd> {
+    let sorted_variables = system.get_system_variables();
+
+    for var in sorted_variables.iter().rev() {
+        let predecessors = system.predecessors_under_variable(var.as_str(), set);
+
+        // Should be equivalent to "predecessors \not\subseteq result".
+        if !predecessors.imp(set).is_true() {
+            let result = predecessors.or(&set);
+            return Some(result);
+        }
+    }
+
+    None
+}
+
+fn build_update_fn<D: SymbolicDomain<u8> + Debug>(sbml_path: &str) -> SmartSystemUpdateFn<D, u8> {
+    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);
+
+    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
+}
+
+impl ComputationStep {
+
+    pub fn new(sbml_path: &str) -> ComputationStep {
+        let system_unary = build_update_fn::<UnaryIntegerDomain>(sbml_path);
+        let system_binary = build_update_fn::<BinaryIntegerDomain<u8>>(sbml_path);
+        let system_gray = build_update_fn::<GrayCodeIntegerDomain<u8>>(sbml_path);
+        //let system_petri_net = build_update_fn::<PetriNetIntegerDomain>(sbml_path);
+        ComputationStep {
+            steps: 0,
+            result_unary: Some(pick_state(&system_unary, &system_unary.unit_vertex_set())),
+            result_binary: Some(pick_state(&system_binary, &system_binary.unit_vertex_set())),
+            result_gray: Some(pick_state(&system_gray, &system_gray.unit_vertex_set())),
+            //result_petri_net: Some(pick_state(&system_petri_net, &system_petri_net.unit_vertex_set())),
+            system_unary,
+            system_binary,
+            system_gray,
+            //system_petri_net,
+        }
+    }
+
+    pub fn perform_step(&mut self) {
+        self.steps += 1;
+        self.result_unary = bwd_step(&self.system_unary, self.result_unary.as_ref().unwrap());
+        self.result_binary = bwd_step(&self.system_binary, self.result_binary.as_ref().unwrap());
+        self.result_gray = bwd_step(&self.system_gray, self.result_gray.as_ref().unwrap());
+        //self.result_petri_net = bwd_step(&self.system_petri_net, self.result_petri_net.as_ref().unwrap());
+    }
+
+    pub fn check_consistency(&self) {
+        let count_unary = self.result_unary.as_ref().map(|it| {
+            count_states(&self.system_unary, &it)
+        });
+        let count_binary = self.result_binary.as_ref().map(|it| {
+            count_states(&self.system_binary, &it)
+        });
+        let count_gray = self.result_gray.as_ref().map(|it| {
+            count_states(&self.system_gray, &it)
+        });
+        /*let count_petri_net = self.result_petri_net.as_ref().map(|it| {
+            count_states(&self.system_petri_net, &it)
+        });*/
+        if count_unary != count_binary || count_binary != count_gray /*|| count_gray != count_petri_net*/ {
+            panic!(
+                "Error at step {}. {:?} <> {:?} <> {:?} <> {:?}",
+                self.steps,
+                count_unary,
+                count_binary,
+                count_gray,
+                "??"//count_petri_net
+            )
+        } else {
+            println!("Step {} successful. Current result state count: {:?}", self.steps, count_unary);
+            println!(
+                " > BDD sizes: {:?} {:?} {:?} {:?}",
+                self.result_unary.as_ref().map(|it| it.size()),
+                self.result_binary.as_ref().map(|it| it.size()),
+                self.result_gray.as_ref().map(|it| it.size()),
+                "??"//self.result_petri_net.as_ref().map(|it| it.size()),
+            );
+        }
+    }
+
+}

reachability_integration_test.py

@@ -0,0 +1,27 @@
+import os
+import subprocess
+
+TIMEOUT = "1m"
+
+os.system("cargo build --release --example test_reachability_bwd")
+if os.path.isdir("./data/test-results"):
+    print("Test results already exist. Won't overwrite existing data.")
+    exit(2)
+
+os.system("mkdir -p ./data/test-results")
+
+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_bwd ./data/test-models/{file} &> ./data/test-results/{name}.bwd_test.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/test-results/{name}.bwd_test.txt` for details.")
+    else:
+        print(f"[PASS] Completed `{file}`.")

src/prototype/reachability.rs

@@ -18,7 +18,11 @@ pub fn reachability_benchmark<D: SymbolicDomain<u8> + Debug>(sbml_path: &str) {
         smart_system_update_fn
     };
 
-    let mut universe = smart_system_update_fn.get_bdd_variable_set().mk_true();
+    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(&smart_system_update_fn, &universe);
         loop {
@@ -27,16 +31,20 @@ pub fn reachability_benchmark<D: SymbolicDomain<u8> + Debug>(sbml_path: &str) {
 
             // 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());
+                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 (card={}%%, size={})", log_percent(&weak_scc, &universe), weak_scc.size());
-
+        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),
+        );
     }
 }
 
@@ -49,21 +57,26 @@ pub fn reach_fwd<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D,
     // 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());
+    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!(" >> (card={}%%, size={})", log_percent(&result, &universe), result.size());
+                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. (card={}%%, size={})", log_percent(&result, &universe), result.size());
+        println!(" >> Done. (states={}, size={})", count_states(system, &result), result.size());
         return result;
     }
 }
@@ -75,21 +88,26 @@ pub fn reach_fwd<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D,
 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());
+    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!(" >> (card={}%%, size={})", log_percent(&result, &universe), result.size());
+                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. (card={}%%, size={})", log_percent(&result, &universe), result.size());
+        println!(" >> Done. (states={}, size={})", count_states(system, &result), result.size());
         return result;
     }
 }
@@ -111,4 +129,15 @@ pub fn pick_state<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D,
 
 pub fn log_percent(set: &Bdd, universe: &Bdd) -> f64 {
     set.cardinality().log2() / universe.cardinality().log2() * 100.0
+}
+
+/// Compute an (approximate) count of state in the given `set` using the encoding of `system`.
+pub fn count_states<D: SymbolicDomain<u8> + Debug>(system: &SmartSystemUpdateFn<D, u8>, set: &Bdd) -> f64 {
+    let symbolic_var_count = system.get_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)
 }