working on hitting sets

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

src/opt/hitting_sets.cpp

@@ -0,0 +1,370 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    hitting_sets.h
+
+Abstract:
+   
+    Hitting set approximations.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-06-06
+
+Notes:
+
+--*/
+#include "vector.h"
+#include "util.h"
+#include "hitting_sets.h"
+#include "simplex.h"
+#include "sparse_matrix_def.h"
+#include "simplex_def.h"
+
+typedef simplex::simplex<simplex::mpz_ext> Simplex;
+typedef simplex::sparse_matrix<simplex::mpz_ext> sparse_matrix;
+
+
+namespace opt {
+
+    struct hitting_sets::imp {
+        typedef unsigned_vector set;
+        volatile bool           m_cancel;
+        rational                m_lower;
+        rational                m_upper;
+        vector<rational>        m_weights;
+        rational                m_max_weight;
+        rational                m_denominator;
+        vector<set>             m_S;
+        svector<lbool>          m_value;
+        unsigned_vector         m_value_trail;
+        unsigned_vector         m_value_lim;
+        vector<unsigned_vector> m_use_list;        
+        unsynch_mpz_manager     m;
+        Simplex                 m_simplex;
+        unsigned                m_weights_var;
+
+        imp():m_cancel(false) {}
+        ~imp() {}
+
+        void add_weight(rational const& w) {
+            SASSERT(w.is_pos());
+            unsigned var = m_weights.size();
+            m_simplex.ensure_var(var);
+            m_simplex.set_lower(var, mpq_inf(mpq(0),mpq(0)));
+            m_simplex.set_upper(var, mpq_inf(mpq(1),mpq(0)));
+            m_weights.push_back(w);
+            m_value.push_back(l_undef);
+            m_use_list.push_back(unsigned_vector());
+            m_max_weight += w;
+        }
+
+        void add_set(unsigned sz, unsigned const* S) {
+            if (sz == 0) {
+                return;
+            }
+            for (unsigned i = 0; i < sz; ++i) {
+                m_use_list[S[i]].push_back(m_S.size());
+            }
+            init_weights();
+            m_S.push_back(unsigned_vector(sz, S));        
+            add_simplex_row(sz, S);            
+        }
+
+        bool compute_lower() {
+            m_lower.reset();
+            return L1() && L2() && L3();            
+        }
+
+        bool compute_upper() {            
+            m_upper = m_max_weight;
+            return U1();
+        }
+
+        rational get_lower() {
+            return m_lower/m_denominator;
+        }
+
+        rational get_upper() {
+            return m_upper/m_denominator;
+        }
+
+        void set_cancel(bool f) {
+            m_cancel = f;
+            m_simplex.set_cancel(f);
+        }
+
+        void collect_statistics(::statistics& st) const {
+            m_simplex.collect_statistics(st);
+        }
+
+        void reset() {
+            m_lower.reset();
+            m_upper = m_max_weight;
+        }
+
+        void init_weights() {
+            if (m_weights_var != 0) {
+                return;
+            }
+            m_weights_var = m_weights.size();
+            unsigned_vector vars;
+            scoped_mpz_vector coeffs(m);
+
+            // normalize weights to integral.
+            rational d(1);
+            for (unsigned i = 0; i < m_weights.size(); ++i) {
+                d = lcm(d, denominator(m_weights[i]));
+            }
+            m_denominator = d;
+            if (!d.is_one()) {
+                for (unsigned i = 0; i < m_weights.size(); ++i) {
+                    m_weights[i] *= d;
+                }
+            }
+            // set up Simplex objective function.
+            for (unsigned i = 0; i < m_weights.size(); ++i) {
+                vars.push_back(i);
+                coeffs.push_back(m_weights[i].to_mpq().numerator());
+            }
+            m_simplex.ensure_var(m_weights_var);
+            vars.push_back(m_weights_var);
+            coeffs.push_back(mpz(-1));
+            m_simplex.add_row(m_weights_var, coeffs.size(), vars.c_ptr(), coeffs.c_ptr());
+        }
+
+        struct scoped_select {
+            imp& s;
+            unsigned sz;
+            scoped_select(imp& s):s(s), sz(s.m_value_trail.size()) {
+            }
+            ~scoped_select() {
+                s.undo_select(sz);
+            }
+        };
+
+        struct value_lt {
+            vector<rational> const& weights;
+            unsigned_vector const& scores;
+            value_lt(vector<rational> const& weights, unsigned_vector const& scores):
+                weights(weights), scores(scores) {}
+            bool operator()(int v1, int v2) const {
+                //     - score1 / w1 < - score2 / w2
+                // <=> 
+                //     score1 / w1 > score2 / w2
+                // <=>
+                //     score1*w2 > score2*w1
+                unsigned score1 = scores[v1];
+                unsigned score2 = scores[v2];
+                rational w1 = weights[v1];
+                rational w2 = weights[v2];
+                return rational(score1)*w2 > rational(score2)*w1;
+            }
+        };
+
+
+        // compute upper bound for hitting set.
+        bool U1() {        
+            rational w(0);
+            scoped_select _sc(*this);
+            
+            // score each variable by the number of
+            // unassigned sets they occur in.
+            unsigned_vector scores;
+            init_scores(scores);
+            
+            //
+            // Sort indices.
+            // The least literals are those where -score/w is minimized.
+            //
+            unsigned_vector indices;
+            for (unsigned i = 0; i < m_value.size(); ++i) {
+                indices.push_back(i);
+            }        
+            value_lt lt(m_weights, scores);
+            while (!m_cancel) {
+                std::sort(indices.begin(), indices.end(), lt);
+                unsigned idx = indices[0];
+                if (scores[idx] == 0) {
+                    break;
+                }
+                update_scores(scores, idx);
+                select(idx);
+                w += m_weights[idx];
+            }
+            if (w < m_upper) {
+                m_upper = w;
+            }
+            return !m_cancel;
+        }
+
+        void init_scores(unsigned_vector & scores) {
+            scores.reset();
+            for (unsigned i = 0; i < m_value.size(); ++i) {
+                scores.push_back(0);
+            }
+            for (unsigned i = 0; i < m_S.size(); ++i) {
+                set const& S = m_S[i];
+                if (!has_selected(S)) {
+                    for (unsigned j = 0; j < S.size(); ++j) {
+                        scores[S[j]]++;
+                    }
+                }
+            }
+        }
+
+        void update_scores(unsigned_vector& scores, unsigned v) {
+            unsigned_vector const& v_uses = m_use_list[v];
+            for (unsigned i = 0; i < v_uses.size(); ++i) {
+                set const& S = m_S[v_uses[i]];
+                if (!has_selected(S)) {
+                    for (unsigned j = 0; j < S.size(); ++j) {
+                        --scores[S[j]];
+                    }
+                }
+            }
+        }
+
+        bool L1() {
+            rational w(0);
+            scoped_select _sc(*this);
+            for (unsigned i = 0; !m_cancel && i < m_S.size(); ++i) {
+                set const& S = m_S[i];
+                SASSERT(!S.empty());
+                if (!has_selected(S)) {
+                    w += m_weights[select_min(S)];                
+                    for (unsigned j = 0; j < S.size(); ++j) {
+                        select(S[j]);
+                    }
+                }
+            }
+            if (m_lower < w) {
+                m_lower = w;
+            }
+            return !m_cancel;
+        }
+
+        bool L2() {
+            rational w(0);
+            scoped_select _sc(*this);
+            int n = 0;
+            for (unsigned i = 0; i < m_S.size(); ++i) {
+                if (!has_selected(m_S[i])) ++n;
+            }
+            unsigned_vector scores;
+            init_scores(scores);
+            unsigned_vector indices;
+            for (unsigned i = 0; i < m_value.size(); ++i) {
+                indices.push_back(i);
+            }        
+            value_lt lt(m_weights, scores);
+
+            std::sort(indices.begin(), indices.end(), lt);
+            for(unsigned i = 0; i < indices.size() && n > 0; ++i) {
+                // deg(c) = score(c)
+                // wt(c) = m_weights[c]
+                unsigned idx = indices[i];
+                if (scores[idx] == 0) {
+                    break;
+                }
+                if (scores[idx] < static_cast<unsigned>(n) || m_weights[idx].is_one()) {
+                    w += m_weights[idx];
+                }
+                else {
+                    w += div((rational(n)*m_weights[idx]), rational(scores[idx]));
+                }
+                n -= scores[idx];
+            }
+            if (m_lower < w) {
+                m_lower = w;
+            }
+            return !m_cancel;
+        }
+
+        bool L3() {        
+            TRACE("simplex", m_simplex.display(tout););
+            VERIFY(l_true == m_simplex.make_feasible());
+            TRACE("simplex", m_simplex.display(tout););
+            VERIFY(l_true == m_simplex.minimize(m_weights_var));
+            mpq_inf const& val = m_simplex.get_value(m_weights_var);
+            unsynch_mpq_inf_manager mg;
+            unsynch_mpq_manager& mq = mg.mpq_manager();
+            scoped_mpq c(mq);
+            mg.ceil(val, c);
+            rational w = rational(c);
+            if (w > m_lower) {
+                m_lower = w;
+            }
+            return true;
+        }    
+
+        void add_simplex_row(unsigned sz, unsigned const* S) {
+            unsigned_vector vars;
+            scoped_mpz_vector coeffs(m);
+            for (unsigned i = 0; i < sz; ++i) {
+                vars.push_back(S[i]);
+                coeffs.push_back(mpz(1));
+            }
+            unsigned base_var = m_S.size() + m_weights.size();
+            m_simplex.ensure_var(base_var);
+            vars.push_back(base_var);
+            coeffs.push_back(mpz(-1));
+            // S - base_var = 0
+            // base_var >= 1
+            m_simplex.set_lower(base_var, mpq_inf(mpq(1),mpq(0)));
+            m_simplex.add_row(base_var, coeffs.size(), vars.c_ptr(), coeffs.c_ptr());
+        }
+
+        void undo_select(unsigned sz) {
+            for (unsigned j = sz; j < m_value_trail.size(); ++j) {
+                m_value[m_value_trail[j]] = l_undef;
+            }
+            m_value_trail.resize(sz);
+        }
+        
+        unsigned select_min(set const& S) {
+            unsigned result = S[0];
+            for (unsigned i = 1; i < S.size(); ++i) {
+                if (m_weights[result] > m_weights[S[i]]) {
+                    result = S[i];
+                }
+            }
+            return result;
+        }
+        
+        
+        lbool selected(unsigned j) const { 
+            return m_value[j];
+        }
+        
+        void select(unsigned j) {
+            m_value[j] = l_true;
+            m_value_trail.push_back(j);
+        }
+        
+        bool has_selected(set const& S) {
+            for (unsigned i = 0; i < S.size(); ++i) {
+                if (l_true == selected(S[i])) {
+                    return true;
+                }
+            }
+            return false;
+        }
+        
+    };
+    hitting_sets::hitting_sets() { m_imp = alloc(imp); }
+    hitting_sets::~hitting_sets() { dealloc(m_imp); }
+    void hitting_sets::add_weight(rational const& w) { m_imp->add_weight(w); }
+    void hitting_sets::add_set(unsigned sz, unsigned const* elems) { m_imp->add_set(sz, elems); }
+    bool hitting_sets::compute_lower() { return m_imp->compute_lower(); }
+    bool hitting_sets::compute_upper() { return m_imp->compute_upper(); }
+    rational hitting_sets::get_lower() { return m_imp->get_lower(); }
+    rational hitting_sets::get_upper() { return m_imp->get_upper(); }
+    void hitting_sets::set_cancel(bool f) { m_imp->set_cancel(f); }
+    void hitting_sets::collect_statistics(::statistics& st) const { m_imp->collect_statistics(st); }
+    void hitting_sets::reset() { m_imp->reset(); }
+
+
+};

