From 3e18cd977b6e8d9729a4aa0f4cfc12710d21c863 Mon Sep 17 00:00:00 2001 From: Gereon Kremer Date: Thu, 30 Jul 2020 23:59:40 +0200 Subject: [PATCH] Cad implementation (#4774) This commit implements the CAD interface added in #4773. --- src/theory/arith/nl/cad/cdcac.cpp | 313 +++++++++++++++++++++++- src/theory/arith/nl/cad/cdcac.h | 6 + src/theory/arith/nl/cad/cdcac_utils.cpp | 6 +- src/theory/arith/nl/cad/cdcac_utils.h | 6 + 4 files changed, 324 insertions(+), 7 deletions(-) diff --git a/src/theory/arith/nl/cad/cdcac.cpp b/src/theory/arith/nl/cad/cdcac.cpp index 18b16bbe8..dbdae3e0b 100644 --- a/src/theory/arith/nl/cad/cdcac.cpp +++ b/src/theory/arith/nl/cad/cdcac.cpp @@ -17,15 +17,37 @@ #include "theory/arith/nl/cad/cdcac.h" +#ifdef CVC4_POLY_IMP + #include "theory/arith/nl/cad/projections.h" #include "theory/arith/nl/cad/variable_ordering.h" +namespace std { +/** Generic streaming operator for std::vector. */ +template +std::ostream& operator<<(std::ostream& os, const std::vector& v) +{ + CVC4::container_to_stream(os, v); + return os; +} +} // namespace std + namespace CVC4 { namespace theory { namespace arith { namespace nl { namespace cad { +namespace { +/** Removed duplicates from a vector. */ +template +void removeDuplicates(std::vector& v) +{ + std::sort(v.begin(), v.end()); + v.erase(std::unique(v.begin(), v.end()), v.end()); +} +} // namespace + CDCAC::CDCAC() {} CDCAC::CDCAC(const std::vector& ordering) @@ -39,7 +61,22 @@ void CDCAC::reset() d_assignment.clear(); } -void CDCAC::computeVariableOrdering() {} +void CDCAC::computeVariableOrdering() +{ + // Actually compute the variable ordering + d_variableOrdering = d_varOrder(d_constraints.getConstraints(), + VariableOrderingStrategy::BROWN); + Trace("cdcac") << "Variable ordering is now " << d_variableOrdering + << std::endl; + + // Write variable ordering back to libpoly. + lp_variable_order_t* vo = poly::Context::get_context().get_variable_order(); + lp_variable_order_clear(vo); + for (const auto& v : d_variableOrdering) + { + lp_variable_order_push(vo, v.get_internal()); + } +} Constraints& CDCAC::getConstraints() { return d_constraints; } const Constraints& CDCAC::getConstraints() const { return d_constraints; } @@ -54,19 +91,124 @@ const std::vector& CDCAC::getVariableOrdering() const std::vector CDCAC::getUnsatIntervals( std::size_t cur_variable) const { - return {}; + std::vector res; + for (const auto& c : d_constraints.getConstraints()) + { + const poly::Polynomial& p = std::get<0>(c); + poly::SignCondition sc = std::get<1>(c); + const Node& n = std::get<2>(c); + + if (main_variable(p) != d_variableOrdering[cur_variable]) + { + // Constraint is in another variable, ignore it. + continue; + } + + Trace("cdcac") << "Infeasible intervals for " << p << " " << sc + << " 0 over " << d_assignment << std::endl; + auto intervals = infeasible_regions(p, d_assignment, sc); + for (const auto& i : intervals) + { + Trace("cdcac") << "-> " << i << std::endl; + std::vector l, u, m, d; + if (!is_minus_infinity(get_lower(i))) l.emplace_back(p); + if (!is_plus_infinity(get_upper(i))) u.emplace_back(p); + m.emplace_back(p); + res.emplace_back(CACInterval{i, l, u, m, d, {n}}); + } + } + cleanIntervals(res); + return res; } std::vector CDCAC::requiredCoefficients( const poly::Polynomial& p) const { - return {}; + std::vector res; + for (long deg = degree(p); deg >= 0; --deg) + { + auto coeff = coefficient(p, deg); + if (lp_polynomial_is_constant(coeff.get_internal())) break; + res.emplace_back(coeff); + if (evaluate_constraint(coeff, d_assignment, poly::SignCondition::NE)) + { + break; + } + } + return res; } std::vector CDCAC::constructCharacterization( std::vector& intervals) { - return {}; + Assert(!intervals.empty()) << "A covering can not be empty"; + Trace("cdcac") << "Constructing characterization now" << std::endl; + std::vector res; + + for (const auto& i : intervals) + { + Trace("cdcac") << "Considering " << i.d_interval << std::endl; + Trace("cdcac") << "-> " << i.d_lowerPolys << " / " << i.d_upperPolys + << " and " << i.d_mainPolys << " / " << i.d_downPolys + << std::endl; + Trace("cdcac") << "-> " << i.d_origins << std::endl; + for (const auto& p : i.d_downPolys) + { + // Add all polynomial from lower levels. + addPolynomial(res, p); + } + for (const auto& p : i.d_mainPolys) + { + Trace("cdcac") << "Discriminant of " << p << " -> " << discriminant(p) + << std::endl; + // Add all discriminants + addPolynomial(res, discriminant(p)); + + for (const auto& q : requiredCoefficients(p)) + { + // Add all required coefficients + Trace("cdcac") << "Coeff of " << p << " -> " << q << std::endl; + addPolynomial(res, q); + } + // TODO(cvc4-projects #210): Only add if p(s \times a) = 0 for some a <= l + for (const auto& q : i.d_lowerPolys) + { + if (p == q) continue; + Trace("cdcac") << "Resultant of " << p << " and " << q << " -> " + << resultant(p, q) << std::endl; + addPolynomial(res, resultant(p, q)); + } + // TODO(cvc4-projects #210): Only add if p(s \times a) = 0 for some a >= u + for (const auto& q : i.d_upperPolys) + { + if (p == q) continue; + Trace("cdcac") << "Resultant of " << p << " and " << q << " -> " + << resultant(p, q) << std::endl; + addPolynomial(res, resultant(p, q)); + } + } + } + + for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i) + { + // Add resultants of consecutive intervals. + cad::makeFinestSquareFreeBasis(intervals[i].d_upperPolys, + intervals[i + 1].d_lowerPolys); + for (const auto& p : intervals[i].d_upperPolys) + { + for (const auto& q : intervals[i + 1].d_lowerPolys) + { + Trace("cdcac") << "Resultant of " << p << " and " << q << " -> " + << resultant(p, q) << std::endl; + addPolynomial(res, resultant(p, q)); + } + } + } + + removeDuplicates(res); + makeFinestSquareFreeBasis(res); + + return res; } CACInterval CDCAC::intervalFromCharacterization( @@ -74,12 +216,169 @@ CACInterval CDCAC::intervalFromCharacterization( std::size_t cur_variable, const poly::Value& sample) { - return {}; + std::vector l; + std::vector u; + std::vector m; + std::vector d; + + for (const auto& p : characterization) + { + // Add polynomials to either main or down + if (main_variable(p) == d_variableOrdering[cur_variable]) + { + m.emplace_back(p); + } + else + { + d.emplace_back(p); + } + } + + // Collect -oo, all roots, oo + std::vector roots; + roots.emplace_back(poly::Value::minus_infty()); + for (const auto& p : m) + { + auto tmp = isolate_real_roots(p, d_assignment); + roots.insert(roots.end(), tmp.begin(), tmp.end()); + } + roots.emplace_back(poly::Value::plus_infty()); + std::sort(roots.begin(), roots.end()); + + // Now find the interval bounds + poly::Value lower; + poly::Value upper; + for (std::size_t i = 0, n = roots.size(); i < n; ++i) + { + if (sample < roots[i]) + { + lower = roots[i - 1]; + upper = roots[i]; + break; + } + if (roots[i] == sample) + { + lower = sample; + upper = sample; + break; + } + } + Assert(!is_none(lower) && !is_none(upper)); + + if (!is_minus_infinity(lower)) + { + // Identify polynomials that have a root at the lower bound + d_assignment.set(d_variableOrdering[cur_variable], lower); + for (const auto& p : m) + { + if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ)) + { + l.emplace_back(p); + } + } + d_assignment.unset(d_variableOrdering[cur_variable]); + } + if (!is_plus_infinity(upper)) + { + // Identify polynomials that have a root at the upper bound + d_assignment.set(d_variableOrdering[cur_variable], upper); + for (const auto& p : m) + { + if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ)) + { + u.emplace_back(p); + } + } + d_assignment.unset(d_variableOrdering[cur_variable]); + } + + if (lower == upper) + { + // construct a point interval + return CACInterval{ + poly::Interval(lower, false, upper, false), l, u, m, d, {}}; + } + else + { + // construct an open interval + Assert(lower < upper); + return CACInterval{ + poly::Interval(lower, true, upper, true), l, u, m, d, {}}; + } } std::vector CDCAC::getUnsatCover(std::size_t cur_variable) { - return {}; + Trace("cdcac") << "Looking for unsat cover for " + << d_variableOrdering[cur_variable] << std::endl; + std::vector intervals = getUnsatIntervals(cur_variable); + + if (Trace.isOn("cdcac")) + { + Trace("cdcac") << "Unsat intervals for " << d_variableOrdering[cur_variable] + << ":" << std::endl; + for (const auto& i : intervals) + { + Trace("cdcac") << "-> " << i.d_interval << " from " << i.d_origins + << std::endl; + } + } + poly::Value sample; + + while (sampleOutside(intervals, sample)) + { + d_assignment.set(d_variableOrdering[cur_variable], sample); + Trace("cdcac") << "Sample: " << d_assignment << std::endl; + if (cur_variable == d_variableOrdering.size() - 1) + { + // We have a full assignment. SAT! + Trace("cdcac") << "Found full assignment: " << d_assignment << std::endl; + return {}; + } + // Recurse to next variable + auto cov = getUnsatCover(cur_variable + 1); + if (cov.empty()) + { + // Found SAT! + Trace("cdcac") << "SAT!" << std::endl; + return {}; + } + Trace("cdcac") << "Refuting Sample: " << d_assignment << std::endl; + auto characterization = constructCharacterization(cov); + Trace("cdcac") << "Characterization: " << characterization << std::endl; + + d_assignment.unset(d_variableOrdering[cur_variable]); + + auto new_interval = + intervalFromCharacterization(characterization, cur_variable, sample); + new_interval.d_origins = collectConstraints(cov); + intervals.emplace_back(new_interval); + + Trace("cdcac") << "Added " << intervals.back().d_interval << std::endl; + Trace("cdcac") << "\tlower: " << intervals.back().d_lowerPolys + << std::endl; + Trace("cdcac") << "\tupper: " << intervals.back().d_upperPolys + << std::endl; + Trace("cdcac") << "\tmain: " << intervals.back().d_mainPolys + << std::endl; + Trace("cdcac") << "\tdown: " << intervals.back().d_downPolys + << std::endl; + Trace("cdcac") << "\torigins: " << intervals.back().d_origins << std::endl; + + // Remove redundant intervals + cleanIntervals(intervals); + } + + if (Trace.isOn("cdcac")) + { + Trace("cdcac") << "Returning intervals for " + << d_variableOrdering[cur_variable] << ":" << std::endl; + for (const auto& i : intervals) + { + Trace("cdcac") << "-> " << i.d_interval << std::endl; + } + } + return intervals; } } // namespace cad @@ -87,3 +386,5 @@ std::vector CDCAC::getUnsatCover(std::size_t cur_variable) } // namespace arith } // namespace theory } // namespace CVC4 + +#endif diff --git a/src/theory/arith/nl/cad/cdcac.h b/src/theory/arith/nl/cad/cdcac.h index 88e260ada..3ea281cec 100644 --- a/src/theory/arith/nl/cad/cdcac.h +++ b/src/theory/arith/nl/cad/cdcac.h @@ -20,6 +20,10 @@ #ifndef CVC4__THEORY__ARITH__NL__CAD__CDCAC_H #define CVC4__THEORY__ARITH__NL__CAD__CDCAC_H +#include "util/real_algebraic_number.h" + +#ifdef CVC4_POLY_IMP + #include #include @@ -136,3 +140,5 @@ class CDCAC } // namespace CVC4 #endif + +#endif diff --git a/src/theory/arith/nl/cad/cdcac_utils.cpp b/src/theory/arith/nl/cad/cdcac_utils.cpp index 02df1a323..a9b4c6128 100644 --- a/src/theory/arith/nl/cad/cdcac_utils.cpp +++ b/src/theory/arith/nl/cad/cdcac_utils.cpp @@ -16,6 +16,8 @@ #include "theory/arith/nl/cad/cdcac_utils.h" +#ifdef CVC4_POLY_IMP + namespace CVC4 { namespace theory { namespace arith { @@ -266,4 +268,6 @@ bool sampleOutside(const std::vector& infeasible, Value& sample) } // namespace nl } // namespace arith } // namespace theory -} // namespace CVC4 \ No newline at end of file +} // namespace CVC4 + +#endif diff --git a/src/theory/arith/nl/cad/cdcac_utils.h b/src/theory/arith/nl/cad/cdcac_utils.h index d5fa70bd9..ed1c0d1c2 100644 --- a/src/theory/arith/nl/cad/cdcac_utils.h +++ b/src/theory/arith/nl/cad/cdcac_utils.h @@ -19,6 +19,10 @@ #ifndef CVC4__THEORY__ARITH__NL__CAD__CDCAC_UTILS_H #define CVC4__THEORY__ARITH__NL__CAD__CDCAC_UTILS_H +#include "util/real_algebraic_number.h" + +#ifdef CVC4_POLY_IMP + #include #include @@ -98,3 +102,5 @@ bool sampleOutside(const std::vector& infeasible, } // namespace CVC4 #endif + +#endif -- 2.30.2