diff options
9 files changed, 176 insertions, 93 deletions
diff --git a/src/theory/arith/nl/transcendental/exponential_solver.cpp b/src/theory/arith/nl/transcendental/exponential_solver.cpp index 7d7d0c23c..482e3bc21 100644 --- a/src/theory/arith/nl/transcendental/exponential_solver.cpp +++ b/src/theory/arith/nl/transcendental/exponential_solver.cpp @@ -170,31 +170,41 @@ void ExponentialSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx) d_data->d_im.addPendingArithLemma(lem, InferenceId::NL_T_TANGENT, nullptr, true); } -void ExponentialSolver::doSecantLemmas( - TNode e, TNode poly_approx, TNode c, TNode poly_approx_c, unsigned d) +void ExponentialSolver::doSecantLemmas(TNode e, + TNode poly_approx, + TNode center, + TNode cval, + unsigned d, + unsigned actual_d) { - d_data->doSecantLemmas( - getSecantBounds(e, c, d), poly_approx, c, poly_approx_c, e, d, 1); + d_data->doSecantLemmas(getSecantBounds(e, center, d), + poly_approx, + center, + cval, + e, + Convexity::CONVEX, + d, + actual_d); } std::pair<Node, Node> ExponentialSolver::getSecantBounds(TNode e, - TNode c, + TNode center, unsigned d) { - std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, c, d); + std::pair<Node, Node> bounds = d_data->getClosestSecantPoints(e, center, d); // Check if we already have neighboring secant points if (bounds.first.isNull()) { // pick c-1 bounds.first = Rewriter::rewrite( - NodeManager::currentNM()->mkNode(Kind::MINUS, c, d_data->d_one)); + NodeManager::currentNM()->mkNode(Kind::MINUS, center, d_data->d_one)); } if (bounds.second.isNull()) { // pick c+1 bounds.second = Rewriter::rewrite( - NodeManager::currentNM()->mkNode(Kind::PLUS, c, d_data->d_one)); + NodeManager::currentNM()->mkNode(Kind::PLUS, center, d_data->d_one)); } return bounds; } diff --git a/src/theory/arith/nl/transcendental/exponential_solver.h b/src/theory/arith/nl/transcendental/exponential_solver.h index b20c23e56..b4d08559a 100644 --- a/src/theory/arith/nl/transcendental/exponential_solver.h +++ b/src/theory/arith/nl/transcendental/exponential_solver.h @@ -87,12 +87,16 @@ class ExponentialSolver void doTangentLemma(TNode e, TNode c, TNode poly_approx); /** Sent secant lemmas around c for e */ - void doSecantLemmas( - TNode e, TNode poly_approx, TNode c, TNode poly_approx_c, unsigned d); + void doSecantLemmas(TNode e, + TNode poly_approx, + TNode center, + TNode cval, + unsigned d, + unsigned actual_d); private: /** Generate bounds for secant lemmas */ - std::pair<Node, Node> getSecantBounds(TNode e, TNode c, unsigned d); + std::pair<Node, Node> getSecantBounds(TNode e, TNode center, unsigned d); /** Holds common state for transcendental solvers */ TranscendentalState* d_data; diff --git a/src/theory/arith/nl/transcendental/sine_solver.cpp b/src/theory/arith/nl/transcendental/sine_solver.cpp index 31fd47503..cba85b80e 100644 --- a/src/theory/arith/nl/transcendental/sine_solver.cpp +++ b/src/theory/arith/nl/transcendental/sine_solver.cpp @@ -266,19 +266,20 @@ void SineSolver::doTangentLemma(TNode e, TNode c, TNode poly_approx, int region) // Figure 3: T( x ) // We use zero slope tangent planes, since the concavity of the Taylor // approximation cannot be easily established. - int concavity = regionToConcavity(region); + Convexity convexity = regionToConvexity(region); int mdir = regionToMonotonicityDir(region); + bool usec = (mdir == 1) == (convexity == Convexity::CONCAVE); Node lem = nm->mkNode( Kind::IMPLIES, nm->mkNode( Kind::AND, - nm->mkNode(Kind::GEQ, - e[0], - mdir == concavity ? Node(c) : regionToLowerBound(region)), - nm->mkNode(Kind::LEQ, - e[0], - mdir != concavity ? Node(c) : regionToUpperBound(region))), - nm->mkNode(concavity == 1 ? Kind::GEQ : Kind::LEQ, e, poly_approx)); + nm->mkNode( + Kind::GEQ, e[0], usec ? Node(c) : regionToLowerBound(region)), + nm->mkNode( + Kind::LEQ, e[0], usec ? Node(c) : regionToUpperBound(region))), + nm->mkNode(convexity == Convexity::CONVEX ? Kind::GEQ : Kind::LEQ, + e, + poly_approx)); Trace("nl-ext-sine") << "*** Tangent plane lemma (pre-rewrite): " << lem << std::endl; @@ -294,6 +295,7 @@ void SineSolver::doSecantLemmas(TNode e, TNode c, TNode poly_approx_c, unsigned d, + unsigned actual_d, int region) { d_data->doSecantLemmas(getSecantBounds(e, c, d, region), @@ -301,8 +303,9 @@ void SineSolver::doSecantLemmas(TNode e, c, poly_approx_c, e, + regionToConvexity(region), d, - regionToConcavity(region)); + actual_d); } std::pair<Node, Node> SineSolver::getSecantBounds(TNode e, diff --git a/src/theory/arith/nl/transcendental/sine_solver.h b/src/theory/arith/nl/transcendental/sine_solver.h index 9f9102a53..e41e6bd4f 100644 --- a/src/theory/arith/nl/transcendental/sine_solver.h +++ b/src/theory/arith/nl/transcendental/sine_solver.h @@ -93,6 +93,7 @@ class SineSolver TNode c, TNode poly_approx_c, unsigned d, + unsigned actual_d, int region); private: @@ -152,15 +153,15 @@ class SineSolver default: return 0; } } - int regionToConcavity(int region) + Convexity regionToConvexity(int region) { switch (region) { case 1: - case 2: return -1; + case 2: return Convexity::CONCAVE; case 3: - case 4: return 1; - default: return 0; + case 4: return Convexity::CONVEX; + default: return Convexity::UNKNOWN; } } diff --git a/src/theory/arith/nl/transcendental/taylor_generator.cpp b/src/theory/arith/nl/transcendental/taylor_generator.cpp index a373339e9..1b7257f8f 100644 --- a/src/theory/arith/nl/transcendental/taylor_generator.cpp +++ b/src/theory/arith/nl/transcendental/taylor_generator.cpp @@ -212,7 +212,7 @@ void TaylorGenerator::getPolynomialApproximationBounds( } } -void TaylorGenerator::getPolynomialApproximationBoundForArg( +unsigned TaylorGenerator::getPolynomialApproximationBoundForArg( Kind k, Node c, unsigned d, std::vector<Node>& pbounds) { getPolynomialApproximationBounds(k, d, pbounds); @@ -252,7 +252,9 @@ void TaylorGenerator::getPolynomialApproximationBoundForArg( getPolynomialApproximationBounds(k, ds, pboundss); pbounds[2] = pboundss[2]; } + return ds; } + return d; } std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlModel& model) diff --git a/src/theory/arith/nl/transcendental/taylor_generator.h b/src/theory/arith/nl/transcendental/taylor_generator.h index c647f7fd2..2dbfcccde 100644 --- a/src/theory/arith/nl/transcendental/taylor_generator.h +++ b/src/theory/arith/nl/transcendental/taylor_generator.h @@ -79,11 +79,13 @@ class TaylorGenerator * polynomials may depend on c. In particular, for P_u+[x] for <k>( c ) where * c>0, we return the P_u+[x] from the function above for the minimum degree * d' >= d such that (1-c^{2*d'+1}/(2*d'+1)!) is positive. + * @return the actual degree of the polynomial approximations (which may be + * larger than d). */ - void getPolynomialApproximationBoundForArg(Kind k, - Node c, - unsigned d, - std::vector<Node>& pbounds); + unsigned getPolynomialApproximationBoundForArg(Kind k, + Node c, + unsigned d, + std::vector<Node>& pbounds); /** get transcendental function model bounds * diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.cpp b/src/theory/arith/nl/transcendental/transcendental_solver.cpp index 8eb69e50b..ad3f49576 100644 --- a/src/theory/arith/nl/transcendental/transcendental_solver.cpp +++ b/src/theory/arith/nl/transcendental/transcendental_solver.cpp @@ -233,7 +233,8 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d) // mapped to for signs of c std::map<int, Node> poly_approx_bounds[2]; std::vector<Node> pbounds; - d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds); + unsigned actual_d = + d_tstate.d_taylor.getPolynomialApproximationBoundForArg(k, c, d, pbounds); poly_approx_bounds[0][1] = pbounds[0]; poly_approx_bounds[0][-1] = pbounds[1]; poly_approx_bounds[1][1] = pbounds[2]; @@ -362,11 +363,12 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d) { if (k == EXPONENTIAL) { - d_expSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d); + d_expSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d, actual_d); } else if (k == Kind::SINE) { - d_sineSlv.doSecantLemmas(tf, poly_approx, c, poly_approx_c, d, region); + d_sineSlv.doSecantLemmas( + tf, poly_approx, c, poly_approx_c, d, actual_d, region); } } return true; diff --git a/src/theory/arith/nl/transcendental/transcendental_state.cpp b/src/theory/arith/nl/transcendental/transcendental_state.cpp index 41ed2c53d..69678c621 100644 --- a/src/theory/arith/nl/transcendental/transcendental_state.cpp +++ b/src/theory/arith/nl/transcendental/transcendental_state.cpp @@ -267,26 +267,26 @@ void TranscendentalState::getCurrentPiBounds() } std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e, - TNode c, + TNode center, unsigned d) { // bounds are the minimum and maximum previous secant points // should not repeat secant points: secant lemmas should suffice to // rule out previous assignment - Assert( - std::find(d_secant_points[e][d].begin(), d_secant_points[e][d].end(), c) - == d_secant_points[e][d].end()); + Assert(std::find( + d_secant_points[e][d].begin(), d_secant_points[e][d].end(), center) + == d_secant_points[e][d].end()); // Insert into the (temporary) vector. We do not update this vector // until we are sure this secant plane lemma has been processed. We do // this by mapping the lemma to a side effect below. std::vector<Node> spoints = d_secant_points[e][d]; - spoints.push_back(c); + spoints.push_back(center); // sort sortByNlModel(spoints.begin(), spoints.end(), &d_model); // get the resulting index of c unsigned index = - std::find(spoints.begin(), spoints.end(), c) - spoints.begin(); + std::find(spoints.begin(), spoints.end(), center) - spoints.begin(); // bounds are the next closest upper/lower bound values return {index > 0 ? spoints[index - 1] : Node(), @@ -294,24 +294,40 @@ std::pair<Node, Node> TranscendentalState::getClosestSecantPoints(TNode e, } Node TranscendentalState::mkSecantPlane( - TNode arg, TNode b, TNode c, TNode approx_b, TNode approx_c) + TNode arg, TNode lower, TNode upper, TNode lval, TNode uval) { NodeManager* nm = NodeManager::currentNM(); // Figure 3: S_l( x ), S_u( x ) for s = 0,1 - Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, b, c)); + Node rcoeff_n = Rewriter::rewrite(nm->mkNode(Kind::MINUS, lower, upper)); Assert(rcoeff_n.isConst()); Rational rcoeff = rcoeff_n.getConst<Rational>(); Assert(rcoeff.sgn() != 0); - return nm->mkNode(Kind::PLUS, - approx_b, - nm->mkNode(Kind::MULT, - nm->mkNode(Kind::MINUS, approx_b, approx_c), - nm->mkConst(rcoeff.inverse()), - nm->mkNode(Kind::MINUS, arg, b))); + Node res = + nm->mkNode(Kind::PLUS, + lval, + nm->mkNode(Kind::MULT, + nm->mkNode(Kind::DIVISION, + nm->mkNode(Kind::MINUS, lval, uval), + nm->mkNode(Kind::MINUS, lower, upper)), + nm->mkNode(Kind::MINUS, arg, lower))); + Trace("nl-trans") << "Creating secant plane for transcendental function of " + << arg << std::endl; + Trace("nl-trans") << "\tfrom ( " << lower << " ; " << lval << " ) to ( " + << upper << " ; " << uval << " )" << std::endl; + Trace("nl-trans") << "\t" << res << std::endl; + Trace("nl-trans") << "\trewritten: " << Rewriter::rewrite(res) << std::endl; + return res; } -NlLemma TranscendentalState::mkSecantLemma( - TNode lower, TNode upper, int concavity, TNode tf, TNode splane) +NlLemma TranscendentalState::mkSecantLemma(TNode lower, + TNode upper, + TNode lapprox, + TNode uapprox, + int csign, + Convexity convexity, + TNode tf, + TNode splane, + unsigned actual_d) { NodeManager* nm = NodeManager::currentNM(); // With respect to Figure 3, this is slightly different. @@ -329,60 +345,73 @@ NlLemma TranscendentalState::mkSecantLemma( Node antec_n = nm->mkNode(Kind::AND, nm->mkNode(Kind::GEQ, tf[0], lower), nm->mkNode(Kind::LEQ, tf[0], upper)); + Trace("nl-trans") << "Bound for secant plane: " << lower << " <= " << tf[0] + << " <= " << upper << std::endl; + Trace("nl-trans") << "\t" << antec_n << std::endl; + // Convex: actual value is below the secant + // Concave: actual value is above the secant Node lem = nm->mkNode( Kind::IMPLIES, antec_n, - nm->mkNode(concavity == 1 ? Kind::LEQ : Kind::GEQ, tf, splane)); - Trace("nl-ext-tftp-debug2") - << "*** Secant plane lemma (pre-rewrite) : " << lem << std::endl; + nm->mkNode( + convexity == Convexity::CONVEX ? Kind::LEQ : Kind::GEQ, tf, splane)); + Trace("nl-trans-lemma") << "*** Secant plane lemma (pre-rewrite) : " << lem + << std::endl; lem = Rewriter::rewrite(lem); - Trace("nl-ext-tftp-lemma") << "*** Secant plane lemma : " << lem << std::endl; + Trace("nl-trans-lemma") << "*** Secant plane lemma : " << lem << std::endl; Assert(d_model.computeAbstractModelValue(lem) == d_false); return NlLemma(lem, LemmaProperty::NONE, nullptr, InferenceId::NL_T_SECANT); } void TranscendentalState::doSecantLemmas(const std::pair<Node, Node>& bounds, TNode poly_approx, - TNode c, - TNode approx_c, + TNode center, + TNode cval, TNode tf, + Convexity convexity, unsigned d, - int concavity) + unsigned actual_d) { - Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds.first - << " ... " << bounds.second << std::endl; - // take the model value of l or u (since may contain PI) - // Make secant from bounds.first to c - Node lval = d_model.computeAbstractModelValue(bounds.first); - Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds.first - << " is " << lval << std::endl; - if (lval != c) + int csign = center.getConst<Rational>().sgn(); + Trace("nl-trans") << "...do secant lemma with center " << center << " val " + << cval << " sign " << csign << std::endl; + Trace("nl-trans") << "...secant bounds are : " << bounds.first << " ... " + << bounds.second << std::endl; + // take the model value of lower (since may contain PI) + // Make secant from bounds.first to center + Node lower = d_model.computeAbstractModelValue(bounds.first); + Trace("nl-trans") << "...model value of bound " << bounds.first << " is " + << lower << std::endl; + if (lower != center) { // Figure 3 : P(l), P(u), for s = 0 - Node approx_l = Rewriter::rewrite( - poly_approx.substitute(d_taylor.getTaylorVariable(), lval)); - Node splane = mkSecantPlane(tf[0], lval, c, approx_l, approx_c); - NlLemma nlem = mkSecantLemma(lval, c, concavity, tf, splane); + Node lval = Rewriter::rewrite( + poly_approx.substitute(d_taylor.getTaylorVariable(), lower)); + Node splane = mkSecantPlane(tf[0], lower, center, lval, cval); + NlLemma nlem = mkSecantLemma( + lower, center, lval, cval, csign, convexity, tf, splane, actual_d); // The side effect says that if lem is added, then we should add the // secant point c for (tf,d). - nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c)); + nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center)); d_im.addPendingArithLemma(nlem, true); } - // Make secant from c to bounds.second - Node uval = d_model.computeAbstractModelValue(bounds.second); - Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds.second - << " is " << uval << std::endl; - if (c != uval) + // take the model value of upper (since may contain PI) + // Make secant from center to bounds.second + Node upper = d_model.computeAbstractModelValue(bounds.second); + Trace("nl-trans") << "...model value of bound " << bounds.second << " is " + << upper << std::endl; + if (center != upper) { // Figure 3 : P(l), P(u), for s = 1 - Node approx_u = Rewriter::rewrite( - poly_approx.substitute(d_taylor.getTaylorVariable(), uval)); - Node splane = mkSecantPlane(tf[0], c, uval, approx_c, approx_u); - NlLemma nlem = mkSecantLemma(c, uval, concavity, tf, splane); + Node uval = Rewriter::rewrite( + poly_approx.substitute(d_taylor.getTaylorVariable(), upper)); + Node splane = mkSecantPlane(tf[0], center, upper, cval, uval); + NlLemma nlem = mkSecantLemma( + center, upper, cval, uval, csign, convexity, tf, splane, actual_d); // The side effect says that if lem is added, then we should add the // secant point c for (tf,d). - nlem.d_secantPoint.push_back(std::make_tuple(tf, d, c)); + nlem.d_secantPoint.push_back(std::make_tuple(tf, d, center)); d_im.addPendingArithLemma(nlem, true); } } diff --git a/src/theory/arith/nl/transcendental/transcendental_state.h b/src/theory/arith/nl/transcendental/transcendental_state.h index 0a4e46eca..e1510c3b3 100644 --- a/src/theory/arith/nl/transcendental/transcendental_state.h +++ b/src/theory/arith/nl/transcendental/transcendental_state.h @@ -27,6 +27,24 @@ namespace nl { namespace transcendental { /** + * This enum indicates whether some function is convex, concave or unknown at + * some point. + */ +enum class Convexity +{ + CONVEX, + CONCAVE, + UNKNOWN +}; +inline std::ostream& operator<<(std::ostream& os, Convexity c) { + switch (c) { + case Convexity::CONVEX: return os << "CONVEX"; + case Convexity::CONCAVE: return os << "CONCAVE"; + default: return os << "UNKNOWN"; + } +} + +/** * Holds common state and utilities for transcendental solvers. * * This includes common lookups and caches as well as generic utilities for @@ -60,20 +78,24 @@ struct TranscendentalState * Get the two closest secant points from the once stored already. * "closest" is determined according to the current model. * @param e The transcendental term (like (exp t)) - * @param c The point currently under consideration (probably the model of t) + * @param center The point currently under consideration (probably the model + * of t) * @param d The taylor degree. */ - std::pair<Node, Node> getClosestSecantPoints(TNode e, TNode c, unsigned d); + std::pair<Node, Node> getClosestSecantPoints(TNode e, + TNode center, + unsigned d); /** - * Construct a secant plane between b and c + * Construct a secant plane as function in arg between lower and upper * @param arg The argument of the transcendental term - * @param b Left secant point - * @param c Right secant point - * @param approx Approximation for b (not yet substituted) - * @param approx_c Approximation for c (already substituted) + * @param lower Left secant point + * @param upper Right secant point + * @param lval Evaluation at lower + * @param uval Evaluation at upper */ - Node mkSecantPlane(TNode arg, TNode b, TNode c, TNode approx, TNode approx_c); + Node mkSecantPlane( + TNode arg, TNode lower, TNode upper, TNode lval, TNode uval); /** * Construct a secant lemma between lower and upper for tf. @@ -83,26 +105,34 @@ struct TranscendentalState * @param tf Current transcendental term * @param splane Secant plane as computed by mkSecantPlane() */ - NlLemma mkSecantLemma( - TNode lower, TNode upper, int concavity, TNode tf, TNode splane); + NlLemma mkSecantLemma(TNode lower, + TNode upper, + TNode lapprox, + TNode uapprox, + int csign, + Convexity convexity, + TNode tf, + TNode splane, + unsigned actual_d); /** * Construct and send secant lemmas (if appropriate) * @param bounds Secant bounds * @param poly_approx Polynomial approximation - * @param c Current point - * @param approx_c Approximation for c + * @param center Current point + * @param cval Evaluation at c * @param tf Current transcendental term * @param d Current taylor degree * @param concavity Concavity in region of c */ void doSecantLemmas(const std::pair<Node, Node>& bounds, TNode poly_approx, - TNode c, - TNode approx_c, + TNode center, + TNode cval, TNode tf, + Convexity convexity, unsigned d, - int concavity); + unsigned actual_d); Node d_true; Node d_false; |