Project: Моделирование труднорешаемых задач

Path: hard-problems-formulations / minimum-exact-cover.ipynb

Views: ²²

Visibility: Unlisted (only visible to those who know the link)

Image: ubuntu2204

Kernel: .venv

MINIMUM EXACT COVER / Покрытие множеств минимальным объемом

In [1]:

import random
from IPython.core.display import SVG
import pyomo.environ as pyo
from pysat.solvers import Solver
from pysat.formula import CNF 
import py_svg_combinatorics as psc
from ipywidgets import widgets, HBox
from collections import Counter
from pprint import pprint
from random import randint
import numpy as np
from IPython.display import IFrame
import IPython
from copy import copy
import os
from pathlib import Path
nbname = ''
try:
    nbname = __vsc_ipynb_file__
except:
    if 'COCALC_JUPYTER_FILENAME' in os.environ:
        nbname = os.environ['COCALC_JUPYTER_FILENAME']    
title_ = Path(nbname).stem.replace('-', '_').title()    
IFrame(f'https://discopal.ispras.ru/index.php?title=Hardprob/{title_}&useskin=cleanmonobook', width=1280, height=300)

Постановка задачи

Представим, что мы имеем конечное множество «U» (universum) и «S» — список подмножеств множества «U».

Покрытием называют семейство $\mathcal{C}\subseteq\mathcal{S}$ наименьшей мощности, объединением которых является $\mathcal{U}$ . В случае постановки вопроса о разрешении на вход подаётся пара $(\mathcal{U},\mathcal{S})$ и целое число k; вопросом является существование покрывающего множества мощности k (или менее).

Оптимизационная версия задачи, '''минимальное покрытие множеств''', задаёт вопрос о минимальном объеме покрытия, измеряемого не в числе множеств из «S» покрывающих «U», а в числе элементов в этих подмножествах

In [2]:

sets = [
        [f"{i:02}" for i in range(1,18)],
        [f"{2*i:02}" for i in range(1,10)]+['01'],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 1)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 2)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 3)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 4)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 5)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 6)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 7)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 8)],
        [f"{i:02}" for i in np.random.choice(np.arange(1,20), 9)],
]

In [3]:

svg = psc.subsets2svg(sets)
SVG(data=svg)

Реализация в Pyomo

In [4]:

def print_solution(m):
    for v in m.component_data_objects(pyo.Var):
        if v.value and v.value > 0:
            print(str(v), v.value)

In [5]:

def get_model(sets):
    sets_ = sets
    if isinstance(sets_, dict):
        sets_ = [val for _, val in sets_.items()]

    model = pyo.ConcreteModel()
    model.S = [set(str(item) for item in sublist) for sublist in sets_]
    model.m = len(model.S)
    model.J = range(model.m)
    model.U = sorted(set([str(item) for sublist in sets_ for item in sublist]))
    model.n = len(model.U)

    #Включаем ли это подмножество.
    model.x = pyo.Var(model.J, domain=pyo.Binary)
    model.set_count = pyo.Objective(expr = sum( len(model.S[j]) * model.x[j] for j in model.J), sense=pyo.minimize)

    @model.Constraint(model.U)
    def каждый_элемент_покрыт_хотябы_одним_множеством(m, u):
        return sum(m.x[j] for j in m.J if u in m.S[j]) >= 1

    return model


m = get_model(sets)
m.U

['01',
 '02',
 '03',
 '04',
 '05',
 '06',
 '07',
 '08',
 '09',
 '10',
 '11',
 '12',
 '13',
 '14',
 '15',
 '16',
 '17',
 '18',
 '19']

In [6]:

m.S

[{'01',
  '02',
  '03',
  '04',
  '05',
  '06',
  '07',
  '08',
  '09',
  '10',
  '11',
  '12',
  '13',
  '14',
  '15',
  '16',
  '17'},
 {'01', '02', '04', '06', '08', '10', '12', '14', '16', '18'},
 {'10'},
 {'11', '14'},
 {'02', '08', '18'},
 {'16', '17', '18', '19'},
 {'01', '07', '08', '11', '17'},
 {'01', '05', '07', '12', '18', '19'},
 {'04', '06', '12', '14', '18', '19'},
 {'02', '04', '05', '09', '10', '15', '18'},
 {'04', '06', '07', '09', '11', '17', '18', '19'}]

