MAXIMUM SET SPLITTING / Максимальное расщепление множеств

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
from collections import defaultdict
from itertools import combinations
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)

def get_random_subsets(N, sizes):
    res = []
    for size_, repeat_ in sizes:
        for _ in range(repeat_):
#             ok = False
#             while not ok:
#                 set_ = set(f"{i:02}" for i in np.random.choice(np.arange(1, N), size_))
#                 if len(set_) == size_:
#                     ok = True
            set_ = set()
            while len(set_) != size_:
                set_ = set(f"{i:02}" for i in np.random.choice(np.arange(1, N), size_))
            res.append(set_)
    return res

N = 20
sets = get_random_subsets(N, [(2, 10), (3, 2), (8, 1)])
#  [(num_of_elements, num_of_subset)]

svg = psc.subsets2svg(sets)
SVG(data=svg)
# Генерация происходит сверху вниз, вы можете наглядно проверить, что программа сгенирировала то, что надо

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):
    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.I = range(model.m)
    model.U = sorted(set([str(item) for sublist in sets_ for item in sublist])) # universe множество всех элементов (отсортированных)
    model.n = len(model.U) # длина подмножества по горизонтали
    model.J = range(model.n)

    # как мы делаем разрез по элементам
    model.x = pyo.Var(model.J, domain=pyo.Binary)

    # пересекает ли подмножество этот разрез
    model.y = pyo.Var(model.I, domain=pyo.Binary)
    model.количество_разрезанных_подмножеств = pyo.Objective(expr = sum( model.y[i] for i in model.I), sense=pyo.maximize)

    # Для корректного разреза мы должны как выбирать, так и не выбирать элемент:
    # Выбираем элемент
    @model.Constraint(model.I)
    def надо_что_то_выбрать_в_множестве_чтобы_считать_его_разрезанным(m, i):
        return sum(m.x[j] for j in m.J if m.U[j] in m.S[i]) >= m.y[i]

    # Не выбираем элемент
    @model.Constraint(model.I)
    def надо_что_то_выбрать_в_множестве_чтобы_считать_его_разрезанным2(m, i):
        return len(m.S[i]) - sum(m.x[j] for j in m.J if m.U[j] in m.S[i]) >= m.y[i]

    return model


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

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

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

selected_sets = [i for i in m.I if m.y[i].value > 0]
selected_sets

SVG(psc.subsets2svg(sets, selected=selected_sets, selected_items=selected_items))

sets

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

from pysat.solvers import Solver

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

print(solver.get_model())

print(list(solver.enum_models()))

def lit2str(lit, i=None):
    '''
    литерал в элемент-строку
    '''
    xn = abs(lit)
    mod = ''
    if i is not None:
        mod = f'^{{{i}}}'
    res = f'x{mod}_{{{xn}}}'
    if lit < 0:
        res = f'\overline{{{res}}}'
#         res = f'¬{{{res}}}'
    return res

# def clause2tex(clause):
#     '''
#     скобка в латех-строку
#     '''
#     res = '∧'.join([lit2str(l) for l in clause])
#     return res
0

def cnf32setsplitting(cnf):
    '''
    Будем считать, что на входе 3SAT формула. 
    Можно улучшить и т.п.
    '''
    n = cnf.nv
    m = len(cnf.clauses)
    
    sets = {}

    # для всех переменных добавляем пары литералов
    # разрезать это множество — выбрать значение.
    for j in range(1, n+1):
        sets[f'x_{{{j}}}'] = [lit2str(j), lit2str(-j)]         

    for i in range(m):
        clause = cnf.clauses[i]
        # Для скобок делает «отрезание от нуля» — специального символа «О»
        clause_set = ['O']
        for lit in clause:
            sets[f'C_{{{i}}}({lit2str(lit)})'] = [lit2str(lit, i), lit2str(-lit)]         
            clause_set.append(lit2str(lit, i))
        sets[f'C_{{{i}}}'] = clause_set
    
    return sets

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

