diff --git a/src/math/lp/lp_bound_propagator.cpp b/src/math/lp/lp_bound_propagator.cpp index 10f155027..bf7171cd8 100644 --- a/src/math/lp/lp_bound_propagator.cpp +++ b/src/math/lp/lp_bound_propagator.cpp @@ -119,7 +119,7 @@ void lp_bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool co } bool lp_bound_propagator::nl_monomial_upper_bound_should_be_taken(unsigned j) const { return (!upper_bound_is_available_for_column(j)) || ( - nl_monomial_upper_bound_is_available(j) && m_nla_solver->get_lower_bound(j) < m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j]); + nl_monomial_upper_bound_is_available(j) && m_nla_solver->get_upper_bound(j) < m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j]); } bool lp_bound_propagator::nl_monomial_lower_bound_should_be_taken(unsigned j) const { @@ -142,7 +142,11 @@ void lp_bound_propagator::explain_implied_bound(implied_bound & ib) { int sign = j_sign * a_sign; if (sign > 0) { if (nl_monomial_upper_bound_should_be_taken(j)) { - SASSERT(false); + TRACE("nla_intervals", tout << "using the monomial upper bound\n";); + svector expl; + m_nla_solver->get_explanation_of_upper_bound_for_monomial(j, expl); + for (auto c : expl) + consume(a, c); } else { const ul_pair & ul = m_lar_solver.m_columns_to_ul_pairs[j]; auto witness = ul.upper_bound_witness(); @@ -151,7 +155,12 @@ void lp_bound_propagator::explain_implied_bound(implied_bound & ib) { } } else { if (nl_monomial_lower_bound_should_be_taken(j)) { - SASSERT(false); + TRACE("nla_intervals", tout << "using the monomial lower bound\n";); + svector expl; + m_nla_solver->get_explanation_of_lower_bound_for_monomial(j, expl); + for (auto c : expl) + consume(a, c); + } else { const ul_pair & ul = m_lar_solver.m_columns_to_ul_pairs[j]; auto witness = ul.lower_bound_witness(); diff --git a/src/math/lp/nla_intervals.cpp b/src/math/lp/nla_intervals.cpp index c15afcf8b..657c4064d 100755 --- a/src/math/lp/nla_intervals.cpp +++ b/src/math/lp/nla_intervals.cpp @@ -18,66 +18,72 @@ bool intervals::check() { return true; } // create a product of interval signs together with the depencies -intervals::interval intervals::mul_signs_with_deps(int sign, const svector& vars) const { - interval a, b; - m_imanager.set(a, mpq(sign)); +intervals::interval intervals::mul_signs_with_deps(const svector& vars) const { + interval a, b, c; + m_imanager.set(a, mpq(1)); for (lpvar v : vars) { - set_var_interval(v, b); + set_var_interval_signs_with_deps(v, b); + if (m_imanager.is_zero(b)) + return b; + interval_deps deps; + m_imanager.mul(a, b, c, deps); + m_imanager.set(a, c); + m_config.add_deps(a, b, deps, a); } return a; } -bool intervals::check(monomial const& m) { - interval a, b, c, d; - m_imanager.set(a, mpq(1)); - set_var_interval(m.var(), d); - if (m_imanager.lower_is_inf(d) && m_imanager.upper_is_inf(d)) { - return true; - } - for (lpvar v : m.vars()) { - // TBD allow for division to get range of a - // m = a*b*c, where m and b*c are bounded, then interval for a is m/b*c - if (m_imanager.lower_is_inf(a) && m_imanager.upper_is_inf(a)) { - return true; - } - // TBD: deal with powers v^n interval instead of multiplying v*v .. * v - set_var_interval(v, b); - interval_deps deps; - m_imanager.mul(a, b, c, deps); - m_imanager.set(a, c); - m_config.add_deps(a, b, deps, a); - } - if (m_imanager.before(a, d)) { - svector cs; - m_dep_manager.linearize(a.m_upper_dep, cs); - m_dep_manager.linearize(d.m_lower_dep, cs); - for (auto ci : cs) { - (void)ci; - SASSERT(false); - //expl.push_justification(ci); - } - // TBD conflict - return false; - } - if (m_imanager.before(d, a)) { - svector cs; - m_dep_manager.linearize(a.m_lower_dep, cs); - m_dep_manager.linearize(d.m_upper_dep, cs); - for (auto ci : cs) { - (void)ci; - SASSERT(false); //expl.push_justification(ci); - } - // TBD conflict - return false; - } - // could also perform bounds propagation: - // a has tighter lower/upper bound than m.var(), - // -> transfer bound to m.var() - // all but one variable has bound - // -> transfer bound to that variable using division - return true; -} +// bool intervals::check(monomial const& m) { +// interval a, b, c, d; +// m_imanager.set(a, mpq(1)); +// set_var_interval(m.var(), d); +// if (m_imanager.lower_is_inf(d) && m_imanager.upper_is_inf(d)) { +// return true; +// } +// for (lpvar v : m.vars()) { +// // TBD allow for division to get range of a +// // m = a*b*c, where m and b*c are bounded, then interval for a is m/b*c +// if (m_imanager.lower_is_inf(a) && m_imanager.upper_is_inf(a)) { +// return true; +// } +// // TBD: deal with powers v^n interval instead of multiplying v*v .. * v +// set_var_interval(v, b); +// interval_deps deps; +// m_imanager.mul(a, b, c, deps); +// m_imanager.set(a, c); +// m_config.add_deps(a, b, deps, a); +// } +// if (m_imanager.before(a, d)) { +// svector cs; +// m_dep_manager.linearize(a.m_upper_dep, cs); +// m_dep_manager.linearize(d.m_lower_dep, cs); +// for (auto ci : cs) { +// (void)ci; +// SASSERT(false); +// //expl.push_justification(ci); +// } +// // TBD conflict +// return false; +// } +// if (m_imanager.before(d, a)) { +// svector cs; +// m_dep_manager.linearize(a.m_lower_dep, cs); +// m_dep_manager.linearize(d.m_upper_dep, cs); +// for (auto ci : cs) { +// (void)ci; +// SASSERT(false); //expl.push_justification(ci); +// } +// // TBD conflict +// return false; +// } +// // could also perform bounds propagation: +// // a has tighter lower/upper bound than m.var(), +// // -> transfer bound to m.var() +// // all but one variable has bound +// // -> transfer bound to that variable using division +// return true; +// } void intervals::set_var_interval(lpvar v, interval& b) const { lp::constraint_index ci; @@ -87,23 +93,19 @@ void intervals::set_var_interval(lpvar v, interval& b) const { m_config.set_lower(b, val); m_config.set_lower_is_open(b, is_strict); m_config.set_lower_is_inf(b, false); - b.m_lower_dep = mk_dep(ci); } else { m_config.set_lower_is_open(b, true); m_config.set_lower_is_inf(b, true); - b.m_lower_dep = nullptr; } if (ls().has_upper_bound(v, ci, val, is_strict)) { m_config.set_upper(b, val); m_config.set_upper_is_open(b, is_strict); m_config.set_upper_is_inf(b, false); - b.m_upper_dep = mk_dep(ci); } else { m_config.set_upper_is_open(b, true); m_config.set_upper_is_inf(b, true); - b.m_upper_dep = nullptr; } } @@ -165,51 +167,63 @@ intervals::ci_dependency *intervals::mk_dep(lp::constraint_index ci) const { return m_dep_manager.mk_leaf(ci); } -bool intervals::check(lp::lar_term const& t) { - // convert term into factors for improved precision - return true; -} - lp::impq intervals::get_upper_bound_of_monomial(lpvar j) const { const monomial& m = m_core->emons()[j]; - interval a = mul(1, m.vars()); - + interval a = mul(m.vars()); + SASSERT(!m_imanager.upper_is_inf(a)); auto r = lp::impq(a.m_upper); if (a.m_upper_open) r.y = -1; + TRACE("nla_intervals", m_core->print_monomial_with_vars(m, tout) << "upper = " << r << "\n";); return r; } lp::impq intervals::get_lower_bound_of_monomial(lpvar j) const { const monomial& m = m_core->emons()[j]; - interval a = mul(1, m.vars()); - + interval a = mul(m.vars()); + SASSERT(!a.m_lower_inf); auto r = lp::impq(a.m_lower); if (a.m_lower_open) r.y = 1; + TRACE("nla_intervals", m_core->print_monomial_with_vars(m, tout) << "lower = " << r << "\n";); return r; } -intervals::interval intervals::mul(int sign, const svector& vars) const { +intervals::interval intervals::mul(const svector& vars) const { interval a; - m_imanager.set(a, mpq(sign)); + m_imanager.set(a, mpq(1)); for (lpvar j : vars) { interval b, c; set_var_interval(j, b); - m_imanager.mul(a, b, c); - if (m_imanager.is_zero(c)) { - TRACE("nla_intervals", tout << "sign = " << sign << "\nproduct = "; - m_core->print_product_with_vars(vars, tout) << "collapsed to zero\n";); - return c; + if (m_imanager.is_zero(b)) { + return b; } + m_imanager.mul(a, b, c); + m_imanager.set(a, c); + } + return a; +} + +intervals::interval intervals::mul_signs(const svector& vars) const { + interval a; + m_imanager.set(a, mpq(1)); + + for (lpvar j : vars) { + interval b, c; + set_var_interval_signs(j, b); + if (m_imanager.is_zero(b)) { + return b; + } + m_imanager.mul(a, b, c); m_imanager.set(a, c); } return a; } bool intervals::product_has_upper_bound(int sign, const svector& vars) const { - interval a = mul(sign, vars); - return !m_imanager.upper_is_inf(a); + interval a = mul_signs(vars); + SASSERT(sign == 1 || sign == -1); + return sign == 1 ? !m_imanager.upper_is_inf(a) : !m_imanager.lower_is_inf(a); } bool intervals::monomial_has_lower_bound(lpvar j) const { @@ -222,6 +236,18 @@ bool intervals::monomial_has_upper_bound(lpvar j) const { return product_has_upper_bound(1, m.vars()); } lp::lar_solver& intervals::ls() { return m_core->m_lar_solver; } + const lp::lar_solver& intervals::ls() const { return m_core->m_lar_solver; } + +void intervals::get_explanation_of_upper_bound_for_monomial(lpvar j, svector& expl) const { + interval a = mul_signs_with_deps(m_core->emons()[j].vars()); + m_dep_manager.linearize(a.m_upper_dep, expl); } +void intervals::get_explanation_of_lower_bound_for_monomial(lpvar j, svector& expl) const{ + interval a = mul_signs_with_deps(m_core->emons()[j].vars()); + m_dep_manager.linearize(a.m_lower_dep, expl); + // return m_intervals.get_explanation_of_lower_bound_for_monomial(j, expl ) +} +} +// instantiate the template template class interval_manager; diff --git a/src/math/lp/nla_intervals.h b/src/math/lp/nla_intervals.h index be06c1733..0766181d8 100755 --- a/src/math/lp/nla_intervals.h +++ b/src/math/lp/nla_intervals.h @@ -78,7 +78,7 @@ namespace nla { }; - void add_deps(interval const& a, interval const& b, interval_deps const& deps, interval& i) { + void add_deps(interval const& a, interval const& b, interval_deps const& deps, interval& i) const { ci_dependency* lo = mk_dependency(a, b, deps.m_lower_deps); ci_dependency* hi = mk_dependency(a, b, deps.m_upper_deps); i.m_lower_dep = lo; @@ -119,7 +119,7 @@ namespace nla { im_config(numeral_manager & m, ci_dependency_manager& d):m_manager(m), m_dep_manager(d) {} private: - ci_dependency* mk_dependency(interval const& a, interval const& b, bound_deps bd) { + ci_dependency* mk_dependency(interval const& a, interval const& b, bound_deps bd) const { ci_dependency* dep = nullptr; if (dep_in_lower1(bd)) { dep = m_dep_manager.mk_join(dep, a.m_lower_dep); @@ -174,7 +174,10 @@ namespace nla { bool product_has_upper_bound(int sign, const svector&) const; lp::impq get_upper_bound_of_monomial(lpvar j) const; lp::impq get_lower_bound_of_monomial(lpvar j) const; - interval mul(int sign, const svector&) const; - interval mul_signs_with_deps(int sign, const svector&) const; + interval mul(const svector&) const; + interval mul_signs(const svector&) const; + interval mul_signs_with_deps(const svector&) const; + void get_explanation_of_upper_bound_for_monomial(lpvar j, svector& expl) const; + void get_explanation_of_lower_bound_for_monomial(lpvar j, svector& expl) const; }; } // end of namespace nla diff --git a/src/math/lp/nla_solver.cpp b/src/math/lp/nla_solver.cpp index 3f58ead75..c9dbd8b06 100644 --- a/src/math/lp/nla_solver.cpp +++ b/src/math/lp/nla_solver.cpp @@ -77,5 +77,10 @@ bool solver::monomial_has_lower_bound(lpvar j) const { bool solver::monomial_has_upper_bound(lpvar j) const { return m_intervals.monomial_has_upper_bound(j); } - +void solver::get_explanation_of_upper_bound_for_monomial(lpvar j, svector& expl) const { + m_intervals.get_explanation_of_upper_bound_for_monomial(j, expl); +} +void solver::get_explanation_of_lower_bound_for_monomial(lpvar j, svector& expl) const{ + m_intervals.get_explanation_of_lower_bound_for_monomial(j, expl); +} } diff --git a/src/math/lp/nla_solver.h b/src/math/lp/nla_solver.h index 8f1fd7622..173ff4dfe 100644 --- a/src/math/lp/nla_solver.h +++ b/src/math/lp/nla_solver.h @@ -51,5 +51,7 @@ public: lp::impq get_upper_bound(lpvar j) const; bool monomial_has_lower_bound(lpvar j) const; bool monomial_has_upper_bound(lpvar j) const; + void get_explanation_of_upper_bound_for_monomial(lpvar j, svector&) const; + void get_explanation_of_lower_bound_for_monomial(lpvar j, svector&) const; }; }