src/opt/hitting_sets.h

@@ -0,0 +1,47 @@
+/*++
+Copyright (c) 2014 Microsoft Corporation
+
+Module Name:
+
+    hitting_sets.h
+
+Abstract:
+   
+    Hitting set approximations.
+
+Author:
+
+    Nikolaj Bjorner (nbjorner) 2014-06-06
+
+Notes:
+
+--*/
+#ifndef _HITTING_SETS_H_
+#define _HITTING_SETS_H_
+
+#include "rational.h"
+#include "statistics.h"
+
+namespace opt {
+
+    class hitting_sets {
+        struct imp;
+        imp* m_imp;
+    public:        
+        hitting_sets();
+        ~hitting_sets();
+        void add_weight(rational const& w);
+        void add_set(unsigned sz, unsigned const* elems);
+        bool compute_lower();
+        bool compute_upper();
+        rational get_lower();
+        rational get_upper();        
+        void set_cancel(bool f);
+        void collect_statistics(::statistics& st) const;
+        void reset();
+    };
+
+
+};
+
+#endif

src/opt/weighted_maxsat.cpp

@@ -38,8 +38,8 @@ Notes:
 #include "cancel_eh.h"
 #include "scoped_timer.h"
 #include "optsmt.h"
+#include "hitting_sets.h"
 
