CoCalc Logo Icon
DocsShareSupport News Sign UpSign In
Views: 22
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Embed | Download | Raw |
Kernel: Python 3 (ipykernel)

MAXIMUM SET PACKING / Максимальная упаковка множеств

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». Задача упаковки задаётся вопросом, есть ли ''k'' подмножеств в списке, которые попарно не пересекаются. Оптимизационная версия задачи, '''максимальная упаковка множеств''', задаёт вопрос о максимальном числе попарно непересекающихся множеств из «S».

https://en.wikipedia.org/wiki/Set_packing

sets = [
        [f"{2*i:02}" for i in range(8)],
        [f"{3*i:02}" for i in range(5)],
        [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)],
]


svg = psc.subsets2svg(sets)
SVG(data=svg)
Image in a Jupyter notebook

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

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)

def get_model(sets):
    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( model.x[j] for j in model.J), sense=pyo.maximize)

    @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, m.S
(['00', '02', '03', '04', '05', '06', '07', '08', '09', '10', '12', '13', '14', '17'], [{'00', '02', '04', '06', '08', '10', '12', '14'}, {'00', '03', '06', '09', '12'}, {'06'}, {'09', '14'}, {'07', '13'}, {'03', '05', '10', '17'}])
m.pprint()
2 Set Declarations x_index : Size=1, Index=None, Ordered=Insertion Key : Dimen : Domain : Size : Members None : 1 : Any : 6 : {0, 1, 2, 3, 4, 5} каждый_элемент_в_одном_множестве_index : Size=1, Index=None, Ordered=Insertion Key : Dimen : Domain : Size : Members None : 1 : Any : 14 : {'00', '02', '03', '04', '05', '06', '07', '08', '09', '10', '12', '13', '14', '17'} 1 Var Declarations x : Size=6, Index=x_index 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 1 Objective Declarations set_count : Size=1, Index=None, Active=True Key : Active : Sense : Expression None : True : maximize : x[0] + x[1] + x[2] + x[3] + x[4] + x[5] 1 Constraint Declarations каждый_элемент_в_одном_множестве : Size=14, Index=каждый_элемент_в_одном_множестве_index, Active=True Key : Lower : Body : Upper : Active 00 : -Inf : x[0] + x[1] : 1.0 : True 02 : -Inf : x[0] : 1.0 : True 03 : -Inf : x[1] + x[5] : 1.0 : True 04 : -Inf : x[0] : 1.0 : True 05 : -Inf : x[5] : 1.0 : True 06 : -Inf : x[0] + x[1] + x[2] : 1.0 : True 07 : -Inf : x[4] : 1.0 : True 08 : -Inf : x[0] : 1.0 : True 09 : -Inf : x[1] + x[3] : 1.0 : True 10 : -Inf : x[0] + x[5] : 1.0 : True 12 : -Inf : x[0] + x[1] : 1.0 : True 13 : -Inf : x[4] : 1.0 : True 14 : -Inf : x[0] + x[3] : 1.0 : True 17 : -Inf : x[5] : 1.0 : True 5 Declarations: x_index x set_count каждый_элемент_в_одном_множестве_index каждый_элемент_в_одном_множестве 
ilp_solver = pyo.SolverFactory('cbc')
ilp_solver.solve(m).write()
print_solution(m)
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 4.0 Upper bound: 4.0 Number of objectives: 1 Number of constraints: 6 Number of variables: 5 Number of binary variables: 6 Number of integer variables: 6 Number of nonzeros: 5 Sense: maximize # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: ok User time: -1.0 System time: 0.0 Wallclock time: 0.01 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.07566499710083008 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0 x[2] 1.0 x[3] 1.0 x[4] 1.0 x[5] 1.0 
selected = [j for j in m.J if m.x[j].value > 0]
selected
[2, 3, 4, 5]
SVG(psc.subsets2svg(sets, selected))
Image in a Jupyter notebook

image.png

image.png

image.png



def sat_to_set_packing(sat_instance):
    # sat_instance is a list of m clauses, each represented as a set of literals
    # e.g. [{1, -2, 3}, {2, 3, -4}, {-1, 4}]
    
    n = max(abs(l) for clause in sat_instance for l in clause)
    m = len(sat_instance)
     
    
    # create a set for each variable's occurrences in the clauses
    sets = []
    for i in range(1, n+1):
        sets.append([j for j in range(m) if i in sat_instance[j] or -i in sat_instance[j]])
    
    # construct the Set Packing instance
    set_packing_instance = [[i+1 for i in s] for s in sets]
    
    return set_packing_instance