sets = cnf32setsplitting(cnf3)
pprint(sets)

SVG(psc.subsets2svg(sets))

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

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

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

splitted_sets = [i for i in m.I if m.y[i].value > 0]
splitted_sets

cnf3.clauses

pprint(sets)

SVG(psc.subsets2svg(sets, selected_items=selected_items, selected=splitted_sets))

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

sets = cnf32setsplitting(cnf3)
pprint(sets)

SVG(psc.subsets2svg(sets))

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

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

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

splitted_sets = [i for i in m.I if m.y[i].value > 0]
splitted_sets

cnf3.clauses

pprint(sets)

SVG(psc.subsets2svg(sets, selected_items=selected_items, selected=splitted_sets))

import re
def max_set_splitting_sol_2_3satsol(m, sets):
    '''
    Конвертация решения этой задачи через ЦЛП, в решение для 3SAT
    '''
    # словарь номер_переменной_в_3SAT → значение (0/1)
    satsol = {}
    
    selected_items = [m.U[j] for j in m.J if m.x[j].value > 0]
    splitted_sets = [i for i in m.I if m.y[i].value > 0]
    O_value = m.y[0].value
    need_negation = O_value
    re_ = re.compile(r"x\_\{(?P<varn>\d+)\}")
    for j in m.J:
        rm = re_.match(m.U[j])
        if rm:
            val = int(m.x[j].value)
            varn = int(rm.group('varn'))
            if need_negation:
                val = not val
            assert varn not in satsol or val == satsol[varn]
            satsol[varn] = val
    return satsol     

cnf.clauses

max_set_splitting_sol_2_3satsol(m, sets)

def solve_cnf_by_set_splitting(cnf):
    sets = cnf32setsplitting(cnf)
    m = get_model(sets)
    ilp_solver.solve(m)   
    sol = max_set_splitting_sol_2_3satsol(m, sets)
    return sol  

cnf = CNF(from_clauses=psc.rand3cnf(7, 4))
cnf.clauses, solve_cnf_by_set_splitting(cnf)

from pyomo.opt import SolverStatus
from pyomo.opt import TerminationCondition

def create_tests(num_of_clauses):
    cnf = CNF(from_clauses=psc.rand3cnf(num_of_clauses, 3))
    test_sets = cnf32setsplitting(cnf)
    test_model = get_model(test_sets)
    test_result_struct = ilp_solver.solve(test_model)
    test_result = (test_result_struct.solver.status == SolverStatus.ok and test_result_struct.solver.termination_condition == TerminationCondition.optimal)
    test_no_result = (test_result_struct.solver.termination_condition == TerminationCondition.infeasible)
    
    py3sat_test_solver = Solver(bootstrap_with=cnf)
    py3sat_test_result = py3sat_test_solver.solve()
    
    if ((py3sat_test_result == True and test_no_result)\
        or\
        (py3sat_test_result == False and test_result)):
        print("SOLUTION FAILED num amount: ", num_of_clauses, "Pyomo: ", test_result, ", test no result: ", test_no_result, ", 3SAT: ", py3sat_test_result)
        return False, test_sets, test_model, py3sat_test_solver.get_model()
    else:
        print("SOLUTION FOUND num amount: ", num_of_clauses, ", answer: ", py3sat_test_result)
        return True, test_sets, test_model, py3sat_test_solver.get_model()

create_tests(5)

for i in range(1, 50):
    result_of_tests = create_tests(i)
    done = result_of_tests[0]
    bad_case_sets = result_of_tests[1]
    model = result_of_tests[2]
    py3sat_model = result_of_tests[3]
    if done != True:
        break

print(bad_case_sets, '\n')

bad_case_selected = [model.U[j] for j in model.J if model.x[j].value > 0]
bad_case_splitted = [i for i in model.I if model.y[i].value > 0]

print("-----\n")
print(bad_case_selected)
print(bad_case_splitted)

SVG(psc.subsets2svg(got_sets, selected_items=result_selected, selected=ssplitted_sets))