diff options
| -rw-r--r-- | src/CMakeLists.txt | 8 | ||||
| -rw-r--r-- | src/theory/arith/nl/cad/constraints.cpp | 92 | ||||
| -rw-r--r-- | src/theory/arith/nl/cad/constraints.h | 104 | ||||
| -rw-r--r-- | src/theory/arith/nl/cad/projections.cpp | 132 | ||||
| -rw-r--r-- | src/theory/arith/nl/cad/projections.h | 80 | ||||
| -rw-r--r-- | src/theory/arith/nl/cad/variable_ordering.cpp | 146 | ||||
| -rw-r--r-- | src/theory/arith/nl/cad/variable_ordering.h | 78 | ||||
| -rw-r--r-- | src/theory/arith/nl/nl_model.cpp | 39 | ||||
| -rw-r--r-- | src/theory/arith/nl/nl_model.h | 73 | ||||
| -rw-r--r-- | src/theory/arith/nl/nonlinear_extension.cpp | 7 | ||||
| -rw-r--r-- | src/theory/arith/nl/nonlinear_extension.h | 5 | ||||
| -rw-r--r-- | src/theory/arith/nl/poly_conversion.cpp | 516 | ||||
| -rw-r--r-- | src/theory/arith/nl/poly_conversion.h | 135 | ||||
| -rw-r--r-- | src/util/poly_util.cpp | 91 | ||||
| -rw-r--r-- | src/util/poly_util.h | 15 |
15 files changed, 1481 insertions, 40 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 306813122..422888e73 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -306,6 +306,12 @@ libcvc4_add_sources( theory/arith/linear_equality.h theory/arith/matrix.cpp theory/arith/matrix.h + theory/arith/nl/cad/constraints.cpp + theory/arith/nl/cad/constraints.h + theory/arith/nl/cad/projections.cpp + theory/arith/nl/cad/projections.h + theory/arith/nl/cad/variable_ordering.cpp + theory/arith/nl/cad/variable_ordering.h theory/arith/nl/iand_solver.cpp theory/arith/nl/iand_solver.h theory/arith/nl/inference.cpp @@ -322,6 +328,8 @@ libcvc4_add_sources( theory/arith/nl/nl_solver.h theory/arith/nl/nonlinear_extension.cpp theory/arith/nl/nonlinear_extension.h + theory/arith/nl/poly_conversion.cpp + theory/arith/nl/poly_conversion.h theory/arith/nl/stats.cpp theory/arith/nl/stats.h theory/arith/nl/transcendental_solver.cpp diff --git a/src/theory/arith/nl/cad/constraints.cpp b/src/theory/arith/nl/cad/constraints.cpp new file mode 100644 index 000000000..b1b9edccc --- /dev/null +++ b/src/theory/arith/nl/cad/constraints.cpp @@ -0,0 +1,92 @@ +/********************* */ +/*! \file constraints.cpp + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Implements a container for CAD constraints. + ** + ** Implements a container for CAD constraints. + **/ + +#include "constraints.h" + +#ifdef CVC4_POLY_IMP + +#include <algorithm> + +#include "../poly_conversion.h" +#include "util/poly_util.h" + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { +namespace cad { + +using namespace poly; + +void Constraints::sort_constraints() +{ + using Tpl = std::tuple<poly::Polynomial, poly::SignCondition, Node>; + std::sort(mConstraints.begin(), + mConstraints.end(), + [](const Tpl& at, const Tpl& bt) { + // Check if a is smaller than b + const poly::Polynomial& a = std::get<0>(at); + const poly::Polynomial& b = std::get<0>(bt); + bool ua = is_univariate(a); + bool ub = is_univariate(b); + if (ua != ub) return ua; + std::size_t tda = poly_utils::totalDegree(a); + std::size_t tdb = poly_utils::totalDegree(b); + if (tda != tdb) return tda < tdb; + return degree(a) < degree(b); + }); + for (auto& c : mConstraints) + { + auto* p = std::get<0>(c).get_internal(); + lp_polynomial_set_external(p); + } +} + +void Constraints::add_constraint(const Polynomial& lhs, + SignCondition sc, + Node n) +{ + mConstraints.emplace_back(lhs, sc, n); + sort_constraints(); +} + +void Constraints::add_constraint(Node n) +{ + auto c = as_poly_constraint(n, mVarMapper); + add_constraint(c.first, c.second, n); + sort_constraints(); +} + +const Constraints::ConstraintVector& Constraints::get_constraints() const +{ + return mConstraints; +} + +void Constraints::reset() { mConstraints.clear(); } + +} // namespace cad +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif diff --git a/src/theory/arith/nl/cad/constraints.h b/src/theory/arith/nl/cad/constraints.h new file mode 100644 index 000000000..eda785d8e --- /dev/null +++ b/src/theory/arith/nl/cad/constraints.h @@ -0,0 +1,104 @@ +/********************* */ +/*! \file constraints.h + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Implements a container for CAD constraints. + ** + ** Implements a container for CAD constraints. + **/ + +#ifndef CVC4__THEORY__ARITH__NL__CAD__CONSTRAINTS_H +#define CVC4__THEORY__ARITH__NL__CAD__CONSTRAINTS_H + +#include "util/real_algebraic_number.h" + +#ifdef CVC4_POLY_IMP + +#include <poly/polyxx.h> + +#include <map> +#include <tuple> +#include <vector> + +#include "../poly_conversion.h" +#include "expr/kind.h" +#include "expr/node_manager_attributes.h" + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { +namespace cad { + +class Constraints +{ + public: + /** Type alias for a list of constraints. */ + using Constraint = std::tuple<poly::Polynomial, poly::SignCondition, Node>; + using ConstraintVector = std::vector<Constraint>; + + private: + /** + * A list of constraints, each comprised of a polynomial and a sign + * condition. + */ + ConstraintVector mConstraints; + + /** + * A mapping from CVC4 variables to poly variables. + */ + VariableMapper mVarMapper; + + void sort_constraints(); + + public: + VariableMapper& var_mapper() { return mVarMapper; } + + /** + * Add a constraint (represented by a polynomial and a sign condition) to the + * list of constraints. + */ + void add_constraint(const poly::Polynomial& lhs, + poly::SignCondition sc, + Node n); + + /** + * Add a constraints (represented by a node) to the list of constraints. + * The given node can either be a negation (NOT) or a suitable relation symbol + * as checked by is_suitable_relation(). + */ + void add_constraint(Node n); + + /** + * Gives the list of added constraints. + */ + const ConstraintVector& get_constraints() const; + + /** + * Remove all constraints. + */ + void reset(); +}; + +} // namespace cad +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif + +#endif diff --git a/src/theory/arith/nl/cad/projections.cpp b/src/theory/arith/nl/cad/projections.cpp new file mode 100644 index 000000000..87d8a3945 --- /dev/null +++ b/src/theory/arith/nl/cad/projections.cpp @@ -0,0 +1,132 @@ +/********************* */ +/*! \file projections.cpp + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Implements utilities for CAD projection operators. + ** + ** Implements utilities for CAD projection operators. + **/ + +#include "projections.h" + +#ifdef CVC4_POLY_IMP + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { +namespace cad { + +using namespace poly; + +void reduce_projection_polynomials(std::vector<Polynomial>& polys) +{ + std::sort(polys.begin(), polys.end()); + auto it = std::unique(polys.begin(), polys.end()); + polys.erase(it, polys.end()); +} + +void add_polynomial(std::vector<Polynomial>& polys, const Polynomial& poly) +{ + for (const auto& p : square_free_factors(poly)) + { + if (is_constant(p)) continue; + polys.emplace_back(p); + } +} + +void add_polynomials(std::vector<Polynomial>& polys, + const std::vector<Polynomial>& p) +{ + for (const auto& q : p) add_polynomial(polys, q); +} + +void make_finest_square_free_basis(std::vector<Polynomial>& polys) +{ + for (std::size_t i = 0, n = polys.size(); i < n; ++i) + { + for (std::size_t j = i + 1; j < n; ++j) + { + Polynomial g = gcd(polys[i], polys[j]); + if (!is_constant(g)) + { + polys[i] = div(polys[i], g); + polys[j] = div(polys[j], g); + polys.emplace_back(g); + } + } + } + auto it = std::remove_if(polys.begin(), polys.end(), [](const Polynomial& p) { + return is_constant(p); + }); + polys.erase(it, polys.end()); + reduce_projection_polynomials(polys); +} + +void make_finest_square_free_basis(std::vector<poly::Polynomial>& lhs, + std::vector<poly::Polynomial>& rhs) +{ + for (std::size_t i = 0, ln = lhs.size(); i < ln; ++i) + { + for (std::size_t j = 0, rn = rhs.size(); j < rn; ++j) + { + if (lhs[i] == rhs[j]) continue; + Polynomial g = gcd(lhs[i], rhs[j]); + if (!is_constant(g)) + { + lhs[i] = div(lhs[i], g); + rhs[j] = div(rhs[j], g); + lhs.emplace_back(g); + rhs.emplace_back(g); + } + } + } + reduce_projection_polynomials(lhs); + reduce_projection_polynomials(rhs); +} + +std::vector<Polynomial> projection_mccallum( + const std::vector<Polynomial>& polys) +{ + std::vector<Polynomial> res; + + for (const auto& p : polys) + { + for (const auto& coeff : coefficients(p)) + { + add_polynomial(res, coeff); + } + add_polynomial(res, discriminant(p)); + } + for (std::size_t i = 0, n = polys.size(); i < n; ++i) + { + for (std::size_t j = i + 1; j < n; ++j) + { + add_polynomial(res, resultant(polys[i], polys[j])); + } + } + + reduce_projection_polynomials(res); + return res; +} + +} // namespace cad +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif diff --git a/src/theory/arith/nl/cad/projections.h b/src/theory/arith/nl/cad/projections.h new file mode 100644 index 000000000..39d2b1221 --- /dev/null +++ b/src/theory/arith/nl/cad/projections.h @@ -0,0 +1,80 @@ +/********************* */ +/*! \file projections.h + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Implements utilities for CAD projection operators. + ** + ** Implements utilities for CAD projection operators. + **/ + +#ifndef CVC4__THEORY__ARITH__NL__CAD_PROJECTIONS_H +#define CVC4__THEORY__ARITH__NL__CAD_PROJECTIONS_H + +#include "util/real_algebraic_number.h" + +#ifdef CVC4_USE_POLY + +#include <poly/polyxx.h> + +#include <algorithm> +#include <iostream> +#include <vector> + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { +namespace cad { + +/** Sort and remove duplicates from the list of polynomials. */ +void reduce_projection_polynomials(std::vector<poly::Polynomial>& polys); + +/** + * Adds a polynomial to the list of projection polynomials. + * Before adding, it factorizes the polynomials and removed constant factors. + */ +void add_polynomial(std::vector<poly::Polynomial>& polys, + const poly::Polynomial& poly); + +/** Adds a list of polynomials using add_polynomial(). */ +void add_polynomials(std::vector<poly::Polynomial>& polys, + const std::vector<poly::Polynomial>& p); + +/** Make a set of polynomials a finest square-free basis. */ +void make_finest_square_free_basis(std::vector<poly::Polynomial>& polys); + +/** + * Ensure that two sets of polynomials are finest square-free basis relative to + * each other. + */ +void make_finest_square_free_basis(std::vector<poly::Polynomial>& lhs, + std::vector<poly::Polynomial>& rhs); + +/** + * Computes McCallum's projection operator. + */ +std::vector<poly::Polynomial> projection_mccallum( + const std::vector<poly::Polynomial>& polys); + +} // namespace cad +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif + +#endif diff --git a/src/theory/arith/nl/cad/variable_ordering.cpp b/src/theory/arith/nl/cad/variable_ordering.cpp new file mode 100644 index 000000000..efb2e6350 --- /dev/null +++ b/src/theory/arith/nl/cad/variable_ordering.cpp @@ -0,0 +1,146 @@ +/********************* */ +/*! \file variable_ordering.cpp + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Implements variable orderings tailored to CAD. + ** + ** Implements variable orderings tailored to CAD. + **/ + +#include "variable_ordering.h" + +#ifdef CVC4_POLY_IMP + +#include "util/poly_util.h" + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { +namespace cad { + +using namespace poly; + +std::vector<poly_utils::VariableInformation> collect_information( + const Constraints::ConstraintVector& polys, bool with_totals) +{ + VariableCollector vc; + for (const auto& c : polys) + { + vc(std::get<0>(c)); + } + std::vector<poly_utils::VariableInformation> res; + for (const auto& v : vc.get_variables()) + { + res.emplace_back(); + res.back().var = v; + for (const auto& c : polys) + { + poly_utils::getVariableInformation(res.back(), std::get<0>(c)); + } + } + if (with_totals) + { + res.emplace_back(); + for (const auto& c : polys) + { + poly_utils::getVariableInformation(res.back(), std::get<0>(c)); + } + } + return res; +} + +std::vector<poly::Variable> get_variables( + const std::vector<poly_utils::VariableInformation>& vi) +{ + std::vector<poly::Variable> res; + for (const auto& v : vi) + { + res.emplace_back(v.var); + } + return res; +} + +std::vector<poly::Variable> sort_byid( + const Constraints::ConstraintVector& polys) +{ + auto vi = collect_information(polys); + std::sort( + vi.begin(), + vi.end(), + [](const poly_utils::VariableInformation& a, + const poly_utils::VariableInformation& b) { return a.var < b.var; }); + return get_variables(vi); +}; + +std::vector<poly::Variable> sort_brown( + const Constraints::ConstraintVector& polys) +{ + auto vi = collect_information(polys); + std::sort(vi.begin(), + vi.end(), + [](const poly_utils::VariableInformation& a, + const poly_utils::VariableInformation& b) { + if (a.max_degree != b.max_degree) + return a.max_degree > b.max_degree; + if (a.max_terms_tdegree != b.max_terms_tdegree) + return a.max_terms_tdegree > b.max_terms_tdegree; + return a.num_terms > b.num_terms; + }); + return get_variables(vi); +}; + +std::vector<poly::Variable> sort_triangular( + const Constraints::ConstraintVector& polys) +{ + auto vi = collect_information(polys); + std::sort(vi.begin(), + vi.end(), + [](const poly_utils::VariableInformation& a, + const poly_utils::VariableInformation& b) { + if (a.max_degree != b.max_degree) + return a.max_degree > b.max_degree; + if (a.max_lc_degree != b.max_lc_degree) + return a.max_lc_degree > b.max_lc_degree; + return a.sum_poly_degree > b.sum_poly_degree; + }); + return get_variables(vi); +}; + +VariableOrdering::VariableOrdering() {} +VariableOrdering::~VariableOrdering() {} + +std::vector<poly::Variable> VariableOrdering::operator()( + const Constraints::ConstraintVector& polys, + VariableOrderingStrategy vos) const +{ + switch (vos) + { + case VariableOrderingStrategy::BYID: return sort_byid(polys); + case VariableOrderingStrategy::BROWN: return sort_brown(polys); + case VariableOrderingStrategy::TRIANGULAR: return sort_triangular(polys); + default: Assert(false) << "Unsupported variable ordering."; + } + return {}; +} + +} // namespace cad +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif diff --git a/src/theory/arith/nl/cad/variable_ordering.h b/src/theory/arith/nl/cad/variable_ordering.h new file mode 100644 index 000000000..3d4bee861 --- /dev/null +++ b/src/theory/arith/nl/cad/variable_ordering.h @@ -0,0 +1,78 @@ +/********************* */ +/*! \file variable_ordering.h + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Implements variable orderings tailored to CAD. + ** + ** Implements variable orderings tailored to CAD. + **/ + +#ifndef CVC4__THEORY__ARITH__NL__CAD__VARIABLE_ORDERING_H +#define CVC4__THEORY__ARITH__NL__CAD__VARIABLE_ORDERING_H + +#include "util/real_algebraic_number.h" + +#ifdef CVC4_POLY_IMP + +#include <poly/polyxx.h> + +#include "constraints.h" +#include "util/poly_util.h" + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { +namespace cad { + +/** Variable orderings for real variables in the context of CAD. */ +enum class VariableOrderingStrategy +{ + /** Dummy ordering by variable ID. */ + BYID, + /** Triangular as of DOI:10.1145/2755996.2756678 */ + TRIANGULAR, + /** Brown as of DOI:10.1145/2755996.2756678 */ + BROWN +}; + +class VariableOrdering +{ + public: + VariableOrdering(); + ~VariableOrdering(); + std::vector<poly::Variable> operator()( + const Constraints::ConstraintVector& polys, + VariableOrderingStrategy vos) const; +}; + +/** + * Retrieves variable information for all variables with the given polynomials. + * If with_totals is set, the last element of the vector contains totals as + * computed by get_variable_information if no variable is specified. + */ +std::vector<poly_utils::VariableInformation> collect_information( + const Constraints::ConstraintVector& polys, bool with_totals = false); + +} // namespace cad +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif + +#endif diff --git a/src/theory/arith/nl/nl_model.cpp b/src/theory/arith/nl/nl_model.cpp index d471bf0f6..daa36f972 100644 --- a/src/theory/arith/nl/nl_model.cpp +++ b/src/theory/arith/nl/nl_model.cpp @@ -57,6 +57,7 @@ void NlModel::resetCheck() d_used_approx = false; d_check_model_solved.clear(); d_check_model_bounds.clear(); + d_check_model_witnesses.clear(); d_check_model_vars.clear(); d_check_model_subs.clear(); } @@ -362,6 +363,9 @@ bool NlModel::addCheckModelSubstitution(TNode v, TNode s) return false; } } + Assert(d_check_model_witnesses.find(v) == d_check_model_witnesses.end()) + << "We tried to add a substitution where we already had a witness term." + << std::endl; std::vector<Node> varsTmp; varsTmp.push_back(v); std::vector<Node> subsTmp; @@ -415,12 +419,34 @@ bool NlModel::addCheckModelBound(TNode v, TNode l, TNode u) return true; } +bool NlModel::addCheckModelWitness(TNode v, TNode w) +{ + Trace("nl-ext-model") << "* check model witness : " << v << " -> " << w + << std::endl; + // should not set a witness for a value that is already set + if (std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) + != d_check_model_vars.end()) + { + Trace("nl-ext-model") << "...ERROR: setting witness for variable that " + "already has a constant value." + << std::endl; + Assert(false); + return false; + } + d_check_model_witnesses.emplace(v, w); + return true; +} + bool NlModel::hasCheckModelAssignment(Node v) const { if (d_check_model_bounds.find(v) != d_check_model_bounds.end()) { return true; } + if (d_check_model_witnesses.find(v) != d_check_model_witnesses.end()) + { + return true; + } return std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) != d_check_model_vars.end(); } @@ -1087,6 +1113,11 @@ bool NlModel::simpleCheckModelMsum(const std::map<Node, Node>& msum, bool pol) } else { + Assert(d_check_model_witnesses.find(vc) + == d_check_model_witnesses.end()) + << "No variable should be assigned a witness term if we get " + "here. " + << vc << " is, though." << std::endl; Trace("nl-ext-cms-debug") << std::endl; Trace("nl-ext-cms") << " failed due to unknown bound for " << vc << std::endl; @@ -1276,7 +1307,8 @@ void NlModel::printModelValue(const char* c, Node n, unsigned prec) const void NlModel::getModelValueRepair( std::map<Node, Node>& arithModel, - std::map<Node, std::pair<Node, Node>>& approximations) + std::map<Node, std::pair<Node, Node>>& approximations, + std::map<Node, Node>& witnesses) { Trace("nl-model") << "NlModel::getModelValueRepair:" << std::endl; @@ -1314,6 +1346,11 @@ void NlModel::getModelValueRepair( Trace("nl-model") << v << " exact approximation is " << l << std::endl; } } + for (const auto& vw : d_check_model_witnesses) + { + Trace("nl-model") << vw.first << " witness is " << vw.second << std::endl; + witnesses.emplace(vw.first, vw.second); + } // Also record the exact values we used. An exact value can be seen as a // special kind approximation of the form (witness x. x = exact_value). // Notice that the above term gets rewritten such that the choice function diff --git a/src/theory/arith/nl/nl_model.h b/src/theory/arith/nl/nl_model.h index fdce446fc..1e6f6c8f3 100644 --- a/src/theory/arith/nl/nl_model.h +++ b/src/theory/arith/nl/nl_model.h @@ -48,21 +48,18 @@ class NlModel public: NlModel(context::Context* c); ~NlModel(); - /** reset - * + /** * This method is called once at the beginning of a last call effort check, * where m is the model of the theory of arithmetic. This method resets the * cache of computed model values. */ void reset(TheoryModel* m, std::map<Node, Node>& arithModel); - /** reset check - * + /** * This method is called when the non-linear arithmetic solver restarts * its computation of lemmas and models during a last call effort check. */ void resetCheck(); - /** compute model value - * + /** * This computes model values for terms based on two semantics, a "concrete" * semantics and an "abstract" semantics. * @@ -88,7 +85,8 @@ class NlModel Node computeAbstractModelValue(Node n); Node computeModelValue(Node n, bool isConcrete); - /** Compare arithmetic terms i and j based an ordering. + /** + * Compare arithmetic terms i and j based an ordering. * * This returns: * -1 if i < j, 1 if i > j, or 0 if i == j @@ -101,7 +99,8 @@ class NlModel * values. */ int compare(Node i, Node j, bool isConcrete, bool isAbsolute); - /** Compare arithmetic terms i and j based an ordering. + /** + * Compare arithmetic terms i and j based an ordering. * * This returns: * -1 if i < j, 1 if i > j, or 0 if i == j @@ -111,8 +110,7 @@ class NlModel int compareValue(Node i, Node j, bool isAbsolute) const; //------------------------------ recording model substitutions and bounds - /** add check model substitution - * + /** * Adds the model substitution v -> s. This applies the substitution * { v -> s } to each term in d_check_model_subs and adds v,s to * d_check_model_vars and d_check_model_subs respectively. @@ -120,24 +118,28 @@ class NlModel * with the current substitution and bounds. */ bool addCheckModelSubstitution(TNode v, TNode s); - /** add check model bound - * + /** * Adds the bound x -> < l, u > to the map above, and records the * approximation ( x, l <= x <= u ) in the model. This method returns false * if the bound is inconsistent with the current model substitution or * bounds. */ bool addCheckModelBound(TNode v, TNode l, TNode u); - /** has check model assignment - * + /** + * Adds a model witness v -> w to the underlying theory model. + * The witness should only contain a single variable v and evaluate to true + * for exactly one value of v. The variable v is then (implicitly, + * declaratively) assigned to this single value that satisfies the witness w. + */ + bool addCheckModelWitness(TNode v, TNode w); + /** * Have we assigned v in the current checkModel(...) call? * * This method returns true if variable v is in the domain of * d_check_model_bounds or if it occurs in d_check_model_vars. */ bool hasCheckModelAssignment(Node v) const; - /** Check model - * + /** * Checks the current model based on solving for equalities, and using error * bounds on the Taylor approximation. * @@ -164,8 +166,7 @@ class NlModel void setUsedApproximate(); /** Did we use an approximation during this check? */ bool usedApproximate() const; - /** Set tautology - * + /** * This states that formula n is a tautology (satisfied in all models). * We call this on internally generated lemmas. This method computes a * set of literals that are implied by n, that are hence tautological @@ -183,14 +184,12 @@ class NlModel void addTautology(Node n); //------------------------------ end recording model substitutions and bounds - /** print model value, for debugging. - * + /** * This prints both the abstract and concrete model values for arithmetic * term n on Trace c with precision prec. */ void printModelValue(const char* c, Node n, unsigned prec = 5) const; - /** get model value repair - * + /** * This gets mappings that indicate how to repair the model generated by the * linear arithmetic solver. This method should be called after a successful * call to checkModel above. @@ -199,11 +198,14 @@ class NlModel * to their (exact) value that was computed during checkModel; the mapping * approximations is updated to store approximate values in the form of a * pair (P, w), where P is a predicate that describes the possible values of - * v and w is a witness point that satisfies this predicate. + * v and w is a witness point that satisfies this predicate; the mapping + * witnesses is filled with witness terms that are satisfied by a single + * value. */ void getModelValueRepair( std::map<Node, Node>& arithModel, - std::map<Node, std::pair<Node, Node>>& approximations); + std::map<Node, std::pair<Node, Node>>& approximations, + std::map<Node, Node>& witnesses); private: /** The current model */ @@ -216,8 +218,7 @@ class NlModel Node getRepresentative(Node n) const; //---------------------------check model - /** solve equality simple - * + /** * This method is used during checkModel(...). It takes as input an * equality eq. If it returns true, then eq is correct-by-construction based * on the information stored in our model representation (see @@ -235,7 +236,8 @@ class NlModel */ bool solveEqualitySimple(Node eq, unsigned d, std::vector<NlLemma>& lemmas); - /** simple check model for transcendental functions for literal + /** + * simple check model for transcendental functions for literal * * This method returns true if literal is true for all interpretations of * transcendental functions within their error bounds (as stored @@ -257,8 +259,7 @@ class NlModel bool simpleCheckModelLit(Node lit); bool simpleCheckModelMsum(const std::map<Node, Node>& msum, bool pol); //---------------------------end check model - /** get approximate sqrt - * + /** * This approximates the square root of positive constant c. If this method * returns true, then l and u are updated to constants such that * l <= sqrt( c ) <= u @@ -282,7 +283,8 @@ class NlModel * sending to TheoryModel during collectModelInfo. */ std::map<Node, Node> d_arithVal; - /** cache of model values + /** + * cache of model values * * Stores the the concrete/abstract model values. This is a cache of the * computeModelValue method. @@ -296,7 +298,8 @@ class NlModel */ std::vector<Node> d_check_model_vars; std::vector<Node> d_check_model_subs; - /** lower and upper bounds for check model + /** + * lower and upper bounds for check model * * For each term t in the domain of this map, if this stores the pair * (c_l, c_u) then the model M is such that c_l <= M( t ) <= c_u. @@ -310,6 +313,14 @@ class NlModel */ std::map<Node, std::pair<Node, Node>> d_check_model_bounds; /** + * witnesses for check model + * + * Stores witnesses for vatiables that define implicit variable assignments. + * For some variable v, we map to a formulas that is true for exactly one + * value of v. + */ + std::map<Node, Node> d_check_model_witnesses; + /** * The map from literals that our model construction solved, to the variable * that was solved for. Examples of such literals are: * (1) Equalities x = t, which we turned into a model substitution x -> t, diff --git a/src/theory/arith/nl/nonlinear_extension.cpp b/src/theory/arith/nl/nonlinear_extension.cpp index 12772f4d2..432c25f27 100644 --- a/src/theory/arith/nl/nonlinear_extension.cpp +++ b/src/theory/arith/nl/nonlinear_extension.cpp @@ -651,6 +651,10 @@ void NonlinearExtension::check(Theory::Effort e) tm->recordApproximation(a.first, a.second.first, a.second.second); } } + for (const auto& vw : d_witnesses) + { + tm->recordApproximation(vw.first, vw.second); + } } } @@ -870,8 +874,9 @@ void NonlinearExtension::interceptModel(std::map<Node, Node>& arithModel) { Trace("nl-ext") << "interceptModel: do model repair" << std::endl; d_approximations.clear(); + d_witnesses.clear(); // modify the model values - d_model.getModelValueRepair(arithModel, d_approximations); + d_model.getModelValueRepair(arithModel, d_approximations, d_witnesses); } } diff --git a/src/theory/arith/nl/nonlinear_extension.h b/src/theory/arith/nl/nonlinear_extension.h index b4e5df976..36875d8a3 100644 --- a/src/theory/arith/nl/nonlinear_extension.h +++ b/src/theory/arith/nl/nonlinear_extension.h @@ -330,6 +330,11 @@ class NonlinearExtension * NlModel::getModelValueRepair. */ std::map<Node, std::pair<Node, Node>> d_approximations; + /** + * The witnesses computed during collectModelInfo. For details, see + * NlModel::getModelValueRepair. + */ + std::map<Node, Node> d_witnesses; /** have we successfully built the model in this SAT context? */ context::CDO<bool> d_builtModel; }; /* class NonlinearExtension */ diff --git a/src/theory/arith/nl/poly_conversion.cpp b/src/theory/arith/nl/poly_conversion.cpp new file mode 100644 index 000000000..8292380f9 --- /dev/null +++ b/src/theory/arith/nl/poly_conversion.cpp @@ -0,0 +1,516 @@ +/********************* */ +/*! \file poly_conversion.cpp + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Utilities for converting to and from LibPoly objects. + ** + ** Utilities for converting to and from LibPoly objects. + **/ + +#include "poly_conversion.h" + +#ifdef CVC4_POLY_IMP + +#include "expr/node.h" +#include "expr/node_manager_attributes.h" +#include "util/integer.h" +#include "util/poly_util.h" +#include "util/rational.h" + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { + +poly::Variable VariableMapper::operator()(const CVC4::Node& n) +{ + auto it = mVarCVCpoly.find(n); + if (it == mVarCVCpoly.end()) + { + std::string name; + if (n.isVar()) + { + if (!n.getAttribute(expr::VarNameAttr(), name)) + { + Trace("poly::conversion") + << "Variable " << n << " has no name, using ID instead." + << std::endl; + name = "v_" + std::to_string(n.getId()); + } + } + else + { + name = "v_" + std::to_string(n.getId()); + } + it = mVarCVCpoly.emplace(n, poly::Variable(name.c_str())).first; + mVarpolyCVC.emplace(it->second, n); + } + return it->second; +} + +CVC4::Node VariableMapper::operator()(const poly::Variable& n) +{ + auto it = mVarpolyCVC.find(n); + Assert(it != mVarpolyCVC.end()) + << "Expect variable " << n << " to be added already."; + return it->second; +} + +CVC4::Node as_cvc_upolynomial(const poly::UPolynomial& p, const CVC4::Node& var) +{ + Trace("poly::conversion") + << "Converting " << p << " over " << var << std::endl; + + std::vector<poly::Integer> coeffs = coefficients(p); + + auto* nm = NodeManager::currentNM(); + + Node res = nm->mkConst(Rational(0)); + Node monomial = nm->mkConst(Rational(1)); + for (std::size_t i = 0, n = coeffs.size(); i < n; ++i) + { + if (!is_zero(coeffs[i])) + { + Node coeff = nm->mkConst(poly_utils::toRational(coeffs[i])); + Node term = nm->mkNode(Kind::MULT, coeff, monomial); + res = nm->mkNode(Kind::PLUS, res, term); + } + monomial = nm->mkNode(Kind::NONLINEAR_MULT, monomial, var); + } + Trace("poly::conversion") << "-> " << res << std::endl; + return res; +} + +poly::UPolynomial as_poly_upolynomial_impl(const CVC4::Node& n, + poly::Integer& denominator, + const CVC4::Node& var) +{ + denominator = poly::Integer(1); + if (n.isVar()) + { + Assert(n == var) << "Unexpected variable: should be " << var << " but is " + << n; + return poly::UPolynomial({0, 1}); + } + switch (n.getKind()) + { + case Kind::CONST_RATIONAL: + { + Rational r = n.getConst<Rational>(); + denominator = poly_utils::toInteger(r.getDenominator()); + return poly::UPolynomial(poly_utils::toInteger(r.getNumerator())); + } + case Kind::PLUS: + { + poly::UPolynomial res; + poly::Integer denom; + for (const auto& child : n) + { + poly::UPolynomial tmp = as_poly_upolynomial_impl(child, denom, var); + /** Normalize denominators + */ + poly::Integer g = gcd(denom, denominator); + res = res * (denom / g) + tmp * (denominator / g); + denominator *= (denom / g); + } + return res; + } + case Kind::MULT: + case Kind::NONLINEAR_MULT: + { + poly::UPolynomial res(denominator); + poly::Integer denom; + for (const auto& child : n) + { + res = res * as_poly_upolynomial_impl(child, denom, var); + denominator *= denom; + } + return res; + } + default: + Warning() << "Unhandled node " << n << " with kind " << n.getKind() + << std::endl; + } + return poly::UPolynomial(); +} + +poly::UPolynomial as_poly_upolynomial(const CVC4::Node& n, + const CVC4::Node& var) +{ + poly::Integer denom; + return as_poly_upolynomial_impl(n, denom, var); +} + +poly::Polynomial as_poly_polynomial_impl(const CVC4::Node& n, + poly::Integer& denominator, + VariableMapper& vm) +{ + denominator = poly::Integer(1); + if (n.isVar()) + { + return poly::Polynomial(vm(n)); + } + switch (n.getKind()) + { + case Kind::CONST_RATIONAL: + { + Rational r = n.getConst<Rational>(); + denominator = poly_utils::toInteger(r.getDenominator()); + return poly::Polynomial(poly_utils::toInteger(r.getNumerator())); + } + case Kind::PLUS: + { + poly::Polynomial res; + poly::Integer denom; + for (const auto& child : n) + { + poly::Polynomial tmp = as_poly_polynomial_impl(child, denom, vm); + /** Normalize denominators + */ + poly::Integer g = gcd(denom, denominator); + res = res * (denom / g) + tmp * (denominator / g); + denominator *= (denom / g); + } + return res; + } + case Kind::MULT: + case Kind::NONLINEAR_MULT: + { + poly::Polynomial res(denominator); + poly::Integer denom; + for (const auto& child : n) + { + res *= as_poly_polynomial_impl(child, denom, vm); + denominator *= denom; + } + return res; + } + default: return poly::Polynomial(vm(n)); + } + return poly::Polynomial(); +} +poly::Polynomial as_poly_polynomial(const CVC4::Node& n, VariableMapper& vm) +{ + poly::Integer denom; + return as_poly_polynomial_impl(n, denom, vm); +} + +poly::SignCondition normalize_kind(CVC4::Kind kind, + bool negated, + poly::Polynomial& lhs) +{ + switch (kind) + { + case Kind::EQUAL: + { + return negated ? poly::SignCondition::NE : poly::SignCondition::EQ; + } + case Kind::LT: + { + if (negated) + { + lhs = -lhs; + return poly::SignCondition::LE; + } + return poly::SignCondition::LT; + } + case Kind::LEQ: + { + if (negated) + { + lhs = -lhs; + return poly::SignCondition::LT; + } + return poly::SignCondition::LE; + } + case Kind::GT: + { + if (negated) + { + return poly::SignCondition::LE; + } + lhs = -lhs; + return poly::SignCondition::LT; + } + case Kind::GEQ: + { + if (negated) + { + return poly::SignCondition::LT; + } + lhs = -lhs; + return poly::SignCondition::LE; + } + default: + Assert(false) << "This function only deals with arithmetic relations."; + return poly::SignCondition::EQ; + } +} + +std::pair<poly::Polynomial, poly::SignCondition> as_poly_constraint( + Node n, VariableMapper& vm) +{ + bool negated = false; + Node origin = n; + if (n.getKind() == Kind::NOT) + { + Assert(n.getNumChildren() == 1) + << "Expect negations to have a single child."; + negated = true; + n = *n.begin(); + } + Assert((n.getKind() == Kind::EQUAL) || (n.getKind() == Kind::GT) + || (n.getKind() == Kind::GEQ) || (n.getKind() == Kind::LT) + || (n.getKind() == Kind::LEQ)) + << "Found a constraint with unsupported relation " << n.getKind(); + + Assert(n.getNumChildren() == 2) + << "Supported relations only have two children."; + auto childit = n.begin(); + poly::Integer ldenom; + poly::Polynomial left = as_poly_polynomial_impl(*childit++, ldenom, vm); + poly::Integer rdenom; + poly::Polynomial right = as_poly_polynomial_impl(*childit++, rdenom, vm); + Assert(childit == n.end()) << "Screwed up iterator handling."; + Assert(ldenom > poly::Integer(0) && rdenom > poly::Integer(0)) + << "Expected denominators to be always positive."; + + poly::Integer g = gcd(ldenom, rdenom); + poly::Polynomial lhs = left * (rdenom / g) - right * (ldenom / g); + poly::SignCondition sc = normalize_kind(n.getKind(), negated, lhs); + return {lhs, sc}; +} + +Node ran_to_node(const RealAlgebraicNumber& ran, const Node& ran_variable) +{ + return ran_to_node(ran.getValue(), ran_variable); +} + +Node ran_to_node(const poly::AlgebraicNumber& an, const Node& ran_variable) +{ + auto* nm = NodeManager::currentNM(); + + const poly::DyadicInterval& di = get_isolating_interval(an); + if (is_point(di)) + { + return nm->mkConst(poly_utils::toRational(get_point(di))); + } + Assert(di.get_internal()->a_open && di.get_internal()->b_open) + << "We assume an open interval here."; + + Node poly = as_cvc_upolynomial(get_defining_polynomial(an), ran_variable); + Node lower = nm->mkConst(poly_utils::toRational(get_lower(di))); + Node upper = nm->mkConst(poly_utils::toRational(get_upper(di))); + + // Construct witness: + return nm->mkNode(Kind::AND, + // poly(var) == 0 + nm->mkNode(Kind::EQUAL, poly, nm->mkConst(Rational(0))), + // lower_bound < var + nm->mkNode(Kind::LT, lower, ran_variable), + // var < upper_bound + nm->mkNode(Kind::LT, ran_variable, upper)); +} + +Node value_to_node(const poly::Value& v, const Node& ran_variable) +{ + Assert(!is_minus_infinity(v)) << "Can not convert minus infinity."; + Assert(!is_none(v)) << "Can not convert none."; + Assert(!is_plus_infinity(v)) << "Can not convert plus infinity."; + + if (is_algebraic_number(v)) + { + return ran_to_node(as_algebraic_number(v), ran_variable); + } + auto* nm = NodeManager::currentNM(); + if (is_dyadic_rational(v)) + { + return nm->mkConst(poly_utils::toRational(as_dyadic_rational(v))); + } + if (is_integer(v)) + { + return nm->mkConst(poly_utils::toRational(as_integer(v))); + } + if (is_rational(v)) + { + return nm->mkConst(poly_utils::toRational(as_rational(v))); + } + Assert(false) << "All cases should be covered."; + return nm->mkConst(Rational(0)); +} + +Node excluding_interval_to_lemma(const Node& variable, + const poly::Interval& interval) +{ + auto* nm = NodeManager::currentNM(); + const auto& lv = poly::get_lower(interval); + const auto& uv = poly::get_upper(interval); + bool li = poly::is_minus_infinity(lv); + bool ui = poly::is_plus_infinity(uv); + if (li && ui) return nm->mkConst(true); + if (poly::is_point(interval)) + { + if (is_algebraic_number(lv)) + { + return nm->mkNode( + Kind::AND, + nm->mkNode( + Kind::GT, variable, nm->mkConst(poly_utils::toRationalBelow(lv))), + nm->mkNode(Kind::LT, + variable, + nm->mkConst(poly_utils::toRationalAbove(lv)))); + } + else + { + return nm->mkNode( + Kind::EQUAL, variable, nm->mkConst(poly_utils::toRationalBelow(lv))); + } + } + if (li) + { + return nm->mkNode( + Kind::GT, variable, nm->mkConst(poly_utils::toRationalAbove(uv))); + } + if (ui) + { + return nm->mkNode( + Kind::LT, variable, nm->mkConst(poly_utils::toRationalBelow(lv))); + } + return nm->mkNode( + Kind::AND, + nm->mkNode( + Kind::GT, variable, nm->mkConst(poly_utils::toRationalAbove(uv))), + nm->mkNode( + Kind::LT, variable, nm->mkConst(poly_utils::toRationalBelow(lv)))); +} + +Maybe<Rational> get_lower_bound(const Node& n) +{ + if (n.getNumChildren() != 2) return Maybe<Rational>(); + if (n.getKind() == Kind::LT) + { + if (!n[0].isConst()) return Maybe<Rational>(); + if (!n[1].isVar()) return Maybe<Rational>(); + return n[0].getConst<Rational>(); + } + else if (n.getKind() == Kind::GT) + { + if (!n[0].isVar()) return Maybe<Rational>(); + if (!n[1].isConst()) return Maybe<Rational>(); + return n[1].getConst<Rational>(); + } + return Maybe<Rational>(); +} +Maybe<Rational> get_upper_bound(const Node& n) +{ + if (n.getNumChildren() != 2) return Maybe<Rational>(); + if (n.getKind() == Kind::LT) + { + if (!n[0].isVar()) return Maybe<Rational>(); + if (!n[1].isConst()) return Maybe<Rational>(); + return n[1].getConst<Rational>(); + } + else if (n.getKind() == Kind::GT) + { + if (!n[0].isConst()) return Maybe<Rational>(); + if (!n[1].isVar()) return Maybe<Rational>(); + return n[0].getConst<Rational>(); + } + return Maybe<Rational>(); +} + +/** Returns indices of appropriate parts of ran encoding. + * Returns (poly equation ; lower bound ; upper bound) + */ +std::tuple<Node, Rational, Rational> detect_ran_encoding(const Node& n) +{ + Assert(n.getKind() == Kind::AND) << "Invalid node structure."; + Assert(n.getNumChildren() == 3) << "Invalid node structure."; + + Node poly_eq; + if (n[0].getKind() == Kind::EQUAL) + poly_eq = n[0]; + else if (n[1].getKind() == Kind::EQUAL) + poly_eq = n[1]; + else if (n[2].getKind() == Kind::EQUAL) + poly_eq = n[2]; + else + Assert(false) << "Could not identify polynomial equation."; + + Node poly; + Assert(poly_eq.getNumChildren() == 2) << "Invalid polynomial equation."; + if (poly_eq[0].isConst()) + { + Assert(poly_eq[0].getConst<Rational>() == Rational(0)) + << "Invalid polynomial equation."; + poly = poly_eq[1]; + } + else if (poly_eq[1].isConst()) + { + Assert(poly_eq[1].getConst<Rational>() == Rational(0)) + << "Invalid polynomial equation."; + poly = poly_eq[0]; + } + else + { + Assert(false) << "Invalid polynomial equation."; + } + + Maybe<Rational> lower = get_lower_bound(n[0]); + if (!lower) lower = get_lower_bound(n[1]); + if (!lower) lower = get_lower_bound(n[2]); + Assert(lower) << "Could not identify lower bound."; + + Maybe<Rational> upper = get_upper_bound(n[0]); + if (!upper) upper = get_upper_bound(n[1]); + if (!upper) upper = get_upper_bound(n[2]); + Assert(upper) << "Could not identify upper bound."; + + return std::tuple<Node, Rational, Rational>( + poly, lower.value(), upper.value()); +} + +poly::AlgebraicNumber node_to_poly_ran(const Node& n, const Node& ran_variable) +{ + // Identify poly, lower and upper + auto encoding = detect_ran_encoding(n); + // Construct polynomial + poly::UPolynomial pol = + as_poly_upolynomial(std::get<0>(encoding), ran_variable); + // Construct algebraic number + return poly_utils::toPolyRanWithRefinement( + std::move(pol), std::get<1>(encoding), std::get<2>(encoding)); +} +RealAlgebraicNumber node_to_ran(const Node& n, const Node& ran_variable) +{ + return RealAlgebraicNumber(node_to_poly_ran(n, ran_variable)); +} + +poly::Value node_to_value(const Node& n, const Node& ran_variable) +{ + if (n.isConst()) + { + return poly_utils::toRational(n.getConst<Rational>()); + } + return node_to_poly_ran(n, ran_variable); +} + +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif diff --git a/src/theory/arith/nl/poly_conversion.h b/src/theory/arith/nl/poly_conversion.h new file mode 100644 index 000000000..73509f83c --- /dev/null +++ b/src/theory/arith/nl/poly_conversion.h @@ -0,0 +1,135 @@ +/********************* */ +/*! \file poly_conversion.h + ** \verbatim + ** Top contributors (to current version): + ** Gereon Kremer + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2020 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.\endverbatim + ** + ** \brief Utilities for converting to and from LibPoly objects. + ** + ** Utilities for converting to and from LibPoly objects. + **/ + +#ifndef CVC4__THEORY__ARITH__NL__POLY_CONVERSION_H +#define CVC4__THEORY__ARITH__NL__POLY_CONVERSION_H + +#include "util/real_algebraic_number.h" + +#ifdef CVC4_POLY_IMP + +#include <poly/polyxx.h> + +#include <iostream> + +#include "expr/node.h" +#include "util/real_algebraic_number.h" + +namespace CVC4 { +namespace theory { +namespace arith { +namespace nl { + +/** Bijective mapping between CVC4 variables and poly variables. */ +struct VariableMapper +{ + /** A mapping from CVC4 variables to poly variables. */ + std::map<CVC4::Node, poly::Variable> mVarCVCpoly; + /** A mapping from poly variables to CVC4 variables. */ + std::map<poly::Variable, CVC4::Node> mVarpolyCVC; + + /** Retrieves the according poly variable. */ + poly::Variable operator()(const CVC4::Node& n); + /** Retrieves the according CVC4 variable. */ + CVC4::Node operator()(const poly::Variable& n); +}; + +/** Convert a poly univariate polynomial to a CVC4::Node. */ +CVC4::Node as_cvc_upolynomial(const poly::UPolynomial& p, + const CVC4::Node& var); + +/** Convert a CVC4::Node to a poly univariate polynomial. */ +poly::UPolynomial as_poly_upolynomial(const CVC4::Node& n, + const CVC4::Node& var); + +/** + * Constructs a polynomial from the given node. + * + * While a Node may contain rationals, a Polynomial does not. + * We therefore also store the denominator of the returned polynomial and + * use it to construct the integer polynomial recursively. + * Once the polynomial has been fully constructed, we can ignore the + * denominator (except for its sign, which is always positive, though). + */ +poly::Polynomial as_poly_polynomial(const CVC4::Node& n, VariableMapper& vm); + +/** + * Constructs a constraints (a polynomial and a sign condition) from the given + * node. + */ +std::pair<poly::Polynomial, poly::SignCondition> as_poly_constraint( + Node n, VariableMapper& vm); + +/** + * Transforms a real algebraic number to a node suitable for putting it into a + * model. The resulting node can be either a constant (suitable for + * addCheckModelSubstitution) or a witness term (suitable for + * addCheckModelWitness). + */ +Node ran_to_node(const RealAlgebraicNumber& ran, const Node& ran_variable); + +Node ran_to_node(const poly::AlgebraicNumber& an, const Node& ran_variable); + +/** + * Transforms a poly::Value to a node. + * The resulting node can be either a constant or a witness term. + */ +Node value_to_node(const poly::Value& v, const Node& ran_variable); + +/** + * Constructs a lemma that excludes a given interval from the feasible values of + * a variable. The resulting lemma has the form + * (OR + * (<= var interval.lower) + * (<= interval.upper var) + * ) + */ +Node excluding_interval_to_lemma(const Node& variable, + const poly::Interval& interval); + +/** + * Transforms a node to a poly::AlgebraicNumber. + * Expects a node of the following form: + * (AND + * (= (polynomial in __z) 0) + * (< CONST __z) + * (< __z CONST) + * ) + */ +poly::AlgebraicNumber node_to_poly_ran(const Node& n, const Node& ran_variable); + +/** Transforms a node to a RealAlgebraicNumber by calling node_to_poly_ran. */ +RealAlgebraicNumber node_to_ran(const Node& n, const Node& ran_variable); + +/** + * Transforms a node to a poly::Value. + */ +poly::Value node_to_value(const Node& n, const Node& ran_variable); + +} // namespace nl +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif + +#endif diff --git a/src/util/poly_util.cpp b/src/util/poly_util.cpp index b4c5d1bf2..251ad7ea3 100644 --- a/src/util/poly_util.cpp +++ b/src/util/poly_util.cpp @@ -86,6 +86,48 @@ Rational toRational(const poly::DyadicRational& dr) { return Rational(toInteger(numerator(dr)), toInteger(denominator(dr))); } +Rational toRationalAbove(const poly::Value& v) +{ + if (is_algebraic_number(v)) + { + return toRational(get_upper_bound(as_algebraic_number(v))); + } + else if (is_dyadic_rational(v)) + { + return toRational(as_dyadic_rational(v)); + } + else if (is_integer(v)) + { + return toRational(as_integer(v)); + } + else if (is_rational(v)) + { + return toRational(as_rational(v)); + } + Assert(false) << "Can not convert " << v << " to rational."; + return Rational(); +} +Rational toRationalBelow(const poly::Value& v) +{ + if (is_algebraic_number(v)) + { + return toRational(get_lower_bound(as_algebraic_number(v))); + } + else if (is_dyadic_rational(v)) + { + return toRational(as_dyadic_rational(v)); + } + else if (is_integer(v)) + { + return toRational(as_integer(v)); + } + else if (is_rational(v)) + { + return toRational(as_rational(v)); + } + Assert(false) << "Can not convert " << v << " to rational."; + return Rational(); +} poly::Integer toInteger(const Integer& i) { @@ -232,6 +274,31 @@ std::size_t totalDegree(const poly::Polynomial& p) return tdeg; } +std::ostream& operator<<(std::ostream& os, const VariableInformation& vi) +{ + if (vi.var == poly::Variable()) + { + os << "Totals: "; + os << "max deg " << vi.max_degree; + os << ", sum term deg " << vi.sum_term_degree; + os << ", sum poly deg " << vi.sum_poly_degree; + os << ", num polys " << vi.num_polynomials; + os << ", num terms " << vi.num_terms; + } + else + { + os << "Info for " << vi.var << ": "; + os << "max deg " << vi.max_degree; + os << ", max lc deg: " << vi.max_lc_degree; + os << ", max term tdeg: " << vi.max_terms_tdegree; + os << ", sum term deg " << vi.sum_term_degree; + os << ", sum poly deg " << vi.sum_poly_degree; + os << ", num polys " << vi.num_polynomials; + os << ", num terms " << vi.num_terms; + } + return os; +} + struct GetVarInfo { VariableInformation* info; @@ -257,13 +324,19 @@ void getVariableInformation(VariableInformation& vi, if (m->p[i].x == info->var) { info->max_degree = std::max(info->max_degree, m->p[i].d); - info->sum_degree += m->p[i].d; - ++info->num_terms; + info->sum_term_degree += m->p[i].d; vardeg = m->p[i].d; } } - if (vardeg > 0) + if (info->var == poly::Variable()) + { + ++info->num_terms; + info->max_degree = std::max(info->max_degree, tdeg); + info->sum_term_degree += tdeg; + } + else if (vardeg > 0) { + ++info->num_terms; if (gvi->cur_var_degree < vardeg) { gvi->cur_lc_degree = tdeg - vardeg; @@ -271,9 +344,19 @@ void getVariableInformation(VariableInformation& vi, info->max_terms_tdegree = std::max(info->max_terms_tdegree, tdeg); } }; - + std::size_t tmp_max_degree = vi.max_degree; + std::size_t tmp_num_terms = vi.num_terms; + vi.max_degree = 0; + vi.num_terms = 0; lp_polynomial_traverse(poly.get_internal(), f, &varinfo); vi.max_lc_degree = std::max(vi.max_lc_degree, varinfo.cur_lc_degree); + if (vi.num_terms > 0) + { + ++vi.num_polynomials; + } + vi.sum_poly_degree += vi.max_degree; + vi.max_degree = std::max(vi.max_degree, tmp_max_degree); + vi.num_terms += tmp_num_terms; } } // namespace poly_utils diff --git a/src/util/poly_util.h b/src/util/poly_util.h index e54652002..0547c4f74 100644 --- a/src/util/poly_util.h +++ b/src/util/poly_util.h @@ -25,7 +25,6 @@ #ifndef CVC4__POLY_UTIL_H #define CVC4__POLY_UTIL_H - #include <vector> #include "maybe.h" @@ -59,6 +58,11 @@ Rational toRational(const poly::Rational& r); /** Converts a poly::DyadicRational to a CVC4::Rational. */ Rational toRational(const poly::DyadicRational& dr); +/** Converts a poly::Value to a CVC4::Rational (that may be a bit above). */ +Rational toRationalAbove(const poly::Value& v); +/** Converts a poly::Value to a CVC4::Rational (that may be a bit below). */ +Rational toRationalBelow(const poly::Value& v); + /** Converts a CVC4::Integer to a poly::Integer. */ poly::Integer toInteger(const Integer& i); /** Converts a vector of CVC4::Integers to a vector of poly::Integers. */ @@ -119,11 +123,16 @@ struct VariableInformation std::size_t max_lc_degree = 0; /** Maximum of total degrees of terms that contain this variable. */ std::size_t max_terms_tdegree = 0; - /** Sum of degrees of this variable. */ - std::size_t sum_degree = 0; + /** Sum of degrees of this variable within all terms. */ + std::size_t sum_term_degree = 0; + /** Sum of degrees of this variable within all polynomials. */ + std::size_t sum_poly_degree = 0; + /** Number of polynomials that contain this variable. */ + std::size_t num_polynomials = 0; /** Number of terms that contain this variable. */ std::size_t num_terms = 0; }; +std::ostream& operator<<(std::ostream& os, const VariableInformation& vi); void getVariableInformation(VariableInformation& vi, const poly::Polynomial& poly); |