cnf3 = CNF(from_clauses=psc.rand3cnf(2, 3))
cnf3.clauses
[(3, -1, 2), (1, -3, 2)]
satsets = sat_to_set_packing(cnf3.clauses) 
satsets
[[1, 2], [1, 2], [1, 2]]
SVG(psc.subsets2svg(satsets))
Image in a Jupyter notebook
from pysat.solvers import Solver

solver = Solver(bootstrap_with=cnf3)
res = solver.solve()
res
True
print(solver.get_model())
[-1, -2, -3] 
print(list(solver.enum_models()))
[[-1, -2, -3], [1, -2, 3], [1, 2, 3], [1, 2, -3], [-1, 2, -3], [-1, 2, 3]] 
from collections import defaultdict
def cnf32setpacking(cnf):
    '''
    Будем считать, что на входе 3SAT формула без повторов литералов в дизьюнкциях и т.п.
    '''
    n = cnf.nv
    m = len(cnf.clauses)
    
    sets = []
    items4var_1 = defaultdict(list)
    items4var_0 = defaultdict(list)
    for i in range(m):
        clause = cnf.clauses[i]
        for lit in clause:
            xn = abs(lit)
            def lit2str(lit):
                xn = abs(lit)
                return f'x_{{{xn}}}^{{{i}}}'

            itp = f'{lit2str(lit)}=1'
            items4var_1[xn].append(itp)
            itn = f'{lit2str(lit)}=0'
            items4var_0[xn].append(itn)

            if lit < 0:
                it = itp
                # Для множеств от скобок делаем наоборот!, чтобы один из вариантов оказался совместимым с множество от переменных!
            else:
                it = itn

            # remain_clause = [l for l in clause if l != lit]
            # print(lit, remain_clause)

            sets.append([f'C_{{{i}}}', it])
            # from itertools import product
            # for bits in product([0, 1], repeat=len(remain_clause)):
            #     sets.append([it, f'{lit2str(remain_clause[0])}={bits[0]}', f'{lit2str(remain_clause[1])}={bits[1]}'])        

    for j in range(1, n+1):
        sets.append([f'x_{{{j}}}'] + items4var_0[j])
        sets.append([f'x_{{{j}}}'] + items4var_1[j])

    return sets


cnf3 = CNF(from_clauses=psc.rand3cnf(1, 3))
cnf3.clauses
[(3, -1, -2)]
sets = cnf32setpacking(cnf3)
pprint(sets)
[['C_{0}', 'x_{3}^{0}=0'], ['C_{0}', 'x_{1}^{0}=1'], ['C_{0}', 'x_{2}^{0}=1'], ['x_{1}', 'x_{1}^{0}=0'], ['x_{1}', 'x_{1}^{0}=1'], ['x_{2}', 'x_{2}^{0}=0'], ['x_{2}', 'x_{2}^{0}=1'], ['x_{3}', 'x_{3}^{0}=0'], ['x_{3}', 'x_{3}^{0}=1']] 
cnf3.clauses
[(3, -1, -2)]
SVG(psc.subsets2svg(sets))
Image in a Jupyter notebook
m = get_model(sets)
m.U, m.S
(['C_{0}', 'x_{1}', 'x_{1}^{0}=0', 'x_{1}^{0}=1', 'x_{2}', 'x_{2}^{0}=0', 'x_{2}^{0}=1', 'x_{3}', 'x_{3}^{0}=0', 'x_{3}^{0}=1'], [{'C_{0}', 'x_{3}^{0}=0'}, {'C_{0}', 'x_{1}^{0}=1'}, {'C_{0}', 'x_{2}^{0}=1'}, {'x_{1}', 'x_{1}^{0}=0'}, {'x_{1}', 'x_{1}^{0}=1'}, {'x_{2}', 'x_{2}^{0}=0'}, {'x_{2}', 'x_{2}^{0}=1'}, {'x_{3}', 'x_{3}^{0}=0'}, {'x_{3}', 'x_{3}^{0}=1'}])
# ilp_solver = pyo.SolverFactory('cbc')
ilp_solver.solve(m).write()
print_solution(m)
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 4.0 Upper bound: 4.0 Number of objectives: 1 Number of constraints: 7 Number of variables: 9 Number of binary variables: 9 Number of integer variables: 9 Number of nonzeros: 9 Sense: maximize # ---------------------------------------------------------- # 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.08013749122619629 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0 x[2] 1.0 x[3] 1.0 x[5] 1.0 x[8] 1.0 
selected = [j for j in m.J if m.x[j].value > 0]
selected
[2, 3, 5, 8]
cnf3.clauses
[(3, -1, -2)]
pprint(sets)
[['C_{0}', 'x_{3}^{0}=0'], ['C_{0}', 'x_{1}^{0}=1'], ['C_{0}', 'x_{2}^{0}=1'], ['x_{1}', 'x_{1}^{0}=0'], ['x_{1}', 'x_{1}^{0}=1'], ['x_{2}', 'x_{2}^{0}=0'], ['x_{2}', 'x_{2}^{0}=1'], ['x_{3}', 'x_{3}^{0}=0'], ['x_{3}', 'x_{3}^{0}=1']] 
SVG(psc.subsets2svg(sets, selected))
Image in a Jupyter notebook


