mirror of
https://git.cs.ou.nl/joshua.moerman/mealy-decompose.git
synced 2025-04-29 17:57:44 +02:00
145 lines
4.8 KiB
Python
145 lines
4.8 KiB
Python
# Copyright 2024-2025 Joshua Moerman, Open Universiteit. All rights reserved
|
||
# SPDX-License-Identifier: EUPL-1.2
|
||
|
||
from pysat.solvers import Solver
|
||
from pysat.card import CardEnc
|
||
from pysat.formula import IDPool
|
||
from pysat.formula import CNF
|
||
|
||
|
||
# Een verzameling X ontbinden in factoren X1 ... Xc, zodat X ⊆ X1 × ... × Xc.
|
||
# Hierbij is c een parameter (het aantal componenten), en ook het aantal
|
||
# elementen van X1 t/m Xc moet vooraf bepaald zijn, dat is 'total_size'.
|
||
|
||
|
||
# Voorbeeld data
|
||
# snel voorbeeld: n = 27, c = 3 en total_size = 9
|
||
# langzaam vb.: n = 151, c = 4 en total_size = 15
|
||
n = 151
|
||
c = 4
|
||
total_size = 15 # als deze te laag is => UNSAT => duurt lang
|
||
|
||
os = [i for i in range(n)] # outputs
|
||
rids = [i for i in range(c)] # components
|
||
|
||
# Optimale decompositie met slechts verzamelingen. Zie ook A007600 (omkering
|
||
# A000792) in de OEIS, dit groeit met ongeveer O(log(n)). Ik weet niet hoeveel
|
||
# c moet(/mag) zijn. Ook relevant: Katona's problem en Edmonds' problem.
|
||
#
|
||
# Tabel: bij elke n de optimale totale grootte (en de minimale c daarbij)
|
||
# n = 1 2 3 4 5 6 ≤9 ≤12 ≤16 ≤18 ≤29 ≤27 ≤36 ≤48 ≤54 ≤64 ≤81 ≤108 ≤144 ≤111
|
||
# -------------------------------------------------------------------------
|
||
# size = 1 2 3 4 5 5 6 7 8 8 9 9 10 11 11 12 12 13 14 15
|
||
# c = 1 1 1 1 1 2 2 2 2 3 2 3 3 3 4 3 4 4 4 5
|
||
|
||
|
||
print('Start encoding')
|
||
vpool = IDPool()
|
||
cnf = CNF()
|
||
|
||
|
||
# Een hulp variabele voor False en True, maakt de andere variabelen eenvoudiger
|
||
def var_const(b):
|
||
return vpool.id(('bool', 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(('var_rel', rid, so1, so2))
|
||
|
||
|
||
# Voor elke relatie, en elke equivalentie-klasse, kiezen we precies 1 element
|
||
# als representant. Deze variabele geeft aan welk element.
|
||
def var_rep(rid, o):
|
||
return vpool.id(('var_rep', rid, o))
|
||
|
||
|
||
# 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')
|
||
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)])
|
||
|
||
# Constraint zodat de relaties samen alle elementen kunnen onderscheiden.
|
||
# (Aka: the bijbehorende quotienten zijn joint-injective.)
|
||
print('- injectivity')
|
||
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])
|
||
|
||
# De constraints die zorgen dat representanten ook echt representanten zijn.
|
||
print('- representatives')
|
||
for rid in rids:
|
||
for xi, xo in enumerate(os):
|
||
# Belangrijkste: een element is een representant, of equivalent met een
|
||
# later element. We forceren hiermee dat de solver representanten moet
|
||
# kiezen (achter aan de lijst).
|
||
cnf.append([var_rep(rid, xo)] + [var_rel(rid, xo, yo) for yo in os[xi + 1 :]])
|
||
for _, yo in enumerate(os[xi + 1 :], xi + 1):
|
||
# xo en yo kunnen niet beide een representant zijn, tenzij ze
|
||
# niet gerelateerd zijn.
|
||
cnf.append([-var_rep(rid, xo), -var_rep(rid, yo), -var_rel(rid, xo, yo)])
|
||
|
||
# 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_rep(rid, xo) for rid in rids for xo in os], total_size, vpool=vpool)
|
||
cnf.extend(cnf_optim)
|
||
|
||
|
||
# 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('sat :-)')
|
||
|
||
# Even omzetten in een makkelijkere data structuur
|
||
print('- get model')
|
||
m = solver.get_model()
|
||
model = {}
|
||
for l in m:
|
||
model[abs(l)] = l > 0
|
||
|
||
# print equivalence classes
|
||
count = 0
|
||
for rid in rids:
|
||
# Eerst verzamelen we de representanten
|
||
reps = []
|
||
for xo in os:
|
||
if model[var_rep(rid, xo)]:
|
||
reps.append(xo)
|
||
count += 1
|
||
|
||
# Dan zoeken we wat bij de representant hoort en printen dat.
|
||
for r in reps:
|
||
rep_class = [r]
|
||
for yo in os:
|
||
if yo != r and model[var_rel(rid, r, yo)]:
|
||
rep_class.append(yo)
|
||
print(f'{rid} => class of {r} = {sorted(rep_class)}')
|
||
|
||
# count moet gelijk zijn aan cost
|
||
print(f'total size = {count}')
|