Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	src/theory/arith/nl/cad/cdcac.cpp	Lines:	228	271	84.1 %
Date:	2021-05-24	Branches:	395	924	42.7 %


/******************************************************************************
 * Top contributors (to current version):
 *   Gereon Kremer, Aina Niemetz
 *
 * This file is part of the cvc5 project.
 *
 * Copyright (c) 2009-2021 by the authors listed in the file AUTHORS
 * in the top-level source directory and their institutional affiliations.
 * All rights reserved.  See the file COPYING in the top-level source
 * directory for licensing information.
 * ****************************************************************************
 *
 * Implements the CDCAC approach as described in
 * https://arxiv.org/pdf/2003.05633.pdf.
 */

#include "theory/arith/nl/cad/cdcac.h"

#ifdef CVC5_POLY_IMP

#include "options/arith_options.h"
#include "theory/arith/nl/cad/projections.h"
#include "theory/arith/nl/cad/variable_ordering.h"
#include "theory/arith/nl/nl_model.h"

namespace std {
/** Generic streaming operator for std::vector. */
template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
{
  cvc5::container_to_stream(os, v);
  return os;
}
}  // namespace std

namespace cvc5 {
namespace theory {
namespace arith {
namespace nl {
namespace cad {

CDCAC::CDCAC(context::Context* ctx,
             ProofNodeManager* pnm,
             const std::vector<poly::Variable>& ordering)
    : d_variableOrdering(ordering)
{
  if (pnm != nullptr)
  {
    d_proof.reset(new CADProofGenerator(ctx, pnm));
  }
}

void CDCAC::reset()
{
  d_constraints.reset();
  d_assignment.clear();
}

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());
  }
}

void CDCAC::retrieveInitialAssignment(NlModel& model, const Node& ran_variable)
{
  if (!options::nlCadUseInitial()) return;
  d_initialAssignment.clear();
  Trace("cdcac") << "Retrieving initial assignment:" << std::endl;
  for (const auto& var : d_variableOrdering)
  {
    Node v = getConstraints().varMapper()(var);
    Node val = model.computeConcreteModelValue(v);
    poly::Value value = node_to_value(val, ran_variable);
    Trace("cdcac") << "\t" << var << " = " << value << std::endl;
    d_initialAssignment.emplace_back(value);
  }
}
Constraints& CDCAC::getConstraints() { return d_constraints; }
const Constraints& CDCAC::getConstraints() const { return d_constraints; }

const poly::Assignment& CDCAC::getModel() const { return d_assignment; }

const std::vector<poly::Variable>& CDCAC::getVariableOrdering() const
{
  return d_variableOrdering;
}

std::vector<CACInterval> CDCAC::getUnsatIntervals(std::size_t cur_variable)
{
  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;
      PolyVector l, u, m, d;
      m.add(p);
      m.pushDownPolys(d, d_variableOrdering[cur_variable]);
      if (!is_minus_infinity(get_lower(i))) l = m;
      if (!is_plus_infinity(get_upper(i))) u = m;
      res.emplace_back(CACInterval{i, l, u, m, d, {n}});
      if (isProofEnabled())
      {
        d_proof->addDirect(
            d_constraints.varMapper()(d_variableOrdering[cur_variable]),
            d_constraints.varMapper(),
            p,
            d_assignment,
            sc,
            i,
            n);
      }
    }
  }
  pruneRedundantIntervals(res);
  return res;
}

bool CDCAC::sampleOutsideWithInitial(const std::vector<CACInterval>& infeasible,
                                     poly::Value& sample,
                                     std::size_t cur_variable)
{
  if (options::nlCadUseInitial() && cur_variable < d_initialAssignment.size())
  {
    const poly::Value& suggested = d_initialAssignment[cur_variable];
    for (const auto& i : infeasible)
    {
      if (poly::contains(i.d_interval, suggested))
      {
        d_initialAssignment.clear();
        return sampleOutside(infeasible, sample);
      }
    }
    Trace("cdcac") << "Using suggested initial value" << std::endl;
    sample = suggested;
    return true;
  }
  return sampleOutside(infeasible, sample);
}

