From 047b8f69db1ab46ad68a2693565066f2a8d40b29 Mon Sep 17 00:00:00 2001 From: Andrew Reynolds Date: Wed, 22 Nov 2017 12:56:13 -0600 Subject: [PATCH] Transcendental tangent planes utilities (#1288) --- src/theory/arith/nonlinear_extension.cpp | 492 ++++++++++++++++----- src/theory/arith/nonlinear_extension.h | 525 +++++++++++++++++++++-- 2 files changed, 858 insertions(+), 159 deletions(-) diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp index 7993ba912..e3e078bbe 100644 --- a/src/theory/arith/nonlinear_extension.cpp +++ b/src/theory/arith/nonlinear_extension.cpp @@ -229,6 +229,12 @@ NonlinearExtension::NonlinearExtension(TheoryArith& containing, d_order_points.push_back(d_neg_one); d_order_points.push_back(d_zero); d_order_points.push_back(d_one); + d_taylor_real_fv = NodeManager::currentNM()->mkBoundVar( + "x", NodeManager::currentNM()->realType()); + d_taylor_real_fv_base = NodeManager::currentNM()->mkBoundVar( + "a", NodeManager::currentNM()->realType()); + d_taylor_real_fv_base_rem = NodeManager::currentNM()->mkBoundVar( + "b", NodeManager::currentNM()->realType()); } NonlinearExtension::~NonlinearExtension() {} @@ -663,31 +669,48 @@ bool NonlinearExtension::getCurrentSubstitution( std::pair NonlinearExtension::isExtfReduced( int effort, Node n, Node on, const std::vector& exp) const { if (n != d_zero) { - return std::make_pair(n.getKind() != kind::NONLINEAR_MULT, Node::null()); + Kind k = n.getKind(); + return std::make_pair(k != kind::NONLINEAR_MULT && !isTranscendentalKind(k), + Node::null()); } Assert(n == d_zero); - Trace("nl-ext-zero-exp") << "Infer zero : " << on << " == " << n << std::endl; - // minimize explanation - const std::set vars(on.begin(), on.end()); - - for (unsigned i = 0; i < exp.size(); i++) { - Trace("nl-ext-zero-exp") << " exp[" << i << "] = " << exp[i] << std::endl; - std::vector eqs; - if (exp[i].getKind() == kind::EQUAL) { - eqs.push_back(exp[i]); - } else if (exp[i].getKind() == kind::AND) { - for (unsigned j = 0; j < exp[i].getNumChildren(); j++) { - if (exp[i][j].getKind() == kind::EQUAL) { - eqs.push_back(exp[i][j]); + if (on.getKind() == kind::NONLINEAR_MULT) + { + Trace("nl-ext-zero-exp") << "Infer zero : " << on << " == " << n + << std::endl; + // minimize explanation if a substitution+rewrite results in zero + const std::set vars(on.begin(), on.end()); + + for (unsigned i = 0, size = exp.size(); i < size; i++) + { + Trace("nl-ext-zero-exp") << " exp[" << i << "] = " << exp[i] + << std::endl; + std::vector eqs; + if (exp[i].getKind() == kind::EQUAL) + { + eqs.push_back(exp[i]); + } + else if (exp[i].getKind() == kind::AND) + { + for (const Node& ec : exp[i]) + { + if (ec.getKind() == kind::EQUAL) + { + eqs.push_back(ec); + } } } - } - for (unsigned j = 0; j < eqs.size(); j++) { - for (unsigned r = 0; r < 2; r++) { - if (eqs[j][r] == d_zero && vars.find(eqs[j][1 - r]) != vars.end()) { - Trace("nl-ext-zero-exp") << "...single exp : " << eqs[j] << std::endl; - return std::make_pair(true, eqs[j]); + for (unsigned j = 0; j < eqs.size(); j++) + { + for (unsigned r = 0; r < 2; r++) + { + if (eqs[j][r] == d_zero && vars.find(eqs[j][1 - r]) != vars.end()) + { + Trace("nl-ext-zero-exp") << "...single exp : " << eqs[j] + << std::endl; + return std::make_pair(true, eqs[j]); + } } } } @@ -716,6 +739,8 @@ Node NonlinearExtension::computeModelValue(Node n, unsigned index) { //Assert( ret.isConst() ); } else if (n.getNumChildren() == 0) { if( n.getKind()==kind::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); @@ -1101,7 +1126,8 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, d_ci_exp.clear(); d_ci_max.clear(); d_tf_rep_map.clear(); - + d_tf_region.clear(); + int lemmas_proc = 0; std::vector lemmas; @@ -1165,6 +1191,8 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, getCurrentPiBounds( lemmas ); } Node shift = NodeManager::currentNM()->mkSkolem( "s", NodeManager::currentNM()->integerType(), "number of shifts" ); + // FIXME : do not introduce shift here, instead needs model-based + // refinement for constant shifts (#1284) Node shift_lem = NodeManager::currentNM()->mkNode( kind::AND, mkValidPhase( y, d_pi ), a[0].eqNode( NodeManager::currentNM()->mkNode( kind::PLUS, y, NodeManager::currentNM()->mkNode( kind::MULT, NodeManager::currentNM()->mkConst( Rational(2) ), shift, d_pi ) ) ), @@ -1184,7 +1212,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, if( itrm!=d_tf_rep_map[a.getKind()].end() ){ //verify they have the same model value if( d_mv[1][a]!=d_mv[1][itrm->second] ){ - //congruence lemma + // if not, add congruence lemma Node cong_lemma = NodeManager::currentNM()->mkNode( kind::IMPLIES, a[0].eqNode( itrm->second[0] ), a.eqNode( itrm->second ) ); lemmas.push_back( cong_lemma ); //Assert( false ); @@ -1608,6 +1636,10 @@ int NonlinearExtension::compare_value(Node i, Node j, } 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); @@ -2654,42 +2686,9 @@ std::vector NonlinearExtension::checkTranscendentalInitialRefine() { std::vector< Node > lemmas; Trace("nl-ext") << "Get initial refinement lemmas for transcendental functions..." << std::endl; for( std::map< Kind, std::map< Node, Node > >::iterator it = d_tf_rep_map.begin(); it != d_tf_rep_map.end(); ++it ){ - std::map< Node, Node > mv_to_term; for( std::map< Node, Node >::iterator itt = it->second.begin(); itt != it->second.end(); ++itt ){ Node t = itt->second; Assert( d_mv[1].find( t )!=d_mv[1].end() ); - Node tv = d_mv[1][t]; - mv_to_term[tv] = t; - // easy model-based bounds (TODO: make these unconditional?) - if( it->first==kind::SINE ){ - Assert( tv.isConst() ); - /* - if( tv.getConst() > d_one.getConst() ){ - lemmas.push_back( NodeManager::currentNM()->mkNode( kind::LEQ, t, d_one ) ); - }else if( tv.getConst() < d_neg_one.getConst() ){ - lemmas.push_back( NodeManager::currentNM()->mkNode( kind::GEQ, t, d_neg_one ) ); - } - */ - /* - if( tv.getConst().sgn()!=0 ){ - //symmetry (model-based) - Node neg_tv = NodeManager::currentNM()->mkConst( -tv.getConst() ); - if( mv_to_term.find( neg_tv )!=mv_to_term.end() ){ - Node sum = NodeManager::currentNM()->mkNode( kind::PLUS, mv_to_term[neg_tv][0], t[0] ); - Node res = computeModelValue( sum, 0 ); - if( !res.isConst() || res.getConst().sgn()!=0 ){ - Node tsum = NodeManager::currentNM()->mkNode( kind::PLUS, mv_to_term[neg_tv], t ); - Node sym_lem = NodeManager::currentNM()->mkNode( kind::IMPLIES, tsum.eqNode( d_zero ), sum.eqNode( d_zero ) ); - lemmas.push_back( sym_lem ); - } - } - } - */ - }else if( it->first==kind::EXPONENTIAL ){ - if( tv.getConst().sgn()==-1 ){ - lemmas.push_back( NodeManager::currentNM()->mkNode( kind::GT, t, d_zero ) ); - } - } //initial refinements if( d_tf_initial_refine.find( t )==d_tf_initial_refine.end() ){ d_tf_initial_refine[t] = true; @@ -2736,15 +2735,25 @@ std::vector NonlinearExtension::checkTranscendentalInitialRefine() { NodeManager::currentNM()->mkNode( kind::GT, t[0], d_pi_neg ), NodeManager::currentNM()->mkNode( kind::GT, t, NodeManager::currentNM()->mkNode( kind::MINUS, d_pi_neg, t[0] ) ) ) ) ); }else if( it->first==kind::EXPONENTIAL ){ - // exp(x)>0 ^ x < 0 <=> exp( x ) < 1 ^ ( x = 0 V exp( x ) > x + 1 ) - lem = NodeManager::currentNM()->mkNode( kind::AND, NodeManager::currentNM()->mkNode( kind::EQUAL, t[0].eqNode(d_zero), t.eqNode(d_one)), - NodeManager::currentNM()->mkNode( kind::EQUAL, - NodeManager::currentNM()->mkNode( kind::LT, t[0], d_zero ), - NodeManager::currentNM()->mkNode( kind::LT, t, d_one ) ), - NodeManager::currentNM()->mkNode( kind::OR, - NodeManager::currentNM()->mkNode( kind::LEQ, t[0], d_zero ), - NodeManager::currentNM()->mkNode( kind::GT, t, NodeManager::currentNM()->mkNode( kind::PLUS, t[0], d_one ) ) ) ); - + // ( exp(x) > 0 ) ^ ( x=0 <=> exp( x ) = 1 ) ^ ( x < 0 <=> exp( x ) < + // 1 ) ^ ( x <= 0 V exp( x ) > x + 1 ) + lem = NodeManager::currentNM()->mkNode( + kind::AND, + NodeManager::currentNM()->mkNode(kind::GT, t, d_zero), + NodeManager::currentNM()->mkNode( + kind::EQUAL, t[0].eqNode(d_zero), t.eqNode(d_one)), + NodeManager::currentNM()->mkNode( + kind::EQUAL, + NodeManager::currentNM()->mkNode(kind::LT, t[0], d_zero), + NodeManager::currentNM()->mkNode(kind::LT, t, d_one)), + NodeManager::currentNM()->mkNode( + kind::OR, + NodeManager::currentNM()->mkNode(kind::LEQ, t[0], d_zero), + NodeManager::currentNM()->mkNode( + kind::GT, + t, + NodeManager::currentNM()->mkNode( + kind::PLUS, t[0], d_one)))); } if( !lem.isNull() ){ lemmas.push_back( lem ); @@ -2794,20 +2803,15 @@ std::vector NonlinearExtension::checkTranscendentalMonotonic() { Assert( d_mv[1].find( t )!=d_mv[1].end() ); Trace("nl-ext-tf-mono") << " f-val : " << d_mv[1][t] << std::endl; } - std::vector< int > mdirs; std::vector< Node > mpoints; std::vector< Node > mpoints_vals; if( it->first==kind::SINE ){ - mdirs.push_back( -1 ); mpoints.push_back( d_pi ); - mdirs.push_back( 1 ); mpoints.push_back( d_pi_2 ); - mdirs.push_back( -1 ); + mpoints.push_back(d_zero); mpoints.push_back( d_pi_neg_2 ); - mdirs.push_back( 0 ); mpoints.push_back( d_pi_neg ); }else if( it->first==kind::EXPONENTIAL ){ - mdirs.push_back( 1 ); mpoints.push_back( Node::null() ); } if( !mpoints.empty() ){ @@ -2837,7 +2841,8 @@ std::vector NonlinearExtension::checkTranscendentalMonotonic() { //increment to the proper monotonicity region bool increment = true; - while( increment && mdir_index NonlinearExtension::checkTranscendentalMonotonic() { } if( increment ){ tval = Node::null(); - monotonic_dir = mdirs[mdir_index]; mono_bounds[1] = mpoints[mdir_index]; mdir_index++; - if( mdir_indexfirst, mdir_index); + if (mdir_index < mpoints.size()) + { mono_bounds[0] = mpoints[mdir_index]; }else{ mono_bounds[0] = Node::null(); } } } - - + // store the concavity region + d_tf_region[s] = mdir_index; + Trace("nl-ext-concavity") << "Transcendental function " << s + << " is in region #" << mdir_index; + Trace("nl-ext-concavity") << ", arg model value = " << sargval + << std::endl; + if( !tval.isNull() ){ Node mono_lem; if( monotonic_dir==1 && sval.getConst() > tval.getConst() ){ @@ -2887,6 +2898,7 @@ std::vector NonlinearExtension::checkTranscendentalMonotonic() { lemmas.push_back( mono_lem ); } } + // store the previous values targ = sarg; targval = sargval; t = s; @@ -2895,53 +2907,303 @@ std::vector NonlinearExtension::checkTranscendentalMonotonic() { } } } - - - - return lemmas; } - - -Node NonlinearExtension::getTaylor( Node tf, Node x, unsigned n, std::vector< Node >& lemmas ) { - Node i_exp_base_term = Rewriter::rewrite( NodeManager::currentNM()->mkNode( kind::MINUS, x, tf[0] ) ); - Node i_exp_base = getFactorSkolem( i_exp_base_term, lemmas ); - Node i_derv = tf; - Node i_fact = d_one; - Node i_exp = d_one; - int i_derv_status = 0; - unsigned counter = 0; - std::vector< Node > sum; - do { - counter++; - if( tf.getKind()==kind::EXPONENTIAL ){ - //unchanged - }else if( tf.getKind()==kind::SINE ){ - if( i_derv_status%2==1 ){ - Node arg = NodeManager::currentNM()->mkNode( kind::MINUS, - NodeManager::currentNM()->mkNode( kind::MULT, - NodeManager::currentNM()->mkConst( Rational(1)/Rational(2) ), - NodeManager::currentNM()->mkNullaryOperator( NodeManager::currentNM()->realType(), kind::PI ) ), - tf[0] ); - i_derv = NodeManager::currentNM()->mkNode( kind::SINE, arg ); + +int NonlinearExtension::regionToMonotonicityDir(Kind k, int region) +{ + if (k == kind::EXPONENTIAL) + { + if (region == 1) + { + return 1; + } + } + else if (k == kind::SINE) + { + if (region == 1 || region == 4) + { + return -1; + } + else if (region == 2 || region == 3) + { + return 1; + } + } + return 0; +} + +int NonlinearExtension::regionToConcavity(Kind k, int region) +{ + if (k == kind::EXPONENTIAL) + { + if (region == 1) + { + return 1; + } + } + else if (k == kind::SINE) + { + if (region == 1 || region == 2) + { + return -1; + } + else if (region == 3 || region == 4) + { + return 1; + } + } + return 0; +} + +Node NonlinearExtension::regionToLowerBound(Kind k, int region) +{ + if (k == kind::SINE) + { + if (region == 1) + { + return d_pi_2; + } + else if (region == 2) + { + return d_zero; + } + else if (region == 3) + { + return d_pi_neg_2; + } + else if (region == 4) + { + return d_pi_neg; + } + } + return Node::null(); +} + +Node NonlinearExtension::regionToUpperBound(Kind k, int region) +{ + if (k == kind::SINE) + { + if (region == 1) + { + return d_pi; + } + else if (region == 2) + { + return d_pi_2; + } + else if (region == 3) + { + return d_zero; + } + else if (region == 4) + { + return d_pi_neg_2; + } + } + return Node::null(); +} + +Node NonlinearExtension::getDerivative(Node n, Node x) +{ + Assert(x.isVar()); + // only handle the cases of the taylor expansion of d + if (n.getKind() == kind::EXPONENTIAL) + { + if (n[0] == x) + { + return n; + } + } + else if (n.getKind() == kind::SINE) + { + if (n[0] == x) + { + Node na = NodeManager::currentNM()->mkNode(kind::MINUS, d_pi_2, n[0]); + Node ret = NodeManager::currentNM()->mkNode(kind::SINE, na); + ret = Rewriter::rewrite(ret); + return ret; + } + } + else if (n.getKind() == kind::PLUS) + { + std::vector dchildren; + for (unsigned i = 0; i < n.getNumChildren(); i++) + { + // PLUS is flattened in rewriter, recursion depth is bounded by 1 + Node dc = getDerivative(n[i], x); + if (dc.isNull()) + { + return dc; }else{ - i_derv = tf; + dchildren.push_back(dc); } - if( i_derv_status>=2 ){ - i_derv = NodeManager::currentNM()->mkNode( kind::MINUS, d_zero, i_derv ); + } + return NodeManager::currentNM()->mkNode(kind::PLUS, dchildren); + } + else if (n.getKind() == kind::MULT) + { + Assert(n[0].isConst()); + Node dc = getDerivative(n[1], x); + if (!dc.isNull()) + { + return NodeManager::currentNM()->mkNode(kind::MULT, n[0], dc); + } + } + else if (n.getKind() == kind::NONLINEAR_MULT) + { + unsigned xcount = 0; + std::vector children; + unsigned xindex = 0; + for (unsigned i = 0, size = n.getNumChildren(); i < size; i++) + { + if (n[i] == x) + { + xcount++; + xindex = i; } - i_derv = Rewriter::rewrite( i_derv ); - i_derv_status = i_derv_status==3 ? 0 : i_derv_status+1; + children.push_back(n[i]); } - Node curr = NodeManager::currentNM()->mkNode( kind::MULT, - NodeManager::currentNM()->mkNode( kind::DIVISION, i_derv, i_fact ), i_exp ); - sum.push_back( curr ); - i_fact = Rewriter::rewrite( NodeManager::currentNM()->mkNode( kind::MULT, NodeManager::currentNM()->mkConst( Rational( counter ) ) ) ); - i_exp = Rewriter::rewrite( NodeManager::currentNM()->mkNode( kind::MULT, i_exp_base, i_exp ) ); - }while( countermkNode( kind::PLUS, sum ); + if (xcount == 0) + { + return d_zero; + } + else + { + children[xindex] = NodeManager::currentNM()->mkConst(Rational(xcount)); + } + return NodeManager::currentNM()->mkNode(kind::MULT, children); + } + else if (n.isVar()) + { + return n == x ? d_one : d_zero; + } + else if (n.isConst()) + { + return d_zero; + } + Trace("nl-ext-debug") << "No derivative computed for " << n; + Trace("nl-ext-debug") << " for d/d{" << x << "}" << std::endl; + return Node::null(); } - + +std::pair NonlinearExtension::getTaylor(TNode fa, unsigned n) +{ + Node fac; // what term we cache for fa + if (fa[0] == d_zero) + { + // optimization : simpler to compute (x-fa[0])^n if we are centered around + // 0. + fac = fa; + } + else + { + // otherwise we use a standard factor a in (x-a)^n + fac = NodeManager::currentNM()->mkNode(fa.getKind(), d_taylor_real_fv_base); + } + Node taylor_rem; + Node taylor_sum; + // check if we have already computed this Taylor series + std::unordered_map::iterator itt = d_taylor_sum[fac].find(n); + if (itt == d_taylor_sum[fac].end()) + { + Node i_exp_base; + if (fa[0] == d_zero) + { + i_exp_base = d_taylor_real_fv; + } + else + { + i_exp_base = Rewriter::rewrite(NodeManager::currentNM()->mkNode( + kind::MINUS, d_taylor_real_fv, d_taylor_real_fv_base)); + } + Node i_derv = fac; + Node i_fact = d_one; + Node i_exp = d_one; + int i_derv_status = 0; + unsigned counter = 0; + std::vector sum; + do + { + counter++; + if (fa.getKind() == kind::EXPONENTIAL) + { + // unchanged + } + else if (fa.getKind() == kind::SINE) + { + if (i_derv_status % 2 == 1) + { + Node arg = NodeManager::currentNM()->mkNode( + kind::PLUS, d_pi_2, d_taylor_real_fv_base); + i_derv = NodeManager::currentNM()->mkNode(kind::SINE, arg); + } + else + { + i_derv = fa; + } + if (i_derv_status >= 2) + { + i_derv = + NodeManager::currentNM()->mkNode(kind::MINUS, d_zero, i_derv); + } + i_derv = Rewriter::rewrite(i_derv); + i_derv_status = i_derv_status == 3 ? 0 : i_derv_status + 1; + } + if (counter == (n + 1)) + { + TNode x = d_taylor_real_fv_base; + i_derv = i_derv.substitute(x, d_taylor_real_fv_base_rem); + } + Node curr = NodeManager::currentNM()->mkNode( + kind::MULT, + NodeManager::currentNM()->mkNode(kind::DIVISION, i_derv, i_fact), + i_exp); + if (counter == (n + 1)) + { + taylor_rem = curr; + } + else + { + sum.push_back(curr); + i_fact = Rewriter::rewrite(NodeManager::currentNM()->mkNode( + kind::MULT, + NodeManager::currentNM()->mkConst(Rational(counter)), + i_fact)); + i_exp = Rewriter::rewrite( + NodeManager::currentNM()->mkNode(kind::MULT, i_exp_base, i_exp)); + } + } while (counter <= n); + taylor_sum = sum.size() == 1 + ? sum[0] + : NodeManager::currentNM()->mkNode(kind::PLUS, sum); + + if (fac[0] != d_taylor_real_fv_base) + { + TNode x = d_taylor_real_fv_base; + taylor_sum = taylor_sum.substitute(x, fac[0]); + } + + // cache + d_taylor_sum[fac][n] = taylor_sum; + d_taylor_rem[fac][n] = taylor_rem; + } + else + { + taylor_sum = itt->second; + Assert(d_taylor_rem[fac].find(n) != d_taylor_rem[fac].end()); + taylor_rem = d_taylor_rem[fac][n]; + } + + // must substitute for the argument if we were using a different lookup + if (fa[0] != fac[0]) + { + TNode x = d_taylor_real_fv_base; + taylor_sum = taylor_sum.substitute(x, fa[0]); + } + return std::pair(taylor_sum, taylor_rem); +} + } // namespace arith } // namespace theory } // namespace CVC4 diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h index dcaa91dc5..7c41fa096 100644 --- a/src/theory/arith/nonlinear_extension.h +++ b/src/theory/arith/nonlinear_extension.h @@ -23,6 +23,7 @@ #include #include #include +#include #include #include @@ -42,25 +43,109 @@ namespace arith { typedef std::map NodeMultiset; +// TODO : refactor/document this class (#1287) +/** Non-linear extension class + * + * This class implements model-based refinement schemes + * for non-linear arithmetic, described in: + * + * - "Invariant Checking of NRA Transition Systems + * via Incremental Reduction to LRA with EUF" by + * Cimatti et al., TACAS 2017. + * + * - Section 5 of "Desiging Theory Solvers with + * Extensions" by Reynolds et al., FroCoS 2017. + * + * - "Satisfiability Modulo Transcendental + * Functions via Incremental Linearization" by Cimatti + * et al., CADE 2017. + * + * It's main functionality is a check(...) method, + * which is called by TheoryArithPrivate either: + * (1) at full effort with no conflicts or lemmas emitted, + * or + * (2) at last call effort. + * In this method, this class calls d_out->lemma(...) + * for valid arithmetic theory lemmas, based on the current set of assertions, + * where d_out is the output channel of TheoryArith. + */ class NonlinearExtension { public: NonlinearExtension(TheoryArith& containing, eq::EqualityEngine* ee); ~NonlinearExtension(); + /** Get current substitution + * + * This function and the one below are + * used for context-dependent + * simplification, see Section 3.1 of + * "Designing Theory Solvers with Extensions" + * by Reynolds et al. FroCoS 2017. + * + * effort : an identifier indicating the stage where + * we are performing context-dependent simplification, + * vars : a set of arithmetic variables. + * + * This function populates subs and exp, such that for 0 <= i < vars.size(): + * ( exp[vars[i]] ) => vars[i] = subs[i] + * where exp[vars[i]] is a set of assertions + * that hold in the current context. We call { vars -> subs } a "derivable + * substituion" (see Reynolds et al. FroCoS 2017). + */ bool getCurrentSubstitution(int effort, const std::vector& vars, std::vector& subs, std::map >& exp); - + /** Is the term n in reduced form? + * + * Used for context-dependent simplification. + * + * effort : an identifier indicating the stage where + * we are performing context-dependent simplification, + * on : the original term that we reduced to n, + * exp : an explanation such that ( exp => on = n ). + * + * We return a pair ( b, exp' ) such that + * if b is true, then: + * n is in reduced form + * if exp' is non-null, then ( exp' => on = n ) + * The second part of the pair is used for constructing + * minimal explanations for context-dependent simplifications. + */ std::pair isExtfReduced(int effort, Node n, Node on, const std::vector& exp) const; + /** Check at effort level e. + * + * This call may result in (possibly multiple) + * calls to d_out->lemma(...) where d_out + * is the output channel of TheoryArith. + */ void check(Theory::Effort e); + /** Does this class need a call to check(...) at last call effort? */ bool needsCheckLastEffort() const { return d_needsLastCall; } + /** 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 + * containing the factors of monomial b. + */ bool isMonomialSubset(Node a, Node b) const; void registerMonomialSubset(Node a, Node b); - private: typedef context::CDHashSet NodeSet; // monomial information (context-independent) @@ -83,13 +168,20 @@ class NonlinearExtension { Kind d_type; }; /* struct ConstraintInfo */ - // Check assertions for consistency in the effort LAST_CALL with a subset of - // the assertions, false_asserts, evaluate to false in the current model. - // Returns the number of lemmas added on the output channel. + /** check last call + * + * Check assertions for consistency in the effort LAST_CALL with a subset of + * the assertions, false_asserts, that evaluate to false in the current model. + * + * xts : the list of (non-reduced) extended terms in the current context. + * + * This method returns the number of lemmas added on the output channel of + * TheoryArith. + */ int checkLastCall(const std::vector& assertions, const std::set& false_asserts, const std::vector& xts); - + //---------------------------------------term utilities static bool isArithKind(Kind k); static Node mkLit(Node a, Node b, int status, int orderType = 0); static Node mkAbs(Node a); @@ -99,53 +191,141 @@ class NonlinearExtension { static Kind transKinds(Kind k1, Kind k2); static bool isTranscendentalKind(Kind k); Node mkMonomialRemFactor(Node n, const NodeMultiset& n_exp_rem) const; + //---------------------------------------end term utilities - // register monomial + /** register monomial */ void registerMonomial(Node n); void setMonomialFactor(Node a, Node b, const NodeMultiset& common); void registerConstraint(Node atom); - // index = 0 : concrete, 1 : abstract + /** 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& vars, NodeMultiset& d_order, unsigned orderType); - // Returns the subset of assertions that evaluate to false in the model. + /** Returns the subset of assertions whose concrete values are + * false in the model. + */ std::set getFalseInModel(const std::vector& assertions); - // status - // 0 : equal - // 1 : greater than or equal - // 2 : greater than - // -X : (less) - // in these functions we are iterating over variables of monomials - // initially : exp => ( oa = a ^ ob = b ) + /** In the following functions, status states a relationship + * between two arithmetic terms, where: + * 0 : equal + * 1 : greater than or equal + * 2 : greater than + * -X : (greater -> less) + * TODO (#1287) make this an enum? + */ + /** compute the sign of a. + * + * Calls to this function are such that : + * exp => ( oa = a ^ a 0 ) + * + * This function iterates over the factors of a, + * where a_index is the index of the factor in a + * we are currently looking at. + * + * This function returns a status, which indicates + * a's relationship to 0. + * We add lemmas to lem of the form given by the + * lemma schema checkSign(...). + */ int compareSign(Node oa, Node a, unsigned a_index, int status, std::vector& exp, std::vector& lem); + /** compare monomials a and b + * + * Initially, a call to this function is such that : + * exp => ( oa = a ^ ob = b ) + * + * This function returns true if we can infer a valid + * arithmetic lemma of the form : + * P => abs( a ) >= abs( b ) + * where P is true and abs( a ) >= abs( b ) is false in the + * current model. + * + * This function is implemented by "processing" factors + * of monomials a and b until an inference of the above + * form can be made. For example, if : + * a = x*x*y and b = z*w + * Assuming we are trying to show abs( a ) >= abs( c ), + * then if abs( M( x ) ) >= abs( M( z ) ) where M is the current model, + * then we can add abs( x ) >= abs( z ) to our explanation, and + * mark one factor of x as processed in a, and + * one factor of z as processed in b. The number of processed factors of a + * and b are stored in a_exp_proc and b_exp_proc respectively. + * + * cmp_infers stores information that is helpful + * in discarding redundant inferences. For example, + * we do not want to infer abs( x ) >= abs( z ) if + * we have already inferred abs( x ) >= abs( y ) and + * abs( y ) >= abs( z ). + * It stores entries of the form (status,t1,t2)->F, + * which indicates that we constructed a lemma F that + * showed t1 t2. + * + * We add lemmas to lem of the form given by the + * lemma schema checkMagnitude(...). + */ bool compareMonomial( Node oa, Node a, NodeMultiset& a_exp_proc, Node ob, Node b, NodeMultiset& b_exp_proc, std::vector& exp, std::vector& lem, std::map > >& cmp_infers); + /** helper function for above + * + * The difference is the inputs a_index and b_index, which are the indices of + * children (factors) in monomials a and b which we are currently looking at. + */ bool compareMonomial( Node oa, Node a, unsigned a_index, NodeMultiset& a_exp_proc, Node ob, Node b, unsigned b_index, NodeMultiset& b_exp_proc, int status, std::vector& exp, std::vector& lem, std::map > >& cmp_infers); + /** Check whether we have already inferred a relationship between monomials + * x and y based on the information in cmp_infers. This computes the + * transitive closure of the relation stored in cmp_infers. + */ bool cmp_holds(Node x, Node y, std::map >& cmp_infers, std::vector& exp, std::map& visited); - + /** Is n entailed with polarity pol in the current context? */ bool isEntailed(Node n, bool pol); - // Potentially sends lem on the output channel if lem has not been sent on the - // output channel in this context. Returns the number of lemmas sent on the - // output channel. + /** flush lemmas + * + * Potentially sends lem on the output channel if lem has not been sent on the + * output channel in this context. Returns the number of lemmas sent on the + * output channel of TheoryArith. + */ int flushLemma(Node lem); - // Potentially sends lemmas to the output channel and clears lemmas. Returns - // the number of lemmas sent to the output channel. + /** Potentially sends lemmas to the output channel and clears lemmas. Returns + * the number of lemmas sent to the output channel. + */ int flushLemmas(std::vector& lemmas); // Returns the NodeMultiset for an existing monomial. @@ -176,12 +356,27 @@ class NonlinearExtension { // literals with Skolems (need not be satisfied by model) NodeSet d_skolem_atoms; - // utilities + /** commonly used terms */ Node d_zero; Node d_one; Node d_neg_one; Node d_true; Node d_false; + /** PI + * + * Note that PI is a (symbolic, non-constant) nullary operator. This is + * because its value cannot be computed exactly. We constraint PI to concrete + * lower and upper bounds stored in d_pi_bound below. + */ + Node d_pi; + /** PI/2 */ + Node d_pi_2; + /** -PI/2 */ + Node d_pi_neg_2; + /** -PI */ + Node d_pi_neg; + /** the concrete lower and upper bounds for PI */ + Node d_pi_bound[2]; // The theory of arithmetic containing this extension. TheoryArith& d_containing; @@ -197,7 +392,11 @@ class NonlinearExtension { std::vector d_constraints; // model values/orderings - // model values + /** cache of model values + * + * Stores the the concrete/abstract model values + * at indices 0 and 1 respectively. + */ std::map d_mv[2]; // ordering, stores variables and 0,1,-1 @@ -208,11 +407,6 @@ class NonlinearExtension { std::map d_trig_base; std::map d_trig_is_base; std::map< Node, bool > d_tf_initial_refine; - Node d_pi; - Node d_pi_2; - Node d_pi_neg_2; - Node d_pi_neg; - Node d_pi_bound[2]; void mkPi(); void getCurrentPiBounds( std::vector< Node >& lemmas ); @@ -236,8 +430,14 @@ private: std::map > > d_ci; std::map > > d_ci_exp; std::map > > d_ci_max; - - //information about transcendental functions + + /** transcendental function representative map + * + * For each transcendental function n = tf( x ), + * this stores ( n.getKind(), r ) -> n + * where r is the current representative of x + * in the equality engine assoiated with this class. + */ std::map< Kind, std::map< Node, Node > > d_tf_rep_map; // factor skolems @@ -246,29 +446,266 @@ private: // tangent plane bounds std::map< Node, std::map< Node, Node > > d_tangent_val_bound[4]; - - Node getTaylor( Node tf, Node x, unsigned n, std::vector< Node >& lemmas ); -private: - // Returns a vector containing a split on whether each term is 0 or not for - // those terms that have not been split on in the current context. + + /** secant points (sorted list) for transcendental functions + * + * This is used for tangent plane refinements for + * transcendental functions. This is the set + * "get-previous-secant-points" in "Satisfiability + * Modulo Transcendental Functions via Incremental + * Linearization" by Cimatti et al., CADE 2017, for + * each transcendental function application. + */ + std::unordered_map, NodeHashFunction> d_secant_points; + + /** get Taylor series of degree n for function fa centered around point fa[0]. + * + * Return value is ( P_{n,f(a)}( x ), R_{n+1,f(a)}( x ) ) where + * the first part of the pair is the Taylor series expansion : + * P_{n,f(a)}( x ) = sum_{i=0}^n (f^i( a )/i!)*(x-a)^i + * and the second part of the pair is the Taylor series remainder : + * R_{n+1,f(a),b}( x ) = (f^{n+1}( b )/(n+1)!)*(x-a)^{n+1} + * + * The above values are cached for each (f,n) for a fixed variable "a". + * To compute the Taylor series for fa, we compute the Taylor series + * for ( fa.getKind(), n ) then substitute { a -> fa[0] } if fa[0]!=0. + * We compute P_{n,f(0)}( x )/R_{n+1,f(0),b}( x ) for ( fa.getKind(), n ) + * if fa[0]=0. + * In the latter case, note we compute the exponential x^{n+1} + * instead of (x-a)^{n+1}, which can be done faster. + */ + std::pair getTaylor(TNode fa, unsigned n); + + /** internal variables used for constructing (cached) versions of the Taylor + * series above. + */ + Node d_taylor_real_fv; // x above + Node d_taylor_real_fv_base; // a above + Node d_taylor_real_fv_base_rem; // b above + + /** cache of sum and remainder terms for getTaylor */ + std::unordered_map, NodeHashFunction> + d_taylor_sum; + std::unordered_map, NodeHashFunction> + d_taylor_rem; + + /** concavity region for transcendental functions + * + * This stores an integer that identifies an interval in + * which the current model value for an argument of an + * application of a transcendental function resides. + * + * For exp( x ): + * region #1 is -infty < x < infty + * For sin( x ): + * region #0 is pi < x < infty (this is an invalid region) + * region #1 is pi/2 < x <= pi + * region #2 is 0 < x <= pi/2 + * region #3 is -pi/2 < x <= 0 + * region #4 is -pi < x <= -pi/2 + * region #5 is -infty < x <= -pi (this is an invalid region) + * All regions not listed above, as well as regions 0 and 5 + * for SINE are "invalid". We only process applications + * of transcendental functions whose arguments have model + * values that reside in valid regions. + */ + std::unordered_map d_tf_region; + /** get monotonicity direction + * + * Returns whether the slope is positive (+1) or negative(-1) + * in region of transcendental function with kind k. + * Returns 0 if region is invalid. + */ + int regionToMonotonicityDir(Kind k, int region); + /** get concavity + * + * Returns whether we are concave (+1) or convex (-1) + * in region of transcendental function with kind k, + * where region is defined above. + * Returns 0 if region is invalid. + */ + int regionToConcavity(Kind k, int region); + /** region to lower bound + * + * Returns the term corresponding to the lower + * bound of the region of transcendental function + * with kind k. Returns Node::null if the region + * is invalid, or there is no lower bound for the + * region. + */ + Node regionToLowerBound(Kind k, int region); + /** region to upper bound + * + * Returns the term corresponding to the upper + * bound of the region of transcendental function + * with kind k. Returns Node::null if the region + * is invalid, or there is no upper bound for the + * region. + */ + Node regionToUpperBound(Kind k, int region); + /** get derivative + * + * Returns d/dx n. Supports cases of n + * for transcendental functions applied to x, + * multiplication, addition, constants and variables. + * Returns Node::null() if derivative is an + * unhandled case. + */ + Node getDerivative(Node n, Node x); + + private: + //-------------------------------------------- lemma schemas + /** check split zero + * + * Returns a set of theory lemmas of the form + * t = 0 V t != 0 + * where t is a term that exists in the current context. + */ std::vector checkSplitZero(); - + + /** check monomial sign + * + * Returns a set of valid theory lemmas, based on a + * lemma schema which ensures that non-linear monomials + * respect sign information based on their facts. + * For more details, see Section 5 of "Design Theory + * Solvers with Extensions" by Reynolds et al., FroCoS 2017, + * Figure 5, this is the schema "Sign". + * + * Examples: + * + * x > 0 ^ y > 0 => x*y > 0 + * x < 0 => x*y*y < 0 + * x = 0 => x*y*z = 0 + */ std::vector checkMonomialSign(); - + + /** check monomial magnitude + * + * Returns a set of valid theory lemmas, based on a + * lemma schema which ensures that comparisons between + * non-linear monomials respect the magnitude of their + * factors. + * For more details, see Section 5 of "Design Theory + * Solvers with Extensions" by Reynolds et al., FroCoS 2017, + * Figure 5, this is the schema "Magnitude". + * + * Examples: + * + * |x|>|y| => |x*z|>|y*z| + * |x|>|y| ^ |z|>|w| ^ |x|>=1 => |x*x*z*u|>|y*w| + */ std::vector checkMonomialMagnitude( unsigned c ); - + + /** check monomial inferred bounds + * + * Returns a set of valid theory lemmas, based on a + * lemma schema that infers new constraints about existing + * terms based on mulitplying both sides of an existing + * constraint by a term. + * For more details, see Section 5 of "Design Theory + * Solvers with Extensions" by Reynolds et al., FroCoS 2017, + * Figure 5, this is the schema "Multiply". + * + * Examples: + * + * x > 0 ^ (y > z + w) => x*y > x*(z+w) + * x < 0 ^ (y > z + w) => x*y < x*(z+w) + * ...where (y > z + w) and x*y are a constraint and term + * that occur in the current context. + */ std::vector checkMonomialInferBounds( std::vector& nt_lemmas, const std::set& false_asserts ); - + + /** check factoring + * + * Returns a set of valid theory lemmas, based on a + * lemma schema that states a relationship betwen monomials + * with common factors that occur in the same constraint. + * + * Examples: + * + * x*z+y*z > t => ( k = x + y ^ k*z > t ) + * ...where k is fresh and x*z + y*z > t is a + * constraint that occurs in the current context. + */ std::vector checkFactoring( const std::set& false_asserts ); - + + /** check monomial infer resolution bounds + * + * Returns a set of valid theory lemmas, based on a + * lemma schema which "resolves" upper bounds + * of one inequality with lower bounds for another. + * This schema is not enabled by default, and can + * be enabled by --nl-ext-rbound. + * + * Examples: + * + * ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t + * ...where s <= x*z and x*y <= t are constraints + * that occur in the current context. + */ std::vector checkMonomialInferResBounds(); - + + /** check tangent planes + * + * Returns a set of valid theory lemmas, based on an + * "incremental linearization" of non-linear monomials. + * This linearization is accomplished by adding constraints + * corresponding to "tangent planes" at the current + * model value of each non-linear monomial. In particular + * consider the definition for constants a,b : + * T_{a,b}( x*y ) = b*x + a*y - a*b. + * The lemmas added by this function are of the form : + * ( ( x>a ^ yb) ) => x*y < T_{a,b}( x*y ) + * ( ( x>a ^ y>b) ^ (x x*y > T_{a,b}( x*y ) + * It is inspired by "Invariant Checking of NRA Transition + * Systems via Incremental Reduction to LRA with EUF" by + * Cimatti et al., TACAS 2017. + * This schema is not terminating in general. + * It is not enabled by default, and can + * be enabled by --nl-ext-tplanes. + * + * Examples: + * + * ( ( x>2 ^ y>5) ^ (x<2 ^ y<5) ) => x*y > 5*x + 2*y - 10 + * ( ( x>2 ^ y<5) ^ (x<2 ^ y>5) ) => x*y < 5*x + 2*y - 10 + */ std::vector checkTangentPlanes(); - + + /** check transcendental initial refine + * + * Returns a set of valid theory lemmas, based on + * simple facts about transcendental functions. + * This mostly follows the initial axioms described in + * Section 4 of "Satisfiability + * Modulo Transcendental Functions via Incremental + * Linearization" by Cimatti et al., CADE 2017. + * + * Examples: + * + * sin( x ) = -sin( -x ) + * ( PI > x > 0 ) => 0 < sin( x ) < 1 + * exp( x )>0 + * x<0 => exp( x )<1 + */ std::vector checkTranscendentalInitialRefine(); + /** check transcendental monotonic + * + * Returns a set of valid theory lemmas, based on a + * lemma scheme that ensures that applications + * of transcendental functions respect monotonicity. + * + * Examples: + * + * x > y => exp( x ) > exp( y ) + * PI/2 > x > y > 0 => sin( x ) > sin( y ) + * PI > x > y > PI/2 => sin( x ) < sin( y ) + */ std::vector checkTranscendentalMonotonic(); + + //-------------------------------------------- end lemma schemas }; /* class NonlinearExtension */ } // namespace arith -- 2.30.2