CoCalc Logo Icon
DocsShareSupport News Sign UpSign In
Views: 16
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2204
Embed | Download | Raw |
Kernel: Python 3 (ipykernel)
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 from IPython.display import display import IPython from copy import copy import os from pyomo_helpers import * import pandas from pathlib import Path from pyomo_helpers import pretty_solution 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)
d = 10 n = 5 # Number of vectors in X M = 4 vectors = np.random.randint(0, 2, size=(d, n))
def visualize(vectors): df = pd.DataFrame(vectors, columns=["x" + str(i) for i in range(vectors.shape[1])]) s = df.style.set_table_styles( [{'selector': 'td', 'props': [ ('border-left', '2px solid black !important'), ('border-right','2px solid black !important') ]}]) return s visualize(vectors)
  x0 x1 x2 x3 x4
0 1 0 1 1 0
1 0 1 1 1 0
2 1 1 1 0 0
3 0 1 0 1 0
4 0 0 1 1 0
5 1 0 1 1 1
6 1 0 1 1 0
7 1 1 1 0 0
8 1 1 0 0 1
9 0 1 0 0 0
def модель(вектора, M): m = pyo.ConcreteModel() m.множества = range(M) m.В = range(вектора.shape[1]) m.D = range(вектора.shape[0]) m.вектора = вектора m.X = pyo.Var(m.В, m.множества, domain=pyo.Binary) m.является_покрытием = pyo.Var(m.множества, domain=pyo.Binary) @m.Objective(sense=pyo.maximize) def количество_покрытий(m, в): return sum(m.является_покрытием[...]) # важно указывать слайс при суммировании! @m.Constraint(m.В) def ровно_в_одном_множестве(m, в): return sum(m.X[в, ...]) == 1 @m.Constraint(m.множества, m.D) def проверка_является_покрытием(m, мн, d): #Если хоть кто-то в мн 0 - то мн не покрытие return (m.вектора @ m.X)[d, мн] >= m.является_покрытием[мн] # @m.Constraint() # def wtf(m): # return m.X[1,0]==1 return m
m = модель(vectors, M) m.pprint()
2 Var Declarations X : Size=20, Index={0, 1, 2, 3, 4}*{0, 1, 2, 3} Key : Lower : Value : Upper : Fixed : Stale : Domain (0, 0) : 0 : None : 1 : False : True : Binary (0, 1) : 0 : None : 1 : False : True : Binary (0, 2) : 0 : None : 1 : False : True : Binary (0, 3) : 0 : None : 1 : False : True : Binary (1, 0) : 0 : None : 1 : False : True : Binary (1, 1) : 0 : None : 1 : False : True : Binary (1, 2) : 0 : None : 1 : False : True : Binary (1, 3) : 0 : None : 1 : False : True : Binary (2, 0) : 0 : None : 1 : False : True : Binary (2, 1) : 0 : None : 1 : False : True : Binary (2, 2) : 0 : None : 1 : False : True : Binary (2, 3) : 0 : None : 1 : False : True : Binary (3, 0) : 0 : None : 1 : False : True : Binary (3, 1) : 0 : None : 1 : False : True : Binary (3, 2) : 0 : None : 1 : False : True : Binary (3, 3) : 0 : None : 1 : False : True : Binary (4, 0) : 0 : None : 1 : False : True : Binary (4, 1) : 0 : None : 1 : False : True : Binary (4, 2) : 0 : None : 1 : False : True : Binary (4, 3) : 0 : None : 1 : False : True : Binary является_покрытием : Size=4, Index={0, 1, 2, 3} 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 1 Objective Declarations количество_покрытий : Size=1, Index=None, Active=True Key : Active : Sense : Expression None : True : maximize : является_покрытием[0] + является_покрытием[1] + является_покрытием[2] + является_покрытием[3] 2 Constraint Declarations проверка_является_покрытием : Size=40, Index={0, 1, 2, 3}*{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}, Active=True Key : Lower : Body : Upper : Active (0, 0) : -Inf : является_покрытием[0] - (X[0,0] + 0*X[1,0] + X[2,0] + X[3,0] + 0*X[4,0]) : 0.0 : True (0, 1) : -Inf : является_покрытием[0] - (0*X[0,0] + X[1,0] + X[2,0] + X[3,0] + 0*X[4,0]) : 0.0 : True (0, 2) : -Inf : является_покрытием[0] - (X[0,0] + X[1,0] + X[2,0] + 0*X[3,0] + 0*X[4,0]) : 0.0 : True (0, 3) : -Inf : является_покрытием[0] - (0*X[0,0] + X[1,0] + 0*X[2,0] + X[3,0] + 0*X[4,0]) : 0.0 : True (0, 4) : -Inf : является_покрытием[0] - (0*X[0,0] + 0*X[1,0] + X[2,0] + X[3,0] + 0*X[4,0]) : 0.0 : True (0, 5) : -Inf : является_покрытием[0] - (X[0,0] + 0*X[1,0] + X[2,0] + X[3,0] + X[4,0]) : 0.0 : True (0, 6) : -Inf : является_покрытием[0] - (X[0,0] + 0*X[1,0] + X[2,0] + X[3,0] + 0*X[4,0]) : 0.0 : True (0, 7) : -Inf : является_покрытием[0] - (X[0,0] + X[1,0] + X[2,0] + 0*X[3,0] + 0*X[4,0]) : 0.0 : True (0, 8) : -Inf : является_покрытием[0] - (X[0,0] + X[1,0] + 0*X[2,0] + 0*X[3,0] + X[4,0]) : 0.0 : True (0, 9) : -Inf : является_покрытием[0] - (0*X[0,0] + X[1,0] + 0*X[2,0] + 0*X[3,0] + 0*X[4,0]) : 0.0 : True (1, 0) : -Inf : является_покрытием[1] - (X[0,1] + 0*X[1,1] + X[2,1] + X[3,1] + 0*X[4,1]) : 0.0 : True (1, 1) : -Inf : является_покрытием[1] - (0*X[0,1] + X[1,1] + X[2,1] + X[3,1] + 0*X[4,1]) : 0.0 : True (1, 2) : -Inf : является_покрытием[1] - (X[0,1] + X[1,1] + X[2,1] + 0*X[3,1] + 0*X[4,1]) : 0.0 : True (1, 3) : -Inf : является_покрытием[1] - (0*X[0,1] + X[1,1] + 0*X[2,1] + X[3,1] + 0*X[4,1]) : 0.0 : True (1, 4) : -Inf : является_покрытием[1] - (0*X[0,1] + 0*X[1,1] + X[2,1] + X[3,1] + 0*X[4,1]) : 0.0 : True (1, 5) : -Inf : является_покрытием[1] - (X[0,1] + 0*X[1,1] + X[2,1] + X[3,1] + X[4,1]) : 0.0 : True (1, 6) : -Inf : является_покрытием[1] - (X[0,1] + 0*X[1,1] + X[2,1] + X[3,1] + 0*X[4,1]) : 0.0 : True (1, 7) : -Inf : является_покрытием[1] - (X[0,1] + X[1,1] + X[2,1] + 0*X[3,1] + 0*X[4,1]) : 0.0 : True (1, 8) : -Inf : является_покрытием[1] - (X[0,1] + X[1,1] + 0*X[2,1] + 0*X[3,1] + X[4,1]) : 0.0 : True (1, 9) : -Inf : является_покрытием[1] - (0*X[0,1] + X[1,1] + 0*X[2,1] + 0*X[3,1] + 0*X[4,1]) : 0.0 : True (2, 0) : -Inf : является_покрытием[2] - (X[0,2] + 0*X[1,2] + X[2,2] + X[3,2] + 0*X[4,2]) : 0.0 : True (2, 1) : -Inf : является_покрытием[2] - (0*X[0,2] + X[1,2] + X[2,2] + X[3,2] + 0*X[4,2]) : 0.0 : True (2, 2) : -Inf : является_покрытием[2] - (X[0,2] + X[1,2] + X[2,2] + 0*X[3,2] + 0*X[4,2]) : 0.0 : True (2, 3) : -Inf : является_покрытием[2] - (0*X[0,2] + X[1,2] + 0*X[2,2] + X[3,2] + 0*X[4,2]) : 0.0 : True (2, 4) : -Inf : является_покрытием[2] - (0*X[0,2] + 0*X[1,2] + X[2,2] + X[3,2] + 0*X[4,2]) : 0.0 : True (2, 5) : -Inf : является_покрытием[2] - (X[0,2] + 0*X[1,2] + X[2,2] + X[3,2] + X[4,2]) : 0.0 : True (2, 6) : -Inf : является_покрытием[2] - (X[0,2] + 0*X[1,2] + X[2,2] + X[3,2] + 0*X[4,2]) : 0.0 : True (2, 7) : -Inf : является_покрытием[2] - (X[0,2] + X[1,2] + X[2,2] + 0*X[3,2] + 0*X[4,2]) : 0.0 : True (2, 8) : -Inf : является_покрытием[2] - (X[0,2] + X[1,2] + 0*X[2,2] + 0*X[3,2] + X[4,2]) : 0.0 : True (2, 9) : -Inf : является_покрытием[2] - (0*X[0,2] + X[1,2] + 0*X[2,2] + 0*X[3,2] + 0*X[4,2]) : 0.0 : True (3, 0) : -Inf : является_покрытием[3] - (X[0,3] + 0*X[1,3] + X[2,3] + X[3,3] + 0*X[4,3]) : 0.0 : True (3, 1) : -Inf : является_покрытием[3] - (0*X[0,3] + X[1,3] + X[2,3] + X[3,3] + 0*X[4,3]) : 0.0 : True (3, 2) : -Inf : является_покрытием[3] - (X[0,3] + X[1,3] + X[2,3] + 0*X[3,3] + 0*X[4,3]) : 0.0 : True (3, 3) : -Inf : является_покрытием[3] - (0*X[0,3] + X[1,3] + 0*X[2,3] + X[3,3] + 0*X[4,3]) : 0.0 : True (3, 4) : -Inf : является_покрытием[3] - (0*X[0,3] + 0*X[1,3] + X[2,3] + X[3,3] + 0*X[4,3]) : 0.0 : True (3, 5) : -Inf : является_покрытием[3] - (X[0,3] + 0*X[1,3] + X[2,3] + X[3,3] + X[4,3]) : 0.0 : True (3, 6) : -Inf : является_покрытием[3] - (X[0,3] + 0*X[1,3] + X[2,3] + X[3,3] + 0*X[4,3]) : 0.0 : True (3, 7) : -Inf : является_покрытием[3] - (X[0,3] + X[1,3] + X[2,3] + 0*X[3,3] + 0*X[4,3]) : 0.0 : True (3, 8) : -Inf : является_покрытием[3] - (X[0,3] + X[1,3] + 0*X[2,3] + 0*X[3,3] + X[4,3]) : 0.0 : True (3, 9) : -Inf : является_покрытием[3] - (0*X[0,3] + X[1,3] + 0*X[2,3] + 0*X[3,3] + 0*X[4,3]) : 0.0 : True ровно_в_одном_множестве : Size=5, Index={0, 1, 2, 3, 4}, Active=True Key : Lower : Body : Upper : Active 0 : 1.0 : X[0,0] + X[0,1] + X[0,2] + X[0,3] : 1.0 : True 1 : 1.0 : X[1,0] + X[1,1] + X[1,2] + X[1,3] : 1.0 : True 2 : 1.0 : X[2,0] + X[2,1] + X[2,2] + X[2,3] : 1.0 : True 3 : 1.0 : X[3,0] + X[3,1] + X[3,2] + X[3,3] : 1.0 : True 4 : 1.0 : X[4,0] + X[4,1] + X[4,2] + X[4,3] : 1.0 : True 5 Declarations: X является_покрытием количество_покрытий ровно_в_одном_множестве проверка_является_покрытием
solver = pyo.SolverFactory('cbc') # solver = pyo.SolverFactory('scip') solver.solve(m).write()
# ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Name: unknown Lower bound: 1.0 Upper bound: 1.0 Number of objectives: 1 Number of constraints: 36 Number of variables: 23 Number of binary variables: 24 Number of integer variables: 24 Number of nonzeros: 4 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.03879690170288086 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0
def solve_d_vectors(vectors, M, solver_name='cbc'): m = модель(vectors, M) solver = pyo.SolverFactory(solver_name) solver.solve(m) cover = np.array([pyo.value(m.является_покрытием[i]) for i in range(M)]) solution = np.array([[int(pyo.value(m.X[i, j])) for j in range(M)] for i in range(vectors.shape[1])]) solution = solution[..., cover > 0] answer = m.количество_покрытий() return answer, solution
def visualize_solution(solution): df = pd.DataFrame(solution, columns=["s" + str(i) for i in range(solution.shape[1])], index=["x" + str(i) for i in range(solution.shape[0])]) return df.style.background_gradient()
visualize_solution(solve_d_vectors(vectors, M)[1])
  s0
