inverse approx

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/math/polysat/fixplex.h

@@ -47,12 +47,17 @@ namespace polysat {
      */
     template<typename Numeral>
     struct interval {
+        bool emp { false };
         Numeral lo { 0 };
         Numeral hi { 0 };
         interval() {}
         interval(Numeral const& l, Numeral const& h): lo(l), hi(h) {}
         static interval free() { return interval(0, 0); }
+        static interval empty() { interval i(0, 0); i.emp = true; return i; }
         bool is_free() const { return lo == hi; }
+        bool is_empty() const { return emp; }
+        bool contains(Numeral const& n) const;
+        interval operator&(interval const& other) const;
         interval operator+(interval const& other) const;
         interval operator-(interval const& other) const;
         interval operator-() const;
@@ -60,7 +65,11 @@ namespace polysat {
         interval operator+(Numeral const& n) const { return interval(lo + n, hi + n); }
         interval operator-(Numeral const& n) const { return interval(lo - n, hi - n); }
         interval& operator+=(interval const& other) { *this = *this + other; return *this; }
-        std::ostream& display(std::ostream& out) const { return out << "[" << pp(lo) << ", " << pp(hi) << "["; }
+        std::ostream& display(std::ostream& out) const { 
+            if (is_empty()) return out << "empty"; 
+            if (is_free()) return out << "free";
+            return out << "[" << pp(lo) << ", " << pp(hi) << "["; 
+        }
     };
 
     template<typename Numeral>
@@ -109,6 +118,10 @@ namespace polysat {
                 m_base2row(0),
                 m_is_base(false)
             {}
+            var_info& operator&=(interval<numeral> const& range) {
+                interval<numeral>::operator=(range);
+                return *this;
+            }
         };
 
         struct row_info {
@@ -203,7 +216,7 @@ namespace polysat {
         void fixed_var_eh(row const& r, var_t x);
         void eq_eh(var_t x, var_t y, row const& r1, row const& r2);
         void propagate_bounds(row const& r);
-        void new_bound(row const& r, var_t x, numeral const& l, numeral const& h);
+        void new_bound(row const& r, var_t x, interval<numeral> const& range);
         void pivot(var_t x_i, var_t x_j, numeral const& b, numeral const& value);
         numeral value2delta(var_t v, numeral const& new_value) const;
         void update_value(var_t v, numeral const& delta);
@@ -211,8 +224,8 @@ namespace polysat {
         bool has_minimal_trailing_zeros(var_t y, numeral const& b);
         var_t select_pivot_core(var_t x, numeral const& new_value, numeral& out_b);
         bool in_bounds(var_t v) const { return in_bounds(v, value(v)); }
-        bool in_bounds(var_t v, numeral const& b) const { return in_bounds(b, lo(v), hi(v)); }
+        bool in_bounds(var_t v, numeral const& b) const { return in_bounds(b, m_vars[v]); }
-        bool in_bounds(numeral const& val, numeral const& lo, numeral const& hi) const;        
+        bool in_bounds(numeral const& val, interval<numeral> const& range) const { return range.contains(val); }
         bool is_free(var_t v) const { return lo(v) == hi(v); }
         bool is_non_free(var_t v) const { return !is_free(v); }
         bool is_fixed(var_t v) const { return lo(v) + 1 == hi(v); }

src/math/polysat/fixplex_def.h

@@ -23,8 +23,24 @@ Author:
 
 namespace polysat {
 
+    template<typename Numeral>
+    bool interval<Numeral>::contains(Numeral const& n) const {
+        if (is_empty())
+            return false;
+        if (is_free())
+            return true;
+        if (lo < hi)
+            return lo <= n && n < hi;
+        else
+            return lo <= n || n < hi;
+    }
+
     template<typename Numeral>
     interval<Numeral> interval<Numeral>::operator+(interval<Numeral> const& other) const {
+        if (is_empty())
+            return *this;
+        if (other.is_empty())
+            return other;
         if (is_free())
             return *this;
         if (other.is_free())
@@ -42,6 +58,8 @@ namespace polysat {
 
     template<typename Numeral>
     interval<Numeral> interval<Numeral>::operator-() const {
+        if (is_empty())
+            return *this;
         if (is_free())
             return *this;
         return interval(1 - hi, 1 - lo);
@@ -49,6 +67,8 @@ namespace polysat {
 
     template<typename Numeral>
     interval<Numeral> interval<Numeral>::operator*(Numeral const& n) const {
+        if (is_empty())
+            return *this;
         if (n == 0)
             return interval(0, 1);
         if (n == 1)
@@ -71,6 +91,26 @@ namespace polysat {
         }
     }
 
+    template<typename Numeral>
+    interval<Numeral> interval<Numeral>::operator&(interval const& other) const {
+        Numeral l, h;
+        if (is_free() || other.is_empty())
+            return other;
+        if (other.is_free() || is_empty())
+            return *this;
+        if (contains(other.lo))
+            l = other.lo;
+        else if (other.contains(lo))
+            l = lo;
+        else interval::empty();
+        if (contains(other.hi - 1))
+            h = other.hi;
+        else if (other.contains(hi - 1))
+            h = hi;
+        return interval(l, h);
+    }
+
+
     template<typename Ext>
     fixplex<Ext>::~fixplex() {
         reset();
@@ -228,14 +268,6 @@ namespace polysat {
         
     }
 
-    template<typename Ext>
-    bool fixplex<Ext>::in_bounds(numeral const& val, numeral const& lo, numeral const& hi) const {
-        if (lo == hi)
-            return true;
-        if (lo < hi)
-            return lo <= val && val < hi;
-        return val < hi || lo <= val;
-    }
 
     /**
      * Attempt to improve assigment to make x feasible.
@@ -267,6 +299,8 @@ namespace polysat {
         
         pivot(x, y, b, new_value);
 
+        std::cout << "pivot v" << y << " := " << new_value << "\n";
+
         // get_offset_eqs(row(base2row(y)));
 
         return l_true;
@@ -304,7 +338,7 @@ namespace polysat {
             if (!has_minimal_trailing_zeros(y, b))
                 continue;
             numeral new_y_value = (row_value - b*value(y) - a*new_value)/b;
-            bool _in_bounds = in_bounds(new_y_value, lo(y), hi(y));
+            bool _in_bounds = in_bounds(y, new_y_value);
             if (!_in_bounds) {
                 if (lo(y) - new_y_value < new_y_value - hi(y))
                     delta_y = new_y_value - lo(y);
@@ -520,6 +554,7 @@ namespace polysat {
         numeral old_value_y = yI.m_value;
         row_x.m_base = y;
         row_x.m_value = row_x.m_value - b*old_value_y + a*new_value;
+        std::cout << "row value " << row_x.m_value << "\n";
         row_x.m_base_coeff = b;
         yI.m_base2row = rx;
         yI.m_is_base = true;
@@ -663,7 +698,15 @@ namespace polysat {
     void fixplex<Ext>::set_base_value(var_t x) {
         SASSERT(is_base(x));
         row r(base2row(x));
-        m_vars[x].m_value = 0 - (row2value(r) / row2base_coeff(r));
+        auto c = row2base_coeff(r);
+        if (c == 1)
+            m_vars[x].m_value = 0 - row2value(r);
+        else if (c + 1 == 0)
+            m_vars[x].m_value = row2value(r);
+        else {
+            // TBD: compute "inverse" of c?
+            m_vars[x].m_value = 0 - row2value(r) / c;           
+        }
         bool was_integral = row2integral(r);
         m_rows[r.id()].m_integral = is_solved(r);
         if (was_integral && !row2integral(r))
@@ -808,29 +851,25 @@ namespace polysat {
 
         if (free_v != null_var) {
             range = (-range) * free_c;
-            if (range.is_free())
+            new_bound(r, free_v, range);
-                return;
-            new_bound(r, free_v, range.lo, range.hi);
             SASSERT(in_bounds(free_v));
             return;
         }
         for (auto const& e : M.row_entries(r)) {
             var_t v = e.var();
             SASSERT(!is_free(v));
-            numeral const& c = e.coeff();
+            auto range1 = range - m_vars[v] * e.coeff();
-            auto range1 = range - m_vars[v] * c;
+            new_bound(r, v, range1);
-            new_bound(r, v, range1.lo, range1.hi);
             // SASSERT(in_bounds(v));
         }
     }
 
     template<typename Ext>
-    void fixplex<Ext>::new_bound(row const& r, var_t x, numeral const& l, numeral const& h) {
+    void fixplex<Ext>::new_bound(row const& r, var_t x, interval<numeral> const& range) {
-        // TBD: 
+        if (range.is_free())
-        std::cout << is_free(x) << " " << l << " " << lo(x) << " " << hi(x) << " " << h << "\n";
+            return;
-        if (!is_free(x) && l <= lo(x) && hi(x) <= h)
+        m_vars[x] &= range;
-            return;            
+        IF_VERBOSE(0, verbose_stream() << "new-bound v" << x << " " << m_vars[x] << "\n");
-        IF_VERBOSE(0, verbose_stream() << "new-bound v" << x << " [" << pp(l) << "," << pp(h) << "[\n");
     }
 
     template<typename Ext>    
@@ -838,7 +877,7 @@ namespace polysat {
         M.display(out);
         for (unsigned i = 0; i < m_vars.size(); ++i) {
             var_info const& vi = m_vars[i];
-            out << "v" << i << " " << pp(value(i)) << " [" << pp(lo(i)) << ", " << pp(hi(i)) << "[ ";
+            out << "v" << i << " " << pp(value(i)) << " " << vi << " ";
             if (vi.m_is_base) out << "b:" << vi.m_base2row << " ";
             out << "\n";
         }
@@ -853,7 +892,7 @@ namespace polysat {
                 out << pp(e.coeff()) << " * ";
             out << "v" << v << " ";
             if (values) 
-                out << pp(value(v)) << " [" << pp(lo(v)) << ", " << pp(hi(v)) << "[ ";            
+                out << pp(value(v)) << " " << m_vars[v];
         }
         return out << "\n";
     }