1 /********************* */
4 ** Top contributors (to current version):
6 ** This file is part of the CVC4 project.
7 ** Copyright (c) 2009-2020 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.\endverbatim
12 ** \brief Implements the CDCAC approach.
14 ** Implements the CDCAC approach as described in
15 ** https://arxiv.org/pdf/2003.05633.pdf.
18 #include "theory/arith/nl/cad/cdcac.h"
22 #include "theory/arith/nl/cad/projections.h"
23 #include "theory/arith/nl/cad/variable_ordering.h"
26 /** Generic streaming operator for std::vector. */
28 std::ostream
& operator<<(std::ostream
& os
, const std::vector
<T
>& v
)
30 CVC4::container_to_stream(os
, v
);
42 /** Removed duplicates from a vector. */
44 void removeDuplicates(std::vector
<T
>& v
)
46 std::sort(v
.begin(), v
.end());
47 v
.erase(std::unique(v
.begin(), v
.end()), v
.end());
53 CDCAC::CDCAC(const std::vector
<poly::Variable
>& ordering
)
54 : d_variableOrdering(ordering
)
60 d_constraints
.reset();
64 void CDCAC::computeVariableOrdering()
66 // Actually compute the variable ordering
67 d_variableOrdering
= d_varOrder(d_constraints
.getConstraints(),
68 VariableOrderingStrategy::BROWN
);
69 Trace("cdcac") << "Variable ordering is now " << d_variableOrdering
72 // Write variable ordering back to libpoly.
73 lp_variable_order_t
* vo
= poly::Context::get_context().get_variable_order();
74 lp_variable_order_clear(vo
);
75 for (const auto& v
: d_variableOrdering
)
77 lp_variable_order_push(vo
, v
.get_internal());
81 void CDCAC::retrieveInitialAssignment(NlModel
& model
, const Node
& ran_variable
)
83 d_initialAssignment
.clear();
84 Trace("cdcac") << "Retrieving initial assignment:" << std::endl
;
85 for (const auto& var
: d_variableOrdering
)
87 Node v
= getConstraints().varMapper()(var
);
88 Node val
= model
.computeConcreteModelValue(v
);
89 poly::Value value
= node_to_value(val
, ran_variable
);
90 Trace("cdcac") << "\t" << var
<< " = " << value
<< std::endl
;
91 d_initialAssignment
.emplace_back(value
);
94 Constraints
& CDCAC::getConstraints() { return d_constraints
; }
95 const Constraints
& CDCAC::getConstraints() const { return d_constraints
; }
97 const poly::Assignment
& CDCAC::getModel() const { return d_assignment
; }
99 const std::vector
<poly::Variable
>& CDCAC::getVariableOrdering() const
101 return d_variableOrdering
;
104 std::vector
<CACInterval
> CDCAC::getUnsatIntervals(
105 std::size_t cur_variable
) const
107 std::vector
<CACInterval
> res
;
108 for (const auto& c
: d_constraints
.getConstraints())
110 const poly::Polynomial
& p
= std::get
<0>(c
);
111 poly::SignCondition sc
= std::get
<1>(c
);
112 const Node
& n
= std::get
<2>(c
);
114 if (main_variable(p
) != d_variableOrdering
[cur_variable
])
116 // Constraint is in another variable, ignore it.
120 Trace("cdcac") << "Infeasible intervals for " << p
<< " " << sc
121 << " 0 over " << d_assignment
<< std::endl
;
122 auto intervals
= infeasible_regions(p
, d_assignment
, sc
);
123 for (const auto& i
: intervals
)
125 Trace("cdcac") << "-> " << i
<< std::endl
;
126 std::vector
<poly::Polynomial
> l
, u
, m
, d
;
127 if (!is_minus_infinity(get_lower(i
))) l
.emplace_back(p
);
128 if (!is_plus_infinity(get_upper(i
))) u
.emplace_back(p
);
130 res
.emplace_back(CACInterval
{i
, l
, u
, m
, d
, {n
}});
137 bool CDCAC::sampleOutsideWithInitial(const std::vector
<CACInterval
>& infeasible
,
139 std::size_t cur_variable
)
141 if (cur_variable
< d_initialAssignment
.size())
143 const poly::Value
& suggested
= d_initialAssignment
[cur_variable
];
144 for (const auto& i
: infeasible
)
146 if (poly::contains(i
.d_interval
, suggested
))
148 d_initialAssignment
.clear();
149 return sampleOutside(infeasible
, sample
);
152 Trace("cdcac") << "Using suggested initial value" << std::endl
;
156 return sampleOutside(infeasible
, sample
);
159 std::vector
<poly::Polynomial
> CDCAC::requiredCoefficients(
160 const poly::Polynomial
& p
) const
162 std::vector
<poly::Polynomial
> res
;
163 for (long deg
= degree(p
); deg
>= 0; --deg
)
165 auto coeff
= coefficient(p
, deg
);
166 if (lp_polynomial_is_constant(coeff
.get_internal())) break;
167 res
.emplace_back(coeff
);
168 if (evaluate_constraint(coeff
, d_assignment
, poly::SignCondition::NE
))
176 std::vector
<poly::Polynomial
> CDCAC::constructCharacterization(
177 std::vector
<CACInterval
>& intervals
)
179 Assert(!intervals
.empty()) << "A covering can not be empty";
180 Trace("cdcac") << "Constructing characterization now" << std::endl
;
181 std::vector
<poly::Polynomial
> res
;
183 for (const auto& i
: intervals
)
185 Trace("cdcac") << "Considering " << i
.d_interval
<< std::endl
;
186 Trace("cdcac") << "-> " << i
.d_lowerPolys
<< " / " << i
.d_upperPolys
187 << " and " << i
.d_mainPolys
<< " / " << i
.d_downPolys
189 Trace("cdcac") << "-> " << i
.d_origins
<< std::endl
;
190 for (const auto& p
: i
.d_downPolys
)
192 // Add all polynomial from lower levels.
193 addPolynomial(res
, p
);
195 for (const auto& p
: i
.d_mainPolys
)
197 Trace("cdcac") << "Discriminant of " << p
<< " -> " << discriminant(p
)
199 // Add all discriminants
200 addPolynomial(res
, discriminant(p
));
202 for (const auto& q
: requiredCoefficients(p
))
204 // Add all required coefficients
205 Trace("cdcac") << "Coeff of " << p
<< " -> " << q
<< std::endl
;
206 addPolynomial(res
, q
);
208 for (const auto& q
: i
.d_lowerPolys
)
210 if (p
== q
) continue;
211 // Check whether p(s \times a) = 0 for some a <= l
212 if (!hasRootBelow(q
, get_lower(i
.d_interval
))) continue;
213 Trace("cdcac") << "Resultant of " << p
<< " and " << q
<< " -> "
214 << resultant(p
, q
) << std::endl
;
215 addPolynomial(res
, resultant(p
, q
));
217 for (const auto& q
: i
.d_upperPolys
)
219 if (p
== q
) continue;
220 // Check whether p(s \times a) = 0 for some a >= u
221 if (!hasRootAbove(q
, get_upper(i
.d_interval
))) continue;
222 Trace("cdcac") << "Resultant of " << p
<< " and " << q
<< " -> "
223 << resultant(p
, q
) << std::endl
;
224 addPolynomial(res
, resultant(p
, q
));
229 for (std::size_t i
= 0, n
= intervals
.size(); i
< n
- 1; ++i
)
231 // Add resultants of consecutive intervals.
232 cad::makeFinestSquareFreeBasis(intervals
[i
].d_upperPolys
,
233 intervals
[i
+ 1].d_lowerPolys
);
234 for (const auto& p
: intervals
[i
].d_upperPolys
)
236 for (const auto& q
: intervals
[i
+ 1].d_lowerPolys
)
238 Trace("cdcac") << "Resultant of " << p
<< " and " << q
<< " -> "
239 << resultant(p
, q
) << std::endl
;
240 addPolynomial(res
, resultant(p
, q
));
245 removeDuplicates(res
);
246 makeFinestSquareFreeBasis(res
);
251 CACInterval
CDCAC::intervalFromCharacterization(
252 const std::vector
<poly::Polynomial
>& characterization
,
253 std::size_t cur_variable
,
254 const poly::Value
& sample
)
256 std::vector
<poly::Polynomial
> l
;
257 std::vector
<poly::Polynomial
> u
;
258 std::vector
<poly::Polynomial
> m
;
259 std::vector
<poly::Polynomial
> d
;
261 for (const auto& p
: characterization
)
263 // Add polynomials to either main or down
264 if (main_variable(p
) == d_variableOrdering
[cur_variable
])
274 // Collect -oo, all roots, oo
275 std::vector
<poly::Value
> roots
;
276 roots
.emplace_back(poly::Value::minus_infty());
277 for (const auto& p
: m
)
279 auto tmp
= isolate_real_roots(p
, d_assignment
);
280 roots
.insert(roots
.end(), tmp
.begin(), tmp
.end());
282 roots
.emplace_back(poly::Value::plus_infty());
283 std::sort(roots
.begin(), roots
.end());
285 // Now find the interval bounds
288 for (std::size_t i
= 0, n
= roots
.size(); i
< n
; ++i
)
290 if (sample
< roots
[i
])
292 lower
= roots
[i
- 1];
296 if (roots
[i
] == sample
)
303 Assert(!is_none(lower
) && !is_none(upper
));
305 if (!is_minus_infinity(lower
))
307 // Identify polynomials that have a root at the lower bound
308 d_assignment
.set(d_variableOrdering
[cur_variable
], lower
);
309 for (const auto& p
: m
)
311 if (evaluate_constraint(p
, d_assignment
, poly::SignCondition::EQ
))
316 d_assignment
.unset(d_variableOrdering
[cur_variable
]);
318 if (!is_plus_infinity(upper
))
320 // Identify polynomials that have a root at the upper bound
321 d_assignment
.set(d_variableOrdering
[cur_variable
], upper
);
322 for (const auto& p
: m
)
324 if (evaluate_constraint(p
, d_assignment
, poly::SignCondition::EQ
))
329 d_assignment
.unset(d_variableOrdering
[cur_variable
]);
334 // construct a point interval
336 poly::Interval(lower
, false, upper
, false), l
, u
, m
, d
, {}};
340 // construct an open interval
341 Assert(lower
< upper
);
343 poly::Interval(lower
, true, upper
, true), l
, u
, m
, d
, {}};
347 std::vector
<CACInterval
> CDCAC::getUnsatCover(std::size_t curVariable
,
348 bool returnFirstInterval
)
350 Trace("cdcac") << "Looking for unsat cover for "
351 << d_variableOrdering
[curVariable
] << std::endl
;
352 std::vector
<CACInterval
> intervals
= getUnsatIntervals(curVariable
);
354 if (Trace
.isOn("cdcac"))
356 Trace("cdcac") << "Unsat intervals for " << d_variableOrdering
[curVariable
]
358 for (const auto& i
: intervals
)
360 Trace("cdcac") << "-> " << i
.d_interval
<< " from " << i
.d_origins
366 while (sampleOutsideWithInitial(intervals
, sample
, curVariable
))
368 if (!checkIntegrality(curVariable
, sample
))
370 // the variable is integral, but the sample is not.
371 Trace("cdcac") << "Used " << sample
<< " for integer variable "
372 << d_variableOrdering
[curVariable
] << std::endl
;
373 auto newInterval
= buildIntegralityInterval(curVariable
, sample
);
374 Trace("cdcac") << "Adding integrality interval " << newInterval
.d_interval
376 intervals
.emplace_back(newInterval
);
377 cleanIntervals(intervals
);
380 d_assignment
.set(d_variableOrdering
[curVariable
], sample
);
381 Trace("cdcac") << "Sample: " << d_assignment
<< std::endl
;
382 if (curVariable
== d_variableOrdering
.size() - 1)
384 // We have a full assignment. SAT!
385 Trace("cdcac") << "Found full assignment: " << d_assignment
<< std::endl
;
388 // Recurse to next variable
389 auto cov
= getUnsatCover(curVariable
+ 1);
393 Trace("cdcac") << "SAT!" << std::endl
;
396 Trace("cdcac") << "Refuting Sample: " << d_assignment
<< std::endl
;
397 auto characterization
= constructCharacterization(cov
);
398 Trace("cdcac") << "Characterization: " << characterization
<< std::endl
;
400 d_assignment
.unset(d_variableOrdering
[curVariable
]);
403 intervalFromCharacterization(characterization
, curVariable
, sample
);
404 newInterval
.d_origins
= collectConstraints(cov
);
405 intervals
.emplace_back(newInterval
);
407 if (returnFirstInterval
)
412 Trace("cdcac") << "Added " << intervals
.back().d_interval
<< std::endl
;
413 Trace("cdcac") << "\tlower: " << intervals
.back().d_lowerPolys
415 Trace("cdcac") << "\tupper: " << intervals
.back().d_upperPolys
417 Trace("cdcac") << "\tmain: " << intervals
.back().d_mainPolys
419 Trace("cdcac") << "\tdown: " << intervals
.back().d_downPolys
421 Trace("cdcac") << "\torigins: " << intervals
.back().d_origins
<< std::endl
;
423 // Remove redundant intervals
424 cleanIntervals(intervals
);
427 if (Trace
.isOn("cdcac"))
429 Trace("cdcac") << "Returning intervals for "
430 << d_variableOrdering
[curVariable
] << ":" << std::endl
;
431 for (const auto& i
: intervals
)
433 Trace("cdcac") << "-> " << i
.d_interval
<< std::endl
;
439 bool CDCAC::checkIntegrality(std::size_t cur_variable
, const poly::Value
& value
)
441 Node var
= d_constraints
.varMapper()(d_variableOrdering
[cur_variable
]);
442 if (var
.getType() != NodeManager::currentNM()->integerType())
444 // variable is not integral
447 return poly::represents_integer(value
);
450 CACInterval
CDCAC::buildIntegralityInterval(std::size_t cur_variable
,
451 const poly::Value
& value
)
453 poly::Variable var
= d_variableOrdering
[cur_variable
];
454 poly::Integer below
= poly::floor(value
);
455 poly::Integer above
= poly::ceil(value
);
456 // construct var \in (below, above)
457 return CACInterval
{poly::Interval(below
, above
),
460 {var
- below
, var
- above
},
465 bool CDCAC::hasRootAbove(const poly::Polynomial
& p
,
466 const poly::Value
& val
) const
468 auto roots
= poly::isolate_real_roots(p
, d_assignment
);
469 return std::any_of(roots
.begin(), roots
.end(), [&val
](const poly::Value
& r
) {
474 bool CDCAC::hasRootBelow(const poly::Polynomial
& p
,
475 const poly::Value
& val
) const
477 auto roots
= poly::isolate_real_roots(p
, d_assignment
);
478 return std::any_of(roots
.begin(), roots
.end(), [&val
](const poly::Value
& r
) {
486 } // namespace theory