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-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.\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 "options/arith_options.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
  27 namespace std {
  28 /** Generic streaming operator for std::vector. */
  29 template <typename T>
  30 std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
  31 {
  32   CVC5::container_to_stream(os, v);
  33   return os;
  34 }
  35 }  // namespace std
  36
  37 namespace CVC5 {
  38 namespace theory {
  39 namespace arith {
  40 namespace nl {
  41 namespace cad {
  42
  43 CDCAC::CDCAC(context::Context* ctx,
  44              ProofNodeManager* pnm,
  45              const std::vector<poly::Variable>& ordering)
  46     : d_variableOrdering(ordering)
  47 {
  48   if (pnm != nullptr)
  49   {
  50     d_proof.reset(new CADProofGenerator(ctx, pnm));
  51   }
  52 }
  53
  54 void CDCAC::reset()
  55 {
  56   d_constraints.reset();
  57   d_assignment.clear();
  58 }
  59
  60 void CDCAC::computeVariableOrdering()
  61 {
  62   // Actually compute the variable ordering
  63   d_variableOrdering = d_varOrder(d_constraints.getConstraints(),
  64                                   VariableOrderingStrategy::BROWN);
  65   Trace("cdcac") << "Variable ordering is now " << d_variableOrdering
  66                  << std::endl;
  67
  68   // Write variable ordering back to libpoly.
  69   lp_variable_order_t* vo = poly::Context::get_context().get_variable_order();
  70   lp_variable_order_clear(vo);
  71   for (const auto& v : d_variableOrdering)
  72   {
  73     lp_variable_order_push(vo, v.get_internal());
  74   }
  75 }
  76
  77 void CDCAC::retrieveInitialAssignment(NlModel& model, const Node& ran_variable)
  78 {
  79   if (!options::nlCadUseInitial()) return;
  80   d_initialAssignment.clear();
  81   Trace("cdcac") << "Retrieving initial assignment:" << std::endl;
  82   for (const auto& var : d_variableOrdering)
  83   {
  84     Node v = getConstraints().varMapper()(var);
  85     Node val = model.computeConcreteModelValue(v);
  86     poly::Value value = node_to_value(val, ran_variable);
  87     Trace("cdcac") << "\t" << var << " = " << value << std::endl;
  88     d_initialAssignment.emplace_back(value);
  89   }
  90 }
  91 Constraints& CDCAC::getConstraints() { return d_constraints; }
  92 const Constraints& CDCAC::getConstraints() const { return d_constraints; }
  93
  94 const poly::Assignment& CDCAC::getModel() const { return d_assignment; }
  95
  96 const std::vector<poly::Variable>& CDCAC::getVariableOrdering() const
  97 {
  98   return d_variableOrdering;
  99 }
 100
 101 std::vector<CACInterval> CDCAC::getUnsatIntervals(std::size_t cur_variable)
 102 {
 103   std::vector<CACInterval> res;
 104   for (const auto& c : d_constraints.getConstraints())
 105   {
 106     const poly::Polynomial& p = std::get<0>(c);
 107     poly::SignCondition sc = std::get<1>(c);
 108     const Node& n = std::get<2>(c);
 109
 110     if (main_variable(p) != d_variableOrdering[cur_variable])
 111     {
 112       // Constraint is in another variable, ignore it.
 113       continue;
 114     }
 115
 116     Trace("cdcac") << "Infeasible intervals for " << p << " " << sc
 117                    << " 0 over " << d_assignment << std::endl;
 118     auto intervals = infeasible_regions(p, d_assignment, sc);
 119     for (const auto& i : intervals)
 120     {
 121       Trace("cdcac") << "-> " << i << std::endl;
 122       PolyVector l, u, m, d;
 123       m.add(p);
 124       m.pushDownPolys(d, d_variableOrdering[cur_variable]);
 125       if (!is_minus_infinity(get_lower(i))) l = m;
 126       if (!is_plus_infinity(get_upper(i))) u = m;
 127       res.emplace_back(CACInterval{i, l, u, m, d, {n}});
 128       if (isProofEnabled())
 129       {
 130         d_proof->addDirect(
 131             d_constraints.varMapper()(d_variableOrdering[cur_variable]),
 132             d_constraints.varMapper(),
 133             p,
 134             d_assignment,
 135             sc,
 136             i,
 137             n);
 138       }
 139     }
 140   }
 141   pruneRedundantIntervals(res);
 142   return res;
 143 }
 144
 145 bool CDCAC::sampleOutsideWithInitial(const std::vector<CACInterval>& infeasible,
 146                                      poly::Value& sample,
 147                                      std::size_t cur_variable)
 148 {
 149   if (options::nlCadUseInitial() && cur_variable < d_initialAssignment.size())
 150   {
 151     const poly::Value& suggested = d_initialAssignment[cur_variable];
 152     for (const auto& i : infeasible)
 153     {
 154       if (poly::contains(i.d_interval, suggested))
 155       {
 156         d_initialAssignment.clear();
 157         return sampleOutside(infeasible, sample);
 158       }
 159     }
 160     Trace("cdcac") << "Using suggested initial value" << std::endl;
 161     sample = suggested;
 162     return true;
 163   }
 164   return sampleOutside(infeasible, sample);
 165 }
 166
 167 PolyVector CDCAC::requiredCoefficients(const poly::Polynomial& p) const
 168 {
 169   PolyVector res;
 170   for (long deg = degree(p); deg >= 0; --deg)
 171   {
 172     auto coeff = coefficient(p, deg);
 173     if (lp_polynomial_is_constant(coeff.get_internal())) break;
 174     res.add(coeff);
 175     if (evaluate_constraint(coeff, d_assignment, poly::SignCondition::NE))
 176     {
 177       break;
 178     }
 179   }
 180   return res;
 181 }
 182
 183 PolyVector CDCAC::constructCharacterization(std::vector<CACInterval>& intervals)
 184 {
 185   Assert(!intervals.empty()) << "A covering can not be empty";
 186   Trace("cdcac") << "Constructing characterization now" << std::endl;
 187   PolyVector res;
 188
 189   for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
 190   {
 191     cad::makeFinestSquareFreeBasis(intervals[i], intervals[i + 1]);
 192   }
 193
 194   for (const auto& i : intervals)
 195   {
 196     Trace("cdcac") << "Considering " << i.d_interval << std::endl;
 197     Trace("cdcac") << "-> " << i.d_lowerPolys << " / " << i.d_upperPolys
 198                    << " and " << i.d_mainPolys << " / " << i.d_downPolys
 199                    << std::endl;
 200     Trace("cdcac") << "-> " << i.d_origins << std::endl;
 201     for (const auto& p : i.d_downPolys)
 202     {
 203       // Add all polynomial from lower levels.
 204       res.add(p);
 205     }
 206     for (const auto& p : i.d_mainPolys)
 207     {
 208       Trace("cdcac") << "Discriminant of " << p << " -> " << discriminant(p)
 209                      << std::endl;
 210       // Add all discriminants
 211       res.add(discriminant(p));
 212
 213       for (const auto& q : requiredCoefficients(p))
 214       {
 215         // Add all required coefficients
 216         Trace("cdcac") << "Coeff of " << p << " -> " << q << std::endl;
 217         res.add(q);
 218       }
 219       for (const auto& q : i.d_lowerPolys)
 220       {
 221         if (p == q) continue;
 222         // Check whether p(s \times a) = 0 for some a <= l
 223         if (!hasRootBelow(q, get_lower(i.d_interval))) continue;
 224         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
 225                        << resultant(p, q) << std::endl;
 226         res.add(resultant(p, q));
 227       }
 228       for (const auto& q : i.d_upperPolys)
 229       {
 230         if (p == q) continue;
 231         // Check whether p(s \times a) = 0 for some a >= u
 232         if (!hasRootAbove(q, get_upper(i.d_interval))) continue;
 233         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
 234                        << resultant(p, q) << std::endl;
 235         res.add(resultant(p, q));
 236       }
 237     }
 238   }
 239
 240   for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
 241   {
 242     // Add resultants of consecutive intervals.
 243     for (const auto& p : intervals[i].d_upperPolys)
 244     {
 245       for (const auto& q : intervals[i + 1].d_lowerPolys)
 246       {
 247         Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
 248                        << resultant(p, q) << std::endl;
 249         res.add(resultant(p, q));
 250       }
 251     }
 252   }
 253
 254   res.reduce();
 255   res.makeFinestSquareFreeBasis();
 256
 257   return res;
 258 }
 259
 260 CACInterval CDCAC::intervalFromCharacterization(
 261     const PolyVector& characterization,
 262     std::size_t cur_variable,
 263     const poly::Value& sample)
 264 {
 265   PolyVector l;
 266   PolyVector u;
 267   PolyVector m;
 268   PolyVector d;
 269
 270   for (const auto& p : characterization)
 271   {
 272     // Add polynomials to main
 273     m.add(p);
 274   }
 275   // Push lower-dimensional polys to down
 276   m.pushDownPolys(d, d_variableOrdering[cur_variable]);
 277
 278   // Collect -oo, all roots, oo
 279   std::vector<poly::Value> roots;
 280   roots.emplace_back(poly::Value::minus_infty());
 281   for (const auto& p : m)
 282   {
 283     auto tmp = isolate_real_roots(p, d_assignment);
 284     roots.insert(roots.end(), tmp.begin(), tmp.end());
 285   }
 286   roots.emplace_back(poly::Value::plus_infty());
 287   std::sort(roots.begin(), roots.end());
 288
 289   // Now find the interval bounds
 290   poly::Value lower;
 291   poly::Value upper;
 292   for (std::size_t i = 0, n = roots.size(); i < n; ++i)
 293   {
 294     if (sample < roots[i])
 295     {
 296       lower = roots[i - 1];
 297       upper = roots[i];
 298       break;
 299     }
 300     if (roots[i] == sample)
 301     {
 302       lower = sample;
 303       upper = sample;
 304       break;
 305     }
 306   }
 307   Assert(!is_none(lower) && !is_none(upper));
 308
 309   if (!is_minus_infinity(lower))
 310   {
 311     // Identify polynomials that have a root at the lower bound
 312     d_assignment.set(d_variableOrdering[cur_variable], lower);
 313     for (const auto& p : m)
 314     {
 315       if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
 316       {
 317         l.add(p, true);
 318       }
 319     }
 320     d_assignment.unset(d_variableOrdering[cur_variable]);
 321   }
 322   if (!is_plus_infinity(upper))
 323   {
 324     // Identify polynomials that have a root at the upper bound
 325     d_assignment.set(d_variableOrdering[cur_variable], upper);
 326     for (const auto& p : m)
 327     {
 328       if (evaluate_constraint(p, d_assignment, poly::SignCondition::EQ))
 329       {
 330         u.add(p, true);
 331       }
 332     }
 333     d_assignment.unset(d_variableOrdering[cur_variable]);
 334   }
 335
 336   if (lower == upper)
 337   {
 338     // construct a point interval
 339     return CACInterval{
 340         poly::Interval(lower, false, upper, false), l, u, m, d, {}};
 341   }
 342   else
 343   {
 344     // construct an open interval
 345     Assert(lower < upper);
 346     return CACInterval{
 347         poly::Interval(lower, true, upper, true), l, u, m, d, {}};
 348   }
 349 }
 350
 351 std::vector<CACInterval> CDCAC::getUnsatCover(std::size_t curVariable,
 352                                               bool returnFirstInterval)
 353 {
 354   if (isProofEnabled())
 355   {
 356     d_proof->startRecursive();
 357   }
 358   Trace("cdcac") << "Looking for unsat cover for "
 359                  << d_variableOrdering[curVariable] << std::endl;
 360   std::vector<CACInterval> intervals = getUnsatIntervals(curVariable);
 361
 362   if (Trace.isOn("cdcac"))
 363   {
 364     Trace("cdcac") << "Unsat intervals for " << d_variableOrdering[curVariable]
 365                    << ":" << std::endl;
 366     for (const auto& i : intervals)
 367     {
 368       Trace("cdcac") << "-> " << i.d_interval << " from " << i.d_origins
 369                      << std::endl;
 370     }
 371   }
 372   poly::Value sample;
 373
 374   while (sampleOutsideWithInitial(intervals, sample, curVariable))
 375   {
 376     if (!checkIntegrality(curVariable, sample))
 377     {
 378       // the variable is integral, but the sample is not.
 379       Trace("cdcac") << "Used " << sample << " for integer variable "
 380                      << d_variableOrdering[curVariable] << std::endl;
 381       auto newInterval = buildIntegralityInterval(curVariable, sample);
 382       Trace("cdcac") << "Adding integrality interval " << newInterval.d_interval
 383                      << std::endl;
 384       intervals.emplace_back(newInterval);
 385       pruneRedundantIntervals(intervals);
 386       continue;
 387     }
 388     d_assignment.set(d_variableOrdering[curVariable], sample);
 389     Trace("cdcac") << "Sample: " << d_assignment << std::endl;
 390     if (curVariable == d_variableOrdering.size() - 1)
 391     {
 392       // We have a full assignment. SAT!
 393       Trace("cdcac") << "Found full assignment: " << d_assignment << std::endl;
 394       return {};
 395     }
 396     if (isProofEnabled())
 397     {
 398       d_proof->startScope();
 399     }
 400     // Recurse to next variable
 401     auto cov = getUnsatCover(curVariable + 1);
 402     if (cov.empty())
 403     {
 404       // Found SAT!
 405       Trace("cdcac") << "SAT!" << std::endl;
 406       return {};
 407     }
 408     Trace("cdcac") << "Refuting Sample: " << d_assignment << std::endl;
 409     auto characterization = constructCharacterization(cov);
 410     Trace("cdcac") << "Characterization: " << characterization << std::endl;
 411
 412     d_assignment.unset(d_variableOrdering[curVariable]);
 413
 414     auto newInterval =
 415         intervalFromCharacterization(characterization, curVariable, sample);
 416     newInterval.d_origins = collectConstraints(cov);
 417     intervals.emplace_back(newInterval);
 418     if (isProofEnabled())
 419     {
 420       auto cell = d_proof->constructCell(
 421           d_constraints.varMapper()(d_variableOrdering[curVariable]),
 422           newInterval,
 423           d_assignment,
 424           sample,
 425           d_constraints.varMapper());
 426       d_proof->endScope(cell);
 427     }
 428
 429     if (returnFirstInterval)
 430     {
 431       return intervals;
 432     }
 433
 434     Trace("cdcac") << "Added " << intervals.back().d_interval << std::endl;
 435     Trace("cdcac") << "\tlower:   " << intervals.back().d_lowerPolys
 436                    << std::endl;
 437     Trace("cdcac") << "\tupper:   " << intervals.back().d_upperPolys
 438                    << std::endl;
 439     Trace("cdcac") << "\tmain:    " << intervals.back().d_mainPolys
 440                    << std::endl;
 441     Trace("cdcac") << "\tdown:    " << intervals.back().d_downPolys
 442                    << std::endl;
 443     Trace("cdcac") << "\torigins: " << intervals.back().d_origins << std::endl;
 444
 445     // Remove redundant intervals
 446     pruneRedundantIntervals(intervals);
 447   }
 448
 449   if (Trace.isOn("cdcac"))
 450   {
 451     Trace("cdcac") << "Returning intervals for "
 452                    << d_variableOrdering[curVariable] << ":" << std::endl;
 453     for (const auto& i : intervals)
 454     {
 455       Trace("cdcac") << "-> " << i.d_interval << std::endl;
 456     }
 457   }
 458   if (isProofEnabled())
 459   {
 460     d_proof->endRecursive();
 461   }
 462   return intervals;
 463 }
 464
 465 void CDCAC::startNewProof()
 466 {
 467   if (isProofEnabled())
 468   {
 469     d_proof->startNewProof();
 470   }
 471 }
 472
 473 ProofGenerator* CDCAC::closeProof(const std::vector<Node>& assertions)
 474 {
 475   if (isProofEnabled())
 476   {
 477     d_proof->endScope(assertions);
 478     return d_proof->getProofGenerator();
 479   }
 480   return nullptr;
 481 }
 482
 483 bool CDCAC::checkIntegrality(std::size_t cur_variable, const poly::Value& value)
 484 {
 485   Node var = d_constraints.varMapper()(d_variableOrdering[cur_variable]);
 486   if (var.getType() != NodeManager::currentNM()->integerType())
 487   {
 488     // variable is not integral
 489     return true;
 490   }
 491   return poly::represents_integer(value);
 492 }
 493
 494 CACInterval CDCAC::buildIntegralityInterval(std::size_t cur_variable,
 495                                             const poly::Value& value)
 496 {
 497   poly::Variable var = d_variableOrdering[cur_variable];
 498   poly::Integer below = poly::floor(value);
 499   poly::Integer above = poly::ceil(value);
 500   // construct var \in (below, above)
 501   return CACInterval{poly::Interval(below, above),
 502                      {var - below},
 503                      {var - above},
 504                      {var - below, var - above},
 505                      {},
 506                      {}};
 507 }
 508
 509 bool CDCAC::hasRootAbove(const poly::Polynomial& p,
 510                          const poly::Value& val) const
 511 {
 512   auto roots = poly::isolate_real_roots(p, d_assignment);
 513   return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
 514     return r >= val;
 515   });
 516 }
 517
 518 bool CDCAC::hasRootBelow(const poly::Polynomial& p,
 519                          const poly::Value& val) const
 520 {
 521   auto roots = poly::isolate_real_roots(p, d_assignment);
 522   return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
 523     return r <= val;
 524   });
 525 }
 526
 527 void CDCAC::pruneRedundantIntervals(std::vector<CACInterval>& intervals)
 528 {
 529   if (isProofEnabled())
 530   {
 531     std::vector<CACInterval> allIntervals = intervals;
 532     cleanIntervals(intervals);
 533     d_proof->pruneChildren([&allIntervals, &intervals](std::size_t i) {
 534       return std::find(intervals.begin(), intervals.end(), allIntervals[i])
 535              != intervals.end();
 536     });
 537   }
 538   else
 539   {
 540     cleanIntervals(intervals);
 541   }
 542 }
 543
 544 }  // namespace cad
 545 }  // namespace nl
 546 }  // namespace arith
 547 }  // namespace theory
 548 }  // namespace CVC5
 549
 550 #endif