-#define USE_SIMPLEX 0
 
 namespace opt {
 
@@ -614,9 +614,7 @@ namespace opt {
         };
 
         scoped_ptr<maxsmt_solver_base> maxs;
-        optsmt                  m_optsmt;     // hitting set optimizer based on simplex.
-        opt_solver              m_solver;
-        unsigned                m_objective;  // index of objective
+        hitting_sets                   m_hs;
         expr_ref_vector         m_aux;        // auxiliary (indicator) variables.
         expr_ref_vector         m_iaux;       // auxiliary integer (indicator) variables.
         expr_ref_vector         m_naux;       // negation of auxiliary variables.
@@ -628,31 +626,29 @@ namespace opt {
         pb_util                 pb;
         arith_util              a;
         stats                   m_stats;
+        bool                    m_at_lower_bound;
 
 
     public:
         hsmax(solver* s, ast_manager& m, maxsmt_solver_base* maxs):
             maxsmt_solver_base(s, m), 
             maxs(maxs), 
-            m_optsmt(m),
-            m_solver(m, m_params, symbol()),
             m_aux(m), 
             m_iaux(m),
             m_naux(m),
             pb(m),
-            a(m) {
+            a(m),
+            m_at_lower_bound(false) {
         }
         virtual ~hsmax() {}
 
         virtual void set_cancel(bool f) { 
             maxsmt_solver_base::set_cancel(f); 
             maxs->set_cancel(f);
-            m_optsmt.set_cancel(f);
         }
 
         virtual void updt_params(params_ref& p) {
             maxsmt_solver_base::updt_params(p);
-            m_solver.updt_params(p);
         }
 
         virtual void collect_statistics(statistics& st) const {
@@ -685,8 +681,8 @@ namespace opt {
                 if (m_cancel) {
                     return l_undef;
                 }
-                lbool core_found = generate_cores(hs);
                 
+                lbool core_found = generate_cores(hs);
                 switch(core_found) {
                 case l_undef: 
                     return l_undef;
@@ -698,20 +694,20 @@ namespace opt {
                         break;
                     case l_false:
                         TRACE("opt", tout << "no more seeds\n";);
-                        m_lower = m_upper;
-                        return l_true;
+                         m_lower = m_upper;
+                         return l_true;
                     case l_undef:
                         return l_undef;
-                    }
-                    break;
+                     }
+                     break;
                 }
-                case l_false:
+                 case l_false:
                     TRACE("opt", tout << "no more cores\n";);
                     m_lower = m_upper;
                     return l_true;
-                }
-            }
-            return l_true;
+                 }
+             }
+             return l_true;
         }
 
     private:
@@ -740,26 +736,16 @@ namespace opt {
                 m_aux_active.push_back(false);
                 m_core_activity.push_back(0);
                 m_aux2index.insert(m_aux.back(), i);
-#if USE_SIMPLEX
-                m_aux2index.insert(m_iaux.back(), i);
-                fml = m.mk_and(a.mk_le(a.mk_numeral(rational::zero(), true), iaux), 
-                               a.mk_le(iaux, a.mk_numeral(rational::one(), true)));
-                rational const& w = m_weights[i];
-                sum.push_back(a.mk_mul(a.mk_numeral(w, w.is_int()), iaux));
-                m_solver.assert_expr(fml);
-#endif
                 if (tt) {
                     m_asms.push_back(m_aux.back());
                     ensure_active(i);
                 }
             }
-#if USE_SIMPLEX
-            obj = a.mk_add(sum.size(), sum.c_ptr());
-            m_objective = m_optsmt.add(obj);
-            m_optsmt.setup(m_solver);
-#else
             maxs->init_soft(m_weights, m_aux);
-#endif
+
+            for (unsigned i = 0; i < m_weights.size(); ++i) {
+                m_hs.add_weight(m_weights[i]);
+            }
             TRACE("opt", print_seed(tout););
         }
 
@@ -895,54 +881,51 @@ namespace opt {
             }
         }
 