In [7]:

m.pprint()

1 Var Declarations
    x : Size=11, Index={0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :     1 : False :  True : Binary
          1 :     0 :  None :     1 : False :  True : Binary
          2 :     0 :  None :     1 : False :  True : Binary
          3 :     0 :  None :     1 : False :  True : Binary
          4 :     0 :  None :     1 : False :  True : Binary
          5 :     0 :  None :     1 : False :  True : Binary
          6 :     0 :  None :     1 : False :  True : Binary
          7 :     0 :  None :     1 : False :  True : Binary
          8 :     0 :  None :     1 : False :  True : Binary
          9 :     0 :  None :     1 : False :  True : Binary
         10 :     0 :  None :     1 : False :  True : Binary

1 Objective Declarations
    set_count : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 17*x[0] + 10*x[1] + x[2] + 2*x[3] + 3*x[4] + 4*x[5] + 5*x[6] + 6*x[7] + 6*x[8] + 7*x[9] + 8*x[10]

1 Constraint Declarations
    каждый_элемент_покрыт_хотябы_одним_множеством : Size=19, Index={01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19}, Active=True
        Key : Lower : Body                                            : Upper : Active
         01 :   1.0 :                       x[0] + x[1] + x[6] + x[7] :  +Inf :   True
         02 :   1.0 :                       x[0] + x[1] + x[4] + x[9] :  +Inf :   True
         03 :   1.0 :                                            x[0] :  +Inf :   True
         04 :   1.0 :               x[0] + x[1] + x[8] + x[9] + x[10] :  +Inf :   True
         05 :   1.0 :                              x[0] + x[7] + x[9] :  +Inf :   True
         06 :   1.0 :                      x[0] + x[1] + x[8] + x[10] :  +Inf :   True
         07 :   1.0 :                      x[0] + x[6] + x[7] + x[10] :  +Inf :   True
         08 :   1.0 :                       x[0] + x[1] + x[4] + x[6] :  +Inf :   True
         09 :   1.0 :                             x[0] + x[9] + x[10] :  +Inf :   True
         10 :   1.0 :                       x[0] + x[1] + x[2] + x[9] :  +Inf :   True
         11 :   1.0 :                      x[0] + x[3] + x[6] + x[10] :  +Inf :   True
         12 :   1.0 :                       x[0] + x[1] + x[7] + x[8] :  +Inf :   True
         13 :   1.0 :                                            x[0] :  +Inf :   True
         14 :   1.0 :                       x[0] + x[1] + x[3] + x[8] :  +Inf :   True
         15 :   1.0 :                                     x[0] + x[9] :  +Inf :   True
         16 :   1.0 :                              x[0] + x[1] + x[5] :  +Inf :   True
         17 :   1.0 :                      x[0] + x[5] + x[6] + x[10] :  +Inf :   True
         18 :   1.0 : x[1] + x[4] + x[5] + x[7] + x[8] + x[9] + x[10] :  +Inf :   True
         19 :   1.0 :                      x[5] + x[7] + x[8] + x[10] :  +Inf :   True

3 Declarations: x set_count каждый_элемент_покрыт_хотябы_одним_множеством

In [8]:


ilp_solver = pyo.SolverFactory('cbc')
ilp_solver.solve(m).write()
print_solution(m)

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 21.0
  Upper bound: 21.0
  Number of objectives: 1
  Number of constraints: 2
  Number of variables: 7
  Number of binary variables: 11
  Number of integer variables: 11
  Number of nonzeros: 7
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.019840240478515625
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x[0] 1.0
x[5] 1.0

In [9]:

print_solution(m)

x[0] 1.0
x[5] 1.0

In [10]:

selected = [j for j in m.J if m.x[j].value > 0]
selected

[0, 5]

In [11]:

SVG(psc.subsets2svg(sets, selected))

In [12]:

cnf3 = CNF(from_clauses=psc.rand3cnf(1, 3))
cnf3.clauses

[(-2, -3, 1)]

In [13]:

from pysat.solvers import Solver

solver = Solver(bootstrap_with=cnf3)
res = solver.solve()
res

True

In [14]:

print(solver.get_model())

[-1, -2, -3]

In [15]:

cnf3.clauses

[(-2, -3, 1)]

In [16]:

from collections import defaultdict
def cnf32setcovering(cnf):
    '''
    Будем считать, что на входе 3SAT формула без повторов литералов в дизьюнкциях и т.п.
    '''
    n = cnf.nv
    m = len(cnf.clauses)
    S = {}

    for j in range(1, n+1):
        S[f'x_{{{j}}}=0'] = [f'x_{{{j}}}']
        S[f'x_{{{j}}}=1'] = [f'x_{{{j}}}']

    for i in range(m):
        U_clause = f'''C_{{{i}}}'''
        clause = cnf.clauses[i]
        for lit in clause:
            j = abs(lit)
            needvar = f'x_{{{j}}}=1'
            if lit < 0:
                needvar = f'x_{{{j}}}=0'
            S[needvar].append(U_clause)
    # Набиваем мусором, чтобы выровнять множества по размеру        
    max_size = max([len(s) for s in S.values()])        
    max_garbage_items = 0
    for s in S:
        add_garbage_items = max_size - len(S[s])
        max_garbage_items = max(max_garbage_items, add_garbage_items)
        for i in range(add_garbage_items):        
            S[s].append(f'G{i+1}')
    S['G'] = [f'G{i}' for i in range(max_garbage_items+1)]
    return S

In [17]:

cnf3 = CNF(from_clauses=psc.rand3cnf(1, 3))
cnf3.clauses

[(2, 3, 1)]

In [18]:

sets = cnf32setcovering(cnf3)
pprint(sets)

{'G': ['G0', 'G1'],
 'x_{1}=0': ['x_{1}', 'G1'],
 'x_{1}=1': ['x_{1}', 'C_{0}'],
 'x_{2}=0': ['x_{2}', 'G1'],
 'x_{2}=1': ['x_{2}', 'C_{0}'],
 'x_{3}=0': ['x_{3}', 'G1'],
 'x_{3}=1': ['x_{3}', 'C_{0}']}

In [19]:

SVG(psc.subsets2svg(sets))

In [20]:

m = get_model(sets)
m.U, m.S

(['C_{0}', 'G0', 'G1', 'x_{1}', 'x_{2}', 'x_{3}'],
 [{'G1', 'x_{1}'},
  {'C_{0}', 'x_{1}'},
  {'G1', 'x_{2}'},
  {'C_{0}', 'x_{2}'},
  {'G1', 'x_{3}'},
  {'C_{0}', 'x_{3}'},
  {'G0', 'G1'}])

In [21]:

# ilp_solver = pyo.SolverFactory('cbc')
ilp_solver.solve(m).write()
print_solution(m)

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 8.0
  Upper bound: 8.0
  Number of objectives: 1
  Number of constraints: 4
  Number of variables: 6
  Number of binary variables: 7
  Number of integer variables: 7
  Number of nonzeros: 6
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.026528120040893555
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x[1] 1.0
x[2] 1.0
x[4] 1.0
x[6] 1.0

In [22]:

m.x.extract_values()

{0: 0.0, 1: 1.0, 2: 1.0, 3: 0.0, 4: 1.0, 5: 0.0, 6: 1.0}

In [23]:

selected = [j for j in m.J if m.x.extract_values()[j] > 0]
selected

[1, 2, 4, 6]

In [24]:

cnf3.clauses

[(2, 3, 1)]

In [25]:

pprint(sets)

{'G': ['G0', 'G1'],
 'x_{1}=0': ['x_{1}', 'G1'],
 'x_{1}=1': ['x_{1}', 'C_{0}'],
 'x_{2}=0': ['x_{2}', 'G1'],
 'x_{2}=1': ['x_{2}', 'C_{0}'],
 'x_{3}=0': ['x_{3}', 'G1'],
 'x_{3}=1': ['x_{3}', 'C_{0}']}

In [26]:

SVG(psc.subsets2svg(sets, selected))

In [27]:

cnf3 = CNF(from_clauses=psc.rand3cnf(2, 3))
cnf3.clauses

[(-2, -3, -1), (-3, 1, 2)]

In [28]:

sets = cnf32setcovering(cnf3)
pprint(sets)

{'G': ['G0', 'G1', 'G2'],
 'x_{1}=0': ['x_{1}', 'C_{0}', 'G1'],
 'x_{1}=1': ['x_{1}', 'C_{1}', 'G1'],
 'x_{2}=0': ['x_{2}', 'C_{0}', 'G1'],
 'x_{2}=1': ['x_{2}', 'C_{1}', 'G1'],
 'x_{3}=0': ['x_{3}', 'C_{0}', 'C_{1}'],
 'x_{3}=1': ['x_{3}', 'G1', 'G2']}

In [29]:

SVG(psc.subsets2svg(sets))

In [30]:

m = get_model(sets)
m.U, m.S

(['C_{0}', 'C_{1}', 'G0', 'G1', 'G2', 'x_{1}', 'x_{2}', 'x_{3}'],
 [{'C_{0}', 'G1', 'x_{1}'},
  {'C_{1}', 'G1', 'x_{1}'},
  {'C_{0}', 'G1', 'x_{2}'},
  {'C_{1}', 'G1', 'x_{2}'},
  {'C_{0}', 'C_{1}', 'x_{3}'},
  {'G1', 'G2', 'x_{3}'},
  {'G0', 'G1', 'G2'}])

In [31]:

ilp_solver.solve(m).write()
print_solution(m)

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 12.0
  Upper bound: 12.0
  Number of objectives: 1
  Number of constraints: 5
  Number of variables: 6
  Number of binary variables: 7
  Number of integer variables: 7
  Number of nonzeros: 6
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.01924443244934082
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
x[1] 1.0
x[2] 1.0
x[5] 1.0
x[6] 1.0

In [32]:

selected = [j for j in m.J if m.x[j].value > 0]
selected

[1, 2, 5, 6]

In [33]:

cnf3.clauses

[(-2, -3, -1), (-3, 1, 2)]

In [34]:

pprint(sets)

{'G': ['G0', 'G1', 'G2'],
 'x_{1}=0': ['x_{1}', 'C_{0}', 'G1'],
 'x_{1}=1': ['x_{1}', 'C_{1}', 'G1'],
 'x_{2}=0': ['x_{2}', 'C_{0}', 'G1'],
 'x_{2}=1': ['x_{2}', 'C_{1}', 'G1'],
 'x_{3}=0': ['x_{3}', 'C_{0}', 'C_{1}'],
 'x_{3}=1': ['x_{3}', 'G1', 'G2']}

In [35]:

SVG(psc.subsets2svg(sets, selected))

In [36]:

m.pprint()
m.U

1 Var Declarations
    x : Size=7, Index={0, 1, 2, 3, 4, 5, 6}
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :   0.0 :     1 : False : False : Binary
          1 :     0 :   1.0 :     1 : False : False : Binary
          2 :     0 :   1.0 :     1 : False : False : Binary
          3 :     0 :   0.0 :     1 : False : False : Binary
          4 :     0 :   0.0 :     1 : False : False : Binary
          5 :     0 :   1.0 :     1 : False : False : Binary
          6 :     0 :   1.0 :     1 : False : False : Binary

1 Objective Declarations
    set_count : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : 3*x[0] + 3*x[1] + 3*x[2] + 3*x[3] + 3*x[4] + 3*x[5] + 3*x[6]

1 Constraint Declarations
    каждый_элемент_покрыт_хотябы_одним_множеством : Size=8, Index={C_{0}, C_{1}, G0, G1, G2, x_{1}, x_{2}, x_{3}}, Active=True
        Key   : Lower : Body                                    : Upper : Active
        C_{0} :   1.0 :                      x[0] + x[2] + x[4] :  +Inf :   True
        C_{1} :   1.0 :                      x[1] + x[3] + x[4] :  +Inf :   True
           G0 :   1.0 :                                    x[6] :  +Inf :   True
           G1 :   1.0 : x[0] + x[1] + x[2] + x[3] + x[5] + x[6] :  +Inf :   True
           G2 :   1.0 :                             x[5] + x[6] :  +Inf :   True
        x_{1} :   1.0 :                             x[0] + x[1] :  +Inf :   True
        x_{2} :   1.0 :                             x[2] + x[3] :  +Inf :   True
        x_{3} :   1.0 :                             x[4] + x[5] :  +Inf :   True

3 Declarations: x set_count каждый_элемент_покрыт_хотябы_одним_множеством

['C_{0}', 'C_{1}', 'G0', 'G1', 'G2', 'x_{1}', 'x_{2}', 'x_{3}']

Генерация теcтов

In [37]:

def solution_of_ilp_to_3sat(m):
    ''' 
    Конвертируем ЦЛП-решение сведенной из 3SAT к MEC задачи, к вопросу, выполнима ли 3CNF-формула?
    '''
    res = sum(m.x.extract_values().values())-1 == (len(m.x)-1)/2
    return res

In [38]:

from pyomo.opt import SolverStatus, TerminationCondition

def generate_tests(clauses_amount):
    test_cnf3 = CNF(from_clauses=psc.rand3cnf(clauses_amount, 3))
    
    ilp_test_sets = cnf32setcovering(test_cnf3)
    ilp_test_model = get_model(ilp_test_sets)
    ilp_test_result_struct = ilp_solver.solve(ilp_test_model)
    ilp_test_result = solution_of_ilp_to_3sat(ilp_test_model)    
#     ilp_test_result = (ilp_test_result_struct.solver.status == SolverStatus.ok and ilp_test_result_struct.solver.termination_condition == TerminationCondition.optimal)
#     ilp_test_no_result = (ilp_test_result_struct.solver.termination_condition == TerminationCondition.infeasible)
    
    py3sat_test_solver = Solver(bootstrap_with=test_cnf3)
    py3sat_test_result = py3sat_test_solver.solve()
    
    if (py3sat_test_result != ilp_test_result):
#         print("For clauses amount = ", clauses_amount, "solutions found: for Pyomo: ", ilp_test_result, ", if it is infeasible: ", ilp_test_no_result, ", for 3-SAT: ", py3sat_test_result)
        print("Fail!", test_cnf3, py3sat_test_result, ilp_test_result)
        raise Exception('WTF!')
#         return False#, ilp_test_sets, ilp_test_model, py3sat_test_solver.get_model()
    else:
#         print("Both solutions are found for clauses amount = ", clauses_amount, ", result is: ", py3sat_test_result)
        return True#, ilp_test_sets, ilp_test_model, py3sat_test_solver.get_model()

In [39]:

#result, got_sets, ilp_model, py3sat_model = generate_tests(1)
for clauses_amount_index in range(2, 20):
    for _ in range(100):
        result, got_sets, ilp_model, py3sat_model = generate_tests(clauses_amount_index)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/data/cocalc/4dea19db-f2e2-4175-bc0b-f84bc59c2366/hard-problems-formulations/minimum-exact-cover.ipynb Cell 43 line 4
      <a href='vscode-notebook-cell://xn--22-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/4dea19db-f2e2-4175-bc0b-f84bc59c2366/hard-problems-formulations/minimum-exact-cover.ipynb#X60sdnNjb2RlLXJlbW90ZQ%3D%3D?line=1'>2</a> for clauses_amount_index in range(2, 20):
      <a href='vscode-notebook-cell://xn--22-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/4dea19db-f2e2-4175-bc0b-f84bc59c2366/hard-problems-formulations/minimum-exact-cover.ipynb#X60sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a>     for _ in range(100):
----> <a href='vscode-notebook-cell://xn--22-6kce2c.xn--80apqgfe.xn--p1ai/var/data/cocalc/4dea19db-f2e2-4175-bc0b-f84bc59c2366/hard-problems-formulations/minimum-exact-cover.ipynb#X60sdnNjb2RlLXJlbW90ZQ%3D%3D?line=3'>4</a>         result, got_sets, ilp_model, py3sat_model = generate_tests(clauses_amount_index)
TypeError: cannot unpack non-iterable bool object

In [ ]:

#result_selected = [j for j in ilp_model.J if ilp_model.x[j].value > 0]
#print(result_selected)

[1, 3, 4, 5, 6]

In [ ]:

#SVG(psc.subsets2svg(got_sets, result_selected))

In [ ]: