From 39e82b212bd2f957805884c24e6414289a7ec987 Mon Sep 17 00:00:00 2001 From: Gereon Kremer Date: Mon, 22 Feb 2021 16:07:40 +0100 Subject: Cleanup in transcendental solver, add ApproximationBounds struct. (#5945) This PR merges some cleanup in the transcendental solver from proof-new. It adds a new struct ApproximationBounds that replaces an opaque std::vector and does some general refactoring in the TaylorGenerator class, removing dead code and using fixed-width integers. --- .../arith/nl/transcendental/taylor_generator.cpp | 242 +++++++-------------- .../arith/nl/transcendental/taylor_generator.h | 61 +++--- .../nl/transcendental/transcendental_solver.cpp | 12 +- 3 files changed, 121 insertions(+), 194 deletions(-) (limited to 'src/theory') diff --git a/src/theory/arith/nl/transcendental/taylor_generator.cpp b/src/theory/arith/nl/transcendental/taylor_generator.cpp index 745cba17a..375f1757e 100644 --- a/src/theory/arith/nl/transcendental/taylor_generator.cpp +++ b/src/theory/arith/nl/transcendental/taylor_generator.cpp @@ -25,214 +25,138 @@ namespace transcendental { TaylorGenerator::TaylorGenerator() : d_nm(NodeManager::currentNM()), - d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType())), - d_taylor_real_fv_base(d_nm->mkBoundVar("a", d_nm->realType())), - d_taylor_real_fv_base_rem(d_nm->mkBoundVar("b", d_nm->realType())) + d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType())) { } TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; } -std::pair TaylorGenerator::getTaylor(TNode fa, std::uint64_t n) +std::pair TaylorGenerator::getTaylor(Kind k, std::uint64_t n) { - NodeManager* nm = NodeManager::currentNM(); - Node d_zero = nm->mkConst(Rational(0)); - Assert(n > 0); - 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 + // check if we have already computed this Taylor series + auto itt = d_taylor_terms[k].find(n); + if (itt != d_taylor_terms[k].end()) { - // otherwise we use a standard factor a in (x-a)^n - fac = nm->mkNode(fa.getKind(), d_taylor_real_fv_base); + return itt->second; } - Node taylor_rem; - Node taylor_sum; - // check if we have already computed this Taylor series - auto itt = s_taylor_sum[fac].find(n); - if (itt == s_taylor_sum[fac].end()) + + NodeManager* nm = NodeManager::currentNM(); + + // the current factorial `counter!` + Integer factorial = 1; + // the current variable power `x^counter` + Node varpow = nm->mkConst(Rational(1)); + std::vector sum; + for (std::uint64_t counter = 1; counter <= n; ++counter) { - Node i_exp_base; - if (fa[0] == d_zero) + if (k == Kind::EXPONENTIAL) { - i_exp_base = d_taylor_real_fv; + // Maclaurin series for exponential: + // \sum_{n=0}^\infty x^n / n! + sum.push_back( + nm->mkNode(Kind::DIVISION, varpow, nm->mkConst(factorial))); } - else + else if (k == Kind::SINE) { - i_exp_base = Rewriter::rewrite( - nm->mkNode(Kind::MINUS, d_taylor_real_fv, d_taylor_real_fv_base)); - } - Node i_derv = fac; - Node i_fact = nm->mkConst(Rational(1)); - Node i_exp = i_fact; - 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) + // Maclaurin series for exponential: + // \sum_{n=0}^\infty (-1)^n / (2n+1)! * x^(2n+1) + if (counter % 2 == 0) { - if (i_derv_status % 2 == 1) - { - Node pi = nm->mkNullaryOperator(nm->realType(), Kind::PI); - Node pi_2 = Rewriter::rewrite(nm->mkNode( - Kind::MULT, pi, nm->mkConst(Rational(1) / Rational(2)))); - - Node arg = nm->mkNode(Kind::PLUS, pi_2, d_taylor_real_fv_base); - i_derv = nm->mkNode(Kind::SINE, arg); - } - else - { - i_derv = fa; - } - if (i_derv_status >= 2) - { - i_derv = nm->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; + int sign = (counter % 4 == 0 ? -1 : 1); + sum.push_back(nm->mkNode(Kind::MULT, + nm->mkNode(Kind::DIVISION, + nm->mkConst(sign), + nm->mkConst(factorial)), + varpow)); } - 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 = nm->mkNode( - Kind::MULT, nm->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( - nm->mkNode(Kind::MULT, nm->mkConst(Rational(counter)), i_fact)); - i_exp = Rewriter::rewrite(nm->mkNode(Kind::MULT, i_exp_base, i_exp)); - } - } while (counter <= n); - taylor_sum = sum.size() == 1 ? sum[0] : nm->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 - s_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]; + factorial *= counter; + varpow = + Rewriter::rewrite(nm->mkNode(Kind::MULT, d_taylor_real_fv, varpow)); } + Node taylor_sum = + Rewriter::rewrite(sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum)); + Node taylor_rem = Rewriter::rewrite( + nm->mkNode(Kind::DIVISION, varpow, nm->mkConst(factorial))); - // 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); + auto res = std::make_pair(taylor_sum, taylor_rem); + + // put result in cache + d_taylor_terms[k][n] = res; + + return res; } void TaylorGenerator::getPolynomialApproximationBounds( - Kind k, unsigned d, std::vector& pbounds) + Kind k, std::uint64_t d, ApproximationBounds& pbounds) { - if (d_poly_bounds[k][d].empty()) + auto it = d_poly_bounds[k].find(d); + if (it == d_poly_bounds[k].end()) { NodeManager* nm = NodeManager::currentNM(); - Node tft = nm->mkNode(k, nm->mkConst(Rational(0))); // n is the Taylor degree we are currently considering - unsigned n = 2 * d; + std::uint64_t n = 2 * d; // n must be even - std::pair taylor = getTaylor(tft, n); - Trace("nl-ext-tftp-debug2") - << "Taylor for " << k << " is : " << taylor.first << std::endl; - Node taylor_sum = Rewriter::rewrite(taylor.first); - Trace("nl-ext-tftp-debug2") - << "Taylor for " << k << " is (post-rewrite) : " << taylor_sum - << std::endl; - Assert(taylor.second.getKind() == Kind::MULT); - Assert(taylor.second.getNumChildren() == 2); - Assert(taylor.second[0].getKind() == Kind::DIVISION); - Trace("nl-ext-tftp-debug2") - << "Taylor remainder for " << k << " is " << taylor.second << std::endl; + std::pair taylor = getTaylor(k, n); + Node taylor_sum = taylor.first; // ru is x^{n+1}/(n+1)! - Node ru = nm->mkNode(Kind::DIVISION, taylor.second[1], taylor.second[0][1]); - ru = Rewriter::rewrite(ru); - Trace("nl-ext-tftp-debug2") - << "Taylor remainder factor is (post-rewrite) : " << ru << std::endl; + Node ru = taylor.second; + Trace("nl-trans") << "Taylor for " << k << " is : " << taylor.first + << std::endl; + Trace("nl-trans") << "Taylor remainder for " << k << " is " << taylor.second + << std::endl; if (k == Kind::EXPONENTIAL) { - pbounds.push_back(taylor_sum); - pbounds.push_back(taylor_sum); - pbounds.push_back(Rewriter::rewrite( + pbounds.d_lower = taylor_sum; + pbounds.d_upperNeg = + Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru)); + pbounds.d_upperPos = Rewriter::rewrite( nm->mkNode(Kind::MULT, taylor_sum, - nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru)))); - pbounds.push_back( - Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru))); + nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru))); } else { Assert(k == Kind::SINE); Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru)); Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru)); - pbounds.push_back(l); - pbounds.push_back(l); - pbounds.push_back(u); - pbounds.push_back(u); + pbounds.d_lower = l; + pbounds.d_upperNeg = u; + pbounds.d_upperPos = u; } - Trace("nl-ext-tf-tplanes") - << "Polynomial approximation for " << k << " is: " << std::endl; - Trace("nl-ext-tf-tplanes") << " Lower (pos): " << pbounds[0] << std::endl; - Trace("nl-ext-tf-tplanes") << " Upper (pos): " << pbounds[2] << std::endl; - Trace("nl-ext-tf-tplanes") << " Lower (neg): " << pbounds[1] << std::endl; - Trace("nl-ext-tf-tplanes") << " Upper (neg): " << pbounds[3] << std::endl; - d_poly_bounds[k][d].insert( - d_poly_bounds[k][d].end(), pbounds.begin(), pbounds.end()); + Trace("nl-trans") << "Polynomial approximation for " << k + << " is: " << std::endl; + Trace("nl-trans") << " Lower: " << pbounds.d_lower << std::endl; + Trace("nl-trans") << " Upper (neg): " << pbounds.d_upperNeg << std::endl; + Trace("nl-trans") << " Upper (pos): " << pbounds.d_upperPos << std::endl; + d_poly_bounds[k].emplace(d, pbounds); } else { - pbounds.insert( - pbounds.end(), d_poly_bounds[k][d].begin(), d_poly_bounds[k][d].end()); + pbounds = it->second; } } -unsigned TaylorGenerator::getPolynomialApproximationBoundForArg( - Kind k, Node c, unsigned d, std::vector& pbounds) +std::uint64_t TaylorGenerator::getPolynomialApproximationBoundForArg( + Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds) { getPolynomialApproximationBounds(k, d, pbounds); + Trace("nl-trans") << "c = " << c << std::endl; Assert(c.isConst()); if (k == Kind::EXPONENTIAL && c.getConst().sgn() == 1) { - NodeManager* nm = NodeManager::currentNM(); - Node tft = nm->mkNode(k, nm->mkConst(Rational(0))); bool success = false; - unsigned ds = d; + std::uint64_t ds = d; TNode ttrf = getTaylorVariable(); TNode tc = c; do { success = true; - unsigned n = 2 * ds; - std::pair taylor = getTaylor(tft, n); + std::uint64_t n = 2 * ds; + std::pair taylor = getTaylor(k, n); // check that 1-c^{n+1}/(n+1)! > 0 - Node ru = - nm->mkNode(Kind::DIVISION, taylor.second[1], taylor.second[0][1]); + Node ru = taylor.second; Node rus = ru.substitute(ttrf, tc); rus = Rewriter::rewrite(rus); Assert(rus.isConst()); @@ -248,16 +172,18 @@ unsigned TaylorGenerator::getPolynomialApproximationBoundForArg( << "*** Increase Taylor bound to " << ds << " > " << d << " for (" << k << " " << c << ")" << std::endl; // must use sound upper bound - std::vector pboundss; + ApproximationBounds pboundss; getPolynomialApproximationBounds(k, ds, pboundss); - pbounds[2] = pboundss[2]; + pbounds.d_upperPos = pboundss.d_upperPos; } return ds; } return d; } -std::pair TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlModel& model) +std::pair TaylorGenerator::getTfModelBounds(Node tf, + std::uint64_t d, + NlModel& model) { // compute the model value of the argument Node c = model.computeAbstractModelValue(tf[0]); @@ -279,7 +205,7 @@ std::pair TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlM } bool isNeg = csign == -1; - std::vector pbounds; + ApproximationBounds pbounds; getPolynomialApproximationBoundForArg(k, c, d, pbounds); std::vector bounds; @@ -287,8 +213,8 @@ std::pair TaylorGenerator::getTfModelBounds(Node tf, unsigned d, NlM TNode tfs = tf[0]; for (unsigned d2 = 0; d2 < 2; d2++) { - int index = d2 == 0 ? (isNeg ? 1 : 0) : (isNeg ? 3 : 2); - Node pab = pbounds[index]; + Node pab = (d2 == 0 ? pbounds.d_lower + : (isNeg ? pbounds.d_upperNeg : pbounds.d_upperPos)); if (!pab.isNull()) { // { x -> M_A(tf[0]) } diff --git a/src/theory/arith/nl/transcendental/taylor_generator.h b/src/theory/arith/nl/transcendental/taylor_generator.h index 261583ca9..ccdc80e0f 100644 --- a/src/theory/arith/nl/transcendental/taylor_generator.h +++ b/src/theory/arith/nl/transcendental/taylor_generator.h @@ -27,6 +27,17 @@ namespace transcendental { class TaylorGenerator { public: + /** Stores the approximation bounds for transcendental functions */ + struct ApproximationBounds + { + /** Lower bound */ + Node d_lower; + /** Upper bound for negative values */ + Node d_upperNeg; + /** Upper bound for positive values */ + Node d_upperPos; + }; + TaylorGenerator(); /** @@ -35,23 +46,17 @@ class TaylorGenerator TNode getTaylorVariable(); /** - * Get Taylor series of degree n for function fa centered around point fa[0]. + * Get Taylor series of degree n for function fa centered around zero. * - * Return value is ( P_{n,f(a)}( x ), R_{n+1,f(a)}( x ) ) where + * Return value is ( P_{n,f(0)}( x ), R_{n+1,f(0)}( 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 + * P_{n,f(0)}( x ) = sum_{i=0}^n (f^i(0)/i!)*x^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} + * R_{n+1,f(0)}( x ) = x^{n+1}/(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. + * The above values are cached for each (f,n). */ - std::pair getTaylor(TNode fa, std::uint64_t n); + std::pair getTaylor(Kind k, std::uint64_t n); /** * polynomial approximation bounds @@ -68,8 +73,8 @@ class TaylorGenerator * given ( c ), use the function below. */ void getPolynomialApproximationBounds(Kind k, - unsigned d, - std::vector& pbounds); + std::uint64_t d, + ApproximationBounds& pbounds); /** * polynomial approximation bounds @@ -82,10 +87,8 @@ class TaylorGenerator * @return the actual degree of the polynomial approximations (which may be * larger than d). */ - unsigned getPolynomialApproximationBoundForArg(Kind k, - Node c, - unsigned d, - std::vector& pbounds); + std::uint64_t getPolynomialApproximationBoundForArg( + Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds); /** get transcendental function model bounds * @@ -93,22 +96,20 @@ class TaylorGenerator * function application tf based on Taylor of degree 2*d, which is dependent * on the model value of its argument. */ - std::pair getTfModelBounds(Node tf, unsigned d, NlModel& model); + std::pair getTfModelBounds(Node tf, + std::uint64_t d, + NlModel& model); private: NodeManager* d_nm; const Node d_taylor_real_fv; - const Node d_taylor_real_fv_base; - const Node d_taylor_real_fv_base_rem; - std::unordered_map, - NodeHashFunction> - s_taylor_sum; - std::unordered_map, - NodeHashFunction> - d_taylor_rem; - std::map>> d_poly_bounds; + + /** + * For every kind (EXP or SINE) and every degree we store the taylor series up + * to this degree and the next term in the series. + */ + std::map>> d_taylor_terms; + std::map> d_poly_bounds; }; } // namespace transcendental diff --git a/src/theory/arith/nl/transcendental/transcendental_solver.cpp b/src/theory/arith/nl/transcendental/transcendental_solver.cpp index 54f430372..ee422ebb6 100644 --- a/src/theory/arith/nl/transcendental/transcendental_solver.cpp +++ b/src/theory/arith/nl/transcendental/transcendental_solver.cpp @@ -267,13 +267,13 @@ bool TranscendentalSolver::checkTfTangentPlanesFun(Node tf, unsigned d) // Figure 3: P_l, P_u // mapped to for signs of c std::map poly_approx_bounds[2]; - std::vector pbounds; - unsigned actual_d = + TaylorGenerator::ApproximationBounds pbounds; + std::uint64_t 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]; - poly_approx_bounds[1][-1] = pbounds[3]; + poly_approx_bounds[0][1] = pbounds.d_lower; + poly_approx_bounds[0][-1] = pbounds.d_lower; + poly_approx_bounds[1][1] = pbounds.d_upperPos; + poly_approx_bounds[1][-1] = pbounds.d_upperNeg; // Figure 3 : v Node v = d_tstate.d_model.computeAbstractModelValue(tf); -- cgit v1.2.3