Iterative

Cargo.toml

@@ -21,3 +21,8 @@ itertools = "0.14.0"
 
 #profile mode
 #debug = true
+
+[[bin]]
+name = "schaeppert"
+path = "src/iterative.rs"
+

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 <AUTOMATON_FILE> 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 <AUTOMATON_FILE> iterate <TMP_DIR>`: 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.
+```
+
+
 

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
+```

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 <n>
+Where <n> 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 <n>")
+        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()

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()