From 9dc32ddc119179ad2bc51ddccc8ea0bf7f34f6aa Mon Sep 17 00:00:00 2001
From: Joshua Moerman <joshua.moerman@ou.nl>
Date: Fri, 17 May 2024 20:58:30 +0200
Subject: [PATCH] New SAT solving script which optimally decomposes an FSM on
 outputs

---
 other/decompose_fsm.py                        | 202 ++++++++++++++++++
 ...ealy.py => decompose_observation_table.py} |   2 +-
 2 files changed, 203 insertions(+), 1 deletion(-)
 create mode 100644 other/decompose_fsm.py
 rename other/{decompose_mealy.py => decompose_observation_table.py} (99%)

diff --git a/other/decompose_fsm.py b/other/decompose_fsm.py
new file mode 100644
index 0000000..f124c3d
--- /dev/null
+++ b/other/decompose_fsm.py
@@ -0,0 +1,202 @@
+from pysat.solvers import Solver
+from pysat.card import CardEnc
+from pysat.formula import IDPool
+from pysat.formula import CNF
+
+### Gebruik:
+# Stap 1: pip3 install python-sat
+# Stap 2: python3 decompose_fsm.py
+
+
+###################################
+# Wat dingetjes over Mealy machines
+
+# Voorbeeld: 2n states, input-alfabet 'a' en 'b', outputs [0...n-1]
+def rick_koenders_machine(N):
+  transition_fun = {((n, 0), 'a') : ((n+1) % N, 0) for n in range(N)}
+  transition_fun |= {(((n+1) % N, 1), 'a') : (n % N, 1) for n in range(N)}
+  transition_fun |= {((n, b), 'b') : (n, not b) for b in [0, 1] for n in range(N)}
+  output_fun = {((n, b), 'a') : n for b in [0, 1] for n in range(N)}
+  output_fun |= {((n, b), 'b') : 0 for b in [0, 1] for n in range(N)}
+  initial_state = (0, 0)
+  inputs = ['a', 'b']
+  outputs = [n for n in range(N)]
+  return {'transition_fun': transition_fun, 'output_fun': output_fun, 'initial_state': initial_state, 'inputs': inputs, 'outputs': outputs}
+
+def print_table(cell, rs, cs):
+  first_col_size = max([len(str(r)) for r in rs])
+  col_size = 1 + max([len(str(c)) for c in cs])
+
+  print(''.rjust(first_col_size), end='')
+  for c in cs:
+    print(str(c).rjust(col_size), end='')
+  print('')
+
+  for r in rs:
+    print(str(r).rjust(first_col_size), end='')
+    for c in cs:
+      print(cell(r, c).rjust(col_size), end='')
+    print('')
+
+
+################
+# Voorbeeld data
+N = 8  # fsm heeft 2N states
+machine = rick_koenders_machine(N)
+
+states = [(n, b) for n in range(N) for b in [0, 1]]
+
+# We zoeken 2 componenten met gezamelijke grootte 6 (minder dan 8)
+# als de de total_size te laag is => UNSAT => duurt lang
+c = 2
+total_size = 10
+
+
+########################
+# Encodering naar logica
+print('Start encoding')
+os = machine['outputs']      # outputs
+rids = [i for i in range(c)] # components
+vpool = IDPool()
+cnf = CNF()
+
+# Een hulp variabele voor False en True, maakt de andere variabelen eenvoudiger
+def var_const(b):
+  return(vpool.id(('const', b)))
+
+cnf.append([var_const(True)])
+cnf.append([-var_const(False)])
+
+# Voor elke relatie en elke twee elementen o1 en o2, is er een variabele die
+# aangeeft of o1 en o2 gerelateerd zijn. Er is 1 variabele voor xRy en yRx, dus
+# symmetrie is al ingebouwd. Reflexiviteit is ook ingebouwd.
+def var_rel(rid, o1, o2):
+  if o1 == o2:
+    return var_const(True)
+
+  [so1, so2] = sorted([o1, o2])
+  return(vpool.id(('rel', rid, so1, so2)))
+
+# De relatie op outputs geeft een relaties op states. Deze relatie moet ook een
+# bisimulatie zijn.
+def var_state_rel(rid, s1, s2):
+  if s1 == s2:
+    return var_const(True)
+
+  [ss1, ss2] = sorted([s1, s2])
+  return(vpool.id(('state_rel', rid, ss1, ss2)))
+
+# Voor elke relatie, en elke equivalentie-klasse, kiezen we precies 1 state
+# als representant. Deze variabele geeft aan welk element.
+def var_state_rep(rid, s):
+  return(vpool.id(('state_rep', rid, s)))
+
+# Contraints zodat de relatie een equivalentie relatie is. We hoeven alleen
+# maar transitiviteit te encoderen, want refl en symm zijn ingebouwd in de var.
+# Dit stukje van het encoderen duurt het langst.
+print('- transitivity (o)')
+for rid in rids:
+  for xo in os:
+    for yo in os:
+      for zo in os:
+        # als xo R yo en yo R zo dan xo R zo
+        cnf.append([-var_rel(rid, xo, yo), -var_rel(rid, yo, zo), var_rel(rid, xo, zo)])
+
+print('- transitivity (s)')
+for rid in rids:
+  for sx in states:
+    for sy in states:
+      for sz in states:
+        # als sx R sy en sy R sz dan sx R sz
+        cnf.append([-var_state_rel(rid, sx, sy), -var_state_rel(rid, sy, sz), var_state_rel(rid, sx, sz)])
+
+# Constraint zodat de relaties samen alle elementen kunnen onderscheiden.
+# (Aka: the bijbehorende quotienten zijn joint-injective.)
+print('- injectivity (o)')
+for xi, xo in enumerate(os):
+  for yo in os[xi+1:]:
+    # Tenminste een rid moet een verschil maken
+    cnf.append([-var_rel(rid, xo, yo) for rid in rids])
+
+# sx ~ sy => for each input: (1) outputs equivalent AND (2) successors related
+# Momenteel hebben we niet de inverse implicatie, is misschien ook niet nodig?
+print('- bisimulation modulo rel')
+for rid in rids:
+  for sx in states:
+    for sy in states:
+      for i in machine['inputs']:
+        # sx ~ sy => output(sx, i) ~ output(sy, i)
+        ox = machine['output_fun'][(sx, i)]
+        oy = machine['output_fun'][(sy, i)]
+        cnf.append([-var_state_rel(rid, sx, sy), var_rel(rid, ox, oy)])
+
+        # sx ~ sy => delta(sx, i) ~ delta(sy, i)
+        tx = machine['transition_fun'][(sx, i)]
+        ty = machine['transition_fun'][(sy, i)]
+        cnf.append([-var_state_rel(rid, sx, sy), var_state_rel(rid, tx, ty)])
+
+# De constraints die zorgen dat representanten ook echt representanten zijn.
+print('- representatives')
+for rid in rids:
+  for ix, sx in enumerate(states):
+    # Belangrijkste: een element is een representant, of equivalent met een
+    # eerder element. We forceren hiermee dat de solver representanten moet
+    # kiezen (voor aan de lijst).
+    cnf.append([var_state_rep(rid, sx)] + [var_state_rel(rid, sx, sy) for sy in states[:ix]] )
+    for sy in states[:ix]:
+      # rx en ry kunnen niet beide een representant zijn, tenzij ze
+      # niet gerelateerd zijn.
+      cnf.append([-var_state_rep(rid, sx), -var_state_rep(rid, sy), -var_state_rel(rid, sx, sy)])
+
+# Tot slot willen we weinig representanten. Dit doen we met een "atmost"
+# formule. Idealiter zoeken we naar de total_size, maar die staat nu vast.
+print('- representatives at most k')
+cnf_optim = CardEnc.atmost([var_state_rep(rid, sx) for rid in rids for sx in states], total_size, vpool=vpool)
+cnf.extend(cnf_optim)
+
+def print_eqrel(rel, xs):
+  print_table(lambda r, c: 'Y' if rel(r, c) else '·', xs, xs)
+
+
+##################################
+# Probleem oplossen met solver :-)
+print('Start solving')
+print('- copying formula')
+with Solver(bootstrap_with=cnf) as solver:
+  print('- actual solve')
+  sat = solver.solve()
+
+  if not sat:
+    print('unsat :-(')
+    exit()
+
+  print(f'sat :-)')
+
+  # Even omzetten in een makkelijkere data structuur
+  print('- get model')
+  m = solver.get_model()
+  model = {}
+  for l in m:
+    if l < 0: model[-l] = False
+    else: model[l] = True
+
+  for rid in rids:
+    print(f'Relation {rid}:')
+    print_eqrel(lambda x, y: model[var_rel(rid, x, y)], os)
+
+  for rid in rids:
+    print(f'State relation {rid}:')
+    print_eqrel(lambda x, y: model[var_state_rel(rid, x, y)], states)
+
+  # print equivalence classes
+  count = 0
+  for rid in rids:
+    print(f'component {rid}')
+    # Eerst verzamelen we de representanten
+    for s in states:
+      if model[var_state_rep(rid, s)]:
+        print(f'- representative state {s}')
+        count += 1
+
+  # count moet gelijk zijn aan cost
+  print(f'total size = {count}')
diff --git a/other/decompose_mealy.py b/other/decompose_observation_table.py
similarity index 99%
rename from other/decompose_mealy.py
rename to other/decompose_observation_table.py
index b3129f9..9fa6d09 100644
--- a/other/decompose_mealy.py
+++ b/other/decompose_observation_table.py
@@ -5,7 +5,7 @@ from pysat.formula import CNF
 
 ### Gebruik:
 # Stap 1: pip3 install python-sat
-# Stap 2: python3 decompose_mealy.py
+# Stap 2: python3 decompose_observation_table.py
 
 
 ###################################