prepare for tuned viable sets

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

src/math/dd/dd_bdd.cpp

@@ -1193,5 +1193,4 @@ namespace dd {
         return r;
     }
 
-
 }

src/math/dd/dd_bdd.h

@@ -259,6 +259,7 @@ namespace dd {
         bool is_constv(bddv const& a);
         rational to_val(bddv const& a);
 
+
         std::ostream& display(std::ostream& out);
         std::ostream& display(std::ostream& out, bdd const& b);
 

src/math/dd/dd_fdd.cpp

@@ -120,4 +120,36 @@ namespace dd {
         return out;
     }
 
+  
+    bool fdd::sup(bdd const& x, bool_vector& lo) {
+        SASSERT(lo.size() == num_bits());
+	//
+	// Assumption: common case is that high-order bits are before lower-order bits also
+	// after re-ordering.
+	//
+
+	// this checks that lo is included in x
+	bdd b = x;
+	while (!b.is_true()) {
+	    unsigned const pos = var2pos(b.var());
+            SASSERT(pos != UINT_MAX && "Unexpected BDD variable");
+	    if (lo[pos] && b.hi().is_false())
+		return false;
+	    if (!lo[pos] && b.lo().is_false())
+		return false;
+	    
+	    if (lo[pos])
+		b = b.hi();
+	    else
+		b = b.lo();
+	}
+	return false;
+    }
+    
+    bool fdd::inf(bdd const& b, bool_vector& hi) {
+        SASSERT(hi.size() == num_bits());
+	return false;
+    }
+
+
 }

src/math/dd/dd_fdd.h

@@ -11,8 +11,8 @@ Abstract:
 
 Author:
 
-    Nikolaj Bjorner (nbjorner) 2021-04-20
     Jakob Rath 2021-04-20
+    Nikolaj Bjorner (nbjorner) 2021-04-20
 
 --*/
 #pragma once
@@ -71,6 +71,23 @@ namespace dd {
 
         /** Like find, but returns hint if it is contained in the BDD. */
         find_t find_hint(bdd b, rational const& hint, rational& out_val) const;
+
+
+
+	/*
+	 * find largest value at lo or above such that bdd b evaluates to true
+	 * at lo and all values between.
+	 * dually, find smallest value below hi that evaluates b to true
+	 * and all values between the value and hi also evaluate b to true.
+	 * \param b - a bdd using variables from this
+	 * \param lo/hi - bound to be traversed.
+	 * \return false if b is false at lo/hi
+	 * \pre variables in b are a subset of variables from the fdd
+	 */ 
+	bool sup(bdd const& b, bool_vector& lo);
+      
+	bool inf(bdd const& b, bool_vector& hi); 
+      
     };
 
 }

src/math/polysat/viable.cpp

@@ -10,6 +10,13 @@ Author:
     Nikolaj Bjorner (nbjorner) 2021-03-19
     Jakob Rath 2021-04-6
 
+Notes:
+
+Use cheap heuristics to narrow viable sets whenever possible.
+If the cheap heuristics fail, compute a BDD representing the viable sets
+and narrow the range using the BDDs that are cached.
+
+
 --*/
 
 #include "math/polysat/viable.h"
