From 1bdb60dc8c415f6edcf4c60136d37bd6ff74697d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20Huvar?= <492849@mail.muni.cz> Date: Mon, 26 Feb 2024 21:11:04 +0100 Subject: [PATCH] Start working on analysis of HCTL with placeholders. --- Cargo.toml | 19 +- src/bin/case_study_tlgl.rs | 23 +- src/bin/small_example.rs | 8 +- src/bin/test_hctl_with_1_hole.rs | 100 +++++++ src/bin/test_hctl_with_2_holes.rs | 41 +++ src/data_processing/data_encoding.rs | 34 ++- src/data_processing/observations.rs | 41 ++- src/hctl_with_holes/mod.rs | 268 ++++++++++++++++++ src/hctl_with_holes/reachability_one_hole.rs | 216 ++++++++++++++ src/hctl_with_holes/reachability_two_holes.rs | 57 ++++ src/hctl_with_holes/utils.rs | 60 ++++ src/inference_attractor_data.rs | 4 +- src/lib.rs | 1 + src/utils.rs | 2 +- 14 files changed, 839 insertions(+), 35 deletions(-) create mode 100644 src/bin/test_hctl_with_1_hole.rs create mode 100644 src/bin/test_hctl_with_2_holes.rs create mode 100644 src/hctl_with_holes/mod.rs create mode 100644 src/hctl_with_holes/reachability_one_hole.rs create mode 100644 src/hctl_with_holes/reachability_two_holes.rs create mode 100644 src/hctl_with_holes/utils.rs diff --git a/Cargo.toml b/Cargo.toml index 644b17e..f0fdccb 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,10 +11,7 @@ keywords = ["inference", "boolean-network", "model-checking", "symbolic", "syste categories = ["science", "simulation"] license = "MIT" exclude = ["benchmark_models", ".github", ".githooks"] -rust-version = "1.72" - -#[profile.release] -#lto = true +rust-version = "1.76" [profile.test] opt-level = 3 @@ -39,10 +36,18 @@ path = "src/bin/case_study_arabidopsis.rs" name = "small-example" path = "src/bin/small_example.rs" +[[bin]] +name = "test-hctl-with-1-hole" +path = "src/bin/test_hctl_with_1_hole.rs" + +[[bin]] +name = "test-hctl-with-2-holes" +path = "src/bin/test_hctl_with_2_holes.rs" + [dependencies] -biodivine-lib-bdd = ">=0.5.7, <1.0.0" -biodivine-lib-param-bn = ">=0.5.1, <1.0.0" -biodivine-hctl-model-checker = { git = "https://github.com/sybila/biodivine-hctl-model-checker.git", branch = "master" } +biodivine-lib-bdd = ">=0.5.11, <0.6.0" +biodivine-lib-param-bn = ">=0.5.9, <0.6.0" +biodivine-hctl-model-checker = { git = "https://github.com/sybila/biodivine-hctl-model-checker.git", rev = "1e893fa" } clap = { version = "4.1.4", features = ["derive"] } termcolor = "1.1.2" rand = "0.8.5" \ No newline at end of file diff --git a/src/bin/case_study_tlgl.rs b/src/bin/case_study_tlgl.rs index 9563aef..df2c702 100644 --- a/src/bin/case_study_tlgl.rs +++ b/src/bin/case_study_tlgl.rs @@ -64,9 +64,10 @@ fn case_study_part_1() { // define data observation and corresponding dynamic property let diseased_attractor = "~Apoptosis_ & S1P & sFas & ~Fas & ~Ceramide_ & ~Caspase & MCL1 & ~BID_ & ~DISC_ & FLIP_ & ~IFNG_ & GPCR_"; let formulae: Vec = vec![mk_formula_attractor(diseased_attractor.to_string())]; + let formulae_slices = formulae.iter().map(|s| s.as_str()).collect(); // apply dynamic constraints - graph = apply_constraints_and_restrict(formulae, graph, "attractor property ensured"); + graph = apply_constraints_and_restrict(formulae_slices, graph, "attractor property ensured"); println!( "{} consistent candidate networks found in total.", graph.mk_unit_colors().approx_cardinality(), // graph has restricted unit colors to satisfying ones @@ -81,8 +82,7 @@ fn case_study_part_1() { println!("Analysing candidate set..."); // compute attractors symbolically - let attrs_all_candidates = - model_check_formula_dirty("!{x}: AG EF {x}".to_string(), &graph).unwrap(); + let attrs_all_candidates = model_check_formula_dirty("!{x}: AG EF {x}", &graph).unwrap(); println!("Attractors for all candidates computed"); println!( "Elapsed time from the start of this computation: {}ms", @@ -92,7 +92,7 @@ fn case_study_part_1() { // check for candidates without attractor for programmed cell death let programmed_cell_death_formula = "Apoptosis_ & ~S1P & ~sFas & ~Fas & ~Ceramide_ & ~Caspase & ~MCL1 & ~BID_ & ~DISC_ & ~FLIP_ & ~CTLA4_ & ~TCR & ~IFNG_ & ~CREB & ~P2 & ~SMAD_ & ~GPCR_ & ~IAP_"; - let pcd = model_check_formula_dirty(programmed_cell_death_formula.to_string(), &graph).unwrap(); + let pcd = model_check_formula_dirty(programmed_cell_death_formula, &graph).unwrap(); let colors_not_pcd = graph .mk_unit_colors() .minus(&attrs_all_candidates.intersect(&pcd).colors()); @@ -108,8 +108,7 @@ fn case_study_part_1() { // check for candidates with unwanted attractor states let unwanted_state_formula = "Apoptosis_ & (S1P | sFas | Fas | Ceramide_ | Caspase | MCL1 | BID_ | DISC_ | FLIP_ | CTLA4_ | TCR | IFNG_ | CREB | P2 | SMAD_ | GPCR_ | IAP_)"; - let unwanted_states = - model_check_formula_dirty(unwanted_state_formula.to_string(), &graph).unwrap(); + let unwanted_states = model_check_formula_dirty(unwanted_state_formula, &graph).unwrap(); let colors_with_unwanted_states = attrs_all_candidates.intersect(&unwanted_states).colors(); println!( "{} candidates have unwanted states in attractors, such as:\n", @@ -167,7 +166,8 @@ fn case_study_part_2(summarize_candidates: bool) { ]; // first ensure attractor existence - graph = apply_constraints_and_restrict(formulae, graph, "attractor property ensured"); + let formulae_slice = formulae.iter().map(|s| s.as_str()).collect(); + graph = apply_constraints_and_restrict(formulae_slice, graph, "attractor property ensured"); println!( "After ensuring both properties regarding attractor presence, {} candidates remain.", graph.mk_unit_colors().approx_cardinality(), @@ -179,7 +179,9 @@ fn case_study_part_2(summarize_candidates: bool) { diseased_attractor.to_string(), ]; let formula = mk_formula_forbid_other_attractors(attr_set); - let inferred_colors = model_check_formula_dirty(formula, &graph).unwrap().colors(); + let inferred_colors = model_check_formula_dirty(&formula, &graph) + .unwrap() + .colors(); println!( "{} consistent candidate networks found in total", inferred_colors.approx_cardinality() @@ -253,16 +255,17 @@ mod tests { mk_formula_fixed_point_specific(healthy_attractor.to_string()), mk_formula_attractor(diseased_attractor.to_string()), ]; + let formulae_slices = formulae.iter().map(|s| s.as_str()).collect(); // first ensure attractor existence - graph = apply_constraints_and_restrict(formulae, graph, "attractor ensured"); + graph = apply_constraints_and_restrict(formulae_slices, graph, "attractor ensured"); // prohibit all other attractors let attr_set = vec![ healthy_attractor.to_string(), diseased_attractor.to_string(), ]; let formula = mk_formula_forbid_other_attractors(attr_set); - let inferred_colors = model_check_formula(formula, &graph).unwrap().colors(); + let inferred_colors = model_check_formula(&formula, &graph).unwrap().colors(); assert_eq!(inferred_colors.approx_cardinality(), 378.); } } diff --git a/src/bin/small_example.rs b/src/bin/small_example.rs index 286fc39..523e48f 100644 --- a/src/bin/small_example.rs +++ b/src/bin/small_example.rs @@ -29,7 +29,7 @@ fn main() { // prior-knowledge dynamic properties only let prior_formula = "3{a}: (3{b}: (3{c}: (@{c}: ((EF {a}) & (EF {b}) & (@{a}: AG EF {a}) & (@{b}: (AG EF {b} & ~ EF {a}))))))"; - let intermediate_result = model_check_formula(prior_formula.to_string(), &graph).unwrap(); + let intermediate_result = model_check_formula(prior_formula, &graph).unwrap(); println!( "After applying prior-knowledge-based dynamic constraints, {} candidates remain.", intermediate_result.colors().approx_cardinality(), @@ -37,7 +37,7 @@ fn main() { // properties regarding both prior knowledge and data let whole_formula = "3{a}: (3{b}: (3{c}: (@{c}: ((EF {a}) & (EF {b}) & (@{a}: AG EF {a}) & (@{b}: (AG EF {b} & ~ EF {a})))))) & (3{x}:@{x}: ~v_1 & ~v_2 & v_3 & AG EF {x}) & (3{y}:@{y}: v_1 & v_2 & ~v_3 & AG EF {y})"; - let result = model_check_formula(whole_formula.to_string(), &graph).unwrap(); + let result = model_check_formula(whole_formula, &graph).unwrap(); let res_color = result.colors(); println!( @@ -80,11 +80,11 @@ mod tests { assert_eq!(graph.mk_unit_colors().approx_cardinality(), 16.); let prior_formula = "3{a}: (3{b}: (3{c}: (@{c}: ((EF {a}) & (EF {b}) & (@{a}: AG EF {a}) & (@{b}: (AG EF {b} & ~ EF {a}))))))"; - let intermediate_result = model_check_formula(prior_formula.to_string(), &graph).unwrap(); + let intermediate_result = model_check_formula(prior_formula, &graph).unwrap(); assert_eq!(intermediate_result.colors().approx_cardinality(), 4.); let whole_formula = "3{a}: (3{b}: (3{c}: (@{c}: ((EF {a}) & (EF {b}) & (@{a}: AG EF {a}) & (@{b}: (AG EF {b} & ~ EF {a})))))) & (3{x}:@{x}: ~v_1 & ~v_2 & v_3 & AG EF {x}) & (3{y}:@{y}: v_1 & v_2 & ~v_3 & AG EF {y})"; - let result = model_check_formula(whole_formula.to_string(), &graph).unwrap(); + let result = model_check_formula(whole_formula, &graph).unwrap(); assert_eq!(result.colors().approx_cardinality(), 1.); } } diff --git a/src/bin/test_hctl_with_1_hole.rs b/src/bin/test_hctl_with_1_hole.rs new file mode 100644 index 0000000..e0dd6f2 --- /dev/null +++ b/src/bin/test_hctl_with_1_hole.rs @@ -0,0 +1,100 @@ +use biodivine_hctl_model_checker::mc_utils::get_extended_symbolic_graph; +use biodivine_hctl_model_checker::model_checking::{ + model_check_formula, model_check_formula_dirty, +}; +use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph; +use biodivine_lib_param_bn::BooleanNetwork; +use boolean_network_sketches::data_processing::data_encoding::encode_observation; +use boolean_network_sketches::data_processing::observations::Observation; +use boolean_network_sketches::hctl_with_holes::reachability_one_hole::{ + general_to_specific_on_the_fly_eval, general_to_specific_precomp_eval, naive_eval, + specific_to_general_eval, +}; +use boolean_network_sketches::hctl_with_holes::utils::encode_obs_weight_pairs; +#[allow(unused_imports)] +use boolean_network_sketches::hctl_with_holes::{MODEL_20V, MODEL_53V}; +use std::time::SystemTime; + +fn main() { + let use_large = true; + + let mut model = MODEL_20V; + let mut start_from_ones = false; + if use_large { + model = MODEL_53V; + start_from_ones = true; + } + + let bn = BooleanNetwork::try_from(model).unwrap(); + let props_vec: Vec<_> = bn + .variables() + .map(|v| bn.get_variable_name(v).clone()) + .collect(); + + // for the first two methods, we must use new STG and sets with symbolic variables (for HCTL model checking) + let stg = get_extended_symbolic_graph(&bn, 1).unwrap(); + + // fixed initial observation (either ones or zeros) + let init_observation = if start_from_ones { + Observation::new_full_true(props_vec.len()) + } else { + Observation::new_full_false(props_vec.len()) + }; + let init_observation_str = encode_observation(&init_observation, &props_vec).unwrap(); + let init_observation_set = + model_check_formula_dirty(init_observation_str.as_str(), &stg).unwrap(); + // variants of the second observation and their weights + let second_observation_weight_pairs = + encode_obs_weight_pairs(&props_vec, !start_from_ones).unwrap(); + + let start = SystemTime::now(); + println!("naive_eval"); + let best_weight = naive_eval( + &stg, + start, + &init_observation_set, + second_observation_weight_pairs.clone(), + ); + println!("TOTAL ELAPSED: {}", start.elapsed().unwrap().as_millis()); + println!("BEST WEIGHT: {:?}", best_weight); + print!("\n\n\n"); + + let start = SystemTime::now(); + println!("specific_to_general_eval"); + let best_weight = specific_to_general_eval( + &stg, + start, + &init_observation_set, + second_observation_weight_pairs.clone(), + ); + println!("TOTAL ELAPSED: {}", start.elapsed().unwrap().as_millis()); + println!("BEST WEIGHT: {:?}", best_weight); + print!("\n\n\n"); + + // for the following methods, we must use new STG and sets without symbolic variables + let stg = SymbolicAsyncGraph::new(&bn).unwrap(); + let init_observation_set = model_check_formula(init_observation_str.as_str(), &stg).unwrap(); + + let start = SystemTime::now(); + println!("general_to_specific_precomp_eval"); + let best_weight = general_to_specific_precomp_eval( + &stg, + start, + &init_observation_set, + second_observation_weight_pairs.clone(), + ); + println!("TOTAL ELAPSED: {}", start.elapsed().unwrap().as_millis()); + println!("BEST WEIGHT: {:?}", best_weight); + print!("\n\n\n"); + + let start = SystemTime::now(); + println!("general_to_specific_on_the_fly_eval"); + let best_weight = general_to_specific_on_the_fly_eval( + &stg, + start, + &init_observation_set, + second_observation_weight_pairs.clone(), + ); + println!("TOTAL ELAPSED: {}", start.elapsed().unwrap().as_millis()); + println!("BEST WEIGHT: {:?}", best_weight); +} diff --git a/src/bin/test_hctl_with_2_holes.rs b/src/bin/test_hctl_with_2_holes.rs new file mode 100644 index 0000000..4bd332c --- /dev/null +++ b/src/bin/test_hctl_with_2_holes.rs @@ -0,0 +1,41 @@ +use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph; +use biodivine_lib_param_bn::BooleanNetwork; +use boolean_network_sketches::hctl_with_holes::reachability_two_holes::two_hole_analysis; +use boolean_network_sketches::hctl_with_holes::utils::encode_obs_weight_pairs; +#[allow(unused_imports)] +use boolean_network_sketches::hctl_with_holes::{MODEL_20V, MODEL_53V}; +use std::time::SystemTime; + +fn main() { + let use_large = false; + + let mut model = MODEL_20V; + let mut start_from_ones = true; + if use_large { + model = MODEL_53V; + start_from_ones = false; + } + + let bn = BooleanNetwork::try_from(model).unwrap(); + let props_vec: Vec<_> = bn + .variables() + .map(|v| bn.get_variable_name(v).clone()) + .collect(); + let stg = SymbolicAsyncGraph::new(&bn).unwrap(); + + let start_time = SystemTime::now(); + + // variants of each observation and their weights (one being ones, the other zeros) + let observation1_weight_pairs = encode_obs_weight_pairs(&props_vec, start_from_ones).unwrap(); + let observation2_weight_pairs = encode_obs_weight_pairs(&props_vec, !start_from_ones).unwrap(); + + let best_weight = two_hole_analysis( + &stg, + observation1_weight_pairs, + observation2_weight_pairs, + start_time, + ); + let time = start_time.elapsed().unwrap().as_millis(); + println!("BEST WEIGHT: {:?}", best_weight); + println!("TOTAL ELAPSED: {time}\n\n"); +} diff --git a/src/data_processing/data_encoding.rs b/src/data_processing/data_encoding.rs index 61d4842..5cc6b6e 100644 --- a/src/data_processing/data_encoding.rs +++ b/src/data_processing/data_encoding.rs @@ -4,7 +4,13 @@ use crate::data_processing::observations::*; /// Encode binarized observation with a formula depicting the corresponding state/sub-space. /// Using binarized values and proposition names, creates a conjunction of literals /// describing that observation. -pub fn encode_observation(observation: &Observation, prop_names: &[String]) -> String { +pub fn encode_observation( + observation: &Observation, + prop_names: &[String], +) -> Result { + if observation.num_values() != prop_names.len() { + return Err("Numbers of observation's values and propositions differs.".to_string()); + } let mut formula = String::new(); formula.push('('); @@ -23,25 +29,25 @@ pub fn encode_observation(observation: &Observation, prop_names: &[String]) -> S formula = formula.strip_suffix(" & ").unwrap().to_string(); } formula.push(')'); - formula + Ok(formula) } /// Encode several observation vectors with conjunction formulae, one by one. fn encode_multiple_observations( observations: &[Observation], prop_names: &[String], -) -> Vec { +) -> Result, String> { observations .iter() .map(|o| encode_observation(o, prop_names)) - .collect() + .collect::, String>>() } /// Encode (ordered) set of observations to a single HCTL formula. The particular formula /// template is chosen depending on the type of data. pub fn encode_observation_list_hctl(observation_list: ObservationList) -> Result { let encoded_observations = - encode_multiple_observations(&observation_list.observations, &observation_list.var_names); + encode_multiple_observations(&observation_list.observations, &observation_list.var_names)?; match observation_list.data_type { ObservationType::Attractor => Ok(mk_formula_attractor_set(encoded_observations)), ObservationType::FixedPoint => Ok(mk_formula_fixed_point_set(encoded_observations)), @@ -70,21 +76,31 @@ mod tests { let observation1 = Observation::try_from_str("001-1".to_string()).unwrap(); let encoded1 = "(~a & ~b & c & e)"; - assert_eq!(encode_observation(&observation1, &prop_names), encoded1); + assert_eq!( + encode_observation(&observation1, &prop_names).unwrap(), + encoded1 + ); let observation2 = Observation::try_from_str("001--".to_string()).unwrap(); let encoded2 = "(~a & ~b & c)"; - assert_eq!(encode_observation(&observation2, &prop_names), encoded2); + assert_eq!( + encode_observation(&observation2, &prop_names).unwrap(), + encoded2 + ); let observation3 = Observation::try_from_str("-----".to_string()).unwrap(); let encoded3 = "(true)"; - assert_eq!(encode_observation(&observation3, &prop_names), encoded3); + assert_eq!( + encode_observation(&observation3, &prop_names).unwrap(), + encoded3 + ); assert_eq!( encode_multiple_observations( &vec![observation1, observation2, observation3], &prop_names - ), + ) + .unwrap(), vec![encoded1, encoded2, encoded3] ); } diff --git a/src/data_processing/observations.rs b/src/data_processing/observations.rs index 5e3f512..af8aee3 100644 --- a/src/data_processing/observations.rs +++ b/src/data_processing/observations.rs @@ -27,12 +27,26 @@ pub struct Observation { } impl Observation { - /// Create observation object from string encoding its values. + /// Create `Observation` object from a vector with values. pub fn new(values: Vec) -> Self { Self { values } } - /// Create observation object from string encoding its values. + /// Create `Observation` encoding a vector of `n` True values. + pub fn new_full_true(n: usize) -> Self { + Self { + values: vec![VarValue::True; n], + } + } + + /// Create `Observation` encoding a vector of `n` False values. + pub fn new_full_false(n: usize) -> Self { + Self { + values: vec![VarValue::False; n], + } + } + + /// Create `Observation` object from string encoding its values. pub fn try_from_str(observation_string: String) -> Result { let mut observation_vec: Vec = Vec::new(); for c in observation_string.chars() { @@ -51,6 +65,29 @@ impl Observation { values: observation_vec, }) } + + pub fn num_values(&self) -> usize { + self.values.len() + } + + pub fn num_unspecified_values(&self) -> usize { + self.values.iter().filter(|&v| *v == VarValue::Any).count() + } + + pub fn num_specified_values(&self) -> usize { + self.values.iter().filter(|&v| *v != VarValue::Any).count() + } + + pub fn num_ones(&self) -> usize { + self.values.iter().filter(|&v| *v == VarValue::True).count() + } + + pub fn num_zeros(&self) -> usize { + self.values + .iter() + .filter(|&v| *v == VarValue::False) + .count() + } } impl fmt::Display for Observation { diff --git a/src/hctl_with_holes/mod.rs b/src/hctl_with_holes/mod.rs new file mode 100644 index 0000000..7a1d540 --- /dev/null +++ b/src/hctl_with_holes/mod.rs @@ -0,0 +1,268 @@ +pub mod reachability_one_hole; +pub mod reachability_two_holes; +pub mod utils; + +pub const MODEL_20V: &str = r" +$v_cMYC:(v_Akt1 | (v_ERa | v_MEK1)) +v_Akt1 -> v_cMYC +v_ERa -> v_cMYC +v_MEK1 -> v_cMYC +$v_p27:(!v_Akt1 & (!v_CDK2 & (!v_CDK4 & (v_ERa & !v_cMYC)))) +v_cMYC -| v_p27 +v_CDK2 -| v_p27 +v_Akt1 -| v_p27 +v_CDK4 -| v_p27 +v_ERa -> v_p27 +$v_CDK2:(v_CycE1 & (!v_p21 & !v_p27)) +v_p27 -| v_CDK2 +v_CycE1 -> v_CDK2 +v_p21 -| v_CDK2 +$v_Akt1:(v_ErbB1 | (v_ErbB1_2 | (v_ErbB1_3 | (v_ErbB2_3 | v_IGF1R)))) +v_ErbB2_3 -> v_Akt1 +v_ErbB1_2 -> v_Akt1 +v_ErbB1_3 -> v_Akt1 +v_IGF1R -> v_Akt1 +v_ErbB1 -> v_Akt1 +$v_CDK4:(v_CycD1 & (!v_p21 & !v_p27)) +v_CycD1 -> v_CDK4 +v_p27 -| v_CDK4 +v_p21 -| v_CDK4 +$v_ERa:(v_Akt1 | v_MEK1) +v_Akt1 -> v_ERa +v_MEK1 -> v_ERa +$v_CycE1:v_cMYC +v_cMYC -> v_CycE1 +$v_p21:(!v_Akt1 & (!v_CDK4 & (v_ERa & !v_cMYC))) +v_cMYC -| v_p21 +v_Akt1 -| v_p21 +v_CDK4 -| v_p21 +v_ERa -> v_p21 +$v_ErbB2_3:(v_ErbB2 & v_ErbB3) +v_ErbB3 -> v_ErbB2_3 +v_ErbB2 -> v_ErbB2_3 +$v_ErbB1_2:(v_ErbB1 & v_ErbB2) +v_ErbB1 -> v_ErbB1_2 +v_ErbB2 -> v_ErbB1_2 +$v_ErbB1_3:(v_ErbB1 & v_ErbB3) +v_ErbB3 -> v_ErbB1_3 +v_ErbB1 -> v_ErbB1_3 +$v_IGF1R:((v_Akt1 & !v_ErbB2_3) | (!v_Akt1 & (v_ERa & !v_ErbB2_3))) +v_ErbB2_3 -| v_IGF1R +v_Akt1 -> v_IGF1R +v_ERa -> v_IGF1R +$v_ErbB1:v_EGF +v_EGF -> v_ErbB1 +$v_MEK1:(v_ErbB1 | (v_ErbB1_2 | (v_ErbB1_3 | (v_ErbB2_3 | v_IGF1R)))) +v_ErbB2_3 -> v_MEK1 +v_ErbB1_2 -> v_MEK1 +v_ErbB1_3 -> v_MEK1 +v_ErbB1 -> v_MEK1 +v_IGF1R -> v_MEK1 +$v_EGF:true +$v_ErbB3:v_EGF +v_EGF -> v_ErbB3 +$v_CycD1:((v_Akt1 & (v_ERa & v_cMYC)) | (!v_Akt1 & (v_ERa & (v_MEK1 & v_cMYC)))) +v_cMYC -> v_CycD1 +v_Akt1 -> v_CycD1 +v_ERa -> v_CycD1 +v_MEK1 -> v_CycD1 +$v_ErbB2:v_EGF +v_EGF -> v_ErbB2 +$v_CDK6:v_CycD1 +v_CycD1 -> v_CDK6 +$v_pRB:(v_CDK4 & v_CDK6) +v_CDK6 -> v_pRB +v_CDK2 ->? v_pRB +v_CDK4 -> v_pRB +"; + +pub const MODEL_53V: &str = r" +$v_DCII_Bacterium:v_DCI_Bacterium +v_DCI_Bacterium -> v_DCII_Bacterium +$v_T0:(v_DCII_Bacterium | v_DCII_TRetortaeformis) +v_DCII_Bacterium -> v_T0 +v_DCII_TRetortaeformis -> v_T0 +$v_DCII_TRetortaeformis:v_DCI_TRetortaeformis +v_DCI_TRetortaeformis -> v_DCII_TRetortaeformis +$v_Th2II_TRetortaeformis:(v_DCII_TRetortaeformis & (!v_IL12II & v_T0)) +v_T0 -> v_Th2II_TRetortaeformis +v_DCII_TRetortaeformis -> v_Th2II_TRetortaeformis +v_IL12II -| v_Th2II_TRetortaeformis +$v_IL12II:((v_DCII_Bacterium & (!v_IL4II & v_T0)) | (!v_DCII_Bacterium & (v_DCII_TRetortaeformis & (!v_IL4II & v_T0)))) +v_T0 -> v_IL12II +v_DCII_Bacterium -> v_IL12II +v_DCII_TRetortaeformis -> v_IL12II +v_IL4II -| v_IL12II +$v_DCI_TRetortaeformis:v_PIC +v_PIC -> v_DCI_TRetortaeformis +$v_IL4II:((v_DCII_Bacterium & (v_EL2 | (!v_IFNgI & (!v_IL12II & (v_T0 | (v_Th2II_Bacterium | v_Th2II_TRetortaeformis)))))) | (!v_DCII_Bacterium & ((v_DCII_TRetortaeformis & (v_EL2 | (!v_IFNgI & (!v_IL12II & (v_T0 | (v_Th2II_Bacterium | v_Th2II_TRetortaeformis)))))) | (!v_DCII_TRetortaeformis & (v_EL2 | (!v_IFNgI & (!v_IL12II & (v_Th2II_Bacterium | v_Th2II_TRetortaeformis)))))))) +v_DCII_Bacterium -> v_IL4II +v_T0 -> v_IL4II +v_Th2II_Bacterium -> v_IL4II +v_DCII_TRetortaeformis -> v_IL4II +v_Th2II_TRetortaeformis -> v_IL4II +v_EL2 -> v_IL4II +v_IL12II -| v_IL4II +v_IFNgI -| v_IL4II +$v_Th2II_Bacterium:(v_DCII_Bacterium & (!v_IL12II & v_T0)) +v_DCII_Bacterium -> v_Th2II_Bacterium +v_T0 -> v_Th2II_Bacterium +v_IL12II -| v_Th2II_Bacterium +$v_EL2:((v_IL13 & v_IL5) | (!v_IL13 & (v_IL5 & v_IgE))) +v_IL13 -> v_EL2 +v_IL5 -> v_EL2 +v_IgE -> v_EL2 +$v_IFNgI:(v_DCI_TRetortaeformis | (v_IFNg_Bacterium | v_Th1I_TRetortaeformis)) +v_Th1I_TRetortaeformis -> v_IFNgI +v_IFNg_Bacterium -> v_IFNgI +v_DCI_TRetortaeformis -> v_IFNgI +$v_IL4I:v_IL4II +v_IL4II -> v_IL4I +$v_BC_TRetortaeformis:(v_BC_TRetortaeformis | v_T0) +v_T0 -> v_BC_TRetortaeformis +v_BC_TRetortaeformis -> v_BC_TRetortaeformis +$v_IS:false +$v_IgA_TRetortaeformis:(v_BC_TRetortaeformis & v_IS) +v_IS -> v_IgA_TRetortaeformis +v_BC_TRetortaeformis -> v_IgA_TRetortaeformis +$v_AP:((v_AgAb_Bacterium & (v_Bb & (v_MPI_Bacterium & v_Th1I_Bacterium))) | (!v_AgAb_Bacterium & (v_Bb & (v_Cb & (v_IgG_Bacterium & (v_MPI_Bacterium & v_Th1I_Bacterium)))))) +v_Cb -> v_AP +v_Th1I_Bacterium -> v_AP +v_IgG_Bacterium -> v_AP +v_AgAb_Bacterium -> v_AP +v_MPI_Bacterium -> v_AP +v_Bb -> v_AP +$v_PH:(v_AP & v_Bb) +v_AP -> v_PH +v_Bb -> v_PH +$v_Bb:(v_Bb & !v_PH) +v_PH -| v_Bb +v_Bb -> v_Bb +$v_EC_Bacterium:v_Bb +v_Bb -> v_EC_Bacterium +$v_Th2I_Bacterium:v_Th2II_Bacterium +v_Th2II_Bacterium -> v_Th2I_Bacterium +$v_IL13:((v_EL & (v_EL2 | (v_IS | (v_Th2I_Bacterium | v_Th2I_TRetortaeformis)))) | (!v_EL & (v_EL2 | (v_Th2I_Bacterium | v_Th2I_TRetortaeformis)))) +v_Th2I_Bacterium -> v_IL13 +v_Th2I_TRetortaeformis -> v_IL13 +v_EL2 -> v_IL13 +v_IS -> v_IL13 +v_EL -> v_IL13 +$v_Th2I_TRetortaeformis:v_Th2II_TRetortaeformis +v_Th2II_TRetortaeformis -> v_Th2I_TRetortaeformis +$v_EL:(!v_EL2 & v_IS) +v_EL2 -| v_EL +v_IS -> v_EL +$v_IL5:(v_EL2 | v_Th2II_TRetortaeformis) +v_Th2II_TRetortaeformis -> v_IL5 +v_EL2 -> v_IL5 +$v_IgE:(v_BC_TRetortaeformis & (v_IL13 | v_IL4II)) +v_IL13 -> v_IgE +v_BC_TRetortaeformis -> v_IgE +v_IL4II -> v_IgE +$v_PIC:((v_AD & (!v_IL10I & !v_IgA_TRetortaeformis)) | (!v_AD & ((v_AP & (!v_IL10I & !v_IgA_TRetortaeformis)) | (!v_AP & ((v_EC_Bacterium & (!v_IL10I & !v_IgA_TRetortaeformis)) | (!v_EC_Bacterium & (v_EC_TRetortaeformis & (!v_IL10I & !v_IgA_TRetortaeformis)))))))) +v_AD -> v_PIC +v_IL10I -| v_PIC +v_IgA_TRetortaeformis -| v_PIC +v_EC_Bacterium -> v_PIC +v_AP -> v_PIC +v_EC_TRetortaeformis -> v_PIC +$v_DCI_Bacterium:(v_Bb & (v_IFNg_Bacterium | v_PIC)) +v_IFNg_Bacterium -> v_DCI_Bacterium +v_PIC -> v_DCI_Bacterium +v_Bb -> v_DCI_Bacterium +$v_TrII:(v_DCII_Bacterium & (v_T0 & v_TTSSII)) +v_DCII_Bacterium -> v_TrII +v_T0 -> v_TrII +v_TTSSII -> v_TrII +$v_TTSSII:v_TTSSI +v_TTSSI -> v_TTSSII +$v_TTSSI:(v_Bb & (!v_IgA_Bacterium & !v_IgG_Bacterium)) +v_IgG_Bacterium -| v_TTSSI +v_IgA_Bacterium -| v_TTSSI +v_Bb -> v_TTSSI +$v_IFNg_Bacterium:(v_DCI_Bacterium | ((v_IL10I_Bacterium & v_MPI_Bacterium) | (!v_IL10I_Bacterium & ((v_IL4I & v_MPI_Bacterium) | (!v_IL4I & (v_MPI_Bacterium | v_Th1I_Bacterium)))))) +v_IL4I -| v_IFNg_Bacterium +v_DCI_Bacterium -> v_IFNg_Bacterium +v_Th1I_Bacterium -> v_IFNg_Bacterium +v_IL10I_Bacterium -| v_IFNg_Bacterium +v_MPI_Bacterium -> v_IFNg_Bacterium +$v_Th1I_Bacterium:v_Th1II_Bacterium +v_Th1II_Bacterium -> v_Th1I_Bacterium +$v_IL10I_Bacterium:(v_MPI_Bacterium | ((v_TTSSI & (v_Th2I_Bacterium | v_TrI_Bacterium)) | (!v_TTSSI & v_TrI_Bacterium))) +v_Th2I_Bacterium -> v_IL10I_Bacterium +v_TrI_Bacterium -> v_IL10I_Bacterium +v_MPI_Bacterium -> v_IL10I_Bacterium +v_TTSSI -> v_IL10I_Bacterium +$v_MPI_Bacterium:(v_Bb & (v_IFNg_Bacterium | v_PIC)) +v_IFNg_Bacterium -> v_MPI_Bacterium +v_PIC -> v_MPI_Bacterium +v_Bb -> v_MPI_Bacterium +$v_AD:((v_AD & (v_IgG & (!v_MPI_Bacterium & !v_NE_TRetortaeformis))) | (!v_AD & (v_IS & (v_IgG & (!v_MPI_Bacterium & !v_NE_TRetortaeformis))))) +v_AD -> v_AD +v_MPI_Bacterium -| v_AD +v_IS -> v_AD +v_IgG -> v_AD +v_NE_TRetortaeformis -| v_AD +$v_IL10I:(v_IL10I_Bacterium | v_Th2I_TRetortaeformis) +v_IL10I_Bacterium -> v_IL10I +v_Th2I_TRetortaeformis -> v_IL10I +$v_EC_TRetortaeformis:(v_AD | v_IS) +v_AD -> v_EC_TRetortaeformis +v_IS -> v_EC_TRetortaeformis +$v_IFNgII:(v_IFNgI | v_IFNg_Bacterium) +v_IFNg_Bacterium -> v_IFNgII +v_IFNgI -> v_IFNgII +$v_Th1I_TRetortaeformis:v_Th1II_TRetortaeformis +v_Th1II_TRetortaeformis -> v_Th1I_TRetortaeformis +$v_IgG_Bacterium:(v_BC_Bacterium | v_IgG_Bacterium) +v_IgG_Bacterium -> v_IgG_Bacterium +v_BC_Bacterium -> v_IgG_Bacterium +$v_BC_Bacterium:(v_BC_Bacterium | v_T0) +v_T0 -> v_BC_Bacterium +v_BC_Bacterium -> v_BC_Bacterium +$v_IgG:v_BC_TRetortaeformis +v_BC_TRetortaeformis -> v_IgG +$v_NE_TRetortaeformis:((v_AD & ((v_IFNgI & ((v_IL10I & v_PIC) | (!v_IL10I & (!v_IL4I | v_PIC)))) | (!v_IFNgI & v_PIC))) | (!v_AD & (v_IFNgI & (!v_IL10I & !v_IL4I)))) +v_AD -> v_NE_TRetortaeformis +v_IL4I -| v_NE_TRetortaeformis +v_IL10I -| v_NE_TRetortaeformis +v_IFNgI -> v_NE_TRetortaeformis +v_PIC -> v_NE_TRetortaeformis +$v_NE_Bacterium:v_PIC +v_PIC -> v_NE_Bacterium +$v_DP:(v_NE_Bacterium & v_TTSSI) +v_NE_Bacterium -> v_DP +v_TTSSI -> v_DP +$v_IgA_Bacterium:((v_BC_Bacterium & v_Bb) | (!v_BC_Bacterium & (v_Bb & v_IgA_Bacterium))) +v_BC_Bacterium -> v_IgA_Bacterium +v_IgA_Bacterium -> v_IgA_Bacterium +v_Bb -> v_IgA_Bacterium +$v_Th1II_Bacterium:(v_DCII_Bacterium & (v_IL12II & v_T0)) +v_DCII_Bacterium -> v_Th1II_Bacterium +v_T0 -> v_Th1II_Bacterium +v_IL12II -> v_Th1II_Bacterium +$v_TrI_Bacterium:v_TrII +v_TrII -> v_TrI_Bacterium +$v_Cb:((v_AgAb_Bacterium & ((v_Bb & (v_IgG_Bacterium | !v_Oag)) | (!v_Bb & v_IgG_Bacterium))) | (!v_AgAb_Bacterium & (v_Bb & !v_Oag))) +v_Oag -| v_Cb +v_IgG_Bacterium -> v_Cb +v_AgAb_Bacterium -> v_Cb +v_Bb -> v_Cb +$v_AgAb_Bacterium:(v_Bb & (v_IgA_Bacterium | v_IgG_Bacterium)) +v_IgG_Bacterium -> v_AgAb_Bacterium +v_IgA_Bacterium -> v_AgAb_Bacterium +v_Bb -> v_AgAb_Bacterium +$v_Oag:v_Bb +v_Bb -> v_Oag +$v_Th1II_TRetortaeformis:(v_DCII_TRetortaeformis & (v_IL12II & v_T0)) +v_T0 -> v_Th1II_TRetortaeformis +v_DCII_TRetortaeformis -> v_Th1II_TRetortaeformis +v_IL12II -> v_Th1II_TRetortaeformis +$v_TEL:(v_EL | v_EL2) +v_EL2 -> v_TEL +v_EL -> v_TEL +$v_TNE:(v_NE_Bacterium | v_NE_TRetortaeformis) +v_NE_Bacterium -> v_TNE +v_NE_TRetortaeformis -> v_TNE +"; diff --git a/src/hctl_with_holes/reachability_one_hole.rs b/src/hctl_with_holes/reachability_one_hole.rs new file mode 100644 index 0000000..e2eace1 --- /dev/null +++ b/src/hctl_with_holes/reachability_one_hole.rs @@ -0,0 +1,216 @@ +use crate::hctl_with_holes::utils::{fwd_saturated, update_weight}; +use biodivine_hctl_model_checker::model_checking::{ + model_check_extended_formula, model_check_extended_formula_dirty, model_check_formula, + model_check_formula_dirty, +}; +use biodivine_lib_param_bn::biodivine_std::traits::Set; +use biodivine_lib_param_bn::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph}; +use std::collections::HashMap; +use std::time::SystemTime; + +/// Naive version of the reachability from given (fixed) initial observation to a second "uncertain" one. +/// Simply runs model checking for all variants of the second observation individually. +pub fn naive_eval( + stg: &SymbolicAsyncGraph, + time_0: SystemTime, + init_observation: &GraphColoredVertices, + second_observation_weight_pairs: Vec<(String, usize)>, +) -> Option { + let mut best_weight = None; + let formula_with_holes = "3{x}: @{x}: (%initial% & EF %second%)".to_string(); + + for (observation_formula, weight) in second_observation_weight_pairs { + let second_observation_set = model_check_formula_dirty(&observation_formula, stg).unwrap(); + let context = HashMap::from([ + ("initial".to_string(), init_observation.clone()), + ("second".to_string(), second_observation_set.clone()), + ]); + + // check if the result is trivial (i.e., some states of the second observation's set are directly + // contained in the set of the initial observation) + if !init_observation + .intersect(&second_observation_set) + .is_empty() + { + println!("trivial success - w: {weight}"); + best_weight = update_weight(best_weight, weight); + continue; + } + + let result = model_check_extended_formula(&formula_with_holes, stg, &context).unwrap(); + let time = time_0.elapsed().unwrap().as_millis(); + if !result.is_empty() { + let cardinality = result.exact_cardinality(); + best_weight = update_weight(best_weight, weight); + println!("success - w: {weight}, time: {time}ms, states: {cardinality}"); + } else { + println!("--- w: {weight}, time: {time}ms"); + } + } + best_weight +} + +/// Slightly better version of reachability from given (fixed) initial observation to a second "uncertain" one. +/// Evaluation goes in the order from the most specific (smallest) set to the most general (largest) set. +/// In each computation,instead of computing full EF again, EF result from the previous step is reused. +pub fn specific_to_general_eval( + stg: &SymbolicAsyncGraph, + time_0: SystemTime, + init_observation: &GraphColoredVertices, + second_observation_weight_pairs: Vec<(String, usize)>, +) -> Option { + let mut best_weight = None; + let formula_with_holes = "3{x}: @{x}: (%initial% & %ef_observation%)".to_string(); + + // remember EF result from the last step (initially empty) + let mut last_ef_result = stg.empty_colored_vertices().clone(); + for (observation_formula, weight) in second_observation_weight_pairs { + let second_observation = model_check_formula_dirty(&observation_formula, stg).unwrap(); + + // check if the result is trivial (i.e., some states of the second observation's set are directly + // contained in the set of the initial observation) + if !init_observation.intersect(&second_observation).is_empty() { + best_weight = update_weight(best_weight, weight); + println!("trivial success - w: {weight}"); + continue; + } + + let extended_observation = second_observation.union(&last_ef_result); + last_ef_result = model_check_extended_formula_dirty( + "EF %extended_observation%", + stg, + &HashMap::from([("extended_observation".to_string(), extended_observation)]), + ) + .unwrap(); + + let context = HashMap::from([ + ("initial".to_string(), init_observation.clone()), + ("ef_observation".to_string(), last_ef_result.clone()), + ]); + + let result = model_check_extended_formula(&formula_with_holes, stg, &context).unwrap(); + let time = time_0.elapsed().unwrap().as_millis(); + if !result.is_empty() { + let cardinality = result.exact_cardinality(); + best_weight = update_weight(best_weight, weight); + println!("success - w: {weight}, time: {time}ms, states: {cardinality}"); + } else { + let last_cardinality = last_ef_result.exact_cardinality(); + println!("--- w: {weight}, time: {time}ms, last_ef: {last_cardinality}"); + } + } + best_weight +} + +/// Better version of reachability from given (fixed) initial observation to a second "uncertain" one. +/// Evaluation goes in the order from the most general (largest) set to the most specific (smallest) set. +/// Starts by pre-computing forward reachability (once) from the initial state, then just iteratively +/// intersects that with the second observation's sets. +pub fn general_to_specific_precomp_eval( + stg: &SymbolicAsyncGraph, + time_0: SystemTime, + init_observation: &GraphColoredVertices, + mut second_observation_weight_pairs: Vec<(String, usize)>, +) -> Option { + let mut best_weight = None; + // fwd reachable set from the initial observation + let fwd_set_from_initial = fwd_saturated(stg, init_observation); + + second_observation_weight_pairs.reverse(); + + for (observation_formula, weight) in second_observation_weight_pairs { + let second_observation = model_check_formula(&observation_formula, stg).unwrap(); + + // check if the result is trivial (i.e., some states of the second observation's set are directly + // contained in the set of the initial observation) + if !init_observation.intersect(&second_observation).is_empty() { + best_weight = update_weight(best_weight, weight); + println!("trivial success - w: {weight}"); + continue; + } + + let result = fwd_set_from_initial.intersect(&second_observation); + let time = time_0.elapsed().unwrap().as_millis(); + if !result.is_empty() { + let cardinality = result.exact_cardinality(); + best_weight = update_weight(best_weight, weight); + println!("success - w: {weight}, time: {time}ms, states: {cardinality}"); + } else { + println!("--- w: {weight}, time: {time}ms"); + } + } + best_weight +} + +/// Better version of reachability from given (fixed) initial observation to a second "uncertain" one. +/// Evaluation in the order from the most general (largest) set to the most specific (smallest) set. +/// Starts computing forward reachability from the initial state, and during this computation +/// checks if then the current reach set intersects with the second observation's sets. +/// +/// This means that we dont have to compute full reachability set to know the result, if lucky. +pub fn general_to_specific_on_the_fly_eval( + stg: &SymbolicAsyncGraph, + time_0: SystemTime, + init_observation: &GraphColoredVertices, + mut second_observation_weight_pairs: Vec<(String, usize)>, +) -> Option { + let mut best_weight = None; + + let mut done = false; + let mut result = init_observation.clone(); + let (mut current_formula, mut current_w) = second_observation_weight_pairs.pop().unwrap(); + let mut current_second_observation = model_check_formula(¤t_formula, stg).unwrap(); + + // saturated reachability + while !done { + done = true; + for var in stg.as_network().unwrap().variables() { + let update = &stg.var_post(var, &result).minus(&result); + if !update.is_empty() { + result = result.union(update); + + // if new reach states found, check for intersection with current set for the second observation + if !result.intersect(¤t_second_observation).is_empty() { + let time = time_0.elapsed().unwrap().as_millis(); + best_weight = update_weight(best_weight, current_w); + + // check if the result is trivial (i.e., some states of the second observation's set are directly + // contained in the set of the initial observation) or not + if !init_observation + .intersect(¤t_second_observation) + .is_empty() + { + println!("trivial success - w: {current_w}, time: {time}ms"); + } else { + println!("success - w: {current_w}, time: {time}ms"); + } + + // update the current set for the second observation for the next one + let current_obs_pair = second_observation_weight_pairs.pop().unwrap(); + current_formula = current_obs_pair.0; + current_w = current_obs_pair.1; + current_second_observation = + model_check_formula(¤t_formula, stg).unwrap(); + } + done = false; + break; + } + } + } + + // check remaining observation's sets (there might have been more than one hits at the last iteration) + second_observation_weight_pairs.reverse(); + for (observation_formula, weight) in second_observation_weight_pairs { + let second_observation = model_check_formula(&observation_formula, stg).unwrap(); + + let inter = result.intersect(&second_observation); + let time = time_0.elapsed().unwrap().as_millis(); + if !inter.is_empty() { + best_weight = update_weight(best_weight, weight); + println!("success - w: {weight}, time: {time}ms"); + } else { + println!("--- w: {weight}, time: {time}ms"); + } + } + best_weight +} diff --git a/src/hctl_with_holes/reachability_two_holes.rs b/src/hctl_with_holes/reachability_two_holes.rs new file mode 100644 index 0000000..b615618 --- /dev/null +++ b/src/hctl_with_holes/reachability_two_holes.rs @@ -0,0 +1,57 @@ +use crate::hctl_with_holes::utils::{fwd_saturated, update_weight}; +use biodivine_hctl_model_checker::model_checking::model_check_formula; +use biodivine_lib_param_bn::biodivine_std::traits::Set; +use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph; +use std::time::SystemTime; + +/// Reachability from an "uncertain" initial observation to another "uncertain" one. +/// Currently, we consider all combinations of the following pattern: +/// 1111 -> 0000, 111* -> 0000, 11** -> 0000, ..., 1111 -> 000*, 1111 -> 00**, ..., 111* -> 000*, ... +pub fn two_hole_analysis( + stg: &SymbolicAsyncGraph, + observation1_weight_pairs: Vec<(String, usize)>, + observation2_weight_pairs: Vec<(String, usize)>, + time_0: SystemTime, +) -> Option { + let mut best_weight = None; + + // the HCTL formula for this would be like "3{x}: @{x}: (%observation1% & EF %observation2%)" + + // Start from the most specific set for the observation 1. Each subsequent set is thus a superset of the set + // from the previous step. This enables us to reuse results of forward reachability from previous steps. + let mut last_fwd_result = stg.empty_colored_vertices().clone(); + for (observation1_formula, weight1) in observation1_weight_pairs { + // Compute saturated reachability from the current set for observation 1. + let raw_observation1 = model_check_formula(&observation1_formula, stg).unwrap(); + let fwd_seed_set = raw_observation1.union(&last_fwd_result); + let fwd_set_from_observation1 = fwd_saturated(stg, &fwd_seed_set); + last_fwd_result = fwd_set_from_observation1.clone(); + + // For all variants of sets for observation 2, just check if there is an intersection with the reachable + // states from observation 1. + for (observation2_formula, weight2) in &observation2_weight_pairs { + let raw_observation2 = model_check_formula(observation2_formula, stg).unwrap(); + + // check if the result is trivial (i.e., some states of the second observation's set are directly + // contained in the set of the initial observation) + if !raw_observation1.intersect(&raw_observation2).is_empty() { + let weight_sum = weight1 + weight2; + best_weight = update_weight(best_weight, weight_sum); + println!("trivial success - w1: {weight1}, w2: {weight2}, w1+w2: {weight_sum}"); + continue; + } + + let result = fwd_set_from_observation1.intersect(&raw_observation2); + let time = time_0.elapsed().unwrap().as_millis(); + if !result.is_empty() { + let cardinality = result.exact_cardinality(); + let weight_sum = weight1 + weight2; + best_weight = update_weight(best_weight, weight_sum); + println!("success - w1: {weight1}, w2: {weight2}, w1+w2: {weight_sum} time: {time}, states: {cardinality}"); + } else { + println!("--- w1: {weight1}, w2: {weight2}, time: {time}"); + } + } + } + best_weight +} diff --git a/src/hctl_with_holes/utils.rs b/src/hctl_with_holes/utils.rs new file mode 100644 index 0000000..8919cd4 --- /dev/null +++ b/src/hctl_with_holes/utils.rs @@ -0,0 +1,60 @@ +use crate::data_processing::data_encoding::encode_observation; +use crate::data_processing::observations::Observation; +use biodivine_lib_param_bn::biodivine_std::traits::Set; +use biodivine_lib_param_bn::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph}; + +/// Compute forward reachable states with saturation-based algorithm. +pub fn fwd_saturated( + graph: &SymbolicAsyncGraph, + initial: &GraphColoredVertices, +) -> GraphColoredVertices { + let mut result = initial.clone(); + let mut done = false; + while !done { + done = true; + for var in graph.as_network().unwrap().variables() { + let update = &graph.var_post(var, &result).minus(&result); + if !update.is_empty() { + result = result.union(update); + done = false; + break; + } + } + } + result +} + +/// Generate list of formulae encoding a sequence of partial observations starting at a full state with +/// all propositions being true (or false), and each successive observation loosens up one variable. +/// For example, for 4 propositions, the sequence is following: +/// 1111, 111*, 11**, 1***, **** +/// Each observation is given a weight equal to the number of free propositions (full state has weight 0). +pub fn encode_obs_weight_pairs( + prop_names: &[String], + use_ones: bool, +) -> Result, String> { + let n = prop_names.len(); + let mut observation_weight_pairs = Vec::new(); + for i in 0..n + 1 { + let mut obs_str = String::new(); + for _ in 0..n - i { + obs_str.push(if use_ones { '1' } else { '0' }); + } + for _ in n - i..n { + obs_str.push('-'); + } + let observation = Observation::try_from_str(obs_str)?; + let observation_formula = encode_observation(&observation, prop_names)?; + observation_weight_pairs.push((observation_formula, i)); + } + Ok(observation_weight_pairs) +} + +/// Update the (potentially None-valued) weight to a new value, if it is better (smaller). +/// If the other weight is the same or worse, the original value is returned. +pub fn update_weight(original: Option, new: usize) -> Option { + match original { + Some(val) if { val < new } => original, + _ => Some(new), + } +} diff --git a/src/inference_attractor_data.rs b/src/inference_attractor_data.rs index a98ee31..886fed9 100644 --- a/src/inference_attractor_data.rs +++ b/src/inference_attractor_data.rs @@ -50,7 +50,7 @@ pub fn perform_inference_with_attractors_specific( }; // compute satisfying colours - inferred_colors = model_check_formula_unsafe_ex(formula, &graph) + inferred_colors = model_check_formula_unsafe_ex(&formula, &graph) .unwrap() .colors(); @@ -77,7 +77,7 @@ pub fn perform_inference_with_attractors_specific( } else { mk_formula_forbid_other_attractors(attr_set) }; - inferred_colors = model_check_formula_unsafe_ex(formula, &graph) + inferred_colors = model_check_formula_unsafe_ex(&formula, &graph) .unwrap() .colors(); } diff --git a/src/lib.rs b/src/lib.rs index 5fee5ef..25b6713 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,5 +2,6 @@ //! with case studies present as binaries. pub mod data_processing; +pub mod hctl_with_holes; pub mod inference_attractor_data; pub mod utils; diff --git a/src/utils.rs b/src/utils.rs index a79000b..c55265b 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -17,7 +17,7 @@ use std::collections::HashMap; /// Apply properties (constraints) given by HCTL `formulae` on the graph's colors. /// Returns a graph with colour space restricted only to the suitable colors. pub fn apply_constraints_and_restrict( - formulae: Vec, + formulae: Vec<&str>, mut graph: SymbolicAsyncGraph, message: &str, ) -> SymbolicAsyncGraph {