From 3aa568dbd217820a625a28f2e34b5547af3f0c4d Mon Sep 17 00:00:00 2001 From: Andrew Reynolds Date: Tue, 1 May 2018 10:47:40 -0500 Subject: [PATCH] Improve tangent planes for transcendental functions (#1832) --- src/options/arith_options.toml | 2 +- src/theory/arith/nonlinear_extension.cpp | 1056 +++++++++++++--------- src/theory/arith/nonlinear_extension.h | 66 +- test/regress/regress1/nl/cos1-tc.smt2 | 2 +- 4 files changed, 699 insertions(+), 427 deletions(-) diff --git a/src/options/arith_options.toml b/src/options/arith_options.toml index 1cccf0884..b2b5f0c31 100644 --- a/src/options/arith_options.toml +++ b/src/options/arith_options.toml @@ -455,7 +455,7 @@ header = "options/arith_options.h" category = "regular" long = "nl-ext-tf-tplanes" type = "bool" - default = "false" + default = "true" read_only = true help = "use non-terminating tangent plane strategy for transcendental functions for non-linear" diff --git a/src/theory/arith/nonlinear_extension.cpp b/src/theory/arith/nonlinear_extension.cpp index a04d39481..c2ebc55b8 100644 --- a/src/theory/arith/nonlinear_extension.cpp +++ b/src/theory/arith/nonlinear_extension.cpp @@ -826,13 +826,13 @@ std::vector NonlinearExtension::checkModel( Node atom = lit.getKind()==NOT ? lit[0] : lit; if( d_skolem_atoms.find( atom )==d_skolem_atoms.end() ){ Node litv = computeModelValue(lit); - Trace("nl-ext-mv") << "M[[ " << lit << " ]] -> " << litv; + Trace("nl-ext-mv-assert") << "M[[ " << lit << " ]] -> " << litv; if (litv != d_true) { - Trace("nl-ext-mv") << " [model-false]" << std::endl; + Trace("nl-ext-mv-assert") << " [model-false]" << std::endl; //Assert(litv == d_false); false_asserts.push_back(lit); } else { - Trace("nl-ext-mv") << std::endl; + Trace("nl-ext-mv-assert") << std::endl; } } } @@ -841,47 +841,176 @@ std::vector NonlinearExtension::checkModel( bool NonlinearExtension::checkModelTf(const std::vector& assertions) { - Trace("nl-ext-tf-check-model") << "check-model : Run" << std::endl; - if (!d_pi.isNull()) + Trace("nl-ext-cm") << "check-model : Run" << std::endl; + d_check_model_vars.clear(); + d_check_model_subs.clear(); + d_check_model_lit.clear(); + d_tf_check_model_bounds.clear(); + + // get the presubstitution + std::vector pvars; + std::vector psubs; + for (std::pair& tb : d_trig_base) { - // add bounds for PI - d_tf_check_model_bounds[d_pi] = - std::pair(d_pi_bound[0], d_pi_bound[1]); + pvars.push_back(tb.first); + psubs.push_back(tb.second); } - for (const std::pair >& tfb : - d_tf_check_model_bounds) + + // initialize representation of assertions + std::vector passertions; + for (const Node& a : assertions) { - Node tf = tfb.first; - Trace("nl-ext-tf-check-model") - << "check-model : satisfied approximate bound : "; - Trace("nl-ext-tf-check-model") << tfb.second.first << " <= " << tf - << " <= " << tfb.second.second << std::endl; + Node pa = a; + if (!pvars.empty()) + { + pa = + pa.substitute(pvars.begin(), pvars.end(), psubs.begin(), psubs.end()); + pa = Rewriter::rewrite(pa); + } + Trace("nl-ext-cm-assert") << "- assert : " << pa << std::endl; + d_check_model_lit[pa] = pa; + passertions.push_back(pa); } + // heuristically, solve for equalities + Trace("nl-ext-cm") << "solving equalities..." << std::endl; + unsigned nassertions_new = passertions.size(); + unsigned curr_index = 0; + unsigned terminate_index = 0; + do + { + Trace("nl-ext-cm-debug") << " indices : " << curr_index << " " + << terminate_index << " " << nassertions_new + << std::endl; + Node lit = passertions[curr_index]; + Node slit = d_check_model_lit[lit]; + Trace("nl-ext-cm-debug") << " process " << lit << std::endl; + // update it based on the current substitution + if (!d_check_model_vars.empty() && !slit.isConst()) + { + // reapply the substitution + slit = slit.substitute(d_check_model_vars.begin(), + d_check_model_vars.end(), + d_check_model_subs.begin(), + d_check_model_subs.end()); + slit = Rewriter::rewrite(slit); + d_check_model_lit[lit] = slit; + Trace("nl-ext-cm-debug") << " ...substituted to " << slit << std::endl; + } + // is it a substitution? + if (slit.getKind() == EQUAL) + { + std::map msum; + if (ArithMSum::getMonomialSumLit(slit, msum)) + { + // find a legal variable to solve for + Node v; + Node slv; + for (std::pair& m : msum) + { + Node mv = m.first; + if (mv.isVar() && ((v.getKind() != SKOLEM && mv.getKind() == SKOLEM) + || slv.isNull())) + { + Node veqc; + if (ArithMSum::isolate(mv, msum, veqc, slv, EQUAL) != 0) + { + Assert(veqc.isNull() && !slv.isNull()); + if (!mv.hasSubterm(v)) + { + v = mv; + } + } + } + } + if (!v.isNull()) + { + Trace("nl-ext-cm") + << " assertion : " << slit + << " can be turned into substitution:" << std::endl; + Trace("nl-ext-cm") << " " << v << " -> " << slv << std::endl; + d_check_model_vars.push_back(v); + d_check_model_subs.push_back(slv); + d_check_model_lit[lit] = d_true; + terminate_index = curr_index; + } + } + } + curr_index++; + if (curr_index == nassertions_new) + { + curr_index = 0; + } + } while (curr_index != terminate_index); + Trace("nl-ext-cm") << "...finished." << std::endl; + + // initialize the check model bounds + for (std::pair >& tfs : d_tf_rep_map) + { + Kind k = tfs.first; + for (std::pair& tfr : tfs.second) + { + // Figure 3 : tf( x ) + Node tf = tfr.second; + Node atf = computeModelValue(tf); + if (k == PI) + { + d_tf_check_model_bounds[atf] = + std::pair(d_pi_bound[0], d_pi_bound[1]); + } + else if (isRefineableTfFun(tf)) + { + d_tf_check_model_bounds[atf] = getTfModelBounds(tf, d_taylor_degree); + } + if (Trace.isOn("nl-ext-cm")) + { + std::map >::iterator it = + d_tf_check_model_bounds.find(atf); + if (it != d_tf_check_model_bounds.end()) + { + Trace("nl-ext-cm") << "check-model : satisfied approximate bound : "; + 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; + } + } + } + } + + std::unordered_set all_assertions; std::vector check_assertions; - for (const Node& a : assertions) + for (const Node& a : passertions) { - Node av = computeModelValue(a); - // simple check - if (!simpleCheckModelTfLit(av)) + Node as = d_check_model_lit[a]; + if (!as.isConst() || !as.getConst()) { - check_assertions.push_back(av); - Trace("nl-ext-tf-check-model") << "check-model : assertion : " << av - << " (from " << a << ")" << std::endl; + Node av = computeModelValue(a); + if (all_assertions.find(av) == all_assertions.end()) + { + all_assertions.insert(av); + // simple check + if (!simpleCheckModelTfLit(av)) + { + check_assertions.push_back(av); + Trace("nl-ext-cm") << "check-model : failed assertion : " << a + << std::endl; + Trace("nl-ext-cm-debug") + << "check-model : failed assertion, value : " << av << std::endl; + } + } } } if (check_assertions.empty()) { - Trace("nl-ext-tf-check-model") << "...simple check succeeded." << std::endl; + Trace("nl-ext-cm") << "...simple check succeeded." << std::endl; return true; } - else - { - Trace("nl-ext-tf-check-model") << "...simple check failed." << std::endl; - // TODO (#1450) check model for general case - return false; - } + + Trace("nl-ext-cm") << "...simple check failed." << std::endl; + // TODO (#1450) check model for general case + return false; } bool NonlinearExtension::simpleCheckModelTfLit(Node lit) @@ -1172,8 +1301,7 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, Node a = xts[i]; computeModelValue(a, 0); computeModelValue(a, 1); - Trace("nl-ext-mv") << " " << a << " -> " << d_mv[1][a] << " [actual: " - << d_mv[0][a] << " ]" << std::endl; + printModelValue("nl-ext-mv", a); //Assert(d_mv[1][a].isConst()); //Assert(d_mv[0][a].isConst()); if (a.getKind() == NONLINEAR_MULT) @@ -1317,7 +1445,31 @@ int NonlinearExtension::checkLastCall(const std::vector& assertions, registerMonomial(v); computeModelValue(v, 0); computeModelValue(v, 1); - Trace("nl-ext-mv") << " " << v << " -> " << d_mv[1][v] << " [actual: " << d_mv[0][v] << " ]" << std::endl; + printModelValue("nl-ext-mv", v); + } + if (Trace.isOn("nl-ext-mv")) + { + Trace("nl-ext-mv") << "Arguments of trancendental functions : " + << std::endl; + for (std::map >::iterator it = + d_tf_rep_map.begin(); + it != d_tf_rep_map.end(); + ++it) + { + Kind k = it->first; + if (k == SINE || k == EXPONENTIAL) + { + for (std::map::iterator itt = it->second.begin(); + itt != it->second.end(); + ++itt) + { + Node v = itt->second[0]; + computeModelValue(v, 0); + computeModelValue(v, 1); + printModelValue("nl-ext-mv", v); + } + } + } } //----------------------------------- possibly split on zero @@ -1638,8 +1790,6 @@ void NonlinearExtension::check(Theory::Effort e) { { d_taylor_degree++; needsRecheck = true; - // clear secant points - d_secant_points.clear(); // increase precision for PI? // Difficult since Taylor series is very slow to converge Trace("nl-ext") << "...increment Taylor degree to " << d_taylor_degree @@ -3234,8 +3384,6 @@ std::vector NonlinearExtension::checkTranscendentalTangentPlanes() std::vector lemmas; Trace("nl-ext") << "Get tangent plane lemmas for transcendental functions..." << std::endl; - - NodeManager* nm = NodeManager::currentNM(); // this implements Figure 3 of "Satisfiaility Modulo Transcendental Functions // via Incremental Linearization" by Cimatti et al for (std::pair >& tfs : d_tf_rep_map) @@ -3248,413 +3396,385 @@ std::vector NonlinearExtension::checkTranscendentalTangentPlanes() // initial approximation is superior. continue; } - Node tft = nm->mkNode(k, d_zero); - Trace("nl-ext-tf-tplanes-debug") << "Taylor variables: " << std::endl; - Trace("nl-ext-tf-tplanes-debug") + Trace("nl-ext-tftp-debug2") << "Taylor variables: " << std::endl; + Trace("nl-ext-tftp-debug2") << " taylor_real_fv : " << d_taylor_real_fv << std::endl; - Trace("nl-ext-tf-tplanes-debug") + Trace("nl-ext-tftp-debug2") << " taylor_real_fv_base : " << d_taylor_real_fv_base << std::endl; - Trace("nl-ext-tf-tplanes-debug") + Trace("nl-ext-tftp-debug2") << " taylor_real_fv_base_rem : " << d_taylor_real_fv_base_rem << std::endl; - Trace("nl-ext-tf-tplanes-debug") << std::endl; + Trace("nl-ext-tftp-debug2") << std::endl; // we substitute into the Taylor sum P_{n,f(0)}( x ) - std::vector taylor_vars; - taylor_vars.push_back(d_taylor_real_fv); - - // Figure 3: P_l, P_u - // mapped to for signs of c - std::map poly_approx_bounds[2]; - // n is the Taylor degree we are currently considering - unsigned n = 2 * d_taylor_degree; - // n must be even - std::pair taylor = getTaylor(tft, n); - Trace("nl-ext-tf-tplanes-debug") << "Taylor for " << k - << " is : " << taylor.first << std::endl; - Node taylor_sum = Rewriter::rewrite(taylor.first); - Trace("nl-ext-tf-tplanes-debug") << "Taylor for " << k - << " is (post-rewrite) : " << taylor_sum - << std::endl; - Assert(taylor.second.getKind() == MULT); - Assert(taylor.second.getNumChildren() == 2); - Assert(taylor.second[0].getKind() == DIVISION); - Trace("nl-ext-tf-tplanes-debug") << "Taylor remainder for " << k << " is " - << taylor.second << std::endl; - // ru is x^{n+1}/(n+1)! - Node ru = nm->mkNode(DIVISION, taylor.second[1], taylor.second[0][1]); - ru = Rewriter::rewrite(ru); - Trace("nl-ext-tf-tplanes-debug") - << "Taylor remainder factor is (post-rewrite) : " << ru << std::endl; - if (k == EXPONENTIAL) - { - poly_approx_bounds[0][1] = taylor_sum; - poly_approx_bounds[0][-1] = taylor_sum; - poly_approx_bounds[1][1] = Rewriter::rewrite( - nm->mkNode(MULT, taylor_sum, nm->mkNode(PLUS, d_one, ru))); - poly_approx_bounds[1][-1] = - Rewriter::rewrite(nm->mkNode(PLUS, taylor_sum, ru)); - } - else - { - Assert(k == SINE); - poly_approx_bounds[0][1] = - Rewriter::rewrite(nm->mkNode(MINUS, taylor_sum, ru)); - poly_approx_bounds[0][-1] = poly_approx_bounds[0][1]; - poly_approx_bounds[1][1] = - Rewriter::rewrite(nm->mkNode(PLUS, taylor_sum, ru)); - poly_approx_bounds[1][-1] = poly_approx_bounds[1][1]; - } - Trace("nl-ext-tf-tplanes") << "Polynomial approximation for " << k - << " is: " << std::endl; - Trace("nl-ext-tf-tplanes") << " Lower (pos): " << poly_approx_bounds[0][1] - << std::endl; - Trace("nl-ext-tf-tplanes") << " Upper (pos): " << poly_approx_bounds[1][1] - << std::endl; - Trace("nl-ext-tf-tplanes") << " Lower (neg): " << poly_approx_bounds[0][-1] - << std::endl; - Trace("nl-ext-tf-tplanes") << " Upper (neg): " << poly_approx_bounds[1][-1] - << std::endl; for (std::pair& tfr : tfs.second) { // Figure 3 : tf( x ) Node tf = tfr.second; - - bool consider = true; - if (k == SINE) + if (isRefineableTfFun(tf)) { - // we do not consider e.g. sin( -1*x ), since considering sin( x ) will - // have the same effect - consider = tf[0].isVar(); - } - int csign; - Node c; - if (consider) - { - // Figure 3 : c - c = computeModelValue(tf[0], 1); - Assert(c.isConst()); - csign = c.getConst().sgn(); - consider = csign != 0; + Trace("nl-ext-tftp") << "Compute tangent planes " << tf << std::endl; + // go until max degree is reached, or we don't meet bound criteria + for (unsigned d = 1; d <= d_taylor_degree; d++) + { + Trace("nl-ext-tftp") << "- run at degree " << d << "..." << std::endl; + unsigned prev = lemmas.size(); + if (!checkTfTangentPlanesFun(tf, d, lemmas)) + { + Trace("nl-ext-tftp") + << "...fail, #lemmas = " << (lemmas.size() - prev) << std::endl; + break; + } + else + { + Trace("nl-ext-tftp") << "...success" << std::endl; + } + } } + } + } - if (consider) - { - Assert(csign == 1 || csign == -1); + return lemmas; +} - // Figure 3 : v - Node v = computeModelValue(tf, 1); +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().sgn(); + if (csign == 0) + { + return false; + } + return true; +} - // check value of tf - Trace("nl-ext-tf-tplanes") << "Process tangent plane refinement for " - << tf << "..." << std::endl; - Trace("nl-ext-tf-tplanes") << " value in model : v = " << v - << std::endl; - Trace("nl-ext-tf-tplanes") << " arg value in model : c = " << c - << std::endl; +bool NonlinearExtension::checkTfTangentPlanesFun(Node tf, + unsigned d, + std::vector& lemmas) +{ + Assert(isRefineableTfFun(tf)); - // compute the concavity - int region = -1; - std::unordered_map::iterator itr = - d_tf_region.find(tf); - if (itr != d_tf_region.end()) + NodeManager* nm = NodeManager::currentNM(); + Kind k = tf.getKind(); + // Figure 3: P_l, P_u + // mapped to for signs of c + std::map poly_approx_bounds[2]; + std::vector pbounds; + getPolynomialApproximationBounds(k, 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]; + + // Figure 3 : c + Node c = computeModelValue(tf[0], 1); + int csign = c.getConst().sgn(); + Assert(csign == 1 || csign == -1); + + // Figure 3 : v + Node v = computeModelValue(tf, 1); + + // check value of tf + Trace("nl-ext-tftp-debug") << "Process tangent plane refinement for " << tf + << ", degree " << d << "..." << std::endl; + Trace("nl-ext-tftp-debug") << " value in model : " << v << std::endl; + Trace("nl-ext-tftp-debug") << " arg value in model : " << c << std::endl; + + std::vector taylor_vars; + taylor_vars.push_back(d_taylor_real_fv); + + // compute the concavity + int region = -1; + std::unordered_map::iterator itr = + d_tf_region.find(tf); + if (itr != d_tf_region.end()) + { + region = itr->second; + Trace("nl-ext-tftp-debug") << " region is : " << region << std::endl; + } + // Figure 3 : conc + int concavity = regionToConcavity(k, itr->second); + Trace("nl-ext-tftp-debug") << " concavity is : " << concavity << std::endl; + if (concavity == 0) + { + return false; + } + // bounds for which we are this concavity + // Figure 3: < l, u > + Node bounds[2]; + if (k == SINE) + { + bounds[0] = regionToLowerBound(k, region); + Assert(!bounds[0].isNull()); + bounds[1] = regionToUpperBound(k, region); + Assert(!bounds[1].isNull()); + } + + // Figure 3: P + Node poly_approx; + + // compute whether this is a tangent refinement or a secant refinement + bool is_tangent = false; + bool is_secant = false; + std::pair mvb = getTfModelBounds(tf, d); + for (unsigned r = 0; r < 2; r++) + { + Node pab = poly_approx_bounds[r][csign]; + Node v_pab = r == 0 ? mvb.first : mvb.second; + if (!v_pab.isNull()) + { + Assert(v_pab.isConst()); + Trace("nl-ext-tftp-debug2") << "...model value of " << pab << " is " + << v_pab << std::endl; + Node comp = nm->mkNode(r == 0 ? LT : GT, v, v_pab); + Trace("nl-ext-tftp-debug2") << "...compare : " << comp << std::endl; + Node compr = Rewriter::rewrite(comp); + Trace("nl-ext-tftp-debug2") << "...got : " << compr << std::endl; + if (compr == d_true) + { + // beyond the bounds + if (r == 0) { - region = itr->second; - Trace("nl-ext-tf-tplanes") << " region is : " << region << std::endl; + poly_approx = poly_approx_bounds[r][csign]; + is_tangent = concavity == 1; + is_secant = concavity == -1; } - // Figure 3 : conc - int concavity = regionToConcavity(k, itr->second); - Trace("nl-ext-tf-tplanes") << " concavity is : " << concavity - << std::endl; - if (concavity != 0) + else { - // bounds for which we are this concavity - // Figure 3: < l, u > - Node bounds[2]; - if (k == SINE) - { - bounds[0] = regionToLowerBound(k, region); - Assert(!bounds[0].isNull()); - bounds[1] = regionToUpperBound(k, region); - Assert(!bounds[1].isNull()); - } - - // Figure 3: P - Node poly_approx; - - // compute whether this is a tangent refinement or a secant refinement - bool is_tangent = false; - bool is_secant = false; - std::map model_values; - for (unsigned d = 0; d < 2; d++) - { - Node pab = poly_approx_bounds[d][csign]; - if (!pab.isNull()) - { - // { x -> tf[0] } - std::vector taylor_subs; - taylor_subs.push_back(tf[0]); - Assert(taylor_vars.size() == taylor_subs.size()); - pab = pab.substitute(taylor_vars.begin(), - taylor_vars.end(), - taylor_subs.begin(), - taylor_subs.end()); - pab = Rewriter::rewrite(pab); - Node v_pab = computeModelValue(pab, 1); - model_values[d] = v_pab; - Assert(v_pab.isConst()); - Trace("nl-ext-tf-tplanes-debug") << "...model value of " << pab - << " is " << v_pab << std::endl; - Node comp = nm->mkNode(d == 0 ? LT : GT, v, v_pab); - Trace("nl-ext-tf-tplanes-debug") << "...compare : " << comp - << std::endl; - Node compr = Rewriter::rewrite(comp); - Trace("nl-ext-tf-tplanes-debug") << "...got : " << compr - << std::endl; - if (compr == d_true) - { - // beyond the bounds - if (d == 0) - { - poly_approx = poly_approx_bounds[d][csign]; - is_tangent = concavity == 1; - is_secant = concavity == -1; - } - else - { - poly_approx = poly_approx_bounds[d][csign]; - is_tangent = concavity == -1; - is_secant = concavity == 1; - } - Trace("nl-ext-tf-tplanes") << "*** Outside boundary point ("; - Trace("nl-ext-tf-tplanes") << (d == 0 ? "low" : "high") << ") "; - Trace("nl-ext-tf-tplanes") << comp << ", will refine..." - << std::endl; - Trace("nl-ext-tf-tplanes") - << " poly_approx = " << poly_approx << std::endl; - Trace("nl-ext-tf-tplanes") << " is_tangent = " << is_tangent - << std::endl; - Trace("nl-ext-tf-tplanes") << " is_secant = " << is_secant - << std::endl; - break; - } - else - { - Trace("nl-ext-tf-tplanes") << " ...within " - << (d == 0 ? "low" : "high") - << " bound : "; - Trace("nl-ext-tf-tplanes") << comp << std::endl; - } - } - } - - // Figure 3: P( c ) - Node poly_approx_c; - if (is_tangent || is_secant) - { - Assert(!poly_approx.isNull()); - std::vector taylor_subs; - taylor_subs.push_back(c); - Assert(taylor_vars.size() == taylor_subs.size()); - poly_approx_c = poly_approx.substitute(taylor_vars.begin(), - taylor_vars.end(), - taylor_subs.begin(), - taylor_subs.end()); - Trace("nl-ext-tf-tplanes-debug") << "...poly appoximation at c is " - << poly_approx_c << std::endl; - } - else - { - // store for check model bounds - Node atf = computeModelValue(tf); - d_tf_check_model_bounds[atf] = - std::pair(model_values[0], model_values[1]); - } + poly_approx = poly_approx_bounds[r][csign]; + is_tangent = concavity == -1; + is_secant = concavity == 1; + } + if (Trace.isOn("nl-ext-tftp")) + { + Trace("nl-ext-tftp") << "*** Outside boundary point ("; + Trace("nl-ext-tftp") << (r == 0 ? "low" : "high") << ") "; + printRationalApprox("nl-ext-tftp", v_pab); + Trace("nl-ext-tftp") << ", will refine..." << std::endl; + Trace("nl-ext-tftp") << " poly_approx = " << poly_approx + << std::endl; + Trace("nl-ext-tftp") << " is_tangent = " << is_tangent + << std::endl; + Trace("nl-ext-tftp") << " is_secant = " << is_secant << std::endl; + } + break; + } + else + { + Trace("nl-ext-tftp") << " ...within " << (r == 0 ? "low" : "high") + << " bound : "; + printRationalApprox("nl-ext-tftp", v_pab); + Trace("nl-ext-tftp") << std::endl; + } + } + } - if (is_tangent) - { - // compute tangent plane - // Figure 3: T( x ) - Node tplane; - Node poly_approx_deriv = - getDerivative(poly_approx, d_taylor_real_fv); - Assert(!poly_approx_deriv.isNull()); - poly_approx_deriv = Rewriter::rewrite(poly_approx_deriv); - Trace("nl-ext-tf-tplanes-debug") << "...derivative of " - << poly_approx << " is " - << poly_approx_deriv << std::endl; - std::vector taylor_subs; - taylor_subs.push_back(c); - Assert(taylor_vars.size() == taylor_subs.size()); - Node poly_approx_c_deriv = - poly_approx_deriv.substitute(taylor_vars.begin(), - taylor_vars.end(), - taylor_subs.begin(), - taylor_subs.end()); - tplane = nm->mkNode( - PLUS, - poly_approx_c, - nm->mkNode( - MULT, poly_approx_c_deriv, nm->mkNode(MINUS, tf[0], c))); - - Node lem = nm->mkNode(concavity == 1 ? GEQ : LEQ, tf, tplane); - std::vector antec; - for (unsigned i = 0; i < 2; i++) - { - if (!bounds[i].isNull()) - { - antec.push_back( - nm->mkNode(i == 0 ? GEQ : LEQ, tf[0], bounds[i])); - } - } - if (!antec.empty()) - { - Node antec_n = - antec.size() == 1 ? antec[0] : nm->mkNode(AND, antec); - lem = nm->mkNode(IMPLIES, antec_n, lem); - } - Trace("nl-ext-tf-tplanes-debug") - << "*** Tangent plane lemma (pre-rewrite): " << lem - << std::endl; - lem = Rewriter::rewrite(lem); - Trace("nl-ext-tf-tplanes") << "*** Tangent plane lemma : " << lem - << std::endl; - // Figure 3 : line 9 - lemmas.push_back(lem); - } - else if (is_secant) - { - // bounds are the minimum and maximum previous secant points - Assert(std::find(d_secant_points[tf].begin(), - d_secant_points[tf].end(), - c) - == d_secant_points[tf].end()); - // insert into the vector - d_secant_points[tf].push_back(c); - // sort - SortNonlinearExtension smv; - smv.d_nla = this; - smv.d_order_type = 0; - std::sort( - d_secant_points[tf].begin(), d_secant_points[tf].end(), smv); - // get the resulting index of c - unsigned index = - std::find( - d_secant_points[tf].begin(), d_secant_points[tf].end(), c) - - d_secant_points[tf].begin(); - // bounds are the next closest upper/lower bound values - if (index > 0) - { - bounds[0] = d_secant_points[tf][index - 1]; - } - else - { - // otherwise, we use the lower boundary point for this concavity - // region - if (k == SINE) - { - Assert(!bounds[0].isNull()); - } - else if (k == EXPONENTIAL) - { - // pick c-1 - bounds[0] = Rewriter::rewrite(nm->mkNode(MINUS, c, d_one)); - } - } - if (index < d_secant_points[tf].size() - 1) - { - bounds[1] = d_secant_points[tf][index + 1]; - } - else - { - // otherwise, we use the upper boundary point for this concavity - // region - if (k == SINE) - { - Assert(!bounds[1].isNull()); - } - else if (k == EXPONENTIAL) - { - // pick c+1 - bounds[1] = Rewriter::rewrite(nm->mkNode(PLUS, c, d_one)); - } - } - Trace("nl-ext-tf-tplanes-debug") - << "...secant bounds are : " << bounds[0] << " ... " - << bounds[1] << std::endl; + // Figure 3: P( c ) + Node poly_approx_c; + if (is_tangent || is_secant) + { + Assert(!poly_approx.isNull()); + std::vector taylor_subs; + taylor_subs.push_back(c); + Assert(taylor_vars.size() == taylor_subs.size()); + poly_approx_c = poly_approx.substitute(taylor_vars.begin(), + taylor_vars.end(), + taylor_subs.begin(), + taylor_subs.end()); + Trace("nl-ext-tftp-debug2") << "...poly approximation at c is " + << poly_approx_c << std::endl; + } + else + { + // we may want to continue getting better bounds + return true; + } - for (unsigned s = 0; s < 2; s++) - { - // compute secant plane - 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); - Trace("nl-ext-tf-tplanes-debug") << "...model value of bound " - << bounds[s] << " is " << b - << std::endl; - Assert(b.isConst()); - if (c != b) - { - // Figure 3 : P(l), P(u), for s = 0,1 - Node poly_approx_b; - std::vector taylor_subs; - taylor_subs.push_back(b); - Assert(taylor_vars.size() == taylor_subs.size()); - poly_approx_b = poly_approx.substitute(taylor_vars.begin(), - taylor_vars.end(), - taylor_subs.begin(), - taylor_subs.end()); - // Figure 3: S_l( x ), S_u( x ) for s = 0,1 - Node splane; - Node rcoeff_n = Rewriter::rewrite(nm->mkNode(MINUS, b, c)); - Assert(rcoeff_n.isConst()); - Rational rcoeff = rcoeff_n.getConst(); - Assert(rcoeff.sgn() != 0); - splane = nm->mkNode( - PLUS, - poly_approx_b, - nm->mkNode(MULT, - nm->mkNode(MINUS, poly_approx_b, poly_approx_c), - nm->mkConst(Rational(1) / rcoeff), - nm->mkNode(MINUS, tf[0], b))); - - Node lem = nm->mkNode(concavity == 1 ? LEQ : GEQ, tf, splane); - // With respect to Figure 3, this is slightly different. - // In particular, we chose b to be the model value of bounds[s], - // which is a constant although bounds[s] may not be (e.g. if it - // contains PI). - // To ensure that c...b does not cross an inflection point, - // we guard with the symbolic version of bounds[s]. - // This leads to lemmas e.g. of this form: - // ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x - - // b ) + P( b ) ) - // where b = (PI/2)^M, the current value of PI/2 in the model. - // This is sound since we are guarded by the symbolic - // representation of PI/2. - Node antec_n = - nm->mkNode(AND, - nm->mkNode(GEQ, tf[0], s == 0 ? bounds[s] : c), - nm->mkNode(LEQ, tf[0], s == 0 ? c : bounds[s])); - lem = nm->mkNode(IMPLIES, antec_n, lem); - Trace("nl-ext-tf-tplanes-debug") - << "*** Secant plane lemma (pre-rewrite) : " << lem - << std::endl; - lem = Rewriter::rewrite(lem); - Trace("nl-ext-tf-tplanes") << "*** Secant plane lemma : " << lem - << std::endl; - // Figure 3 : line 22 - lemmas.push_back(lem); - } - } - } - } + if (is_tangent) + { + // compute tangent plane + // Figure 3: T( x ) + Node tplane; + Node poly_approx_deriv = getDerivative(poly_approx, d_taylor_real_fv); + Assert(!poly_approx_deriv.isNull()); + poly_approx_deriv = Rewriter::rewrite(poly_approx_deriv); + Trace("nl-ext-tftp-debug2") << "...derivative of " << poly_approx << " is " + << poly_approx_deriv << std::endl; + std::vector taylor_subs; + taylor_subs.push_back(c); + Assert(taylor_vars.size() == taylor_subs.size()); + Node poly_approx_c_deriv = poly_approx_deriv.substitute(taylor_vars.begin(), + taylor_vars.end(), + taylor_subs.begin(), + taylor_subs.end()); + tplane = nm->mkNode( + PLUS, + poly_approx_c, + nm->mkNode(MULT, poly_approx_c_deriv, nm->mkNode(MINUS, tf[0], c))); + + Node lem = nm->mkNode(concavity == 1 ? GEQ : LEQ, tf, tplane); + std::vector antec; + for (unsigned i = 0; i < 2; i++) + { + if (!bounds[i].isNull()) + { + Node ant = nm->mkNode(i == 0 ? GEQ : LEQ, tf[0], bounds[i]); + antec.push_back(ant); } } + if (!antec.empty()) + { + Node antec_n = antec.size() == 1 ? antec[0] : nm->mkNode(AND, antec); + lem = nm->mkNode(IMPLIES, antec_n, lem); + } + Trace("nl-ext-tftp-debug2") + << "*** Tangent plane lemma (pre-rewrite): " << lem << std::endl; + lem = Rewriter::rewrite(lem); + Trace("nl-ext-tftp-lemma") << "*** Tangent plane lemma : " << lem + << std::endl; + Assert(computeModelValue(lem, 1) == d_false); + // Figure 3 : line 9 + lemmas.push_back(lem); } + else if (is_secant) + { + // 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[tf][d].begin(), d_secant_points[tf][d].end(), c) + == d_secant_points[tf][d].end()); + // insert into the vector + d_secant_points[tf][d].push_back(c); + // sort + SortNonlinearExtension smv; + smv.d_nla = this; + smv.d_order_type = 0; + std::sort( + d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), smv); + // get the resulting index of c + unsigned index = + std::find( + d_secant_points[tf][d].begin(), d_secant_points[tf][d].end(), c) + - d_secant_points[tf][d].begin(); + // bounds are the next closest upper/lower bound values + if (index > 0) + { + bounds[0] = d_secant_points[tf][d][index - 1]; + } + else + { + // otherwise, we use the lower boundary point for this concavity + // region + if (k == SINE) + { + Assert(!bounds[0].isNull()); + } + else if (k == EXPONENTIAL) + { + // pick c-1 + bounds[0] = Rewriter::rewrite(nm->mkNode(MINUS, c, d_one)); + } + } + if (index < d_secant_points[tf][d].size() - 1) + { + bounds[1] = d_secant_points[tf][d][index + 1]; + } + else + { + // otherwise, we use the upper boundary point for this concavity + // region + if (k == SINE) + { + Assert(!bounds[1].isNull()); + } + else if (k == EXPONENTIAL) + { + // pick c+1 + bounds[1] = Rewriter::rewrite(nm->mkNode(PLUS, c, d_one)); + } + } + Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds[0] + << " ... " << bounds[1] << std::endl; - return lemmas; + for (unsigned s = 0; s < 2; s++) + { + // compute secant plane + 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); + Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds[s] + << " is " << b << std::endl; + Assert(b.isConst()); + if (c != b) + { + // Figure 3 : P(l), P(u), for s = 0,1 + Node poly_approx_b; + std::vector taylor_subs; + taylor_subs.push_back(b); + Assert(taylor_vars.size() == taylor_subs.size()); + poly_approx_b = poly_approx.substitute(taylor_vars.begin(), + taylor_vars.end(), + taylor_subs.begin(), + taylor_subs.end()); + // Figure 3: S_l( x ), S_u( x ) for s = 0,1 + Node splane; + Node rcoeff_n = Rewriter::rewrite(nm->mkNode(MINUS, b, c)); + Assert(rcoeff_n.isConst()); + Rational rcoeff = rcoeff_n.getConst(); + Assert(rcoeff.sgn() != 0); + splane = nm->mkNode( + PLUS, + poly_approx_b, + nm->mkNode(MULT, + nm->mkNode(MINUS, poly_approx_b, poly_approx_c), + nm->mkConst(Rational(1) / rcoeff), + nm->mkNode(MINUS, tf[0], b))); + + Node lem = nm->mkNode(concavity == 1 ? LEQ : GEQ, tf, splane); + // With respect to Figure 3, this is slightly different. + // In particular, we chose b to be the model value of bounds[s], + // which is a constant although bounds[s] may not be (e.g. if it + // contains PI). + // To ensure that c...b does not cross an inflection point, + // we guard with the symbolic version of bounds[s]. + // This leads to lemmas e.g. of this form: + // ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x - + // b ) + P( b ) ) + // where b = (PI/2)^M, the current value of PI/2 in the model. + // This is sound since we are guarded by the symbolic + // representation of PI/2. + Node antec_n = + nm->mkNode(AND, + nm->mkNode(GEQ, tf[0], s == 0 ? bounds[s] : c), + nm->mkNode(LEQ, tf[0], s == 0 ? c : bounds[s])); + lem = nm->mkNode(IMPLIES, antec_n, lem); + Trace("nl-ext-tftp-debug2") + << "*** Secant plane lemma (pre-rewrite) : " << lem << std::endl; + lem = Rewriter::rewrite(lem); + Trace("nl-ext-tftp-lemma") << "*** Secant plane lemma : " << lem + << std::endl; + // Figure 3 : line 22 + lemmas.push_back(lem); + Assert(computeModelValue(lem, 1) == d_false); + } + } + } + return false; } int NonlinearExtension::regionToMonotonicityDir(Kind k, int region) @@ -3834,13 +3954,13 @@ Node NonlinearExtension::getDerivative(Node n, Node x) return Node::null(); } -std::pair NonlinearExtension::getTaylor(TNode fa, unsigned n) +std::pair NonlinearExtension::getTaylor(Node fa, unsigned n) { + 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. + // optimization : simpler to compute (x-fa[0])^n if we are centered around 0 fac = fa; } else @@ -3949,6 +4069,102 @@ std::pair NonlinearExtension::getTaylor(TNode fa, unsigned n) return std::pair(taylor_sum, taylor_rem); } +void NonlinearExtension::getPolynomialApproximationBounds( + Kind k, unsigned d, std::vector& pbounds) +{ + if (d_poly_bounds[k][d].empty()) + { + NodeManager* nm = NodeManager::currentNM(); + Node tft = nm->mkNode(k, d_zero); + // n is the Taylor degree we are currently considering + unsigned 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() == MULT); + Assert(taylor.second.getNumChildren() == 2); + Assert(taylor.second[0].getKind() == DIVISION); + Trace("nl-ext-tftp-debug2") << "Taylor remainder for " << k << " is " + << taylor.second << std::endl; + // ru is x^{n+1}/(n+1)! + Node ru = nm->mkNode(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; + if (k == EXPONENTIAL) + { + pbounds.push_back(taylor_sum); + pbounds.push_back(taylor_sum); + pbounds.push_back(Rewriter::rewrite( + nm->mkNode(MULT, taylor_sum, nm->mkNode(PLUS, d_one, ru)))); + pbounds.push_back(Rewriter::rewrite(nm->mkNode(PLUS, taylor_sum, ru))); + } + else + { + Assert(k == SINE); + Node l = Rewriter::rewrite(nm->mkNode(MINUS, taylor_sum, ru)); + Node u = Rewriter::rewrite(nm->mkNode(PLUS, taylor_sum, ru)); + pbounds.push_back(l); + pbounds.push_back(l); + pbounds.push_back(u); + pbounds.push_back(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()); + } + else + { + pbounds.insert( + pbounds.end(), d_poly_bounds[k][d].begin(), d_poly_bounds[k][d].end()); + } +} + +std::pair NonlinearExtension::getTfModelBounds(Node tf, unsigned d) +{ + // compute the model value of the argument + Node c = computeModelValue(tf[0], 1); + Assert(c.isConst()); + int csign = c.getConst().sgn(); + Assert(csign != 0); + bool isNeg = csign == -1; + + std::vector pbounds; + getPolynomialApproximationBounds(tf.getKind(), d, pbounds); + + std::vector bounds; + TNode tfv = d_taylor_real_fv; + TNode tfs = tf[0]; + for (unsigned d = 0; d < 2; d++) + { + int index = d == 0 ? (isNeg ? 1 : 0) : (isNeg ? 3 : 2); + Node pab = pbounds[index]; + if (!pab.isNull()) + { + // { x -> tf[0] } + pab = pab.substitute(tfv, tfs); + pab = Rewriter::rewrite(pab); + Node v_pab = computeModelValue(pab, 1); + bounds.push_back(v_pab); + } + else + { + bounds.push_back(Node::null()); + } + } + return std::pair(bounds[0], bounds[1]); +} + } // namespace arith } // namespace theory } // namespace CVC4 diff --git a/src/theory/arith/nonlinear_extension.h b/src/theory/arith/nonlinear_extension.h index 20c239ea7..00792a047 100644 --- a/src/theory/arith/nonlinear_extension.h +++ b/src/theory/arith/nonlinear_extension.h @@ -457,6 +457,22 @@ class NonlinearExtension { std::map > d_c_info_maxm; std::vector d_constraints; + // per last-call effort + + /** + * 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 d_check_model_vars; + std::vector d_check_model_subs; + /** + * Map from all literals appearing in the current set of assertions to their + * rewritten form under the substitution given by d_check_model_solve_form. + */ + std::map d_check_model_lit; + // model values/orderings /** cache of model values * @@ -534,9 +550,13 @@ class NonlinearExtension { * "get-previous-secant-points" in "Satisfiability * Modulo Transcendental Functions via Incremental * Linearization" by Cimatti et al., CADE 2017, for - * each transcendental function application. + * each transcendental function application. We store this set for each + * Taylor degree. */ - std::unordered_map, NodeHashFunction> d_secant_points; + std::unordered_map >, + NodeHashFunction> + d_secant_points; /** get Taylor series of degree n for function fa centered around point fa[0]. * @@ -554,7 +574,7 @@ class NonlinearExtension { * 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); + std::pair getTaylor(Node fa, unsigned n); /** internal variables used for constructing (cached) versions of the Taylor * series above. @@ -568,7 +588,6 @@ class NonlinearExtension { d_taylor_sum; std::unordered_map, NodeHashFunction> d_taylor_rem; - /** taylor degree * * Indicates that the degree of the polynomials in the Taylor approximation of @@ -577,7 +596,33 @@ class NonlinearExtension { * if the option options::nlExtTfIncPrecision() is enabled. */ unsigned d_taylor_degree; - + /** polynomial approximation bounds + * + * This adds P_l+, P_l-, P_u+, P_u- to pbounds, where these are polynomial + * approximations of the Taylor series of ( 0 ) for degree 2*d where + * k is SINE or EXPONENTIAL. + * These correspond to P_l and P_u from Figure 3 of Cimatti et al., CADE 2017, + * for positive/negative (+/-) values of the argument of ( 0 ). + */ + void getPolynomialApproximationBounds(Kind k, + unsigned d, + std::vector& pbounds); + /** cache of the above function */ + std::map > > d_poly_bounds; + /** get transcendental function model bounds + * + * This returns the current lower and upper bounds of transcendental + * 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); + /** 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 a lower/upper approximation of the constant r within the given * level of precision. In other words, this returns a constant c' such that @@ -871,6 +916,17 @@ class NonlinearExtension { * such that c1 ~= .277 and c2 ~= 2.032. */ std::vector checkTranscendentalTangentPlanes(); + /** check transcendental function refinement for tf + * + * This method is called by the above method for each refineable + * transcendental function (see isRefineableTfFun) that occurs in an + * assertion in the current context. + * + * This runs Figure 3 of Cimatti et al., CADE 2017 for transcendental + * function application tf for Taylor degree d. It may add a secant or + * tangent plane lemma to lems. + */ + bool checkTfTangentPlanesFun(Node tf, unsigned d, std::vector& lems); //-------------------------------------------- end lemma schemas }; /* class NonlinearExtension */ diff --git a/test/regress/regress1/nl/cos1-tc.smt2 b/test/regress/regress1/nl/cos1-tc.smt2 index 7ddae1453..f910f0c58 100644 --- a/test/regress/regress1/nl/cos1-tc.smt2 +++ b/test/regress/regress1/nl/cos1-tc.smt2 @@ -1,4 +1,4 @@ -; COMMAND-LINE: --nl-ext --no-nl-ext-tf-inc-prec +; COMMAND-LINE: --nl-ext --no-nl-ext-tf-tplanes --no-nl-ext-tf-inc-prec ; EXPECT: unknown (set-logic UFNRA) (declare-fun f (Real) Real) -- 2.30.2