x0 0
x1 1
x2 1
x3 0
x4 0

Сведение Sat к d-vector covering

Идея

Подробнее посмотрим на условие, что k-е подмножество покрывает i-ю координату вектора: jVijSjk>0\sum_{j}{V_{ij} S_{jk} > 0},

где

  • VV-матрица векторов,

  • SS - матрица множеств.

Или иначе j(Vij=Sjk=1)\lor_{j}{(V_{ij} = S_{jk} = 1)}. Это дизъюнкция.

Значит размерность d будет отвечать за число дизъюнкций в КНФ. А размерность n — число переменных в формуле.

Рассмотрим на примере следующей КНФ:

(x1x2x3)(x1x2x4)(x1 \lor \overline{x2} \lor x3) \land (\overline{x1} \lor x2 \lor \overline{x4})

closurex1x1x1\overline{x1}x2x2x2\overline{x2}x3x3x3\overline{x3}x4x4x4\overline{x4}
x1x2x3x1 \lor \overline{x2} \lor x310011000
x1x2x4\overline{x1} \lor x2 \lor \overline{x4}01100001

Тогда покрывающее множество векторов - решение КНФ

Как не взять всё?

В текущем подходе есть проблема - можно взять всё. Например одновременно переменную и её отрицание.

Решение: будем разбивать переменные в два множества. Их условно можно назвать True и False. КНФ решаема тогда и только тогда когда все переменные можно непротиворечиво разбить на два этих множества.

