diff options
| -rw-r--r-- | src/CMakeLists.txt | 2 | ||||
| -rw-r--r-- | src/theory/arith/nl_model.cpp | 1297 | ||||
| -rw-r--r-- | src/theory/arith/nl_model.h | 308 | ||||
| -rw-r--r-- | src/theory/arith/nonlinear_extension.cpp | 1487 | ||||
| -rw-r--r-- | src/theory/arith/nonlinear_extension.h | 167 |
5 files changed, 1805 insertions, 1456 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f2ccbd765..77db3b4db 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -300,6 +300,8 @@ libcvc4_add_sources( theory/arith/linear_equality.h theory/arith/matrix.cpp theory/arith/matrix.h + theory/arith/nl_model.cpp + theory/arith/nl_model.h theory/arith/nonlinear_extension.cpp theory/arith/nonlinear_extension.h theory/arith/normal_form.cpp diff --git a/src/theory/arith/nl_model.cpp b/src/theory/arith/nl_model.cpp new file mode 100644 index 000000000..571fbda6c --- /dev/null +++ b/src/theory/arith/nl_model.cpp @@ -0,0 +1,1297 @@ +/********************* */ +/*! \file nl_model.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 Model object for the non-linear extension class + **/ + +#include "theory/arith/nl_model.h" + +#include "expr/node_algorithm.h" +#include "theory/arith/arith_msum.h" +#include "theory/arith/arith_utilities.h" +#include "theory/rewriter.h" + +using namespace CVC4::kind; + +namespace CVC4 { +namespace theory { +namespace arith { + +NlModel::NlModel(context::Context* c) : d_used_approx(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_two = NodeManager::currentNM()->mkConst(Rational(2)); +} + +NlModel::~NlModel() {} + +void NlModel::reset(TheoryModel* m) +{ + d_model = m; + d_mv[0].clear(); + d_mv[1].clear(); +} + +void NlModel::resetCheck() +{ + d_used_approx = false; + d_check_model_solved.clear(); + d_check_model_bounds.clear(); + d_check_model_vars.clear(); + d_check_model_subs.clear(); +} + +Node NlModel::computeConcreteModelValue(Node n) +{ + return computeModelValue(n, true); +} + +Node NlModel::computeAbstractModelValue(Node n) +{ + return computeModelValue(n, false); +} + +Node NlModel::computeModelValue(Node n, bool isConcrete) +{ + unsigned index = isConcrete ? 0 : 1; + std::map<Node, Node>::iterator it = d_mv[index].find(n); + if (it != d_mv[index].end()) + { + return it->second; + } + Trace("nl-ext-mv-debug") << "computeModelValue " << n << ", index=" << index + << std::endl; + Node ret; + if (n.isConst()) + { + ret = n; + } + else if (index == 1 + && (n.getKind() == NONLINEAR_MULT + || isTranscendentalKind(n.getKind()))) + { + if (hasTerm(n)) + { + // use model value for abstraction + ret = getRepresentative(n); + } + else + { + // abstraction does not exist, use model value + ret = getValueInternal(n); + } + } + else if (n.getNumChildren() == 0) + { + if (n.getKind() == PI) + { + // we are interested in the exact value of PI, which cannot be computed. + // hence, we return PI itself when asked for the concrete value. + ret = n; + } + else + { + ret = getValueInternal(n); + } + } + else + { + // otherwise, compute true value + std::vector<Node> children; + if (n.getMetaKind() == metakind::PARAMETERIZED) + { + children.push_back(n.getOperator()); + } + for (unsigned i = 0; i < n.getNumChildren(); i++) + { + Node mc = computeModelValue(n[i], isConcrete); + children.push_back(mc); + } + ret = NodeManager::currentNM()->mkNode(n.getKind(), children); + if (n.getKind() == APPLY_UF) + { + ret = getValueInternal(ret); + } + else + { + ret = Rewriter::rewrite(ret); + } + } + Trace("nl-ext-mv-debug") << "computed " << (index == 0 ? "M" : "M_A") << "[" + << n << "] = " << ret << std::endl; + d_mv[index][n] = ret; + return ret; +} + +Node NlModel::getValueInternal(Node n) const +{ + return d_model->getValue(n); + /* + std::map< Node, Node >::const_iterator it = d_arithVal.find(n); + if (it!=d_arithVal.end()) + { + return it->second; + } + return Node::null(); + */ +} + +bool NlModel::hasTerm(Node n) const +{ + return d_model->hasTerm(n); + // return d_arithVal.find(n)!=d_arithVal.end(); +} + +Node NlModel::getRepresentative(Node n) const +{ + return d_model->getRepresentative(n); + /* + std::map< Node, Node >::const_iterator it = d_arithVal.find(n); + if (it!=d_arithVal.end()) + { + std::map< Node, Node >::const_iterator itr = d_valToRep.find(it->second); + if (itr != d_valToRep.end()) + { + return itr->second; + } + Assert(false); + } + return Node::null(); + */ +} + +int NlModel::compare(Node i, Node j, bool isConcrete, bool isAbsolute) +{ + Node ci = computeModelValue(i, isConcrete); + Node cj = computeModelValue(j, isConcrete); + if (ci.isConst()) + { + if (cj.isConst()) + { + return compareValue(ci, cj, isAbsolute); + } + return 1; + } + return cj.isConst() ? -1 : 0; +} + +int NlModel::compareValue(Node i, Node j, bool isAbsolute) const +{ + Assert(i.isConst() && j.isConst()); + int ret; + if (i == j) + { + ret = 0; + } + else if (!isAbsolute) + { + ret = i.getConst<Rational>() < j.getConst<Rational>() ? 1 : -1; + } + else + { + ret = (i.getConst<Rational>().abs() == j.getConst<Rational>().abs() + ? 0 + : (i.getConst<Rational>().abs() < j.getConst<Rational>().abs() + ? 1 + : -1)); + } + return ret; +} + +bool NlModel::checkModel(const std::vector<Node>& assertions, + const std::vector<Node>& false_asserts, + unsigned d, + std::vector<Node>& lemmas, + std::vector<Node>& gs) +{ + Trace("nl-ext-cm-debug") << " solve for equalities..." << std::endl; + for (const Node& atom : false_asserts) + { + // see if it corresponds to a univariate polynomial equation of degree two + if (atom.getKind() == EQUAL) + { + if (!solveEqualitySimple(atom, d, lemmas)) + { + // no chance we will satisfy this equality + Trace("nl-ext-cm") << "...check-model : failed to solve equality : " + << atom << std::endl; + } + } + } + + // all remaining variables are constrained to their exact model values + Trace("nl-ext-cm-debug") << " set exact bounds for remaining variables..." + << std::endl; + std::unordered_set<TNode, TNodeHashFunction> visited; + std::vector<TNode> visit; + TNode cur; + for (const Node& a : assertions) + { + visit.push_back(a); + do + { + cur = visit.back(); + visit.pop_back(); + if (visited.find(cur) == visited.end()) + { + visited.insert(cur); + if (cur.getType().isReal() && !cur.isConst()) + { + Kind k = cur.getKind(); + if (k != MULT && k != PLUS && k != NONLINEAR_MULT + && !isTranscendentalKind(k)) + { + // if we have not set an approximate bound for it + if (!hasCheckModelAssignment(cur)) + { + // set its exact model value in the substitution + Node curv = computeConcreteModelValue(cur); + Trace("nl-ext-cm") + << "check-model-bound : exact : " << cur << " = "; + printRationalApprox("nl-ext-cm", curv); + Trace("nl-ext-cm") << std::endl; + bool ret = addCheckModelSubstitution(cur, curv); + AlwaysAssert(ret); + } + } + } + for (const Node& cn : cur) + { + visit.push_back(cn); + } + } + } while (!visit.empty()); + } + + Trace("nl-ext-cm-debug") << " check assertions..." << std::endl; + std::vector<Node> check_assertions; + for (const Node& a : assertions) + { + if (d_check_model_solved.find(a) == d_check_model_solved.end()) + { + Node av = a; + // apply the substitution to a + if (!d_check_model_vars.empty()) + { + av = av.substitute(d_check_model_vars.begin(), + d_check_model_vars.end(), + d_check_model_subs.begin(), + d_check_model_subs.end()); + av = Rewriter::rewrite(av); + } + // simple check literal + if (!simpleCheckModelLit(av)) + { + Trace("nl-ext-cm") << "...check-model : assertion failed : " << a + << std::endl; + check_assertions.push_back(av); + Trace("nl-ext-cm-debug") + << "...check-model : failed assertion, value : " << av << std::endl; + } + } + } + + if (!check_assertions.empty()) + { + Trace("nl-ext-cm") << "...simple check failed." << std::endl; + // TODO (#1450) check model for general case + return false; + } + Trace("nl-ext-cm") << "...simple check succeeded!" << std::endl; + + // must assert and re-check if produce models is true + if (options::produceModels()) + { + NodeManager* nm = NodeManager::currentNM(); + // model guard whose semantics is "the model we constructed holds" + Node mg = nm->mkSkolem("model", nm->booleanType()); + gs.push_back(mg); + // assert the constructed model as assertions + for (const std::pair<const Node, std::pair<Node, Node> > cb : + d_check_model_bounds) + { + Node l = cb.second.first; + Node u = cb.second.second; + Node v = cb.first; + Node pred = nm->mkNode(AND, nm->mkNode(GEQ, v, l), nm->mkNode(GEQ, u, v)); + pred = nm->mkNode(OR, mg.negate(), pred); + lemmas.push_back(pred); + } + } + return true; +} + +bool NlModel::addCheckModelSubstitution(TNode v, TNode s) +{ + // should not substitute the same variable twice + Trace("nl-ext-model") << "* check model substitution : " << v << " -> " << s + << std::endl; + // should not set exact bound more than once + if (std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) + != d_check_model_vars.end()) + { + Trace("nl-ext-model") << "...ERROR: already has value." << std::endl; + // this should never happen since substitutions should be applied eagerly + Assert(false); + return false; + } + // if we previously had an approximate bound, the exact bound should be in its + // range + std::map<Node, std::pair<Node, Node> >::iterator itb = + d_check_model_bounds.find(v); + if (itb != d_check_model_bounds.end()) + { + if (s.getConst<Rational>() >= itb->second.first.getConst<Rational>() + || s.getConst<Rational>() <= itb->second.second.getConst<Rational>()) + { + Trace("nl-ext-model") + << "...ERROR: already has bound which is out of range." << std::endl; + return false; + } + } + for (unsigned i = 0, size = d_check_model_subs.size(); i < size; i++) + { + Node ms = d_check_model_subs[i]; + Node mss = ms.substitute(v, s); + if (mss != ms) + { + mss = Rewriter::rewrite(mss); + } + d_check_model_subs[i] = mss; + } + d_check_model_vars.push_back(v); + d_check_model_subs.push_back(s); + return true; +} + +bool NlModel::addCheckModelBound(TNode v, TNode l, TNode u) +{ + Trace("nl-ext-model") << "* check model bound : " << v << " -> [" << l << " " + << u << "]" << std::endl; + if (l == u) + { + // bound is exact, can add as substitution + return addCheckModelSubstitution(v, l); + } + // should not set a bound for a value that is exact + if (std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) + != d_check_model_vars.end()) + { + Trace("nl-ext-model") + << "...ERROR: setting bound for variable that already has exact value." + << std::endl; + Assert(false); + return false; + } + Assert(l.isConst()); + Assert(u.isConst()); + Assert(l.getConst<Rational>() <= u.getConst<Rational>()); + d_check_model_bounds[v] = std::pair<Node, Node>(l, u); + if (Trace.isOn("nl-ext-cm")) + { + Trace("nl-ext-cm") << "check-model-bound : approximate : "; + printRationalApprox("nl-ext-cm", l); + Trace("nl-ext-cm") << " <= " << v << " <= "; + printRationalApprox("nl-ext-cm", u); + Trace("nl-ext-cm") << std::endl; + } + return true; +} + +bool NlModel::hasCheckModelAssignment(Node v) const +{ + if (d_check_model_bounds.find(v) != d_check_model_bounds.end()) + { + return true; + } + return std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) + != d_check_model_vars.end(); +} + +void NlModel::setUsedApproximate() { d_used_approx = true; } + +bool NlModel::usedApproximate() const { return d_used_approx; } + +bool NlModel::solveEqualitySimple(Node eq, + unsigned d, + std::vector<Node>& lemmas) +{ + Node seq = eq; + if (!d_check_model_vars.empty()) + { + seq = eq.substitute(d_check_model_vars.begin(), + d_check_model_vars.end(), + d_check_model_subs.begin(), + d_check_model_subs.end()); + seq = Rewriter::rewrite(seq); + if (seq.isConst()) + { + if (seq.getConst<bool>()) + { + d_check_model_solved[eq] = Node::null(); + return true; + } + return false; + } + } + Trace("nl-ext-cms") << "simple solve equality " << seq << "..." << std::endl; + Assert(seq.getKind() == EQUAL); + std::map<Node, Node> msum; + if (!ArithMSum::getMonomialSumLit(seq, msum)) + { + Trace("nl-ext-cms") << "...fail, could not determine monomial sum." + << std::endl; + return false; + } + bool is_valid = true; + // the variable we will solve a quadratic equation for + Node var; + Node a = d_zero; + Node b = d_zero; + Node c = d_zero; + NodeManager* nm = NodeManager::currentNM(); + // the list of variables that occur as a monomial in msum, and whose value + // is so far unconstrained in the model. + std::unordered_set<Node, NodeHashFunction> unc_vars; + // the list of variables that occur as a factor in a monomial, and whose + // value is so far unconstrained in the model. + std::unordered_set<Node, NodeHashFunction> unc_vars_factor; + for (std::pair<const Node, Node>& m : msum) + { + Node v = m.first; + Node coeff = m.second.isNull() ? d_one : m.second; + if (v.isNull()) + { + c = coeff; + } + else if (v.getKind() == NONLINEAR_MULT) + { + if (v.getNumChildren() == 2 && v[0].isVar() && v[0] == v[1] + && (var.isNull() || var == v[0])) + { + // may solve quadratic + a = coeff; + var = v[0]; + } + else + { + is_valid = false; + Trace("nl-ext-cms-debug") + << "...invalid due to non-linear monomial " << v << std::endl; + // may wish to set an exact bound for a factor and repeat + for (const Node& vc : v) + { + unc_vars_factor.insert(vc); + } + } + } + else if (!v.isVar() || (!var.isNull() && var != v)) + { + Trace("nl-ext-cms-debug") + << "...invalid due to factor " << v << std::endl; + // cannot solve multivariate + if (is_valid) + { + is_valid = false; + // if b is non-zero, then var is also an unconstrained variable + if (b != d_zero) + { + unc_vars.insert(var); + unc_vars_factor.insert(var); + } + } + // if v is unconstrained, we may turn this equality into a substitution + unc_vars.insert(v); + unc_vars_factor.insert(v); + } + else + { + // set the variable to solve for + b = coeff; + var = v; + } + } + if (!is_valid) + { + // see if we can solve for a variable? + for (const Node& uv : unc_vars) + { + Trace("nl-ext-cm-debug") << "check subs var : " << uv << std::endl; + // cannot already have a bound + if (uv.isVar() && !hasCheckModelAssignment(uv)) + { + Node slv; + Node veqc; + if (ArithMSum::isolate(uv, msum, veqc, slv, EQUAL) != 0) + { + Assert(!slv.isNull()); + // currently do not support substitution-with-coefficients + if (veqc.isNull() && !expr::hasSubterm(slv, uv)) + { + Trace("nl-ext-cm") + << "check-model-subs : " << uv << " -> " << slv << std::endl; + bool ret = addCheckModelSubstitution(uv, slv); + if (ret) + { + Trace("nl-ext-cms") << "...success, model substitution " << uv + << " -> " << slv << std::endl; + d_check_model_solved[eq] = uv; + } + return ret; + } + } + } + } + // see if we can assign a variable to a constant + for (const Node& uvf : unc_vars_factor) + { + Trace("nl-ext-cm-debug") << "check set var : " << uvf << std::endl; + // cannot already have a bound + if (uvf.isVar() && !hasCheckModelAssignment(uvf)) + { + Node uvfv = computeConcreteModelValue(uvf); + Trace("nl-ext-cm") << "check-model-bound : exact : " << uvf << " = "; + printRationalApprox("nl-ext-cm", uvfv); + Trace("nl-ext-cm") << std::endl; + bool ret = addCheckModelSubstitution(uvf, uvfv); + // recurse + return ret ? solveEqualitySimple(eq, d, lemmas) : false; + } + } + Trace("nl-ext-cms") << "...fail due to constrained invalid terms." + << std::endl; + return false; + } + else if (var.isNull() || var.getType().isInteger()) + { + // cannot solve quadratic equations for integer variables + Trace("nl-ext-cms") << "...fail due to variable to solve for." << std::endl; + return false; + } + + // we are linear, it is simple + if (a == d_zero) + { + if (b == d_zero) + { + Trace("nl-ext-cms") << "...fail due to zero a/b." << std::endl; + Assert(false); + return false; + } + Node val = nm->mkConst(-c.getConst<Rational>() / b.getConst<Rational>()); + Trace("nl-ext-cm") << "check-model-bound : exact : " << var << " = "; + printRationalApprox("nl-ext-cm", val); + Trace("nl-ext-cm") << std::endl; + bool ret = addCheckModelSubstitution(var, val); + if (ret) + { + Trace("nl-ext-cms") << "...success, solved linear." << std::endl; + d_check_model_solved[eq] = var; + } + return ret; + } + Trace("nl-ext-quad") << "Solve quadratic : " << seq << std::endl; + Trace("nl-ext-quad") << " a : " << a << std::endl; + Trace("nl-ext-quad") << " b : " << b << std::endl; + Trace("nl-ext-quad") << " c : " << c << std::endl; + Node two_a = nm->mkNode(MULT, d_two, a); + two_a = Rewriter::rewrite(two_a); + Node sqrt_val = nm->mkNode( + MINUS, nm->mkNode(MULT, b, b), nm->mkNode(MULT, d_two, two_a, c)); + sqrt_val = Rewriter::rewrite(sqrt_val); + Trace("nl-ext-quad") << "Will approximate sqrt " << sqrt_val << std::endl; + Assert(sqrt_val.isConst()); + // if it is negative, then we are in conflict + if (sqrt_val.getConst<Rational>().sgn() == -1) + { + Node conf = seq.negate(); + Trace("nl-ext-lemma") << "NlModel::Lemma : quadratic no root : " << conf + << std::endl; + lemmas.push_back(conf); + Trace("nl-ext-cms") << "...fail due to negative discriminant." << std::endl; + return false; + } + if (hasCheckModelAssignment(var)) + { + Trace("nl-ext-cms") << "...fail due to bounds on variable to solve for." + << std::endl; + // two quadratic equations for same variable, give up + return false; + } + // approximate the square root of sqrt_val + Node l, u; + if (!getApproximateSqrt(sqrt_val, l, u, 15 + d)) + { + Trace("nl-ext-cms") << "...fail, could not approximate sqrt." << std::endl; + return false; + } + d_used_approx = true; + Trace("nl-ext-quad") << "...got " << l << " <= sqrt(" << sqrt_val + << ") <= " << u << std::endl; + Node negb = nm->mkConst(-b.getConst<Rational>()); + Node coeffa = nm->mkConst(Rational(1) / two_a.getConst<Rational>()); + // two possible bound regions + Node bounds[2][2]; + Node diff_bound[2]; + Node m_var = computeConcreteModelValue(var); + Assert(m_var.isConst()); + for (unsigned r = 0; r < 2; r++) + { + for (unsigned b = 0; b < 2; b++) + { + Node val = b == 0 ? l : u; + // (-b +- approx_sqrt( b^2 - 4ac ))/2a + Node approx = nm->mkNode( + MULT, coeffa, nm->mkNode(r == 0 ? MINUS : PLUS, negb, val)); + approx = Rewriter::rewrite(approx); + bounds[r][b] = approx; + Assert(approx.isConst()); + } + if (bounds[r][0].getConst<Rational>() > bounds[r][1].getConst<Rational>()) + { + // ensure bound is (lower, upper) + Node tmp = bounds[r][0]; + bounds[r][0] = bounds[r][1]; + bounds[r][1] = tmp; + } + Node diff = + nm->mkNode(MINUS, + m_var, + nm->mkNode(MULT, + nm->mkConst(Rational(1) / Rational(2)), + nm->mkNode(PLUS, bounds[r][0], bounds[r][1]))); + Trace("nl-ext-cm-debug") << "Bound option #" << r << " : "; + printRationalApprox("nl-ext-cm-debug", bounds[r][0]); + Trace("nl-ext-cm-debug") << "..."; + printRationalApprox("nl-ext-cm-debug", bounds[r][1]); + Trace("nl-ext-cm-debug") << std::endl; + diff = Rewriter::rewrite(diff); + Assert(diff.isConst()); + diff = nm->mkConst(diff.getConst<Rational>().abs()); + diff_bound[r] = diff; + Trace("nl-ext-cm-debug") << "...diff from model value ("; + printRationalApprox("nl-ext-cm-debug", m_var); + Trace("nl-ext-cm-debug") << ") is "; + printRationalApprox("nl-ext-cm-debug", diff_bound[r]); + Trace("nl-ext-cm-debug") << std::endl; + } + // take the one that var is closer to in the model + Node cmp = nm->mkNode(GEQ, diff_bound[0], diff_bound[1]); + cmp = Rewriter::rewrite(cmp); + Assert(cmp.isConst()); + unsigned r_use_index = cmp == d_true ? 1 : 0; + Trace("nl-ext-cm") << "check-model-bound : approximate (sqrt) : "; + printRationalApprox("nl-ext-cm", bounds[r_use_index][0]); + Trace("nl-ext-cm") << " <= " << var << " <= "; + printRationalApprox("nl-ext-cm", bounds[r_use_index][1]); + Trace("nl-ext-cm") << std::endl; + bool ret = + addCheckModelBound(var, bounds[r_use_index][0], bounds[r_use_index][1]); + if (ret) + { + d_check_model_solved[eq] = var; + Trace("nl-ext-cms") << "...success, solved quadratic." << std::endl; + } + return ret; +} + +bool NlModel::simpleCheckModelLit(Node lit) +{ + Trace("nl-ext-cms") << "*** Simple check-model lit for " << lit << "..." + << std::endl; + if (lit.isConst()) + { + Trace("nl-ext-cms") << " return constant." << std::endl; + return lit.getConst<bool>(); + } + NodeManager* nm = NodeManager::currentNM(); + bool pol = lit.getKind() != kind::NOT; + Node atom = lit.getKind() == kind::NOT ? lit[0] : lit; + + if (atom.getKind() == EQUAL) + { + // x = a is ( x >= a ^ x <= a ) + for (unsigned i = 0; i < 2; i++) + { + Node lit = nm->mkNode(GEQ, atom[i], atom[1 - i]); + if (!pol) + { + lit = lit.negate(); + } + lit = Rewriter::rewrite(lit); + bool success = simpleCheckModelLit(lit); + if (success != pol) + { + // false != true -> one conjunct of equality is false, we fail + // true != false -> one disjunct of disequality is true, we succeed + return success; + } + } + // both checks passed and polarity is true, or both checks failed and + // polarity is false + return pol; + } + else if (atom.getKind() != GEQ) + { + Trace("nl-ext-cms") << " failed due to unknown literal." << std::endl; + return false; + } + // get the monomial sum + std::map<Node, Node> msum; + if (!ArithMSum::getMonomialSumLit(atom, msum)) + { + Trace("nl-ext-cms") << " failed due to get msum." << std::endl; + return false; + } + // simple interval analysis + if (simpleCheckModelMsum(msum, pol)) + { + return true; + } + // can also try reasoning about univariate quadratic equations + Trace("nl-ext-cms-debug") + << "* Try univariate quadratic analysis..." << std::endl; + std::vector<Node> vs_invalid; + std::unordered_set<Node, NodeHashFunction> vs; + std::map<Node, Node> v_a; + std::map<Node, Node> v_b; + // get coefficients... + for (std::pair<const Node, Node>& m : msum) + { + Node v = m.first; + if (!v.isNull()) + { + if (v.isVar()) + { + v_b[v] = m.second.isNull() ? d_one : m.second; + vs.insert(v); + } + else if (v.getKind() == NONLINEAR_MULT && v.getNumChildren() == 2 + && v[0] == v[1] && v[0].isVar()) + { + v_a[v[0]] = m.second.isNull() ? d_one : m.second; + vs.insert(v[0]); + } + else + { + vs_invalid.push_back(v); + } + } + } + // solve the valid variables... + Node invalid_vsum = vs_invalid.empty() ? d_zero + : (vs_invalid.size() == 1 + ? vs_invalid[0] + : nm->mkNode(PLUS, vs_invalid)); + // substitution to try + std::vector<Node> qvars; + std::vector<Node> qsubs; + for (const Node& v : vs) + { + // is it a valid variable? + std::map<Node, std::pair<Node, Node> >::iterator bit = + d_check_model_bounds.find(v); + if (!expr::hasSubterm(invalid_vsum, v) && bit != d_check_model_bounds.end()) + { + std::map<Node, Node>::iterator it = v_a.find(v); + if (it != v_a.end()) + { + Node a = it->second; + Assert(a.isConst()); + int asgn = a.getConst<Rational>().sgn(); + Assert(asgn != 0); + Node t = nm->mkNode(MULT, a, v, v); + Node b = d_zero; + it = v_b.find(v); + if (it != v_b.end()) + { + b = it->second; + t = nm->mkNode(PLUS, t, nm->mkNode(MULT, b, v)); + } + t = Rewriter::rewrite(t); + Trace("nl-ext-cms-debug") << "Trying to find min/max for quadratic " + << t << "..." << std::endl; + Trace("nl-ext-cms-debug") << " a = " << a << std::endl; + Trace("nl-ext-cms-debug") << " b = " << b << std::endl; + // find maximal/minimal value on the interval + Node apex = nm->mkNode( + DIVISION, nm->mkNode(UMINUS, b), nm->mkNode(MULT, d_two, a)); + apex = Rewriter::rewrite(apex); + Assert(apex.isConst()); + // for lower, upper, whether we are greater than the apex + bool cmp[2]; + Node boundn[2]; + for (unsigned r = 0; r < 2; r++) + { + boundn[r] = r == 0 ? bit->second.first : bit->second.second; + Node cmpn = nm->mkNode(GT, boundn[r], apex); + cmpn = Rewriter::rewrite(cmpn); + Assert(cmpn.isConst()); + cmp[r] = cmpn.getConst<bool>(); + } + Trace("nl-ext-cms-debug") << " apex " << apex << std::endl; + Trace("nl-ext-cms-debug") + << " lower " << boundn[0] << ", cmp: " << cmp[0] << std::endl; + Trace("nl-ext-cms-debug") + << " upper " << boundn[1] << ", cmp: " << cmp[1] << std::endl; + Assert(boundn[0].getConst<Rational>() + <= boundn[1].getConst<Rational>()); + Node s; + qvars.push_back(v); + if (cmp[0] != cmp[1]) + { + Assert(!cmp[0] && cmp[1]); + // does the sign match the bound? + if ((asgn == 1) == pol) + { + // the apex is the max/min value + s = apex; + Trace("nl-ext-cms-debug") << " ...set to apex." << std::endl; + } + else + { + // it is one of the endpoints, plug in and compare + Node tcmpn[2]; + for (unsigned r = 0; r < 2; r++) + { + qsubs.push_back(boundn[r]); + Node ts = t.substitute( + qvars.begin(), qvars.end(), qsubs.begin(), qsubs.end()); + tcmpn[r] = Rewriter::rewrite(ts); + qsubs.pop_back(); + } + Node tcmp = nm->mkNode(LT, tcmpn[0], tcmpn[1]); + Trace("nl-ext-cms-debug") + << " ...both sides of apex, compare " << tcmp << std::endl; + tcmp = Rewriter::rewrite(tcmp); + Assert(tcmp.isConst()); + unsigned bindex_use = (tcmp.getConst<bool>() == pol) ? 1 : 0; + Trace("nl-ext-cms-debug") + << " ...set to " << (bindex_use == 1 ? "upper" : "lower") + << std::endl; + s = boundn[bindex_use]; + } + } + else + { + // both to one side of the apex + // we figure out which bound to use (lower or upper) based on + // three factors: + // (1) whether a's sign is positive, + // (2) whether we are greater than the apex of the parabola, + // (3) the polarity of the constraint, i.e. >= or <=. + // there are 8 cases of these factors, which we test here. + unsigned bindex_use = (((asgn == 1) == cmp[0]) == pol) ? 0 : 1; + Trace("nl-ext-cms-debug") + << " ...set to " << (bindex_use == 1 ? "upper" : "lower") + << std::endl; + s = boundn[bindex_use]; + } + Assert(!s.isNull()); + qsubs.push_back(s); + Trace("nl-ext-cms") << "* set bound based on quadratic : " << v + << " -> " << s << std::endl; + } + } + } + if (!qvars.empty()) + { + Assert(qvars.size() == qsubs.size()); + Node slit = + lit.substitute(qvars.begin(), qvars.end(), qsubs.begin(), qsubs.end()); + slit = Rewriter::rewrite(slit); + return simpleCheckModelLit(slit); + } + return false; +} + +bool NlModel::simpleCheckModelMsum(const std::map<Node, Node>& msum, bool pol) +{ + Trace("nl-ext-cms-debug") << "* Try simple interval analysis..." << std::endl; + NodeManager* nm = NodeManager::currentNM(); + // map from transcendental functions to whether they were set to lower + // bound + bool simpleSuccess = true; + std::map<Node, bool> set_bound; + std::vector<Node> sum_bound; + for (const std::pair<const Node, Node>& m : msum) + { + Node v = m.first; + if (v.isNull()) + { + sum_bound.push_back(m.second.isNull() ? d_one : m.second); + } + else + { + Trace("nl-ext-cms-debug") << "- monomial : " << v << std::endl; + // --- whether we should set a lower bound for this monomial + bool set_lower = + (m.second.isNull() || m.second.getConst<Rational>().sgn() == 1) + == pol; + Trace("nl-ext-cms-debug") + << "set bound to " << (set_lower ? "lower" : "upper") << std::endl; + + // --- Collect variables and factors in v + std::vector<Node> vars; + std::vector<unsigned> factors; + if (v.getKind() == NONLINEAR_MULT) + { + unsigned last_start = 0; + for (unsigned i = 0, nchildren = v.getNumChildren(); i < nchildren; i++) + { + // are we at the end? + if (i + 1 == nchildren || v[i + 1] != v[i]) + { + unsigned vfact = 1 + (i - last_start); + last_start = (i + 1); + vars.push_back(v[i]); + factors.push_back(vfact); + } + } + } + else + { + vars.push_back(v); + factors.push_back(1); + } + + // --- Get the lower and upper bounds and sign information. + // Whether we have an (odd) number of negative factors in vars, apart + // from the variable at choose_index. + bool has_neg_factor = false; + int choose_index = -1; + std::vector<Node> ls; + std::vector<Node> us; + // the relevant sign information for variables with odd exponents: + // 1: both signs of the interval of this variable are positive, + // -1: both signs of the interval of this variable are negative. + std::vector<int> signs; + Trace("nl-ext-cms-debug") << "get sign information..." << std::endl; + for (unsigned i = 0, size = vars.size(); i < size; i++) + { + Node vc = vars[i]; + unsigned vcfact = factors[i]; + if (Trace.isOn("nl-ext-cms-debug")) + { + Trace("nl-ext-cms-debug") << "-- " << vc; + if (vcfact > 1) + { + Trace("nl-ext-cms-debug") << "^" << vcfact; + } + Trace("nl-ext-cms-debug") << " "; + } + std::map<Node, std::pair<Node, Node> >::iterator bit = + d_check_model_bounds.find(vc); + // if there is a model bound for this term + if (bit != d_check_model_bounds.end()) + { + Node l = bit->second.first; + Node u = bit->second.second; + ls.push_back(l); + us.push_back(u); + int vsign = 0; + if (vcfact % 2 == 1) + { + vsign = 1; + int lsgn = l.getConst<Rational>().sgn(); + int usgn = u.getConst<Rational>().sgn(); + Trace("nl-ext-cms-debug") + << "bound_sign(" << lsgn << "," << usgn << ") "; + if (lsgn == -1) + { + if (usgn < 1) + { + // must have a negative factor + has_neg_factor = !has_neg_factor; + vsign = -1; + } + else if (choose_index == -1) + { + // set the choose index to this + choose_index = i; + vsign = 0; + } + else + { + // ambiguous, can't determine the bound + Trace("nl-ext-cms") + << " failed due to ambiguious monomial." << std::endl; + return false; + } + } + } + Trace("nl-ext-cms-debug") << " -> " << vsign << std::endl; + signs.push_back(vsign); + } + else + { + Trace("nl-ext-cms-debug") << std::endl; + Trace("nl-ext-cms") + << " failed due to unknown bound for " << vc << std::endl; + // should either assign a model bound or eliminate the variable + // via substitution + Assert(false); + return false; + } + } + // whether we will try to minimize/maximize (-1/1) the absolute value + int setAbs = (set_lower == has_neg_factor) ? 1 : -1; + Trace("nl-ext-cms-debug") + << "set absolute value to " << (setAbs == 1 ? "maximal" : "minimal") + << std::endl; + + std::vector<Node> vbs; + Trace("nl-ext-cms-debug") << "set bounds..." << std::endl; + for (unsigned i = 0, size = vars.size(); i < size; i++) + { + Node vc = vars[i]; + unsigned vcfact = factors[i]; + Node l = ls[i]; + Node u = us[i]; + bool vc_set_lower; + int vcsign = signs[i]; + Trace("nl-ext-cms-debug") + << "Bounds for " << vc << " : " << l << ", " << u + << ", sign : " << vcsign << ", factor : " << vcfact << std::endl; + if (l == u) + { + // by convention, always say it is lower if they are the same + vc_set_lower = true; + Trace("nl-ext-cms-debug") + << "..." << vc << " equal bound, set to lower" << std::endl; + } + else + { + if (vcfact % 2 == 0) + { + // minimize or maximize its absolute value + Rational la = l.getConst<Rational>().abs(); + Rational ua = u.getConst<Rational>().abs(); + if (la == ua) + { + // by convention, always say it is lower if abs are the same + vc_set_lower = true; + Trace("nl-ext-cms-debug") + << "..." << vc << " equal abs, set to lower" << std::endl; + } + else + { + vc_set_lower = (la > ua) == (setAbs == 1); + } + } + else if (signs[i] == 0) + { + // we choose this index to match the overall set_lower + vc_set_lower = set_lower; + } + else + { + vc_set_lower = (signs[i] != setAbs); + } + Trace("nl-ext-cms-debug") + << "..." << vc << " set to " << (vc_set_lower ? "lower" : "upper") + << std::endl; + } + // check whether this is a conflicting bound + std::map<Node, bool>::iterator itsb = set_bound.find(vc); + if (itsb == set_bound.end()) + { + set_bound[vc] = vc_set_lower; + } + else if (itsb->second != vc_set_lower) + { + Trace("nl-ext-cms") + << " failed due to conflicting bound for " << vc << std::endl; + return false; + } + // must over/under approximate based on vc_set_lower, computed above + Node vb = vc_set_lower ? l : u; + for (unsigned i = 0; i < vcfact; i++) + { + vbs.push_back(vb); + } + } + if (!simpleSuccess) + { + break; + } + Node vbound = vbs.size() == 1 ? vbs[0] : nm->mkNode(MULT, vbs); + sum_bound.push_back(ArithMSum::mkCoeffTerm(m.second, vbound)); + } + } + // if the exact bound was computed via simple analysis above + // make the bound + Node bound; + if (sum_bound.size() > 1) + { + bound = nm->mkNode(kind::PLUS, sum_bound); + } + else if (sum_bound.size() == 1) + { + bound = sum_bound[0]; + } + else + { + bound = d_zero; + } + // make the comparison + Node comp = nm->mkNode(kind::GEQ, bound, d_zero); + if (!pol) + { + comp = comp.negate(); + } + Trace("nl-ext-cms") << " comparison is : " << comp << std::endl; + comp = Rewriter::rewrite(comp); + Assert(comp.isConst()); + Trace("nl-ext-cms") << " returned : " << comp << std::endl; + return comp == d_true; +} + +bool NlModel::isRefineableTfFun(Node tf) +{ + Assert(tf.getKind() == SINE || tf.getKind() == EXPONENTIAL); + if (tf.getKind() == SINE) + { + // we do not consider e.g. sin( -1*x ), since considering sin( x ) will + // have the same effect. We also do not consider sin(x+y) since this is + // handled by introducing a fresh variable (see the map d_tr_base in + // NonlinearExtension). + if (!tf[0].isVar()) + { + return false; + } + } + // Figure 3 : c + Node c = computeAbstractModelValue(tf[0]); + Assert(c.isConst()); + int csign = c.getConst<Rational>().sgn(); + if (csign == 0) + { + return false; + } + return true; +} + +bool NlModel::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>(); + + 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; +} + +void NlModel::recordApproximations() +{ + // Record the approximations we used. This code calls the + // recordApproximation method of the model, which overrides the model + // values for variables that we solved for, using techniques specific to + // this class. + NodeManager* nm = NodeManager::currentNM(); + for (const std::pair<const Node, std::pair<Node, Node> >& cb : + d_check_model_bounds) + { + Node l = cb.second.first; + Node u = cb.second.second; + Node pred; + Node v = cb.first; + if (l != u) + { + pred = nm->mkNode(AND, nm->mkNode(GEQ, v, l), nm->mkNode(GEQ, u, v)); + } + else if (!d_model->areEqual(v, l)) + { + // only record if value was not equal already + pred = v.eqNode(l); + } + if (!pred.isNull()) + { + pred = Rewriter::rewrite(pred); + d_model->recordApproximation(v, pred); + } + } + // Also record the exact values we used. An exact value can be seen as a + // special kind approximation of the form (choice x. x = exact_value). + // Notice that the above term gets rewritten such that the choice function + // is eliminated. + for (size_t i = 0, num = d_check_model_vars.size(); i < num; i++) + { + Node v = d_check_model_vars[i]; + Node s = d_check_model_subs[i]; + if (!d_model->areEqual(v, s)) + { + Node pred = v.eqNode(s); + pred = Rewriter::rewrite(pred); + d_model->recordApproximation(v, pred); + } + } +} +void NlModel::printModelValue(const char* c, Node n, unsigned prec) const +{ + if (Trace.isOn(c)) + { + Trace(c) << " " << n << " -> "; + for (unsigned i = 0; i < 2; i++) + { + std::map<Node, Node>::const_iterator it = d_mv[1 - i].find(n); + Assert(it != d_mv[1 - i].end()); + if (it->second.isConst()) + { + printRationalApprox(c, it->second, prec); + } + else + { + Trace(c) << "?"; // it->second; + } + Trace(c) << (i == 0 ? " [actual: " : " ]"); + } + Trace(c) << std::endl; + } +} +} // namespace arith +} // namespace theory +} // namespace CVC4 diff --git a/src/theory/arith/nl_model.h b/src/theory/arith/nl_model.h new file mode 100644 index 000000000..19341cc5f --- /dev/null +++ b/src/theory/arith/nl_model.h @@ -0,0 +1,308 @@ +/********************* */ +/*! \file nl_model.h + ** \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 Model object for the non-linear extension class + **/ + +#ifndef CVC4__THEORY__ARITH__NL_MODEL_H +#define CVC4__THEORY__ARITH__NL_MODEL_H + +#include <map> +#include <unordered_map> +#include <vector> + +#include "context/cdo.h" +#include "context/context.h" +#include "expr/kind.h" +#include "expr/node.h" +#include "theory/theory_model.h" + +namespace CVC4 { +namespace theory { +namespace arith { + +class NonlinearExtension; + +/** Non-linear model finder + * + * This class is responsible for all queries related to the (candidate) model + * that is being processed by the non-linear arithmetic solver. It further + * implements techniques for finding modifications to the current candidate + * model in the case it can determine that a model exists. These include + * techniques based on solving (quadratic) equations and bound analysis. + */ +class NlModel +{ + friend class NonlinearExtension; + + public: + NlModel(context::Context* c); + ~NlModel(); + /** reset + * + * This method is called once at the beginning of a last call effort check, + * where m is the model of the theory of arithmetic. This method resets the + * cache of computed model values. + */ + void reset(TheoryModel* m); + /** reset check + * + * This method is called when the non-linear arithmetic solver restarts + * its computation of lemmas and models during a last call effort check. + */ + void resetCheck(); + /** compute model value + * + * This computes model values for terms based on two semantics, a "concrete" + * semantics and an "abstract" semantics. + * + * if isConcrete is true, this means compute the value of n based on its + * children recursively. (we call this its "concrete" value) + * if isConcrete is false, this means lookup the value of n in the model. + * (we call this its "abstract" value) + * In other words, !isConcrete treats multiplication terms and transcendental + * function applications as variables, whereas isConcrete computes their + * actual values based on the semantics of multiplication. This is a key + * distinction used in the model-based refinement scheme in Cimatti et al. + * TACAS 2017. + * + * For example, if M( a ) = 2, M( b ) = 3, M( a*b ) = 5, i.e. the variable + * for a*b has been assigned a value 5 by the linear solver, then : + * + * computeModelValue( a*b, true ) = + * computeModelValue( a, true )*computeModelValue( b, true ) = 2*3 = 6 + * whereas: + * computeModelValue( a*b, false ) = 5 + */ + Node computeConcreteModelValue(Node n); + Node computeAbstractModelValue(Node n); + Node computeModelValue(Node n, bool isConcrete); + + /** Compare arithmetic terms i and j based an ordering. + * + * This returns: + * -1 if i < j, 1 if i > j, or 0 if i == j + * + * If isConcrete is true, we consider the concrete model values of i and j, + * otherwise, we consider their abstract model values. For definitions of + * concrete vs abstract model values, see NlModel::computeModelValue. + * + * If isAbsolute is true, we compare the absolute value of thee above + * values. + */ + int compare(Node i, Node j, bool isConcrete, bool isAbsolute); + /** Compare arithmetic terms i and j based an ordering. + * + * This returns: + * -1 if i < j, 1 if i > j, or 0 if i == j + * + * If isAbsolute is true, we compare the absolute value of i and j + */ + int compareValue(Node i, Node j, bool isAbsolute) const; + + //------------------------------ recording model substitutions and bounds + /** add check model substitution + * + * Adds the model substitution v -> s. This applies the substitution + * { v -> s } to each term in d_check_model_subs and adds v,s to + * d_check_model_vars and d_check_model_subs respectively. + * If this method returns false, then the substitution v -> s is inconsistent + * with the current substitution and bounds. + */ + bool addCheckModelSubstitution(TNode v, TNode s); + /** add check model bound + * + * Adds the bound x -> < l, u > to the map above, and records the + * approximation ( x, l <= x <= u ) in the model. This method returns false + * if the bound is inconsistent with the current model substitution or + * bounds. + */ + bool addCheckModelBound(TNode v, TNode l, TNode u); + /** has check model assignment + * + * Have we assigned v in the current checkModel(...) call? + * + * This method returns true if variable v is in the domain of + * d_check_model_bounds or if it occurs in d_check_model_vars. + */ + bool hasCheckModelAssignment(Node v) const; + /** Check model + * + * Checks the current model based on solving for equalities, and using error + * bounds on the Taylor approximation. + * + * If this function returns true, then all assertions in the input argument + * "assertions" are satisfied for all interpretations of variables within + * their computed bounds (as stored in d_check_model_bounds). + * + * For details, see Section 3 of Cimatti et al CADE 2017 under the heading + * "Detecting Satisfiable Formulas". + * + * d is a degree indicating how precise our computations are. + */ + bool checkModel(const std::vector<Node>& assertions, + const std::vector<Node>& false_asserts, + unsigned d, + std::vector<Node>& lemmas, + std::vector<Node>& gs); + /** + * Set that we have used an approximation during this check. This flag is + * reset on a call to resetCheck. It is set when we use reasoning that + * is limited by a degree of precision we are using. In other words, if we + * used an approximation, then we maybe could still establish a lemma or + * determine the input is SAT if we increased our precision. + */ + void setUsedApproximate(); + /** Did we use an approximation during this check? */ + bool usedApproximate() const; + //------------------------------ end recording model substitutions and bounds + + /** print model value, for debugging. + * + * This prints both the abstract and concrete model values for arithmetic + * term n on Trace c with precision prec. + */ + void printModelValue(const char* c, Node n, unsigned prec = 5) const; + /** record approximations in the current model + * + * This makes necessary calls that notify the model of any approximations + * that were used by this solver. + */ + void recordApproximations(); + + private: + /** The current model */ + TheoryModel* d_model; + /** Get the model value of n from the model object above */ + Node getValueInternal(Node n) const; + /** Does the equality engine of the model have term n? */ + bool hasTerm(Node n) const; + /** Get the representative of n in the model */ + Node getRepresentative(Node n) const; + + //---------------------------check model + /** solve equality simple + * + * This method is used during checkModel(...). It takes as input an + * equality eq. If it returns true, then eq is correct-by-construction based + * on the information stored in our model representation (see + * d_check_model_vars, d_check_model_subs, d_check_model_bounds), and eq + * is added to d_check_model_solved. The equality eq may involve any + * number of variables, and monomials of arbitrary degree. If this method + * returns false, then we did not show that the equality was true in the + * model. This method uses incomplete techniques based on interval + * analysis and quadratic equation solving. + * + * If it can be shown that the equality must be false in the current + * model, then we may add a lemma to lemmas explaining why this is the case. + * For instance, if eq reduces to a univariate quadratic equation with no + * root, we send a conflict clause of the form a*x^2 + b*x + c != 0. + */ + bool solveEqualitySimple(Node eq, unsigned d, std::vector<Node>& lemmas); + + /** simple check model for transcendental functions for literal + * + * This method returns true if literal is true for all interpretations of + * transcendental functions within their error bounds (as stored + * in d_check_model_bounds). This is determined by a simple under/over + * approximation of the value of sum of (linear) monomials. For example, + * if we determine that .8 < sin( 1 ) < .9, this function will return + * true for literals like: + * 2.0*sin( 1 ) > 1.5 + * -1.0*sin( 1 ) < -0.79 + * -1.0*sin( 1 ) > -0.91 + * sin( 1 )*sin( 1 ) + sin( 1 ) > 0.0 + * It will return false for literals like: + * sin( 1 ) > 0.85 + * It will also return false for literals like: + * -0.3*sin( 1 )*sin( 2 ) + sin( 2 ) > .7 + * sin( sin( 1 ) ) > .5 + * since the bounds on these terms cannot quickly be determined. + */ + bool simpleCheckModelLit(Node lit); + bool simpleCheckModelMsum(const std::map<Node, Node>& msum, bool pol); + //---------------------------end check model + + /** is refinable transcendental function + * + * A transcendental function application is not refineable if its current + * model value is zero, or if it is an application of SINE applied + * to a non-variable. + */ + bool isRefineableTfFun(Node tf); + + /** 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; + + /** commonly used terms */ + Node d_zero; + Node d_one; + Node d_two; + Node d_true; + Node d_false; + Node d_null; + /** cache of model values + * + * Stores the the concrete/abstract model values. This is a cache of the + * computeModelValue method. + */ + std::map<Node, Node> d_mv[2]; + /** + * A substitution from variables that appear in assertions to a solved form + * term. These vectors are ordered in the form: + * x_1 -> t_1 ... x_n -> t_n + * where x_i is not in the free variables of t_j for j>=i. + */ + std::vector<Node> d_check_model_vars; + std::vector<Node> d_check_model_subs; + /** lower and upper bounds for check model + * + * For each term t in the domain of this map, if this stores the pair + * (c_l, c_u) then the model M is such that c_l <= M( t ) <= c_u. + * + * We add terms whose value is approximated in the model to this map, which + * includes: + * (1) applications of transcendental functions, whose value is approximated + * by the Taylor series, + * (2) variables we have solved quadratic equations for, whose value + * involves approximations of square roots. + */ + std::map<Node, std::pair<Node, Node> > d_check_model_bounds; + /** + * The map from literals that our model construction solved, to the variable + * that was solved for. Examples of such literals are: + * (1) Equalities x = t, which we turned into a model substitution x -> t, + * where x not in FV( t ), and + * (2) Equalities a*x*x + b*x + c = 0, which we turned into a model bound + * -b+s*sqrt(b*b-4*a*c)/2a - E <= x <= -b+s*sqrt(b*b-4*a*c)/2a + E. + * + * These literals are exempt from check-model, since they are satisfied by + * definition of our model construction. + */ + std::unordered_map<Node, Node, NodeHashFunction> d_check_model_solved; + /** did we use an approximation on this call to last-call effort? */ + bool d_used_approx; +}; /* class NlModel */ + +} // namespace arith +} // namespace theory +} // namespace CVC4 + +#endif /* CVC4__THEORY__ARITH__NONLINEAR_EXTENSION_H */ diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp index 839520c0b..e14f07a4f 100644 --- a/src/theory/arith/nonlinear_extension.cpp +++ b/src/theory/arith/nonlinear_extension.cpp @@ -113,19 +113,46 @@ void debugPrintBound(const char* c, Node coeff, Node x, Kind type, Node rhs) { Trace(c) << t << " " << type << " " << rhs; } -struct SortNonlinearExtension { - SortNonlinearExtension() - : d_nla(NULL), d_order_type(0), d_reverse_order(false) {} - NonlinearExtension* d_nla; - unsigned d_order_type; +struct SortNlModel +{ + SortNlModel() + : d_nlm(nullptr), + d_isConcrete(true), + d_isAbsolute(false), + d_reverse_order(false) + { + } + /** pointer to the model */ + NlModel* d_nlm; + /** are we comparing concrete model values? */ + bool d_isConcrete; + /** are we comparing absolute values? */ + bool d_isAbsolute; + /** are we in reverse order? */ bool d_reverse_order; + /** the comparison */ bool operator()(Node i, Node j) { - int cv = d_nla->compare(i, j, d_order_type); + int cv = d_nlm->compare(i, j, d_isConcrete, d_isAbsolute); if (cv == 0) { return i < j; - } else { - return d_reverse_order ? cv < 0 : cv > 0; } + return d_reverse_order ? cv < 0 : cv > 0; + } +}; +struct SortNonlinearDegree +{ + SortNonlinearDegree(NodeMultiset& m) : d_mdegree(m) {} + /** pointer to the non-linear extension */ + NodeMultiset& d_mdegree; + /** + * Sorts by degree of the monomials, where lower degree monomials come + * first. + */ + bool operator()(Node i, Node j) + { + unsigned i_count = getCount(d_mdegree, i); + unsigned j_count = getCount(d_mdegree, j); + return i_count == j_count ? (i < j) : (i_count < j_count ? true : false); } }; @@ -156,13 +183,14 @@ bool hasNewMonomials(Node n, const std::vector<Node>& existing) { NonlinearExtension::NonlinearExtension(TheoryArith& containing, eq::EqualityEngine* ee) - : d_builtModel(containing.getSatContext(), false), - d_lemmas(containing.getUserContext()), + : d_lemmas(containing.getUserContext()), d_zero_split(containing.getUserContext()), d_skolem_atoms(containing.getUserContext()), d_containing(containing), d_ee(ee), - d_needsLastCall(false) + d_needsLastCall(false), + d_model(containing.getSatContext()), + d_builtModel(containing.getSatContext(), false) { d_true = NodeManager::currentNM()->mkConst(true); d_false = NodeManager::currentNM()->mkConst(false); @@ -180,7 +208,6 @@ NonlinearExtension::NonlinearExtension(TheoryArith& containing, d_taylor_real_fv_base_rem = NodeManager::currentNM()->mkBoundVar( "b", NodeManager::currentNM()->realType()); d_taylor_degree = options::nlExtTfTaylorDegree(); - d_used_approx = false; } NonlinearExtension::~NonlinearExtension() {} @@ -318,76 +345,6 @@ std::pair<bool, Node> NonlinearExtension::isExtfReduced( return std::make_pair(true, Node::null()); } -Node NonlinearExtension::computeModelValue(Node n, unsigned index) { - std::map<Node, Node>::iterator it = d_mv[index].find(n); - if (it != d_mv[index].end()) { - return it->second; - } else { - Trace("nl-ext-mv-debug") << "computeModelValue " << n << ", index=" << index << std::endl; - Node ret; - if (n.isConst()) { - ret = n; - } - else if (index == 1 && (n.getKind() == NONLINEAR_MULT - || isTranscendentalKind(n.getKind()))) - { - if (d_containing.getValuation().getModel()->hasTerm(n)) { - // use model value for abstraction - ret = d_containing.getValuation().getModel()->getRepresentative(n); - } else { - // abstraction does not exist, use model value - //ret = computeModelValue(n, 0); - ret = d_containing.getValuation().getModel()->getValue(n); - } - //Assert( ret.isConst() ); - } else if (n.getNumChildren() == 0) { - if (n.getKind() == PI) - { - // we are interested in the exact value of PI, which cannot be computed. - // hence, we return PI itself when asked for the concrete value. - ret = n; - }else{ - ret = d_containing.getValuation().getModel()->getValue(n); - } - } else { - // otherwise, compute true value - std::vector<Node> children; - if (n.getMetaKind() == metakind::PARAMETERIZED) - { - children.push_back(n.getOperator()); - } - for (unsigned i = 0; i < n.getNumChildren(); i++) { - Node mc = computeModelValue(n[i], index); - children.push_back(mc); - } - ret = NodeManager::currentNM()->mkNode(n.getKind(), children); - if (n.getKind() == APPLY_UF) - { - ret = d_containing.getValuation().getModel()->getValue(ret); - }else{ - ret = Rewriter::rewrite(ret); - } - /* - if (!ret.isConst()) { - Trace("nl-ext-debug") << "...got non-constant : " << ret << " for " - << n << ", ask model directly." << std::endl; - ret = d_containing.getValuation().getModel()->getValue(ret); - } - */ - } - //if (ret.getType().isReal() && !isArithKind(n.getKind())) { - // Trace("nl-ext-mv-debug") << ( index==0 ? "M" : "M_A" ) << "[ " << n - // << " ] -> " << ret << std::endl; - //may involve transcendental functions - //Assert(ret.isConst()); - //} - Trace("nl-ext-mv-debug") << "computed " << (index == 0 ? "M" : "M_A") << "[" - << n << "] = " << ret << std::endl; - d_mv[index][n] = ret; - return ret; - } -} - void NonlinearExtension::registerMonomial(Node n) { if (!IsInVector(d_monomials, n)) { d_monomials.push_back(n); @@ -497,12 +454,16 @@ bool NonlinearExtension::isArithKind(Kind k) { return k == PLUS || k == MULT || k == NONLINEAR_MULT; } -Node NonlinearExtension::mkLit(Node a, Node b, int status, int orderType) { +Node NonlinearExtension::mkLit(Node a, Node b, int status, bool isAbsolute) +{ if (status == 0) { Node a_eq_b = a.eqNode(b); - if (orderType == 0) { + if (!isAbsolute) + { return a_eq_b; - } else { + } + else + { // return mkAbs( a ).eqNode( mkAbs( b ) ); Node negate_b = NodeManager::currentNM()->mkNode(UMINUS, b); return a_eq_b.orNode(a.eqNode(negate_b)); @@ -513,9 +474,12 @@ Node NonlinearExtension::mkLit(Node a, Node b, int status, int orderType) { Assert(status == 1 || status == 2); NodeManager* nm = NodeManager::currentNM(); Kind greater_op = status == 1 ? GEQ : GT; - if (orderType == 0) { + if (!isAbsolute) + { return nm->mkNode(greater_op, a, b); - } else { + } + else + { // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) ); Node zero = mkRationalNode(0); Node a_is_nonnegative = nm->mkNode(GEQ, a, zero); @@ -751,7 +715,7 @@ std::vector<Node> NonlinearExtension::checkModelEval( Node lit = assertions[i]; Node atom = lit.getKind()==NOT ? lit[0] : lit; if( d_skolem_atoms.find( atom )==d_skolem_atoms.end() ){ - Node litv = computeModelValue(lit); + Node litv = d_model.computeConcreteModelValue(lit); Trace("nl-ext-mv-assert") << "M[[ " << lit << " ]] -> " << litv; if (litv != d_true) { Trace("nl-ext-mv-assert") << " [model-false]" << std::endl; @@ -769,10 +733,6 @@ bool NonlinearExtension::checkModel(const std::vector<Node>& assertions, const std::vector<Node>& false_asserts) { Trace("nl-ext-cm") << "--- check-model ---" << std::endl; - d_check_model_solved.clear(); - d_check_model_bounds.clear(); - d_check_model_vars.clear(); - d_check_model_subs.clear(); // get the presubstitution Trace("nl-ext-cm-debug") << " apply pre-substitution..." << std::endl; @@ -811,16 +771,16 @@ bool NonlinearExtension::checkModel(const std::vector<Node>& assertions, { bool success = true; // tf is Figure 3 : tf( x ) - Node atf = computeModelValue(tf, 0); + Node atf = d_model.computeConcreteModelValue(tf); if (k == PI) { - success = addCheckModelBound(atf, d_pi_bound[0], d_pi_bound[1]); + success = d_model.addCheckModelBound(atf, d_pi_bound[0], d_pi_bound[1]); } - else if (isRefineableTfFun(tf)) + else if (d_model.isRefineableTfFun(tf)) { - d_used_approx = true; + d_model.setUsedApproximate(); std::pair<Node, Node> bounds = getTfModelBounds(tf, d_taylor_degree); - success = addCheckModelBound(atf, bounds.first, bounds.second); + success = d_model.addCheckModelBound(atf, bounds.first, bounds.second); } if (!success) { @@ -829,952 +789,28 @@ bool NonlinearExtension::checkModel(const std::vector<Node>& assertions, << std::endl; return false; } - if (Trace.isOn("nl-ext-cm")) - { - std::map<Node, std::pair<Node, Node> >::iterator it = - d_check_model_bounds.find(tf); - if (it != d_check_model_bounds.end()) - { - Trace("nl-ext-cm") << "check-model-bound : approximate (taylor) : "; - printRationalApprox("nl-ext-cm", it->second.first); - Trace("nl-ext-cm") << " <= " << tf << " <= "; - printRationalApprox("nl-ext-cm", it->second.second); - Trace("nl-ext-cm") << std::endl; - } - } - } - } - - Trace("nl-ext-cm-debug") << " solve for equalities..." << std::endl; - for (const Node& atom : false_asserts) - { - // see if it corresponds to a univariate polynomial equation of degree two - if (atom.getKind() == EQUAL) - { - if (!solveEqualitySimple(atom)) - { - // no chance we will satisfy this equality - Trace("nl-ext-cm") << "...check-model : failed to solve equality : " - << atom << std::endl; - } - } - } - - // all remaining variables are constrained to their exact model values - Trace("nl-ext-cm-debug") << " set exact bounds for remaining variables..." - << std::endl; - std::unordered_set<TNode, TNodeHashFunction> visited; - std::vector<TNode> visit; - TNode cur; - for (const Node& a : passertions) - { - visit.push_back(a); - do - { - cur = visit.back(); - visit.pop_back(); - if (visited.find(cur) == visited.end()) - { - visited.insert(cur); - if (cur.getType().isReal() && !cur.isConst()) - { - Kind k = cur.getKind(); - if (k != MULT && k != PLUS && k != NONLINEAR_MULT - && !isTranscendentalKind(k)) - { - // if we have not set an approximate bound for it - if (!hasCheckModelAssignment(cur)) - { - // set its exact model value in the substitution - Node curv = computeModelValue(cur); - Trace("nl-ext-cm") - << "check-model-bound : exact : " << cur << " = "; - printRationalApprox("nl-ext-cm", curv); - Trace("nl-ext-cm") << std::endl; - bool ret = addCheckModelSubstitution(cur, curv); - AlwaysAssert(ret); - } - } - } - for (const Node& cn : cur) - { - visit.push_back(cn); - } - } - } while (!visit.empty()); - } - - Trace("nl-ext-cm-debug") << " check assertions..." << std::endl; - std::vector<Node> check_assertions; - for (const Node& a : passertions) - { - if (d_check_model_solved.find(a) == d_check_model_solved.end()) - { - Node av = a; - // apply the substitution to a - if (!d_check_model_vars.empty()) - { - av = av.substitute(d_check_model_vars.begin(), - d_check_model_vars.end(), - d_check_model_subs.begin(), - d_check_model_subs.end()); - av = Rewriter::rewrite(av); - } - // simple check literal - if (!simpleCheckModelLit(av)) - { - Trace("nl-ext-cm") << "...check-model : assertion failed : " << a - << std::endl; - check_assertions.push_back(av); - Trace("nl-ext-cm-debug") - << "...check-model : failed assertion, value : " << av << std::endl; - } } } - - if (!check_assertions.empty()) - { - Trace("nl-ext-cm") << "...simple check failed." << std::endl; - // TODO (#1450) check model for general case - return false; - } - Trace("nl-ext-cm") << "...simple check succeeded!" << std::endl; - - // must assert and re-check if produce models is true - if (options::produceModels()) + std::vector<Node> lemmas; + std::vector<Node> gs; + bool ret = d_model.checkModel( + passertions, false_asserts, d_taylor_degree, lemmas, gs); + for (Node& mg : gs) { - NodeManager* nm = NodeManager::currentNM(); - // model guard whose semantics is "the model we constructed holds" - Node mg = nm->mkSkolem("model", nm->booleanType()); mg = Rewriter::rewrite(mg); mg = d_containing.getValuation().ensureLiteral(mg); d_containing.getOutputChannel().requirePhase(mg, true); - // assert the constructed model as assertions - for (const std::pair<const Node, std::pair<Node, Node> > cb : - d_check_model_bounds) - { - Node l = cb.second.first; - Node u = cb.second.second; - Node v = cb.first; - Node pred = nm->mkNode(AND, nm->mkNode(GEQ, v, l), nm->mkNode(GEQ, u, v)); - pred = nm->mkNode(OR, mg.negate(), pred); - Trace("nl-ext-lemma-model") << "Assert : " << pred << std::endl; - d_containing.getOutputChannel().lemma(pred); - } d_builtModel = true; } - return true; -} - -bool NonlinearExtension::addCheckModelSubstitution(TNode v, TNode s) -{ - // should not substitute the same variable twice - Trace("nl-ext-model") << "* check model substitution : " << v << " -> " << s << std::endl; - // should not set exact bound more than once - if(std::find(d_check_model_vars.begin(),d_check_model_vars.end(),v)!=d_check_model_vars.end()) - { - Trace("nl-ext-model") << "...ERROR: already has value." << std::endl; - // this should never happen since substitutions should be applied eagerly - Assert(false); - return false; - } - // if we previously had an approximate bound, the exact bound should be in its - // range - std::map<Node, std::pair<Node, Node> >::iterator itb = - d_check_model_bounds.find(v); - if (itb != d_check_model_bounds.end()) - { - if (s.getConst<Rational>() >= itb->second.first.getConst<Rational>() - || s.getConst<Rational>() <= itb->second.second.getConst<Rational>()) - { - Trace("nl-ext-model") - << "...ERROR: already has bound which is out of range." << std::endl; - return false; - } - } - for (unsigned i = 0, size = d_check_model_subs.size(); i < size; i++) + for (Node& lem : lemmas) { - Node ms = d_check_model_subs[i]; - Node mss = ms.substitute(v, s); - if (mss != ms) - { - mss = Rewriter::rewrite(mss); - } - d_check_model_subs[i] = mss; - } - d_check_model_vars.push_back(v); - d_check_model_subs.push_back(s); - return true; -} - -bool NonlinearExtension::addCheckModelBound(TNode v, TNode l, TNode u) -{ - Trace("nl-ext-model") << "* check model bound : " << v << " -> [" << l << " " << u << "]" << std::endl; - if( l==u ) - { - // bound is exact, can add as substitution - return addCheckModelSubstitution(v, l); - } - // should not set a bound for a value that is exact - if (std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) - != d_check_model_vars.end()) - { - Trace("nl-ext-model") - << "...ERROR: setting bound for variable that already has exact value." - << std::endl; - Assert(false); - return false; - } - Assert(l.isConst()); - Assert(u.isConst()); - Assert(l.getConst<Rational>() <= u.getConst<Rational>()); - d_check_model_bounds[v] = std::pair<Node, Node>(l, u); - return true; -} - -bool NonlinearExtension::hasCheckModelAssignment(Node v) const -{ - if (d_check_model_bounds.find(v) != d_check_model_bounds.end()) - { - return true; - } - return std::find(d_check_model_vars.begin(), d_check_model_vars.end(), v) - != d_check_model_vars.end(); -} - -bool NonlinearExtension::solveEqualitySimple(Node eq) -{ - Node seq = eq; - if (!d_check_model_vars.empty()) - { - seq = eq.substitute(d_check_model_vars.begin(), - d_check_model_vars.end(), - d_check_model_subs.begin(), - d_check_model_subs.end()); - seq = Rewriter::rewrite(seq); - if (seq.isConst()) - { - if (seq.getConst<bool>()) - { - d_check_model_solved[eq] = Node::null(); - return true; - } - return false; - } - } - Trace("nl-ext-cms") << "simple solve equality " << seq << "..." << std::endl; - Assert(seq.getKind() == EQUAL); - std::map<Node, Node> msum; - if (!ArithMSum::getMonomialSumLit(seq, msum)) - { - Trace("nl-ext-cms") << "...fail, could not determine monomial sum." - << std::endl; - return false; - } - bool is_valid = true; - // the variable we will solve a quadratic equation for - Node var; - Node a = d_zero; - Node b = d_zero; - Node c = d_zero; - NodeManager* nm = NodeManager::currentNM(); - // the list of variables that occur as a monomial in msum, and whose value - // is so far unconstrained in the model. - std::unordered_set<Node, NodeHashFunction> unc_vars; - // the list of variables that occur as a factor in a monomial, and whose - // value is so far unconstrained in the model. - std::unordered_set<Node, NodeHashFunction> unc_vars_factor; - for (std::pair<const Node, Node>& m : msum) - { - Node v = m.first; - Node coeff = m.second.isNull() ? d_one : m.second; - if (v.isNull()) - { - c = coeff; - } - else if (v.getKind() == NONLINEAR_MULT) - { - if (v.getNumChildren() == 2 && v[0].isVar() && v[0] == v[1] - && (var.isNull() || var == v[0])) - { - // may solve quadratic - a = coeff; - var = v[0]; - } - else - { - is_valid = false; - Trace("nl-ext-cms-debug") - << "...invalid due to non-linear monomial " << v << std::endl; - // may wish to set an exact bound for a factor and repeat - for (const Node& vc : v) - { - unc_vars_factor.insert(vc); - } - } - } - else if (!v.isVar() || (!var.isNull() && var != v)) - { - Trace("nl-ext-cms-debug") - << "...invalid due to factor " << v << std::endl; - // cannot solve multivariate - if (is_valid) - { - is_valid = false; - // if b is non-zero, then var is also an unconstrained variable - if (b != d_zero) - { - unc_vars.insert(var); - unc_vars_factor.insert(var); - } - } - // if v is unconstrained, we may turn this equality into a substitution - unc_vars.insert(v); - unc_vars_factor.insert(v); - } - else - { - // set the variable to solve for - b = coeff; - var = v; - } - } - if (!is_valid) - { - // see if we can solve for a variable? - for (const Node& uv : unc_vars) - { - Trace("nl-ext-cm-debug") << "check subs var : " << uv << std::endl; - // cannot already have a bound - if (uv.isVar() && !hasCheckModelAssignment(uv)) - { - Node slv; - Node veqc; - if (ArithMSum::isolate(uv, msum, veqc, slv, EQUAL) != 0) - { - Assert(!slv.isNull()); - // currently do not support substitution-with-coefficients - if (veqc.isNull() && !expr::hasSubterm(slv, uv)) - { - Trace("nl-ext-cm") - << "check-model-subs : " << uv << " -> " << slv << std::endl; - bool ret = addCheckModelSubstitution(uv, slv); - if (ret) - { - Trace("nl-ext-cms") << "...success, model substitution " << uv - << " -> " << slv << std::endl; - d_check_model_solved[eq] = uv; - } - return ret; - } - } - } - } - // see if we can assign a variable to a constant - for (const Node& uvf : unc_vars_factor) - { - Trace("nl-ext-cm-debug") << "check set var : " << uvf << std::endl; - // cannot already have a bound - if (uvf.isVar() && !hasCheckModelAssignment(uvf)) - { - Node uvfv = computeModelValue(uvf); - Trace("nl-ext-cm") << "check-model-bound : exact : " << uvf << " = "; - printRationalApprox("nl-ext-cm", uvfv); - Trace("nl-ext-cm") << std::endl; - bool ret = addCheckModelSubstitution(uvf, uvfv); - // recurse - return ret ? solveEqualitySimple(eq) : false; - } - } - Trace("nl-ext-cms") << "...fail due to constrained invalid terms." - << std::endl; - return false; - } - else if (var.isNull() || var.getType().isInteger()) - { - // cannot solve quadratic equations for integer variables - Trace("nl-ext-cms") << "...fail due to variable to solve for." << std::endl; - return false; - } - - // we are linear, it is simple - if (a == d_zero) - { - if (b == d_zero) - { - Trace("nl-ext-cms") << "...fail due to zero a/b." << std::endl; - Assert(false); - return false; - } - Node val = nm->mkConst(-c.getConst<Rational>() / b.getConst<Rational>()); - Trace("nl-ext-cm") << "check-model-bound : exact : " << var << " = "; - printRationalApprox("nl-ext-cm", val); - Trace("nl-ext-cm") << std::endl; - bool ret = addCheckModelSubstitution(var, val); - if (ret) - { - Trace("nl-ext-cms") << "...success, solved linear." << std::endl; - d_check_model_solved[eq] = var; - } - return ret; - } - Trace("nl-ext-quad") << "Solve quadratic : " << seq << std::endl; - Trace("nl-ext-quad") << " a : " << a << std::endl; - Trace("nl-ext-quad") << " b : " << b << std::endl; - Trace("nl-ext-quad") << " c : " << c << std::endl; - Node two_a = nm->mkNode(MULT, d_two, a); - two_a = Rewriter::rewrite(two_a); - Node sqrt_val = nm->mkNode( - MINUS, nm->mkNode(MULT, b, b), nm->mkNode(MULT, d_two, two_a, c)); - sqrt_val = Rewriter::rewrite(sqrt_val); - Trace("nl-ext-quad") << "Will approximate sqrt " << sqrt_val << std::endl; - Assert(sqrt_val.isConst()); - // if it is negative, then we are in conflict - if (sqrt_val.getConst<Rational>().sgn() == -1) - { - Node conf = seq.negate(); - Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : quadratic no root : " - << conf << std::endl; - d_containing.getOutputChannel().lemma(conf); - Trace("nl-ext-cms") << "...fail due to negative discriminant." << std::endl; - return false; - } - if (hasCheckModelAssignment(var)) - { - Trace("nl-ext-cms") << "...fail due to bounds on variable to solve for." - << std::endl; - // two quadratic equations for same variable, give up - return false; - } - // approximate the square root of sqrt_val - Node l, u; - if (!getApproximateSqrt(sqrt_val, l, u, 15 + d_taylor_degree)) - { - Trace("nl-ext-cms") << "...fail, could not approximate sqrt." << std::endl; - return false; - } - d_used_approx = true; - Trace("nl-ext-quad") << "...got " << l << " <= sqrt(" << sqrt_val - << ") <= " << u << std::endl; - Node negb = nm->mkConst(-b.getConst<Rational>()); - Node coeffa = nm->mkConst(Rational(1) / two_a.getConst<Rational>()); - // two possible bound regions - Node bounds[2][2]; - Node diff_bound[2]; - Node m_var = computeModelValue(var, 0); - Assert(m_var.isConst()); - for (unsigned r = 0; r < 2; r++) - { - for (unsigned b = 0; b < 2; b++) - { - Node val = b == 0 ? l : u; - // (-b +- approx_sqrt( b^2 - 4ac ))/2a - Node approx = nm->mkNode( - MULT, coeffa, nm->mkNode(r == 0 ? MINUS : PLUS, negb, val)); - approx = Rewriter::rewrite(approx); - bounds[r][b] = approx; - Assert(approx.isConst()); - } - if (bounds[r][0].getConst<Rational>() > bounds[r][1].getConst<Rational>()) - { - // ensure bound is (lower, upper) - Node tmp = bounds[r][0]; - bounds[r][0] = bounds[r][1]; - bounds[r][1] = tmp; - } - Node diff = - nm->mkNode(MINUS, - m_var, - nm->mkNode(MULT, - nm->mkConst(Rational(1) / Rational(2)), - nm->mkNode(PLUS, bounds[r][0], bounds[r][1]))); - Trace("nl-ext-cm-debug") << "Bound option #" << r << " : "; - printRationalApprox("nl-ext-cm-debug", bounds[r][0]); - Trace("nl-ext-cm-debug") << "..."; - printRationalApprox("nl-ext-cm-debug", bounds[r][1]); - Trace("nl-ext-cm-debug") << std::endl; - diff = Rewriter::rewrite(diff); - Assert(diff.isConst()); - diff = nm->mkConst(diff.getConst<Rational>().abs()); - diff_bound[r] = diff; - Trace("nl-ext-cm-debug") << "...diff from model value ("; - printRationalApprox("nl-ext-cm-debug", m_var); - Trace("nl-ext-cm-debug") << ") is "; - printRationalApprox("nl-ext-cm-debug", diff_bound[r]); - Trace("nl-ext-cm-debug") << std::endl; - } - // take the one that var is closer to in the model - Node cmp = nm->mkNode(GEQ, diff_bound[0], diff_bound[1]); - cmp = Rewriter::rewrite(cmp); - Assert(cmp.isConst()); - unsigned r_use_index = cmp == d_true ? 1 : 0; - Trace("nl-ext-cm") << "check-model-bound : approximate (sqrt) : "; - printRationalApprox("nl-ext-cm", bounds[r_use_index][0]); - Trace("nl-ext-cm") << " <= " << var << " <= "; - printRationalApprox("nl-ext-cm", bounds[r_use_index][1]); - Trace("nl-ext-cm") << std::endl; - bool ret = - addCheckModelBound(var, bounds[r_use_index][0], bounds[r_use_index][1]); - if (ret) - { - d_check_model_solved[eq] = var; - Trace("nl-ext-cms") << "...success, solved quadratic." << std::endl; + Trace("nl-ext-lemma-model") + << "Lemma from check model : " << lem << std::endl; + d_containing.getOutputChannel().lemma(lem); } return ret; } -bool NonlinearExtension::simpleCheckModelLit(Node lit) -{ - Trace("nl-ext-cms") << "*** Simple check-model lit for " << lit << "..." - << std::endl; - if (lit.isConst()) - { - Trace("nl-ext-cms") << " return constant." << std::endl; - return lit.getConst<bool>(); - } - NodeManager* nm = NodeManager::currentNM(); - bool pol = lit.getKind() != kind::NOT; - Node atom = lit.getKind() == kind::NOT ? lit[0] : lit; - - if (atom.getKind() == EQUAL) - { - // x = a is ( x >= a ^ x <= a ) - for (unsigned i = 0; i < 2; i++) - { - Node lit = nm->mkNode(GEQ, atom[i], atom[1 - i]); - if (!pol) - { - lit = lit.negate(); - } - lit = Rewriter::rewrite(lit); - bool success = simpleCheckModelLit(lit); - if (success != pol) - { - // false != true -> one conjunct of equality is false, we fail - // true != false -> one disjunct of disequality is true, we succeed - return success; - } - } - // both checks passed and polarity is true, or both checks failed and - // polarity is false - return pol; - } - else if (atom.getKind() != GEQ) - { - Trace("nl-ext-cms") << " failed due to unknown literal." << std::endl; - return false; - } - // get the monomial sum - std::map<Node, Node> msum; - if (!ArithMSum::getMonomialSumLit(atom, msum)) - { - Trace("nl-ext-cms") << " failed due to get msum." << std::endl; - return false; - } - // simple interval analysis - if (simpleCheckModelMsum(msum, pol)) - { - return true; - } - // can also try reasoning about univariate quadratic equations - Trace("nl-ext-cms-debug") - << "* Try univariate quadratic analysis..." << std::endl; - std::vector<Node> vs_invalid; - std::unordered_set<Node, NodeHashFunction> vs; - std::map<Node, Node> v_a; - std::map<Node, Node> v_b; - // get coefficients... - for (std::pair<const Node, Node>& m : msum) - { - Node v = m.first; - if (!v.isNull()) - { - if (v.isVar()) - { - v_b[v] = m.second.isNull() ? d_one : m.second; - vs.insert(v); - } - else if (v.getKind() == NONLINEAR_MULT && v.getNumChildren() == 2 - && v[0] == v[1] && v[0].isVar()) - { - v_a[v[0]] = m.second.isNull() ? d_one : m.second; - vs.insert(v[0]); - } - else - { - vs_invalid.push_back(v); - } - } - } - // solve the valid variables... - Node invalid_vsum = vs_invalid.empty() ? d_zero - : (vs_invalid.size() == 1 - ? vs_invalid[0] - : nm->mkNode(PLUS, vs_invalid)); - // substitution to try - std::vector<Node> qvars; - std::vector<Node> qsubs; - for (const Node& v : vs) - { - // is it a valid variable? - std::map<Node, std::pair<Node, Node> >::iterator bit = - d_check_model_bounds.find(v); - if (!expr::hasSubterm(invalid_vsum, v) && bit != d_check_model_bounds.end()) - { - std::map<Node, Node>::iterator it = v_a.find(v); - if (it != v_a.end()) - { - Node a = it->second; - Assert(a.isConst()); - int asgn = a.getConst<Rational>().sgn(); - Assert(asgn != 0); - Node t = nm->mkNode(MULT, a, v, v); - Node b = d_zero; - it = v_b.find(v); - if (it != v_b.end()) - { - b = it->second; - t = nm->mkNode(PLUS, t, nm->mkNode(MULT, b, v)); - } - t = Rewriter::rewrite(t); - Trace("nl-ext-cms-debug") << "Trying to find min/max for quadratic " - << t << "..." << std::endl; - Trace("nl-ext-cms-debug") << " a = " << a << std::endl; - Trace("nl-ext-cms-debug") << " b = " << b << std::endl; - // find maximal/minimal value on the interval - Node apex = nm->mkNode( - DIVISION, nm->mkNode(UMINUS, b), nm->mkNode(MULT, d_two, a)); - apex = Rewriter::rewrite(apex); - Assert(apex.isConst()); - // for lower, upper, whether we are greater than the apex - bool cmp[2]; - Node boundn[2]; - for (unsigned r = 0; r < 2; r++) - { - boundn[r] = r == 0 ? bit->second.first : bit->second.second; - Node cmpn = nm->mkNode(GT, boundn[r], apex); - cmpn = Rewriter::rewrite(cmpn); - Assert(cmpn.isConst()); - cmp[r] = cmpn.getConst<bool>(); - } - Trace("nl-ext-cms-debug") << " apex " << apex << std::endl; - Trace("nl-ext-cms-debug") - << " lower " << boundn[0] << ", cmp: " << cmp[0] << std::endl; - Trace("nl-ext-cms-debug") - << " upper " << boundn[1] << ", cmp: " << cmp[1] << std::endl; - Assert(boundn[0].getConst<Rational>() - <= boundn[1].getConst<Rational>()); - Node s; - qvars.push_back(v); - if (cmp[0] != cmp[1]) - { - Assert(!cmp[0] && cmp[1]); - // does the sign match the bound? - if ((asgn == 1) == pol) - { - // the apex is the max/min value - s = apex; - Trace("nl-ext-cms-debug") << " ...set to apex." << std::endl; - } - else - { - // it is one of the endpoints, plug in and compare - Node tcmpn[2]; - for (unsigned r = 0; r < 2; r++) - { - qsubs.push_back(boundn[r]); - Node ts = t.substitute( - qvars.begin(), qvars.end(), qsubs.begin(), qsubs.end()); - tcmpn[r] = Rewriter::rewrite(ts); - qsubs.pop_back(); - } - Node tcmp = nm->mkNode(LT, tcmpn[0], tcmpn[1]); - Trace("nl-ext-cms-debug") - << " ...both sides of apex, compare " << tcmp << std::endl; - tcmp = Rewriter::rewrite(tcmp); - Assert(tcmp.isConst()); - unsigned bindex_use = (tcmp.getConst<bool>() == pol) ? 1 : 0; - Trace("nl-ext-cms-debug") - << " ...set to " << (bindex_use == 1 ? "upper" : "lower") - << std::endl; - s = boundn[bindex_use]; - } - } - else - { - // both to one side of the apex - // we figure out which bound to use (lower or upper) based on - // three factors: - // (1) whether a's sign is positive, - // (2) whether we are greater than the apex of the parabola, - // (3) the polarity of the constraint, i.e. >= or <=. - // there are 8 cases of these factors, which we test here. - unsigned bindex_use = (((asgn == 1) == cmp[0]) == pol) ? 0 : 1; - Trace("nl-ext-cms-debug") - << " ...set to " << (bindex_use == 1 ? "upper" : "lower") - << std::endl; - s = boundn[bindex_use]; - } - Assert(!s.isNull()); - qsubs.push_back(s); - Trace("nl-ext-cms") << "* set bound based on quadratic : " << v - << " -> " << s << std::endl; - } - } - } - if (!qvars.empty()) - { - Assert(qvars.size() == qsubs.size()); - Node slit = - lit.substitute(qvars.begin(), qvars.end(), qsubs.begin(), qsubs.end()); - slit = Rewriter::rewrite(slit); - return simpleCheckModelLit(slit); - } - return false; -} - -bool NonlinearExtension::simpleCheckModelMsum(const std::map<Node, Node>& msum, - bool pol) -{ - Trace("nl-ext-cms-debug") << "* Try simple interval analysis..." << std::endl; - NodeManager* nm = NodeManager::currentNM(); - // map from transcendental functions to whether they were set to lower - // bound - bool simpleSuccess = true; - std::map<Node, bool> set_bound; - std::vector<Node> sum_bound; - for (const std::pair<const Node, Node>& m : msum) - { - Node v = m.first; - if (v.isNull()) - { - sum_bound.push_back(m.second.isNull() ? d_one : m.second); - } - else - { - Trace("nl-ext-cms-debug") << "- monomial : " << v << std::endl; - // --- whether we should set a lower bound for this monomial - bool set_lower = - (m.second.isNull() || m.second.getConst<Rational>().sgn() == 1) - == pol; - Trace("nl-ext-cms-debug") - << "set bound to " << (set_lower ? "lower" : "upper") << std::endl; - - // --- Collect variables and factors in v - std::vector<Node> vars; - std::vector<unsigned> factors; - if (v.getKind() == NONLINEAR_MULT) - { - unsigned last_start = 0; - for (unsigned i = 0, nchildren = v.getNumChildren(); i < nchildren; i++) - { - // are we at the end? - if (i + 1 == nchildren || v[i + 1] != v[i]) - { - unsigned vfact = 1 + (i - last_start); - last_start = (i + 1); - vars.push_back(v[i]); - factors.push_back(vfact); - } - } - } - else - { - vars.push_back(v); - factors.push_back(1); - } - - // --- Get the lower and upper bounds and sign information. - // Whether we have an (odd) number of negative factors in vars, apart - // from the variable at choose_index. - bool has_neg_factor = false; - int choose_index = -1; - std::vector<Node> ls; - std::vector<Node> us; - // the relevant sign information for variables with odd exponents: - // 1: both signs of the interval of this variable are positive, - // -1: both signs of the interval of this variable are negative. - std::vector<int> signs; - Trace("nl-ext-cms-debug") << "get sign information..." << std::endl; - for (unsigned i = 0, size = vars.size(); i < size; i++) - { - Node vc = vars[i]; - unsigned vcfact = factors[i]; - if (Trace.isOn("nl-ext-cms-debug")) - { - Trace("nl-ext-cms-debug") << "-- " << vc; - if (vcfact > 1) - { - Trace("nl-ext-cms-debug") << "^" << vcfact; - } - Trace("nl-ext-cms-debug") << " "; - } - std::map<Node, std::pair<Node, Node> >::iterator bit = - d_check_model_bounds.find(vc); - // if there is a model bound for this term - if (bit != d_check_model_bounds.end()) - { - Node l = bit->second.first; - Node u = bit->second.second; - ls.push_back(l); - us.push_back(u); - int vsign = 0; - if (vcfact % 2 == 1) - { - vsign = 1; - int lsgn = l.getConst<Rational>().sgn(); - int usgn = u.getConst<Rational>().sgn(); - Trace("nl-ext-cms-debug") - << "bound_sign(" << lsgn << "," << usgn << ") "; - if (lsgn == -1) - { - if (usgn < 1) - { - // must have a negative factor - has_neg_factor = !has_neg_factor; - vsign = -1; - } - else if (choose_index == -1) - { - // set the choose index to this - choose_index = i; - vsign = 0; - } - else - { - // ambiguous, can't determine the bound - Trace("nl-ext-cms") - << " failed due to ambiguious monomial." << std::endl; - return false; - } - } - } - Trace("nl-ext-cms-debug") << " -> " << vsign << std::endl; - signs.push_back(vsign); - } - else - { - Trace("nl-ext-cms-debug") << std::endl; - Trace("nl-ext-cms") - << " failed due to unknown bound for " << vc << std::endl; - // should either assign a model bound or eliminate the variable - // via substitution - Assert(false); - return false; - } - } - // whether we will try to minimize/maximize (-1/1) the absolute value - int setAbs = (set_lower == has_neg_factor) ? 1 : -1; - Trace("nl-ext-cms-debug") - << "set absolute value to " << (setAbs == 1 ? "maximal" : "minimal") - << std::endl; - - std::vector<Node> vbs; - Trace("nl-ext-cms-debug") << "set bounds..." << std::endl; - for (unsigned i = 0, size = vars.size(); i < size; i++) - { - Node vc = vars[i]; - unsigned vcfact = factors[i]; - Node l = ls[i]; - Node u = us[i]; - bool vc_set_lower; - int vcsign = signs[i]; - Trace("nl-ext-cms-debug") - << "Bounds for " << vc << " : " << l << ", " << u - << ", sign : " << vcsign << ", factor : " << vcfact << std::endl; - if (l == u) - { - // by convention, always say it is lower if they are the same - vc_set_lower = true; - Trace("nl-ext-cms-debug") - << "..." << vc << " equal bound, set to lower" << std::endl; - } - else - { - if (vcfact % 2 == 0) - { - // minimize or maximize its absolute value - Rational la = l.getConst<Rational>().abs(); - Rational ua = u.getConst<Rational>().abs(); - if (la == ua) - { - // by convention, always say it is lower if abs are the same - vc_set_lower = true; - Trace("nl-ext-cms-debug") - << "..." << vc << " equal abs, set to lower" << std::endl; - } - else - { - vc_set_lower = (la > ua) == (setAbs == 1); - } - } - else if (signs[i] == 0) - { - // we choose this index to match the overall set_lower - vc_set_lower = set_lower; - } - else - { - vc_set_lower = (signs[i] != setAbs); - } - Trace("nl-ext-cms-debug") - << "..." << vc << " set to " << (vc_set_lower ? "lower" : "upper") - << std::endl; - } - // check whether this is a conflicting bound - std::map<Node, bool>::iterator itsb = set_bound.find(vc); - if (itsb == set_bound.end()) - { - set_bound[vc] = vc_set_lower; - } - else if (itsb->second != vc_set_lower) - { - Trace("nl-ext-cms") - << " failed due to conflicting bound for " << vc << std::endl; - return false; - } - // must over/under approximate based on vc_set_lower, computed above - Node vb = vc_set_lower ? l : u; - for (unsigned i = 0; i < vcfact; i++) - { - vbs.push_back(vb); - } - } - if (!simpleSuccess) - { - break; - } - Node vbound = vbs.size() == 1 ? vbs[0] : nm->mkNode(MULT, vbs); - sum_bound.push_back(ArithMSum::mkCoeffTerm(m.second, vbound)); - } - } - // if the exact bound was computed via simple analysis above - // make the bound - Node bound; - if (sum_bound.size() > 1) - { - bound = nm->mkNode(kind::PLUS, sum_bound); - } - else if (sum_bound.size() == 1) - { - bound = sum_bound[0]; - } - else - { - bound = d_zero; - } - // make the comparison - Node comp = nm->mkNode(kind::GEQ, bound, d_zero); - if (!pol) - { - comp = comp.negate(); - } - Trace("nl-ext-cms") << " comparison is : " << comp << std::endl; - comp = Rewriter::rewrite(comp); - Assert(comp.isConst()); - Trace("nl-ext-cms") << " returned : " << comp << std::endl; - return comp == d_true; -} std::vector<Node> NonlinearExtension::checkSplitZero() { std::vector<Node> lemmas; @@ -1849,11 +885,9 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, for (unsigned i = 0, xsize = xts.size(); i < xsize; i++) { Node a = xts[i]; - computeModelValue(a, 0); - computeModelValue(a, 1); - printModelValue("nl-ext-mv", a); - //Assert(d_mv[1][a].isConst()); - //Assert(d_mv[0][a].isConst()); + d_model.computeConcreteModelValue(a); + d_model.computeAbstractModelValue(a); + d_model.printModelValue("nl-ext-mv", a); Kind ak = a.getKind(); if (ak == NONLINEAR_MULT) { @@ -1868,25 +902,12 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, if (!IsInVector(d_ms_vars, itvl->second[k])) { d_ms_vars.push_back(itvl->second[k]); } - Node mvk = computeModelValue( itvl->second[k], 1 ); + 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) - std::map< Node, std::map< Node, unsigned > >::iterator itme = - d_m_exp.find( a ); Assert( itme!=d_m_exp.end() ); for( std::map< - Node, unsigned >::iterator itme2 = itme->second.begin(); itme2 != - itme->second.end(); ++itme2 ){ Node v = itme->first; Assert( - d_mv[0].find( v )!=d_mv[0].end() ); Node mvv = d_mv[0][ v ]; if( - mvv==d_one || mvv==d_neg_one ){ ms_proc[ a ] = true; - Trace("nl-ext-mv") - << "...mark " << a << " reduced since has 1 factor." << std::endl; - break; - } - } - */ + // mark processed if has a "one" factor (will look at reduced monomial)? } else if (a.getNumChildren() > 0) { @@ -1932,7 +953,9 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, { // apply congruence to pairs of terms that are disequal and congruent Assert(aa.getNumChildren() == a.getNumChildren()); - if (d_mv[1][a] != d_mv[1][aa]) + Node mvaa = d_model.computeAbstractModelValue(a); + Node mvaaa = d_model.computeAbstractModelValue(aa); + if (mvaa != mvaaa) { std::vector<Node> exp; for (unsigned j = 0, size = a.getNumChildren(); j < size; j++) @@ -2029,8 +1052,8 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, registerMonomial(d_one); for (unsigned j = 0; j < d_order_points.size(); j++) { Node c = d_order_points[j]; - computeModelValue(c, 0); - computeModelValue(c, 1); + d_model.computeConcreteModelValue(c); + d_model.computeAbstractModelValue(c); } // register variables @@ -2038,9 +1061,9 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, for (unsigned i = 0; i < d_ms_vars.size(); i++) { Node v = d_ms_vars[i]; registerMonomial(v); - computeModelValue(v, 0); - computeModelValue(v, 1); - printModelValue("nl-ext-mv", v); + d_model.computeConcreteModelValue(v); + d_model.computeAbstractModelValue(v); + d_model.printModelValue("nl-ext-mv", v); } if (Trace.isOn("nl-ext-mv")) { @@ -2054,9 +1077,9 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, for (const Node& tf : tfl.second) { Node v = tf[0]; - computeModelValue(v, 0); - computeModelValue(v, 1); - printModelValue("nl-ext-mv", v); + d_model.computeConcreteModelValue(v); + d_model.computeAbstractModelValue(v); + d_model.printModelValue("nl-ext-mv", v); } } } @@ -2099,14 +1122,15 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, //-----------------------------------lemmas based on magnitude of non-zero monomials Trace("nl-ext-proc") << "Assign order ids..." << std::endl; - unsigned r = 3; - assignOrderIds(d_ms_vars, d_order_vars, r); + // 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; - SortNonlinearExtension smv; - smv.d_nla = this; - smv.d_order_type = r; + 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); @@ -2126,10 +1150,7 @@ int NonlinearExtension::checkLastCall(const std::vector<Node>& assertions, // sort monomials by degree Trace("nl-ext-proc") << "Sort monomials by degree..." << std::endl; - SortNonlinearExtension snlad; - snlad.d_nla = this; - snlad.d_order_type = 4; - snlad.d_reverse_order = false; + 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()); @@ -2239,8 +1260,8 @@ void NonlinearExtension::check(Theory::Effort e) { getAssertions(assertions); // reset cached information - d_mv[0].clear(); - d_mv[1].clear(); + TheoryModel* tm = d_containing.getValuation().getModel(); + d_model.reset(tm); Trace("nl-ext-mv-assert") << "Getting model values... check for [model-false]" << std::endl; @@ -2279,9 +1300,9 @@ void NonlinearExtension::check(Theory::Effort e) { { TNode shared_term = *its; // compute its value in the model, and its evaluation in the model - Node stv0 = computeModelValue(shared_term, 0); - Node stv1 = computeModelValue(shared_term, 1); - printModelValue("nl-ext-mv", shared_term); + Node stv0 = d_model.computeConcreteModelValue(shared_term); + Node stv1 = d_model.computeAbstractModelValue(shared_term); + d_model.printModelValue("nl-ext-mv", shared_term); if (stv0 != stv1) { num_shared_wrong_value++; @@ -2311,7 +1332,7 @@ void NonlinearExtension::check(Theory::Effort e) { bool needsRecheck; do { - d_used_approx = false; + d_model.resetCheck(); needsRecheck = false; Assert(e == Theory::EFFORT_LAST_CALL); // complete_status: @@ -2392,10 +1413,9 @@ void NonlinearExtension::check(Theory::Effort e) { } // we are incomplete - if (options::nlExtIncPrecision() && d_used_approx) + if (options::nlExtIncPrecision() && d_model.usedApproximate()) { d_taylor_degree++; - d_used_approx = false; needsRecheck = true; // increase precision for PI? // Difficult since Taylor series is very slow to converge @@ -2423,49 +1443,8 @@ void NonlinearExtension::check(Theory::Effort e) { // don't need to build the model yet return; } - // Record the approximations we used. This code calls the - // recordApproximation method of the model, which overrides the model - // values for variables that we solved for, using techniques specific to - // this class. - NodeManager* nm = NodeManager::currentNM(); - TheoryModel* m = d_containing.getValuation().getModel(); - for (const std::pair<const Node, std::pair<Node, Node> >& cb : - d_check_model_bounds) - { - Node l = cb.second.first; - Node u = cb.second.second; - Node pred; - Node v = cb.first; - if (l != u) - { - pred = nm->mkNode(AND, nm->mkNode(GEQ, v, l), nm->mkNode(GEQ, u, v)); - } - else if (!m->areEqual(v, l)) - { - // only record if value was not equal already - pred = v.eqNode(l); - } - if (!pred.isNull()) - { - pred = Rewriter::rewrite(pred); - m->recordApproximation(v, pred); - } - } - // Also record the exact values we used. An exact value can be seen as a - // special kind approximation of the form (choice x. x = exact_value). - // Notice that the above term gets rewritten such that the choice function - // is eliminated. - for (size_t i = 0, num = d_check_model_vars.size(); i < num; i++) - { - Node v = d_check_model_vars[i]; - Node s = d_check_model_subs[i]; - if (!m->areEqual(v, s)) - { - Node pred = v.eqNode(s); - pred = Rewriter::rewrite(pred); - m->recordApproximation(v, pred); - } - } + // record approximations in the model + d_model.recordApproximations(); return; } } @@ -2477,21 +1456,24 @@ void NonlinearExtension::presolve() void NonlinearExtension::assignOrderIds(std::vector<Node>& vars, NodeMultiset& order, - unsigned orderType) { - SortNonlinearExtension smv; - smv.d_nla = this; - smv.d_order_type = orderType; + 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 = (orderType == 0 || orderType == 2) ? 0 : 1; + unsigned order_index = isConcrete ? 0 : 1; Node prev; for (unsigned j = 0; j < vars.size(); j++) { Node x = vars[j]; - Node v = get_compare_value(x, orderType); + 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) @@ -2504,12 +1486,13 @@ void NonlinearExtension::assignOrderIds(std::vector<Node>& vars, do { success = false; if (order_index < d_order_points.size()) { - Node vv = get_compare_value(d_order_points[order_index], orderType); - if (compare_value(v, vv, orderType) <= 0) { + 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_" << orderType << "[" << d_order_points[order_index] - << "] = " << counter << std::endl; + Trace("nl-ext-mvo") << "O[" << d_order_points[order_index] + << "] = " << counter << std::endl; order[d_order_points[order_index]] = counter; prev = vv; order_index++; @@ -2518,54 +1501,23 @@ void NonlinearExtension::assignOrderIds(std::vector<Node>& vars, } } while (success); } - if (prev.isNull() || compare_value(v, prev, orderType) != 0) { + if (prev.isNull() || d_model.compareValue(v, prev, isAbsolute) != 0) + { counter++; } - Trace("nl-ext-mvo") << "O_" << orderType << "[" << x << "] = " << counter - << std::endl; + 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_" << orderType << "[" - << d_order_points[order_index] << "] = " << counter - << std::endl; + Trace("nl-ext-mvo") << "O[" << d_order_points[order_index] + << "] = " << counter << std::endl; order[d_order_points[order_index]] = counter; order_index++; } } -int NonlinearExtension::compare(Node i, Node j, unsigned orderType) const { - Assert(orderType >= 0); - if (orderType <= 3) { - Node ci = get_compare_value(i, orderType); - Node cj = get_compare_value(j, orderType); - if( ci.isConst() ){ - if( cj.isConst() ){ - return compare_value(ci, cj, orderType); - }else{ - return 1; - } - }else{ - if( cj.isConst() ){ - return -1; - }else{ - return 0; - } - } - // minimal degree - } else if (orderType == 4) { - unsigned i_count = getCount(d_m_degree, i); - unsigned j_count = getCount(d_m_degree, j); - return i_count == j_count ? 0 : (i_count < j_count ? 1 : -1); - } else { - return 0; - } -} - - - void NonlinearExtension::mkPi(){ if( d_pi.isNull() ){ d_pi = NodeManager::currentNM()->mkNullaryOperator( @@ -2641,80 +1593,16 @@ bool NonlinearExtension::getApproximateSqrt(Node c, return true; } -void NonlinearExtension::printModelValue(const char* c, - Node n, - unsigned prec) const -{ - if (Trace.isOn(c)) - { - Trace(c) << " " << n << " -> "; - for (unsigned i = 0; i < 2; i++) - { - std::map<Node, Node>::const_iterator it = d_mv[1 - i].find(n); - Assert(it != d_mv[1 - i].end()); - if (it->second.isConst()) - { - printRationalApprox(c, it->second, prec); - } - else - { - Trace(c) << "?"; // it->second; - } - Trace(c) << (i == 0 ? " [actual: " : " ]"); - } - Trace(c) << std::endl; - } -} - -int NonlinearExtension::compare_value(Node i, Node j, - unsigned orderType) const { - Assert(orderType >= 0 && orderType <= 3); - Assert(i.isConst() && j.isConst()); - Trace("nl-ext-debug") << "compare value " << i << " " << j - << ", o = " << orderType << std::endl; - int ret; - if (i == j) { - ret = 0; - } else if (orderType == 0 || orderType == 2) { - ret = i.getConst<Rational>() < j.getConst<Rational>() ? 1 : -1; - } else { - Trace("nl-ext-debug") << "vals : " << i.getConst<Rational>() << " " - << j.getConst<Rational>() << std::endl; - Trace("nl-ext-debug") << i.getConst<Rational>().abs() << " " - << j.getConst<Rational>().abs() << std::endl; - ret = (i.getConst<Rational>().abs() == j.getConst<Rational>().abs() - ? 0 - : (i.getConst<Rational>().abs() < j.getConst<Rational>().abs() - ? 1 - : -1)); - } - Trace("nl-ext-debug") << "...return " << ret << std::endl; - return ret; -} - -Node NonlinearExtension::get_compare_value(Node i, unsigned orderType) const { - if (i.isConst()) - { - return i; - } - Trace("nl-ext-debug") << "Compare variable " << i << " " << orderType - << std::endl; - Assert(orderType >= 0 && orderType <= 3); - unsigned mindex = orderType <= 1 ? 0 : 1; - std::map<Node, Node>::const_iterator iti = d_mv[mindex].find(i); - Assert(iti != d_mv[mindex].end()); - return iti->second; -} - // 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<Node>& exp, std::vector<Node>& lem) { Trace("nl-ext-debug") << "Process " << a << " at index " << a_index << ", status is " << status << std::endl; - Assert(d_mv[1].find(oa) != d_mv[1].end()); + Node mvaoa = d_model.computeAbstractModelValue(oa); if (a_index == d_m_vlist[a].size()) { - if (d_mv[1][oa].getConst<Rational>().sgn() != status) { + if (mvaoa.getConst<Rational>().sgn() != status) + { Node lemma = safeConstructNary(AND, exp).impNode(mkLit(oa, d_zero, status * 2)); lem.push_back(lemma); @@ -2725,13 +1613,13 @@ int NonlinearExtension::compareSign(Node oa, Node a, unsigned a_index, Node av = d_m_vlist[a][a_index]; unsigned aexp = d_m_exp[a][av]; // take current sign in model - Assert(d_mv[1].find(av) != d_mv[0].end()); - Assert(d_mv[1][av].isConst()); - int sgn = d_mv[1][av].getConst<Rational>().sgn(); + Node mvaav = d_model.computeAbstractModelValue(av); + int sgn = mvaav.getConst<Rational>().sgn(); Trace("nl-ext-debug") << "Process var " << av << "^" << aexp << ", model sign = " << sgn << std::endl; if (sgn == 0) { - if (d_mv[1][oa].getConst<Rational>().sgn() != 0) { + if (mvaoa.getConst<Rational>().sgn() != 0) + { Node lemma = av.eqNode(d_zero).impNode(oa.eqNode(d_zero)); lem.push_back(lemma); } @@ -2806,8 +1694,8 @@ bool NonlinearExtension::compareMonomial( << " " << 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 abstract values - int modelStatus = compare(oa, ob, 3) * -2; + // 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; @@ -2821,7 +1709,7 @@ bool NonlinearExtension::compareMonomial( } } Node clem = NodeManager::currentNM()->mkNode( - IMPLIES, safeConstructNary(AND, exp), mkLit(oa, ob, status, 1)); + 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; @@ -2868,7 +1756,7 @@ bool NonlinearExtension::compareMonomial( 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, 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); @@ -2881,7 +1769,7 @@ bool NonlinearExtension::compareMonomial( 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, 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); @@ -2896,7 +1784,7 @@ bool NonlinearExtension::compareMonomial( 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, 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); @@ -2907,7 +1795,7 @@ bool NonlinearExtension::compareMonomial( 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, 1)); + 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); @@ -2919,7 +1807,7 @@ bool NonlinearExtension::compareMonomial( 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, 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); @@ -2988,7 +1876,12 @@ std::vector<Node> NonlinearExtension::checkMonomialSign() { Node a = d_ms[j]; if (d_ms_proc.find(a) == d_ms_proc.end()) { std::vector<Node> exp; - Trace("nl-ext-debug") << " process " << a << ", mv=" << d_mv[0][a] << "..." << std::endl; + 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) { @@ -3032,8 +1925,9 @@ std::vector<Node> NonlinearExtension::checkMonomialMagnitude( unsigned c ) { // compare magnitude against variables for (unsigned k = 0; k < d_ms_vars.size(); k++) { Node v = d_ms_vars[k]; - Assert(d_mv[0].find(v) != d_mv[0].end()); - if( d_mv[0][v].isConst() ){ + Node mvcv = d_model.computeConcreteModelValue(v); + if (mvcv.isConst()) + { std::vector<Node> exp; NodeMultiset a_exp_proc; NodeMultiset b_exp_proc; @@ -3172,8 +2066,8 @@ std::vector<Node> NonlinearExtension::checkTangentPlanes() { dproc[a][b] = true; Trace("nl-ext-tplanes") << " decomposable into : " << a << " * " << b << std::endl; - Node a_v_c = computeModelValue(a, 1); - Node b_v_c = computeModelValue(b, 1); + 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 ); @@ -3320,7 +2214,7 @@ std::vector<Node> NonlinearExtension::checkMonomialInferBounds( // model Node lhs = ArithMSum::mkCoeffTerm(coeff, x); Node query = NodeManager::currentNM()->mkNode(GT, lhs, rhs); - Node query_mv = computeModelValue(query, 1); + Node query_mv = d_model.computeAbstractModelValue(query); if (query_mv == d_true) { exp = query; type = GT; @@ -3429,7 +2323,7 @@ std::vector<Node> NonlinearExtension::checkMonomialInferBounds( Node mult = d_m_contain_mult[x][y]; // x <k> t => m*x <k'> t where y = m*x // get the sign of mult - Node mmv = computeModelValue(mult); + Node mmv = d_model.computeConcreteModelValue(mult); Trace("nl-ext-bound-debug2") << "Model value of " << mult << " is " << mmv << std::endl; if(mmv.isConst()){ @@ -3454,7 +2348,7 @@ std::vector<Node> NonlinearExtension::checkMonomialInferBounds( Trace("nl-ext-bound-debug2") << " ...rewritten : " << infer << std::endl; // check whether it is false in model for abstraction - Node infer_mv = computeModelValue(infer, 1); + Node infer_mv = d_model.computeAbstractModelValue(infer); Trace("nl-ext-bound-debug") << " ...infer model value is " << infer_mv << std::endl; @@ -3509,7 +2403,7 @@ std::vector<Node> NonlinearExtension::checkFactoring( { bool polarity = lit.getKind() != NOT; Node atom = lit.getKind() == NOT ? lit[0] : lit; - Node litv = computeModelValue(lit); + Node litv = d_model.computeConcreteModelValue(lit); bool considerLit = false; if( d_skolem_atoms.find(atom) != d_skolem_atoms.end() ) { @@ -3649,11 +2543,11 @@ std::vector<Node> NonlinearExtension::checkMonomialInferResBounds() { if (ita != d_mono_diff[a].end()) { std::map<Node, Node>::iterator itb = d_mono_diff[b].find(a); Assert(itb != d_mono_diff[b].end()); - Node mv_a = computeModelValue(ita->second, 1); + Node mv_a = d_model.computeAbstractModelValue(ita->second); Assert(mv_a.isConst()); int mv_a_sgn = mv_a.getConst<Rational>().sgn(); Assert(mv_a_sgn != 0); - Node mv_b = computeModelValue(itb->second, 1); + Node mv_b = d_model.computeAbstractModelValue(itb->second); Assert(mv_b.isConst()); int mv_b_sgn = mv_b.getConst<Rational>().sgn(); Assert(mv_b_sgn != 0); @@ -3737,7 +2631,8 @@ std::vector<Node> NonlinearExtension::checkMonomialInferResBounds() { { Node conc = NodeManager::currentNM()->mkNode( jk, rhs_a_res, rhs_b_res); - Node conc_mv = computeModelValue(conc, 1); + Node conc_mv = + d_model.computeAbstractModelValue(conc); if (conc_mv == d_false) { Node rblem = NodeManager::currentNM()->mkNode( IMPLIES, @@ -3781,7 +2676,6 @@ std::vector<Node> NonlinearExtension::checkTranscendentalInitialRefine() { Kind k = tfl.first; for (const Node& t : tfl.second) { - Assert(d_mv[1].find(t) != d_mv[1].end()); //initial refinements if( d_tf_initial_refine.find( t )==d_tf_initial_refine.end() ){ d_tf_initial_refine[t] = true; @@ -3895,9 +2789,8 @@ std::vector<Node> NonlinearExtension::checkTranscendentalMonotonic() { for (const Node& tf : tfl.second) { Node a = tf[0]; - computeModelValue(a, 1); - Assert(d_mv[1].find(a) != d_mv[1].end()); - if (d_mv[1][a].isConst()) + Node mvaa = d_model.computeAbstractModelValue(a); + if (mvaa.isConst()) { Trace("nl-ext-tf-mono-debug") << "...tf term : " << a << std::endl; sorted_tf_args[k].push_back(a); @@ -3906,11 +2799,11 @@ std::vector<Node> NonlinearExtension::checkTranscendentalMonotonic() { } } } - - SortNonlinearExtension smv; - smv.d_nla = this; + + SortNlModel smv; + smv.d_nlm = &d_model; //sort by concrete values - smv.d_order_type = 0; + smv.d_isConcrete = true; smv.d_reverse_order = true; for (std::pair<const Kind, std::vector<Node> >& tfl : d_f_map) { @@ -3921,11 +2814,12 @@ std::vector<Node> NonlinearExtension::checkTranscendentalMonotonic() { for (unsigned i = 0; i < sorted_tf_args[k].size(); i++) { Node targ = sorted_tf_args[k][i]; - Assert(d_mv[1].find(targ) != d_mv[1].end()); - Trace("nl-ext-tf-mono") << " " << targ << " -> " << d_mv[1][targ] << std::endl; + Node mvatarg = d_model.computeAbstractModelValue(targ); + Trace("nl-ext-tf-mono") + << " " << targ << " -> " << mvatarg << std::endl; Node t = tf_arg_to_term[k][targ]; - Assert(d_mv[1].find(t) != d_mv[1].end()); - Trace("nl-ext-tf-mono") << " f-val : " << d_mv[1][t] << std::endl; + Node mvat = d_model.computeAbstractModelValue(t); + Trace("nl-ext-tf-mono") << " f-val : " << mvat << std::endl; } std::vector< Node > mpoints; std::vector< Node > mpoints_vals; @@ -3946,7 +2840,7 @@ std::vector<Node> NonlinearExtension::checkTranscendentalMonotonic() { for( unsigned i=0; i<mpoints.size(); i++ ){ Node mpv; if( !mpoints[i].isNull() ){ - mpv = computeModelValue( mpoints[i], 1 ); + mpv = d_model.computeAbstractModelValue(mpoints[i]); Assert(mpv.isConst()); } mpoints_vals.push_back( mpv ); @@ -3959,12 +2853,10 @@ std::vector<Node> NonlinearExtension::checkTranscendentalMonotonic() { for (unsigned i = 0, size = sorted_tf_args[k].size(); i < size; i++) { Node sarg = sorted_tf_args[k][i]; - Assert(d_mv[1].find(sarg) != d_mv[1].end()); - Node sargval = d_mv[1][sarg]; + Node sargval = d_model.computeAbstractModelValue(sarg); Assert(sargval.isConst()); Node s = tf_arg_to_term[k][ sarg ]; - Assert(d_mv[1].find(s) != d_mv[1].end()); - Node sval = d_mv[1][s]; + Node sval = d_model.computeAbstractModelValue(s); Assert(sval.isConst()); //increment to the proper monotonicity region @@ -4074,7 +2966,7 @@ std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes() for (const Node& tf : tfs.second) { // tf is Figure 3 : tf( x ) - if (isRefineableTfFun(tf)) + if (d_model.isRefineableTfFun(tf)) { Trace("nl-ext-tftp") << "Compute tangent planes " << tf << std::endl; // go until max degree is reached, or we don't meet bound criteria @@ -4100,34 +2992,11 @@ std::vector<Node> NonlinearExtension::checkTranscendentalTangentPlanes() return lemmas; } -bool NonlinearExtension::isRefineableTfFun(Node tf) -{ - Assert(tf.getKind() == SINE || tf.getKind() == EXPONENTIAL); - if (tf.getKind() == SINE) - { - // we do not consider e.g. sin( -1*x ), since considering sin( x ) will - // have the same effect - if (!tf[0].isVar()) - { - return false; - } - } - // Figure 3 : c - Node c = computeModelValue(tf[0], 1); - Assert(c.isConst()); - int csign = c.getConst<Rational>().sgn(); - if (csign == 0) - { - return false; - } - return true; -} - bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, unsigned d, std::vector<Node>& lemmas) { - Assert(isRefineableTfFun(tf)); + Assert(d_model.isRefineableTfFun(tf)); NodeManager* nm = NodeManager::currentNM(); Kind k = tf.getKind(); @@ -4142,12 +3011,12 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, poly_approx_bounds[1][-1] = pbounds[3]; // Figure 3 : c - Node c = computeModelValue(tf[0], 1); + Node c = d_model.computeAbstractModelValue(tf[0]); int csign = c.getConst<Rational>().sgn(); Assert(csign == 1 || csign == -1); // Figure 3 : v - Node v = computeModelValue(tf, 1); + Node v = d_model.computeAbstractModelValue(tf); // check value of tf Trace("nl-ext-tftp-debug") << "Process tangent plane refinement for " << tf @@ -4298,7 +3167,7 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, lem = Rewriter::rewrite(lem); Trace("nl-ext-tftp-lemma") << "*** Tangent plane lemma : " << lem << std::endl; - Assert(computeModelValue(lem, 1) == d_false); + Assert(d_model.computeAbstractModelValue(lem) == d_false); // Figure 3 : line 9 lemmas.push_back(lem); } @@ -4313,9 +3182,9 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, // insert into the vector d_secant_points[tf][d].push_back(c); // sort - SortNonlinearExtension smv; - smv.d_nla = this; - smv.d_order_type = 0; + SortNlModel smv; + smv.d_nlm = &d_model; + smv.d_isConcrete = true; std::sort( d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), smv); // get the resulting index of c @@ -4369,7 +3238,7 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, Assert(!poly_approx.isNull()); Assert(!bounds[s].isNull()); // take the model value of l or u (since may contain PI) - Node b = computeModelValue(bounds[s], 1); + Node b = d_model.computeAbstractModelValue(bounds[s]); Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds[s] << " is " << b << std::endl; Assert(b.isConst()); @@ -4425,7 +3294,7 @@ bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, << std::endl; // Figure 3 : line 22 lemmas.push_back(lem); - Assert(computeModelValue(lem, 1) == d_false); + Assert(d_model.computeAbstractModelValue(lem) == d_false); } } } @@ -4830,7 +3699,7 @@ void NonlinearExtension::getPolynomialApproximationBoundForArg( std::pair<Node, Node> NonlinearExtension::getTfModelBounds(Node tf, unsigned d) { // compute the model value of the argument - Node c = computeModelValue(tf[0], 1); + Node c = d_model.computeAbstractModelValue(tf[0]); Assert(c.isConst()); int csign = c.getConst<Rational>().sgn(); Assert(csign != 0); @@ -4851,7 +3720,7 @@ std::pair<Node, Node> NonlinearExtension::getTfModelBounds(Node tf, unsigned d) // { x -> tf[0] } pab = pab.substitute(tfv, tfs); pab = Rewriter::rewrite(pab); - Node v_pab = computeModelValue(pab, 1); + Node v_pab = d_model.computeAbstractModelValue(pab); bounds.push_back(v_pab); } else diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h index 7bed514cd..b76877414 100644 --- a/src/theory/arith/nonlinear_extension.h +++ b/src/theory/arith/nonlinear_extension.h @@ -34,6 +34,7 @@ #include "context/context.h" #include "expr/kind.h" #include "expr/node.h" +#include "theory/arith/nl_model.h" #include "theory/arith/theory_arith.h" #include "theory/uf/equality_engine.h" @@ -129,23 +130,6 @@ class NonlinearExtension { * on the output channel of TheoryArith in this function. */ void presolve(); - /** Compare arithmetic terms i and j based an ordering. - * - * orderType = 0 : compare concrete model values - * orderType = 1 : compare abstract model values - * orderType = 2 : compare abs of concrete model values - * orderType = 3 : compare abs of abstract model values - * TODO (#1287) make this an enum? - * - * For definitions of concrete vs abstract model values, - * see computeModelValue below. - */ - int compare(Node i, Node j, unsigned orderType) const; - /** Compare constant rationals i and j based an ordering. - * orderType is the same as above. - */ - int compare_value(Node i, Node j, unsigned orderType) const; - private: /** returns true if the multiset containing the * factors of monomial a is a subset of the multiset @@ -191,7 +175,7 @@ class NonlinearExtension { const std::vector<Node>& xts); //---------------------------------------term utilities static bool isArithKind(Kind k); - static Node mkLit(Node a, Node b, int status, int orderType = 0); + static Node mkLit(Node a, Node b, int status, bool isAbsolute = false); static Node mkAbs(Node a); static Node mkValidPhase(Node a, Node pi); static Node mkBounded( Node l, Node a, Node u ); @@ -203,35 +187,11 @@ class NonlinearExtension { void setMonomialFactor(Node a, Node b, const NodeMultiset& common); void registerConstraint(Node atom); - /** compute model value - * - * This computes model values for terms based on two semantics, a "concrete" - * semantics and an "abstract" semantics. - * - * index = 0 means compute the value of n based on its children recursively. - * (we call this its "concrete" value) - * index = 1 means lookup the value of n in the model. - * (we call this its "abstract" value) - * In other words, index = 1 treats multiplication terms and transcendental - * function applications as variables, whereas index = 0 computes their - * actual values. This is a key distinction used in the model-based - * refinement scheme in Cimatti et al. TACAS 2017. - * - * For example, if M( a ) = 2, M( b ) = 3, M( a * b ) = 5, then : - * - * computeModelValue( a*b, 0 ) = - * computeModelValue( a, 0 )*computeModelValue( b, 0 ) = 2*3 = 6 - * whereas: - * computeModelValue( a*b, 1 ) = 5 - */ - Node computeModelValue(Node n, unsigned index = 0); - /** returns the Node corresponding to the value of i in the - * type of order orderType, which is one of values - * described above ::compare(...). - */ - Node get_compare_value(Node i, unsigned orderType) const; - void assignOrderIds(std::vector<Node>& vars, NodeMultiset& d_order, - unsigned orderType); + /** assign order ids */ + void assignOrderIds(std::vector<Node>& vars, + NodeMultiset& d_order, + bool isConcrete, + bool isAbsolute); /** get assertions * @@ -269,98 +229,6 @@ class NonlinearExtension { */ bool checkModel(const std::vector<Node>& assertions, const std::vector<Node>& false_asserts); - - /** solve equality simple - * - * This method is used during checkModel(...). It takes as input an - * equality eq. If it returns true, then eq is correct-by-construction based - * on the information stored in our model representation (see - * d_check_model_vars, d_check_model_subs, d_check_model_bounds), and eq - * is added to d_check_model_solved. - */ - bool solveEqualitySimple(Node eq); - - /** simple check model for transcendental functions for literal - * - * This method returns true if literal is true for all interpretations of - * transcendental functions within their error bounds (as stored - * in d_check_model_bounds). This is determined by a simple under/over - * approximation of the value of sum of (linear) monomials. For example, - * if we determine that .8 < sin( 1 ) < .9, this function will return - * true for literals like: - * 2.0*sin( 1 ) > 1.5 - * -1.0*sin( 1 ) < -0.79 - * -1.0*sin( 1 ) > -0.91 - * sin( 1 )*sin( 1 ) + sin( 1 ) > 0.0 - * It will return false for literals like: - * sin( 1 ) > 0.85 - * It will also return false for literals like: - * -0.3*sin( 1 )*sin( 2 ) + sin( 2 ) > .7 - * sin( sin( 1 ) ) > .5 - * since the bounds on these terms cannot quickly be determined. - */ - bool simpleCheckModelLit(Node lit); - bool simpleCheckModelMsum(const std::map<Node, Node>& msum, bool pol); - /** - * A substitution from variables that appear in assertions to a solved form - * term. These vectors are ordered in the form: - * x_1 -> t_1 ... x_n -> t_n - * where x_i is not in the free variables of t_j for j>=i. - */ - std::vector<Node> d_check_model_vars; - std::vector<Node> d_check_model_subs; - /** add check model substitution - * - * Adds the model substitution v -> s. This applies the substitution - * { v -> s } to each term in d_check_model_subs and adds v,s to - * d_check_model_vars and d_check_model_subs respectively. - * If this method returns false, then the substitution v -> s is inconsistent - * with the current substitution and bounds. - */ - bool addCheckModelSubstitution(TNode v, TNode s); - /** lower and upper bounds for check model - * - * For each term t in the domain of this map, if this stores the pair - * (c_l, c_u) then the model M is such that c_l <= M( t ) <= c_u. - * - * We add terms whose value is approximated in the model to this map, which - * includes: - * (1) applications of transcendental functions, whose value is approximated - * by the Taylor series, - * (2) variables we have solved quadratic equations for, whose value - * involves approximations of square roots. - */ - std::map<Node, std::pair<Node, Node> > d_check_model_bounds; - /** add check model bound - * - * Adds the bound x -> < l, u > to the map above, and records the - * approximation ( x, l <= x <= u ) in the model. This method returns false - * if the bound is inconsistent with the current model substitution or - * bounds. - */ - bool addCheckModelBound(TNode v, TNode l, TNode u); - /** - * The map from literals that our model construction solved, to the variable - * that was solved for. Examples of such literals are: - * (1) Equalities x = t, which we turned into a model substitution x -> t, - * where x not in FV( t ), and - * (2) Equalities a*x*x + b*x + c = 0, which we turned into a model bound - * -b+s*sqrt(b*b-4*a*c)/2a - E <= x <= -b+s*sqrt(b*b-4*a*c)/2a + E. - * - * These literals are exempt from check-model, since they are satisfied by - * definition of our model construction. - */ - std::unordered_map<Node, Node, NodeHashFunction> d_check_model_solved; - /** has check model assignment - * - * Have we assigned v in the current checkModel(...) call? - * - * This method returns true if variable v is in the domain of - * d_check_model_bounds or if it occurs in d_check_model_vars. - */ - bool hasCheckModelAssignment(Node v) const; - /** have we successfully built the model in this SAT context? */ - context::CDO<bool> d_builtModel; //---------------------------end check model /** In the following functions, status states a relationship @@ -533,30 +401,35 @@ class NonlinearExtension { // per last-call effort // model values/orderings - /** cache of model values + + /** The non-linear model object * - * Stores the the concrete/abstract model values - * at indices 0 and 1 respectively. + * This class is responsible for computing model values for arithmetic terms + * and for establishing when we are able to answer "SAT". */ - std::map<Node, Node> d_mv[2]; + NlModel d_model; + /** have we successfully built the model in this SAT context? */ + context::CDO<bool> d_builtModel; // ordering, stores variables and 0,1,-1 std::map<Node, unsigned> d_order_vars; std::vector<Node> d_order_points; //transcendental functions + /** + * Maps arguments of SINE applications to a fresh skolem. This is used for + * ensuring that the argument of SINE we process are on the interval + * [-pi .. pi]. + */ std::map<Node, Node> d_tr_base; + /** Stores skolems in the range of the above map */ std::map<Node, bool> d_tr_is_base; std::map< Node, bool > d_tf_initial_refine; /** the list of lemmas we are waiting to flush until after check model */ std::vector<Node> d_waiting_lemmas; - /** did we use an approximation on this call to last-call effort? */ - bool d_used_approx; void mkPi(); void getCurrentPiBounds( std::vector< Node >& lemmas ); - /** print model value */ - void printModelValue(const char* c, Node n, unsigned prec = 5) const; private: //per last-call effort check |