throttle intervals

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

src/math/lp/CMakeLists.txt

@@ -28,6 +28,7 @@ z3_add_component(lp
     nla_basics_lemmas.cpp
     nla_common.cpp
     nla_core.cpp
+    nla_intervals.cpp
     nla_monotone_lemmas.cpp
     nla_order_lemmas.cpp
     nla_solver.cpp

src/math/lp/nla_basics_lemmas.cpp

@@ -250,8 +250,8 @@ bool basics::basic_lemma_for_mon_zero(const monomial& rm, const factorization& f
 bool basics::basic_lemma(bool derived) {
     if (basic_sign_lemma(derived))
         return true;
-    if (derived)
-        return false; // c().m_intervals.get_lemmas();
+    if (derived) 
+        return c().m_intervals.get_lemmas();
     const auto& mon_inds_to_ref = c().m_to_refine;
     TRACE("nla_solver", tout << "mon_inds_to_ref = "; print_vector(mon_inds_to_ref, tout););
     unsigned start = c().random();

src/math/lp/nla_core.cpp

@@ -25,6 +25,7 @@ core::core(lp::lar_solver& s, reslimit & lim) :
     m_evars(),
     m_lar_solver(s),
     m_tangents(this),
+    m_intervals(this, lim),
     m_basics(this),
     m_order(this),
     m_monotone(this),
@@ -136,11 +137,13 @@ void core::add_monomial(lpvar v, unsigned sz, lpvar const* vs) {
 void core::push() {
     TRACE("nla_solver",);
     m_emons.push();
+    m_intervals.push();
 }
 
      
 void core::pop(unsigned n) {
     TRACE("nla_solver", tout << "n = " << n << "\n";);
+    m_intervals.pop(n);
     m_emons.pop(n);
     SASSERT(elists_are_consistent(false));
 }
@@ -955,6 +958,7 @@ void core::clear() {
 void core::init_search() {
     clear();
     init_vars_equivalence();
+    m_intervals.init();
     SASSERT(elists_are_consistent(false));
 }
 

src/math/lp/nla_core.h

@@ -27,7 +27,7 @@
 #include "math/lp/nla_monotone_lemmas.h"
 #include "math/lp/emonomials.h"
 #include "math/lp/nla_settings.h"
-// #include "math/lp/nla_intervals.h"
+#include "math/lp/nla_intervals.h"
 namespace nla {
 
 template <typename A, typename B>
@@ -83,6 +83,7 @@ public:
     vector<lemma> *          m_lemma_vec;
     svector<lpvar>           m_to_refine;
     tangents                 m_tangents;
+    intervals                m_intervals;
     basics                   m_basics;
     order                    m_order;
     monotone                 m_monotone;

src/math/lp/nla_intervals.cpp

@@ -0,0 +1,298 @@
+#include "math/lp/nla_core.h"
+#include "math/interval/interval_def.h"
+#include "math/lp/nla_intervals.h"
+
+namespace nla {
+
+bool intervals::get_lemmas() {
+    m_region.reset();
+    bool ret = false;
+    for (auto const& k : c().m_to_refine) {
+        if (get_lemma(c().emons()[k])) {
+            ret = true;
+        }
+        if (c().done())
+            break;
+    }
+    return ret;
+}
+// create a product of interval signs together with the depencies
+intervals::interval intervals::mul_signs_with_deps(const svector<lpvar>& vars) const {
+    interval a, b, c;
+    m_imanager.set(a, mpq(1));
+    for (lpvar v : vars) {
+        set_var_interval_signs_with_deps(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.is_zero(a))
+            return a;
+    }
+    return a;    
+}
+
+void intervals::get_lemma_for_zero_interval(monomial const& m) {
+    if (val(m).is_zero()) return;
+    interval signs_a = mul_signs_with_deps(m.vars());
+    add_empty_lemma();
+    svector<lp::constraint_index> expl;
+    m_dep_manager.linearize(signs_a.m_lower_dep, expl); 
+    _().current_expl().add_expl(expl);
+    mk_ineq(m.var(), llc::EQ);
+    TRACE("nla_solver", _().print_lemma(tout); );
+}
+
+bool intervals::get_lemma_for_lower(const monomial& m, const interval& a) {
+    if (m_vars_pushed_up[m.var()] > 10)
+        return false;
+    lp::impq lb(rational(a.m_lower));
+    if (m_config.lower_is_open(a))
+        lb.y = 1;
+
+    lp::impq v(val(m.var()));        
+    if (v < lb) {
+        m_vars_pushed_up[m.var()] = m_vars_pushed_up[m.var()] + 1;
+        interval signs_a = mul_signs_with_deps(m.vars());
+        add_empty_lemma();
+        svector<lp::constraint_index> expl;
+        m_dep_manager.linearize(signs_a.m_lower_dep, expl); 
+        _().current_expl().add_expl(expl);
+        llc cmp = m_config.lower_is_open(a)? llc::GT: llc::GE;
+        mk_ineq(m.var(), cmp,  lb.x);
+        TRACE("nla_solver", _().print_lemma(tout); );
+        return true;
+    }
+    return false;
+}
+
+bool intervals::get_lemma_for_upper(const monomial& m, const interval& a) {
+    if (m_vars_pushed_down[m.var()] > 10)
+        return false;
+    lp::impq ub(rational(a.m_upper));
+    if (m_config.upper_is_open(a))
+        ub.y = 1;
+    lp::impq v(val(m.var()));
+    if (v > ub) {
+        m_vars_pushed_down[m.var()] = m_vars_pushed_down[m.var()] + 1;
+        interval signs_a = mul_signs_with_deps(m.vars());
+        add_empty_lemma();
+        svector<lp::constraint_index> expl;
+        m_dep_manager.linearize(signs_a.m_upper_dep, expl); 
+        _().current_expl().add_expl(expl);
+        llc cmp = m_config.upper_is_open(a)? llc::LT: llc::LE;
+        mk_ineq(m.var(), cmp,  ub.x);
+        TRACE("nla_solver", _().print_lemma(tout); );
+        return true;
+    }
+    return false;
+}
+bool intervals::get_lemma(monomial const& m) {
+    interval  b, c, d;
+    interval a = mul(m.vars());
+    if (m_config.is_zero(a)) {
+        get_lemma_for_zero_interval(m);
+        return true;
+    }
+    if (!m_imanager.lower_is_inf(a)) {
+        return get_lemma_for_lower(m, a);
+    }
+
+    if (!m_imanager.upper_is_inf(a)) {
+        return get_lemma_for_upper(m, a);
+    }
+    return false;
+} 
+void intervals::set_var_interval(lpvar v, interval& b) const {
+    lp::constraint_index ci;
+    rational val;
+    bool is_strict;
+    if (ls().has_lower_bound(v, ci, val, is_strict)) {            
+        m_config.set_lower(b, val);
+        m_config.set_lower_is_open(b, is_strict);
+        m_config.set_lower_is_inf(b, false);
+    }
+    else {
+        m_config.set_lower_is_open(b, true);
+        m_config.set_lower_is_inf(b, true);
+    }
+    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);
+    }
+    else {
+        m_config.set_upper_is_open(b, true);
+        m_config.set_upper_is_inf(b, true);
+    }
+}
+
+rational sign(const rational& v) { return v.is_zero()? v : (rational(v.is_pos()? 1 : -1)); }
+
+void intervals::set_var_interval_signs(lpvar v, interval& b) const {
+    lp::constraint_index ci;
+    rational val;
+    bool is_strict;
+    if (ls().has_lower_bound(v, ci, val, is_strict)) {            
+        m_config.set_lower(b, sign(val));
+        m_config.set_lower_is_open(b, is_strict);
+        m_config.set_lower_is_inf(b, false);
+    }
+    else {
+        m_config.set_lower_is_open(b, true);
+        m_config.set_lower_is_inf(b, true);
+    }
+    if (ls().has_upper_bound(v, ci, val, is_strict)) {
+        m_config.set_upper(b, sign(val));
+        m_config.set_upper_is_open(b, is_strict);
+        m_config.set_upper_is_inf(b, false);
+    }
+    else {
+        m_config.set_upper_is_open(b, true);
+        m_config.set_upper_is_inf(b, true);
+    }
+}
+
+void intervals::set_var_interval_signs_with_deps(lpvar v, interval& b) const {
+    lp::constraint_index ci;
+    rational val;
+    bool is_strict;
+    if (ls().has_lower_bound(v, ci, val, is_strict)) {            
+        m_config.set_lower(b, sign(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, sign(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;
+    }
+}
+    
+intervals::ci_dependency *intervals::mk_dep(lp::constraint_index ci) const {
+    return m_dep_manager.mk_leaf(ci);
+}
+    
+lp::impq intervals::get_upper_bound_of_monomial(lpvar j) const {
+    const monomial& m = m_core->emons()[j];
+    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_detail", 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(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_detail", m_core->print_monomial_with_vars(m, tout) << "lower = " << r << "\n";);
+    return r;
+}
+
+intervals::interval intervals::mul(const svector<lpvar>& vars) const {
+    interval a;
+    m_imanager.set(a, mpq(1));
+    
+    for (lpvar j : vars) {
+        interval b, c;
+        set_var_interval(j, b);
+        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<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_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 {
+    const monomial& m = m_core->emons()[j];
+    return product_has_upper_bound(-1, m.vars());
+}
+
+bool intervals::monomial_has_upper_bound(lpvar j) const {
+    const monomial& m = m_core->emons()[j];
+    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; }
+
+std::ostream& intervals::print_explanations(const svector<lp::constraint_index> &expl , std::ostream& out) const {
+    out << "interv expl:\n ";
+    for (auto ci : expl)
+        m_core->m_lar_solver.print_constraint_indices_only(ci, out);
+    return out;
+}
+
+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);
+    TRACE("nla_intervals", print_explanations(expl, tout););
+}
+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);    
+    TRACE("nla_intervals", print_explanations(expl, tout););
+    //    return m_intervals.get_explanation_of_lower_bound_for_monomial(j, expl ) 
+}
+void intervals::push() {
+    m_vars_pushed_up.push();
+    m_vars_pushed_down.push();
+}
+void intervals::pop(unsigned k) {
+    m_vars_pushed_up.pop(k);
+    m_vars_pushed_down.pop(k);
+}
+
+void intervals::init() {
+    SASSERT(m_vars_pushed_down.size() == m_vars_pushed_up.size());
+    unsigned n = c().m_lar_solver.number_of_vars();
+    while (m_vars_pushed_up.size() < n) {
+        m_vars_pushed_up.push_back(0);
+        m_vars_pushed_down.push_back(0);
+    }
+}
+}
+
+// instantiate the template
+template class interval_manager<nla::intervals::im_config>;

src/math/lp/nla_intervals.h

@@ -0,0 +1,189 @@
+/*++
+  Copyright (c) 2017 Microsoft Corporation
+
+  Module Name:
+
+  <name>
+
+  Abstract:
+
+  <abstract>
+
+  Author:
+  Nikolaj Bjorner (nbjorner)
+  Lev Nachmanson (levnach)
+
+  Revision History:
+
+
+  --*/
+#pragma once
+#include "util/dependency.h"
+#include "util/small_object_allocator.h"
+#include "math/lp/nla_common.h"
+#include "math/lp/lar_solver.h"
+#include "math/interval/interval.h"
+
+
+namespace nla {
+class core;
+
+class intervals : common {
+    // fields to throttle the propagation on intervals
+    lp::stacked_vector<unsigned> m_vars_pushed_up;
+    lp::stacked_vector<unsigned> m_vars_pushed_down;
+
+    class ci_value_manager {
+    public:
+        void inc_ref(lp::constraint_index const & v) {
+        }
+
+        void dec_ref(lp::constraint_index const & v) {
+        }
+    };
+
+    struct ci_dependency_config {
+        typedef ci_value_manager        value_manager;
+        typedef small_object_allocator  allocator;
+        static const bool ref_count = false;
+        typedef lp::constraint_index value;
+    };
+        
+    typedef dependency_manager<ci_dependency_config> ci_dependency_manager;
+
+    typedef ci_dependency_manager::dependency ci_dependency;
+
+    class im_config {
+        unsynch_mpq_manager&        m_manager;
+        ci_dependency_manager&       m_dep_manager;
+
+    public:
+        typedef unsynch_mpq_manager numeral_manager;
+            
+
+        struct interval {
+            interval():
+                m_lower(), m_upper(), 
+                m_lower_open(1), m_upper_open(1), 
+                m_lower_inf(1), m_upper_inf(1), 
+                m_lower_dep(nullptr), m_upper_dep(nullptr) {}
+            mpq   m_lower;
+            mpq   m_upper;
+            unsigned  m_lower_open:1;
+            unsigned  m_upper_open:1;
+            unsigned  m_lower_inf:1;
+            unsigned  m_upper_inf:1;
+            ci_dependency * m_lower_dep; // justification for the lower bound
+            ci_dependency * m_upper_dep; // justification for the upper bound
+        };
+
+
+        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;
+            i.m_upper_dep = hi;
+        }
+            
+        // Should be NOOPs for precise mpq types.
+        // For imprecise types (e.g., floats) it should set the rounding mode.
+        void round_to_minus_inf() {}
+        void round_to_plus_inf() {}
+        void set_rounding(bool to_plus_inf) {}
+            
+        // Getters
+        mpq const & lower(interval const & a) const { return a.m_lower; }
+        mpq const & upper(interval const & a) const { return a.m_upper; }
+        mpq & lower(interval & a) { return a.m_lower; }
+        mpq & upper(interval & a) { return a.m_upper; }
+        bool lower_is_open(interval const & a) const { return a.m_lower_open; }
+        bool upper_is_open(interval const & a) const { return a.m_upper_open; }
+        bool lower_is_inf(interval const & a) const { return a.m_lower_inf; }
+        bool upper_is_inf(interval const & a) const { return a.m_upper_inf; }
+        bool is_zero(interval const & a) const {
+            return unsynch_mpq_manager::is_zero(a.m_lower)
+                && unsynch_mpq_manager::is_zero(a.m_upper); }
+            
+        // Setters
+        void set_lower(interval & a, mpq const & n) const { m_manager.set(a.m_lower, n); }
+        void set_upper(interval & a, mpq const & n) const { m_manager.set(a.m_upper, n); }
+        void set_lower(interval & a, rational const & n) const { set_lower(a, n.to_mpq()); }
+        void set_upper(interval & a, rational const & n) const { set_upper(a, n.to_mpq()); }
+        void set_lower_is_open(interval & a, bool v) const { a.m_lower_open = v; }
+        void set_upper_is_open(interval & a, bool v) const { a.m_upper_open = v; }
+        void set_lower_is_inf(interval & a, bool v) const { a.m_lower_inf = v; }
+        void set_upper_is_inf(interval & a, bool v) const { a.m_upper_inf = v; }
+            
+        // Reference to numeral manager
+        numeral_manager & m() const { return m_manager; }
+            
+        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) const {
+            ci_dependency* dep = nullptr;
+            if (dep_in_lower1(bd)) {
+                dep = m_dep_manager.mk_join(dep, a.m_lower_dep);
+            }
+            if (dep_in_lower2(bd)) {
+                dep = m_dep_manager.mk_join(dep, b.m_lower_dep);
+            }
+            if (dep_in_upper1(bd)) {
+                dep = m_dep_manager.mk_join(dep, a.m_upper_dep);
+            }
+            if (dep_in_upper2(bd)) {
+                dep = m_dep_manager.mk_join(dep, b.m_upper_dep);
+            }
+            return dep;
+        }
+    };
+
+    small_object_allocator      m_alloc;
+    ci_value_manager            m_val_manager;
+    unsynch_mpq_manager         m_num_manager;
+    mutable ci_dependency_manager       m_dep_manager;
+    im_config                   m_config;
+    mutable interval_manager<im_config> m_imanager;
+    region                      m_region;
+
+    typedef interval_manager<im_config>::interval interval;
+
+    void set_var_interval(lpvar v, interval & b) const;
+
+    void set_var_interval_signs(lpvar v, interval & b) const;
+
+    void set_var_interval_signs_with_deps(lpvar v, interval & b) const;
+
+    ci_dependency* mk_dep(lp::constraint_index ci) const;
+
+    lp::lar_solver& ls();
+    const lp::lar_solver& ls() const;
+public:
+    intervals(core* c, reslimit& lim) :
+        common(c),
+        m_alloc("intervals"),
+        m_dep_manager(m_val_manager, m_alloc),
+        m_config(m_num_manager, m_dep_manager),
+        m_imanager(lim, im_config(m_num_manager, m_dep_manager))
+    {}
+    bool get_lemmas();
+    bool get_lemma(monomial const& m);
+    void get_lemma_for_zero_interval(monomial const& m);
+    bool get_lemma_for_lower(monomial const& m, const interval& );
+    bool get_lemma_for_upper(monomial const& m, const interval &);
+    bool monomial_has_lower_bound(lpvar j) const;
+    bool monomial_has_upper_bound(lpvar j) const;
+    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(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;
+    std::ostream& print_explanations(const svector<lp::constraint_index> &, std::ostream&) const;
+    void push();
+    void pop(unsigned k);
+    void init();
+    
+};
+} // end of namespace nla

src/math/lp/nla_solver.cpp

@@ -24,7 +24,7 @@
 #include "math/lp/var_eqs.h"
 #include "math/lp/factorization.h"
 #include "math/lp/nla_solver.h"
-//#include "math/lp/nla_intervals.h"
+#include "math/lp/nla_intervals.h"
 namespace nla {
 
 void solver::add_monomial(lpvar v, unsigned sz, lpvar const* vs) {

src/math/lp/nla_solver.h

@@ -26,7 +26,7 @@ Revision History:
 #include "math/lp/lar_solver.h"
 #include "math/lp/monomial.h"
 #include "math/lp/nla_core.h"
-// #include "math/lp/nla_intervals.h"
+#include "math/lp/nla_intervals.h"
 
 namespace nla {
 

src/math/lp/nla_tangent_lemmas.h

@@ -23,28 +23,28 @@
 #include "math/lp/nla_common.h"
 
 namespace nla {
-    class core;
+class core;
 
-    struct point {
-        rational x;
-        rational y;
-        point(const rational& a, const rational& b): x(a), y(b) {}
-        point() {}
-        inline point& operator*=(rational a) {
-            x *= a;
-            y *= a;
-            return *this;
-        } 
-        inline point operator+(const point& b) const {
-            return point(x + b.x, y + b.y);
-        } 
+struct point {
+    rational x;
+    rational y;
+    point(const rational& a, const rational& b): x(a), y(b) {}
+    point() {}
+    inline point& operator*=(rational a) {
+        x *= a;
+        y *= a;
+        return *this;
+    } 
+    inline point operator+(const point& b) const {
+        return point(x + b.x, y + b.y);
+    } 
 
-        inline point operator-(const point& b) const {
-            return point(x - b.x, y - b.y);
-        }
-    };
+    inline point operator-(const point& b) const {
+        return point(x - b.x, y - b.y);
+    }
+};
 
-    inline std::ostream& operator<<(std::ostream& out, point const& a) { return out << "(" << a.x <<  ", " << a.y << ")"; }
+inline std::ostream& operator<<(std::ostream& out, point const& a) { return out << "(" << a.x <<  ", " << a.y << ")"; }
 
 
 class tangents : common {