Осталось только форсировать непротиворечивое разбиение. Для этого добавим тавтологические дизъюнкции в КНФ xnxnxn \lor \overline{xn}

closurex1x1x1\overline{x1}x2x2x2\overline{x2}x3x3x3\overline{x3}x4x4x4\overline{x4}
x1x2x3x1 \lor \overline{x2} \lor x310011000
x1x2x4\overline{x1} \lor x2 \lor \overline{x4}01100001
x1x1x1 \lor \overline{x1}11000000
x2x2x2 \lor \overline{x2}00110000
x3x3x3 \lor \overline{x3}00001100
x4x4x4 \lor \overline{x4}00000011

Теперь чтобы получить 2 покрывающих множества одно из них обязано содержать переменную а второе её отрицание.

Если вдруг в одном множестве окажется и переменная и отрицание, то второе не сможет покрыть всё.

Например пусть разбиение следующее {x1,x1,x2,x3,x4},{x2,x3,x4}\{x1, \overline{x1}, x2, x3, x4\}, \{\overline{x2}, \overline{x3}, \overline{x4}\}

closurex1x1x1\overline{x1}x2x2x3x3x4x4sum
x1x2x3x1 \lor \overline{x2} \lor x3100102
x1x2x4\overline{x1} \lor x2 \lor \overline{x4}011002
x1x1x1 \lor \overline{x1}110002
x2x2x2 \lor \overline{x2}001001
x3x3x3 \lor \overline{x3}000101
x4x4x4 \lor \overline{x4}000011
closurex2\overline{x2}x3\overline{x3}x4\overline{x4}sum
x1x2x3x1 \lor \overline{x2} \lor x31001
x1x2x4\overline{x1} \lor x2 \lor \overline{x4}0011
x1x1x1 \lor \overline{x1}0000
x2x2x2 \lor \overline{x2}1001
x3x3x3 \lor \overline{x3}0101
x4x4x4 \lor \overline{x4}0011