PolyVector CDCAC::requiredCoefficients(const poly::Polynomial& p) const
{
  PolyVector res;
  for (long deg = degree(p); deg >= 0; --deg)
  {
    auto coeff = coefficient(p, deg);
    if (lp_polynomial_is_constant(coeff.get_internal())) break;
    res.add(coeff);
    if (evaluate_constraint(coeff, d_assignment, poly::SignCondition::NE))
    {
      break;
    }
  }
  return res;
}

PolyVector CDCAC::constructCharacterization(std::vector<CACInterval>& intervals)
{
  Assert(!intervals.empty()) << "A covering can not be empty";
  Trace("cdcac") << "Constructing characterization now" << std::endl;
  PolyVector res;

  for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
  {
    cad::makeFinestSquareFreeBasis(intervals[i], intervals[i + 1]);
  }

  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.
      res.add(p);
    }
    for (const auto& p : i.d_mainPolys)
    {
      Trace("cdcac") << "Discriminant of " << p << " -> " << discriminant(p)
                     << std::endl;
      // Add all discriminants
      res.add(discriminant(p));

      for (const auto& q : requiredCoefficients(p))
      {
        // Add all required coefficients
        Trace("cdcac") << "Coeff of " << p << " -> " << q << std::endl;
        res.add(q);
      }
      for (const auto& q : i.d_lowerPolys)
      {
        if (p == q) continue;
        // Check whether p(s \times a) = 0 for some a <= l
        if (!hasRootBelow(q, get_lower(i.d_interval))) continue;
        Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
                       << resultant(p, q) << std::endl;
        res.add(resultant(p, q));
      }
      for (const auto& q : i.d_upperPolys)
      {
        if (p == q) continue;
        // Check whether p(s \times a) = 0 for some a >= u
        if (!hasRootAbove(q, get_upper(i.d_interval))) continue;
        Trace("cdcac") << "Resultant of " << p << " and " << q << " -> "
                       << resultant(p, q) << std::endl;
        res.add(resultant(p, q));
      }
    }
  }

  for (std::size_t i = 0, n = intervals.size(); i < n - 1; ++i)
  {
    // Add resultants of consecutive intervals.
    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;
        res.add(resultant(p, q));
      }
    }
  }

  res.reduce();
  res.makeFinestSquareFreeBasis();

  return res;
}