+        lbool next_seed(ptr_vector<expr>& hs, lbool core_found) {
+            
+            if (core_found == l_false && m_at_lower_bound) {
+                return l_true;
+            }            
+            lbool is_sat = next_seed();
+            switch(is_sat) {
+            case l_true:
+                seed2hs(false, hs);
+                return m_at_lower_bound?l_true:l_false;
+            case l_false:
+                TRACE("opt", tout << "no more seeds\n";);
+                return l_true;
+            case l_undef:
+                return l_undef;
+            }
+            return l_undef;
+        }
+
         //
         // retrieve the next seed that satisfies state of maxs.
         // state of maxs must be satisfiable before optimization is called.
         // 
-        //
         // find a satisfying assignment to maxs state, that 
         // minimizes objective function.
         // 
         lbool next_seed() {
             scoped_stopwatch _sw(m_stats.m_aux_sat_time);
             TRACE("opt", tout << "\n";);
-#if USE_SIMPLEX
-            m_solver.display(std::cout);
-            lbool is_sat = m_optsmt.lex(m_objective);
-            if (is_sat == l_true) {
-                model_ref mdl;
-                m_optsmt.get_model(mdl);
-                
-                for (unsigned i = 0; i < num_soft(); ++i) {
-                    if (is_active(i)) {
-                        m_seed[i] = is_one(mdl, m_iaux[i].get());
-                    }
-                    else {
-                        m_seed[i] = false;
-                    }
-                }
-                print_seed(std::cout);
-                TRACE("opt", print_seed(tout););
-            }
-            
-#else
 
             // min c_i*(not x_i) for x_i are soft clauses.
             // max c_i*x_i for x_i are soft clauses
             
             lbool is_sat = l_true;
+            m_at_lower_bound = false;
+            expr_ref fml(m);
             if (m_lower.is_pos()) {
-                expr_ref fml(m);
                 solver::scoped_push _scope(maxs->s());
                 fml = pb.mk_le(num_soft(), m_weights.c_ptr(), m_naux.c_ptr(), m_lower);
                 maxs->add_hard(fml);
-                //fml = pb.mk_ge(num_soft(), m_weights.c_ptr(), m_naux.c_ptr(), m_lower);
-                //maxs->add_hard(fml);
-                std::cout << fml << "\n";
                 is_sat = maxs->s().check_sat(0,0);
                 if (is_sat == l_true) {
                     maxs->set_model();
                     extract_seed();
+                    m_at_lower_bound = true;
                     return l_true;
                 }
             }
@@ -951,17 +934,32 @@ namespace opt {
                 maxs->set_model();
             }
             else {
+                m_at_lower_bound = true;
                 return is_sat;
             }
-            is_sat = (*maxs)();            
-                                
+            is_sat = (*maxs)();
+
             if (is_sat == l_true) {
                 extract_seed();
             }
-#endif
             return is_sat;
         }
 
