src/theory/arith/nl/cad/cdcac.cpp

   1 /*********************                                                        */
   2 /*! \file cdcac.cpp
   3  ** \verbatim
   4  ** Top contributors (to current version):
   5  **   Gereon Kremer
   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
  11  **
  12  ** \brief Implements the CDCAC approach.
  13  **
  14  ** Implements the CDCAC approach as described in
  15  ** https://arxiv.org/pdf/2003.05633.pdf.
  16  **/
  17
  18 #include "theory/arith/nl/cad/cdcac.h"
  19
  20 #ifdef CVC4_POLY_IMP
  21
  22 #include "theory/arith/nl/cad/projections.h"
  23 #include "theory/arith/nl/cad/variable_ordering.h"
  24
  25 namespace std {
  26 /** Generic streaming operator for std::vector. */
  27 template <typename T>
  28 std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
  29 {
  30   CVC4::container_to_stream(os, v);
  31   return os;
  32 }
  33 }  // namespace std
  34
  35 namespace CVC4 {
  36 namespace theory {
  37 namespace arith {
  38 namespace nl {
  39 namespace cad {
  40
  41 namespace {
  42 /** Removed duplicates from a vector. */
  43 template <typename T>
  44 void removeDuplicates(std::vector<T>& v)
  45 {
  46   std::sort(v.begin(), v.end());
  47   v.erase(std::unique(v.begin(), v.end()), v.end());
  48 }
  49 }  // namespace
  50
  51 CDCAC::CDCAC() {}
  52
  53 CDCAC::CDCAC(const std::vector<poly::Variable>& ordering)
  54     : d_variableOrdering(ordering)
  55 {
  56 }
  57
  58 void CDCAC::reset()
  59 {
  60   d_constraints.reset();
  61   d_assignment.clear();
  62 }
  63
  64 void CDCAC::computeVariableOrdering()
  65 {
  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
  70                  << std::endl;
  71
  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)
  76   {
  77     lp_variable_order_push(vo, v.get_internal());
  78   }
  79 }
  80
  81 void CDCAC::retrieveInitialAssignment(NlModel& model, const Node& ran_variable)
  82 {
  83   d_initialAssignment.clear();
  84   Trace("cdcac") << "Retrieving initial assignment:" << std::endl;
  85   for (const auto& var : d_variableOrdering)
  86   {
  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);
  92   }
  93 }
  94 Constraints& CDCAC::getConstraints() { return d_constraints; }
  95 const Constraints& CDCAC::getConstraints() const { return d_constraints; }
  96
  97 const poly::Assignment& CDCAC::getModel() const { return d_assignment; }
  98
  99 const std::vector<poly::Variable>& CDCAC::getVariableOrdering() const
 100 {
 101   return d_variableOrdering;
 102 }
 103
 104 std::vector<CACInterval> CDCAC::getUnsatIntervals(
 105     std::size_t cur_variable) const
 106 {
 107   std::vector<CACInterval> res;
 108   for (const auto& c : d_constraints.getConstraints())
 109   {
 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);
 113
 114     if (main_variable(p) != d_variableOrdering[cur_variable])
 115     {
 116       // Constraint is in another variable, ignore it.
 117       continue;
 118     }
 119
 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)
 124     {
 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);
 129       m.emplace_back(p);
 130       res.emplace_back(CACInterval{i, l, u, m, d, {n}});
 131     }
 132   }
 133   cleanIntervals(res);
 134   return res;
 135 }
 136
 137 bool CDCAC::sampleOutsideWithInitial(const std::vector<CACInterval>& infeasible,
 138                                      poly::Value& sample,
 139                                      std::size_t cur_variable)
 140 {
 141   if (cur_variable < d_initialAssignment.size())
 142   {
 143     const poly::Value& suggested = d_initialAssignment[cur_variable];
 144     for (const auto& i : infeasible)
 145     {
 146       if (poly::contains(i.d_interval, suggested))
 147       {
 148         d_initialAssignment.clear();
 149         return sampleOutside(infeasible, sample);
 150       }
 151     }
 152     Trace("cdcac") << "Using suggested initial value" << std::endl;
 153     sample = suggested;
 154     return true;
 155   }
 156   return sampleOutside(infeasible, sample);
 157 }
 158
 159 std::vector<poly::Polynomial> CDCAC::requiredCoefficients(
 160     const poly::Polynomial& p) const
 161 {
 162   std::vector<poly::Polynomial> res;
 163   for (long deg = degree(p); deg >= 0; --deg)
 164   {
 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))
 169     {
 170       break;
 171     }
 172   }
 173   return res;
 174 }
 175
 176 std::vector<poly::Polynomial> CDCAC::constructCharacterization(
 177     std::vector<CACInterval>& intervals)
 178 {
 179   Assert(!intervals.empty()) << "A covering can not be empty";
 180   Trace("cdcac") << "Constructing characterization now" << std::endl;
 181   std::vector<poly::Polynomial> res;
 182
 183   for (const auto& i : intervals)
 184   {
 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
 188                    << std::endl;
 189     Trace("cdcac") << "-> " << i.d_origins << std::endl;
 190     for (const auto& p : i.d_downPolys)
 191     {
 192       // Add all polynomial from lower levels.
 193       addPolynomial(res, p);
 194     }
 195     for (const auto& p : i.d_mainPolys)
 196     {
 197       Trace("cdcac") << "Discriminant of " << p << " -> " << discriminant(p)
 198                      << std::endl;
 199       // Add all discriminants
 200       addPolynomial(res, discriminant(p));
 201
 202       for (const auto& q : requiredCoefficients(p))
 203       {
 204         // Add all required coefficients
 205         Trace("cdcac") << "Coeff of " << p << " -> " << q << std::endl;
 206         addPolynomial(res, q);
 207       }
 208       for (const auto& q : i.d_lowerPolys)
 209       {
 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));
 216       }
 217       for (const auto& q : i.d_upperPolys)
 218       {
 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));
 225       }
 226     }
 227   }
 228
 229   for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
 230   {
 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)
 235     {
 236       for (const auto& q : intervals[i + 1].d_lowerPolys)
 237       {
 238         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
 239                        << resultant(p, q) << std::endl;
 240         addPolynomial(res, resultant(p, q));
 241       }
 242     }
 243   }
 244
 245   removeDuplicates(res);
 246   makeFinestSquareFreeBasis(res);
 247
 248   return res;
 249 }
 250
 251 CACInterval CDCAC::intervalFromCharacterization(
 252     const std::vector<poly::Polynomial>& characterization,
 253     std::size_t cur_variable,
 254     const poly::Value& sample)
 255 {
 256   std::vector<poly::Polynomial> l;
 257   std::vector<poly::Polynomial> u;
 258   std::vector<poly::Polynomial> m;
 259   std::vector<poly::Polynomial> d;
 260
 261   for (const auto& p : characterization)
 262   {
 263     // Add polynomials to either main or down
 264     if (main_variable(p) == d_variableOrdering[cur_variable])
 265     {
 266       m.emplace_back(p);
 267     }
 268     else
 269     {
 270       d.emplace_back(p);
 271     }
 272   }
 273
 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)
 278   {
 279     auto tmp = isolate_real_roots(p, d_assignment);
 280     roots.insert(roots.end(), tmp.begin(), tmp.end());
 281   }
 282   roots.emplace_back(poly::Value::plus_infty());
 283   std::sort(roots.begin(), roots.end());
 284
 285   // Now find the interval bounds
 286   poly::Value lower;
 287   poly::Value upper;
 288   for (std::size_t i = 0, n = roots.size(); i < n; ++i)
 289   {
 290     if (sample < roots[i])
 291     {
 292       lower = roots[i - 1];
 293       upper = roots[i];
 294       break;
 295     }
 296     if (roots[i] == sample)
 297     {
 298       lower = sample;
 299       upper = sample;
 300       break;
 301     }
 302   }
 303   Assert(!is_none(lower) && !is_none(upper));
 304
 305   if (!is_minus_infinity(lower))
 306   {
 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)
 310     {
 311       if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
 312       {
 313         l.emplace_back(p);
 314       }
 315     }
 316     d_assignment.unset(d_variableOrdering[cur_variable]);
 317   }
 318   if (!is_plus_infinity(upper))
 319   {
 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)
 323     {
 324       if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
 325       {
 326         u.emplace_back(p);
 327       }
 328     }
 329     d_assignment.unset(d_variableOrdering[cur_variable]);
 330   }
 331
 332   if (lower == upper)
 333   {
 334     // construct a point interval
 335     return CACInterval{
 336         poly::Interval(lower, false, upper, false), l, u, m, d, {}};
 337   }
 338   else
 339   {
 340     // construct an open interval
 341     Assert(lower < upper);
 342     return CACInterval{
 343         poly::Interval(lower, true, upper, true), l, u, m, d, {}};
 344   }
 345 }
 346
 347 std::vector<CACInterval> CDCAC::getUnsatCover(std::size_t curVariable,
 348                                               bool returnFirstInterval)
 349 {
 350   Trace("cdcac") << "Looking for unsat cover for "
 351                  << d_variableOrdering[curVariable] << std::endl;
 352   std::vector<CACInterval> intervals = getUnsatIntervals(curVariable);
 353
 354   if (Trace.isOn("cdcac"))
 355   {
 356     Trace("cdcac") << "Unsat intervals for " << d_variableOrdering[curVariable]
 357                    << ":" << std::endl;
 358     for (const auto& i : intervals)
 359     {
 360       Trace("cdcac") << "-> " << i.d_interval << " from " << i.d_origins
 361                      << std::endl;
 362     }
 363   }
 364   poly::Value sample;
 365
 366   while (sampleOutsideWithInitial(intervals, sample, curVariable))
 367   {
 368     if (!checkIntegrality(curVariable, sample))
 369     {
 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
 375                      << std::endl;
 376       intervals.emplace_back(newInterval);
 377       cleanIntervals(intervals);
 378       continue;
 379     }
 380     d_assignment.set(d_variableOrdering[curVariable], sample);
 381     Trace("cdcac") << "Sample: " << d_assignment << std::endl;
 382     if (curVariable == d_variableOrdering.size() - 1)
 383     {
 384       // We have a full assignment. SAT!
 385       Trace("cdcac") << "Found full assignment: " << d_assignment << std::endl;
 386       return {};
 387     }
 388     // Recurse to next variable
 389     auto cov = getUnsatCover(curVariable + 1);
 390     if (cov.empty())
 391     {
 392       // Found SAT!
 393       Trace("cdcac") << "SAT!" << std::endl;
 394       return {};
 395     }
 396     Trace("cdcac") << "Refuting Sample: " << d_assignment << std::endl;
 397     auto characterization = constructCharacterization(cov);
 398     Trace("cdcac") << "Characterization: " << characterization << std::endl;
 399
 400     d_assignment.unset(d_variableOrdering[curVariable]);
 401
 402     auto newInterval =
 403         intervalFromCharacterization(characterization, curVariable, sample);
 404     newInterval.d_origins = collectConstraints(cov);
 405     intervals.emplace_back(newInterval);
 406
 407     if (returnFirstInterval)
 408     {
 409       return intervals;
 410     }
 411
 412     Trace("cdcac") << "Added " << intervals.back().d_interval << std::endl;
 413     Trace("cdcac") << "\tlower:   " << intervals.back().d_lowerPolys
 414                    << std::endl;
 415     Trace("cdcac") << "\tupper:   " << intervals.back().d_upperPolys
 416                    << std::endl;
 417     Trace("cdcac") << "\tmain:    " << intervals.back().d_mainPolys
 418                    << std::endl;
 419     Trace("cdcac") << "\tdown:    " << intervals.back().d_downPolys
 420                    << std::endl;
 421     Trace("cdcac") << "\torigins: " << intervals.back().d_origins << std::endl;
 422
 423     // Remove redundant intervals
 424     cleanIntervals(intervals);
 425   }
 426
 427   if (Trace.isOn("cdcac"))
 428   {
 429     Trace("cdcac") << "Returning intervals for "
 430                    << d_variableOrdering[curVariable] << ":" << std::endl;
 431     for (const auto& i : intervals)
 432     {
 433       Trace("cdcac") << "-> " << i.d_interval << std::endl;
 434     }
 435   }
 436   return intervals;
 437 }
 438
 439 bool CDCAC::checkIntegrality(std::size_t cur_variable, const poly::Value& value)
 440 {
 441   Node var = d_constraints.varMapper()(d_variableOrdering[cur_variable]);
 442   if (var.getType() != NodeManager::currentNM()->integerType())
 443   {
 444     // variable is not integral
 445     return true;
 446   }
 447   return poly::represents_integer(value);
 448 }
 449
 450 CACInterval CDCAC::buildIntegralityInterval(std::size_t cur_variable,
 451                                             const poly::Value& value)
 452 {
 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),
 458                      {var - below},
 459                      {var - above},
 460                      {var - below, var - above},
 461                      {},
 462                      {}};
 463 }
 464
 465 bool CDCAC::hasRootAbove(const poly::Polynomial& p,
 466                          const poly::Value& val) const
 467 {
 468   auto roots = poly::isolate_real_roots(p, d_assignment);
 469   return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
 470     return r >= val;
 471   });
 472 }
 473
 474 bool CDCAC::hasRootBelow(const poly::Polynomial& p,
 475                          const poly::Value& val) const
 476 {
 477   auto roots = poly::isolate_real_roots(p, d_assignment);
 478   return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
 479     return r <= val;
 480   });
 481 }
 482
 483 }  // namespace cad
 484 }  // namespace nl
 485 }  // namespace arith
 486 }  // namespace theory
 487 }  // namespace CVC4
 488
 489 #endif