From 04a8bf833bad57329a4e83b3c5aafb7537de885d Mon Sep 17 00:00:00 2001 From: Andrew Reynolds Date: Fri, 10 Apr 2020 19:15:18 -0500 Subject: [PATCH] Split the non-linear solver (#4219) This splits the "non-linear solver" from "NonlinearExtension". The non-linear solver is the module that implements the inference schemas whereas NonlinearExtension is the glue code that manages the solver(s) for non-linear. This also involves moving utilities from the non-linear solver to their own file. --- src/CMakeLists.txt | 6 + src/theory/arith/nl_constraint.cpp | 123 ++ src/theory/arith/nl_constraint.h | 86 ++ src/theory/arith/nl_lemma_utils.h | 4 +- src/theory/arith/nl_monomial.cpp | 334 +++++ src/theory/arith/nl_monomial.h | 147 ++ src/theory/arith/nl_solver.cpp | 1585 +++++++++++++++++++++ src/theory/arith/nl_solver.h | 368 +++++ src/theory/arith/nonlinear_extension.cpp | 1654 +--------------------- src/theory/arith/nonlinear_extension.h | 363 +---- src/theory/arith/transcendental_solver.h | 9 + 11 files changed, 2681 insertions(+), 1998 deletions(-) create mode 100644 src/theory/arith/nl_constraint.cpp create mode 100644 src/theory/arith/nl_constraint.h create mode 100644 src/theory/arith/nl_monomial.cpp create mode 100644 src/theory/arith/nl_monomial.h create mode 100644 src/theory/arith/nl_solver.cpp create mode 100644 src/theory/arith/nl_solver.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dec5a7f16..5284aa12c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -298,10 +298,16 @@ libcvc4_add_sources( theory/arith/linear_equality.h theory/arith/matrix.cpp theory/arith/matrix.h + theory/arith/nl_constraint.cpp + theory/arith/nl_constraint.h theory/arith/nl_lemma_utils.cpp theory/arith/nl_lemma_utils.h theory/arith/nl_model.cpp theory/arith/nl_model.h + theory/arith/nl_monomial.cpp + theory/arith/nl_monomial.h + theory/arith/nl_solver.cpp + theory/arith/nl_solver.h theory/arith/nonlinear_extension.cpp theory/arith/nonlinear_extension.h theory/arith/normal_form.cpp diff --git a/src/theory/arith/nl_constraint.cpp b/src/theory/arith/nl_constraint.cpp new file mode 100644 index 000000000..8fb4535ea --- /dev/null +++ b/src/theory/arith/nl_constraint.cpp @@ -0,0 +1,123 @@ +/********************* */ +/*! \file nl_constraint.cpp + ** \verbatim + ** Top contributors (to current version): + ** Andrew Reynolds + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2019 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 Implementation of utilities for non-linear constraints + **/ + +#include "theory/arith/nl_constraint.h" + +#include "theory/arith/arith_msum.h" +#include "theory/arith/arith_utilities.h" + +using namespace CVC4::kind; + +namespace CVC4 { +namespace theory { +namespace arith { + +ConstraintDb::ConstraintDb(MonomialDb& mdb) : d_mdb(mdb) {} + +void ConstraintDb::registerConstraint(Node atom) +{ + if (std::find(d_constraints.begin(), d_constraints.end(), atom) + != d_constraints.end()) + { + return; + } + d_constraints.push_back(atom); + Trace("nl-ext-debug") << "Register constraint : " << atom << std::endl; + std::map msum; + if (ArithMSum::getMonomialSumLit(atom, msum)) + { + Trace("nl-ext-debug") << "got monomial sum: " << std::endl; + if (Trace.isOn("nl-ext-debug")) + { + ArithMSum::debugPrintMonomialSum(msum, "nl-ext-debug"); + } + unsigned max_degree = 0; + std::vector all_m; + std::vector max_deg_m; + for (std::map::iterator itm = msum.begin(); itm != msum.end(); + ++itm) + { + if (!itm->first.isNull()) + { + all_m.push_back(itm->first); + d_mdb.registerMonomial(itm->first); + Trace("nl-ext-debug2") + << "...process monomial " << itm->first << std::endl; + unsigned d = d_mdb.getDegree(itm->first); + if (d > max_degree) + { + max_degree = d; + max_deg_m.clear(); + } + if (d >= max_degree) + { + max_deg_m.push_back(itm->first); + } + } + } + // isolate for each maximal degree monomial + for (unsigned i = 0; i < all_m.size(); i++) + { + Node m = all_m[i]; + Node rhs, coeff; + int res = ArithMSum::isolate(m, msum, coeff, rhs, atom.getKind()); + if (res != 0) + { + Kind type = atom.getKind(); + if (res == -1) + { + type = reverseRelationKind(type); + } + Trace("nl-ext-constraint") << "Constraint : " << atom << " <=> "; + if (!coeff.isNull()) + { + Trace("nl-ext-constraint") << coeff << " * "; + } + Trace("nl-ext-constraint") + << m << " " << type << " " << rhs << std::endl; + ConstraintInfo& ci = d_c_info[atom][m]; + ci.d_rhs = rhs; + ci.d_coeff = coeff; + ci.d_type = type; + } + } + for (unsigned i = 0; i < max_deg_m.size(); i++) + { + Node m = max_deg_m[i]; + d_c_info_maxm[atom][m] = true; + } + } + else + { + Trace("nl-ext-debug") << "...failed to get monomial sum." << std::endl; + } +} + +const std::map >& +ConstraintDb::getConstraints() +{ + return d_c_info; +} + +bool ConstraintDb::isMaximal(Node atom, Node x) const +{ + std::map >::const_iterator itcm = + d_c_info_maxm.find(atom); + Assert(itcm != d_c_info_maxm.end()); + return itcm->second.find(x) != itcm->second.end(); +} + +} // namespace arith +} // namespace theory +} // namespace CVC4 diff --git a/src/theory/arith/nl_constraint.h b/src/theory/arith/nl_constraint.h new file mode 100644 index 000000000..faa3cc812 --- /dev/null +++ b/src/theory/arith/nl_constraint.h @@ -0,0 +1,86 @@ +/********************* */ +/*! \file nl_constraint.h + ** \verbatim + ** Top contributors (to current version): + ** Andrew Reynolds, Tim King + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2019 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 non-linear constraints + **/ + +#ifndef CVC4__THEORY__ARITH__NL_CONSTRAINT_H +#define CVC4__THEORY__ARITH__NL_CONSTRAINT_H + +#include +#include + +#include "expr/kind.h" +#include "expr/node.h" +#include "theory/arith/nl_monomial.h" + +namespace CVC4 { +namespace theory { +namespace arith { + +/** constraint information + * + * The struct ( d_rhs, d_coeff, d_type ) represents that a literal is of the + * form (d_coeff * x) d_rhs. + */ +struct ConstraintInfo +{ + public: + /** The term on the right hand side of the constraint */ + Node d_rhs; + /** The coefficent */ + Node d_coeff; + /** The type (relation) of the constraint */ + Kind d_type; +}; /* struct ConstraintInfo */ + +/** A database for constraints */ +class ConstraintDb +{ + public: + ConstraintDb(MonomialDb& mdb); + ~ConstraintDb() {} + /** register constraint + * + * This ensures that atom is in the domain of the constraints maintained by + * this database. + */ + void registerConstraint(Node atom); + /** get constraints + * + * Returns a map m such that whenever + * m[lit][x] = ( r, coeff, k ), then + * ( lit <=> (coeff * x) r ) + */ + const std::map >& getConstraints(); + /** Returns true if m is of maximal degree in atom + * + * For example, for atom x^2 + x*y + y >=0, the monomials x^2 and x*y + * are of maximal degree (2). + */ + bool isMaximal(Node atom, Node m) const; + + private: + /** Reference to a monomial database */ + MonomialDb& d_mdb; + /** List of all constraints */ + std::vector d_constraints; + /** Is maximal degree */ + std::map > d_c_info_maxm; + /** Constraint information */ + std::map > d_c_info; +}; + +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif /* CVC4__THEORY__ARITH__NL_SOLVER_H */ diff --git a/src/theory/arith/nl_lemma_utils.h b/src/theory/arith/nl_lemma_utils.h index 9ad9f2ca5..bd729dad9 100644 --- a/src/theory/arith/nl_lemma_utils.h +++ b/src/theory/arith/nl_lemma_utils.h @@ -71,9 +71,9 @@ struct SortNlModel struct SortNonlinearDegree { - SortNonlinearDegree(std::map& m) : d_mdegree(m) {} + SortNonlinearDegree(const std::map& m) : d_mdegree(m) {} /** pointer to the non-linear extension */ - std::map& d_mdegree; + const std::map& d_mdegree; /** Get the degree of n in d_mdegree */ unsigned getDegree(Node n) const; /** diff --git a/src/theory/arith/nl_monomial.cpp b/src/theory/arith/nl_monomial.cpp new file mode 100644 index 000000000..be472217d --- /dev/null +++ b/src/theory/arith/nl_monomial.cpp @@ -0,0 +1,334 @@ +/********************* */ +/*! \file nl_monomial.cpp + ** \verbatim + ** Top contributors (to current version): + ** Andrew Reynolds + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2019 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 Implementation of utilities for monomials + **/ + +#include "theory/arith/nl_monomial.h" + +#include "theory/arith/arith_utilities.h" +#include "theory/arith/nl_lemma_utils.h" +#include "theory/rewriter.h" + +using namespace CVC4::kind; + +namespace CVC4 { +namespace theory { +namespace arith { + +// Returns a[key] if key is in a or value otherwise. +unsigned getCountWithDefault(const NodeMultiset& a, Node key, unsigned value) +{ + NodeMultiset::const_iterator it = a.find(key); + return (it == a.end()) ? value : it->second; +} +// Given two multisets return the multiset difference a \ b. +NodeMultiset diffMultiset(const NodeMultiset& a, const NodeMultiset& b) +{ + NodeMultiset difference; + for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) + { + Node key = it_a->first; + const unsigned a_value = it_a->second; + const unsigned b_value = getCountWithDefault(b, key, 0); + if (a_value > b_value) + { + difference[key] = a_value - b_value; + } + } + return difference; +} + +// Return a vector containing a[key] repetitions of key in a multiset a. +std::vector ExpandMultiset(const NodeMultiset& a) +{ + std::vector expansion; + for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) + { + expansion.insert(expansion.end(), it_a->second, it_a->first); + } + return expansion; +} + +// status 0 : n equal, -1 : n superset, 1 : n subset +void MonomialIndex::addTerm(Node n, + const std::vector& reps, + MonomialDb* nla, + int status, + unsigned argIndex) +{ + if (status == 0) + { + if (argIndex == reps.size()) + { + d_monos.push_back(n); + } + else + { + d_data[reps[argIndex]].addTerm(n, reps, nla, status, argIndex + 1); + } + } + for (std::map::iterator it = d_data.begin(); + it != d_data.end(); + ++it) + { + if (status != 0 || argIndex == reps.size() || it->first != reps[argIndex]) + { + // if we do not contain this variable, then if we were a superset, + // fail (-2), otherwise we are subset. if we do contain this + // variable, then if we were equal, we are superset since variables + // are ordered, otherwise we remain the same. + int new_status = + std::find(reps.begin(), reps.end(), it->first) == reps.end() + ? (status >= 0 ? 1 : -2) + : (status == 0 ? -1 : status); + if (new_status != -2) + { + it->second.addTerm(n, reps, nla, new_status, argIndex); + } + } + } + // compare for subsets + for (unsigned i = 0; i < d_monos.size(); i++) + { + Node m = d_monos[i]; + if (m != n) + { + // we are superset if we are equal and haven't traversed all variables + int cstatus = status == 0 ? (argIndex == reps.size() ? 0 : -1) : status; + Trace("nl-ext-mindex-debug") << " compare " << n << " and " << m + << ", status = " << cstatus << std::endl; + if (cstatus <= 0 && nla->isMonomialSubset(m, n)) + { + nla->registerMonomialSubset(m, n); + Trace("nl-ext-mindex-debug") << "...success" << std::endl; + } + else if (cstatus >= 0 && nla->isMonomialSubset(n, m)) + { + nla->registerMonomialSubset(n, m); + Trace("nl-ext-mindex-debug") << "...success (rev)" << std::endl; + } + } + } +} + +MonomialDb::MonomialDb() +{ + d_one = NodeManager::currentNM()->mkConst(Rational(1)); +} + +void MonomialDb::registerMonomial(Node n) +{ + if (std::find(d_monomials.begin(), d_monomials.end(), n) != d_monomials.end()) + { + return; + } + d_monomials.push_back(n); + Trace("nl-ext-debug") << "Register monomial : " << n << std::endl; + Kind k = n.getKind(); + if (k == NONLINEAR_MULT) + { + // get exponent count + unsigned nchild = n.getNumChildren(); + for (unsigned i = 0; i < nchild; i++) + { + d_m_exp[n][n[i]]++; + if (i == 0 || n[i] != n[i - 1]) + { + d_m_vlist[n].push_back(n[i]); + } + } + d_m_degree[n] = nchild; + } + else if (n == d_one) + { + d_m_exp[n].clear(); + d_m_vlist[n].clear(); + d_m_degree[n] = 0; + } + else + { + Assert(k != PLUS && k != MULT); + d_m_exp[n][n] = 1; + d_m_vlist[n].push_back(n); + d_m_degree[n] = 1; + } + std::sort(d_m_vlist[n].begin(), d_m_vlist[n].end()); + Trace("nl-ext-mindex") << "Add monomial to index : " << n << std::endl; + d_m_index.addTerm(n, d_m_vlist[n], this); +} + +void MonomialDb::registerMonomialSubset(Node a, Node b) +{ + Assert(isMonomialSubset(a, b)); + + const NodeMultiset& a_exponent_map = getMonomialExponentMap(a); + const NodeMultiset& b_exponent_map = getMonomialExponentMap(b); + + std::vector diff_children = + ExpandMultiset(diffMultiset(b_exponent_map, a_exponent_map)); + Assert(!diff_children.empty()); + + d_m_contain_parent[a].push_back(b); + d_m_contain_children[b].push_back(a); + + Node mult_term = safeConstructNary(MULT, diff_children); + Node nlmult_term = safeConstructNary(NONLINEAR_MULT, diff_children); + d_m_contain_mult[a][b] = mult_term; + d_m_contain_umult[a][b] = nlmult_term; + Trace("nl-ext-mindex") << "..." << a << " is a subset of " << b + << ", difference is " << mult_term << std::endl; +} + +bool MonomialDb::isMonomialSubset(Node am, Node bm) const +{ + const NodeMultiset& a = getMonomialExponentMap(am); + const NodeMultiset& b = getMonomialExponentMap(bm); + for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) + { + Node key = it_a->first; + const unsigned a_value = it_a->second; + const unsigned b_value = getCountWithDefault(b, key, 0); + if (a_value > b_value) + { + return false; + } + } + return true; +} + +const NodeMultiset& MonomialDb::getMonomialExponentMap(Node monomial) const +{ + MonomialExponentMap::const_iterator it = d_m_exp.find(monomial); + Assert(it != d_m_exp.end()); + return it->second; +} + +unsigned MonomialDb::getExponent(Node monomial, Node v) const +{ + MonomialExponentMap::const_iterator it = d_m_exp.find(monomial); + if (it == d_m_exp.end()) + { + return 0; + } + std::map::const_iterator itv = it->second.find(v); + if (itv == it->second.end()) + { + return 0; + } + return itv->second; +} + +const std::vector& MonomialDb::getVariableList(Node monomial) const +{ + std::map >::const_iterator itvl = + d_m_vlist.find(monomial); + Assert(itvl != d_m_vlist.end()); + return itvl->second; +} + +unsigned MonomialDb::getDegree(Node monomial) const +{ + std::map::const_iterator it = d_m_degree.find(monomial); + Assert(it != d_m_degree.end()); + return it->second; +} + +void MonomialDb::sortByDegree(std::vector& ms) const +{ + SortNonlinearDegree snlad(d_m_degree); + std::sort(ms.begin(), ms.end(), snlad); +} + +void MonomialDb::sortVariablesByModel(std::vector& ms, NlModel& m) +{ + SortNlModel smv; + smv.d_nlm = &m; + smv.d_isConcrete = false; + smv.d_isAbsolute = true; + smv.d_reverse_order = true; + for (const Node& msc : ms) + { + std::sort(d_m_vlist[msc].begin(), d_m_vlist[msc].end(), smv); + } +} + +const std::map >& MonomialDb::getContainsChildrenMap() +{ + return d_m_contain_children; +} + +const std::map >& MonomialDb::getContainsParentMap() +{ + return d_m_contain_parent; +} + +Node MonomialDb::getContainsDiff(Node a, Node b) const +{ + std::map >::const_iterator it = + d_m_contain_mult.find(a); + if (it == d_m_contain_umult.end()) + { + return Node::null(); + } + std::map::const_iterator it2 = it->second.find(b); + if (it2 == it->second.end()) + { + return Node::null(); + } + return it2->second; +} + +Node MonomialDb::getContainsDiffNl(Node a, Node b) const +{ + std::map >::const_iterator it = + d_m_contain_umult.find(a); + if (it == d_m_contain_umult.end()) + { + return Node::null(); + } + std::map::const_iterator it2 = it->second.find(b); + if (it2 == it->second.end()) + { + return Node::null(); + } + return it2->second; +} + +Node MonomialDb::mkMonomialRemFactor(Node n, + const NodeMultiset& n_exp_rem) const +{ + std::vector children; + const NodeMultiset& exponent_map = getMonomialExponentMap(n); + for (NodeMultiset::const_iterator itme2 = exponent_map.begin(); + itme2 != exponent_map.end(); + ++itme2) + { + Node v = itme2->first; + unsigned inc = itme2->second; + Trace("nl-ext-mono-factor") + << "..." << inc << " factors of " << v << std::endl; + unsigned count_in_n_exp_rem = getCountWithDefault(n_exp_rem, v, 0); + Assert(count_in_n_exp_rem <= inc); + inc -= count_in_n_exp_rem; + Trace("nl-ext-mono-factor") + << "......rem, now " << inc << " factors of " << v << std::endl; + children.insert(children.end(), inc, v); + } + Node ret = safeConstructNary(MULT, children); + ret = Rewriter::rewrite(ret); + Trace("nl-ext-mono-factor") << "...return : " << ret << std::endl; + return ret; +} + +} // namespace arith +} // namespace theory +} // namespace CVC4 diff --git a/src/theory/arith/nl_monomial.h b/src/theory/arith/nl_monomial.h new file mode 100644 index 000000000..81665c4d9 --- /dev/null +++ b/src/theory/arith/nl_monomial.h @@ -0,0 +1,147 @@ +/********************* */ +/*! \file nl_monomial.h + ** \verbatim + ** Top contributors (to current version): + ** Andrew Reynolds, Tim King + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2019 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 monomials + **/ + +#ifndef CVC4__THEORY__ARITH__NL_MONOMIAL_H +#define CVC4__THEORY__ARITH__NL_MONOMIAL_H + +#include +#include + +#include "expr/node.h" + +namespace CVC4 { +namespace theory { +namespace arith { + +class MonomialDb; +class NlModel; + +typedef std::map NodeMultiset; +typedef std::map MonomialExponentMap; + +/** An index data structure for node multisets (monomials) */ +class MonomialIndex +{ + public: + /** + * Add term to this trie. The argument status indicates what the status + * of n is with respect to the current node in the trie, where: + * 0 : n is equal, -1 : n is superset, 1 : n is subset + * of the node described by the current path in the trie. + */ + void addTerm(Node n, + const std::vector& reps, + MonomialDb* nla, + int status = 0, + unsigned argIndex = 0); + + private: + /** The children of this node */ + std::map d_data; + /** The monomials at this node */ + std::vector d_monos; +}; /* class MonomialIndex */ + +/** Context-independent database for monomial information */ +class MonomialDb +{ + public: + MonomialDb(); + ~MonomialDb() {} + /** register monomial */ + void registerMonomial(Node n); + /** + * Register monomial subset. This method is called when we infer that b is + * a subset of monomial a, e.g. x*y^2 is a subset of x^3*y^2*z. + */ + void registerMonomialSubset(Node a, Node b); + /** + * returns true if the multiset containing the + * factors of monomial a is a subset of the multiset + * containing the factors of monomial b. + */ + bool isMonomialSubset(Node a, Node b) const; + /** Returns the NodeMultiset for a registered monomial. */ + const NodeMultiset& getMonomialExponentMap(Node monomial) const; + /** Returns the exponent of variable v in the given monomial */ + unsigned getExponent(Node monomial, Node v) const; + /** Get the list of unique variables is the monomial */ + const std::vector& getVariableList(Node monomial) const; + /** Get degree of monomial, e.g. the degree of x^2*y^2 = 4 */ + unsigned getDegree(Node monomial) const; + /** Sort monomials in ms by their degree + * + * Updates ms so that degree(ms[i]) <= degree(ms[j]) for i <= j. + */ + void sortByDegree(std::vector& ms) const; + /** Sort the variable lists based on model values + * + * This updates the variable lists of monomials in ms based on the absolute + * value of their current model values in m. + * + * In other words, for each i, getVariableList(ms[i]) returns + * v1, ..., vn where |m(v1)| <= ... <= |m(vn)| after this method is invoked. + */ + void sortVariablesByModel(std::vector& ms, NlModel& m); + /** Get monomial contains children map + * + * This maps monomials to other monomials that are contained in them, e.g. + * x^2 * y may map to { x, x^2, y } if these three terms exist have been + * registered to this class. + */ + const std::map >& getContainsChildrenMap(); + /** Get monomial contains parent map, reverse of above */ + const std::map >& getContainsParentMap(); + /** + * Get contains difference. Return the difference of a and b or null if it + * does not exist. In other words, this returns a term equivalent to a/b + * that does not contain division. + */ + Node getContainsDiff(Node a, Node b) const; + /** + * Get contains difference non-linear. Same as above, but stores terms of kind + * NONLINEAR_MULT instead of MULT. + */ + Node getContainsDiffNl(Node a, Node b) const; + /** Make monomial remainder factor */ + Node mkMonomialRemFactor(Node n, const NodeMultiset& n_exp_rem) const; + + private: + /** commonly used terms */ + Node d_one; + /** list of all monomials */ + std::vector d_monomials; + /** Map from monomials to var^index. */ + MonomialExponentMap d_m_exp; + /** + * Mapping from monomials to the list of variables that occur in it. For + * example, x*x*y*z -> { x, y, z }. + */ + std::map > d_m_vlist; + /** Degree information */ + std::map d_m_degree; + /** monomial index, by sorted variables */ + MonomialIndex d_m_index; + /** containment ordering */ + std::map > d_m_contain_children; + std::map > d_m_contain_parent; + std::map > d_m_contain_mult; + std::map > d_m_contain_umult; +}; + +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif /* CVC4__THEORY__ARITH__NL_MONOMIAL_H */ diff --git a/src/theory/arith/nl_solver.cpp b/src/theory/arith/nl_solver.cpp new file mode 100644 index 000000000..893d3dbd7 --- /dev/null +++ b/src/theory/arith/nl_solver.cpp @@ -0,0 +1,1585 @@ +/********************* */ +/*! \file nl_solver.cpp + ** \verbatim + ** Top contributors (to current version): + ** Andrew Reynolds + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2019 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 Implementation of non-linear solver + **/ + +#include "theory/arith/nl_solver.h" + +#include "options/arith_options.h" +#include "theory/arith/arith_msum.h" +#include "theory/arith/arith_utilities.h" +#include "theory/arith/theory_arith.h" +#include "theory/theory_model.h" + +using namespace CVC4::kind; + +namespace CVC4 { +namespace theory { +namespace arith { + +void debugPrintBound(const char* c, Node coeff, Node x, Kind type, Node rhs) +{ + Node t = ArithMSum::mkCoeffTerm(coeff, x); + Trace(c) << t << " " << type << " " << rhs; +} + +bool hasNewMonomials(Node n, const std::vector& existing) +{ + std::set visited; + + std::vector worklist; + worklist.push_back(n); + while (!worklist.empty()) + { + Node current = worklist.back(); + worklist.pop_back(); + if (visited.find(current) == visited.end()) + { + visited.insert(current); + if (current.getKind() == NONLINEAR_MULT) + { + if (std::find(existing.begin(), existing.end(), current) + == existing.end()) + { + return true; + } + } + else + { + worklist.insert(worklist.end(), current.begin(), current.end()); + } + } + } + return false; +} + +NlSolver::NlSolver(TheoryArith& containing, NlModel& model) + : d_containing(containing), + d_model(model), + d_cdb(d_mdb), + d_zero_split(containing.getUserContext()) +{ + NodeManager* nm = NodeManager::currentNM(); + d_true = nm->mkConst(true); + d_false = nm->mkConst(false); + d_zero = nm->mkConst(Rational(0)); + d_one = nm->mkConst(Rational(1)); + d_neg_one = nm->mkConst(Rational(-1)); + d_order_points.push_back(d_neg_one); + d_order_points.push_back(d_zero); + d_order_points.push_back(d_one); +} + +NlSolver::~NlSolver() {} + +void NlSolver::initLastCall(const std::vector& assertions, + const std::vector& false_asserts, + const std::vector& xts) +{ + d_ms_vars.clear(); + d_ms_proc.clear(); + d_ms.clear(); + d_mterms.clear(); + d_m_nconst_factor.clear(); + d_tplane_refine.clear(); + d_ci.clear(); + d_ci_exp.clear(); + d_ci_max.clear(); + + Trace("nl-ext-mv") << "Extended terms : " << std::endl; + // for computing congruence + std::map argTrie; + for (unsigned i = 0, xsize = xts.size(); i < xsize; i++) + { + Node a = xts[i]; + d_model.computeConcreteModelValue(a); + d_model.computeAbstractModelValue(a); + d_model.printModelValue("nl-ext-mv", a); + Kind ak = a.getKind(); + if (ak == NONLINEAR_MULT) + { + d_ms.push_back(a); + + // context-independent registration + d_mdb.registerMonomial(a); + + const std::vector& varList = d_mdb.getVariableList(a); + for (const Node& v : varList) + { + if (std::find(d_ms_vars.begin(), d_ms_vars.end(), v) == d_ms_vars.end()) + { + d_ms_vars.push_back(v); + } + Node mvk = d_model.computeAbstractModelValue(v); + if (!mvk.isConst()) + { + d_m_nconst_factor[a] = true; + } + } + // mark processed if has a "one" factor (will look at reduced monomial)? + } + } + + // register constants + d_mdb.registerMonomial(d_one); + for (unsigned j = 0; j < d_order_points.size(); j++) + { + Node c = d_order_points[j]; + d_model.computeConcreteModelValue(c); + d_model.computeAbstractModelValue(c); + } + + // register variables + Trace("nl-ext-mv") << "Variables in monomials : " << std::endl; + for (unsigned i = 0; i < d_ms_vars.size(); i++) + { + Node v = d_ms_vars[i]; + d_mdb.registerMonomial(v); + d_model.computeConcreteModelValue(v); + d_model.computeAbstractModelValue(v); + d_model.printModelValue("nl-ext-mv", v); + } + + Trace("nl-ext") << "We have " << d_ms.size() << " monomials." << std::endl; +} + +void NlSolver::setMonomialFactor(Node a, Node b, const NodeMultiset& common) +{ + // Could not tell if this was being inserted intentionally or not. + std::map& mono_diff_a = d_mono_diff[a]; + if (mono_diff_a.find(b) == mono_diff_a.end()) + { + Trace("nl-ext-mono-factor") + << "Set monomial factor for " << a << "/" << b << std::endl; + mono_diff_a[b] = d_mdb.mkMonomialRemFactor(a, common); + } +} + +std::vector NlSolver::checkSplitZero() +{ + std::vector lemmas; + for (unsigned i = 0; i < d_ms_vars.size(); i++) + { + Node v = d_ms_vars[i]; + if (d_zero_split.insert(v)) + { + Node eq = v.eqNode(d_zero); + eq = Rewriter::rewrite(eq); + Node literal = d_containing.getValuation().ensureLiteral(eq); + d_containing.getOutputChannel().requirePhase(literal, true); + lemmas.push_back(literal.orNode(literal.negate())); + } + } + return lemmas; +} + +void NlSolver::assignOrderIds(std::vector& vars, + NodeMultiset& order, + bool isConcrete, + bool isAbsolute) +{ + SortNlModel smv; + smv.d_nlm = &d_model; + smv.d_isConcrete = isConcrete; + smv.d_isAbsolute = isAbsolute; + smv.d_reverse_order = false; + std::sort(vars.begin(), vars.end(), smv); + + order.clear(); + // assign ordering id's + unsigned counter = 0; + unsigned order_index = isConcrete ? 0 : 1; + Node prev; + for (unsigned j = 0; j < vars.size(); j++) + { + Node x = vars[j]; + Node v = d_model.computeModelValue(x, isConcrete); + if (!v.isConst()) + { + Trace("nl-ext-mvo") << "..do not assign order to " << x << " : " << v + << std::endl; + // don't assign for non-constant values (transcendental function apps) + break; + } + Trace("nl-ext-mvo") << " order " << x << " : " << v << std::endl; + if (v != prev) + { + // builtin points + bool success; + do + { + success = false; + if (order_index < d_order_points.size()) + { + Node vv = d_model.computeModelValue(d_order_points[order_index], + isConcrete); + if (d_model.compareValue(v, vv, isAbsolute) <= 0) + { + counter++; + Trace("nl-ext-mvo") << "O[" << d_order_points[order_index] + << "] = " << counter << std::endl; + order[d_order_points[order_index]] = counter; + prev = vv; + order_index++; + success = true; + } + } + } while (success); + } + if (prev.isNull() || d_model.compareValue(v, prev, isAbsolute) != 0) + { + counter++; + } + Trace("nl-ext-mvo") << "O[" << x << "] = " << counter << std::endl; + order[x] = counter; + prev = v; + } + while (order_index < d_order_points.size()) + { + counter++; + Trace("nl-ext-mvo") << "O[" << d_order_points[order_index] + << "] = " << counter << std::endl; + order[d_order_points[order_index]] = counter; + order_index++; + } +} + +// show a <> 0 by inequalities between variables in monomial a w.r.t 0 +int NlSolver::compareSign(Node oa, + Node a, + unsigned a_index, + int status, + std::vector& exp, + std::vector& lem) +{ + Trace("nl-ext-debug") << "Process " << a << " at index " << a_index + << ", status is " << status << std::endl; + NodeManager* nm = NodeManager::currentNM(); + Node mvaoa = d_model.computeAbstractModelValue(oa); + const std::vector& vla = d_mdb.getVariableList(a); + if (a_index == vla.size()) + { + if (mvaoa.getConst().sgn() != status) + { + Node lemma = + safeConstructNary(AND, exp).impNode(mkLit(oa, d_zero, status * 2)); + lem.push_back(lemma); + } + return status; + } + Assert(a_index < vla.size()); + Node av = vla[a_index]; + unsigned aexp = d_mdb.getExponent(a, av); + // take current sign in model + Node mvaav = d_model.computeAbstractModelValue(av); + int sgn = mvaav.getConst().sgn(); + Trace("nl-ext-debug") << "Process var " << av << "^" << aexp + << ", model sign = " << sgn << std::endl; + if (sgn == 0) + { + if (mvaoa.getConst().sgn() != 0) + { + Node lemma = av.eqNode(d_zero).impNode(oa.eqNode(d_zero)); + lem.push_back(lemma); + } + return 0; + } + if (aexp % 2 == 0) + { + exp.push_back(av.eqNode(d_zero).negate()); + return compareSign(oa, a, a_index + 1, status, exp, lem); + } + exp.push_back(nm->mkNode(sgn == 1 ? GT : LT, av, d_zero)); + return compareSign(oa, a, a_index + 1, status * sgn, exp, lem); +} + +bool NlSolver::compareMonomial( + Node oa, + Node a, + NodeMultiset& a_exp_proc, + Node ob, + Node b, + NodeMultiset& b_exp_proc, + std::vector& exp, + std::vector& lem, + std::map > >& cmp_infers) +{ + Trace("nl-ext-comp-debug") + << "Check |" << a << "| >= |" << b << "|" << std::endl; + unsigned pexp_size = exp.size(); + if (compareMonomial( + oa, a, 0, a_exp_proc, ob, b, 0, b_exp_proc, 0, exp, lem, cmp_infers)) + { + return true; + } + exp.resize(pexp_size); + Trace("nl-ext-comp-debug") + << "Check |" << b << "| >= |" << a << "|" << std::endl; + if (compareMonomial( + ob, b, 0, b_exp_proc, oa, a, 0, a_exp_proc, 0, exp, lem, cmp_infers)) + { + return true; + } + return false; +} + +Node NlSolver::mkLit(Node a, Node b, int status, bool isAbsolute) +{ + if (status == 0) + { + Node a_eq_b = a.eqNode(b); + if (!isAbsolute) + { + return a_eq_b; + } + Node negate_b = NodeManager::currentNM()->mkNode(UMINUS, b); + return a_eq_b.orNode(a.eqNode(negate_b)); + } + else if (status < 0) + { + return mkLit(b, a, -status); + } + Assert(status == 1 || status == 2); + NodeManager* nm = NodeManager::currentNM(); + Kind greater_op = status == 1 ? GEQ : GT; + if (!isAbsolute) + { + return nm->mkNode(greater_op, a, b); + } + // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) ); + Node zero = mkRationalNode(0); + Node a_is_nonnegative = nm->mkNode(GEQ, a, zero); + Node b_is_nonnegative = nm->mkNode(GEQ, b, zero); + Node negate_a = nm->mkNode(UMINUS, a); + Node negate_b = nm->mkNode(UMINUS, b); + return a_is_nonnegative.iteNode( + b_is_nonnegative.iteNode(nm->mkNode(greater_op, a, b), + nm->mkNode(greater_op, a, negate_b)), + b_is_nonnegative.iteNode(nm->mkNode(greater_op, negate_a, b), + nm->mkNode(greater_op, negate_a, negate_b))); +} + +bool NlSolver::cmp_holds(Node x, + Node y, + std::map >& cmp_infers, + std::vector& exp, + std::map& visited) +{ + if (x == y) + { + return true; + } + else if (visited.find(x) != visited.end()) + { + return false; + } + visited[x] = true; + std::map >::iterator it = cmp_infers.find(x); + if (it != cmp_infers.end()) + { + for (std::map::iterator itc = it->second.begin(); + itc != it->second.end(); + ++itc) + { + exp.push_back(itc->second); + if (cmp_holds(itc->first, y, cmp_infers, exp, visited)) + { + return true; + } + exp.pop_back(); + } + } + return false; +} + +// trying to show a ( >, = ) b by inequalities between variables in +// monomials a,b +bool NlSolver::compareMonomial( + Node oa, + Node a, + unsigned a_index, + NodeMultiset& a_exp_proc, + Node ob, + Node b, + unsigned b_index, + NodeMultiset& b_exp_proc, + int status, + std::vector& exp, + std::vector& lem, + std::map > >& cmp_infers) +{ + Trace("nl-ext-comp-debug") + << "compareMonomial " << oa << " and " << ob << ", indices = " << a_index + << " " << b_index << std::endl; + Assert(status == 0 || status == 2); + const std::vector& vla = d_mdb.getVariableList(a); + const std::vector& vlb = d_mdb.getVariableList(b); + if (a_index == vla.size() && b_index == vlb.size()) + { + // finished, compare absolute value of abstract model values + int modelStatus = d_model.compare(oa, ob, false, true) * -2; + Trace("nl-ext-comp") << "...finished comparison with " << oa << " <" + << status << "> " << ob + << ", model status = " << modelStatus << std::endl; + if (status != modelStatus) + { + Trace("nl-ext-comp-infer") + << "infer : " << oa << " <" << status << "> " << ob << std::endl; + if (status == 2) + { + // must state that all variables are non-zero + for (unsigned j = 0; j < vla.size(); j++) + { + exp.push_back(vla[j].eqNode(d_zero).negate()); + } + } + NodeManager* nm = NodeManager::currentNM(); + Node clem = nm->mkNode( + IMPLIES, safeConstructNary(AND, exp), mkLit(oa, ob, status, true)); + Trace("nl-ext-comp-lemma") << "comparison lemma : " << clem << std::endl; + lem.push_back(clem); + cmp_infers[status][oa][ob] = clem; + } + return true; + } + // get a/b variable information + Node av; + unsigned aexp = 0; + unsigned avo = 0; + if (a_index < vla.size()) + { + av = vla[a_index]; + unsigned aexpTotal = d_mdb.getExponent(a, av); + Assert(a_exp_proc[av] <= aexpTotal); + aexp = aexpTotal - a_exp_proc[av]; + if (aexp == 0) + { + return compareMonomial(oa, + a, + a_index + 1, + a_exp_proc, + ob, + b, + b_index, + b_exp_proc, + status, + exp, + lem, + cmp_infers); + } + Assert(d_order_vars.find(av) != d_order_vars.end()); + avo = d_order_vars[av]; + } + Node bv; + unsigned bexp = 0; + unsigned bvo = 0; + if (b_index < vlb.size()) + { + bv = vlb[b_index]; + unsigned bexpTotal = d_mdb.getExponent(b, bv); + Assert(b_exp_proc[bv] <= bexpTotal); + bexp = bexpTotal - b_exp_proc[bv]; + if (bexp == 0) + { + return compareMonomial(oa, + a, + a_index, + a_exp_proc, + ob, + b, + b_index + 1, + b_exp_proc, + status, + exp, + lem, + cmp_infers); + } + Assert(d_order_vars.find(bv) != d_order_vars.end()); + bvo = d_order_vars[bv]; + } + // get "one" information + Assert(d_order_vars.find(d_one) != d_order_vars.end()); + unsigned ovo = d_order_vars[d_one]; + Trace("nl-ext-comp-debug") << "....vars : " << av << "^" << aexp << " " << bv + << "^" << bexp << std::endl; + + //--- cases + if (av.isNull()) + { + if (bvo <= ovo) + { + Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl; + // can multiply b by <=1 + exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true)); + return compareMonomial(oa, + a, + a_index, + a_exp_proc, + ob, + b, + b_index + 1, + b_exp_proc, + bvo == ovo ? status : 2, + exp, + lem, + cmp_infers); + } + Trace("nl-ext-comp-debug") + << "...failure, unmatched |b|>1 component." << std::endl; + return false; + } + else if (bv.isNull()) + { + if (avo >= ovo) + { + Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl; + // can multiply a by >=1 + exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true)); + return compareMonomial(oa, + a, + a_index + 1, + a_exp_proc, + ob, + b, + b_index, + b_exp_proc, + avo == ovo ? status : 2, + exp, + lem, + cmp_infers); + } + Trace("nl-ext-comp-debug") + << "...failure, unmatched |a|<1 component." << std::endl; + return false; + } + Assert(!av.isNull() && !bv.isNull()); + if (avo >= bvo) + { + if (bvo < ovo && avo >= ovo) + { + Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl; + // do avo>=1 instead + exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true)); + return compareMonomial(oa, + a, + a_index + 1, + a_exp_proc, + ob, + b, + b_index, + b_exp_proc, + avo == ovo ? status : 2, + exp, + lem, + cmp_infers); + } + unsigned min_exp = aexp > bexp ? bexp : aexp; + a_exp_proc[av] += min_exp; + b_exp_proc[bv] += min_exp; + Trace("nl-ext-comp-debug") << "...take leading " << min_exp << " from " + << av << " and " << bv << std::endl; + exp.push_back(mkLit(av, bv, avo == bvo ? 0 : 2, true)); + bool ret = compareMonomial(oa, + a, + a_index, + a_exp_proc, + ob, + b, + b_index, + b_exp_proc, + avo == bvo ? status : 2, + exp, + lem, + cmp_infers); + a_exp_proc[av] -= min_exp; + b_exp_proc[bv] -= min_exp; + return ret; + } + if (bvo <= ovo) + { + Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl; + // try multiply b <= 1 + exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true)); + return compareMonomial(oa, + a, + a_index, + a_exp_proc, + ob, + b, + b_index + 1, + b_exp_proc, + bvo == ovo ? status : 2, + exp, + lem, + cmp_infers); + } + Trace("nl-ext-comp-debug") + << "...failure, leading |b|>|a|>1 component." << std::endl; + return false; +} + +std::vector NlSolver::checkMonomialSign() +{ + std::vector lemmas; + std::map signs; + Trace("nl-ext") << "Get monomial sign lemmas..." << std::endl; + for (unsigned j = 0; j < d_ms.size(); j++) + { + Node a = d_ms[j]; + if (d_ms_proc.find(a) == d_ms_proc.end()) + { + std::vector exp; + if (Trace.isOn("nl-ext-debug")) + { + Node cmva = d_model.computeConcreteModelValue(a); + Trace("nl-ext-debug") + << " process " << a << ", mv=" << cmva << "..." << std::endl; + } + if (d_m_nconst_factor.find(a) == d_m_nconst_factor.end()) + { + signs[a] = compareSign(a, a, 0, 1, exp, lemmas); + if (signs[a] == 0) + { + d_ms_proc[a] = true; + Trace("nl-ext-debug") + << "...mark " << a << " reduced since its value is 0." + << std::endl; + } + } + else + { + Trace("nl-ext-debug") + << "...can't conclude sign lemma for " << a + << " since model value of a factor is non-constant." << std::endl; + } + } + } + return lemmas; +} + +std::vector NlSolver::checkMonomialMagnitude(unsigned c) +{ + // ensure information is setup + if (c == 0) + { + // sort by absolute values of abstract model values + assignOrderIds(d_ms_vars, d_order_vars, false, true); + + // sort individual variable lists + Trace("nl-ext-proc") << "Assign order var lists..." << std::endl; + d_mdb.sortVariablesByModel(d_ms, d_model); + } + + unsigned r = 1; + std::vector lemmas; + // if (x,y,L) in cmp_infers, then x > y inferred as conclusion of L + // in lemmas + std::map > > cmp_infers; + Trace("nl-ext") << "Get monomial comparison lemmas (order=" << r + << ", compare=" << c << ")..." << std::endl; + for (unsigned j = 0; j < d_ms.size(); j++) + { + Node a = d_ms[j]; + if (d_ms_proc.find(a) == d_ms_proc.end() + && d_m_nconst_factor.find(a) == d_m_nconst_factor.end()) + { + if (c == 0) + { + // compare magnitude against 1 + std::vector exp; + NodeMultiset a_exp_proc; + NodeMultiset b_exp_proc; + compareMonomial(a, + a, + a_exp_proc, + d_one, + d_one, + b_exp_proc, + exp, + lemmas, + cmp_infers); + } + else + { + const NodeMultiset& mea = d_mdb.getMonomialExponentMap(a); + if (c == 1) + { + // could compare not just against containing variables? + // compare magnitude against variables + for (unsigned k = 0; k < d_ms_vars.size(); k++) + { + Node v = d_ms_vars[k]; + Node mvcv = d_model.computeConcreteModelValue(v); + if (mvcv.isConst()) + { + std::vector exp; + NodeMultiset a_exp_proc; + NodeMultiset b_exp_proc; + if (mea.find(v) != mea.end()) + { + a_exp_proc[v] = 1; + b_exp_proc[v] = 1; + setMonomialFactor(a, v, a_exp_proc); + setMonomialFactor(v, a, b_exp_proc); + compareMonomial(a, + a, + a_exp_proc, + v, + v, + b_exp_proc, + exp, + lemmas, + cmp_infers); + } + } + } + } + else + { + // compare magnitude against other non-linear monomials + for (unsigned k = (j + 1); k < d_ms.size(); k++) + { + Node b = d_ms[k]; + //(signs[a]==signs[b])==(r==0) + if (d_ms_proc.find(b) == d_ms_proc.end() + && d_m_nconst_factor.find(b) == d_m_nconst_factor.end()) + { + const NodeMultiset& meb = d_mdb.getMonomialExponentMap(b); + + std::vector exp; + // take common factors of monomials, set minimum of + // common exponents as processed + NodeMultiset a_exp_proc; + NodeMultiset b_exp_proc; + for (NodeMultiset::const_iterator itmea2 = mea.begin(); + itmea2 != mea.end(); + ++itmea2) + { + NodeMultiset::const_iterator itmeb2 = meb.find(itmea2->first); + if (itmeb2 != meb.end()) + { + unsigned min_exp = itmea2->second > itmeb2->second + ? itmeb2->second + : itmea2->second; + a_exp_proc[itmea2->first] = min_exp; + b_exp_proc[itmea2->first] = min_exp; + Trace("nl-ext-comp") << "Common exponent : " << itmea2->first + << " : " << min_exp << std::endl; + } + } + if (!a_exp_proc.empty()) + { + setMonomialFactor(a, b, a_exp_proc); + setMonomialFactor(b, a, b_exp_proc); + } + /* + if( !a_exp_proc.empty() ){ + //reduction based on common exponents a > 0 => ( a * b + <> a * c <=> b <> c ), a < 0 => ( a * b <> a * c <=> b + !<> c ) ? }else{ compareMonomial( a, a, a_exp_proc, b, + b, b_exp_proc, exp, lemmas ); + } + */ + compareMonomial( + a, a, a_exp_proc, b, b, b_exp_proc, exp, lemmas, cmp_infers); + } + } + } + } + } + } + // remove redundant lemmas, e.g. if a > b, b > c, a > c were + // inferred, discard lemma with conclusion a > c + Trace("nl-ext-comp") << "Compute redundancies for " << lemmas.size() + << " lemmas." << std::endl; + // naive + std::vector r_lemmas; + for (std::map > >::iterator itb = + cmp_infers.begin(); + itb != cmp_infers.end(); + ++itb) + { + for (std::map >::iterator itc = + itb->second.begin(); + itc != itb->second.end(); + ++itc) + { + for (std::map::iterator itc2 = itc->second.begin(); + itc2 != itc->second.end(); + ++itc2) + { + std::map visited; + for (std::map::iterator itc3 = itc->second.begin(); + itc3 != itc->second.end(); + ++itc3) + { + if (itc3->first != itc2->first) + { + std::vector exp; + if (cmp_holds(itc3->first, itc2->first, itb->second, exp, visited)) + { + r_lemmas.push_back(itc2->second); + Trace("nl-ext-comp") + << "...inference of " << itc->first << " > " << itc2->first + << " was redundant." << std::endl; + break; + } + } + } + } + } + } + std::vector nr_lemmas; + for (unsigned i = 0; i < lemmas.size(); i++) + { + if (std::find(r_lemmas.begin(), r_lemmas.end(), lemmas[i]) + == r_lemmas.end()) + { + nr_lemmas.push_back(lemmas[i]); + } + } + // could only take maximal lower/minimial lower bounds? + + Trace("nl-ext-comp") << nr_lemmas.size() << " / " << lemmas.size() + << " were non-redundant." << std::endl; + return nr_lemmas; +} + +std::vector NlSolver::checkTangentPlanes() +{ + std::vector lemmas; + Trace("nl-ext") << "Get monomial tangent plane lemmas..." << std::endl; + NodeManager* nm = NodeManager::currentNM(); + const std::map >& ccMap = + d_mdb.getContainsChildrenMap(); + unsigned kstart = d_ms_vars.size(); + for (unsigned k = kstart; k < d_mterms.size(); k++) + { + Node t = d_mterms[k]; + // if this term requires a refinement + if (d_tplane_refine.find(t) == d_tplane_refine.end()) + { + continue; + } + Trace("nl-ext-tplanes") + << "Look at monomial requiring refinement : " << t << std::endl; + // get a decomposition + std::map >::const_iterator it = ccMap.find(t); + if (it == ccMap.end()) + { + continue; + } + std::map > dproc; + for (unsigned j = 0; j < it->second.size(); j++) + { + Node tc = it->second[j]; + if (tc != d_one) + { + Node tc_diff = d_mdb.getContainsDiffNl(tc, t); + Assert(!tc_diff.isNull()); + Node a = tc < tc_diff ? tc : tc_diff; + Node b = tc < tc_diff ? tc_diff : tc; + if (dproc[a].find(b) == dproc[a].end()) + { + dproc[a][b] = true; + Trace("nl-ext-tplanes") + << " decomposable into : " << a << " * " << b << std::endl; + Node a_v_c = d_model.computeAbstractModelValue(a); + Node b_v_c = d_model.computeAbstractModelValue(b); + // points we will add tangent planes for + std::vector pts[2]; + pts[0].push_back(a_v_c); + pts[1].push_back(b_v_c); + // if previously refined + bool prevRefine = d_tangent_val_bound[0][a].find(b) + != d_tangent_val_bound[0][a].end(); + // a_min, a_max, b_min, b_max + for (unsigned p = 0; p < 4; p++) + { + Node curr_v = p <= 1 ? a_v_c : b_v_c; + if (prevRefine) + { + Node pt_v = d_tangent_val_bound[p][a][b]; + Assert(!pt_v.isNull()); + if (curr_v != pt_v) + { + Node do_extend = + nm->mkNode((p == 1 || p == 3) ? GT : LT, curr_v, pt_v); + do_extend = Rewriter::rewrite(do_extend); + if (do_extend == d_true) + { + for (unsigned q = 0; q < 2; q++) + { + pts[p <= 1 ? 0 : 1].push_back(curr_v); + pts[p <= 1 ? 1 : 0].push_back( + d_tangent_val_bound[p <= 1 ? 2 + q : q][a][b]); + } + } + } + } + else + { + d_tangent_val_bound[p][a][b] = curr_v; + } + } + + for (unsigned p = 0; p < pts[0].size(); p++) + { + Node a_v = pts[0][p]; + Node b_v = pts[1][p]; + + // tangent plane + Node tplane = nm->mkNode( + MINUS, + nm->mkNode( + PLUS, nm->mkNode(MULT, b_v, a), nm->mkNode(MULT, a_v, b)), + nm->mkNode(MULT, a_v, b_v)); + for (unsigned d = 0; d < 4; d++) + { + Node aa = nm->mkNode(d == 0 || d == 3 ? GEQ : LEQ, a, a_v); + Node ab = nm->mkNode(d == 1 || d == 3 ? GEQ : LEQ, b, b_v); + Node conc = nm->mkNode(d <= 1 ? LEQ : GEQ, t, tplane); + Node tlem = nm->mkNode(OR, aa.negate(), ab.negate(), conc); + Trace("nl-ext-tplanes") + << "Tangent plane lemma : " << tlem << std::endl; + lemmas.push_back(tlem); + } + + // tangent plane reverse implication + + // t <= tplane -> ( (a <= a_v ^ b >= b_v) v + // (a >= a_v ^ b <= b_v) ). + // in clause form, the above becomes + // t <= tplane -> a <= a_v v b <= b_v. + // t <= tplane -> b >= b_v v a >= a_v. + Node a_leq_av = nm->mkNode(LEQ, a, a_v); + Node b_leq_bv = nm->mkNode(LEQ, b, b_v); + Node a_geq_av = nm->mkNode(GEQ, a, a_v); + Node b_geq_bv = nm->mkNode(GEQ, b, b_v); + + Node t_leq_tplane = nm->mkNode(LEQ, t, tplane); + Node a_leq_av_or_b_leq_bv = nm->mkNode(OR, a_leq_av, b_leq_bv); + Node b_geq_bv_or_a_geq_av = nm->mkNode(OR, b_geq_bv, a_geq_av); + Node ub_reverse1 = + nm->mkNode(OR, t_leq_tplane.negate(), a_leq_av_or_b_leq_bv); + Trace("nl-ext-tplanes") + << "Tangent plane lemma (reverse) : " << ub_reverse1 + << std::endl; + lemmas.push_back(ub_reverse1); + Node ub_reverse2 = + nm->mkNode(OR, t_leq_tplane.negate(), b_geq_bv_or_a_geq_av); + Trace("nl-ext-tplanes") + << "Tangent plane lemma (reverse) : " << ub_reverse2 + << std::endl; + lemmas.push_back(ub_reverse2); + + // t >= tplane -> ( (a <= a_v ^ b <= b_v) v + // (a >= a_v ^ b >= b_v) ). + // in clause form, the above becomes + // t >= tplane -> a <= a_v v b >= b_v. + // t >= tplane -> b >= b_v v a <= a_v + Node t_geq_tplane = nm->mkNode(GEQ, t, tplane); + Node a_leq_av_or_b_geq_bv = nm->mkNode(OR, a_leq_av, b_geq_bv); + Node a_geq_av_or_b_leq_bv = nm->mkNode(OR, a_geq_av, b_leq_bv); + Node lb_reverse1 = + nm->mkNode(OR, t_geq_tplane.negate(), a_leq_av_or_b_geq_bv); + Trace("nl-ext-tplanes") + << "Tangent plane lemma (reverse) : " << lb_reverse1 + << std::endl; + lemmas.push_back(lb_reverse1); + Node lb_reverse2 = + nm->mkNode(OR, t_geq_tplane.negate(), a_geq_av_or_b_leq_bv); + Trace("nl-ext-tplanes") + << "Tangent plane lemma (reverse) : " << lb_reverse2 + << std::endl; + lemmas.push_back(lb_reverse2); + } + } + } + } + } + Trace("nl-ext") << "...trying " << lemmas.size() << " tangent plane lemmas..." + << std::endl; + return lemmas; +} + +std::vector NlSolver::checkMonomialInferBounds( + std::vector& nt_lemmas, + const std::vector& asserts, + const std::vector& false_asserts) +{ + // sort monomials by degree + Trace("nl-ext-proc") << "Sort monomials by degree..." << std::endl; + d_mdb.sortByDegree(d_ms); + // all monomials + d_mterms.insert(d_mterms.end(), d_ms_vars.begin(), d_ms_vars.end()); + d_mterms.insert(d_mterms.end(), d_ms.begin(), d_ms.end()); + + const std::map >& cim = + d_cdb.getConstraints(); + + std::vector lemmas; + NodeManager* nm = NodeManager::currentNM(); + // register constraints + Trace("nl-ext-debug") << "Register bound constraints..." << std::endl; + for (const Node& lit : asserts) + { + bool polarity = lit.getKind() != NOT; + Node atom = lit.getKind() == NOT ? lit[0] : lit; + d_cdb.registerConstraint(atom); + bool is_false_lit = + std::find(false_asserts.begin(), false_asserts.end(), lit) + != false_asserts.end(); + // add information about bounds to variables + std::map >::const_iterator itc = + cim.find(atom); + if (itc == cim.end()) + { + continue; + } + for (const std::pair& itcc : itc->second) + { + Node x = itcc.first; + Node coeff = itcc.second.d_coeff; + Node rhs = itcc.second.d_rhs; + Kind type = itcc.second.d_type; + Node exp = lit; + if (!polarity) + { + // reverse + if (type == EQUAL) + { + // we will take the strict inequality in the direction of the + // model + Node lhs = ArithMSum::mkCoeffTerm(coeff, x); + Node query = nm->mkNode(GT, lhs, rhs); + Node query_mv = d_model.computeAbstractModelValue(query); + if (query_mv == d_true) + { + exp = query; + type = GT; + } + else + { + Assert(query_mv == d_false); + exp = nm->mkNode(LT, lhs, rhs); + type = LT; + } + } + else + { + type = negateKind(type); + } + } + // add to status if maximal degree + d_ci_max[x][coeff][rhs] = d_cdb.isMaximal(atom, x); + if (Trace.isOn("nl-ext-bound-debug2")) + { + Node t = ArithMSum::mkCoeffTerm(coeff, x); + Trace("nl-ext-bound-debug2") << "Add Bound: " << t << " " << type << " " + << rhs << " by " << exp << std::endl; + } + bool updated = true; + std::map::iterator its = d_ci[x][coeff].find(rhs); + if (its == d_ci[x][coeff].end()) + { + d_ci[x][coeff][rhs] = type; + d_ci_exp[x][coeff][rhs] = exp; + } + else if (type != its->second) + { + Trace("nl-ext-bound-debug2") + << "Joining kinds : " << type << " " << its->second << std::endl; + Kind jk = joinKinds(type, its->second); + if (jk == UNDEFINED_KIND) + { + updated = false; + } + else if (jk != its->second) + { + if (jk == type) + { + d_ci[x][coeff][rhs] = type; + d_ci_exp[x][coeff][rhs] = exp; + } + else + { + d_ci[x][coeff][rhs] = jk; + d_ci_exp[x][coeff][rhs] = + nm->mkNode(AND, d_ci_exp[x][coeff][rhs], exp); + } + } + else + { + updated = false; + } + } + if (Trace.isOn("nl-ext-bound")) + { + if (updated) + { + Trace("nl-ext-bound") << "Bound: "; + debugPrintBound("nl-ext-bound", coeff, x, d_ci[x][coeff][rhs], rhs); + Trace("nl-ext-bound") << " by " << d_ci_exp[x][coeff][rhs]; + if (d_ci_max[x][coeff][rhs]) + { + Trace("nl-ext-bound") << ", is max degree"; + } + Trace("nl-ext-bound") << std::endl; + } + } + // compute if bound is not satisfied, and store what is required + // for a possible refinement + if (options::nlExtTangentPlanes()) + { + if (is_false_lit) + { + d_tplane_refine.insert(x); + } + } + } + } + // reflexive constraints + Node null_coeff; + for (unsigned j = 0; j < d_mterms.size(); j++) + { + Node n = d_mterms[j]; + d_ci[n][null_coeff][n] = EQUAL; + d_ci_exp[n][null_coeff][n] = d_true; + d_ci_max[n][null_coeff][n] = false; + } + + Trace("nl-ext") << "Get inferred bound lemmas..." << std::endl; + const std::map >& cpMap = + d_mdb.getContainsParentMap(); + for (unsigned k = 0; k < d_mterms.size(); k++) + { + Node x = d_mterms[k]; + Trace("nl-ext-bound-debug") + << "Process bounds for " << x << " : " << std::endl; + std::map >::const_iterator itm = cpMap.find(x); + if (itm == cpMap.end()) + { + Trace("nl-ext-bound-debug") << "...has no parent monomials." << std::endl; + continue; + } + Trace("nl-ext-bound-debug") + << "...has " << itm->second.size() << " parent monomials." << std::endl; + // check derived bounds + std::map > >::iterator itc = + d_ci.find(x); + if (itc == d_ci.end()) + { + continue; + } + for (std::map >::iterator itcc = + itc->second.begin(); + itcc != itc->second.end(); + ++itcc) + { + Node coeff = itcc->first; + Node t = ArithMSum::mkCoeffTerm(coeff, x); + for (std::map::iterator itcr = itcc->second.begin(); + itcr != itcc->second.end(); + ++itcr) + { + Node rhs = itcr->first; + // only consider this bound if maximal degree + if (!d_ci_max[x][coeff][rhs]) + { + continue; + } + Kind type = itcr->second; + for (unsigned j = 0; j < itm->second.size(); j++) + { + Node y = itm->second[j]; + Node mult = d_mdb.getContainsDiff(x, y); + // x t => m*x t where y = m*x + // get the sign of mult + Node mmv = d_model.computeConcreteModelValue(mult); + Trace("nl-ext-bound-debug2") + << "Model value of " << mult << " is " << mmv << std::endl; + if (!mmv.isConst()) + { + Trace("nl-ext-bound-debug") + << " ...coefficient " << mult + << " is non-constant (probably transcendental)." << std::endl; + continue; + } + int mmv_sign = mmv.getConst().sgn(); + Trace("nl-ext-bound-debug2") + << " sign of " << mmv << " is " << mmv_sign << std::endl; + if (mmv_sign == 0) + { + Trace("nl-ext-bound-debug") + << " ...coefficient " << mult << " is zero." << std::endl; + continue; + } + Trace("nl-ext-bound-debug") + << " from " << x << " * " << mult << " = " << y << " and " << t + << " " << type << " " << rhs << ", infer : " << std::endl; + Kind infer_type = mmv_sign == -1 ? reverseRelationKind(type) : type; + Node infer_lhs = nm->mkNode(MULT, mult, t); + Node infer_rhs = nm->mkNode(MULT, mult, rhs); + Node infer = nm->mkNode(infer_type, infer_lhs, infer_rhs); + Trace("nl-ext-bound-debug") << " " << infer << std::endl; + infer = Rewriter::rewrite(infer); + Trace("nl-ext-bound-debug2") + << " ...rewritten : " << infer << std::endl; + // check whether it is false in model for abstraction + Node infer_mv = d_model.computeAbstractModelValue(infer); + Trace("nl-ext-bound-debug") + << " ...infer model value is " << infer_mv << std::endl; + if (infer_mv == d_false) + { + Node exp = + nm->mkNode(AND, + nm->mkNode(mmv_sign == 1 ? GT : LT, mult, d_zero), + d_ci_exp[x][coeff][rhs]); + Node iblem = nm->mkNode(IMPLIES, exp, infer); + Node pr_iblem = iblem; + iblem = Rewriter::rewrite(iblem); + bool introNewTerms = hasNewMonomials(iblem, d_ms); + Trace("nl-ext-bound-lemma") + << "*** Bound inference lemma : " << iblem + << " (pre-rewrite : " << pr_iblem << ")" << std::endl; + // Trace("nl-ext-bound-lemma") << " intro new + // monomials = " << introNewTerms << std::endl; + if (!introNewTerms) + { + lemmas.push_back(iblem); + } + else + { + nt_lemmas.push_back(iblem); + } + } + } + } + } + } + return lemmas; +} + +std::vector NlSolver::checkFactoring( + const std::vector& asserts, const std::vector& false_asserts) +{ + std::vector lemmas; + NodeManager* nm = NodeManager::currentNM(); + Trace("nl-ext") << "Get factoring lemmas..." << std::endl; + for (const Node& lit : asserts) + { + bool polarity = lit.getKind() != NOT; + Node atom = lit.getKind() == NOT ? lit[0] : lit; + Node litv = d_model.computeConcreteModelValue(lit); + bool considerLit = false; + // Only consider literals that are in false_asserts. + considerLit = std::find(false_asserts.begin(), false_asserts.end(), lit) + != false_asserts.end(); + + if (considerLit) + { + std::map msum; + if (ArithMSum::getMonomialSumLit(atom, msum)) + { + Trace("nl-ext-factor") << "Factoring for literal " << lit + << ", monomial sum is : " << std::endl; + if (Trace.isOn("nl-ext-factor")) + { + ArithMSum::debugPrintMonomialSum(msum, "nl-ext-factor"); + } + std::map > factor_to_mono; + std::map > factor_to_mono_orig; + for (std::map::iterator itm = msum.begin(); + itm != msum.end(); + ++itm) + { + if (!itm->first.isNull()) + { + if (itm->first.getKind() == NONLINEAR_MULT) + { + std::vector children; + for (unsigned i = 0; i < itm->first.getNumChildren(); i++) + { + children.push_back(itm->first[i]); + } + std::map processed; + for (unsigned i = 0; i < itm->first.getNumChildren(); i++) + { + if (processed.find(itm->first[i]) == processed.end()) + { + processed[itm->first[i]] = true; + children[i] = d_one; + if (!itm->second.isNull()) + { + children.push_back(itm->second); + } + Node val = nm->mkNode(MULT, children); + if (!itm->second.isNull()) + { + children.pop_back(); + } + children[i] = itm->first[i]; + val = Rewriter::rewrite(val); + factor_to_mono[itm->first[i]].push_back(val); + factor_to_mono_orig[itm->first[i]].push_back(itm->first); + } + } + } + } + } + for (std::map >::iterator itf = + factor_to_mono.begin(); + itf != factor_to_mono.end(); + ++itf) + { + Node x = itf->first; + if (itf->second.size() == 1) + { + std::map::iterator itm = msum.find(x); + if (itm != msum.end()) + { + itf->second.push_back(itm->second.isNull() ? d_one : itm->second); + factor_to_mono_orig[x].push_back(x); + } + } + if (itf->second.size() <= 1) + { + continue; + } + Node sum = nm->mkNode(PLUS, itf->second); + sum = Rewriter::rewrite(sum); + Trace("nl-ext-factor") + << "* Factored sum for " << x << " : " << sum << std::endl; + Node kf = getFactorSkolem(sum, lemmas); + std::vector poly; + poly.push_back(nm->mkNode(MULT, x, kf)); + std::map >::iterator itfo = + factor_to_mono_orig.find(x); + Assert(itfo != factor_to_mono_orig.end()); + for (std::map::iterator itm = msum.begin(); + itm != msum.end(); + ++itm) + { + if (std::find(itfo->second.begin(), itfo->second.end(), itm->first) + == itfo->second.end()) + { + poly.push_back(ArithMSum::mkCoeffTerm( + itm->second, itm->first.isNull() ? d_one : itm->first)); + } + } + Node polyn = poly.size() == 1 ? poly[0] : nm->mkNode(PLUS, poly); + Trace("nl-ext-factor") + << "...factored polynomial : " << polyn << std::endl; + Node conc_lit = nm->mkNode(atom.getKind(), polyn, d_zero); + conc_lit = Rewriter::rewrite(conc_lit); + if (!polarity) + { + conc_lit = conc_lit.negate(); + } + + std::vector lemma_disj; + lemma_disj.push_back(lit.negate()); + lemma_disj.push_back(conc_lit); + Node flem = nm->mkNode(OR, lemma_disj); + Trace("nl-ext-factor") << "...lemma is " << flem << std::endl; + lemmas.push_back(flem); + } + } + } + } + return lemmas; +} + +Node NlSolver::getFactorSkolem(Node n, std::vector& lemmas) +{ + std::map::iterator itf = d_factor_skolem.find(n); + if (itf == d_factor_skolem.end()) + { + NodeManager* nm = NodeManager::currentNM(); + Node k = nm->mkSkolem("kf", n.getType()); + Node k_eq = Rewriter::rewrite(k.eqNode(n)); + lemmas.push_back(k_eq); + d_factor_skolem[n] = k; + return k; + } + return itf->second; +} + +std::vector NlSolver::checkMonomialInferResBounds() +{ + std::vector lemmas; + NodeManager* nm = NodeManager::currentNM(); + Trace("nl-ext") << "Get monomial resolution inferred bound lemmas..." + << std::endl; + size_t nmterms = d_mterms.size(); + for (unsigned j = 0; j < nmterms; j++) + { + Node a = d_mterms[j]; + std::map > >::iterator itca = + d_ci.find(a); + if (itca == d_ci.end()) + { + continue; + } + for (unsigned k = (j + 1); k < nmterms; k++) + { + Node b = d_mterms[k]; + std::map > >::iterator itcb = + d_ci.find(b); + if (itcb == d_ci.end()) + { + continue; + } + Trace("nl-ext-rbound-debug") << "resolution inferences : compare " << a + << " and " << b << std::endl; + // if they have common factors + std::map::iterator ita = d_mono_diff[a].find(b); + if (ita == d_mono_diff[a].end()) + { + continue; + } + Trace("nl-ext-rbound") << "Get resolution inferences for [a] " << a + << " vs [b] " << b << std::endl; + std::map::iterator itb = d_mono_diff[b].find(a); + Assert(itb != d_mono_diff[b].end()); + Node mv_a = d_model.computeAbstractModelValue(ita->second); + Assert(mv_a.isConst()); + int mv_a_sgn = mv_a.getConst().sgn(); + if (mv_a_sgn == 0) + { + // we don't compare monomials whose current model value is zero + continue; + } + Node mv_b = d_model.computeAbstractModelValue(itb->second); + Assert(mv_b.isConst()); + int mv_b_sgn = mv_b.getConst().sgn(); + if (mv_b_sgn == 0) + { + // we don't compare monomials whose current model value is zero + continue; + } + Trace("nl-ext-rbound") << " [a] factor is " << ita->second + << ", sign in model = " << mv_a_sgn << std::endl; + Trace("nl-ext-rbound") << " [b] factor is " << itb->second + << ", sign in model = " << mv_b_sgn << std::endl; + + std::vector exp; + // bounds of a + for (std::map >::iterator itcac = + itca->second.begin(); + itcac != itca->second.end(); + ++itcac) + { + Node coeff_a = itcac->first; + for (std::map::iterator itcar = itcac->second.begin(); + itcar != itcac->second.end(); + ++itcar) + { + Node rhs_a = itcar->first; + Node rhs_a_res_base = nm->mkNode(MULT, itb->second, rhs_a); + rhs_a_res_base = Rewriter::rewrite(rhs_a_res_base); + if (hasNewMonomials(rhs_a_res_base, d_ms)) + { + continue; + } + Kind type_a = itcar->second; + exp.push_back(d_ci_exp[a][coeff_a][rhs_a]); + + // bounds of b + for (std::map >::iterator itcbc = + itcb->second.begin(); + itcbc != itcb->second.end(); + ++itcbc) + { + Node coeff_b = itcbc->first; + Node rhs_a_res = ArithMSum::mkCoeffTerm(coeff_b, rhs_a_res_base); + for (std::map::iterator itcbr = itcbc->second.begin(); + itcbr != itcbc->second.end(); + ++itcbr) + { + Node rhs_b = itcbr->first; + Node rhs_b_res = nm->mkNode(MULT, ita->second, rhs_b); + rhs_b_res = ArithMSum::mkCoeffTerm(coeff_a, rhs_b_res); + rhs_b_res = Rewriter::rewrite(rhs_b_res); + if (hasNewMonomials(rhs_b_res, d_ms)) + { + continue; + } + Kind type_b = itcbr->second; + exp.push_back(d_ci_exp[b][coeff_b][rhs_b]); + if (Trace.isOn("nl-ext-rbound")) + { + Trace("nl-ext-rbound") << "* try bounds : "; + debugPrintBound("nl-ext-rbound", coeff_a, a, type_a, rhs_a); + Trace("nl-ext-rbound") << std::endl; + Trace("nl-ext-rbound") << " "; + debugPrintBound("nl-ext-rbound", coeff_b, b, type_b, rhs_b); + Trace("nl-ext-rbound") << std::endl; + } + Kind types[2]; + for (unsigned r = 0; r < 2; r++) + { + Node pivot_factor = r == 0 ? itb->second : ita->second; + int pivot_factor_sign = r == 0 ? mv_b_sgn : mv_a_sgn; + types[r] = r == 0 ? type_a : type_b; + if (pivot_factor_sign == (r == 0 ? 1 : -1)) + { + types[r] = reverseRelationKind(types[r]); + } + if (pivot_factor_sign == 1) + { + exp.push_back(nm->mkNode(GT, pivot_factor, d_zero)); + } + else + { + exp.push_back(nm->mkNode(LT, pivot_factor, d_zero)); + } + } + Kind jk = transKinds(types[0], types[1]); + Trace("nl-ext-rbound-debug") + << "trans kind : " << types[0] << " + " << types[1] << " = " + << jk << std::endl; + if (jk != UNDEFINED_KIND) + { + Node conc = nm->mkNode(jk, rhs_a_res, rhs_b_res); + Node conc_mv = d_model.computeAbstractModelValue(conc); + if (conc_mv == d_false) + { + Node rblem = nm->mkNode(IMPLIES, nm->mkNode(AND, exp), conc); + Trace("nl-ext-rbound-lemma-debug") + << "Resolution bound lemma " + "(pre-rewrite) " + ": " + << rblem << std::endl; + rblem = Rewriter::rewrite(rblem); + Trace("nl-ext-rbound-lemma") + << "Resolution bound lemma : " << rblem << std::endl; + lemmas.push_back(rblem); + } + } + exp.pop_back(); + exp.pop_back(); + exp.pop_back(); + } + } + exp.pop_back(); + } + } + } + } + return lemmas; +} + +} // namespace arith +} // namespace theory +} // namespace CVC4 diff --git a/src/theory/arith/nl_solver.h b/src/theory/arith/nl_solver.h new file mode 100644 index 000000000..73ca97f00 --- /dev/null +++ b/src/theory/arith/nl_solver.h @@ -0,0 +1,368 @@ +/********************* */ +/*! \file nl_solver.h + ** \verbatim + ** Top contributors (to current version): + ** Andrew Reynolds, Tim King + ** This file is part of the CVC4 project. + ** Copyright (c) 2009-2019 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 Solver for standard non-linear constraints + **/ + +#ifndef CVC4__THEORY__ARITH__NL_SOLVER_H +#define CVC4__THEORY__ARITH__NL_SOLVER_H + +#include +#include +#include +#include + +#include "context/cdhashset.h" +#include "context/cdinsert_hashmap.h" +#include "context/cdlist.h" +#include "context/cdqueue.h" +#include "context/context.h" +#include "expr/kind.h" +#include "expr/node.h" +#include "theory/arith/nl_constraint.h" +#include "theory/arith/nl_lemma_utils.h" +#include "theory/arith/nl_model.h" +#include "theory/arith/nl_monomial.h" +#include "theory/arith/theory_arith.h" + +namespace CVC4 { +namespace theory { +namespace arith { + +typedef std::map NodeMultiset; + +/** Non-linear solver class + * + * This class implements model-based refinement schemes + * for non-linear arithmetic, described in: + * + * - "Invariant Checking of NRA Transition Systems + * via Incremental Reduction to LRA with EUF" by + * Cimatti et al., TACAS 2017. + * + * - Section 5 of "Desiging Theory Solvers with + * Extensions" by Reynolds et al., FroCoS 2017. + */ +class NlSolver +{ + typedef std::map MonomialExponentMap; + typedef context::CDHashSet NodeSet; + + public: + NlSolver(TheoryArith& containing, NlModel& model); + ~NlSolver(); + + /** init last call + * + * This is called at the beginning of last call effort check, where + * assertions are the set of assertions belonging to arithmetic, + * false_asserts is the subset of assertions that are false in the current + * model, and xts is the set of extended function terms that are active in + * the current context. + */ + void initLastCall(const std::vector& assertions, + const std::vector& false_asserts, + const std::vector& xts); + //-------------------------------------------- lemma schemas + /** check split zero + * + * Returns a set of theory lemmas of the form + * t = 0 V t != 0 + * where t is a term that exists in the current context. + */ + std::vector checkSplitZero(); + + /** check monomial sign + * + * Returns a set of valid theory lemmas, based on a + * lemma schema which ensures that non-linear monomials + * respect sign information based on their facts. + * For more details, see Section 5 of "Design Theory + * Solvers with Extensions" by Reynolds et al., FroCoS 2017, + * Figure 5, this is the schema "Sign". + * + * Examples: + * + * x > 0 ^ y > 0 => x*y > 0 + * x < 0 => x*y*y < 0 + * x = 0 => x*y*z = 0 + */ + std::vector checkMonomialSign(); + + /** check monomial magnitude + * + * Returns a set of valid theory lemmas, based on a + * lemma schema which ensures that comparisons between + * non-linear monomials respect the magnitude of their + * factors. + * For more details, see Section 5 of "Design Theory + * Solvers with Extensions" by Reynolds et al., FroCoS 2017, + * Figure 5, this is the schema "Magnitude". + * + * Examples: + * + * |x|>|y| => |x*z|>|y*z| + * |x|>|y| ^ |z|>|w| ^ |x|>=1 => |x*x*z*u|>|y*w| + * + * Argument c indicates the class of inferences to perform for the + * (non-linear) monomials in the vector d_ms. 0 : compare non-linear monomials + * against 1, 1 : compare non-linear monomials against variables, 2 : compare + * non-linear monomials against other non-linear monomials. + */ + std::vector checkMonomialMagnitude(unsigned c); + + /** check monomial inferred bounds + * + * Returns a set of valid theory lemmas, based on a + * lemma schema that infers new constraints about existing + * terms based on mulitplying both sides of an existing + * constraint by a term. + * For more details, see Section 5 of "Design Theory + * Solvers with Extensions" by Reynolds et al., FroCoS 2017, + * Figure 5, this is the schema "Multiply". + * + * Examples: + * + * x > 0 ^ (y > z + w) => x*y > x*(z+w) + * x < 0 ^ (y > z + w) => x*y < x*(z+w) + * ...where (y > z + w) and x*y are a constraint and term + * that occur in the current context. + */ + std::vector checkMonomialInferBounds( + std::vector& nt_lemmas, + const std::vector& asserts, + const std::vector& false_asserts); + + /** check factoring + * + * Returns a set of valid theory lemmas, based on a + * lemma schema that states a relationship betwen monomials + * with common factors that occur in the same constraint. + * + * Examples: + * + * x*z+y*z > t => ( k = x + y ^ k*z > t ) + * ...where k is fresh and x*z + y*z > t is a + * constraint that occurs in the current context. + */ + std::vector checkFactoring(const std::vector& asserts, + const std::vector& false_asserts); + + /** check monomial infer resolution bounds + * + * Returns a set of valid theory lemmas, based on a + * lemma schema which "resolves" upper bounds + * of one inequality with lower bounds for another. + * This schema is not enabled by default, and can + * be enabled by --nl-ext-rbound. + * + * Examples: + * + * ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t + * ...where s <= x*z and x*y <= t are constraints + * that occur in the current context. + */ + std::vector checkMonomialInferResBounds(); + + /** check tangent planes + * + * Returns a set of valid theory lemmas, based on an + * "incremental linearization" of non-linear monomials. + * This linearization is accomplished by adding constraints + * corresponding to "tangent planes" at the current + * model value of each non-linear monomial. In particular + * consider the definition for constants a,b : + * T_{a,b}( x*y ) = b*x + a*y - a*b. + * The lemmas added by this function are of the form : + * ( ( x>a ^ yb) ) => x*y < T_{a,b}( x*y ) + * ( ( x>a ^ y>b) ^ (x x*y > T_{a,b}( x*y ) + * It is inspired by "Invariant Checking of NRA Transition + * Systems via Incremental Reduction to LRA with EUF" by + * Cimatti et al., TACAS 2017. + * This schema is not terminating in general. + * It is not enabled by default, and can + * be enabled by --nl-ext-tplanes. + * + * Examples: + * + * ( ( x>2 ^ y>5) ^ (x<2 ^ y<5) ) => x*y > 5*x + 2*y - 10 + * ( ( x>2 ^ y<5) ^ (x<2 ^ y>5) ) => x*y < 5*x + 2*y - 10 + */ + std::vector checkTangentPlanes(); + + //-------------------------------------------- end lemma schemas + private: + // The theory of arithmetic containing this extension. + TheoryArith& d_containing; + /** Reference to the non-linear model object */ + NlModel& d_model; + /** commonly used terms */ + Node d_zero; + Node d_one; + Node d_neg_one; + Node d_two; + Node d_true; + Node d_false; + /** Context-independent database of monomial information */ + MonomialDb d_mdb; + /** Context-independent database of constraint information */ + ConstraintDb d_cdb; + + // ( x*y, x*z, y ) for each pair of monomials ( x*y, x*z ) with common factors + std::map > d_mono_diff; + + /** cache of terms t for which we have added the lemma ( t = 0 V t != 0 ). */ + NodeSet d_zero_split; + + // ordering, stores variables and 0,1,-1 + std::map d_order_vars; + std::vector d_order_points; + + // information about monomials + std::vector d_ms; + std::vector d_ms_vars; + std::map d_ms_proc; + std::vector d_mterms; + + // list of monomials with factors whose model value is non-constant in model + // e.g. y*cos( x ) + std::map d_m_nconst_factor; + /** the set of monomials we should apply tangent planes to */ + std::unordered_set d_tplane_refine; + /** maps nodes to their factor skolems */ + std::map d_factor_skolem; + /** tangent plane bounds */ + std::map > d_tangent_val_bound[4]; + // term -> coeff -> rhs -> ( status, exp, b ), + // where we have that : exp => ( coeff * term rhs ) + // b is true if degree( term ) >= degree( rhs ) + std::map > > d_ci; + std::map > > d_ci_exp; + std::map > > d_ci_max; + + /** Make literal */ + static Node mkLit(Node a, Node b, int status, bool isAbsolute = false); + /** register monomial */ + void setMonomialFactor(Node a, Node b, const NodeMultiset& common); + /** assign order ids */ + void assignOrderIds(std::vector& vars, + NodeMultiset& d_order, + bool isConcrete, + bool isAbsolute); + + /** Check whether we have already inferred a relationship between monomials + * x and y based on the information in cmp_infers. This computes the + * transitive closure of the relation stored in cmp_infers. + */ + bool cmp_holds(Node x, + Node y, + std::map >& cmp_infers, + std::vector& exp, + std::map& visited); + /** In the following functions, status states a relationship + * between two arithmetic terms, where: + * 0 : equal + * 1 : greater than or equal + * 2 : greater than + * -X : (greater -> less) + * TODO (#1287) make this an enum? + */ + /** compute the sign of a. + * + * Calls to this function are such that : + * exp => ( oa = a ^ a 0 ) + * + * This function iterates over the factors of a, + * where a_index is the index of the factor in a + * we are currently looking at. + * + * This function returns a status, which indicates + * a's relationship to 0. + * We add lemmas to lem of the form given by the + * lemma schema checkSign(...). + */ + int compareSign(Node oa, + Node a, + unsigned a_index, + int status, + std::vector& exp, + std::vector& lem); + /** compare monomials a and b + * + * Initially, a call to this function is such that : + * exp => ( oa = a ^ ob = b ) + * + * This function returns true if we can infer a valid + * arithmetic lemma of the form : + * P => abs( a ) >= abs( b ) + * where P is true and abs( a ) >= abs( b ) is false in the + * current model. + * + * This function is implemented by "processing" factors + * of monomials a and b until an inference of the above + * form can be made. For example, if : + * a = x*x*y and b = z*w + * Assuming we are trying to show abs( a ) >= abs( c ), + * then if abs( M( x ) ) >= abs( M( z ) ) where M is the current model, + * then we can add abs( x ) >= abs( z ) to our explanation, and + * mark one factor of x as processed in a, and + * one factor of z as processed in b. The number of processed factors of a + * and b are stored in a_exp_proc and b_exp_proc respectively. + * + * cmp_infers stores information that is helpful + * in discarding redundant inferences. For example, + * we do not want to infer abs( x ) >= abs( z ) if + * we have already inferred abs( x ) >= abs( y ) and + * abs( y ) >= abs( z ). + * It stores entries of the form (status,t1,t2)->F, + * which indicates that we constructed a lemma F that + * showed t1 t2. + * + * We add lemmas to lem of the form given by the + * lemma schema checkMagnitude(...). + */ + bool compareMonomial( + Node oa, + Node a, + NodeMultiset& a_exp_proc, + Node ob, + Node b, + NodeMultiset& b_exp_proc, + std::vector& exp, + std::vector& lem, + std::map > >& cmp_infers); + /** helper function for above + * + * The difference is the inputs a_index and b_index, which are the indices of + * children (factors) in monomials a and b which we are currently looking at. + */ + bool compareMonomial( + Node oa, + Node a, + unsigned a_index, + NodeMultiset& a_exp_proc, + Node ob, + Node b, + unsigned b_index, + NodeMultiset& b_exp_proc, + int status, + std::vector& exp, + std::vector& lem, + std::map > >& cmp_infers); + /** Get factor skolem for n, add resulting lemmas to lemmas */ + Node getFactorSkolem(Node n, std::vector& lemmas); +}; /* class NlSolver */ + +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif /* CVC4__THEORY__ARITH__NL_SOLVER_H */ diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp index 61a8b18b0..b638d8a59 100644 --- a/src/theory/arith/nonlinear_extension.cpp +++ b/src/theory/arith/nonlinear_extension.cpp @@ -17,17 +17,10 @@ #include "theory/arith/nonlinear_extension.h" -#include -#include - -#include "expr/node_algorithm.h" -#include "expr/node_builder.h" #include "options/arith_options.h" -#include "theory/arith/arith_msum.h" #include "theory/arith/arith_utilities.h" #include "theory/arith/theory_arith.h" #include "theory/ext_theory.h" -#include "theory/quantifiers/quant_util.h" #include "theory/theory_model.h" using namespace CVC4::kind; @@ -36,181 +29,25 @@ namespace CVC4 { namespace theory { namespace arith { -namespace { - -// Return true if a collection c contains an elem k. Compatible with map and set -// containers. -template -bool Contains(const Container& c, const Key& k) { - return c.find(k) != c.end(); -} - -// Inserts value into the set/map c if the value was not present there. Returns -// true if the value was inserted. -template -bool InsertIfNotPresent(Container* c, const Value& value) { - return (c->insert(value)).second; -} - -// Returns true if a vector c contains an elem t. -template -bool IsInVector(const std::vector& c, const T& t) { - return std::find(c.begin(), c.end(), t) != c.end(); -} - -// Returns the a[key] and assertion fails in debug mode. -inline unsigned getCount(const NodeMultiset& a, Node key) { - NodeMultiset::const_iterator it = a.find(key); - Assert(it != a.end()); - return it->second; -} - -// Returns a[key] if key is in a or value otherwise. -unsigned getCountWithDefault(const NodeMultiset& a, Node key, unsigned value) { - NodeMultiset::const_iterator it = a.find(key); - return (it == a.end()) ? value : it->second; -} - -// Returns true if for any key then a[key] <= b[key] where the value for any key -// not present is interpreted as 0. -bool isSubset(const NodeMultiset& a, const NodeMultiset& b) { - for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) { - Node key = it_a->first; - const unsigned a_value = it_a->second; - const unsigned b_value = getCountWithDefault(b, key, 0); - if (a_value > b_value) { - return false; - } - } - return true; -} - -// Given two multisets return the multiset difference a \ b. -NodeMultiset diffMultiset(const NodeMultiset& a, const NodeMultiset& b) { - NodeMultiset difference; - for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) { - Node key = it_a->first; - const unsigned a_value = it_a->second; - const unsigned b_value = getCountWithDefault(b, key, 0); - if (a_value > b_value) { - difference[key] = a_value - b_value; - } - } - return difference; -} - -// Return a vector containing a[key] repetitions of key in a multiset a. -std::vector ExpandMultiset(const NodeMultiset& a) { - std::vector expansion; - for (NodeMultiset::const_iterator it_a = a.begin(); it_a != a.end(); ++it_a) { - expansion.insert(expansion.end(), it_a->second, it_a->first); - } - return expansion; -} - -void debugPrintBound(const char* c, Node coeff, Node x, Kind type, Node rhs) { - Node t = ArithMSum::mkCoeffTerm(coeff, x); - Trace(c) << t << " " << type << " " << rhs; -} - -bool hasNewMonomials(Node n, const std::vector& existing) { - std::set visited; - - std::vector worklist; - worklist.push_back(n); - while (!worklist.empty()) { - Node current = worklist.back(); - worklist.pop_back(); - if (!Contains(visited, current)) { - visited.insert(current); - if (current.getKind() == NONLINEAR_MULT) - { - if (!IsInVector(existing, current)) { - return true; - } - } else { - worklist.insert(worklist.end(), current.begin(), current.end()); - } - } - } - return false; -} - -} // namespace - NonlinearExtension::NonlinearExtension(TheoryArith& containing, eq::EqualityEngine* ee) : d_lemmas(containing.getUserContext()), - d_zero_split(containing.getUserContext()), d_containing(containing), d_ee(ee), d_needsLastCall(false), d_model(containing.getSatContext()), d_trSlv(d_model), + d_nlSlv(containing, d_model), d_builtModel(containing.getSatContext(), false) { d_true = NodeManager::currentNM()->mkConst(true); - d_false = NodeManager::currentNM()->mkConst(false); d_zero = NodeManager::currentNM()->mkConst(Rational(0)); d_one = NodeManager::currentNM()->mkConst(Rational(1)); d_neg_one = NodeManager::currentNM()->mkConst(Rational(-1)); - d_two = NodeManager::currentNM()->mkConst(Rational(2)); - d_order_points.push_back(d_neg_one); - d_order_points.push_back(d_zero); - d_order_points.push_back(d_one); } NonlinearExtension::~NonlinearExtension() {} -// Returns a reference to either map[key] if it exists in the map -// or to a default value otherwise. -// -// Warning: sped_cial care must be taken if value is a temporary object. -template -const Value& FindWithDefault(const MapType& map, const Key& key, - const Value& value) { - typename MapType::const_iterator it = map.find(key); - if (it == map.end()) { - return value; - } - return it->second; -} - -const NodeMultiset& NonlinearExtension::getMonomialExponentMap( - Node monomial) const { - MonomialExponentMap::const_iterator it = d_m_exp.find(monomial); - Assert(it != d_m_exp.end()); - return it->second; -} - -bool NonlinearExtension::isMonomialSubset(Node a, Node b) const { - const NodeMultiset& a_exponent_map = getMonomialExponentMap(a); - const NodeMultiset& b_exponent_map = getMonomialExponentMap(b); - - return isSubset(a_exponent_map, b_exponent_map); -} - -void NonlinearExtension::registerMonomialSubset(Node a, Node b) { - Assert(isMonomialSubset(a, b)); - - const NodeMultiset& a_exponent_map = getMonomialExponentMap(a); - const NodeMultiset& b_exponent_map = getMonomialExponentMap(b); - - std::vector diff_children = - ExpandMultiset(diffMultiset(b_exponent_map, a_exponent_map)); - Assert(!diff_children.empty()); - - d_m_contain_parent[a].push_back(b); - d_m_contain_children[b].push_back(a); - - Node mult_term = safeConstructNary(MULT, diff_children); - Node nlmult_term = safeConstructNary(NONLINEAR_MULT, diff_children); - d_m_contain_mult[a][b] = mult_term; - d_m_contain_umult[a][b] = nlmult_term; - Trace("nl-ext-mindex") << "..." << a << " is a subset of " << b - << ", difference is " << mult_term << std::endl; -} - bool NonlinearExtension::getCurrentSubstitution( int effort, const std::vector& vars, std::vector& subs, std::map >& exp) { @@ -239,8 +76,6 @@ bool NonlinearExtension::getCurrentSubstitution( // return true if the substitution is non-trivial return retVal; - - // d_containing.getValuation().getModel()->getRepresentative( n ); } std::pair NonlinearExtension::isExtfReduced( @@ -295,194 +130,6 @@ std::pair NonlinearExtension::isExtfReduced( return std::make_pair(true, Node::null()); } -void NonlinearExtension::registerMonomial(Node n) { - if (!IsInVector(d_monomials, n)) { - d_monomials.push_back(n); - Trace("nl-ext-debug") << "Register monomial : " << n << std::endl; - if (n.getKind() == NONLINEAR_MULT) - { - // get exponent count - for (unsigned k = 0; k < n.getNumChildren(); k++) { - d_m_exp[n][n[k]]++; - if (k == 0 || n[k] != n[k - 1]) { - d_m_vlist[n].push_back(n[k]); - } - } - d_m_degree[n] = n.getNumChildren(); - } else if (n == d_one) { - d_m_exp[n].clear(); - d_m_vlist[n].clear(); - d_m_degree[n] = 0; - } else { - Assert(!isArithKind(n.getKind())); - d_m_exp[n][n] = 1; - d_m_vlist[n].push_back(n); - d_m_degree[n] = 1; - } - // TODO: sort necessary here? - std::sort(d_m_vlist[n].begin(), d_m_vlist[n].end()); - Trace("nl-ext-mindex") << "Add monomial to index : " << n << std::endl; - d_m_index.addTerm(n, d_m_vlist[n], this); - } -} - -void NonlinearExtension::setMonomialFactor(Node a, Node b, - const NodeMultiset& common) { - // Could not tell if this was being inserted intentionally or not. - std::map& mono_diff_a = d_mono_diff[a]; - if (!Contains(mono_diff_a, b)) { - Trace("nl-ext-mono-factor") - << "Set monomial factor for " << a << "/" << b << std::endl; - mono_diff_a[b] = mkMonomialRemFactor(a, common); - } -} - -void NonlinearExtension::registerConstraint(Node atom) { - if (!IsInVector(d_constraints, atom)) { - d_constraints.push_back(atom); - Trace("nl-ext-debug") << "Register constraint : " << atom << std::endl; - std::map msum; - if (ArithMSum::getMonomialSumLit(atom, msum)) - { - Trace("nl-ext-debug") << "got monomial sum: " << std::endl; - if (Trace.isOn("nl-ext-debug")) { - ArithMSum::debugPrintMonomialSum(msum, "nl-ext-debug"); - } - unsigned max_degree = 0; - std::vector all_m; - std::vector max_deg_m; - for (std::map::iterator itm = msum.begin(); itm != msum.end(); - ++itm) { - if (!itm->first.isNull()) { - all_m.push_back(itm->first); - registerMonomial(itm->first); - Trace("nl-ext-debug2") - << "...process monomial " << itm->first << std::endl; - Assert(d_m_degree.find(itm->first) != d_m_degree.end()); - unsigned d = d_m_degree[itm->first]; - if (d > max_degree) { - max_degree = d; - max_deg_m.clear(); - } - if (d >= max_degree) { - max_deg_m.push_back(itm->first); - } - } - } - // isolate for each maximal degree monomial - for (unsigned i = 0; i < all_m.size(); i++) { - Node m = all_m[i]; - Node rhs, coeff; - int res = ArithMSum::isolate(m, msum, coeff, rhs, atom.getKind()); - if (res != 0) { - Kind type = atom.getKind(); - if (res == -1) { - type = reverseRelationKind(type); - } - Trace("nl-ext-constraint") << "Constraint : " << atom << " <=> "; - if (!coeff.isNull()) { - Trace("nl-ext-constraint") << coeff << " * "; - } - Trace("nl-ext-constraint") - << m << " " << type << " " << rhs << std::endl; - d_c_info[atom][m].d_rhs = rhs; - d_c_info[atom][m].d_coeff = coeff; - d_c_info[atom][m].d_type = type; - } - } - for (unsigned i = 0; i < max_deg_m.size(); i++) { - Node m = max_deg_m[i]; - d_c_info_maxm[atom][m] = true; - } - } else { - Trace("nl-ext-debug") << "...failed to get monomial sum." << std::endl; - } - } -} - -bool NonlinearExtension::isArithKind(Kind k) { - return k == PLUS || k == MULT || k == NONLINEAR_MULT; -} - -Node NonlinearExtension::mkLit(Node a, Node b, int status, bool isAbsolute) -{ - if (status == 0) { - Node a_eq_b = a.eqNode(b); - if (!isAbsolute) - { - return a_eq_b; - } - else - { - // return mkAbs( a ).eqNode( mkAbs( b ) ); - Node negate_b = NodeManager::currentNM()->mkNode(UMINUS, b); - return a_eq_b.orNode(a.eqNode(negate_b)); - } - } else if (status < 0) { - return mkLit(b, a, -status); - } else { - Assert(status == 1 || status == 2); - NodeManager* nm = NodeManager::currentNM(); - Kind greater_op = status == 1 ? GEQ : GT; - if (!isAbsolute) - { - return nm->mkNode(greater_op, a, b); - } - else - { - // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) ); - Node zero = mkRationalNode(0); - Node a_is_nonnegative = nm->mkNode(GEQ, a, zero); - Node b_is_nonnegative = nm->mkNode(GEQ, b, zero); - Node negate_a = nm->mkNode(UMINUS, a); - Node negate_b = nm->mkNode(UMINUS, b); - return a_is_nonnegative.iteNode( - b_is_nonnegative.iteNode(nm->mkNode(greater_op, a, b), - nm->mkNode(greater_op, a, negate_b)), - b_is_nonnegative.iteNode(nm->mkNode(greater_op, negate_a, b), - nm->mkNode(greater_op, negate_a, negate_b))); - } - } -} - -Node NonlinearExtension::mkAbs(Node a) { - if (a.isConst()) { - return mkRationalNode(a.getConst().abs()); - } else { - NodeManager* nm = NodeManager::currentNM(); - Node a_is_nonnegative = nm->mkNode(GEQ, a, mkRationalNode(0)); - return a_is_nonnegative.iteNode(a, nm->mkNode(UMINUS, a)); - } -} - -Node NonlinearExtension::mkValidPhase(Node a, Node pi) { - return mkBounded( - NodeManager::currentNM()->mkNode(MULT, mkRationalNode(-1), pi), a, pi); -} - -Node NonlinearExtension::mkMonomialRemFactor( - Node n, const NodeMultiset& n_exp_rem) const { - std::vector children; - const NodeMultiset& exponent_map = getMonomialExponentMap(n); - for (NodeMultiset::const_iterator itme2 = exponent_map.begin(); - itme2 != exponent_map.end(); ++itme2) { - Node v = itme2->first; - unsigned inc = itme2->second; - Trace("nl-ext-mono-factor") - << "..." << inc << " factors of " << v << std::endl; - unsigned count_in_n_exp_rem = getCountWithDefault(n_exp_rem, v, 0); - Assert(count_in_n_exp_rem <= inc); - inc -= count_in_n_exp_rem; - Trace("nl-ext-mono-factor") - << "......rem, now " << inc << " factors of " << v << std::endl; - children.insert(children.end(), inc, v); - } - Node ret = safeConstructNary(MULT, children); - ret = Rewriter::rewrite(ret); - Trace("nl-ext-mono-factor") << "...return : " << ret << std::endl; - return ret; -} - void NonlinearExtension::sendLemmas(const std::vector& out, bool preprocess, std::map& lemSE) @@ -732,22 +379,6 @@ bool NonlinearExtension::checkModel(const std::vector& assertions, return ret; } - -std::vector NonlinearExtension::checkSplitZero() { - std::vector lemmas; - for (unsigned i = 0; i < d_ms_vars.size(); i++) { - Node v = d_ms_vars[i]; - if (d_zero_split.insert(v)) { - Node eq = v.eqNode(d_zero); - eq = Rewriter::rewrite(eq); - Node literal = d_containing.getValuation().ensureLiteral(eq); - d_containing.getOutputChannel().requirePhase(literal, true); - lemmas.push_back(literal.orNode(literal.negate())); - } - } - return lemmas; -} - int NonlinearExtension::checkLastCall(const std::vector& assertions, const std::vector& false_asserts, const std::vector& xts, @@ -756,48 +387,8 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, std::vector& wlems, std::map& lemSE) { - d_ms_vars.clear(); - d_ms_proc.clear(); - d_ms.clear(); - d_mterms.clear(); - d_m_nconst_factor.clear(); - d_tplane_refine.clear(); - d_ci.clear(); - d_ci_exp.clear(); - d_ci_max.clear(); - - Trace("nl-ext-mv") << "Extended terms : " << std::endl; - // for computing congruence - std::map argTrie; - for (unsigned i = 0, xsize = xts.size(); i < xsize; i++) - { - Node a = xts[i]; - d_model.computeConcreteModelValue(a); - d_model.computeAbstractModelValue(a); - d_model.printModelValue("nl-ext-mv", a); - Kind ak = a.getKind(); - if (ak == NONLINEAR_MULT) - { - d_ms.push_back( a ); - - //context-independent registration - registerMonomial(a); - - std::map >::iterator itvl = d_m_vlist.find(a); - Assert(itvl != d_m_vlist.end()); - for (unsigned k = 0; k < itvl->second.size(); k++) { - if (!IsInVector(d_ms_vars, itvl->second[k])) { - d_ms_vars.push_back(itvl->second[k]); - } - Node mvk = d_model.computeAbstractModelValue(itvl->second[k]); - if( !mvk.isConst() ){ - d_m_nconst_factor[a] = true; - } - } - // mark processed if has a "one" factor (will look at reduced monomial)? - } - } - + // initialize the non-linear solver + d_nlSlv.initLastCall(assertions, false_asserts, xts); // initialize the trancendental function solver std::vector lemmas; d_trSlv.initLastCall(assertions, false_asserts, xts, lemmas, lemsPp); @@ -810,30 +401,11 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, << " new lemmas during registration." << std::endl; return lems.size() + lemsPp.size(); } - Trace("nl-ext") << "We have " << d_ms.size() << " monomials." << std::endl; - - // register constants - registerMonomial(d_one); - for (unsigned j = 0; j < d_order_points.size(); j++) { - Node c = d_order_points[j]; - d_model.computeConcreteModelValue(c); - d_model.computeAbstractModelValue(c); - } - - // register variables - Trace("nl-ext-mv") << "Variables in monomials : " << std::endl; - for (unsigned i = 0; i < d_ms_vars.size(); i++) { - Node v = d_ms_vars[i]; - registerMonomial(v); - d_model.computeConcreteModelValue(v); - d_model.computeAbstractModelValue(v); - d_model.printModelValue("nl-ext-mv", v); - } //----------------------------------- possibly split on zero if (options::nlExtSplitZero()) { Trace("nl-ext") << "Get zero split lemmas..." << std::endl; - lemmas = checkSplitZero(); + lemmas = d_nlSlv.checkSplitZero(); filterLemmas(lemmas, lems); if (!lems.empty()) { @@ -854,7 +426,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, } //-----------------------------------lemmas based on sign (comparison to zero) - lemmas = checkMonomialSign(); + lemmas = d_nlSlv.checkMonomialSign(); filterLemmas(lemmas, lems); if (!lems.empty()) { @@ -874,23 +446,9 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, } //-----------------------------------lemmas based on magnitude of non-zero monomials - Trace("nl-ext-proc") << "Assign order ids..." << std::endl; - // sort by absolute values of abstract model values - assignOrderIds(d_ms_vars, d_order_vars, false, true); - - // sort individual variable lists - Trace("nl-ext-proc") << "Assign order var lists..." << std::endl; - SortNlModel smv; - smv.d_nlm = &d_model; - smv.d_isConcrete = false; - smv.d_isAbsolute = true; - smv.d_reverse_order = true; - for (unsigned j = 0; j < d_ms.size(); j++) { - std::sort(d_m_vlist[d_ms[j]].begin(), d_m_vlist[d_ms[j]].end(), smv); - } for (unsigned c = 0; c < 3; c++) { // c is effort level - lemmas = checkMonomialMagnitude( c ); + lemmas = d_nlSlv.checkMonomialMagnitude(c); unsigned nlem = lemmas.size(); filterLemmas(lemmas, lems); if (!lems.empty()) @@ -902,18 +460,11 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, } } - // sort monomials by degree - Trace("nl-ext-proc") << "Sort monomials by degree..." << std::endl; - SortNonlinearDegree snlad(d_m_degree); - std::sort(d_ms.begin(), d_ms.end(), snlad); - // all monomials - d_mterms.insert(d_mterms.end(), d_ms_vars.begin(), d_ms_vars.end()); - d_mterms.insert(d_mterms.end(), d_ms.begin(), d_ms.end()); - //-----------------------------------inferred bounds lemmas // e.g. x >= t => y*x >= y*t std::vector< Node > nt_lemmas; - lemmas = checkMonomialInferBounds(nt_lemmas, assertions, false_asserts); + lemmas = + d_nlSlv.checkMonomialInferBounds(nt_lemmas, assertions, false_asserts); // Trace("nl-ext") << "Bound lemmas : " << lemmas.size() << ", " << // nt_lemmas.size() << std::endl; prioritize lemmas that do not // introduce new monomials @@ -921,7 +472,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, if (options::nlExtTangentPlanes() && options::nlExtTangentPlanesInterleave()) { - lemmas = checkTangentPlanes(); + lemmas = d_nlSlv.checkTangentPlanes(); filterLemmas(lemmas, lems); } @@ -944,7 +495,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, //------------------------------------factoring lemmas // x*y + x*z >= t => exists k. k = y + z ^ x*k >= t if( options::nlExtFactor() ){ - lemmas = checkFactoring(assertions, false_asserts); + lemmas = d_nlSlv.checkFactoring(assertions, false_asserts); filterLemmas(lemmas, lems); if (!lems.empty()) { @@ -957,7 +508,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, //------------------------------------resolution bound inferences // e.g. ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t if (options::nlExtResBound()) { - lemmas = checkMonomialInferResBounds(); + lemmas = d_nlSlv.checkMonomialInferResBounds(); filterLemmas(lemmas, lems); if (!lems.empty()) { @@ -970,7 +521,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, //------------------------------------tangent planes if (options::nlExtTangentPlanes() && !options::nlExtTangentPlanesInterleave()) { - lemmas = checkTangentPlanes(); + lemmas = d_nlSlv.checkTangentPlanes(); filterLemmas(lemmas, wlems); } if (options::nlExtTfTangentPlanes()) @@ -1258,1187 +809,6 @@ void NonlinearExtension::presolve() Trace("nl-ext") << "NonlinearExtension::presolve" << std::endl; } -void NonlinearExtension::assignOrderIds(std::vector& vars, - NodeMultiset& order, - bool isConcrete, - bool isAbsolute) -{ - SortNlModel smv; - smv.d_nlm = &d_model; - smv.d_isConcrete = isConcrete; - smv.d_isAbsolute = isAbsolute; - smv.d_reverse_order = false; - std::sort(vars.begin(), vars.end(), smv); - - order.clear(); - // assign ordering id's - unsigned counter = 0; - unsigned order_index = isConcrete ? 0 : 1; - Node prev; - for (unsigned j = 0; j < vars.size(); j++) { - Node x = vars[j]; - Node v = d_model.computeModelValue(x, isConcrete); - if( !v.isConst() ){ - Trace("nl-ext-mvo") << "..do not assign order to " << x << " : " << v << std::endl; - //don't assign for non-constant values (transcendental function apps) - break; - } - Trace("nl-ext-mvo") << " order " << x << " : " << v << std::endl; - if (v != prev) { - // builtin points - bool success; - do { - success = false; - if (order_index < d_order_points.size()) { - Node vv = d_model.computeModelValue(d_order_points[order_index], - isConcrete); - if (d_model.compareValue(v, vv, isAbsolute) <= 0) - { - counter++; - Trace("nl-ext-mvo") << "O[" << d_order_points[order_index] - << "] = " << counter << std::endl; - order[d_order_points[order_index]] = counter; - prev = vv; - order_index++; - success = true; - } - } - } while (success); - } - if (prev.isNull() || d_model.compareValue(v, prev, isAbsolute) != 0) - { - counter++; - } - Trace("nl-ext-mvo") << "O[" << x << "] = " << counter << std::endl; - order[x] = counter; - prev = v; - } - while (order_index < d_order_points.size()) { - counter++; - Trace("nl-ext-mvo") << "O[" << d_order_points[order_index] - << "] = " << counter << std::endl; - order[d_order_points[order_index]] = counter; - order_index++; - } -} - -bool NonlinearExtension::getApproximateSqrt(Node c, - Node& l, - Node& u, - unsigned iter) const -{ - Assert(c.isConst()); - if (c == d_one || c == d_zero) - { - l = c; - u = c; - return true; - } - Rational rc = c.getConst(); - - Rational rl = rc < Rational(1) ? rc : Rational(1); - Rational ru = rc < Rational(1) ? Rational(1) : rc; - unsigned count = 0; - Rational half = Rational(1) / Rational(2); - while (count < iter) - { - Rational curr = half * (rl + ru); - Rational curr_sq = curr * curr; - if (curr_sq == rc) - { - rl = curr_sq; - ru = curr_sq; - break; - } - else if (curr_sq < rc) - { - rl = curr; - } - else - { - ru = curr; - } - count++; - } - - NodeManager* nm = NodeManager::currentNM(); - l = nm->mkConst(rl); - u = nm->mkConst(ru); - return true; -} - -// show a <> 0 by inequalities between variables in monomial a w.r.t 0 -int NonlinearExtension::compareSign(Node oa, Node a, unsigned a_index, - int status, std::vector& exp, - std::vector& lem) { - Trace("nl-ext-debug") << "Process " << a << " at index " << a_index - << ", status is " << status << std::endl; - Node mvaoa = d_model.computeAbstractModelValue(oa); - if (a_index == d_m_vlist[a].size()) { - if (mvaoa.getConst().sgn() != status) - { - Node lemma = - safeConstructNary(AND, exp).impNode(mkLit(oa, d_zero, status * 2)); - lem.push_back(lemma); - } - return status; - } else { - Assert(a_index < d_m_vlist[a].size()); - Node av = d_m_vlist[a][a_index]; - unsigned aexp = d_m_exp[a][av]; - // take current sign in model - Node mvaav = d_model.computeAbstractModelValue(av); - int sgn = mvaav.getConst().sgn(); - Trace("nl-ext-debug") << "Process var " << av << "^" << aexp - << ", model sign = " << sgn << std::endl; - if (sgn == 0) { - if (mvaoa.getConst().sgn() != 0) - { - Node lemma = av.eqNode(d_zero).impNode(oa.eqNode(d_zero)); - lem.push_back(lemma); - } - return 0; - } else { - if (aexp % 2 == 0) { - exp.push_back(av.eqNode(d_zero).negate()); - return compareSign(oa, a, a_index + 1, status, exp, lem); - } else { - exp.push_back( - NodeManager::currentNM()->mkNode(sgn == 1 ? GT : LT, av, d_zero)); - return compareSign(oa, a, a_index + 1, status * sgn, exp, lem); - } - } - } -} - -bool NonlinearExtension::compareMonomial( - Node oa, Node a, NodeMultiset& a_exp_proc, Node ob, Node b, - NodeMultiset& b_exp_proc, std::vector& exp, std::vector& lem, - std::map > >& cmp_infers) { - Trace("nl-ext-comp-debug") - << "Check |" << a << "| >= |" << b << "|" << std::endl; - unsigned pexp_size = exp.size(); - if (compareMonomial(oa, a, 0, a_exp_proc, ob, b, 0, b_exp_proc, 0, exp, lem, - cmp_infers)) { - return true; - } else { - exp.resize(pexp_size); - Trace("nl-ext-comp-debug") - << "Check |" << b << "| >= |" << a << "|" << std::endl; - if (compareMonomial(ob, b, 0, b_exp_proc, oa, a, 0, a_exp_proc, 0, exp, lem, - cmp_infers)) { - return true; - } else { - return false; - } - } -} - -bool NonlinearExtension::cmp_holds( - Node x, Node y, std::map >& cmp_infers, - std::vector& exp, std::map& visited) { - if (x == y) { - return true; - } else if (visited.find(x) == visited.end()) { - visited[x] = true; - std::map >::iterator it = cmp_infers.find(x); - if (it != cmp_infers.end()) { - for (std::map::iterator itc = it->second.begin(); - itc != it->second.end(); ++itc) { - exp.push_back(itc->second); - if (cmp_holds(itc->first, y, cmp_infers, exp, visited)) { - return true; - } - exp.pop_back(); - } - } - } - return false; -} - -// trying to show a ( >, = ) b by inequalities between variables in -// monomials a,b -bool NonlinearExtension::compareMonomial( - Node oa, Node a, unsigned a_index, NodeMultiset& a_exp_proc, Node ob, - Node b, unsigned b_index, NodeMultiset& b_exp_proc, int status, - std::vector& exp, std::vector& lem, - std::map > >& cmp_infers) { - Trace("nl-ext-comp-debug") - << "compareMonomial " << oa << " and " << ob << ", indices = " << a_index - << " " << b_index << std::endl; - Assert(status == 0 || status == 2); - if (a_index == d_m_vlist[a].size() && b_index == d_m_vlist[b].size()) { - // finished, compare absolute value of abstract model values - int modelStatus = d_model.compare(oa, ob, false, true) * -2; - Trace("nl-ext-comp") << "...finished comparison with " << oa << " <" - << status << "> " << ob - << ", model status = " << modelStatus << std::endl; - if (status != modelStatus) { - Trace("nl-ext-comp-infer") - << "infer : " << oa << " <" << status << "> " << ob << std::endl; - if (status == 2) { - // must state that all variables are non-zero - for (unsigned j = 0; j < d_m_vlist[a].size(); j++) { - exp.push_back(d_m_vlist[a][j].eqNode(d_zero).negate()); - } - } - Node clem = NodeManager::currentNM()->mkNode( - IMPLIES, safeConstructNary(AND, exp), mkLit(oa, ob, status, true)); - Trace("nl-ext-comp-lemma") << "comparison lemma : " << clem << std::endl; - lem.push_back(clem); - cmp_infers[status][oa][ob] = clem; - } - return true; - } else { - // get a/b variable information - Node av; - unsigned aexp = 0; - unsigned avo = 0; - if (a_index < d_m_vlist[a].size()) { - av = d_m_vlist[a][a_index]; - Assert(a_exp_proc[av] <= d_m_exp[a][av]); - aexp = d_m_exp[a][av] - a_exp_proc[av]; - if (aexp == 0) { - return compareMonomial(oa, a, a_index + 1, a_exp_proc, ob, b, b_index, - b_exp_proc, status, exp, lem, cmp_infers); - } - Assert(d_order_vars.find(av) != d_order_vars.end()); - avo = d_order_vars[av]; - } - Node bv; - unsigned bexp = 0; - unsigned bvo = 0; - if (b_index < d_m_vlist[b].size()) { - bv = d_m_vlist[b][b_index]; - Assert(b_exp_proc[bv] <= d_m_exp[b][bv]); - bexp = d_m_exp[b][bv] - b_exp_proc[bv]; - if (bexp == 0) { - return compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index + 1, - b_exp_proc, status, exp, lem, cmp_infers); - } - Assert(d_order_vars.find(bv) != d_order_vars.end()); - bvo = d_order_vars[bv]; - } - // get "one" information - Assert(d_order_vars.find(d_one) != d_order_vars.end()); - unsigned ovo = d_order_vars[d_one]; - Trace("nl-ext-comp-debug") << "....vars : " << av << "^" << aexp << " " - << bv << "^" << bexp << std::endl; - - //--- cases - if (av.isNull()) { - if (bvo <= ovo) { - Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl; - // can multiply b by <=1 - exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true)); - return compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index + 1, - b_exp_proc, bvo == ovo ? status : 2, exp, lem, - cmp_infers); - } else { - Trace("nl-ext-comp-debug") - << "...failure, unmatched |b|>1 component." << std::endl; - return false; - } - } else if (bv.isNull()) { - if (avo >= ovo) { - Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl; - // can multiply a by >=1 - exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true)); - return compareMonomial(oa, a, a_index + 1, a_exp_proc, ob, b, b_index, - b_exp_proc, avo == ovo ? status : 2, exp, lem, - cmp_infers); - } else { - Trace("nl-ext-comp-debug") - << "...failure, unmatched |a|<1 component." << std::endl; - return false; - } - } else { - Assert(!av.isNull() && !bv.isNull()); - if (avo >= bvo) { - if (bvo < ovo && avo >= ovo) { - Trace("nl-ext-comp-debug") << "...take leading " << av << std::endl; - // do avo>=1 instead - exp.push_back(mkLit(av, d_one, avo == ovo ? 0 : 2, true)); - return compareMonomial(oa, a, a_index + 1, a_exp_proc, ob, b, b_index, - b_exp_proc, avo == ovo ? status : 2, exp, lem, - cmp_infers); - } else { - unsigned min_exp = aexp > bexp ? bexp : aexp; - a_exp_proc[av] += min_exp; - b_exp_proc[bv] += min_exp; - Trace("nl-ext-comp-debug") - << "...take leading " << min_exp << " from " << av << " and " - << bv << std::endl; - exp.push_back(mkLit(av, bv, avo == bvo ? 0 : 2, true)); - bool ret = compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index, - b_exp_proc, avo == bvo ? status : 2, exp, - lem, cmp_infers); - a_exp_proc[av] -= min_exp; - b_exp_proc[bv] -= min_exp; - return ret; - } - } else { - if (bvo <= ovo) { - Trace("nl-ext-comp-debug") << "...take leading " << bv << std::endl; - // try multiply b <= 1 - exp.push_back(mkLit(d_one, bv, bvo == ovo ? 0 : 2, true)); - return compareMonomial(oa, a, a_index, a_exp_proc, ob, b, b_index + 1, - b_exp_proc, bvo == ovo ? status : 2, exp, lem, - cmp_infers); - } else { - Trace("nl-ext-comp-debug") - << "...failure, leading |b|>|a|>1 component." << std::endl; - return false; - } - } - } - } -} - -// status 0 : n equal, -1 : n superset, 1 : n subset -void NonlinearExtension::MonomialIndex::addTerm(Node n, - const std::vector& reps, - NonlinearExtension* nla, - int status, unsigned argIndex) { - if (status == 0) { - if (argIndex == reps.size()) { - d_monos.push_back(n); - } else { - d_data[reps[argIndex]].addTerm(n, reps, nla, status, argIndex + 1); - } - } - for (std::map::iterator it = d_data.begin(); - it != d_data.end(); ++it) { - if (status != 0 || argIndex == reps.size() || it->first != reps[argIndex]) { - // if we do not contain this variable, then if we were a superset, - // fail (-2), otherwise we are subset. if we do contain this - // variable, then if we were equal, we are superset since variables - // are ordered, otherwise we remain the same. - int new_status = - std::find(reps.begin(), reps.end(), it->first) == reps.end() - ? (status >= 0 ? 1 : -2) - : (status == 0 ? -1 : status); - if (new_status != -2) { - it->second.addTerm(n, reps, nla, new_status, argIndex); - } - } - } - // compare for subsets - for (unsigned i = 0; i < d_monos.size(); i++) { - Node m = d_monos[i]; - if (m != n) { - // we are superset if we are equal and haven't traversed all variables - int cstatus = status == 0 ? (argIndex == reps.size() ? 0 : -1) : status; - Trace("nl-ext-mindex-debug") << " compare " << n << " and " << m - << ", status = " << cstatus << std::endl; - if (cstatus <= 0 && nla->isMonomialSubset(m, n)) { - nla->registerMonomialSubset(m, n); - Trace("nl-ext-mindex-debug") << "...success" << std::endl; - } else if (cstatus >= 0 && nla->isMonomialSubset(n, m)) { - nla->registerMonomialSubset(n, m); - Trace("nl-ext-mindex-debug") << "...success (rev)" << std::endl; - } - } - } -} - -std::vector NonlinearExtension::checkMonomialSign() { - std::vector lemmas; - std::map signs; - Trace("nl-ext") << "Get monomial sign lemmas..." << std::endl; - for (unsigned j = 0; j < d_ms.size(); j++) { - Node a = d_ms[j]; - if (d_ms_proc.find(a) == d_ms_proc.end()) { - std::vector exp; - if (Trace.isOn("nl-ext-debug")) - { - Node cmva = d_model.computeConcreteModelValue(a); - Trace("nl-ext-debug") - << " process " << a << ", mv=" << cmva << "..." << std::endl; - } - if( d_m_nconst_factor.find( a )==d_m_nconst_factor.end() ){ - signs[a] = compareSign(a, a, 0, 1, exp, lemmas); - if (signs[a] == 0) { - d_ms_proc[a] = true; - Trace("nl-ext-debug") << "...mark " << a - << " reduced since its value is 0." << std::endl; - } - }else{ - Trace("nl-ext-debug") << "...can't conclude sign lemma for " << a << " since model value of a factor is non-constant." << std::endl; - //TODO - } - } - } - return lemmas; -} - -std::vector NonlinearExtension::checkMonomialMagnitude( unsigned c ) { - unsigned r = 1; - std::vector lemmas; -// if (x,y,L) in cmp_infers, then x > y inferred as conclusion of L - // in lemmas - std::map > > cmp_infers; - Trace("nl-ext") << "Get monomial comparison lemmas (order=" << r - << ", compare=" << c << ")..." << std::endl; - for (unsigned j = 0; j < d_ms.size(); j++) { - Node a = d_ms[j]; - if (d_ms_proc.find(a) == d_ms_proc.end() && - d_m_nconst_factor.find( a )==d_m_nconst_factor.end()) { - if (c == 0) { - // compare magnitude against 1 - std::vector exp; - NodeMultiset a_exp_proc; - NodeMultiset b_exp_proc; - compareMonomial(a, a, a_exp_proc, d_one, d_one, b_exp_proc, exp, - lemmas, cmp_infers); - } else { - std::map::iterator itmea = d_m_exp.find(a); - Assert(itmea != d_m_exp.end()); - if (c == 1) { - // TODO : not just against containing variables? - // compare magnitude against variables - for (unsigned k = 0; k < d_ms_vars.size(); k++) { - Node v = d_ms_vars[k]; - Node mvcv = d_model.computeConcreteModelValue(v); - if (mvcv.isConst()) - { - std::vector exp; - NodeMultiset a_exp_proc; - NodeMultiset b_exp_proc; - if (itmea->second.find(v) != itmea->second.end()) { - a_exp_proc[v] = 1; - b_exp_proc[v] = 1; - setMonomialFactor(a, v, a_exp_proc); - setMonomialFactor(v, a, b_exp_proc); - compareMonomial(a, a, a_exp_proc, v, v, b_exp_proc, exp, - lemmas, cmp_infers); - } - } - } - } else { - // compare magnitude against other non-linear monomials - for (unsigned k = (j + 1); k < d_ms.size(); k++) { - Node b = d_ms[k]; - //(signs[a]==signs[b])==(r==0) - if (d_ms_proc.find(b) == d_ms_proc.end() && - d_m_nconst_factor.find( b )==d_m_nconst_factor.end()) { - std::map::iterator itmeb = - d_m_exp.find(b); - Assert(itmeb != d_m_exp.end()); - - std::vector exp; - // take common factors of monomials, set minimum of - // common exponents as processed - NodeMultiset a_exp_proc; - NodeMultiset b_exp_proc; - for (NodeMultiset::iterator itmea2 = itmea->second.begin(); - itmea2 != itmea->second.end(); ++itmea2) { - NodeMultiset::iterator itmeb2 = - itmeb->second.find(itmea2->first); - if (itmeb2 != itmeb->second.end()) { - unsigned min_exp = itmea2->second > itmeb2->second - ? itmeb2->second - : itmea2->second; - a_exp_proc[itmea2->first] = min_exp; - b_exp_proc[itmea2->first] = min_exp; - Trace("nl-ext-comp") - << "Common exponent : " << itmea2->first << " : " - << min_exp << std::endl; - } - } - if (!a_exp_proc.empty()) { - setMonomialFactor(a, b, a_exp_proc); - setMonomialFactor(b, a, b_exp_proc); - } - /* - if( !a_exp_proc.empty() ){ - //reduction based on common exponents a > 0 => ( a * b - <> a * c <=> b <> c ), a < 0 => ( a * b <> a * c <=> b - !<> c ) ? }else{ compareMonomial( a, a, a_exp_proc, b, - b, b_exp_proc, exp, lemmas ); - } - */ - compareMonomial(a, a, a_exp_proc, b, b, b_exp_proc, exp, - lemmas, cmp_infers); - } - } - } - } - } - } - // remove redundant lemmas, e.g. if a > b, b > c, a > c were - // inferred, discard lemma with conclusion a > c - Trace("nl-ext-comp") << "Compute redundancies for " << lemmas.size() - << " lemmas." << std::endl; - // naive - std::vector r_lemmas; - for (std::map > >::iterator itb = - cmp_infers.begin(); - itb != cmp_infers.end(); ++itb) { - for (std::map >::iterator itc = - itb->second.begin(); - itc != itb->second.end(); ++itc) { - for (std::map::iterator itc2 = itc->second.begin(); - itc2 != itc->second.end(); ++itc2) { - std::map visited; - for (std::map::iterator itc3 = itc->second.begin(); - itc3 != itc->second.end(); ++itc3) { - if (itc3->first != itc2->first) { - std::vector exp; - if (cmp_holds(itc3->first, itc2->first, itb->second, exp, - visited)) { - r_lemmas.push_back(itc2->second); - Trace("nl-ext-comp") - << "...inference of " << itc->first << " > " - << itc2->first << " was redundant." << std::endl; - break; - } - } - } - } - } - } - std::vector nr_lemmas; - for (unsigned i = 0; i < lemmas.size(); i++) { - if (std::find(r_lemmas.begin(), r_lemmas.end(), lemmas[i]) == - r_lemmas.end()) { - nr_lemmas.push_back(lemmas[i]); - } - } - // TODO: only take maximal lower/minimial lower bounds? - - Trace("nl-ext-comp") << nr_lemmas.size() << " / " << lemmas.size() - << " were non-redundant." << std::endl; - return nr_lemmas; -} - -std::vector NonlinearExtension::checkTangentPlanes() { - std::vector< Node > lemmas; - Trace("nl-ext") << "Get monomial tangent plane lemmas..." << std::endl; - NodeManager* nm = NodeManager::currentNM(); - unsigned kstart = d_ms_vars.size(); - for (unsigned k = kstart; k < d_mterms.size(); k++) { - Node t = d_mterms[k]; - // if this term requires a refinement - if (d_tplane_refine.find(t) != d_tplane_refine.end()) - { - Trace("nl-ext-tplanes") - << "Look at monomial requiring refinement : " << t << std::endl; - // get a decomposition - std::map >::iterator it = - d_m_contain_children.find(t); - if (it != d_m_contain_children.end()) { - std::map > dproc; - for (unsigned j = 0; j < it->second.size(); j++) { - Node tc = it->second[j]; - if (tc != d_one) { - Node tc_diff = d_m_contain_umult[tc][t]; - Assert(!tc_diff.isNull()); - Node a = tc < tc_diff ? tc : tc_diff; - Node b = tc < tc_diff ? tc_diff : tc; - if (dproc[a].find(b) == dproc[a].end()) { - dproc[a][b] = true; - Trace("nl-ext-tplanes") - << " decomposable into : " << a << " * " << b << std::endl; - Node a_v_c = d_model.computeAbstractModelValue(a); - Node b_v_c = d_model.computeAbstractModelValue(b); - // points we will add tangent planes for - std::vector< Node > pts[2]; - pts[0].push_back( a_v_c ); - pts[1].push_back( b_v_c ); - // if previously refined - bool prevRefine = d_tangent_val_bound[0][a].find( b )!=d_tangent_val_bound[0][a].end(); - // a_min, a_max, b_min, b_max - for( unsigned p=0; p<4; p++ ){ - Node curr_v = p<=1 ? a_v_c : b_v_c; - if( prevRefine ){ - Node pt_v = d_tangent_val_bound[p][a][b]; - Assert(!pt_v.isNull()); - if( curr_v!=pt_v ){ - Node do_extend = - nm->mkNode((p == 1 || p == 3) ? GT : LT, curr_v, pt_v); - do_extend = Rewriter::rewrite( do_extend ); - if( do_extend==d_true ){ - for( unsigned q=0; q<2; q++ ){ - pts[ p<=1 ? 0 : 1 ].push_back( curr_v ); - pts[ p<=1 ? 1 : 0 ].push_back( d_tangent_val_bound[ p<=1 ? 2+q : q ][a][b] ); - } - } - } - }else{ - d_tangent_val_bound[p][a][b] = curr_v; - } - } - - for( unsigned p=0; pmkNode(MINUS, - nm->mkNode(PLUS, - nm->mkNode(MULT, b_v, a), - nm->mkNode(MULT, a_v, b)), - nm->mkNode(MULT, a_v, b_v)); - for (unsigned d = 0; d < 4; d++) { - Node aa = nm->mkNode(d == 0 || d == 3 ? GEQ : LEQ, a, a_v); - Node ab = nm->mkNode(d == 1 || d == 3 ? GEQ : LEQ, b, b_v); - Node conc = nm->mkNode(d <= 1 ? LEQ : GEQ, t, tplane); - Node tlem = nm->mkNode(OR, aa.negate(), ab.negate(), conc); - Trace("nl-ext-tplanes") - << "Tangent plane lemma : " << tlem << std::endl; - lemmas.push_back(tlem); - } - - // tangent plane reverse implication - - // t <= tplane -> ( (a <= a_v ^ b >= b_v) v - // (a >= a_v ^ b <= b_v) ). - // in clause form, the above becomes - // t <= tplane -> a <= a_v v b <= b_v. - // t <= tplane -> b >= b_v v a >= a_v. - Node a_leq_av = nm->mkNode(LEQ, a, a_v); - Node b_leq_bv = nm->mkNode(LEQ, b, b_v); - Node a_geq_av = nm->mkNode(GEQ, a, a_v); - Node b_geq_bv = nm->mkNode(GEQ, b, b_v); - - Node t_leq_tplane = nm->mkNode(LEQ, t, tplane); - Node a_leq_av_or_b_leq_bv = nm->mkNode(OR, a_leq_av, b_leq_bv); - Node b_geq_bv_or_a_geq_av = nm->mkNode(OR, b_geq_bv, a_geq_av); - Node ub_reverse1 = - nm->mkNode(OR, t_leq_tplane.negate(), a_leq_av_or_b_leq_bv); - Trace("nl-ext-tplanes") - << "Tangent plane lemma (reverse) : " << ub_reverse1 - << std::endl; - lemmas.push_back(ub_reverse1); - Node ub_reverse2 = - nm->mkNode(OR, t_leq_tplane.negate(), b_geq_bv_or_a_geq_av); - Trace("nl-ext-tplanes") - << "Tangent plane lemma (reverse) : " << ub_reverse2 - << std::endl; - lemmas.push_back(ub_reverse2); - - // t >= tplane -> ( (a <= a_v ^ b <= b_v) v - // (a >= a_v ^ b >= b_v) ). - // in clause form, the above becomes - // t >= tplane -> a <= a_v v b >= b_v. - // t >= tplane -> b >= b_v v a <= a_v - Node t_geq_tplane = nm->mkNode(GEQ, t, tplane); - Node a_leq_av_or_b_geq_bv = nm->mkNode(OR, a_leq_av, b_geq_bv); - Node a_geq_av_or_b_leq_bv = nm->mkNode(OR, a_geq_av, b_leq_bv); - Node lb_reverse1 = - nm->mkNode(OR, t_geq_tplane.negate(), a_leq_av_or_b_geq_bv); - Trace("nl-ext-tplanes") - << "Tangent plane lemma (reverse) : " << lb_reverse1 - << std::endl; - lemmas.push_back(lb_reverse1); - Node lb_reverse2 = - nm->mkNode(OR, t_geq_tplane.negate(), a_geq_av_or_b_leq_bv); - Trace("nl-ext-tplanes") - << "Tangent plane lemma (reverse) : " << lb_reverse2 - << std::endl; - lemmas.push_back(lb_reverse2); - } - } - } - } - } - } - } - Trace("nl-ext") << "...trying " << lemmas.size() - << " tangent plane lemmas..." << std::endl; - return lemmas; -} - -std::vector NonlinearExtension::checkMonomialInferBounds( - std::vector& nt_lemmas, - const std::vector& asserts, - const std::vector& false_asserts) -{ - std::vector< Node > lemmas; - // register constraints - Trace("nl-ext-debug") << "Register bound constraints..." << std::endl; - for (const Node& lit : asserts) - { - bool polarity = lit.getKind() != NOT; - Node atom = lit.getKind() == NOT ? lit[0] : lit; - registerConstraint(atom); - bool is_false_lit = - std::find(false_asserts.begin(), false_asserts.end(), lit) - != false_asserts.end(); - // add information about bounds to variables - std::map >::iterator itc = - d_c_info.find(atom); - std::map >::iterator itcm = - d_c_info_maxm.find(atom); - if (itc != d_c_info.end()) { - Assert(itcm != d_c_info_maxm.end()); - for (std::map::iterator itcc = itc->second.begin(); - itcc != itc->second.end(); ++itcc) { - Node x = itcc->first; - Node coeff = itcc->second.d_coeff; - Node rhs = itcc->second.d_rhs; - Kind type = itcc->second.d_type; - Node exp = lit; - if (!polarity) { - // reverse - if (type == EQUAL) - { - // we will take the strict inequality in the direction of the - // model - Node lhs = ArithMSum::mkCoeffTerm(coeff, x); - Node query = NodeManager::currentNM()->mkNode(GT, lhs, rhs); - Node query_mv = d_model.computeAbstractModelValue(query); - if (query_mv == d_true) { - exp = query; - type = GT; - } else { - Assert(query_mv == d_false); - exp = NodeManager::currentNM()->mkNode(LT, lhs, rhs); - type = LT; - } - } else { - type = negateKind(type); - } - } - // add to status if maximal degree - d_ci_max[x][coeff][rhs] = itcm->second.find(x) != itcm->second.end(); - if (Trace.isOn("nl-ext-bound-debug2")) { - Node t = ArithMSum::mkCoeffTerm(coeff, x); - Trace("nl-ext-bound-debug2") - << "Add Bound: " << t << " " << type << " " << rhs << " by " - << exp << std::endl; - } - bool updated = true; - std::map::iterator its = d_ci[x][coeff].find(rhs); - if (its == d_ci[x][coeff].end()) { - d_ci[x][coeff][rhs] = type; - d_ci_exp[x][coeff][rhs] = exp; - } else if (type != its->second) { - Trace("nl-ext-bound-debug2") - << "Joining kinds : " << type << " " << its->second << std::endl; - Kind jk = joinKinds(type, its->second); - if (jk == UNDEFINED_KIND) - { - updated = false; - } else if (jk != its->second) { - if (jk == type) { - d_ci[x][coeff][rhs] = type; - d_ci_exp[x][coeff][rhs] = exp; - } else { - d_ci[x][coeff][rhs] = jk; - d_ci_exp[x][coeff][rhs] = NodeManager::currentNM()->mkNode( - AND, d_ci_exp[x][coeff][rhs], exp); - } - } else { - updated = false; - } - } - if (Trace.isOn("nl-ext-bound")) { - if (updated) { - Trace("nl-ext-bound") << "Bound: "; - debugPrintBound("nl-ext-bound", coeff, x, d_ci[x][coeff][rhs], rhs); - Trace("nl-ext-bound") << " by " << d_ci_exp[x][coeff][rhs]; - if (d_ci_max[x][coeff][rhs]) { - Trace("nl-ext-bound") << ", is max degree"; - } - Trace("nl-ext-bound") << std::endl; - } - } - // compute if bound is not satisfied, and store what is required - // for a possible refinement - if (options::nlExtTangentPlanes()) { - if (is_false_lit) { - d_tplane_refine.insert(x); - } - } - } - } - } - // reflexive constraints - Node null_coeff; - for (unsigned j = 0; j < d_mterms.size(); j++) { - Node n = d_mterms[j]; - d_ci[n][null_coeff][n] = EQUAL; - d_ci_exp[n][null_coeff][n] = d_true; - d_ci_max[n][null_coeff][n] = false; - } - - Trace("nl-ext") << "Get inferred bound lemmas..." << std::endl; - - for (unsigned k = 0; k < d_mterms.size(); k++) { - Node x = d_mterms[k]; - Trace("nl-ext-bound-debug") - << "Process bounds for " << x << " : " << std::endl; - std::map >::iterator itm = - d_m_contain_parent.find(x); - if (itm != d_m_contain_parent.end()) { - Trace("nl-ext-bound-debug") << "...has " << itm->second.size() - << " parent monomials." << std::endl; - // check derived bounds - std::map > >::iterator itc = - d_ci.find(x); - if (itc != d_ci.end()) { - for (std::map >::iterator itcc = - itc->second.begin(); - itcc != itc->second.end(); ++itcc) { - Node coeff = itcc->first; - Node t = ArithMSum::mkCoeffTerm(coeff, x); - for (std::map::iterator itcr = itcc->second.begin(); - itcr != itcc->second.end(); ++itcr) { - Node rhs = itcr->first; - // only consider this bound if maximal degree - if (d_ci_max[x][coeff][rhs]) { - Kind type = itcr->second; - for (unsigned j = 0; j < itm->second.size(); j++) { - Node y = itm->second[j]; - Assert(d_m_contain_mult[x].find(y) - != d_m_contain_mult[x].end()); - Node mult = d_m_contain_mult[x][y]; - // x t => m*x t where y = m*x - // get the sign of mult - Node mmv = d_model.computeConcreteModelValue(mult); - Trace("nl-ext-bound-debug2") - << "Model value of " << mult << " is " << mmv << std::endl; - if(mmv.isConst()){ - int mmv_sign = mmv.getConst().sgn(); - Trace("nl-ext-bound-debug2") - << " sign of " << mmv << " is " << mmv_sign << std::endl; - if (mmv_sign != 0) { - Trace("nl-ext-bound-debug") - << " from " << x << " * " << mult << " = " << y - << " and " << t << " " << type << " " << rhs - << ", infer : " << std::endl; - Kind infer_type = - mmv_sign == -1 ? reverseRelationKind(type) : type; - Node infer_lhs = - NodeManager::currentNM()->mkNode(MULT, mult, t); - Node infer_rhs = - NodeManager::currentNM()->mkNode(MULT, mult, rhs); - Node infer = NodeManager::currentNM()->mkNode( - infer_type, infer_lhs, infer_rhs); - Trace("nl-ext-bound-debug") << " " << infer << std::endl; - infer = Rewriter::rewrite(infer); - Trace("nl-ext-bound-debug2") - << " ...rewritten : " << infer << std::endl; - // check whether it is false in model for abstraction - Node infer_mv = d_model.computeAbstractModelValue(infer); - Trace("nl-ext-bound-debug") - << " ...infer model value is " << infer_mv - << std::endl; - if (infer_mv == d_false) { - Node exp = NodeManager::currentNM()->mkNode( - AND, - NodeManager::currentNM()->mkNode( - mmv_sign == 1 ? GT : LT, mult, d_zero), - d_ci_exp[x][coeff][rhs]); - Node iblem = - NodeManager::currentNM()->mkNode(IMPLIES, exp, infer); - Node pr_iblem = iblem; - iblem = Rewriter::rewrite(iblem); - bool introNewTerms = hasNewMonomials(iblem, d_ms); - Trace("nl-ext-bound-lemma") - << "*** Bound inference lemma : " << iblem - << " (pre-rewrite : " << pr_iblem << ")" << std::endl; - // Trace("nl-ext-bound-lemma") << " intro new - // monomials = " << introNewTerms << std::endl; - if (!introNewTerms) { - lemmas.push_back(iblem); - } else { - nt_lemmas.push_back(iblem); - } - } - } else { - Trace("nl-ext-bound-debug") << " ...coefficient " << mult - << " is zero." << std::endl; - } - }else{ - Trace("nl-ext-bound-debug") << " ...coefficient " << mult - << " is non-constant (probably transcendental)." << std::endl; - } - } - } - } - } - } - } else { - Trace("nl-ext-bound-debug") << "...has no parent monomials." << std::endl; - } - } - return lemmas; -} - -std::vector NonlinearExtension::checkFactoring( - const std::vector& asserts, const std::vector& false_asserts) -{ - std::vector< Node > lemmas; - Trace("nl-ext") << "Get factoring lemmas..." << std::endl; - for (const Node& lit : asserts) - { - bool polarity = lit.getKind() != NOT; - Node atom = lit.getKind() == NOT ? lit[0] : lit; - Node litv = d_model.computeConcreteModelValue(lit); - bool considerLit = false; - // Only consider literals that are in false_asserts. - considerLit = std::find(false_asserts.begin(), false_asserts.end(), lit) - != false_asserts.end(); - - if (considerLit) - { - std::map msum; - if (ArithMSum::getMonomialSumLit(atom, msum)) - { - Trace("nl-ext-factor") << "Factoring for literal " << lit << ", monomial sum is : " << std::endl; - if (Trace.isOn("nl-ext-factor")) { - ArithMSum::debugPrintMonomialSum(msum, "nl-ext-factor"); - } - std::map< Node, std::vector< Node > > factor_to_mono; - std::map< Node, std::vector< Node > > factor_to_mono_orig; - for( std::map::iterator itm = msum.begin(); itm != msum.end(); ++itm ){ - if( !itm->first.isNull() ){ - if( itm->first.getKind()==NONLINEAR_MULT ){ - std::vector< Node > children; - for( unsigned i=0; ifirst.getNumChildren(); i++ ){ - children.push_back( itm->first[i] ); - } - std::map< Node, bool > processed; - for( unsigned i=0; ifirst.getNumChildren(); i++ ){ - if( processed.find( itm->first[i] )==processed.end() ){ - processed[itm->first[i]] = true; - children[i] = d_one; - if( !itm->second.isNull() ){ - children.push_back( itm->second ); - } - Node val = NodeManager::currentNM()->mkNode(MULT, children); - if( !itm->second.isNull() ){ - children.pop_back(); - } - children[i] = itm->first[i]; - val = Rewriter::rewrite( val ); - factor_to_mono[itm->first[i]].push_back( val ); - factor_to_mono_orig[itm->first[i]].push_back( itm->first ); - } - } - } - } - } - for( std::map< Node, std::vector< Node > >::iterator itf = factor_to_mono.begin(); itf != factor_to_mono.end(); ++itf ){ - Node x = itf->first; - if (itf->second.size() == 1) - { - std::map::iterator itm = msum.find(x); - if (itm != msum.end()) - { - itf->second.push_back(itm->second.isNull() ? d_one : itm->second); - factor_to_mono_orig[x].push_back(x); - } - } - if( itf->second.size()>1 ){ - Node sum = NodeManager::currentNM()->mkNode(PLUS, itf->second); - sum = Rewriter::rewrite( sum ); - Trace("nl-ext-factor") - << "* Factored sum for " << x << " : " << sum << std::endl; - Node kf = getFactorSkolem(sum, lemmas); - std::vector< Node > poly; - poly.push_back(NodeManager::currentNM()->mkNode(MULT, x, kf)); - std::map >::iterator itfo = - factor_to_mono_orig.find(x); - Assert(itfo != factor_to_mono_orig.end()); - for( std::map::iterator itm = msum.begin(); itm != msum.end(); ++itm ){ - if( std::find( itfo->second.begin(), itfo->second.end(), itm->first )==itfo->second.end() ){ - poly.push_back(ArithMSum::mkCoeffTerm( - itm->second, itm->first.isNull() ? d_one : itm->first)); - } - } - Node polyn = poly.size() == 1 - ? poly[0] - : NodeManager::currentNM()->mkNode(PLUS, poly); - Trace("nl-ext-factor") << "...factored polynomial : " << polyn << std::endl; - Node conc_lit = NodeManager::currentNM()->mkNode( atom.getKind(), polyn, d_zero ); - conc_lit = Rewriter::rewrite( conc_lit ); - if( !polarity ){ - conc_lit = conc_lit.negate(); - } - - std::vector< Node > lemma_disj; - lemma_disj.push_back( lit.negate() ); - lemma_disj.push_back( conc_lit ); - Node flem = NodeManager::currentNM()->mkNode(OR, lemma_disj); - Trace("nl-ext-factor") << "...lemma is " << flem << std::endl; - lemmas.push_back( flem ); - } - } - } - } - } - return lemmas; -} - -Node NonlinearExtension::getFactorSkolem( Node n, std::vector< Node >& lemmas ) { - std::map< Node, Node >::iterator itf = d_factor_skolem.find( n ); - if( itf==d_factor_skolem.end() ){ - Node k = NodeManager::currentNM()->mkSkolem( "kf", n.getType() ); - Node k_eq = Rewriter::rewrite( k.eqNode( n ) ); - lemmas.push_back( k_eq ); - d_factor_skolem[n] = k; - return k; - }else{ - return itf->second; - } -} - -std::vector NonlinearExtension::checkMonomialInferResBounds() { - std::vector< Node > lemmas; - Trace("nl-ext") << "Get monomial resolution inferred bound lemmas..." << std::endl; - for (unsigned j = 0; j < d_mterms.size(); j++) { - Node a = d_mterms[j]; - std::map > >::iterator itca = - d_ci.find(a); - if (itca != d_ci.end()) { - for (unsigned k = (j + 1); k < d_mterms.size(); k++) { - Node b = d_mterms[k]; - std::map > >::iterator - itcb = d_ci.find(b); - if (itcb != d_ci.end()) { - Trace("nl-ext-rbound-debug") << "resolution inferences : compare " - << a << " and " << b << std::endl; - // if they have common factors - std::map::iterator ita = d_mono_diff[a].find(b); - if (ita != d_mono_diff[a].end()) { - Trace("nl-ext-rbound") << "Get resolution inferences for [a] " << a - << " vs [b] " << b << std::endl; - std::map::iterator itb = d_mono_diff[b].find(a); - Assert(itb != d_mono_diff[b].end()); - Node mv_a = d_model.computeAbstractModelValue(ita->second); - Assert(mv_a.isConst()); - int mv_a_sgn = mv_a.getConst().sgn(); - if (mv_a_sgn == 0) - { - // we don't compare monomials whose current model value is zero - continue; - } - Node mv_b = d_model.computeAbstractModelValue(itb->second); - Assert(mv_b.isConst()); - int mv_b_sgn = mv_b.getConst().sgn(); - if (mv_b_sgn == 0) - { - // we don't compare monomials whose current model value is zero - continue; - } - Trace("nl-ext-rbound") - << " [a] factor is " << ita->second - << ", sign in model = " << mv_a_sgn << std::endl; - Trace("nl-ext-rbound") - << " [b] factor is " << itb->second - << ", sign in model = " << mv_b_sgn << std::endl; - - std::vector exp; - // bounds of a - for (std::map >::iterator itcac = - itca->second.begin(); - itcac != itca->second.end(); ++itcac) { - Node coeff_a = itcac->first; - for (std::map::iterator itcar = - itcac->second.begin(); - itcar != itcac->second.end(); ++itcar) { - Node rhs_a = itcar->first; - Node rhs_a_res_base = - NodeManager::currentNM()->mkNode(MULT, itb->second, rhs_a); - rhs_a_res_base = Rewriter::rewrite(rhs_a_res_base); - if (!hasNewMonomials(rhs_a_res_base, d_ms)) { - Kind type_a = itcar->second; - exp.push_back(d_ci_exp[a][coeff_a][rhs_a]); - - // bounds of b - for (std::map >::iterator itcbc = - itcb->second.begin(); - itcbc != itcb->second.end(); ++itcbc) { - Node coeff_b = itcbc->first; - Node rhs_a_res = - ArithMSum::mkCoeffTerm(coeff_b, rhs_a_res_base); - for (std::map::iterator itcbr = - itcbc->second.begin(); - itcbr != itcbc->second.end(); ++itcbr) { - Node rhs_b = itcbr->first; - Node rhs_b_res = NodeManager::currentNM()->mkNode( - MULT, ita->second, rhs_b); - rhs_b_res = ArithMSum::mkCoeffTerm(coeff_a, rhs_b_res); - rhs_b_res = Rewriter::rewrite(rhs_b_res); - if (!hasNewMonomials(rhs_b_res, d_ms)) { - Kind type_b = itcbr->second; - exp.push_back(d_ci_exp[b][coeff_b][rhs_b]); - if (Trace.isOn("nl-ext-rbound")) { - Trace("nl-ext-rbound") << "* try bounds : "; - debugPrintBound("nl-ext-rbound", coeff_a, a, type_a, - rhs_a); - Trace("nl-ext-rbound") << std::endl; - Trace("nl-ext-rbound") << " "; - debugPrintBound("nl-ext-rbound", coeff_b, b, type_b, - rhs_b); - Trace("nl-ext-rbound") << std::endl; - } - Kind types[2]; - for (unsigned r = 0; r < 2; r++) { - Node pivot_factor = - r == 0 ? itb->second : ita->second; - int pivot_factor_sign = - r == 0 ? mv_b_sgn : mv_a_sgn; - types[r] = r == 0 ? type_a : type_b; - if (pivot_factor_sign == (r == 0 ? 1 : -1)) { - types[r] = reverseRelationKind(types[r]); - } - if (pivot_factor_sign == 1) { - exp.push_back(NodeManager::currentNM()->mkNode( - GT, pivot_factor, d_zero)); - } else { - exp.push_back(NodeManager::currentNM()->mkNode( - LT, pivot_factor, d_zero)); - } - } - Kind jk = transKinds(types[0], types[1]); - Trace("nl-ext-rbound-debug") - << "trans kind : " << types[0] << " + " - << types[1] << " = " << jk << std::endl; - if (jk != UNDEFINED_KIND) - { - Node conc = NodeManager::currentNM()->mkNode( - jk, rhs_a_res, rhs_b_res); - Node conc_mv = - d_model.computeAbstractModelValue(conc); - if (conc_mv == d_false) { - Node rblem = NodeManager::currentNM()->mkNode( - IMPLIES, - NodeManager::currentNM()->mkNode(AND, exp), - conc); - Trace("nl-ext-rbound-lemma-debug") - << "Resolution bound lemma " - "(pre-rewrite) " - ": " - << rblem << std::endl; - rblem = Rewriter::rewrite(rblem); - Trace("nl-ext-rbound-lemma") - << "Resolution bound lemma : " << rblem - << std::endl; - lemmas.push_back(rblem); - } - } - exp.pop_back(); - exp.pop_back(); - exp.pop_back(); - } - } - } - exp.pop_back(); - } - } - } - } - } - } - } - } - return lemmas; -} } // namespace arith } // namespace theory diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h index bfc713b12..5aa2070a6 100644 --- a/src/theory/arith/nonlinear_extension.h +++ b/src/theory/arith/nonlinear_extension.h @@ -19,23 +19,15 @@ #define CVC4__THEORY__ARITH__NONLINEAR_EXTENSION_H #include - #include -#include -#include -#include -#include #include -#include "context/cdhashset.h" -#include "context/cdinsert_hashmap.h" #include "context/cdlist.h" -#include "context/cdqueue.h" -#include "context/context.h" #include "expr/kind.h" #include "expr/node.h" #include "theory/arith/nl_lemma_utils.h" #include "theory/arith/nl_model.h" +#include "theory/arith/nl_solver.h" #include "theory/arith/theory_arith.h" #include "theory/arith/transcendental_solver.h" #include "theory/uf/equality_engine.h" @@ -44,9 +36,6 @@ namespace CVC4 { namespace theory { namespace arith { -typedef std::map NodeMultiset; - -// TODO : refactor/document this class (#1287) /** Non-linear extension class * * This class implements model-based refinement schemes @@ -72,6 +61,8 @@ typedef std::map NodeMultiset; * where d_out is the output channel of TheoryArith. */ class NonlinearExtension { + typedef context::CDHashSet NodeSet; + public: NonlinearExtension(TheoryArith& containing, eq::EqualityEngine* ee); ~NonlinearExtension(); @@ -187,34 +178,7 @@ class NonlinearExtension { bool modelBasedRefinement(std::vector& mlems, std::vector& mlemsPp, std::map& lemSE); - /** returns true if the multiset containing the - * factors of monomial a is a subset of the multiset - * containing the factors of monomial b. - */ - bool isMonomialSubset(Node a, Node b) const; - void registerMonomialSubset(Node a, Node b); - - typedef context::CDHashSet NodeSet; - - // monomial information (context-independent) - class MonomialIndex { - public: - // status 0 : n equal, -1 : n superset, 1 : n subset - void addTerm(Node n, const std::vector& reps, NonlinearExtension* nla, - int status = 0, unsigned argIndex = 0); - - private: - std::map d_data; - std::vector d_monos; - }; /* class MonomialIndex */ - // constraint information (context-independent) - struct ConstraintInfo { - public: - Node d_rhs; - Node d_coeff; - Kind d_type; - }; /* struct ConstraintInfo */ /** check last call * @@ -249,24 +213,6 @@ class NonlinearExtension { std::vector& lemsPp, std::vector& wlems, std::map& lemSE); - //---------------------------------------term utilities - static bool isArithKind(Kind k); - static Node mkLit(Node a, Node b, int status, bool isAbsolute = false); - static Node mkAbs(Node a); - static Node mkValidPhase(Node a, Node pi); - Node mkMonomialRemFactor(Node n, const NodeMultiset& n_exp_rem) const; - //---------------------------------------end term utilities - - /** register monomial */ - void registerMonomial(Node n); - void setMonomialFactor(Node a, Node b, const NodeMultiset& common); - - void registerConstraint(Node atom); - /** assign order ids */ - void assignOrderIds(std::vector& vars, - NodeMultiset& d_order, - bool isConcrete, - bool isAbsolute); /** get assertions * @@ -312,85 +258,6 @@ class NonlinearExtension { std::vector& gs); //---------------------------end check model - /** In the following functions, status states a relationship - * between two arithmetic terms, where: - * 0 : equal - * 1 : greater than or equal - * 2 : greater than - * -X : (greater -> less) - * TODO (#1287) make this an enum? - */ - /** compute the sign of a. - * - * Calls to this function are such that : - * exp => ( oa = a ^ a 0 ) - * - * This function iterates over the factors of a, - * where a_index is the index of the factor in a - * we are currently looking at. - * - * This function returns a status, which indicates - * a's relationship to 0. - * We add lemmas to lem of the form given by the - * lemma schema checkSign(...). - */ - int compareSign(Node oa, Node a, unsigned a_index, int status, - std::vector& exp, std::vector& lem); - /** compare monomials a and b - * - * Initially, a call to this function is such that : - * exp => ( oa = a ^ ob = b ) - * - * This function returns true if we can infer a valid - * arithmetic lemma of the form : - * P => abs( a ) >= abs( b ) - * where P is true and abs( a ) >= abs( b ) is false in the - * current model. - * - * This function is implemented by "processing" factors - * of monomials a and b until an inference of the above - * form can be made. For example, if : - * a = x*x*y and b = z*w - * Assuming we are trying to show abs( a ) >= abs( c ), - * then if abs( M( x ) ) >= abs( M( z ) ) where M is the current model, - * then we can add abs( x ) >= abs( z ) to our explanation, and - * mark one factor of x as processed in a, and - * one factor of z as processed in b. The number of processed factors of a - * and b are stored in a_exp_proc and b_exp_proc respectively. - * - * cmp_infers stores information that is helpful - * in discarding redundant inferences. For example, - * we do not want to infer abs( x ) >= abs( z ) if - * we have already inferred abs( x ) >= abs( y ) and - * abs( y ) >= abs( z ). - * It stores entries of the form (status,t1,t2)->F, - * which indicates that we constructed a lemma F that - * showed t1 t2. - * - * We add lemmas to lem of the form given by the - * lemma schema checkMagnitude(...). - */ - bool compareMonomial( - Node oa, Node a, NodeMultiset& a_exp_proc, Node ob, Node b, - NodeMultiset& b_exp_proc, std::vector& exp, std::vector& lem, - std::map > >& cmp_infers); - /** helper function for above - * - * The difference is the inputs a_index and b_index, which are the indices of - * children (factors) in monomials a and b which we are currently looking at. - */ - bool compareMonomial( - Node oa, Node a, unsigned a_index, NodeMultiset& a_exp_proc, Node ob, - Node b, unsigned b_index, NodeMultiset& b_exp_proc, int status, - std::vector& exp, std::vector& lem, - std::map > >& cmp_infers); - /** Check whether we have already inferred a relationship between monomials - * x and y based on the information in cmp_infers. This computes the - * transitive closure of the relation stored in cmp_infers. - */ - bool cmp_holds(Node x, Node y, - std::map >& cmp_infers, - std::vector& exp, std::map& visited); /** Is n entailed with polarity pol in the current context? */ bool isEntailed(Node n, bool pol); @@ -412,61 +279,19 @@ class NonlinearExtension { /** Process side effect se */ void processSideEffect(const NlLemmaSideEffect& se); - // Returns the NodeMultiset for an existing monomial. - const NodeMultiset& getMonomialExponentMap(Node monomial) const; - - // Map from monomials to var^index. - typedef std::map MonomialExponentMap; - MonomialExponentMap d_m_exp; - - /** - * Mapping from monomials to the list of variables that occur in it. For - * example, x*x*y*z -> { x, y, z }. - */ - std::map > d_m_vlist; - NodeMultiset d_m_degree; - // monomial index, by sorted variables - MonomialIndex d_m_index; - // list of all monomials - std::vector d_monomials; - // containment ordering - std::map > d_m_contain_children; - std::map > d_m_contain_parent; - std::map > d_m_contain_mult; - std::map > d_m_contain_umult; - // ( x*y, x*z, y ) for each pair of monomials ( x*y, x*z ) with common factors - std::map > d_mono_diff; - /** cache of all lemmas sent on the output channel (user-context-dependent) */ NodeSet d_lemmas; - /** cache of terms t for which we have added the lemma ( t = 0 V t != 0 ). */ - NodeSet d_zero_split; - /** commonly used terms */ Node d_zero; Node d_one; Node d_neg_one; - Node d_two; Node d_true; - Node d_false; - // The theory of arithmetic containing this extension. TheoryArith& d_containing; - // pointer to used equality engine eq::EqualityEngine* d_ee; // needs last call effort bool d_needsLastCall; - - // if d_c_info[lit][x] = ( r, coeff, k ), then ( lit <=> (coeff * x) r ) - std::map > d_c_info; - std::map > d_c_info_maxm; - std::vector d_constraints; - - // per last-call effort - - // model values/orderings - /** The non-linear model object * * This class is responsible for computing model values for arithmetic terms @@ -479,6 +304,12 @@ class NonlinearExtension { * transcendental functions. */ TranscendentalSolver d_trSlv; + /** The nonlinear extension object + * + * This is the subsolver responsible for running the procedure for + * constraints involving nonlinear mulitplication. + */ + NlSolver d_nlSlv; /** * The lemmas we computed during collectModelInfo. We store two vectors of * lemmas to be sent out on the output channel of TheoryArith. The first @@ -495,182 +326,6 @@ class NonlinearExtension { std::map> d_approximations; /** have we successfully built the model in this SAT context? */ context::CDO d_builtModel; - - // ordering, stores variables and 0,1,-1 - std::map d_order_vars; - std::vector d_order_points; - - - private: - //per last-call effort check - - //information about monomials - std::vector< Node > d_ms; - std::vector< Node > d_ms_vars; - std::map d_ms_proc; - std::vector d_mterms; - - //list of monomials with factors whose model value is non-constant in model - // e.g. y*cos( x ) - std::map d_m_nconst_factor; - /** the set of monomials we should apply tangent planes to */ - std::unordered_set d_tplane_refine; - // term -> coeff -> rhs -> ( status, exp, b ), - // where we have that : exp => ( coeff * term rhs ) - // b is true if degree( term ) >= degree( rhs ) - std::map > > d_ci; - std::map > > d_ci_exp; - std::map > > d_ci_max; - - // factor skolems - std::map< Node, Node > d_factor_skolem; - Node getFactorSkolem( Node n, std::vector< Node >& lemmas ); - - // tangent plane bounds - std::map< Node, std::map< Node, Node > > d_tangent_val_bound[4]; - - /** 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 - * The argument iter is the number of iterations in the binary search to - * perform. By default, this is set to 15, which is usually enough to be - * precise in the majority of simple cases, whereas not prohibitively - * expensive to compute. - */ - bool getApproximateSqrt(Node c, Node& l, Node& u, unsigned iter = 15) const; - - private: - //-------------------------------------------- lemma schemas - /** check split zero - * - * Returns a set of theory lemmas of the form - * t = 0 V t != 0 - * where t is a term that exists in the current context. - */ - std::vector checkSplitZero(); - - /** check monomial sign - * - * Returns a set of valid theory lemmas, based on a - * lemma schema which ensures that non-linear monomials - * respect sign information based on their facts. - * For more details, see Section 5 of "Design Theory - * Solvers with Extensions" by Reynolds et al., FroCoS 2017, - * Figure 5, this is the schema "Sign". - * - * Examples: - * - * x > 0 ^ y > 0 => x*y > 0 - * x < 0 => x*y*y < 0 - * x = 0 => x*y*z = 0 - */ - std::vector checkMonomialSign(); - - /** check monomial magnitude - * - * Returns a set of valid theory lemmas, based on a - * lemma schema which ensures that comparisons between - * non-linear monomials respect the magnitude of their - * factors. - * For more details, see Section 5 of "Design Theory - * Solvers with Extensions" by Reynolds et al., FroCoS 2017, - * Figure 5, this is the schema "Magnitude". - * - * Examples: - * - * |x|>|y| => |x*z|>|y*z| - * |x|>|y| ^ |z|>|w| ^ |x|>=1 => |x*x*z*u|>|y*w| - * - * Argument c indicates the class of inferences to perform for the (non-linear) - * monomials in the vector d_ms. - * 0 : compare non-linear monomials against 1, - * 1 : compare non-linear monomials against variables, - * 2 : compare non-linear monomials against other non-linear monomials. - */ - std::vector checkMonomialMagnitude( unsigned c ); - - /** check monomial inferred bounds - * - * Returns a set of valid theory lemmas, based on a - * lemma schema that infers new constraints about existing - * terms based on mulitplying both sides of an existing - * constraint by a term. - * For more details, see Section 5 of "Design Theory - * Solvers with Extensions" by Reynolds et al., FroCoS 2017, - * Figure 5, this is the schema "Multiply". - * - * Examples: - * - * x > 0 ^ (y > z + w) => x*y > x*(z+w) - * x < 0 ^ (y > z + w) => x*y < x*(z+w) - * ...where (y > z + w) and x*y are a constraint and term - * that occur in the current context. - */ - std::vector checkMonomialInferBounds( - std::vector& nt_lemmas, - const std::vector& asserts, - const std::vector& false_asserts); - - /** check factoring - * - * Returns a set of valid theory lemmas, based on a - * lemma schema that states a relationship betwen monomials - * with common factors that occur in the same constraint. - * - * Examples: - * - * x*z+y*z > t => ( k = x + y ^ k*z > t ) - * ...where k is fresh and x*z + y*z > t is a - * constraint that occurs in the current context. - */ - std::vector checkFactoring(const std::vector& asserts, - const std::vector& false_asserts); - - /** check monomial infer resolution bounds - * - * Returns a set of valid theory lemmas, based on a - * lemma schema which "resolves" upper bounds - * of one inequality with lower bounds for another. - * This schema is not enabled by default, and can - * be enabled by --nl-ext-rbound. - * - * Examples: - * - * ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t - * ...where s <= x*z and x*y <= t are constraints - * that occur in the current context. - */ - std::vector checkMonomialInferResBounds(); - - /** check tangent planes - * - * Returns a set of valid theory lemmas, based on an - * "incremental linearization" of non-linear monomials. - * This linearization is accomplished by adding constraints - * corresponding to "tangent planes" at the current - * model value of each non-linear monomial. In particular - * consider the definition for constants a,b : - * T_{a,b}( x*y ) = b*x + a*y - a*b. - * The lemmas added by this function are of the form : - * ( ( x>a ^ yb) ) => x*y < T_{a,b}( x*y ) - * ( ( x>a ^ y>b) ^ (x x*y > T_{a,b}( x*y ) - * It is inspired by "Invariant Checking of NRA Transition - * Systems via Incremental Reduction to LRA with EUF" by - * Cimatti et al., TACAS 2017. - * This schema is not terminating in general. - * It is not enabled by default, and can - * be enabled by --nl-ext-tplanes. - * - * Examples: - * - * ( ( x>2 ^ y>5) ^ (x<2 ^ y<5) ) => x*y > 5*x + 2*y - 10 - * ( ( x>2 ^ y<5) ^ (x<2 ^ y>5) ) => x*y < 5*x + 2*y - 10 - */ - std::vector checkTangentPlanes(); - - //-------------------------------------------- end lemma schemas }; /* class NonlinearExtension */ } // namespace arith diff --git a/src/theory/arith/transcendental_solver.h b/src/theory/arith/transcendental_solver.h index 8e539bbf0..88de49af3 100644 --- a/src/theory/arith/transcendental_solver.h +++ b/src/theory/arith/transcendental_solver.h @@ -47,6 +47,15 @@ class TranscendentalSolver ~TranscendentalSolver(); /** init last call + * + * This is called at the beginning of last call effort check, where + * assertions are the set of assertions belonging to arithmetic, + * false_asserts is the subset of assertions that are false in the current + * model, and xts is the set of extended function terms that are active in + * the current context. + * + * This call may add lemmas to lems/lemsPp based on registering term + * information (for example, purification of sine terms). */ void initLastCall(const std::vector& assertions, const std::vector& false_asserts, -- 2.30.2