@@ -24,13 +31,10 @@ namespace polysat {
     dd::find_t viable_set::find_hint(rational const& d, rational& val) const {
         if (is_empty())
             return dd::find_t::empty;         
-        // 
+        if (contains(d))
-        // ignore d since with a single interval, 
+            val = d;
-        // backtracking does not ensure that non-boundary bounds are removed.
+        else 
-        // viable_set could have multiple intervals to support arbitrary partitions
+            val = lo;
-        // this is similar to interval_set or might even be an instance of it.
-        // 
-        val = lo;
         if (lo + 1 == hi || hi == 0 && is_max(lo))
             return dd::find_t::singleton;            
         return dd::find_t::multiple;
@@ -40,7 +44,7 @@ namespace polysat {
         return a + 1 == rational::power_of_two(m_num_bits);
     }
     
-    bool viable_set::is_singleton(rational& val) const {
+    bool viable_set::is_singleton() const {
         return !is_empty() && (lo + 1 == hi || (hi == 0 && is_max(lo)));
     }
 
@@ -60,9 +64,9 @@ namespace polysat {
                 return;
             if (a == lo && a + 1 == hi)
                 set_empty();
-            if (a == lo && hi == 0 && is_max(a))
+            else if (a == lo && hi == 0 && is_max(a))
                 set_empty();
-            else if (a == lo)
+            else if (a == lo && !is_max(a))
                 lo = a + 1;
             else if (a + 1 == hi)
                 hi = a;
@@ -80,14 +84,16 @@ namespace polysat {
             else {
                 rational a_inv;
                 VERIFY(a.mult_inverse(m_num_bits, a_inv));
-                intersect_eq(mod(a_inv * -b, rational::power_of_two(m_num_bits)), is_positive);
+                intersect_eq(mod(a_inv * -b, p2()), is_positive);
             }
         }
+	    else
+	        std::cout << "intersect " << a << "*x " << " == " << b << " " << is_positive << "\n";
     }
 
-    void viable_set::intersect_ule(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
+    bool viable_set::intersect_ule(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
         // x <= 0
-        if (a == 1 && b == 0 && c == 0 && d == 0) 
+        if (a.is_odd() && b == 0 && c == 0 && d == 0)
             intersect_eq(b, is_positive);
         else if (a == 1 && b == 0 && c == 0) {
             // x <= d
@@ -108,6 +114,31 @@ namespace polysat {
             else
                 set_hi(b - 1);
         }
+        else if (c == 0 && d == 0) {
+            // ax + b <= 0
+            // or ax + b > 0
+            intersect_eq(a, b, is_positive);
+        }
+        else
+            return false;
+
+        return true;
+    }
+
+    void viable_set::intersect_ule(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive, unsigned& budget) {
+        auto eval = [&](rational const& x) {
+            return is_positive == mod(a * x + b, p2()) <= mod(c * x + d, p2());
+        };        
+        while (budget > 0 && !eval(lo) && !is_max(lo) && !is_empty()) {
+            --budget;
+            lo += 1;
+            set_lo(lo);
+        }
+        while (budget > 0 && hi > 0 && !eval(hi - 1) && !is_empty()) {
+            --budget;
+            hi = hi - 1;
+            set_hi(hi);
+        }
     }
 
     void viable_set::set_hi(rational const& d) {
@@ -136,11 +167,8 @@ namespace polysat {
 #endif
 
     viable::viable(solver& s):
-        s(s)
+        s(s),
-#if !NEW_VIABLE
-        ,
         m_bdd(1000)
-#endif
     {}
 
     void viable::push_viable(pvar v) {
@@ -191,7 +219,35 @@ namespace polysat {
     void viable::intersect_ule(pvar v, rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive) {
 #if NEW_VIABLE
         push_viable(v);
-        m_viable[v].intersect_ule(a, b, c, d, is_positive);
+        if (!m_viable[v].intersect_ule(a, b, c, d, is_positive)) {
+            unsigned budget = 10;
+            m_viable[v].intersect_ule(a, b, c, d, is_positive, budget);
+            if (budget == 0) {
+                std::cout << "miss: " << a << " " << b << " " << c << " " << d << " " << is_positive << "\n";
+                unsigned sz = var2bits(v).num_bits();
+                bdd le = m_bdd.mk_true();
+                ineq_entry entry0(sz, a, b, c, d, le);
+                ineq_entry* other = nullptr;
+                if (!m_ineq_cache.find(&entry0, other)) {
+                    std::cout << "ADD-to-cache\n";
+                    bddv const& x = var2bits(v).var();
+                    le = ((a * x) + b) <= ((c * x) + d);
+                    other = alloc(ineq_entry, sz, a, b, c, d, le);
+                    m_ineq_cache.insert(other);
+                }
+                le = is_positive ? other->repr : !other->repr;
+                other->m_activity++;
+                // le(lo) is false: find min x >= lo, such that le(x) is false, le(x+1) is true
+                // le(hi) is false: find max x =< hi, such that le(x) is false, le(x-1) is true
+                
+                // bdd x_is_lo = m.mk_num(sz, m_viable[v].lo) == var2bits(v).var();
+                // bdd lo_is_sat = x_is_lo && lo;
+                // if (lo_is_sat.is_false())
+                //     find_min_above(lo, le);
+
+            }
+
+        }
         if (m_viable[v].is_empty())
             s.set_conflict(v);
 #else
@@ -230,6 +286,7 @@ namespace polysat {
 #if NEW_VIABLE
         push_viable(v);
         m_viable[v].set_ne(val);
+        std::cout << " v" << v << " != " << val << "\n";
         if (m_viable[v].is_empty())
             s.set_conflict(v);
 #else
@@ -257,7 +314,6 @@ namespace polysat {
 #endif
     }
 
-#if !NEW_VIABLE
     dd::fdd const& viable::sz2bits(unsigned sz) {
         m_bits.reserve(sz + 1);
         auto* bits = m_bits[sz];
@@ -267,7 +323,6 @@ namespace polysat {
         }
         return *bits;
     }
-#endif
 
 #if POLYSAT_LOGGING_ENABLED
     void viable::log() {
@@ -291,9 +346,7 @@ namespace polysat {
     }
 #endif
 
-#if !NEW_VIABLE
     dd::fdd const& viable::var2bits(pvar v) { return sz2bits(s.size(v)); }
-#endif
 
 
 }

src/math/polysat/viable.h

@@ -41,29 +41,54 @@ namespace polysat {
     //
     class viable_set : public mod_interval<rational> {
         unsigned m_num_bits;
+        rational p2() const { return rational::power_of_two(m_num_bits); }
         bool is_max(rational const& a) const;
         void set_lo(rational const& lo);
         void set_hi(rational const& hi);
         void intersect_eq(rational const& a, bool is_positive);
     public:
         viable_set(unsigned num_bits): m_num_bits(num_bits) {}
-        bool is_singleton(rational& val) const; 
+        bool is_singleton() const; 
         dd::find_t find_hint(rational const& c, rational& val) const;
         void set_ne(rational const& a) { intersect_eq(a, false); }
         void intersect_eq(rational const& a, rational const& b, bool is_positive);
-        void intersect_ule(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive);
+        bool intersect_ule(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive);
+        void intersect_ule(rational const& a, rational const& b, rational const& c, rational const& d, bool is_positive, unsigned& budget);
     };
 #endif
 
     class viable {
         solver& s;
-#if NEW_VIABLE
-        vector<viable_set>           m_viable; // future representation of viable values
-        vector<std::pair<pvar, viable_set>> m_viable_trail;
-#else
         typedef dd::bdd bdd;
+        typedef dd::fdd fdd;
         dd::bdd_manager              m_bdd;
         scoped_ptr_vector<dd::fdd>   m_bits;
+#if NEW_VIABLE
+        struct ineq_entry {
+            unsigned m_num_bits;
+            rational a, b, c, d;
+            bdd      repr;
+            unsigned m_activity = 0;
+            ineq_entry(unsigned n, rational const& a, rational const& b, rational const& c, rational const& d, bdd& f) :
+                m_num_bits(n), a(a), b(b), c(c), d(d), repr(f) {}
+
+            struct hash {
+                unsigned operator()(ineq_entry const* e) const {
+                    return mk_mix(e->a.hash(), e->b.hash(), mk_mix(e->c.hash(), e->d.hash(), e->m_num_bits));
+                }
+            };
+            struct eq {
+                bool operator()(ineq_entry const* x, ineq_entry const* y) const {
+                    return x->a == y->a && x->b == y->b && x->c == y->c && x->d == y->d && x->m_num_bits == y->m_num_bits;
+                }
+            };
+        };
+        vector<viable_set>                                       m_viable; 
+        vector<std::pair<pvar, viable_set>>                      m_viable_trail;
+        hashtable<ineq_entry*, ineq_entry::hash, ineq_entry::eq> m_ineq_cache;
+#else
+
+        
         vector<bdd>                  m_viable;   // set of viable values.
         vector<std::pair<pvar, bdd>> m_viable_trail;
 
@@ -72,11 +97,11 @@ namespace polysat {
          * Register all values that are not contained in vals as non-viable.
          */
         void intersect_viable(pvar v, bdd vals);
-
+#endif
         dd::bdd_manager& get_bdd() { return m_bdd; }
         dd::fdd const& sz2bits(unsigned sz);
         dd::fdd const& var2bits(pvar v);
-#endif
+
 
     public:
         viable(solver& s);