Skip to content

Commit

Permalink
Start working on analysis of HCTL with placeholders.
Browse files Browse the repository at this point in the history
  • Loading branch information
ondrej33 committed Feb 26, 2024
1 parent 65668a5 commit 1bdb60d
Show file tree
Hide file tree
Showing 14 changed files with 839 additions and 35 deletions.
19 changes: 12 additions & 7 deletions Cargo.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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"
23 changes: 13 additions & 10 deletions src/bin/case_study_tlgl.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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<String> = 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
Expand All @@ -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",
Expand All @@ -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());
Expand All @@ -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",
Expand Down Expand Up @@ -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(),
Expand All @@ -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()
Expand Down Expand Up @@ -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.);
}
}
8 changes: 4 additions & 4 deletions src/bin/small_example.rs
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,15 @@ 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(),
);

// 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!(
Expand Down Expand Up @@ -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.);
}
}
100 changes: 100 additions & 0 deletions src/bin/test_hctl_with_1_hole.rs
Original file line number Diff line number Diff line change
@@ -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);
}
41 changes: 41 additions & 0 deletions src/bin/test_hctl_with_2_holes.rs
Original file line number Diff line number Diff line change
@@ -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");
}
34 changes: 25 additions & 9 deletions src/data_processing/data_encoding.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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<String, String> {
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('(');

Expand All @@ -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<String> {
) -> Result<Vec<String>, String> {
observations
.iter()
.map(|o| encode_observation(o, prop_names))
.collect()
.collect::<Result<Vec<String>, 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<String, String> {
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)),
Expand Down Expand Up @@ -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]
);
}
Expand Down
Loading

0 comments on commit 1bdb60d

Please sign in to comment.