get bounds from interval multiplication

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

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<lp::constraint_index> 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<lp::constraint_index> 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();

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<lpvar>& vars) const {
+intervals::interval intervals::mul_signs_with_deps(const svector<lpvar>& vars) const {
-    interval a, b;
+    interval a, b, c;
-    m_imanager.set(a, mpq(sign));
+    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) {
+// bool intervals::check(monomial const& m) {
-    interval a, b, c, d;
+//     interval a, b, c, d;
-    m_imanager.set(a, mpq(1));
+//     m_imanager.set(a, mpq(1));
-    set_var_interval(m.var(), d);
+//     set_var_interval(m.var(), d);
-    if (m_imanager.lower_is_inf(d) && m_imanager.upper_is_inf(d)) {
+//     if (m_imanager.lower_is_inf(d) && m_imanager.upper_is_inf(d)) {
-        return true;
+//         return true;
-    }
+//     }
-    for (lpvar v : m.vars()) {
+//     for (lpvar v : m.vars()) {
-        // TBD allow for division to get range of a
+//         // 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
+//         // 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)) {
+//         if (m_imanager.lower_is_inf(a) && m_imanager.upper_is_inf(a)) {
-            return true;
+//             return true;
-        }
+//         }
-        // TBD: deal with powers v^n interval instead of multiplying v*v .. * v
+//         // TBD: deal with powers v^n interval instead of multiplying v*v .. * v
-        set_var_interval(v, b);
+//         set_var_interval(v, b);
-        interval_deps deps;
+//         interval_deps deps;
-        m_imanager.mul(a, b, c, deps);
+//         m_imanager.mul(a, b, c, deps);
-        m_imanager.set(a, c);
+//         m_imanager.set(a, c);
-        m_config.add_deps(a, b, deps, a);
+//         m_config.add_deps(a, b, deps, a);
-    }
+//     }
-    if (m_imanager.before(a, d)) {
+//     if (m_imanager.before(a, d)) {
-        svector<lp::constraint_index> cs;
+//         svector<lp::constraint_index> cs;
-        m_dep_manager.linearize(a.m_upper_dep, cs);
+//         m_dep_manager.linearize(a.m_upper_dep, cs);
-        m_dep_manager.linearize(d.m_lower_dep, cs);
+//         m_dep_manager.linearize(d.m_lower_dep, cs);
-        for (auto ci : cs) {
+//         for (auto ci : cs) {
-            (void)ci;
+//             (void)ci;
-            SASSERT(false);
+//             SASSERT(false);
-            //expl.push_justification(ci);
+//             //expl.push_justification(ci);
-        }
+//         }
-        // TBD conflict
+//         // TBD conflict
-        return false;
+//         return false;
-    }
+//     }
-    if (m_imanager.before(d, a)) {
+//     if (m_imanager.before(d, a)) {
-        svector<lp::constraint_index> cs;
+//         svector<lp::constraint_index> cs;
-        m_dep_manager.linearize(a.m_lower_dep, cs);
+//         m_dep_manager.linearize(a.m_lower_dep, cs);
-        m_dep_manager.linearize(d.m_upper_dep, cs);
+//         m_dep_manager.linearize(d.m_upper_dep, cs);
-        for (auto ci : cs) {
+//         for (auto ci : cs) {
-            (void)ci;
+//             (void)ci;
-            SASSERT(false); //expl.push_justification(ci);
+//             SASSERT(false); //expl.push_justification(ci);
-        }
+//         }
-        // TBD conflict
+//         // TBD conflict
-        return false;
+//         return false;
-    }
+//     }
-    // could also perform bounds propagation:
+//     // could also perform bounds propagation:
-    // a has tighter lower/upper bound than m.var(), 
+//     // a has tighter lower/upper bound than m.var(), 
-    // -> transfer bound to m.var()
+//     // -> transfer bound to m.var()
-    // all but one variable has bound
+//     // all but one variable has bound
-    // -> transfer bound to that variable using division
+//     // -> transfer bound to that variable using division
-    return true;
+//     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<lpvar>& vars) const {
+intervals::interval intervals::mul(const svector<lpvar>& 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(b)) {
-        if (m_imanager.is_zero(c)) {
+            return b;
-            TRACE("nla_intervals", tout << "sign = " << sign << "\nproduct = ";
-                  m_core->print_product_with_vars(vars, tout) << "collapsed to zero\n";);
-            return c;
         }
+        m_imanager.mul(a, b, c);
+        m_imanager.set(a, c);
+    }
+    return a;
+}
+
+intervals::interval intervals::mul_signs(const svector<lpvar>& 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<lpvar>& vars) const {
-    interval a = mul(sign, vars);
+    interval a = mul_signs(vars);
-    return !m_imanager.upper_is_inf(a);
+    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<lp::constraint_index>& 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<lp::constraint_index>& 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<nla::intervals::im_config>;

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<lpvar>&) 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<lpvar>&) const;
+        interval mul(const svector<lpvar>&) const;
-        interval mul_signs_with_deps(int sign, const svector<lpvar>&) const;
+        interval mul_signs(const svector<lpvar>&) const;
+        interval mul_signs_with_deps(const svector<lpvar>&) const;
+        void get_explanation_of_upper_bound_for_monomial(lpvar j, svector<lp::constraint_index>& expl) const;
+        void get_explanation_of_lower_bound_for_monomial(lpvar j, svector<lp::constraint_index>& expl) const;
     };
 } // end of namespace nla

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<lp::constraint_index>& 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<lp::constraint_index>& expl) const{
+    m_intervals.get_explanation_of_lower_bound_for_monomial(j, expl);
+}
 }

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<lp::constraint_index>&) const;
+    void get_explanation_of_lower_bound_for_monomial(lpvar j, svector<lp::constraint_index>&) const;
 };
 }