From a6eaa95690f3bc9a0848f647f88a5e28bef0586a Mon Sep 17 00:00:00 2001 From: Samuel Pastva Date: Sat, 12 Dec 2020 00:15:56 +0100 Subject: [PATCH] Improvements to algorithm and attractor visualisation. --- Cargo.toml | 2 +- src/bin/benchmark_validator.rs | 87 ++++++--- src/bin/experiment.rs | 21 ++- src/main.rs | 263 ++++++++++++++++----------- src/scc/algo_effectively_constant.rs | 223 ++++++++++++++++++++--- src/scc/algo_par_utils.rs | 13 +- src/scc/algo_symbolic_components.rs | 240 +++++++++++++++++++++--- src/scc/algo_symbolic_reach.rs | 48 +++-- src/scc/impl_classifier.rs | 94 +++++++--- src/scc/mod.rs | 7 +- 10 files changed, 755 insertions(+), 243 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index ccfb12b..573a0d7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -17,7 +17,7 @@ debug = true [dependencies] biodivine-lib-std = { git = "https://github.com/sybila/biodivine-lib-std.git" } -biodivine-lib-param-bn = { git = "https://github.com/sybila/biodivine-lib-param-bn.git", rev = "62db744ff0b217522081db4f9edacceee698daae" } +biodivine-lib-param-bn = { git = "https://github.com/sybila/biodivine-lib-param-bn.git", rev = "f9dd65e584f81debe9e17ea2c6a439070fdaf6a9" } biodivine-lib-bdd = "0.1.0" rocket = "0.4.6" rayon = "1.5.0" diff --git a/src/bin/benchmark_validator.rs b/src/bin/benchmark_validator.rs index d640e51..fdc3db3 100644 --- a/src/bin/benchmark_validator.rs +++ b/src/bin/benchmark_validator.rs @@ -1,7 +1,7 @@ -use biodivine_lib_param_bn::{BooleanNetwork, FnUpdate}; -use std::convert::TryFrom; use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph; +use biodivine_lib_param_bn::{BooleanNetwork, FnUpdate}; use std::collections::HashSet; +use std::convert::TryFrom; fn main() { //let args: Vec = std::env::args().collect(); @@ -24,8 +24,11 @@ fn main() { let sbml_model_path = bench_dir.path().join("model.sbml"); let sbml_model = if !sbml_model_path.exists() { errors += 1; - eprintln!("ERROR: Missing model.sbml in {}.", bench_dir.path().display()); - continue + eprintln!( + "ERROR: Missing model.sbml in {}.", + bench_dir.path().display() + ); + continue; } else { // Check that the sbml model is readable: let model_string = std::fs::read_to_string(sbml_model_path).unwrap(); @@ -33,17 +36,23 @@ fn main() { match model { Err(err) => { errors += 1; - eprintln!("ERROR: Invalid SBML model in {}.", bench_dir.path().display()); + eprintln!( + "ERROR: Invalid SBML model in {}.", + bench_dir.path().display() + ); eprintln!("\t\t{}", err); - continue + continue; } - Ok((model, _)) => model + Ok((model, _)) => model, } }; let aeon_model_path = bench_dir.path().join("model.aeon"); if !aeon_model_path.exists() { errors += 1; - eprintln!("ERROR: Missing model.aeon in {}.", bench_dir.path().display()); + eprintln!( + "ERROR: Missing model.aeon in {}.", + bench_dir.path().display() + ); } else { // Check that the aeon model is valid: let model_string = std::fs::read_to_string(aeon_model_path.clone()).unwrap(); @@ -51,27 +60,48 @@ fn main() { match model { Ok(mut model) => { // Check that basic properties match SBML model. But note that variables can be re-ordered... - let mut models_match = model.graph().num_vars() == sbml_model.graph().num_vars(); + let mut models_match = + model.graph().num_vars() == sbml_model.graph().num_vars(); if model.graph().num_vars() != sbml_model.graph().num_vars() { - eprintln!("{} != {}", model.graph().num_vars(), sbml_model.graph().num_vars()); + eprintln!( + "{} != {}", + model.graph().num_vars(), + sbml_model.graph().num_vars() + ); } for v in model.graph().variable_ids() { - let regulators_in_model: HashSet = model.graph().regulators(v).into_iter() + let regulators_in_model: HashSet = model + .graph() + .regulators(v) + .into_iter() .map(|r| model.graph().get_variable(r).get_name().clone()) .collect(); - let regulators_in_sbml_model: HashSet = sbml_model.graph().regulators( - sbml_model.graph().find_variable(model.graph().get_variable(v).get_name()).unwrap() - ).into_iter() + let regulators_in_sbml_model: HashSet = sbml_model + .graph() + .regulators( + sbml_model + .graph() + .find_variable(model.graph().get_variable(v).get_name()) + .unwrap(), + ) + .into_iter() .map(|r| sbml_model.graph().get_variable(r).get_name().clone()) .collect(); if regulators_in_model != regulators_in_sbml_model { - eprintln!("{:?} != {:?}", regulators_in_model, regulators_in_sbml_model); + eprintln!( + "{:?} != {:?}", + regulators_in_model, regulators_in_sbml_model + ); } - models_match = models_match && regulators_in_model == regulators_in_sbml_model; + models_match = + models_match && regulators_in_model == regulators_in_sbml_model; } if !models_match { errors += 1; - eprintln!("ERROR: SBML and AEON model are different in {}.", bench_dir.path().display()); + eprintln!( + "ERROR: SBML and AEON model are different in {}.", + bench_dir.path().display() + ); } // Check that all update functions are set (for non-parametrized model anyway). let mut model_ok = true; @@ -84,7 +114,10 @@ fn main() { } if !model_ok { errors += 1; - eprintln!("ERROR: Model in {} contains unconstrained variables.", bench_dir.path().display()); + eprintln!( + "ERROR: Model in {} contains unconstrained variables.", + bench_dir.path().display() + ); if auto_fix { std::fs::write(aeon_model_path, model.to_string()).unwrap(); } else { @@ -97,12 +130,19 @@ fn main() { Ok(graph) => { if graph.unit_colors().cardinality() != 1.0 { errors += 1; - eprintln!("ERROR: Default model has {} colors in {}.", graph.unit_colors().cardinality(), bench_dir.path().display()); + eprintln!( + "ERROR: Default model has {} colors in {}.", + graph.unit_colors().cardinality(), + bench_dir.path().display() + ); } } Err(err) => { errors += 1; - eprintln!("ERROR: Cannot build graph from model in {}.", bench_dir.path().display()); + eprintln!( + "ERROR: Cannot build graph from model in {}.", + bench_dir.path().display() + ); eprintln!("{}", err); } } @@ -110,7 +150,10 @@ fn main() { } Err(err) => { errors += 1; - eprintln!("ERROR: Invalid AEON model in {}.", bench_dir.path().display()); + eprintln!( + "ERROR: Invalid AEON model in {}.", + bench_dir.path().display() + ); eprintln!("\t\t{}", err); } } @@ -118,4 +161,4 @@ fn main() { println!("OK: {}", bench_dir.path().display()); } println!("TOTAL ERRORS: {}", errors); -} \ No newline at end of file +} diff --git a/src/bin/experiment.rs b/src/bin/experiment.rs index 13bf137..b226567 100644 --- a/src/bin/experiment.rs +++ b/src/bin/experiment.rs @@ -1,4 +1,4 @@ -use biodivine_aeon_server::scc::algo_symbolic_components::components; +use biodivine_aeon_server::scc::algo_symbolic_components::{components, components_2}; use biodivine_aeon_server::scc::{Classifier, ProgressTracker}; use biodivine_lib_param_bn::symbolic_async_graph::SymbolicAsyncGraph; use biodivine_lib_param_bn::BooleanNetwork; @@ -32,18 +32,19 @@ fn main() { "Admissible parametrisations: {}", graph.unit_colors().cardinality() ); - println!( - "State space: {}", - graph.unit_vertices().cardinality() - ); + println!("State space: {}", graph.unit_vertices().cardinality()); let classifier = Classifier::new(&graph); let progress = ProgressTracker::new(&graph); - components(&graph, &progress, &AtomicBool::new(false), |component| { - println!("Found attractor..."); - println!("Remaining: {}", progress.get_percent_string()); - classifier.add_component(component, &graph); - }); + components_2( + &graph, + /*&progress, &AtomicBool::new(false), */ + |component| { + println!("Found attractor..."); + println!("Remaining: {}", progress.get_percent_string()); + classifier.add_component(component, &graph); + }, + ); classifier.print(); diff --git a/src/main.rs b/src/main.rs index 6309af6..da20bbb 100644 --- a/src/main.rs +++ b/src/main.rs @@ -21,16 +21,16 @@ use std::convert::TryFrom; mod test_main; use biodivine_lib_param_bn::symbolic_async_graph::{GraphColors, SymbolicAsyncGraph}; +use biodivine_lib_std::collections::bitvectors::{ArrayBitVector, BitVector}; use rocket::config::Environment; use rocket::{Config, Data}; -use std::collections::HashMap; +use std::cmp::max; +use std::collections::{HashMap, HashSet}; use std::io::Read; use std::sync::atomic::{AtomicBool, Ordering}; -use std::sync::{Arc, Mutex, RwLock}; +use std::sync::{Arc, RwLock}; use std::thread::JoinHandle; use std::time::{Duration, SystemTime, UNIX_EPOCH}; -use biodivine_lib_std::collections::bitvectors::{ArrayBitVector, BitVector}; -use std::cmp::max; /// Computation keeps all information struct Computation { @@ -75,8 +75,8 @@ fn max_parameter_cardinality(function: &FnUpdate) -> usize { FnUpdate::Not(inner) => max_parameter_cardinality(inner), FnUpdate::Binary(_, left, right) => max( max_parameter_cardinality(left), - max_parameter_cardinality(right) - ) + max_parameter_cardinality(right), + ), }; } @@ -92,29 +92,38 @@ fn check_update_function(data: Data) -> BackendResponse { let lock = CHECK_UPDATE_FUNCTION_LOCK.clone(); let mut lock = lock.write().unwrap(); let start = SystemTime::now(); - let graph = BooleanNetwork::try_from(model_string.as_str()).and_then(|model| { - let mut max_size = 0; - for v in model.graph().variable_ids() { - if let Some(update_function) = model.get_update_function(v) { - max_size = max(max_size, max_parameter_cardinality(update_function)); + let graph = BooleanNetwork::try_from(model_string.as_str()) + .and_then(|model| { + let mut max_size = 0; + for v in model.graph().variable_ids() { + if let Some(update_function) = model.get_update_function(v) { + max_size = max(max_size, max_parameter_cardinality(update_function)); + } else { + max_size = max(max_size, model.graph().regulators(v).len()) + } + } + if max_size <= 4 { + println!( + "Start partial function analysis. {} variables and complexity {}.", + model.graph().num_vars(), + max_size + ); + SymbolicAsyncGraph::new(model) } else { - max_size = max(max_size, model.graph().regulators(v).len()) + Err("Function too large for on-the-fly analysis.".to_string()) } - } - if max_size <= 4 { - println!("Start partial function analysis. {} variables and complexity {}.", model.graph().num_vars(), max_size); - SymbolicAsyncGraph::new(model) - } else { - Err("Function too large for on-the-fly analysis.".to_string()) - } - }).map(|g| g.unit_colors().cardinality()); - println!("Elapsed: {}, result {:?}", start.elapsed().unwrap().as_millis(), graph); + }) + .map(|g| g.unit_colors().cardinality()); + println!( + "Elapsed: {}, result {:?}", + start.elapsed().unwrap().as_millis(), + graph + ); (*lock) = !(*lock); match graph { - Ok(cardinality) => BackendResponse::ok(&format!( - "{{\"cardinality\":\"{}\"}}", - cardinality - )), + Ok(cardinality) => { + BackendResponse::ok(&format!("{{\"cardinality\":\"{}\"}}", cardinality)) + } Err(error) => BackendResponse::err(&error), } } @@ -324,99 +333,132 @@ fn get_attractors(class_str: String) -> BackendResponse { if let Some(has_class) = try_get_class_params(classifier, &class) { if let Some(class) = has_class { if let Some(graph) = &cmp.graph { - let witness: BooleanNetwork = graph.make_witness(&class); - let witness_str = format!("{}", witness); - let mut all_variables = vec![]; - for id in witness.graph().variable_ids() { - all_variables.push(format!( - "{:?}", - witness.graph().get_variable(id).to_string() + let witness_colour = class.pick_color(graph).unwrap(); + let witness_network: BooleanNetwork = graph.make_witness(&class); + let witness_graph = + SymbolicAsyncGraph::new(witness_network.clone()).unwrap(); + let witness_str = witness_network.to_string(); + let witness_attractors = classifier.attractors(&witness_colour, graph); + let variable_name_strings = + witness_network.graph().variable_ids().map(|id| { + format!( + "\"{}\"", + witness_network.graph().get_variable(id).get_name() + ) + }); + + let mut all_attractors: Vec<( + Behaviour, + Vec<(ArrayBitVector, ArrayBitVector)>, + HashSet, + )> = Vec::new(); + + // Note that the choice of graph/witness_graph is not arbitrary. + // The attractor set is from the original graph, but source_set/target_set + // are based on the witness_graph. This means they have different number + // of BDD variables inside! + for (attractor, behaviour) in witness_attractors.iter() { + println!( + "Attractor {:?} state count: {}", + behaviour, + attractor.cardinality() + ); + if attractor.cardinality() >= 2_000.0 { + return BackendResponse::err(&format!("Attractor has {} states. Visualisation size limit exceeded.", attractor.cardinality())); + } + let mut attractor_graph: Vec<(ArrayBitVector, ArrayBitVector)> = + Vec::new(); + let mut not_fixed_vars: HashSet = HashSet::new(); + if *behaviour == Behaviour::Stability { + // This is a sink - no edges + assert_eq!(attractor.states(&graph).count(), 1); + let sink: ArrayBitVector = + attractor.states(&graph).next().unwrap(); + attractor_graph.push((sink.clone(), sink)); + for i in 0..witness_network.graph().num_vars() { + // In sink, we mark everything as "not-fixed" because we want to just display it normally. + not_fixed_vars.insert(i); + } + } else { + for source in attractor.states(&graph) { + let source_set = witness_graph.vertex(&source); + let mut target_set = witness_graph.empty_vertices().clone(); + for v in witness_graph.network().graph().variable_ids() { + let post = witness_graph.any_post(v, &source_set); + if !post.is_empty() { + not_fixed_vars.insert(v.into()); + target_set = target_set.union(&post); + } + } + + for target in target_set + .state_projection(&witness_graph) + .states(&witness_graph) + { + attractor_graph.push((source.clone(), target)); + } + } + } + + all_attractors.push(( + behaviour.clone(), + attractor_graph, + not_fixed_vars, )); } - let witness_graph = SymbolicAsyncGraph::new(witness); - match witness_graph { - Ok(witness_graph) => { - let all_attractors: Mutex)>> = Mutex::new(vec![]); - - // This is a witness network, so there is exactly ONE parametrisation... - biodivine_aeon_server::scc::algo_symbolic_components::components( - &witness_graph, - &ProgressTracker::new(&witness_graph), - &AtomicBool::new(false), - |component| { - let mut attractor_graph: Vec<(ArrayBitVector, ArrayBitVector)> = vec![]; - - let behaviour = Classifier::classify_component(&component, &witness_graph); - assert_eq!(behaviour.len(), 1); - let behaviour = behaviour.into_iter().next().unwrap().0; - - for source in component.state_projection(&witness_graph).states(&witness_graph) { - let source_set = witness_graph.vertex(&source); - let mut target_set = witness_graph.empty_vertices().clone(); - for v in witness_graph.network().graph().variable_ids() { - let post = witness_graph.any_post(v, &source_set); - target_set = target_set.union(&post); - } - - let mut is_sink = true; - for target in target_set.state_projection(&witness_graph).states(&witness_graph) { - is_sink = false; - attractor_graph.push((source.clone(), target)); - } - - if is_sink { - attractor_graph.push((source.clone(), source.clone())); - } - } - all_attractors - .lock() - .unwrap() - .push((behaviour, attractor_graph)); - }, - ); - // now the data is stored in `all_attractors`, just convert it to json: - let mut json = String::new(); - - for (i, (behavior, graph)) in - all_attractors.into_inner().unwrap().iter().enumerate() - { - if i != 0 { - json += ","; - } // json? no trailing commas for you - json += &format!( - "{{\"class\":\"{:?}\", \"graph\":[", - behavior - ); - let mut edge_count = 0; - for (j, edge) in graph.iter().enumerate() { - let from: String = format!("{:?}", edge.0.values()); - let to: String = format!("{:?}", edge.1.values()); - if j != 0 { - json += "," + // now the data is stored in `all_attractors`, just convert it to json: + let mut json = String::new(); + + for (i, (behavior, graph, not_fixed)) in + all_attractors.iter().enumerate() + { + if i != 0 { + json += ","; + } // json? no trailing commas for you + json += &format!("{{\"class\":\"{:?}\", \"graph\":[", behavior); + let mut edge_count = 0; + for (j, edge) in graph.iter().enumerate() { + fn state_to_binary( + state: &ArrayBitVector, + not_fixed: &HashSet, + ) -> String { + let mut result = String::new(); + for i in 0..state.len() { + if not_fixed.contains(&i) { + result.push(if state.get(i) { '1' } else { '0' }); + } else { + result.push(if state.get(i) { + '⊤' + } else { + '⊥' + }); } - json += &format!("[\"{}\", \"{}\"]", from, to); - edge_count += 1; } - json += &format!("], \"edges\":{}}}", edge_count); + return result; } - json = "{ \"attractors\":[".to_owned() - + &json - + "], \"variables\":["; - for (i, var) in all_variables.iter().enumerate() { - if i != 0 { - json += ","; - } - json += var; + let from: String = state_to_binary(&edge.0, not_fixed); + let to: String = state_to_binary(&edge.1, not_fixed); + if j != 0 { + json += "," } - json += &format!( - "], \"model\":{}", - &object! { "model" => witness_str }.to_string() - ); - BackendResponse::ok(&(json + "}")) + json += &format!("[\"{}\", \"{}\"]", from, to); + edge_count += 1; + } + json += &format!("], \"edges\":{}}}", edge_count); + } + json = "{ \"attractors\":[".to_owned() + &json + "], \"variables\":["; + for (i, var) in variable_name_strings.enumerate() { + if i != 0 { + json += ","; } - Err(error_msg) => BackendResponse::err(&error_msg), + json += var.as_str(); } + json += &format!( + "], \"model\":{}", + &object! { "model" => witness_str }.to_string() + ); + BackendResponse::ok(&(json + "}")) } else { return BackendResponse::err(&"No results available.".to_string()); } @@ -747,7 +789,8 @@ impl BackendResponse { message: object! { "status" => false, "message" => message.replace("\n", "
").clone(), - }.to_string(), + } + .to_string(), }; } } diff --git a/src/scc/algo_effectively_constant.rs b/src/scc/algo_effectively_constant.rs index a547e4f..7971525 100644 --- a/src/scc/algo_effectively_constant.rs +++ b/src/scc/algo_effectively_constant.rs @@ -1,7 +1,15 @@ -use biodivine_lib_param_bn::symbolic_async_graph::{SymbolicAsyncGraph, GraphColoredVertices}; +use crate::scc::algo_reach::reach; +use biodivine_lib_param_bn::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph}; use biodivine_lib_param_bn::VariableId; +use std::cmp::min; +use std::collections::HashSet; +use std::env::var; +use std::fs::read; -pub fn remove_effectively_constant_states_2(graph: &SymbolicAsyncGraph, set: GraphColoredVertices) -> GraphColoredVertices { +pub fn remove_effectively_constant_states_2( + graph: &SymbolicAsyncGraph, + set: GraphColoredVertices, +) -> GraphColoredVertices { println!("Remove effectively constant states."); let mut universe = set; for variable in graph.network().graph().variable_ids() { @@ -10,7 +18,11 @@ pub fn remove_effectively_constant_states_2(graph: &SymbolicAsyncGraph, set: Gra return universe; } -fn cut_variable(graph: &SymbolicAsyncGraph, variable: VariableId, universe: GraphColoredVertices) -> GraphColoredVertices { +fn cut_variable( + graph: &SymbolicAsyncGraph, + variable: VariableId, + universe: GraphColoredVertices, +) -> GraphColoredVertices { let mut scc_candidate = universe.clone(); let mut not_attractor = graph.empty_vertices().clone(); @@ -19,7 +31,8 @@ fn cut_variable(graph: &SymbolicAsyncGraph, variable: VariableId, universe: Grap loop { let vertices_where_var_can_jump = graph.has_any_post(variable, &scc_candidate); let vertices_where_var_jumped = graph.any_post(variable, &scc_candidate); - let reachable_after_jump = reach_fwd_excluding(graph, &vertices_where_var_jumped, &scc_candidate, variable); + let reachable_after_jump = + reach_fwd_excluding(graph, &vertices_where_var_jumped, &scc_candidate, variable); let can_jump_again = reachable_after_jump.intersect(&vertices_where_var_can_jump); let will_never_jump_again = vertices_where_var_can_jump.minus(&can_jump_again); scc_candidate = reachable_after_jump; @@ -27,12 +40,21 @@ fn cut_variable(graph: &SymbolicAsyncGraph, variable: VariableId, universe: Grap break; } not_attractor = not_attractor.union(&will_never_jump_again); - println!("{:?} will never jump again: {}", variable, will_never_jump_again.cardinality()); + println!( + "{:?} will never jump again: {}", + variable, + will_never_jump_again.cardinality() + ); } // Now finish the SCCs let vertices_where_var_can_jump = graph.has_any_post(variable, &scc_candidate); - let scc = reach_bwd_excluding(graph, &vertices_where_var_can_jump, &scc_candidate, variable); + let scc = reach_bwd_excluding( + graph, + &vertices_where_var_can_jump, + &scc_candidate, + variable, + ); let not_scc = graph.unit_vertices().minus(&scc); for other_variable in graph.network().graph().variable_ids() { let can_jump_out = graph.any_pre(other_variable, ¬_scc).intersect(&scc); @@ -43,7 +65,11 @@ fn cut_variable(graph: &SymbolicAsyncGraph, variable: VariableId, universe: Grap } let to_remove = reach_bwd_excluding(graph, ¬_attractor, &universe, variable); - println!("Eliminated: {}/{}", to_remove.cardinality(), universe.cardinality()); + println!( + "Eliminated: {}/{}", + to_remove.cardinality(), + universe.cardinality() + ); return universe.minus(&to_remove); } @@ -55,39 +81,183 @@ fn cut_variable(graph: &SymbolicAsyncGraph, variable: VariableId, universe: Grap /// /// Note that this does not mean the variable has to strictly always jump - that is why we need the /// backward reachability to detect states that can actually achieve irreversible jump. -pub fn remove_effectively_constant_states(graph: &SymbolicAsyncGraph, set: GraphColoredVertices) -> GraphColoredVertices { +pub fn remove_effectively_constant_states( + graph: &SymbolicAsyncGraph, + set: GraphColoredVertices, +) -> (GraphColoredVertices, Vec) { println!("Remove effectively constant states."); let original_size = set.cardinality(); let mut universe = set; let mut stop = false; + let mut variables: Vec = graph.network().graph().variable_ids().collect(); while !stop { stop = true; let mut to_remove = graph.empty_vertices().clone(); - for variable in graph.network().graph().variable_ids() { - let vertices_where_var_can_jump = graph.has_post(variable, &universe); - let vertices_where_var_jumped = graph.post(variable, &universe, &universe); - let reachable_after_jump = reach_fwd_excluding(graph, &vertices_where_var_jumped, &universe, variable); + for variable in &variables { + let active_variables: Vec = variables + .iter() + .filter(|it| *it != variable) + .cloned() + .collect(); + let vertices_where_var_can_jump = graph.has_post(*variable, &universe); + let vertices_where_var_jumped = graph.post(*variable, &universe, &universe); + let reachable_after_jump = reach_saturated_fwd_excluding( + graph, + &vertices_where_var_jumped, + &universe, + &active_variables, + ); let can_jump_again = reachable_after_jump.intersect(&vertices_where_var_can_jump); let will_never_jump_again = vertices_where_var_can_jump.minus(&can_jump_again); - if !will_never_jump_again.is_empty() { - println!("{:?} will never jump again: {}", variable, will_never_jump_again.cardinality()); - } if !will_never_jump_again.is_empty() { stop = false; - let to_remove_for_var = reach_bwd_excluding(graph, &will_never_jump_again, &universe, variable); + let to_remove_for_var = reach_saturated_bwd_excluding( + graph, + &will_never_jump_again, + &universe, + &active_variables, + ); to_remove = to_remove.union(&to_remove_for_var); - //universe = universe.minus(&to_remove); - println!("Eliminated {}/{} {:+e}%", to_remove.cardinality(), universe.cardinality(), (to_remove.cardinality()/universe.cardinality()) * 100.0); + //universe = universe.minus(&to_remove); THIS IS A BAD IDEA... + println!( + "{:?} will never jump again: {}", + variable, + will_never_jump_again.cardinality() + ); + println!( + "Eliminated {}/{} {:+e}%", + to_remove.cardinality(), + universe.cardinality(), + (to_remove.cardinality() / universe.cardinality()) * 100.0 + ); } } universe = universe.minus(&to_remove); - println!("Universe now has {} nodes.", universe.clone().into_bdd().size()); + let original_vars = variables.len(); + variables = variables + .into_iter() + .filter(|v| !graph.has_post(*v, &universe).is_empty()) + .collect(); + println!( + "Universe now has {} nodes. Eliminated {} variables.", + universe.clone().into_bdd().size(), + original_vars - variables.len() + ); } - println!("Removed {}/{} {:+e}%; {} nodes.", universe.cardinality(), original_size, (universe.cardinality()/original_size) * 100.0, universe.clone().into_bdd().size()); - return universe; + + println!("Final active variables: {}", variables.len()); + println!( + "Removed {}/{} {:+e}%; {} nodes.", + universe.cardinality(), + original_size, + (universe.cardinality() / original_size) * 100.0, + universe.clone().into_bdd().size() + ); + + for v in &variables { + let vertices_where_var_can_jump = graph.has_post(*v, &universe); + let reachable_before_jump = reach_saturated_bwd_excluding( + graph, + &vertices_where_var_can_jump, + &universe, + &variables, + ); + let reachable_after_jump = reach_saturated_fwd_excluding( + graph, + &vertices_where_var_can_jump, + &universe, + &variables, + ); + let components = reachable_before_jump.intersect(&reachable_after_jump); + let below = reachable_after_jump.minus(&components); + let can_reach_below = + reach_saturated_bwd_excluding(graph, &below, &universe, &variables).minus(&below); + println!( + "({:?}) Below: {} Can reach below: {}", + v, + below.cardinality(), + can_reach_below.cardinality() + ); + universe = universe.minus(&can_reach_below); + } + + println!("Final active variables: {}", variables.len()); + println!( + "Removed {}/{} {:+e}%; {} nodes.", + universe.cardinality(), + original_size, + (universe.cardinality() / original_size) * 100.0, + universe.clone().into_bdd().size() + ); + return (universe, variables); +} + +pub fn reach_saturated_fwd_excluding( + graph: &SymbolicAsyncGraph, + initial: &GraphColoredVertices, + guard: &GraphColoredVertices, + variables: &Vec, +) -> GraphColoredVertices { + if variables.is_empty() { + return initial.clone(); + } + let mut result = initial.clone(); + let last_variable = variables.len() - 1; + let mut active_variable = last_variable; + loop { + let variable = variables[active_variable]; + let post = graph.post(variable, &result, guard).minus(&result); + result = result.union(&post); + + if !post.is_empty() { + active_variable = last_variable; + } else { + if active_variable == 0 { + break; + } else { + active_variable -= 1; + } + } + } + return result; +} + +pub fn reach_saturated_bwd_excluding( + graph: &SymbolicAsyncGraph, + initial: &GraphColoredVertices, + guard: &GraphColoredVertices, + variables: &Vec, +) -> GraphColoredVertices { + if variables.is_empty() { + return initial.clone(); + } + let mut result = initial.clone(); + let last_variable = variables.len() - 1; + let mut active_variable = last_variable; + loop { + let variable = variables[active_variable]; + let post = graph.pre(variable, &result, guard).minus(&result); + result = result.union(&post); + + if !post.is_empty() { + active_variable = last_variable; + } else { + if active_variable == 0 { + break; + } else { + active_variable -= 1; + } + } + } + return result; } -pub fn reach_fwd_excluding(graph: &SymbolicAsyncGraph, initial: &GraphColoredVertices, guard: &GraphColoredVertices, exclude: VariableId) -> GraphColoredVertices { +pub fn reach_fwd_excluding( + graph: &SymbolicAsyncGraph, + initial: &GraphColoredVertices, + guard: &GraphColoredVertices, + exclude: VariableId, +) -> GraphColoredVertices { let mut result = initial.clone(); //println!("Reach fwd excluding {:?}...", exclude); loop { @@ -123,7 +293,12 @@ pub fn reach_fwd_excluding(graph: &SymbolicAsyncGraph, initial: &GraphColoredVer return result; } -pub fn reach_bwd_excluding(graph: &SymbolicAsyncGraph, initial: &GraphColoredVertices, guard: &GraphColoredVertices, exclude: VariableId) -> GraphColoredVertices { +pub fn reach_bwd_excluding( + graph: &SymbolicAsyncGraph, + initial: &GraphColoredVertices, + guard: &GraphColoredVertices, + exclude: VariableId, +) -> GraphColoredVertices { let mut result = initial.clone(); //println!("Reach bwd excluding {:?}...", exclude); loop { @@ -156,4 +331,4 @@ pub fn reach_bwd_excluding(graph: &SymbolicAsyncGraph, initial: &GraphColoredVer } return result; -} \ No newline at end of file +} diff --git a/src/scc/algo_par_utils.rs b/src/scc/algo_par_utils.rs index 65d4863..e6fc819 100644 --- a/src/scc/algo_par_utils.rs +++ b/src/scc/algo_par_utils.rs @@ -1,14 +1,18 @@ use rayon::prelude::*; // TODO: Something, something, FUTURES! -pub fn par_fold(items: Vec, action: F) -> T where F: Fn(&T, &T) -> T + Sync { +pub fn par_fold(items: Vec, action: F) -> T +where + F: Fn(&T, &T) -> T + Sync, +{ if items.is_empty() { panic!("Empty parallel fold"); } else if items.len() == 1 { return items[0].clone(); } else { let data: &[T] = &items; - let joined: Vec = data.par_iter() + let joined: Vec = data + .par_iter() .chunks(2) .map(|chunk| { if chunk.len() == 2 { @@ -16,7 +20,8 @@ pub fn par_fold(items: Vec, action: F) -> T where } else { chunk[0].clone() } - }).collect(); + }) + .collect(); return par_fold(joined, action); } -} \ No newline at end of file +} diff --git a/src/scc/algo_symbolic_components.rs b/src/scc/algo_symbolic_components.rs index c0608ae..53c2f26 100644 --- a/src/scc/algo_symbolic_components.rs +++ b/src/scc/algo_symbolic_components.rs @@ -1,14 +1,21 @@ +use crate::scc::algo_effectively_constant::{ + reach_saturated_bwd_excluding, reach_saturated_fwd_excluding, + remove_effectively_constant_states, +}; use crate::scc::algo_symbolic_reach::{guarded_reach_bwd, guarded_reach_fwd}; use crate::scc::ProgressTracker; use biodivine_lib_param_bn::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph}; +use biodivine_lib_param_bn::VariableId; use biodivine_lib_std::param_graph::Params; -use std::option::Option::Some; -use std::sync::atomic::{AtomicBool, Ordering}; use std::io; use std::io::Write; -use crate::scc::algo_effectively_constant::{remove_effectively_constant_states}; +use std::option::Option::Some; +use std::sync::atomic::{AtomicBool, Ordering}; -pub fn prune_sources(graph: &SymbolicAsyncGraph, set: GraphColoredVertices) -> GraphColoredVertices { +pub fn prune_sources( + graph: &SymbolicAsyncGraph, + set: GraphColoredVertices, +) -> GraphColoredVertices { let start = set.cardinality(); let mut result = set; println!("Pruning sources..."); @@ -20,28 +27,32 @@ pub fn prune_sources(graph: &SymbolicAsyncGraph, set: GraphColoredVertices) -> G let jumps_false_to_true = graph.has_any_post(v, &result.intersect(&v_false)); let jumps_true_to_false = graph.has_any_post(v, &result.intersect(&v_true)); /* - Now we are looking for colors that only jump one way: for these, we can fix the value - because they will never flip back... - */ - let colors_false_to_true = jumps_false_to_true.color_projection(); - let colors_true_to_false = jumps_true_to_false.color_projection(); + Now we are looking for colors that only jump one way: for these, we can fix the value + because they will never flip back... + */ + let colors_false_to_true = jumps_false_to_true.color_projection(graph); + let colors_true_to_false = jumps_true_to_false.color_projection(graph); let colors_that_jump_both_ways = colors_true_to_false.intersect(&colors_true_to_false); - let colors_that_jump_only_one_way = (colors_false_to_true.union(&colors_true_to_false)).minus(&colors_that_jump_both_ways); + let colors_that_jump_only_one_way = (colors_false_to_true.union(&colors_true_to_false)) + .minus(&colors_that_jump_both_ways); if !colors_that_jump_only_one_way.is_empty() { found_something = true; - let one_way_colors_false_to_true = colors_false_to_true.intersect(&colors_that_jump_only_one_way); - let one_way_colors_true_to_false = colors_true_to_false.intersect(&colors_that_jump_only_one_way); + let one_way_colors_false_to_true = + colors_false_to_true.intersect(&colors_that_jump_only_one_way); + let one_way_colors_true_to_false = + colors_true_to_false.intersect(&colors_that_jump_only_one_way); result = result.minus(&v_true.intersect_colors(&one_way_colors_true_to_false)); result = result.minus(&v_false.intersect_colors(&one_way_colors_false_to_true)); } println!("{:?} -> {}", v, colors_that_jump_only_one_way.cardinality()); } - println!("Keeping {}/{} ({:+e}%, nodes result({}))", - result.cardinality(), - start, - (result.cardinality()/start) * 100.0, - result.clone().into_bdd().size() + println!( + "Keeping {}/{} ({:+e}%, nodes result({}))", + result.cardinality(), + start, + (result.cardinality() / start) * 100.0, + result.clone().into_bdd().size() ); if !found_something { return result; @@ -49,14 +60,159 @@ pub fn prune_sources(graph: &SymbolicAsyncGraph, set: GraphColoredVertices) -> G } } +pub fn components_2(graph: &SymbolicAsyncGraph, on_component: F) +where + F: Fn(GraphColoredVertices) -> () + Send + Sync, +{ + let mut universe = graph.unit_vertices().clone(); + + /*let mut can_be_sink = graph.unit_vertices().clone(); // intentionally use all vertices + //panic!(""); + for variable in graph.network().graph().variable_ids() { + print!("{:?}...", variable); + io::stdout().flush().unwrap(); + let has_successor = &graph.has_any_post(variable, graph.unit_vertices()); + can_be_sink = can_be_sink.minus(has_successor); + } + println!(); + + let mut is_sink = can_be_sink.clone(); + + let sinks = is_sink.clone(); + // Now we have to report sinks, but we have to satisfy that every reported set has only one component: + while !is_sink.is_empty() { + let to_report = is_sink.pivots(graph); + is_sink = is_sink.minus(&to_report); + on_component(to_report); + } + + println!("Sinks detected: {}", sinks.cardinality()); + + let variables: Vec = graph.network().graph().variable_ids().collect(); + universe = universe.minus(&reach_saturated_bwd_excluding(graph, &sinks, graph.unit_vertices(), &variables)); + if universe.is_empty() { + println!("Everything is in sink!"); + return; + }*/ + + let (constrained, variables) = remove_effectively_constant_states(graph, universe); + universe = constrained.clone(); + + let mut attractor_candidates = graph.empty_vertices().clone(); + while !universe.is_empty() { + println!("Pick a new pivot branch..."); + let mut pivot = universe.pivots(graph); + //while !pivot.is_empty() { + let bwd_pivot = reach_saturated_bwd_excluding(graph, &pivot, &universe, &variables); + let component_with_pivot = + reach_saturated_fwd_excluding(graph, &pivot, &bwd_pivot, &variables); + let after_component = post(graph, &component_with_pivot).minus(&component_with_pivot); + let is_candidate = component_with_pivot + .color_projection(graph) + .minus(&after_component.color_projection(graph)); + if !is_candidate.is_empty() { + //attractor_candidates = attractor_candidates.union(&component_with_pivot.intersect_colors(&is_candidate)); + on_component(component_with_pivot.intersect_colors(&is_candidate)); + } + pivot = after_component; + universe = universe.minus(&bwd_pivot); + println!( + "Attractor candidates: {} Universe: {}", + attractor_candidates.cardinality(), + universe.cardinality() + ); + //} + } + println!( + "Attractor candidates: {}/{}", + attractor_candidates.cardinality(), + constrained.cardinality() + ); +} + +/* +/// Compute the SCCs that are induced by the vertices in `initial`. +/// +/// Uses lockstep, so in terms of symbolic steps, this should be optimal. +fn components_with(graph: &SymbolicAsyncGraph, initial: &GraphColoredVertices) -> GraphColoredVertices { + println!("Lockstep component search."); + let mut still_running = initial.color_projection(graph); + let mut fwd = initial.clone(); + let mut bwd = initial.clone(); + let mut fwd_frontier = initial.clone(); + let mut bwd_frontier = initial.clone(); + let mut fwd_unfinished = graph.empty_vertices().clone(); + let mut bwd_unfinished = graph.empty_vertices().clone(); + // First, find fwd/bwd templates using lockstep. + while !still_running.is_empty() { + println!("(1) Still running {}; {} {} {} {}", still_running.cardinality(), fwd.as_bdd().size(), bwd.as_bdd().size(), fwd_unfinished.as_bdd().size(), bwd_unfinished.as_bdd().size()); + let post = post(graph, &fwd_frontier).minus(&fwd); + let pre = pre(graph, &bwd_frontier).minus(&bwd); + // Compute colours that FINISHED in this iteration. + let fwd_finished_now = fwd_frontier.color_projection(graph).minus(&post.color_projection(graph)); + let bwd_finished_now = bwd_frontier.color_projection(graph).minus(&pre.color_projection(graph)); + let bwd_finished_now = bwd_finished_now.minus(&fwd_finished_now); // Resolve ties. + // Remove BWD-FINISHED colours from FWD frontier and move all vertices with these colours to unfinished. + fwd_frontier = post.minus_colors(&bwd_finished_now); + fwd_unfinished = fwd_unfinished.union(&fwd.intersect_colors(&bwd_finished_now)); + fwd = fwd.minus_colors(&bwd_finished_now).union(&fwd_frontier); + // Remove FWD-FINISHED colours from BWD frontier and move all vertices with these colours to unfinished. + bwd_frontier = pre.minus_colors(&fwd_finished_now); + bwd_unfinished = bwd_unfinished.union(&bwd.intersect_colors(&fwd_finished_now)); + bwd = bwd.minus_colors(&fwd_finished_now).union(&bwd_frontier); + // Mark finished colours as done. + still_running = still_running.minus(&fwd_finished_now.union(&bwd_finished_now)); + } + // Now fwd and bwd contain VALID forward/backward sets and are colour disjoint. + // We have to continue computing unfinished within these valid bounds. + assert!(fwd.color_projection(graph).intersect(&bwd.color_projection(graph)).is_empty()); + assert!(initial.color_projection(graph).minus(&fwd.color_projection(graph)).minus(&bwd.color_projection(graph)).is_empty()); + + loop { + let post = post(graph, &fwd_unfinished).minus(&fwd_unfinished).intersect(&bwd); + let pre = pre(graph, &bwd_unfinished).minus(&bwd_unfinished).intersect(&fwd); + fwd_unfinished = fwd_unfinished.union(&post); + bwd_unfinished = bwd_unfinished.union(&pre); + println!("(2) Still running {}", post.cardinality() + pre.cardinality()); + if post.is_empty() && pre.is_empty() { + break; + } + } + // At this point, fwd_unfinished and bwd_unfinished are completely done within the bounds of bwd/fwd. + // We can therefore join them back + fwd = fwd.union(&fwd_unfinished); + bwd = bwd.union(&bwd_unfinished); + + // The final result is the intersection of the two. + return fwd.intersect(&bwd); +}*/ + +fn post(graph: &SymbolicAsyncGraph, initial: &GraphColoredVertices) -> GraphColoredVertices { + let mut result = graph.empty_vertices().clone(); + for v in graph.network().graph().variable_ids() { + result = result.union(&graph.any_post(v, initial)); + } + return result; +} + +fn pre(graph: &SymbolicAsyncGraph, initial: &GraphColoredVertices) -> GraphColoredVertices { + let mut result = graph.empty_vertices().clone(); + for v in graph.network().graph().variable_ids() { + result = result.union(&graph.any_pre(v, initial)); + } + return result; +} + pub fn components( graph: &SymbolicAsyncGraph, - progress: &ProgressTracker, - cancelled: &AtomicBool, + _progress: &ProgressTracker, + _cancelled: &AtomicBool, on_component: F, ) where F: Fn(GraphColoredVertices) -> () + Send + Sync, { + components_2(graph, on_component); + /* crossbeam::thread::scope(|scope| { println!("Detect eventually stable..."); // TODO: This is not correct, because for parametrisations can_move will never be empty... @@ -116,7 +272,7 @@ pub fn components( let sinks = is_sink.clone(); // Now we have to report sinks, but we have to satisfy that every reported set has only one component: while !is_sink.is_empty() { - let to_report = is_sink.pivots(); + let to_report = is_sink.pivots(graph); is_sink = is_sink.minus(&to_report); on_component(to_report); } @@ -132,7 +288,7 @@ pub fn components( .collect(); let has_successors = par_fold(has_successors, |a, b| a.union(b));*/ - let not_constant = remove_effectively_constant_states(graph, graph.unit_vertices().clone()); + let (not_constant, _) = remove_effectively_constant_states(graph, graph.unit_vertices().clone()); println!("Not constant: {}/{}", not_constant.cardinality(), graph.unit_vertices().cardinality()); if cancelled.load(Ordering::SeqCst) { @@ -170,14 +326,16 @@ pub fn components( let remaining: f64 = queue.iter().map(|(a, _)| *a).sum::() + universe_cardinality; progress.update_remaining(remaining); println!("Look for pivots..."); - let pivots = universe.pivots(); - println!("Pivots cardinality: {}", pivots.cardinality()); + let pivots = universe.pivots(graph); let backward = guarded_reach_bwd(&graph, &pivots, &universe, cancelled, progress); + let pivots = improve_pivots(graph, pivots.clone(), &(universe.minus(&backward).union(&pivots))); + println!("Pivots cardinality: {}", pivots.cardinality()); + let backward = guarded_reach_bwd(&graph, &backward.union(&pivots), &universe, cancelled, progress); let component_with_pivots = guarded_reach_fwd(&graph, &pivots, &backward, cancelled, progress); - let mut is_terminal = component_with_pivots.color_projection(); + let mut is_terminal = component_with_pivots.color_projection(graph); for v in graph.network().graph().variable_ids() { - let successors = graph.any_post(v, &component_with_pivots).minus(&component_with_pivots).color_projection(); + let successors = graph.any_post(v, &component_with_pivots).minus(&component_with_pivots).color_projection(graph); if !successors.is_empty() { is_terminal = is_terminal.minus(&successors); } @@ -240,5 +398,35 @@ pub fn components( println!("Main component loop done..."); }) - .unwrap(); + .unwrap();*/ } + +/*fn improve_pivots(graph: &SymbolicAsyncGraph, pivots: GraphColoredVertices, universe: &GraphColoredVertices) -> GraphColoredVertices { + let mut discovered = pivots.clone(); + let mut final_pivots = graph.empty_vertices().clone(); + let mut frontier = pivots.clone(); + println!("Pivot optimisation..."); + while !frontier.is_empty() { + if discovered.as_bdd().size() > 10_000 { + println!("{}/{} ({:+e}%, nodes result({}))", + discovered.cardinality(), + universe.cardinality(), + (discovered.cardinality() / universe.cardinality()) * 100.0, + discovered.as_bdd().size() + ); + } + let mut new_successor = graph.empty_vertices().clone(); + for v in graph.network().graph().variable_ids() { + new_successor = new_successor.union(&graph.any_post(v, &frontier)); + } + new_successor = new_successor.minus(&discovered); + let stopped_in_this_step = frontier.color_projection(graph).minus(&new_successor.color_projection(graph)); + if !stopped_in_this_step.is_empty() { + // For colours that stopped in this iteration + final_pivots = final_pivots.union(&frontier.intersect_colors(&stopped_in_this_step).pivots(graph)) + } + frontier = new_successor; + discovered = discovered.union(&frontier); + } + return final_pivots; +}*/ diff --git a/src/scc/algo_symbolic_reach.rs b/src/scc/algo_symbolic_reach.rs index b226cab..7519460 100644 --- a/src/scc/algo_symbolic_reach.rs +++ b/src/scc/algo_symbolic_reach.rs @@ -1,8 +1,8 @@ use crate::scc::ProgressTracker; use biodivine_lib_param_bn::symbolic_async_graph::{GraphColoredVertices, SymbolicAsyncGraph}; -use std::sync::atomic::{AtomicBool, Ordering}; use std::io; use std::io::Write; +use std::sync::atomic::{AtomicBool, Ordering}; pub fn guarded_reach_fwd( graph: &SymbolicAsyncGraph, @@ -57,12 +57,15 @@ pub fn guarded_reach_fwd( progress.update_last_wave(result.cardinality()); - println!("{}/{} ({:+e}%, nodes result({}))", - result.cardinality(), - guard.cardinality(), - (result.cardinality()/guard.cardinality()) * 100.0, - result.clone().into_bdd().size() - ); + if result.as_bdd().size() > 10_000 { + println!( + "{}/{} ({:+e}%, nodes result({}))", + result.cardinality(), + guard.cardinality(), + (result.cardinality() / guard.cardinality()) * 100.0, + result.clone().into_bdd().size() + ); + } let mut successors = graph.empty_vertices().clone(); for variable in graph.network().graph().variable_ids() { io::stdout().flush().unwrap(); @@ -75,8 +78,12 @@ pub fn guarded_reach_fwd( println!("{:?} -> {}", variable, s.into_bdd().size()); }*/ } - print!(" || {}", successors.clone().into_bdd().size()); - println!(); + if result.as_bdd().size() > 10_000 { + println!( + "Successor node count: {}", + successors.clone().into_bdd().size() + ); + } successors = successors.minus(&result); if successors.is_empty() { break; @@ -156,12 +163,15 @@ pub fn guarded_reach_bwd( progress.update_last_wave(result.cardinality()); - println!("{}/{} ({:+e}%, nodes result({}))", - result.cardinality(), - guard.cardinality(), - (result.cardinality()/guard.cardinality()) * 100.0, - result.clone().into_bdd().size() - ); + if result.as_bdd().size() > 10_000 { + println!( + "{}/{} ({:+e}%, nodes result({}))", + result.cardinality(), + guard.cardinality(), + (result.cardinality() / guard.cardinality()) * 100.0, + result.clone().into_bdd().size() + ); + } let mut predecessors = graph.empty_vertices().clone(); for variable in graph.network().graph().variable_ids() { io::stdout().flush().unwrap(); @@ -174,8 +184,12 @@ pub fn guarded_reach_bwd( println!("{:?} -> {}", variable, s.into_bdd().size()); }*/ } - print!(" || {}", predecessors.clone().into_bdd().size()); - println!(); + if result.as_bdd().size() > 10_000 { + print!( + "Predecessor node count {}", + predecessors.clone().into_bdd().size() + ); + } predecessors = predecessors.minus(&result); if predecessors.is_empty() { break; diff --git a/src/scc/impl_classifier.rs b/src/scc/impl_classifier.rs index bb7cd33..e9c5f2a 100644 --- a/src/scc/impl_classifier.rs +++ b/src/scc/impl_classifier.rs @@ -1,6 +1,6 @@ use super::{Behaviour, Class, Classifier}; use biodivine_lib_param_bn::symbolic_async_graph::{ - GraphColoredVertices, GraphColors, SymbolicAsyncGraph, + GraphColoredVertices, GraphColors, GraphVertices, SymbolicAsyncGraph, }; use biodivine_lib_std::param_graph::Params; use std::collections::HashMap; @@ -12,6 +12,7 @@ impl Classifier { map.insert(Class::new_empty(), graph.unit_colors().clone()); return Classifier { classes: Mutex::new(map), + attractors: Mutex::new(Vec::new()), }; } @@ -50,13 +51,16 @@ impl Classifier { } /// Static function to classify just one component and immediately obtain results. - pub fn classify_component(component: &GraphColoredVertices, graph: &SymbolicAsyncGraph) -> HashMap { + pub fn classify_component( + component: &GraphColoredVertices, + graph: &SymbolicAsyncGraph, + ) -> HashMap { let classifier = Classifier::new(graph); classifier.add_component(component.clone(), graph); let mut result: HashMap = HashMap::new(); for (class, colors) in classifier.export_result() { if class.0.is_empty() { - continue // This is an empty class - those colors were not in the attractor. + continue; // This is an empty class - those colors were not in the attractor. } else if class.0.len() > 1 { unreachable!("Multiple behaviours in one component."); } else { @@ -66,33 +70,71 @@ impl Classifier { return result; } + /// Find attractor of the given witness colour. The argument set must be a singleton. + pub fn attractors( + &self, + witness_colour: &GraphColors, + graph: &SymbolicAsyncGraph, + ) -> Vec<(GraphVertices, Behaviour)> { + if witness_colour.cardinality() != 1.0 { + eprintln!("WARNING: Computing attractor witnesses for non-singleton set. (This may be just a floating point error in large models)."); + } + let mut result = Vec::new(); + let attractors = self.attractors.lock().unwrap(); + for (attractor, behaviour) in attractors.iter() { + let attractor_states = attractor.intersect_colors(witness_colour); + if attractor_states.is_empty() { + continue; + } + let attractor_states = attractor_states.state_projection(graph); + let attractor_behaviour = behaviour + .iter() + .find(|(_, c)| witness_colour.is_subset(c)) + .unwrap() + .0 + .clone(); + result.push((attractor_states, attractor_behaviour)); + } + return result; + } + // TODO: Parallelism pub fn add_component(&self, component: GraphColoredVertices, graph: &SymbolicAsyncGraph) { - let without_sinks = self.filter_sinks(component, graph); - let not_sink_params = without_sinks.color_projection(); - if not_sink_params.is_empty() { - return; + let mut component_classification = HashMap::new(); + let without_sinks = self.filter_sinks(component.clone(), graph); + let not_sink_params = without_sinks.color_projection(graph); + let sink_params = component.color_projection(graph).minus(¬_sink_params); + if !sink_params.is_empty() { + component_classification.insert(Behaviour::Stability, sink_params); } - let mut disorder = graph.empty_colors().clone(); - for variable in graph.network().graph().variable_ids() { - let found_first_successor = &graph.has_any_post(variable, &without_sinks); - for next_variable in graph.network().graph().variable_ids() { - if next_variable == variable { - continue; + if !not_sink_params.is_empty() { + let mut disorder = graph.empty_colors().clone(); + for variable in graph.network().graph().variable_ids() { + let found_first_successor = &graph.has_any_post(variable, &without_sinks); + for next_variable in graph.network().graph().variable_ids() { + if next_variable == variable { + continue; + } + let found_second_successor = + &graph.has_any_post(next_variable, &found_first_successor); + disorder = disorder.union(&found_second_successor.color_projection(graph)); } - let found_second_successor = - &graph.has_any_post(next_variable, &found_first_successor); - disorder = disorder.union(&found_second_successor.color_projection()); + } + let cycle = without_sinks.color_projection(graph).minus(&disorder); + if !cycle.is_empty() { + println!("Found cycle: {}", cycle.cardinality()); + component_classification.insert(Behaviour::Oscillation, cycle.clone()); + self.push(Behaviour::Oscillation, cycle); + } + if !disorder.is_empty() { + println!("Found disorder: {}", disorder.cardinality()); + component_classification.insert(Behaviour::Disorder, disorder.clone()); + self.push(Behaviour::Disorder, disorder); } } - let cycle = without_sinks.color_projection().minus(&disorder); - if !cycle.is_empty() { - println!("Found cycle: {}", cycle.cardinality()); - self.push(Behaviour::Oscillation, cycle); - } - if !disorder.is_empty() { - println!("Found disorder: {}", disorder.cardinality()); - self.push(Behaviour::Disorder, disorder); + { + let mut attractors = self.attractors.lock().unwrap(); + attractors.push((component, component_classification)); } } @@ -149,8 +191,8 @@ impl Classifier { } } let is_sink = component - .color_projection() - .minus(&is_not_sink.color_projection()); + .color_projection(graph) + .minus(&is_not_sink.color_projection(graph)); if !is_sink.is_empty() { self.push(Behaviour::Stability, is_sink); } diff --git a/src/scc/mod.rs b/src/scc/mod.rs index 646813c..05afb59 100644 --- a/src/scc/mod.rs +++ b/src/scc/mod.rs @@ -1,17 +1,17 @@ use biodivine_lib_param_bn::bdd_params::BddParams; -use biodivine_lib_param_bn::symbolic_async_graph::GraphColors; +use biodivine_lib_param_bn::symbolic_async_graph::{GraphColoredVertices, GraphColors}; use std::cmp::Ordering; use std::collections::HashMap; use std::sync::Mutex; use std::vec::IntoIter; pub mod algo_components; +pub mod algo_effectively_constant; pub mod algo_par_reach; +pub mod algo_par_utils; pub mod algo_reach; pub mod algo_symbolic_components; pub mod algo_symbolic_reach; -pub mod algo_par_utils; -pub mod algo_effectively_constant; mod impl_class; mod impl_classifier; mod impl_old_classifier; @@ -61,6 +61,7 @@ impl PartialOrd for Class { pub struct Classifier { //graph: &'a AsyncGraph, classes: Mutex>, + attractors: Mutex)>>, } pub struct OldClassifier {