Cad implementation (#4774)

This commit implements the CAD interface added in #4773.

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 <typename T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T>& 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 <typename T>
+void removeDuplicates(std::vector<T>& 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<poly::Variable>& 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<poly::Variable>& CDCAC::getVariableOrdering() const
 std::vector<CACInterval> CDCAC::getUnsatIntervals(
     std::size_t cur_variable) const
 {
-  return {};
+  std::vector<CACInterval> 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<poly::Polynomial> 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<poly::Polynomial> CDCAC::requiredCoefficients(
     const poly::Polynomial& p) const
 {
-  return {};
+  std::vector<poly::Polynomial> 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<poly::Polynomial> CDCAC::constructCharacterization(
     std::vector<CACInterval>& intervals)
 {
-  return {};
+  Assert(!intervals.empty()) << "A covering can not be empty";
+  Trace("cdcac") << "Constructing characterization now" << std::endl;
+  std::vector<poly::Polynomial> 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<poly::Polynomial> l;
+  std::vector<poly::Polynomial> u;
+  std::vector<poly::Polynomial> m;
+  std::vector<poly::Polynomial> 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<poly::Value> 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<CACInterval> CDCAC::getUnsatCover(std::size_t cur_variable)
 {
-  return {};
+  Trace("cdcac") << "Looking for unsat cover for "
+                 << d_variableOrdering[cur_variable] << std::endl;
+  std::vector<CACInterval> 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<CACInterval> 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 <poly/polyxx.h>
 
 #include <vector>
@@ -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<CACInterval>& 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 <poly/polyxx.h>
 
 #include <vector>
@@ -98,3 +102,5 @@ bool sampleOutside(const std::vector<CACInterval>& infeasible,
 }  // namespace CVC4
 
 #endif
+
+#endif