test_lra_theory.py 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440
  1. from sympy.core.numbers import Rational, I, oo
  2. from sympy.core.relational import Eq
  3. from sympy.core.symbol import symbols
  4. from sympy.core.singleton import S
  5. from sympy.matrices.dense import Matrix
  6. from sympy.matrices.dense import randMatrix
  7. from sympy.assumptions.ask import Q
  8. from sympy.logic.boolalg import And
  9. from sympy.abc import x, y, z
  10. from sympy.assumptions.cnf import CNF, EncodedCNF
  11. from sympy.functions.elementary.trigonometric import cos
  12. from sympy.external import import_module
  13. from sympy.logic.algorithms.lra_theory import LRASolver, UnhandledInput, LRARational, HANDLE_NEGATION
  14. from sympy.core.random import random, choice, randint
  15. from sympy.core.sympify import sympify
  16. from sympy.ntheory.generate import randprime
  17. from sympy.core.relational import StrictLessThan, StrictGreaterThan
  18. import itertools
  19. from sympy.testing.pytest import raises, XFAIL, skip
  20. def make_random_problem(num_variables=2, num_constraints=2, sparsity=.1, rational=True,
  21. disable_strict = False, disable_nonstrict=False, disable_equality=False):
  22. def rand(sparsity=sparsity):
  23. if random() < sparsity:
  24. return sympify(0)
  25. if rational:
  26. int1, int2 = [randprime(0, 50) for _ in range(2)]
  27. return Rational(int1, int2) * choice([-1, 1])
  28. else:
  29. return randint(1, 10) * choice([-1, 1])
  30. variables = symbols('x1:%s' % (num_variables + 1))
  31. constraints = []
  32. for _ in range(num_constraints):
  33. lhs, rhs = sum(rand() * x for x in variables), rand(sparsity=0) # sparsity=0 bc of bug with smtlib_code
  34. options = []
  35. if not disable_equality:
  36. options += [Eq(lhs, rhs)]
  37. if not disable_nonstrict:
  38. options += [lhs <= rhs, lhs >= rhs]
  39. if not disable_strict:
  40. options += [lhs < rhs, lhs > rhs]
  41. constraints.append(choice(options))
  42. return constraints
  43. def check_if_satisfiable_with_z3(constraints):
  44. from sympy.external.importtools import import_module
  45. from sympy.printing.smtlib import smtlib_code
  46. from sympy.logic.boolalg import And
  47. boolean_formula = And(*constraints)
  48. z3 = import_module("z3")
  49. if z3:
  50. smtlib_string = smtlib_code(boolean_formula)
  51. s = z3.Solver()
  52. s.from_string(smtlib_string)
  53. res = str(s.check())
  54. if res == 'sat':
  55. return True
  56. elif res == 'unsat':
  57. return False
  58. else:
  59. raise ValueError(f"z3 was not able to check the satisfiability of {boolean_formula}")
  60. def find_rational_assignment(constr, assignment, iter=20):
  61. eps = sympify(1)
  62. for _ in range(iter):
  63. assign = {key: val[0] + val[1]*eps for key, val in assignment.items()}
  64. try:
  65. for cons in constr:
  66. assert cons.subs(assign) == True
  67. return assign
  68. except AssertionError:
  69. eps = eps/2
  70. return None
  71. def boolean_formula_to_encoded_cnf(bf):
  72. cnf = CNF.from_prop(bf)
  73. enc = EncodedCNF()
  74. enc.from_cnf(cnf)
  75. return enc
  76. def test_from_encoded_cnf():
  77. s1, s2 = symbols("s1 s2")
  78. # Test preprocessing
  79. # Example is from section 3 of paper.
  80. phi = (x >= 0) & ((x + y <= 2) | (x + 2 * y - z >= 6)) & (Eq(x + y, 2) | (x + 2 * y - z > 4))
  81. enc = boolean_formula_to_encoded_cnf(phi)
  82. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  83. assert lra.A.shape == (2, 5)
  84. assert str(lra.slack) == '[_s1, _s2]'
  85. assert str(lra.nonslack) == '[x, y, z]'
  86. assert lra.A == Matrix([[ 1, 1, 0, -1, 0],
  87. [-1, -2, 1, 0, -1]])
  88. assert {(str(b.var), b.bound, b.upper, b.equality, b.strict) for b in lra.enc_to_boundary.values()} == {('_s1', 2, None, True, False),
  89. ('_s1', 2, True, False, False),
  90. ('_s2', -4, True, False, True),
  91. ('_s2', -6, True, False, False),
  92. ('x', 0, False, False, False)}
  93. def test_problem():
  94. from sympy.logic.algorithms.lra_theory import LRASolver
  95. from sympy.assumptions.cnf import CNF, EncodedCNF
  96. cons = [-2 * x - 2 * y >= 7, -9 * y >= 7, -6 * y >= 5]
  97. cnf = CNF().from_prop(And(*cons))
  98. enc = EncodedCNF()
  99. enc.from_cnf(cnf)
  100. lra, _ = LRASolver.from_encoded_cnf(enc)
  101. lra.assert_lit(1)
  102. lra.assert_lit(2)
  103. lra.assert_lit(3)
  104. is_sat, assignment = lra.check()
  105. assert is_sat is True
  106. def test_random_problems():
  107. z3 = import_module("z3")
  108. if z3 is None:
  109. skip("z3 is not installed")
  110. special_cases = []; x1, x2, x3 = symbols("x1 x2 x3")
  111. special_cases.append([x1 - 3 * x2 <= -5, 6 * x1 + 4 * x2 <= 0, -7 * x1 + 3 * x2 <= 3])
  112. special_cases.append([-3 * x1 >= 3, Eq(4 * x1, -1)])
  113. special_cases.append([-4 * x1 < 4, 6 * x1 <= -6])
  114. special_cases.append([-3 * x2 >= 7, 6 * x1 <= -5, -3 * x2 <= -4])
  115. special_cases.append([x + y >= 2, x + y <= 1])
  116. special_cases.append([x >= 0, x + y <= 2, x + 2 * y - z >= 6]) # from paper example
  117. special_cases.append([-2 * x1 - 2 * x2 >= 7, -9 * x1 >= 7, -6 * x1 >= 5])
  118. special_cases.append([2 * x1 > -3, -9 * x1 < -6, 9 * x1 <= 6])
  119. special_cases.append([-2*x1 < -4, 9*x1 > -9])
  120. special_cases.append([-6*x1 >= -1, -8*x1 + x2 >= 5, -8*x1 + 7*x2 < 4, x1 > 7])
  121. special_cases.append([Eq(x1, 2), Eq(5*x1, -2), Eq(-7*x2, -6), Eq(9*x1 + 10*x2, 9)])
  122. special_cases.append([Eq(3*x1, 6), Eq(x1 - 8*x2, -9), Eq(-7*x1 + 5*x2, 3), Eq(3*x2, 7)])
  123. special_cases.append([-4*x1 < 4, 6*x1 <= -6])
  124. special_cases.append([-3*x1 + 8*x2 >= -8, -10*x2 > 9, 8*x1 - 4*x2 < 8, 10*x1 - 9*x2 >= -9])
  125. special_cases.append([x1 + 5*x2 >= -6, 9*x1 - 3*x2 >= -9, 6*x1 + 6*x2 < -10, -3*x1 + 3*x2 < -7])
  126. special_cases.append([-9*x1 < 7, -5*x1 - 7*x2 < -1, 3*x1 + 7*x2 > 1, -6*x1 - 6*x2 > 9])
  127. special_cases.append([9*x1 - 6*x2 >= -7, 9*x1 + 4*x2 < -8, -7*x2 <= 1, 10*x2 <= -7])
  128. feasible_count = 0
  129. for i in range(50):
  130. if i % 8 == 0:
  131. constraints = make_random_problem(num_variables=1, num_constraints=2, rational=False)
  132. elif i % 8 == 1:
  133. constraints = make_random_problem(num_variables=2, num_constraints=4, rational=False, disable_equality=True,
  134. disable_nonstrict=True)
  135. elif i % 8 == 2:
  136. constraints = make_random_problem(num_variables=2, num_constraints=4, rational=False, disable_strict=True)
  137. elif i % 8 == 3:
  138. constraints = make_random_problem(num_variables=3, num_constraints=12, rational=False)
  139. else:
  140. constraints = make_random_problem(num_variables=3, num_constraints=6, rational=False)
  141. if i < len(special_cases):
  142. constraints = special_cases[i]
  143. if False in constraints or True in constraints:
  144. continue
  145. phi = And(*constraints)
  146. if phi == False:
  147. continue
  148. cnf = CNF.from_prop(phi); enc = EncodedCNF()
  149. enc.from_cnf(cnf)
  150. assert all(0 not in clause for clause in enc.data)
  151. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  152. s_subs = lra.s_subs
  153. lra.run_checks = True
  154. s_subs_rev = {value: key for key, value in s_subs.items()}
  155. lits = {lit for clause in enc.data for lit in clause}
  156. bounds = [(lra.enc_to_boundary[l], l) for l in lits if l in lra.enc_to_boundary]
  157. bounds = sorted(bounds, key=lambda x: (str(x[0].var), x[0].bound, str(x[0].upper))) # to remove nondeterminism
  158. for b, l in bounds:
  159. if lra.result and lra.result[0] == False:
  160. break
  161. lra.assert_lit(l)
  162. feasible = lra.check()
  163. if feasible[0] == True:
  164. feasible_count += 1
  165. assert check_if_satisfiable_with_z3(constraints) is True
  166. cons_funcs = [cons.func for cons in constraints]
  167. assignment = feasible[1]
  168. assignment = {key.var : value for key, value in assignment.items()}
  169. if not (StrictLessThan in cons_funcs or StrictGreaterThan in cons_funcs):
  170. assignment = {key: value[0] for key, value in assignment.items()}
  171. for cons in constraints:
  172. assert cons.subs(assignment) == True
  173. else:
  174. rat_assignment = find_rational_assignment(constraints, assignment)
  175. assert rat_assignment is not None
  176. else:
  177. assert check_if_satisfiable_with_z3(constraints) is False
  178. conflict = feasible[1]
  179. assert len(conflict) >= 2
  180. conflict = {lra.enc_to_boundary[-l].get_inequality() for l in conflict}
  181. conflict = {clause.subs(s_subs_rev) for clause in conflict}
  182. assert check_if_satisfiable_with_z3(conflict) is False
  183. # check that conflict clause is probably minimal
  184. for subset in itertools.combinations(conflict, len(conflict)-1):
  185. assert check_if_satisfiable_with_z3(subset) is True
  186. @XFAIL
  187. def test_pos_neg_zero():
  188. bf = Q.positive(x) & Q.negative(x) & Q.zero(y)
  189. enc = boolean_formula_to_encoded_cnf(bf)
  190. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  191. for lit in enc.encoding.values():
  192. if lra.assert_lit(lit) is not None:
  193. break
  194. assert len(lra.enc_to_boundary) == 3
  195. assert lra.check()[0] == False
  196. bf = Q.positive(x) & Q.lt(x, -1)
  197. enc = boolean_formula_to_encoded_cnf(bf)
  198. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  199. for lit in enc.encoding.values():
  200. if lra.assert_lit(lit) is not None:
  201. break
  202. assert len(lra.enc_to_boundary) == 2
  203. assert lra.check()[0] == False
  204. bf = Q.positive(x) & Q.zero(x)
  205. enc = boolean_formula_to_encoded_cnf(bf)
  206. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  207. for lit in enc.encoding.values():
  208. if lra.assert_lit(lit) is not None:
  209. break
  210. assert len(lra.enc_to_boundary) == 2
  211. assert lra.check()[0] == False
  212. bf = Q.positive(x) & Q.zero(y)
  213. enc = boolean_formula_to_encoded_cnf(bf)
  214. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  215. for lit in enc.encoding.values():
  216. if lra.assert_lit(lit) is not None:
  217. break
  218. assert len(lra.enc_to_boundary) == 2
  219. assert lra.check()[0] == True
  220. @XFAIL
  221. def test_pos_neg_infinite():
  222. bf = Q.positive_infinite(x) & Q.lt(x, 10000000) & Q.positive_infinite(y)
  223. enc = boolean_formula_to_encoded_cnf(bf)
  224. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  225. for lit in enc.encoding.values():
  226. if lra.assert_lit(lit) is not None:
  227. break
  228. assert len(lra.enc_to_boundary) == 3
  229. assert lra.check()[0] == False
  230. bf = Q.positive_infinite(x) & Q.gt(x, 10000000) & Q.positive_infinite(y)
  231. enc = boolean_formula_to_encoded_cnf(bf)
  232. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  233. for lit in enc.encoding.values():
  234. if lra.assert_lit(lit) is not None:
  235. break
  236. assert len(lra.enc_to_boundary) == 3
  237. assert lra.check()[0] == True
  238. bf = Q.positive_infinite(x) & Q.negative_infinite(x)
  239. enc = boolean_formula_to_encoded_cnf(bf)
  240. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  241. for lit in enc.encoding.values():
  242. if lra.assert_lit(lit) is not None:
  243. break
  244. assert len(lra.enc_to_boundary) == 2
  245. assert lra.check()[0] == False
  246. def test_binrel_evaluation():
  247. bf = Q.gt(3, 2)
  248. enc = boolean_formula_to_encoded_cnf(bf)
  249. lra, conflicts = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  250. assert len(lra.enc_to_boundary) == 0
  251. assert conflicts == [[1]]
  252. bf = Q.lt(3, 2)
  253. enc = boolean_formula_to_encoded_cnf(bf)
  254. lra, conflicts = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  255. assert len(lra.enc_to_boundary) == 0
  256. assert conflicts == [[-1]]
  257. def test_negation():
  258. assert HANDLE_NEGATION is True
  259. bf = Q.gt(x, 1) & ~Q.gt(x, 0)
  260. enc = boolean_formula_to_encoded_cnf(bf)
  261. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  262. for clause in enc.data:
  263. for lit in clause:
  264. lra.assert_lit(lit)
  265. assert len(lra.enc_to_boundary) == 2
  266. assert lra.check()[0] == False
  267. assert sorted(lra.check()[1]) in [[-1, 2], [-2, 1]]
  268. bf = ~Q.gt(x, 1) & ~Q.lt(x, 0)
  269. enc = boolean_formula_to_encoded_cnf(bf)
  270. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  271. for clause in enc.data:
  272. for lit in clause:
  273. lra.assert_lit(lit)
  274. assert len(lra.enc_to_boundary) == 2
  275. assert lra.check()[0] == True
  276. bf = ~Q.gt(x, 0) & ~Q.lt(x, 1)
  277. enc = boolean_formula_to_encoded_cnf(bf)
  278. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  279. for clause in enc.data:
  280. for lit in clause:
  281. lra.assert_lit(lit)
  282. assert len(lra.enc_to_boundary) == 2
  283. assert lra.check()[0] == False
  284. bf = ~Q.gt(x, 0) & ~Q.le(x, 0)
  285. enc = boolean_formula_to_encoded_cnf(bf)
  286. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  287. for clause in enc.data:
  288. for lit in clause:
  289. lra.assert_lit(lit)
  290. assert len(lra.enc_to_boundary) == 2
  291. assert lra.check()[0] == False
  292. bf = ~Q.le(x+y, 2) & ~Q.ge(x-y, 2) & ~Q.ge(y, 0)
  293. enc = boolean_formula_to_encoded_cnf(bf)
  294. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  295. for clause in enc.data:
  296. for lit in clause:
  297. lra.assert_lit(lit)
  298. assert len(lra.enc_to_boundary) == 3
  299. assert lra.check()[0] == False
  300. assert len(lra.check()[1]) == 3
  301. assert all(i > 0 for i in lra.check()[1])
  302. def test_unhandled_input():
  303. nan = S.NaN
  304. bf = Q.gt(3, nan) & Q.gt(x, nan)
  305. enc = boolean_formula_to_encoded_cnf(bf)
  306. raises(ValueError, lambda: LRASolver.from_encoded_cnf(enc, testing_mode=True))
  307. bf = Q.gt(3, I) & Q.gt(x, I)
  308. enc = boolean_formula_to_encoded_cnf(bf)
  309. raises(UnhandledInput, lambda: LRASolver.from_encoded_cnf(enc, testing_mode=True))
  310. bf = Q.gt(3, float("inf")) & Q.gt(x, float("inf"))
  311. enc = boolean_formula_to_encoded_cnf(bf)
  312. raises(UnhandledInput, lambda: LRASolver.from_encoded_cnf(enc, testing_mode=True))
  313. bf = Q.gt(3, oo) & Q.gt(x, oo)
  314. enc = boolean_formula_to_encoded_cnf(bf)
  315. raises(UnhandledInput, lambda: LRASolver.from_encoded_cnf(enc, testing_mode=True))
  316. # test non-linearity
  317. bf = Q.gt(x**2 + x, 2)
  318. enc = boolean_formula_to_encoded_cnf(bf)
  319. raises(UnhandledInput, lambda: LRASolver.from_encoded_cnf(enc, testing_mode=True))
  320. bf = Q.gt(cos(x) + x, 2)
  321. enc = boolean_formula_to_encoded_cnf(bf)
  322. raises(UnhandledInput, lambda: LRASolver.from_encoded_cnf(enc, testing_mode=True))
  323. @XFAIL
  324. def test_infinite_strict_inequalities():
  325. # Extensive testing of the interaction between strict inequalities
  326. # and constraints containing infinity is needed because
  327. # the paper's rule for strict inequalities don't work when
  328. # infinite numbers are allowed. Using the paper's rules you
  329. # can end up with situations where oo + delta > oo is considered
  330. # True when oo + delta should be equal to oo.
  331. # See https://math.stackexchange.com/questions/4757069/can-this-method-of-converting-strict-inequalities-to-equisatisfiable-nonstrict-i
  332. bf = (-x - y >= -float("inf")) & (x > 0) & (y >= float("inf"))
  333. enc = boolean_formula_to_encoded_cnf(bf)
  334. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  335. for lit in sorted(enc.encoding.values()):
  336. if lra.assert_lit(lit) is not None:
  337. break
  338. assert len(lra.enc_to_boundary) == 3
  339. assert lra.check()[0] == True
  340. def test_pivot():
  341. for _ in range(10):
  342. m = randMatrix(5)
  343. rref = m.rref()
  344. for _ in range(5):
  345. i, j = randint(0, 4), randint(0, 4)
  346. if m[i, j] != 0:
  347. assert LRASolver._pivot(m, i, j).rref() == rref
  348. def test_reset_bounds():
  349. bf = Q.ge(x, 1) & Q.lt(x, 1)
  350. enc = boolean_formula_to_encoded_cnf(bf)
  351. lra, _ = LRASolver.from_encoded_cnf(enc, testing_mode=True)
  352. for clause in enc.data:
  353. for lit in clause:
  354. lra.assert_lit(lit)
  355. assert len(lra.enc_to_boundary) == 2
  356. assert lra.check()[0] == False
  357. lra.reset_bounds()
  358. assert lra.check()[0] == True
  359. for var in lra.all_var:
  360. assert var.upper == LRARational(float("inf"), 0)
  361. assert var.upper_from_eq == False
  362. assert var.upper_from_neg == False
  363. assert var.lower == LRARational(-float("inf"), 0)
  364. assert var.lower_from_eq == False
  365. assert var.lower_from_neg == False
  366. assert var.assign == LRARational(0, 0)
  367. assert var.var is not None
  368. assert var.col_idx is not None
  369. def test_empty_cnf():
  370. cnf = CNF()
  371. enc = EncodedCNF()
  372. enc.from_cnf(cnf)
  373. lra, conflict = LRASolver.from_encoded_cnf(enc)
  374. assert len(conflict) == 0
  375. assert lra.check() == (True, {})