Как видно, второе множество не является покрывающим

Как найти второе решение

У текущего подхода есть ещё одна проблема. В общем случае неверно, что если (x1,x2,)(x1, x2, \ldots) решение SAT то (x1,x2,)(\overline{x1}, \overline{x2}, \ldots) тоже решение. Но надо как то сделать так, чтобы отрицание было покрывающем множеством. Для этого введем новую переменную (назовем её ff). И объявим её ложной.

То есть сделаем ещё одно эквивалентное преобразование: (x1x2x3)(x1x2x3False)(x1x2x3f)(x1 \lor \overline{x2} \lor x3) \Rightarrow (x1 \lor \overline{x2} \lor x3 \lor False) \Rightarrow (x1 \lor \overline{x2} \lor x3 \lor f) при f=Falsef = False Проделаем так с каждой оригинальной скобкой в КНФ:

closurex1x1x1\overline{x1}x2x2x2\overline{x2}x3x3x3\overline{x3}x4x4x4\overline{x4}f
x1x2x3x1 \lor \overline{x2} \lor x3100110001
x1x2x4\overline{x1} \lor x2 \lor \overline{x4}011000011
x1x1x1 \lor \overline{x1}110000000
x2x2x2 \lor \overline{x2}001100000
x3x3x3 \lor \overline{x3}000011000
x4x4x4 \lor \overline{x4}000000110