cnf3 = CNF(from_clauses=psc.rand3cnf(2, 3))
cnf3.clauses
[(2, 3, -1), (2, -1, 3)]
sets = cnf32setpacking(cnf3)
pprint(sets)
[['C_{0}', 'x_{2}^{0}=0'], ['C_{0}', 'x_{3}^{0}=0'], ['C_{0}', 'x_{1}^{0}=1'], ['C_{1}', 'x_{2}^{1}=0'], ['C_{1}', 'x_{1}^{1}=1'], ['C_{1}', 'x_{3}^{1}=0'], ['x_{1}', 'x_{1}^{0}=0', 'x_{1}^{1}=0'], ['x_{1}', 'x_{1}^{0}=1', 'x_{1}^{1}=1'], ['x_{2}', 'x_{2}^{0}=0', 'x_{2}^{1}=0'], ['x_{2}', 'x_{2}^{0}=1', 'x_{2}^{1}=1'], ['x_{3}', 'x_{3}^{0}=0', 'x_{3}^{1}=0'], ['x_{3}', 'x_{3}^{0}=1', 'x_{3}^{1}=1']] 
SVG(psc.subsets2svg(sets))
Image in a Jupyter notebook
m = get_model(sets)
m.U, m.S
(['C_{0}', 'C_{1}', 'x_{1}', 'x_{1}^{0}=0', 'x_{1}^{0}=1', 'x_{1}^{1}=0', 'x_{1}^{1}=1', 'x_{2}', 'x_{2}^{0}=0', 'x_{2}^{0}=1', 'x_{2}^{1}=0', 'x_{2}^{1}=1', 'x_{3}', 'x_{3}^{0}=0', 'x_{3}^{0}=1', 'x_{3}^{1}=0', 'x_{3}^{1}=1'], [{'C_{0}', 'x_{2}^{0}=0'}, {'C_{0}', 'x_{3}^{0}=0'}, {'C_{0}', 'x_{1}^{0}=1'}, {'C_{1}', 'x_{2}^{1}=0'}, {'C_{1}', 'x_{1}^{1}=1'}, {'C_{1}', 'x_{3}^{1}=0'}, {'x_{1}', 'x_{1}^{0}=0', 'x_{1}^{1}=0'}, {'x_{1}', 'x_{1}^{0}=1', 'x_{1}^{1}=1'}, {'x_{2}', 'x_{2}^{0}=0', 'x_{2}^{1}=0'}, {'x_{2}', 'x_{2}^{0}=1', 'x_{2}^{1}=1'}, {'x_{3}', 'x_{3}^{0}=0', 'x_{3}^{1}=0'}, {'x_{3}', 'x_{3}^{0}=1', 'x_{3}^{1}=1'}])
ilp_solver.solve(m).write()
print_solution(m)
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 5.0 Upper bound: 5.0 Number of objectives: 1 Number of constraints: 11 Number of variables: 12 Number of binary variables: 12 Number of integer variables: 12 Number of nonzeros: 12 Sense: maximize # ---------------------------------------------------------- # 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.03926873207092285 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0 x[2] 1.0 x[5] 1.0 x[6] 1.0 x[9] 1.0 x[11] 1.0 
selected = [j for j in m.J if m.x[j].value > 0]
selected
[2, 5, 6, 9, 11]
cnf3.clauses
[(2, 3, -1), (2, -1, 3)]
pprint(sets)
[['C_{0}', 'x_{2}^{0}=0'], ['C_{0}', 'x_{3}^{0}=0'], ['C_{0}', 'x_{1}^{0}=1'], ['C_{1}', 'x_{2}^{1}=0'], ['C_{1}', 'x_{1}^{1}=1'], ['C_{1}', 'x_{3}^{1}=0'], ['x_{1}', 'x_{1}^{0}=0', 'x_{1}^{1}=0'], ['x_{1}', 'x_{1}^{0}=1', 'x_{1}^{1}=1'], ['x_{2}', 'x_{2}^{0}=0', 'x_{2}^{1}=0'], ['x_{2}', 'x_{2}^{0}=1', 'x_{2}^{1}=1'], ['x_{3}', 'x_{3}^{0}=0', 'x_{3}^{1}=0'], ['x_{3}', 'x_{3}^{0}=1', 'x_{3}^{1}=1']] 
cnf3.clauses
[(2, 3, -1), (2, -1, 3)]
SVG(psc.subsets2svg(sets, selected))
Image in a Jupyter notebook