CACInterval CDCAC::intervalFromCharacterization(
    const PolyVector& characterization,
    std::size_t cur_variable,
    const poly::Value& sample)
{
  PolyVector l;
  PolyVector u;
  PolyVector m;
  PolyVector d;

  for (const auto& p : characterization)
  {
    // Add polynomials to main
    m.add(p);
  }
  // Push lower-dimensional polys to down
  m.pushDownPolys(d, d_variableOrdering[cur_variable]);

  // 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.add(p, true);
      }
    }
    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.add(p, true);
      }
    }
    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 curVariable,
                                              bool returnFirstInterval)
{
  if (isProofEnabled())
  {
    d_proof->startRecursive();
  }
  Trace("cdcac") << "Looking for unsat cover for "
                 << d_variableOrdering[curVariable] << std::endl;
  std::vector<CACInterval> intervals = getUnsatIntervals(curVariable);

  if (Trace.isOn("cdcac"))
  {
    Trace("cdcac") << "Unsat intervals for " << d_variableOrdering[curVariable]
                   << ":" << std::endl;
    for (const auto& i : intervals)
    {
      Trace("cdcac") << "-> " << i.d_interval << " from " << i.d_origins
                     << std::endl;
    }
  }
  poly::Value sample;

  while (sampleOutsideWithInitial(intervals, sample, curVariable))
  {
    if (!checkIntegrality(curVariable, sample))
    {
      // the variable is integral, but the sample is not.
      Trace("cdcac") << "Used " << sample << " for integer variable "
                     << d_variableOrdering[curVariable] << std::endl;
      auto newInterval = buildIntegralityInterval(curVariable, sample);
      Trace("cdcac") << "Adding integrality interval " << newInterval.d_interval
                     << std::endl;
      intervals.emplace_back(newInterval);
      pruneRedundantIntervals(intervals);
      continue;
    }
    d_assignment.set(d_variableOrdering[curVariable], sample);
    Trace("cdcac") << "Sample: " << d_assignment << std::endl;
    if (curVariable == d_variableOrdering.size() - 1)
    {
      // We have a full assignment. SAT!
      Trace("cdcac") << "Found full assignment: " << d_assignment << std::endl;
      return {};
    }
    if (isProofEnabled())
    {
      d_proof->startScope();
    }
    // Recurse to next variable
    auto cov = getUnsatCover(curVariable + 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[curVariable]);

    auto newInterval =
        intervalFromCharacterization(characterization, curVariable, sample);
    newInterval.d_origins = collectConstraints(cov);
    intervals.emplace_back(newInterval);
    if (isProofEnabled())
    {
      auto cell = d_proof->constructCell(
          d_constraints.varMapper()(d_variableOrdering[curVariable]),
          newInterval,
          d_assignment,
          sample,
          d_constraints.varMapper());
      d_proof->endScope(cell);
    }

    if (returnFirstInterval)
    {
      return intervals;
    }

    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
    pruneRedundantIntervals(intervals);
  }

  if (Trace.isOn("cdcac"))
  {
    Trace("cdcac") << "Returning intervals for "
                   << d_variableOrdering[curVariable] << ":" << std::endl;
    for (const auto& i : intervals)
    {
      Trace("cdcac") << "-> " << i.d_interval << std::endl;
    }
  }
  if (isProofEnabled())
  {
    d_proof->endRecursive();
  }
  return intervals;
}

void CDCAC::startNewProof()
{
  if (isProofEnabled())
  {
    d_proof->startNewProof();
  }
}

ProofGenerator* CDCAC::closeProof(const std::vector<Node>& assertions)
{
  if (isProofEnabled())
  {
    d_proof->endScope(assertions);
    return d_proof->getProofGenerator();
  }
  return nullptr;
}

bool CDCAC::checkIntegrality(std::size_t cur_variable, const poly::Value& value)
{
  Node var = d_constraints.varMapper()(d_variableOrdering[cur_variable]);
  if (var.getType() != NodeManager::currentNM()->integerType())
  {
    // variable is not integral
    return true;
  }
  return poly::represents_integer(value);
}

CACInterval CDCAC::buildIntegralityInterval(std::size_t cur_variable,
                                            const poly::Value& value)
{
  poly::Variable var = d_variableOrdering[cur_variable];
  poly::Integer below = poly::floor(value);
  poly::Integer above = poly::ceil(value);
  // construct var \in (below, above)
  return CACInterval{poly::Interval(below, above),
                     {var - below},
                     {var - above},
                     {var - below, var - above},
                     {},
                     {}};
}

bool CDCAC::hasRootAbove(const poly::Polynomial& p,
                         const poly::Value& val) const
{
  auto roots = poly::isolate_real_roots(p, d_assignment);
  return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
    return r >= val;
  });
}

bool CDCAC::hasRootBelow(const poly::Polynomial& p,
                         const poly::Value& val) const
{
  auto roots = poly::isolate_real_roots(p, d_assignment);
  return std::any_of(roots.begin(), roots.end(), [&val](const poly::Value& r) {
    return r <= val;
  });
}

void CDCAC::pruneRedundantIntervals(std::vector<CACInterval>& intervals)
{
  if (isProofEnabled())
  {
    std::vector<CACInterval> allIntervals = intervals;
    cleanIntervals(intervals);
    d_proof->pruneChildren([&allIntervals, &intervals](std::size_t i) {
      return std::find(intervals.begin(), intervals.end(), allIntervals[i])
             != intervals.end();
    });
  }
  else
  {
    cleanIntervals(intervals);
  }
}

}  // namespace cad
}  // namespace nl
}  // namespace arith
}  // namespace theory
}  // namespace cvc5

#endif


Generated by: GCOVR (Version 3.2)

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