Branch data Line data Source code
1 : : /******************************************************************************
2 : : * This file is part of the cvc5 project.
3 : : *
4 : : * Copyright (c) 2009-2026 by the authors listed in the file AUTHORS
5 : : * in the top-level source directory and their institutional affiliations.
6 : : * All rights reserved. See the file COPYING in the top-level source
7 : : * directory for licensing information.
8 : : * ****************************************************************************
9 : : *
10 : : * Implements the CDCAC approach as described in
11 : : * https://arxiv.org/pdf/2003.05633.pdf.
12 : : */
13 : :
14 : : #include "theory/arith/nl/coverings/cdcac.h"
15 : :
16 : : #ifdef CVC5_POLY_IMP
17 : :
18 : : #include "options/arith_options.h"
19 : : #include "theory/arith/nl/coverings/lazard_evaluation.h"
20 : : #include "theory/arith/nl/coverings/projections.h"
21 : : #include "theory/arith/nl/coverings/variable_ordering.h"
22 : : #include "theory/arith/nl/nl_model.h"
23 : : #include "theory/rewriter.h"
24 : : #include "util/resource_manager.h"
25 : :
26 : : using namespace cvc5::internal::kind;
27 : :
28 : : namespace std {
29 : : /** Generic streaming operator for std::vector. */
30 : : template <typename T>
31 : 0 : std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
32 : : {
33 : 0 : cvc5::internal::container_to_stream(os, v);
34 : 0 : return os;
35 : : }
36 : : } // namespace std
37 : :
38 : : namespace cvc5::internal {
39 : : namespace theory {
40 : : namespace arith {
41 : : namespace nl {
42 : : namespace coverings {
43 : :
44 : 34056 : CDCAC::CDCAC(Env& env, const std::vector<poly::Variable>& ordering)
45 : : : EnvObj(env),
46 : 34056 : d_assignment(nodeManager()->getPolyContext()),
47 : 34056 : d_constraints(nodeManager()->getPolyContext()),
48 : 34056 : d_variableOrdering(ordering),
49 : 68112 : d_varOrder(nodeManager()->getPolyContext())
50 : : {
51 [ + + ]: 34056 : if (d_env.isTheoryProofProducing())
52 : : {
53 : 7532 : d_proof.reset(new CoveringsProofGenerator(env, userContext()));
54 : : }
55 : 34056 : }
56 : :
57 : 274 : void CDCAC::reset()
58 : : {
59 : 274 : d_constraints.reset();
60 : 274 : d_assignment.clear();
61 : 274 : d_nextIntervalId = 1;
62 : 274 : }
63 : :
64 : 0 : poly::detail::variable_printer CDCAC::get_stream_variable(
65 : : const poly::Variable& v)
66 : : {
67 : 0 : const poly::Context& ctx = nodeManager()->getPolyContext();
68 : 0 : return stream_variable(ctx, v);
69 : : }
70 : :
71 : 0 : std::vector<poly::detail::variable_printer> CDCAC::get_stream_variables(
72 : : const std::vector<poly::Variable>& vars)
73 : : {
74 : 0 : std::vector<poly::detail::variable_printer> result;
75 : 0 : result.reserve(vars.size());
76 : :
77 : 0 : const poly::Context& ctx = nodeManager()->getPolyContext();
78 : :
79 [ - - ]: 0 : for (const poly::Variable& v : vars)
80 : : {
81 : 0 : result.emplace_back(stream_variable(ctx, v));
82 : : }
83 : :
84 : 0 : return result;
85 : 0 : }
86 : :
87 : 279 : void CDCAC::computeVariableOrdering()
88 : : {
89 : : // Actually compute the variable ordering
90 : 558 : d_variableOrdering = d_varOrder(d_constraints.getConstraints(),
91 : 279 : VariableOrderingStrategy::BROWN);
92 [ + - ]: 558 : Trace("cdcac") << "Variable ordering is now "
93 [ - + ][ - - ]: 279 : << get_stream_variables(d_variableOrdering) << std::endl;
94 : :
95 : : // Write variable ordering back to libpoly.
96 : : lp_variable_order_t* vo =
97 : 279 : nodeManager()->getPolyContext().get_variable_order();
98 : 279 : lp_variable_order_clear(vo);
99 [ + + ]: 1313 : for (const auto& v : d_variableOrdering)
100 : : {
101 : 1034 : lp_variable_order_push(vo, v.get_internal());
102 : : }
103 : 279 : }
104 : :
105 : 272 : void CDCAC::retrieveInitialAssignment(NlModel& model, const Node& ran_variable)
106 : : {
107 [ + - ]: 272 : if (options().arith.nlCovLinearModel == options::nlCovLinearModelMode::NONE)
108 : 272 : return;
109 : 0 : d_initialAssignment.clear();
110 [ - - ]: 0 : Trace("cdcac") << "Retrieving initial assignment:" << std::endl;
111 [ - - ]: 0 : for (const auto& var : d_variableOrdering)
112 : : {
113 : 0 : Node v = getConstraints().varMapper()(var);
114 : 0 : Node val = model.computeConcreteModelValue(v);
115 : 0 : poly::Value value = node_to_value(val, ran_variable);
116 [ - - ]: 0 : Trace("cdcac") << "\t" << get_stream_variable(var) << " = " << value
117 : 0 : << std::endl;
118 : 0 : d_initialAssignment.emplace_back(value);
119 : 0 : }
120 : : }
121 : 6412 : Constraints& CDCAC::getConstraints() { return d_constraints; }
122 : 0 : const Constraints& CDCAC::getConstraints() const { return d_constraints; }
123 : :
124 : 384 : const poly::Assignment& CDCAC::getModel() const { return d_assignment; }
125 : :
126 : 128 : const std::vector<poly::Variable>& CDCAC::getVariableOrdering() const
127 : : {
128 : 128 : return d_variableOrdering;
129 : : }
130 : :
131 : 2665 : std::vector<CACInterval> CDCAC::getUnsatIntervals(std::size_t cur_variable)
132 : : {
133 : 2665 : std::vector<CACInterval> res;
134 : 2665 : LazardEvaluation le(statisticsRegistry(), nodeManager()->getPolyContext());
135 : 2665 : prepareRootIsolation(le, cur_variable);
136 [ + + ]: 47566 : for (const auto& c : d_constraints.getConstraints())
137 : : {
138 : 44901 : const poly::Polynomial& p = std::get<0>(c);
139 : 44901 : poly::SignCondition sc = std::get<1>(c);
140 : 44901 : const Node& n = std::get<2>(c);
141 : :
142 [ + + ]: 44901 : if (main_variable(p) != d_variableOrdering[cur_variable])
143 : : {
144 : : // Constraint is in another variable, ignore it.
145 : 33129 : continue;
146 : : }
147 : :
148 [ + - ]: 23544 : Trace("cdcac") << "Infeasible intervals for " << p << " " << sc
149 : 11772 : << " 0 over " << d_assignment << std::endl;
150 : 11772 : std::vector<poly::Interval> intervals;
151 [ + + ]: 11772 : if (options().arith.nlCovLifting == options::nlCovLiftingMode::LAZARD)
152 : : {
153 : 572 : intervals = le.infeasibleRegions(p, sc);
154 [ - + ]: 572 : if (TraceIsOn("cdcac"))
155 : : {
156 : 0 : auto reference = poly::infeasible_regions(p, d_assignment, sc);
157 [ - - ]: 0 : Trace("cdcac") << "Lazard: " << intervals << std::endl;
158 [ - - ]: 0 : Trace("cdcac") << "Regular: " << reference << std::endl;
159 : 0 : }
160 : : }
161 : : else
162 : : {
163 : 11200 : intervals = poly::infeasible_regions(p, d_assignment, sc);
164 : : }
165 [ + + ]: 24485 : for (const auto& i : intervals)
166 : : {
167 [ + - ]: 12713 : Trace("cdcac") << "-> " << i << std::endl;
168 : 12713 : PolyVector l, u, m, d;
169 : 12713 : m.add(p);
170 : 12713 : m.pushDownPolys(d, d_variableOrdering[cur_variable]);
171 [ + + ]: 12713 : if (!is_minus_infinity(get_lower(i))) l = m;
172 [ + + ]: 12713 : if (!is_plus_infinity(get_upper(i))) u = m;
173 : 25426 : res.emplace_back(CACInterval{d_nextIntervalId++, i, l, u, m, d, {n}});
174 [ + + ]: 12713 : if (isProofEnabled())
175 : : {
176 : 9584 : d_proof->addDirect(
177 : 4792 : d_constraints.varMapper()(d_variableOrdering[cur_variable]),
178 : : d_constraints.varMapper(),
179 : : p,
180 : 2396 : d_assignment,
181 : : i,
182 : : n,
183 : 2396 : res.back().d_id);
184 : : }
185 : 12713 : }
186 : 11772 : }
187 : 2665 : pruneRedundantIntervals(res);
188 : 5330 : return res;
189 : 2665 : }
190 : :
191 : 4853 : bool CDCAC::sampleOutsideWithInitial(const std::vector<CACInterval>& infeasible,
192 : : poly::Value& sample,
193 : : std::size_t cur_variable)
194 : : {
195 : 4853 : if (options().arith.nlCovLinearModel != options::nlCovLinearModelMode::NONE
196 [ - + ][ - - ]: 4853 : && cur_variable < d_initialAssignment.size())
[ - + ]
197 : : {
198 : 0 : const poly::Value& suggested = d_initialAssignment[cur_variable];
199 [ - - ]: 0 : for (const auto& i : infeasible)
200 : : {
201 [ - - ]: 0 : if (poly::contains(i.d_interval, suggested))
202 : : {
203 : 0 : if (options().arith.nlCovLinearModel
204 [ - - ]: 0 : == options::nlCovLinearModelMode::INITIAL)
205 : : {
206 : 0 : d_initialAssignment.clear();
207 : : }
208 : 0 : return sampleOutside(infeasible, sample);
209 : : }
210 : : }
211 [ - - ]: 0 : Trace("cdcac") << "Using suggested initial value" << std::endl;
212 : 0 : sample = suggested;
213 : 0 : return true;
214 : : }
215 : 4853 : return sampleOutside(infeasible, sample);
216 : : }
217 : :
218 : : namespace {
219 : :
220 : : /**
221 : : * This method follows the projection operator as detailed in algorithm 6 of
222 : : * 10.1016/j.jlamp.2020.100633, which mostly follows the projection operator due
223 : : * to McCallum. It uses all coefficients until one is either constant or does
224 : : * not vanish over the current assignment.
225 : : */
226 : 4064 : PolyVector requiredCoefficientsOriginal(const poly::Polynomial& p,
227 : : const poly::Assignment& assignment)
228 : : {
229 : 4064 : PolyVector res;
230 [ + - ]: 4066 : for (long deg = degree(p); deg >= 0; --deg)
231 : : {
232 : 4066 : auto coeff = coefficient(p, deg);
233 [ - + ][ - + ]: 4066 : Assert(poly::is_constant(coeff)
[ - - ]
234 : : == lp_polynomial_is_constant(coeff.get_internal()));
235 [ + + ]: 4066 : if (poly::is_constant(coeff)) break;
236 : 516 : res.add(coeff);
237 [ + + ]: 516 : if (evaluate_constraint(coeff, assignment, poly::SignCondition::NE))
238 : : {
239 : 514 : break;
240 : : }
241 [ + + ]: 4066 : }
242 : 4064 : return res;
243 : 0 : }
244 : :
245 : : /**
246 : : * This method follows the original projection operator due to Lazard from
247 : : * section 3 of 10.1007/978-1-4612-2628-4_29. It uses the leading and trailing
248 : : * coefficient, and makes some limited efforts to avoid them in certain cases:
249 : : * We avoid the leading coefficient if it is constant. We avoid the trailing
250 : : * coefficient if the leading coefficient does not vanish over the current
251 : : * assignment and the trailing coefficient is not constant.
252 : : */
253 : 0 : PolyVector requiredCoefficientsLazard(const poly::Polynomial& p,
254 : : const poly::Assignment& assignment)
255 : : {
256 : 0 : PolyVector res;
257 : 0 : auto lc = poly::leading_coefficient(p);
258 [ - - ]: 0 : if (poly::is_constant(lc)) return res;
259 : 0 : res.add(lc);
260 [ - - ]: 0 : if (evaluate_constraint(lc, assignment, poly::SignCondition::NE)) return res;
261 : 0 : auto tc = poly::coefficient(p, 0);
262 [ - - ]: 0 : if (poly::is_constant(tc)) return res;
263 : 0 : res.add(tc);
264 : 0 : return res;
265 : 0 : }
266 : :
267 : : /**
268 : : * This method follows the enhancements from 10.1007/978-3-030-60026-6_8 for the
269 : : * projection operator due to Lazard, more specifically Section 3.3 and
270 : : * Definition 4. In essence, we can skip the trailing coefficient if we can show
271 : : * that p is not nullified by any n-1 dimensional point. The statement in the
272 : : * paper is slightly more general, but there is no simple way to have such a
273 : : * procedure T here. We simply try to show that T(f) = {} by using the extended
274 : : * rewriter to simplify (and (= f_i 0)) (f_i being the coefficients of f) to
275 : : * false.
276 : : */
277 : 0 : PolyVector requiredCoefficientsLazardModified(
278 : : NodeManager* nm,
279 : : const poly::Polynomial& p,
280 : : const poly::Assignment& assignment,
281 : : VariableMapper& vm,
282 : : Rewriter* rewriter)
283 : : {
284 : 0 : PolyVector res;
285 : 0 : auto lc = poly::leading_coefficient(p);
286 : : // if leading coefficient is constant
287 [ - - ]: 0 : if (poly::is_constant(lc)) return res;
288 : : // add leading coefficient
289 : 0 : res.add(lc);
290 : 0 : auto tc = poly::coefficient(p, 0);
291 : : // if trailing coefficient is constant
292 [ - - ]: 0 : if (poly::is_constant(tc)) return res;
293 : : // if leading coefficient does not vanish over the current assignment
294 [ - - ]: 0 : if (evaluate_constraint(lc, assignment, poly::SignCondition::NE)) return res;
295 : :
296 : : // construct phi := (and (= p_i 0)) with p_i the coefficients of p
297 : 0 : std::vector<Node> conditions;
298 : 0 : auto zero = nm->mkConstReal(Rational(0));
299 [ - - ]: 0 : for (const auto& coeff : poly::coefficients(p))
300 : : {
301 : 0 : conditions.emplace_back(NodeManager::mkNode(
302 : 0 : Kind::EQUAL, nl::as_cvc_polynomial(nm, coeff, vm), zero));
303 : 0 : }
304 : : // if phi is false (i.e. p can not vanish)
305 : 0 : Node rewritten = rewriter->extendedRewrite(nm->mkAnd(conditions));
306 [ - - ]: 0 : if (rewritten.isConst())
307 : : {
308 : 0 : Assert(rewritten.getKind() == Kind::CONST_BOOLEAN);
309 : 0 : Assert(!rewritten.getConst<bool>());
310 : 0 : return res;
311 : : }
312 : : // otherwise add trailing coefficient as well
313 : 0 : res.add(tc);
314 : 0 : return res;
315 : 0 : }
316 : :
317 : : } // namespace
318 : :
319 : 4064 : PolyVector CDCAC::requiredCoefficients(const poly::Polynomial& p)
320 : : {
321 [ - + ]: 4064 : if (TraceIsOn("cdcac::projection"))
322 : : {
323 [ - - ]: 0 : Trace("cdcac::projection")
324 : 0 : << "Poly: " << p << " over " << d_assignment << std::endl;
325 [ - - ]: 0 : Trace("cdcac::projection")
326 : 0 : << "Lazard: " << requiredCoefficientsLazard(p, d_assignment)
327 : 0 : << std::endl;
328 [ - - ]: 0 : Trace("cdcac::projection")
329 : 0 : << "LMod: "
330 : 0 : << requiredCoefficientsLazardModified(d_env.getNodeManager(),
331 : : p,
332 : 0 : d_assignment,
333 : : d_constraints.varMapper(),
334 : 0 : d_env.getRewriter())
335 : 0 : << std::endl;
336 [ - - ]: 0 : Trace("cdcac::projection")
337 : 0 : << "Original: " << requiredCoefficientsOriginal(p, d_assignment)
338 : 0 : << std::endl;
339 : : }
340 [ + - ][ - - ]: 4064 : switch (options().arith.nlCovProjection)
341 : : {
342 : 4064 : case options::nlCovProjectionMode::MCCALLUM:
343 : 4064 : return requiredCoefficientsOriginal(p, d_assignment);
344 : 0 : case options::nlCovProjectionMode::LAZARD:
345 : 0 : return requiredCoefficientsLazard(p, d_assignment);
346 : 0 : case options::nlCovProjectionMode::LAZARDMOD:
347 : 0 : return requiredCoefficientsLazardModified(d_env.getNodeManager(),
348 : : p,
349 : 0 : d_assignment,
350 : : d_constraints.varMapper(),
351 : 0 : d_env.getRewriter());
352 : 0 : default:
353 : 0 : DebugUnhandled();
354 : : return requiredCoefficientsOriginal(p, d_assignment);
355 : : }
356 : : }
357 : :
358 : 2102 : PolyVector CDCAC::constructCharacterization(std::vector<CACInterval>& intervals)
359 : : {
360 [ - + ][ - + ]: 2102 : Assert(!intervals.empty()) << "A covering can not be empty";
[ - - ]
361 [ + - ]: 2102 : Trace("cdcac") << "Constructing characterization now" << std::endl;
362 : 2102 : PolyVector res;
363 : :
364 [ + + ]: 3831 : for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
365 : : {
366 : 1729 : coverings::makeFinestSquareFreeBasis(intervals[i], intervals[i + 1]);
367 : : }
368 : :
369 [ + + ]: 5933 : for (const auto& i : intervals)
370 : : {
371 [ + - ]: 3831 : Trace("cdcac") << "Considering " << i.d_interval << std::endl;
372 [ + - ]: 7662 : Trace("cdcac") << "-> " << i.d_lowerPolys << " / " << i.d_upperPolys
373 : 0 : << " and " << i.d_mainPolys << " / " << i.d_downPolys
374 : 3831 : << std::endl;
375 [ + - ]: 3831 : Trace("cdcac") << "-> " << i.d_origins << std::endl;
376 [ + + ]: 17048 : for (const auto& p : i.d_downPolys)
377 : : {
378 : : // Add all polynomial from lower levels.
379 : 13217 : res.add(p);
380 : : }
381 [ + + ]: 7895 : for (const auto& p : i.d_mainPolys)
382 : : {
383 [ + - ]: 8128 : Trace("cdcac::projection")
384 [ - + ][ - - ]: 4064 : << "Discriminant of " << p << " -> " << discriminant(p) << std::endl;
385 : : // Add all discriminants
386 : 4064 : res.add(discriminant(p));
387 : :
388 : : // Add pairwise resultants
389 [ + + ]: 8828 : for (const auto& q : i.d_mainPolys)
390 : : {
391 : : // avoid symmetric duplicates
392 [ + + ]: 4764 : if (p >= q) continue;
393 : 350 : res.add(resultant(p, q));
394 : : }
395 : :
396 [ + + ]: 4932 : for (const auto& q : requiredCoefficients(p))
397 : : {
398 : : // Add all required coefficients
399 [ + - ]: 1736 : Trace("cdcac::projection")
400 : 868 : << "Coeff of " << p << " -> " << q << std::endl;
401 : 868 : res.add(q);
402 : 4064 : }
403 [ + + ]: 6027 : for (const auto& q : i.d_lowerPolys)
404 : : {
405 [ + + ]: 1963 : if (p == q) continue;
406 : : // Check whether p(s \times a) = 0 for some a <= l
407 [ - + ]: 234 : if (!hasRootBelow(q, get_lower(i.d_interval))) continue;
408 [ + - ]: 468 : Trace("cdcac::projection") << "Resultant of " << p << " and " << q
409 [ - + ][ - - ]: 234 : << " -> " << resultant(p, q) << std::endl;
410 : 234 : res.add(resultant(p, q));
411 : : }
412 [ + + ]: 6057 : for (const auto& q : i.d_upperPolys)
413 : : {
414 [ + + ]: 1993 : if (p == q) continue;
415 : : // Check whether p(s \times a) = 0 for some a >= u
416 [ - + ]: 260 : if (!hasRootAbove(q, get_upper(i.d_interval))) continue;
417 [ + - ]: 520 : Trace("cdcac::projection") << "Resultant of " << p << " and " << q
418 [ - + ][ - - ]: 260 : << " -> " << resultant(p, q) << std::endl;
419 : 260 : res.add(resultant(p, q));
420 : : }
421 : : }
422 : : }
423 : :
424 [ + + ]: 3831 : for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
425 : : {
426 : : // Add resultants of consecutive intervals.
427 [ + + ]: 3462 : for (const auto& p : intervals[i].d_upperPolys)
428 : : {
429 [ + + ]: 3466 : for (const auto& q : intervals[i + 1].d_lowerPolys)
430 : : {
431 [ + - ]: 3466 : Trace("cdcac::projection") << "Resultant of " << p << " and " << q
432 [ - + ][ - - ]: 1733 : << " -> " << resultant(p, q) << std::endl;
433 : 1733 : res.add(resultant(p, q));
434 : : }
435 : : }
436 : : }
437 : :
438 : 2102 : res.reduce();
439 : 2102 : res.makeFinestSquareFreeBasis();
440 : :
441 : 2102 : return res;
442 : 0 : }
443 : :
444 : 2102 : CACInterval CDCAC::intervalFromCharacterization(
445 : : const PolyVector& characterization,
446 : : std::size_t cur_variable,
447 : : const poly::Value& sample)
448 : : {
449 : 2102 : PolyVector l;
450 : 2102 : PolyVector u;
451 : 2102 : PolyVector m;
452 : 2102 : PolyVector d;
453 : :
454 [ + + ]: 11804 : for (const auto& p : characterization)
455 : : {
456 : : // Add polynomials to main
457 : 9702 : m.add(p);
458 : : }
459 : : // Push lower-dimensional polys to down
460 : 2102 : m.pushDownPolys(d, d_variableOrdering[cur_variable]);
461 : :
462 : : // Collect -oo, all roots, oo
463 : :
464 : 2102 : LazardEvaluation le(statisticsRegistry(), nodeManager()->getPolyContext());
465 : 2102 : prepareRootIsolation(le, cur_variable);
466 : 2102 : std::vector<poly::Value> roots;
467 : 2102 : roots.emplace_back(poly::Value::minus_infty());
468 [ + + ]: 4606 : for (const auto& p : m)
469 : : {
470 [ + - ]: 5008 : Trace("cdcac") << "Isolating real roots of " << p << " over "
471 : 2504 : << d_assignment << std::endl;
472 : :
473 : 2504 : auto tmp = isolateRealRoots(le, p);
474 : 2504 : roots.insert(roots.end(), tmp.begin(), tmp.end());
475 : 2504 : }
476 : 2102 : roots.emplace_back(poly::Value::plus_infty());
477 : 2102 : std::sort(roots.begin(), roots.end());
478 : :
479 : : // Now find the interval bounds
480 : 2102 : poly::Value lower;
481 : 2102 : poly::Value upper;
482 [ + - ]: 5173 : for (std::size_t i = 0, n = roots.size(); i < n; ++i)
483 : : {
484 [ + + ]: 5173 : if (sample < roots[i])
485 : : {
486 : 1544 : lower = roots[i - 1];
487 : 1544 : upper = roots[i];
488 : 1544 : break;
489 : : }
490 [ + + ]: 3629 : if (roots[i] == sample)
491 : : {
492 : 558 : lower = sample;
493 : 558 : upper = sample;
494 : 558 : break;
495 : : }
496 : : }
497 [ + - ][ + - ]: 2102 : Assert(!is_none(lower) && !is_none(upper));
[ - + ][ - + ]
[ - - ]
498 : :
499 [ + + ]: 2102 : if (!is_minus_infinity(lower))
500 : : {
501 : : // Identify polynomials that have a root at the lower bound
502 : 1359 : d_assignment.set(d_variableOrdering[cur_variable], lower);
503 [ + + ]: 3109 : for (const auto& p : m)
504 : : {
505 [ + - ]: 3500 : Trace("cdcac") << "Evaluating " << p << " = 0 over " << d_assignment
506 : 1750 : << std::endl;
507 [ + + ]: 1750 : if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
508 : : {
509 : 1359 : l.add(p, true);
510 : : }
511 : : }
512 : 1359 : d_assignment.unset(d_variableOrdering[cur_variable]);
513 : : }
514 [ + + ]: 2102 : if (!is_plus_infinity(upper))
515 : : {
516 : : // Identify polynomials that have a root at the upper bound
517 : 1332 : d_assignment.set(d_variableOrdering[cur_variable], upper);
518 [ + + ]: 3100 : for (const auto& p : m)
519 : : {
520 [ + - ]: 3536 : Trace("cdcac") << "Evaluating " << p << " = 0 over " << d_assignment
521 : 1768 : << std::endl;
522 [ + + ]: 1768 : if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
523 : : {
524 : 1336 : u.add(p, true);
525 : : }
526 : : }
527 : 1332 : d_assignment.unset(d_variableOrdering[cur_variable]);
528 : : }
529 : :
530 [ + + ]: 2102 : if (lower == upper)
531 : : {
532 : : // construct a point interval
533 : 558 : return CACInterval{d_nextIntervalId++,
534 : : poly::Interval(lower, false, upper, false),
535 : : l,
536 : : u,
537 : : m,
538 : : d,
539 : 558 : {}};
540 : : }
541 : : else
542 : : {
543 : : // construct an open interval
544 [ - + ][ - + ]: 1544 : Assert(lower < upper);
[ - - ]
545 : 1544 : return CACInterval{d_nextIntervalId++,
546 : : poly::Interval(lower, true, upper, true),
547 : : l,
548 : : u,
549 : : m,
550 : : d,
551 : 1544 : {}};
552 : : }
553 : 2102 : }
554 : :
555 : 2665 : std::vector<CACInterval> CDCAC::getUnsatCoverImpl(std::size_t curVariable,
556 : : bool returnFirstInterval)
557 : : {
558 : 2665 : d_env.getResourceManager()->spendResource(Resource::ArithNlCoveringStep);
559 [ + - ]: 5330 : Trace("cdcac") << "Looking for unsat cover for "
560 : 2665 : << get_stream_variable(d_variableOrdering[curVariable])
561 : 2665 : << std::endl;
562 : 2665 : std::vector<CACInterval> intervals = getUnsatIntervals(curVariable);
563 : :
564 [ - + ]: 2665 : if (TraceIsOn("cdcac"))
565 : : {
566 [ - - ]: 0 : Trace("cdcac") << "Unsat intervals for "
567 : 0 : << get_stream_variable(d_variableOrdering[curVariable])
568 : 0 : << ":" << std::endl;
569 [ - - ]: 0 : for (const auto& i : intervals)
570 : : {
571 [ - - ]: 0 : Trace("cdcac") << "-> " << i.d_interval << " from " << i.d_origins
572 : 0 : << std::endl;
573 : : }
574 : : }
575 : 2665 : poly::Value sample;
576 : :
577 [ + + ]: 4853 : while (sampleOutsideWithInitial(intervals, sample, curVariable))
578 : : {
579 [ + + ]: 2614 : if (!checkIntegrality(curVariable, sample))
580 : : {
581 : : // the variable is integral, but the sample is not.
582 [ + - ]: 172 : Trace("cdcac") << "Used " << sample << " for integer variable "
583 : 86 : << get_stream_variable(d_variableOrdering[curVariable])
584 : 86 : << std::endl;
585 : 86 : auto newInterval = buildIntegralityInterval(curVariable, sample);
586 [ + - ]: 172 : Trace("cdcac") << "Adding integrality interval " << newInterval.d_interval
587 : 86 : << std::endl;
588 : 86 : intervals.emplace_back(newInterval);
589 : 86 : pruneRedundantIntervals(intervals);
590 : 86 : continue;
591 : 86 : }
592 : 2528 : d_assignment.set(d_variableOrdering[curVariable], sample);
593 [ + - ]: 2528 : Trace("cdcac") << "Sample: " << d_assignment << std::endl;
594 [ + + ]: 2528 : if (curVariable == d_variableOrdering.size() - 1)
595 : : {
596 : : // We have a full assignment. SAT!
597 [ + - ]: 142 : Trace("cdcac") << "Found full assignment: " << d_assignment << std::endl;
598 : 426 : return {};
599 : : }
600 [ + + ]: 2386 : if (isProofEnabled())
601 : : {
602 : 1101 : d_proof->startScope();
603 : 1101 : d_proof->startRecursive();
604 : : }
605 : : // Recurse to next variable
606 : 2386 : auto cov = getUnsatCoverImpl(curVariable + 1);
607 [ + + ]: 2386 : if (cov.empty())
608 : : {
609 : : // Found SAT!
610 [ + - ]: 284 : Trace("cdcac") << "SAT!" << std::endl;
611 : 284 : return {};
612 : : }
613 [ + - ]: 2102 : Trace("cdcac") << "Refuting Sample: " << d_assignment << std::endl;
614 : 2102 : auto characterization = constructCharacterization(cov);
615 [ + - ]: 2102 : Trace("cdcac") << "Characterization: " << characterization << std::endl;
616 : :
617 : 2102 : d_assignment.unset(d_variableOrdering[curVariable]);
618 : :
619 [ + - ]: 2102 : Trace("cdcac") << "Building interval..." << std::endl;
620 : : auto newInterval =
621 : 2102 : intervalFromCharacterization(characterization, curVariable, sample);
622 [ + - ]: 2102 : Trace("cdcac") << "New interval: " << newInterval.d_interval << std::endl;
623 : 2102 : newInterval.d_origins = collectConstraints(cov);
624 : 2102 : intervals.emplace_back(newInterval);
625 [ + + ]: 2102 : if (isProofEnabled())
626 : : {
627 : 1087 : d_proof->endRecursive(newInterval.d_id);
628 : : auto cell = d_proof->constructCell(
629 : 1087 : d_constraints.varMapper()(d_variableOrdering[curVariable]),
630 : : newInterval,
631 : 1087 : d_assignment,
632 : : sample,
633 : 2174 : d_constraints.varMapper());
634 : 1087 : d_proof->endScope(cell);
635 : 1087 : }
636 : :
637 [ - + ]: 2102 : if (returnFirstInterval)
638 : : {
639 : 0 : return intervals;
640 : : }
641 : :
642 [ + - ]: 2102 : Trace("cdcac") << "Added " << intervals.back().d_interval << std::endl;
643 [ + - ]: 4204 : Trace("cdcac") << "\tlower: " << intervals.back().d_lowerPolys
644 : 2102 : << std::endl;
645 [ + - ]: 4204 : Trace("cdcac") << "\tupper: " << intervals.back().d_upperPolys
646 : 2102 : << std::endl;
647 [ + - ]: 4204 : Trace("cdcac") << "\tmain: " << intervals.back().d_mainPolys
648 : 2102 : << std::endl;
649 [ + - ]: 4204 : Trace("cdcac") << "\tdown: " << intervals.back().d_downPolys
650 : 2102 : << std::endl;
651 [ + - ]: 2102 : Trace("cdcac") << "\torigins: " << intervals.back().d_origins << std::endl;
652 : :
653 : : // Remove redundant intervals
654 : 2102 : pruneRedundantIntervals(intervals);
655 [ + - ][ + - ]: 2386 : }
[ + + ]
656 : :
657 [ - + ]: 2239 : if (TraceIsOn("cdcac"))
658 : : {
659 [ - - ]: 0 : Trace("cdcac") << "Returning intervals for "
660 : 0 : << get_stream_variable(d_variableOrdering[curVariable])
661 : 0 : << ":" << std::endl;
662 [ - - ]: 0 : for (const auto& i : intervals)
663 : : {
664 [ - - ]: 0 : Trace("cdcac") << "-> " << i.d_interval << std::endl;
665 : : }
666 : : }
667 : 2239 : return intervals;
668 : 2665 : }
669 : :
670 : 279 : std::vector<CACInterval> CDCAC::getUnsatCover(bool returnFirstInterval)
671 : : {
672 [ + + ]: 279 : if (isProofEnabled())
673 : : {
674 : 71 : d_proof->startRecursive();
675 : : }
676 : 279 : auto res = getUnsatCoverImpl(0, returnFirstInterval);
677 [ + + ]: 279 : if (isProofEnabled())
678 : : {
679 : 71 : d_proof->endRecursive(0);
680 : : }
681 : 279 : return res;
682 : 0 : }
683 : :
684 : 273 : void CDCAC::startNewProof()
685 : : {
686 [ + + ]: 273 : if (isProofEnabled())
687 : : {
688 : 71 : d_proof->startNewProof();
689 : : }
690 : 273 : }
691 : :
692 : 136 : ProofGenerator* CDCAC::closeProof(const std::vector<Node>& assertions)
693 : : {
694 [ + + ]: 136 : if (isProofEnabled())
695 : : {
696 : 61 : d_proof->endScope(assertions);
697 : 61 : return d_proof->getProofGenerator();
698 : : }
699 : 75 : return nullptr;
700 : : }
701 : :
702 : 2614 : bool CDCAC::checkIntegrality(std::size_t cur_variable, const poly::Value& value)
703 : : {
704 : 2614 : Node var = d_constraints.varMapper()(d_variableOrdering[cur_variable]);
705 [ + + ]: 2614 : if (!var.getType().isInteger())
706 : : {
707 : : // variable is not integral
708 : 2376 : return true;
709 : : }
710 : 238 : return poly::represents_integer(value);
711 : 2614 : }
712 : :
713 : 86 : CACInterval CDCAC::buildIntegralityInterval(std::size_t cur_variable,
714 : : const poly::Value& value)
715 : : {
716 : 86 : poly::Variable var = d_variableOrdering[cur_variable];
717 : 86 : poly::Integer below = poly::floor(value);
718 : 86 : poly::Integer above = poly::ceil(value);
719 : 86 : const poly::Context& ctx = nodeManager()->getPolyContext();
720 : 86 : poly::Polynomial pvar = poly::Polynomial(ctx, var);
721 : 86 : poly::Polynomial pbelow = poly::Polynomial(ctx, below);
722 : 86 : poly::Polynomial pabove = poly::Polynomial(ctx, above);
723 : : // construct var \in (below, above)
724 : 86 : return CACInterval{d_nextIntervalId++,
725 : : poly::Interval(below, above),
726 : : {pvar - pbelow},
727 : : {pvar - pabove},
728 : : {pvar - pbelow, pvar - pabove},
729 : : {},
730 : 516 : {}};
731 : 86 : }
732 : :
733 : 260 : bool CDCAC::hasRootAbove(const poly::Polynomial& p,
734 : : const poly::Value& val) const
735 : : {
736 : 260 : auto roots = poly::isolate_real_roots(p, d_assignment);
737 : 260 : return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
738 : 276 : return r >= val;
739 : 520 : });
740 : 260 : }
741 : :
742 : 234 : bool CDCAC::hasRootBelow(const poly::Polynomial& p,
743 : : const poly::Value& val) const
744 : : {
745 : 234 : auto roots = poly::isolate_real_roots(p, d_assignment);
746 : 234 : return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
747 : 234 : return r <= val;
748 : 468 : });
749 : 234 : }
750 : :
751 : 4853 : void CDCAC::pruneRedundantIntervals(std::vector<CACInterval>& intervals)
752 : : {
753 : 4853 : cleanIntervals(intervals);
754 [ - + ]: 4853 : if (options().arith.nlCovPrune)
755 : : {
756 [ - - ]: 0 : if (TraceIsOn("cdcac"))
757 : : {
758 : 0 : auto copy = intervals;
759 : 0 : removeRedundantIntervals(intervals);
760 [ - - ]: 0 : if (copy.size() != intervals.size())
761 : : {
762 [ - - ]: 0 : Trace("cdcac") << "Before pruning:";
763 : 0 : for (const auto& i : copy) Trace("cdcac") << " " << i.d_interval;
764 [ - - ]: 0 : Trace("cdcac") << std::endl;
765 [ - - ]: 0 : Trace("cdcac") << "After pruning: ";
766 : 0 : for (const auto& i : intervals) Trace("cdcac") << " " << i.d_interval;
767 [ - - ]: 0 : Trace("cdcac") << std::endl;
768 : : }
769 : 0 : }
770 : : else
771 : : {
772 : 0 : removeRedundantIntervals(intervals);
773 : : }
774 : : }
775 [ + + ]: 4853 : if (isProofEnabled())
776 : : {
777 : 2259 : d_proof->pruneChildren([&intervals](std::size_t id) {
778 [ - + ]: 4627 : if (id == 0) return false;
779 : 4627 : return std::find_if(intervals.begin(),
780 : : intervals.end(),
781 : 6901 : [id](const CACInterval& i) { return i.d_id == id; })
782 : 9254 : == intervals.end();
783 : : });
784 : : }
785 : 4853 : }
786 : :
787 : 4767 : void CDCAC::prepareRootIsolation(LazardEvaluation& le,
788 : : size_t cur_variable) const
789 : : {
790 [ + + ]: 4767 : if (options().arith.nlCovLifting == options::nlCovLiftingMode::LAZARD)
791 : : {
792 [ + + ]: 1062 : for (size_t vid = 0; vid < cur_variable; ++vid)
793 : : {
794 : 836 : const auto& val = d_assignment.get(d_variableOrdering[vid]);
795 : 836 : le.add(d_variableOrdering[vid], val);
796 : : }
797 : 226 : le.addFreeVariable(d_variableOrdering[cur_variable]);
798 : : }
799 : 4767 : }
800 : :
801 : 2504 : std::vector<poly::Value> CDCAC::isolateRealRoots(
802 : : LazardEvaluation& le, const poly::Polynomial& p) const
803 : : {
804 [ + + ]: 2504 : if (options().arith.nlCovLifting == options::nlCovLiftingMode::LAZARD)
805 : : {
806 : 130 : return le.isolateRealRoots(p);
807 : : }
808 : 2374 : return poly::isolate_real_roots(p, d_assignment);
809 : : }
810 : :
811 : : } // namespace coverings
812 : : } // namespace nl
813 : : } // namespace arith
814 : : } // namespace theory
815 : : } // namespace cvc5::internal
816 : :
817 : : #endif
|
