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