diff --git a/Cargo.toml b/Cargo.toml index fd1b5c0..7e12c97 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -21,3 +21,8 @@ itertools = "0.14.0" #profile mode #debug = true + +[[bin]] +name = "schaeppert" +path = "src/iterative.rs" + diff --git a/README.md b/README.md index abe13b6..b47c809 100644 --- a/README.md +++ b/README.md @@ -134,4 +134,30 @@ shepherd -f examples/example1.tikz -o report.tex && pdflatex report.tex The states of the NFA can be automatically reordered in order to make the generated reports more readable. Either topologically (`-s topological`) or alphabetically (`-s alphabetical`). +## Comparison with iterating prism + +The binary `schaeppert` implements an alternative approach for the random population control problem by iteratively constructing an explicit n-fold product MDP, +and stopping once the combined target is not almost-surely reachable. +This is of course only a semi-decision procedure, as it may not terminate on controllable instances. + +MDP solving is done using [prism model checker](https://github.com/prismmodelchecker/prism) and we **assume that the `prism` binary is available in your PATH**. +The product MDPs are represented in the PRISM language using the [parallel composition](https://www.prismmodelchecker.org/manual/ThePRISMLanguage/ProcessAlgebraOperators). +Accordingly, the size of the representation is linear in the number of NFA transitions and $n$. + +The binary includes two subcommands: + +- `schaeppert convert [n]`: This will read an automaton from the given file and print the n-fold product MDP to standard out. The argument `n` is optional and defaults to 1. +- `schaeppert iterate `: This will read an automaton from the given file and iteratively creates prism-models and calls `prism` on them. + + +```console +./target/release/schaeppert -f dot examples/bottleneck-2-free.dot iterate tmp/ +n=1 -> 1.000 +n=2 -> 1.000 +n=3 -> 0.500 +The value is less than 1.0, stopping the search. +The 3-fold power of this NFA is not controllable. +``` + + diff --git a/examples/guess-empty-n/README.md b/examples/guess-empty-n/README.md new file mode 100644 index 0000000..3e4cf29 --- /dev/null +++ b/examples/guess-empty-n/README.md @@ -0,0 +1,77 @@ + +# A family of uncontrollable NFAs + +Each automaton in this family is a bottleneck of capacity n, facilitated by $O(n)$ states and actions. +The idea is to let all processes go to one of $n+1$ distinct nodes; afterwards, one has to play an action $aj$ +identifying a node $j$ that is not occupied. This move wins, any other loses. + +![nfa for N=2](nfa-2.png) + +All of these models are **not** congrollable for arbitrary number of processes, but the $n$-fold product of the size-$m$ model is controllable iff $n\le m$. + +## Generate model files + +Starting in shepherd's main directory, run the python script to generate 10 models. + +```console +cd examples/guess-empty-n +./generate_models.py 10 + +ls +nfa-10.dot nfa-2.dot nfa-4.dot nfa-6.dot nfa-8.dot generate_models.py +nfa-1.dot nfa-3.dot nfa-5.dot nfa-7.dot nfa-9.dot +``` + +## Solve them using prism + +```console +time ./target/release/schaeppert -f dot examples/guess-empty-n/nfa-5.dot iterate tmp/ +n=1 -> 1.000 +n=2 -> 1.000 +n=3 -> 1.000 +n=4 -> 1.000 +n=5 -> 0.962 +The value is less than 1.0, stopping the search. +The 5-fold power of this NFA is not controllable. + +./target/release/schaeppert -f dot examples/guess-empty-n/nfa-5.dot 5.82s user 0.33s system 216% cpu 2.847 total +``` + + +## Solve them using shepherd + + +```console + time ./target/release/shepherd -f dot examples/guess-empty-n/nfa-5.dot + +Maximal winning strategy; +Answer: + NO (uncontrollable) + +States: ( 3 , 2 , 0 , T , 4 , 1 , 5 ) + Play action 'a' in the downward-closure of + ( _ , _ , 4 , _ , _ , _ , _ ) + + +Play action 'a1' in the downward-closure of + ( ω , ω , _ , _ , ω , _ , ω ) + + +Play action 'a2' in the downward-closure of + ( ω , _ , _ , _ , ω , ω , ω ) + + +Play action 'a3' in the downward-closure of + ( _ , ω , _ , _ , ω , ω , ω ) + + +Play action 'a4' in the downward-closure of + ( ω , ω , _ , _ , _ , ω , ω ) + + +Play action 'a5' in the downward-closure of + ( ω , ω , _ , _ , ω , ω , _ ) + + +./target/release/shepherd -f dot examples/guess-empty-n/nfa-5.dot 0.37s user 0.06s system 175% cpu 0.248 total +``` diff --git a/examples/guess-empty-n/generate_models.py b/examples/guess-empty-n/generate_models.py new file mode 100755 index 0000000..cc03e1c --- /dev/null +++ b/examples/guess-empty-n/generate_models.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 +""" +generate_models.py + +This script generates DOT files representing a family of NFA (nondeterministic finite automaton) models. +For each integer i from 1 to n (inclusive), it creates a file named 'nfa-i.dot' describing an NFA +with a specific structure: +- States 0 through i, plus a target state T. +- Initial transitions from state 0 to each state 1..i+1 labeled 'a'. +- For each state j in 1..i+1, transitions from j to T labeled 'a{k}' for every k in 1..i+1 where k != j. + +Usage: + python generate_models.py +Where is a positive integer. + + +Note: All of these models are not congrollable for arbitrary number of processes, but the n-fold product of the size-m model is controllable by playing action "a", then "aj" if node "j" if none of the n processes reside in node j. +""" +import sys +def generate_model_file(n): + filename = f"nfa-{n}.dot" + with open(filename, "w") as f: + f.write('digraph NFA {\n') + f.write(' rankdir=TB;\n') + f.write(' init [label=" ",shape=none,height=0,width=0];\n') + # States 0..N + for i in range(n + 2): + f.write(f' {i} [label="{i}", shape=circle];\n') + # Target state T + f.write(' T [label="T", shape=doublecircle];\n\n') + # Initial arrow + f.write(' init -> 0;\n') + # Transitions from 0 to 1..N + for i in range(1, n + 2): + f.write(f' 0 -> {i} [label="a"];\n') + # For every two nodes 1 <= i, j <= N with i != j, add j -> T [label="i"] + for j in range(1, n + 2): + for i in range(1, n + 2): + if i != j: + f.write(f' {j} -> T [label="a{i}"];\n') + f.write('}\n') + +def main(): + if len(sys.argv) != 2: + print("Usage: python generate_models.py ") + sys.exit(1) + try: + n = int(sys.argv[1]) + if n < 1: + raise ValueError + except ValueError: + print("n must be a positive integer") + sys.exit(1) + for i in range(1, n + 1): + generate_model_file(i) + +if __name__ == "__main__": + main() diff --git a/examples/guess-empty-n/nfa-2.png b/examples/guess-empty-n/nfa-2.png new file mode 100644 index 0000000..a756590 Binary files /dev/null and b/examples/guess-empty-n/nfa-2.png differ diff --git a/examples/guess-empty-n/run_benchmarks.py b/examples/guess-empty-n/run_benchmarks.py new file mode 100755 index 0000000..f2543df --- /dev/null +++ b/examples/guess-empty-n/run_benchmarks.py @@ -0,0 +1,101 @@ +#!/usr/bin/env python3 +import os +import shutil +import subprocess +import sys +import time +import csv +import matplotlib.pyplot as plt + +# Constants +COMMANDS = [ + "schaeppert -vv -f dot nfa-{n}.dot iterate tmp/", + # "schaeppert -f dot nfa-{n}.dot iterate tmp/", + #"shepherd -vv -f dot nfa-{n}.dot", + "shepherd -f dot nfa-{n}.dot", +] +TIMEOUT = 60 # seconds +TMP_DIR = "tmp" +LOGFILE_TEMPLATE = "log_{cmd_idx}_n{n}.txt" + +def ensure_empty_tmp(): + if os.path.exists(TMP_DIR): + shutil.rmtree(TMP_DIR) + os.makedirs(TMP_DIR) + +def run_command(cmd, logfile, timeout): + start = time.time() + try: + proc = subprocess.run( + cmd, + shell=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + timeout=timeout + ) + elapsed = time.time() - start + with open(logfile, 'wb') as f: + f.write(b'=== STDOUT ===\n') + f.write(proc.stdout) + f.write(b'\n=== STDERR ===\n') + f.write(proc.stderr) + return elapsed + except subprocess.TimeoutExpired as e: + with open(logfile, 'wb') as f: + f.write(b'=== STDOUT ===\n') + f.write((e.stdout or b'')) + f.write(b'\n=== STDERR ===\n') + f.write((e.stderr or b'')) + f.write(b'\n=== TIMEOUT ===\n') + return None + +def main(): + if len(sys.argv) != 2: + print(f"Usage: {sys.argv[0]} N", file=sys.stderr) + sys.exit(1) + try: + N = int(sys.argv[1]) + except ValueError: + print("N must be an integer", file=sys.stderr) + sys.exit(1) + + ensure_empty_tmp() + results = [] + for n in range(1, N + 1): + row = {'n': n} + for cmd_idx, cmd_template in enumerate(COMMANDS): + cmd = cmd_template.format(n=n) + logfile = LOGFILE_TEMPLATE.format(cmd_idx=cmd_idx, n=n) + print(f"Running command {cmd_idx+1} for n={n}: {cmd_template}") + t = run_command(cmd, logfile, TIMEOUT) + if t is None: + break + row[f'cmd{cmd_idx+1}_time'] = t + else: + results.append(row) + continue + break + + # Write CSV to stdout + fieldnames = ['n'] + [f'cmd{i+1}_time' for i in range(len(COMMANDS))] + writer = csv.DictWriter(sys.stdout, fieldnames=fieldnames) + writer.writeheader() + for row in results: + writer.writerow(row) + + # Plotting + ns = [row['n'] for row in results] + for cmd_idx, cmd_template in enumerate(COMMANDS): + times = [row.get(f'cmd{cmd_idx+1}_time') for row in results] + plt.plot(ns, times, label=cmd_template) + plt.xlabel('n') + plt.ylabel('Wallclock time (s)') + plt.legend() + plt.title('Command Running Times') + plt.grid(True) + plt.tight_layout() + plt.savefig("benchmark_results.png") + plt.show() + +if __name__ == '__main__': + main() diff --git a/src/iterative.rs b/src/iterative.rs new file mode 100644 index 0000000..3d4ac82 --- /dev/null +++ b/src/iterative.rs @@ -0,0 +1,265 @@ +use clap::{Parser, Subcommand}; +use log::{debug, info}; +use std::path::{PathBuf,Path}; + +use shepherd::nfa; +mod logging; +use regex::Regex; +use std::fs::File; +use std::io; +use std::io::{Read, Write}; +use std::process::{Command, Stdio}; + +const PRISM_CMD: &str = "prism"; + +#[derive(Parser)] +#[command(version, about, long_about = None)] +pub struct Args { + #[command(subcommand)] + command: Commands, + + #[arg(value_name = "AUTOMATON_FILE", help = "Path to the input")] + pub filename: String, + + #[arg( + short = 'f', + long = "from", + value_enum, + default_value = "tikz", + help = "The input format" + )] + pub input_format: nfa::InputFormat, + + #[arg( + short = 'v', + long = "verbose", + action = clap::ArgAction::Count, + help = "Increase verbosity level" + )] + pub verbosity: u8, + + #[arg( + long, + short = 'l', + value_name = "LOG_FILE", + help = "Optional path to the log file. Defaults to stdout if not specified." + )] + pub log_output: Option, +} + +#[derive(Subcommand)] +enum Commands { + Convert { + /// generate the n-fold product + //#[arg(default_value_t = 2020)] + power : Option, + }, + Iterate { + /// where to write genrated files + tmp_dir: PathBuf, + }, +} + +pub fn main() { + // parse CLI arguments + let cli = Args::parse(); + + // set up logging + logging::setup_logger(cli.verbosity, cli.log_output); + + // parse the input file + let mut nfa = nfa::Nfa::load_from_file( + &cli.filename, + &cli.input_format, + &nfa::StateOrdering::Alphabetical, + ); + + // print the input automaton + info!("{}", nfa); + + // compute the solution + if !nfa.is_complete() { + info!("The automaton is not complete. Completing it..."); + + match nfa.add_state("SINK") { + Ok(sink) => { + info!("Added sink state"); + nfa.complete(Some(sink)); + } + Err(e) => { + info!("Error adding sink state: {}", e); + return; // TODO: handle this error properly + } + } + } + // print the complete automaton again + info!("{}", nfa); + + match &cli.command { + Commands::Iterate { tmp_dir } => { + info!("Iterating prism on the given NFA. Writing files to {tmp_dir:?}."); + let res = iter_prism(&nfa, tmp_dir); + match res { + Ok(stopped_at) =>{ + println!("The {stopped_at}-fold power of this NFA is not controllable."); + } + Err(msg) => { + eprintln!("Error when calling prism:\n {msg}"); + } + } + } + Commands::Convert { power } => { + info!("Converting given NFA to a prism model."); + println!("{}", nfa_to_prism(&nfa, power.unwrap_or(1))) + } + } +} + +fn write_string_to_file(content: &str, file_path: &PathBuf) -> io::Result<()> { + let mut file = File::create(file_path)?; + file.write_all(content.as_bytes())?; + Ok(()) +} + +fn call_prism(args: &[&str]) -> Result { + let mut child = Command::new(PRISM_CMD) + .args(args) + .stdout(Stdio::piped()) + .spawn() + .expect("Failed to call prism"); + + child.wait().expect("failed to wait on child"); + + let stdout = child.stdout.expect("Failed to capture stdout"); + + // Using stdout to read the output from the child process + let mut output = String::new(); + io::BufReader::new(stdout) + .read_to_string(&mut output) + .expect("Failed to read from stdout"); + info!( + "PRISM OUTPUT\n---------------\n{}\n-----------------", + output + ); + + // Compile the regular expression once + let re = Regex::new(r"Value in the initial state: (\d+\.\d+)").unwrap(); + + for line in output.lines() { + if let Some(captures) = re.captures(line) { + if let Some(value) = captures.get(1) { + return value.as_str().parse::().map_err(|_| ()); + } + } + } + Err(()) +} + +fn nfa_to_prism(nfa: &nfa::Nfa, n: u16) -> String { + let mut prism_input = String::new(); + prism_input.push_str("mdp\n\n"); + + // module M1 will be our NFA. + prism_input.push_str("module M1\n"); + + // we assume that there is only one initial state; get it. + let initial: nfa::State = nfa.initial_states().iter().cloned().next().unwrap(); + + // define states string for prism + prism_input.push_str(&format!( + "s1 : [0..{}] init {initial};\n", + nfa.nb_states() - 1 + )); + + // define transitions + for (act, am) in nfa.get_edges().iter() { + // for every alphabet letter + for src in 0..am.dim() { + // for all states + let succs = am.get_successors(src); // get successors + // prism requires explicit floating point numbers to represent distributions. + // here we represent a uniform dist among successors. + let prob = 1.0 / succs.len() as f64; + let update = succs + .iter() + .map(|trg| format!("{}:(s1'={})", prob, trg)) + .collect::>() + .join(" + "); + prism_input.push_str(&format!("[{act}] s1={} -> {};\n", src, update)); + } + } + prism_input.push_str("endmodule\n\n"); + + // Add a copy of the MDP for every power up to n + for i in 2..=n { + prism_input.push_str(&format!("module M{i} = M1 [s1=s{i}, s{i}=s1] endmodule\n")); + } + + // define a label representing global reachability target: + // every component is in one of its final states. + let mut final_line = String::from("\nlabel \"final\" = "); + + let mut conj = Vec::new(); + for i in 1..=n { + conj.push( + nfa.final_states() + .iter() + .map(|f| format!("s{i}={f}")) + .collect::>() + .join("| "), + ); + } + final_line.push_str( + &conj + .iter() + .map(|f| format!("( {f} )")) + .collect::>() + .join(" & "), + ); + final_line.push_str(";\n"); + prism_input.push_str(&final_line); + + // define the global system as the product of all n many copies. + // This uses prisms syntax for parallel composition. + prism_input.push_str("\nsystem\n"); + let prod_string = (1..=n) + .map(|i| format!("M{i}")) + .collect::>() + .join(" || "); + prism_input.push_str(&prod_string); + prism_input.push_str("\nendsystem\n"); + + prism_input +} + +fn iter_prism(nfa: &nfa::Nfa, tmp_dir :&Path) -> Result { + let mut i: u16 = 1; + loop { + // create prism input string + let prism_model = nfa_to_prism(nfa, i); + debug!("{}", prism_model); + + // write prism input to file + let prism_input_path = tmp_dir.join(format!("model-{}.pm", i)); + if let Err(e) = write_string_to_file(&prism_model, &prism_input_path) { + return Err(format!("Error writing to file: {}", e)); + } + + info!("Wrote prism input to file: {}", prism_input_path.display()); + + let value = call_prism(&[ + "-pf", + "Pmax=? [ F \"final\" ]", + &prism_input_path.to_string_lossy(), + ]) + .unwrap(); + println!("n={} -> {:.3}", i, value); + + if value < 1.0 { + println!("The value is less than 1.0, stopping the search."); + break; + } + i += 1; + } + Ok(i) +} diff --git a/src/nfa.rs b/src/nfa.rs index 8df8bb4..b46ecb9 100644 --- a/src/nfa.rs +++ b/src/nfa.rs @@ -320,6 +320,19 @@ impl Nfa { }); } + /// add a new state with given label + /// returns the index of the new state or an Error if the state already exists + pub fn add_state(&mut self, label: &str) -> Result{ + let label = String::from(label); + match self.states.contains(&label) { + true => Err("state with label exists"), + false => { + self.states.push(label); + Ok(self.states.len() - 1) + } + } + } + fn check_state(&self, q: State) { assert!(q < self.nb_states(), "State {} is not in the NFA", q) }