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