Итого:

  • Если нашлось разбиение на 2 множества, то ни одно из них не содержит одновременно переменную и её отрицание

  • Множество не содержащее ff является решением КНФ.

cnf_vars = 5 cnf_clauses = 20 cnf3 = CNF(from_clauses=psc.rand3cnf(cnf_clauses, cnf_vars)) clauses = cnf3.clauses
def solve_sat(cnf_vars, cnf_clauses, clauses): vectors = np.zeros((cnf_clauses + cnf_vars, 2*cnf_vars + 1)) for i in range(cnf_clauses): for j in clauses[i]: vectors[i][j + cnf_vars] = 1 vectors[i][cnf_vars] = 1 # Contradict for i in range(1, cnf_vars+1): vectors[cnf_clauses + i - 1][cnf_vars - i] = 1 vectors[cnf_clauses + i - 1][cnf_vars + i] = 1 sol = solve_d_vectors(vectors, 2) return sol[0] == 2
solve_sat(cnf_vars, cnf_clauses, clauses)
True
cnf_vars = 5 cnf_clauses = 1 cnf3 = CNF(from_clauses=psc.rand3cnf(cnf_clauses, cnf_vars)) clauses = cnf3.clauses
clauses
[[-3, 2, 1]]
solve_sat(cnf_vars, cnf_clauses, clauses)
True
def test_sat_iter(N, cnf_vars, cnf_clauses): mistakes = 0 for _ in range(N): cnf3 = CNF(from_clauses=psc.rand3cnf(cnf_clauses, cnf_vars)) clauses = cnf3.clauses sol_this = solve_sat(cnf_vars, cnf_clauses, clauses) py3sat_test_solver = Solver(bootstrap_with=cnf3) sol_other = py3sat_test_solver.solve() if sol_this != sol_other: mistakes += 1 return mistakes
def stress_test(): for i in range(100): if test_sat_iter(20, 4, 20) > 5: print("Failed") return if i % 10 == 0: print("Passed ", i)
stress_test()
Passed 0 Passed 10 Passed 20 Passed 30 Passed 40 Passed 50 Passed 60 Passed 70 Passed 80 Passed 90