+#if 0
+            if (!m_hs.compute_upper()) {
+                return l_undef;
+            }
+            solver::scoped_push _scope(maxs->s());
+            fml = pb.mk_le(num_soft(), m_weights.c_ptr(), m_naux.c_ptr(), m_hs.get_upper());
+            IF_VERBOSE(0, verbose_stream() << "upper: " << m_hs.get_upper() << " " << m_upper << "\n";);
+            maxs->add_hard(fml);
+            TRACE("opt", tout << "checking with upper bound: " << m_hs.get_upper() << "\n";);
+            is_sat = maxs->s().check_sat(0,0);
+            std::cout << is_sat << "\n";
+            
+            // TBD: uper bound estimate does not include the negative constraints.
+#endif
+
         void extract_seed() {
             model_ref mdl;
             maxs->get_model(mdl);
@@ -1140,16 +1138,6 @@ namespace opt {
             }
             expr_ref_vector fmls(m);
             expr_ref fml(m);
-#if USE_SIMPLEX
-            for (unsigned i = 0; i < num_soft(); ++i) {
-                if (!indices.contains(i)) {
-                    fmls.push_back(m_iaux[i].get());
-                }
-            }
-            fml = a.mk_ge(a.mk_add(fmls.size(), fmls.c_ptr()), a.mk_numeral(rational::one(), true));
-            m_solver.assert_expr(fml);
-
-#else
             for (unsigned i = 0; i < num_soft(); ++i) {
                 if (!indices.contains(i)) {
                     fmls.push_back(m_aux[i].get());
@@ -1158,7 +1146,6 @@ namespace opt {
             fml = m.mk_or(fmls.size(), fmls.c_ptr());
             maxs->add_hard(fml);
             set_upper();
-#endif
             TRACE("opt", tout << fml << "\n";);
         }
 
@@ -1174,25 +1161,17 @@ namespace opt {
         void block_up() {
             expr_ref_vector fmls(m);
             expr_ref fml(m);
-#if USE_SIMPLEX
-            for (unsigned i = 0; i < m_asms.size(); ++i) {
-                unsigned index = m_aux2index.find(m_asms[i]);
-                m_core_activity[index]++;
-                fmls.push_back(m_iaux[index].get());
-            }
-            fml = a.mk_lt(a.mk_add(fmls.size(), fmls.c_ptr()), a.mk_numeral(rational(fmls.size()), true));
-            TRACE("opt", tout << fml << "\n";);
-            m_solver.assert_expr(fml);
-#else
+            unsigned_vector indices;
             for (unsigned i = 0; i < m_asms.size(); ++i) {
                 unsigned index = m_aux2index.find(m_asms[i]);
                 fmls.push_back(m.mk_not(m_asms[i]));
                 m_core_activity[index]++;
+                indices.push_back(index);
             }
             fml = m.mk_or(fmls.size(), fmls.c_ptr());
             TRACE("opt", tout << fml << "\n";);
+            m_hs.add_set(indices.size(), indices.c_ptr());
             maxs->add_hard(fml);            
-#endif
         }
 
 
@@ -1224,7 +1203,6 @@ namespace opt {
             rational r;
             expr_ref val(m);
             VERIFY(mdl->eval(e, val));
-            std::cout << mk_pp(e, m) << " |-> " << val << "\n";
             return a.is_numeral(val, r) && r.is_one();
         }
 
@@ -1684,10 +1662,12 @@ namespace opt {
                 m_maxsmt = alloc(bcd2, s.get(), m);
             }
             else if (m_engine == symbol("hsmax")) {
+                //m_params.set_bool("pb.enable_simplex", true);
                 ref<opt_solver> s0 = alloc(opt_solver, m, m_params, symbol());
                 s0->check_sat(0,0);                
                 maxsmt_solver_base* s2 = alloc(pbmax, s0.get(), m); // , s0->get_context());
                 s2->set_converter(s0->mc_ref().get());
+                
                 m_maxsmt = alloc(hsmax, s.get(), m, s2);
             }
             // NB: this is experimental one-round version of SLS