add init_model, global m_delta, get_value, get_ivalue to push model maintainance into lar_solver #4740

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

src/math/lp/lar_solver.cpp

@@ -14,7 +14,8 @@ lp_settings & lar_solver::settings() { return m_settings;}
 
     lp_settings const& lar_solver::settings() const { return m_settings; }
 
-void clear() {lp_assert(false); // not implemented
+    void clear() {
+        lp_assert(false); // not implemented
     }
 
 
@@ -110,7 +111,8 @@ bool lar_solver::implied_bound_is_correctly_explained(implied_bound const & be,
                 kind = static_cast<lconstraint_kind>(-kind);
             }
             rs_of_evidence /= ratio;
-    } else {
+        }
+        else {
             lar_term const& t = get_term(be.m_j);
             auto first_coeff = t.begin();
             unsigned j = (*first_coeff).column();
@@ -407,7 +409,8 @@ bool lar_solver::move_non_basic_column_to_bounds(unsigned j, bool force_change)
                 else
                     set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
                 return true;
-        } else if (force_change && m_settings.random_next() % 3 == 0) {
+            }
+            else if (force_change && m_settings.random_next() % 3 == 0) {
                 set_value_for_nbasic_column(j,
                     at_l ? lcs.m_r_upper_bounds()[j] : lcs.m_r_lower_bounds()[j]);
                 return true;
@@ -510,7 +513,8 @@ lp_status lar_solver::maximize_term(unsigned j_or_term,
         auto backup = m_mpq_lar_core_solver.m_r_x;
         if (m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis()) {
             prev_value = term.apply(m_mpq_lar_core_solver.m_r_x);
-    } else {
+        }
+        else {
             m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = false;
             if (solve() != lp_status::OPTIMAL)
                 return lp_status::UNBOUNDED;
@@ -593,7 +597,8 @@ void lar_solver::register_monoid_in_map(std::unordered_map<var_index, mpq> & coe
         auto it = coeffs.find(j);
         if (it == coeffs.end()) {
             coeffs[j] = a;
-    } else {
+        }
+        else {
             it->second += a;
         }
     }
@@ -606,7 +611,8 @@ void lar_solver::substitute_terms_in_linear_expression(const vector<std::pair<mp
             unsigned j = t.second;
             if (!tv::is_term(j)) {
                 register_monoid_in_map(coeffs, t.first, j);
-        } else {
+            }
+            else {
                 const lar_term& term = *m_terms[tv::unmask_term(t.second)];
 
                 for (auto p : term) {
@@ -696,7 +702,8 @@ void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j,
                     (m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj) ? "feas" : "inf") << std::endl;);
 
             }
-    } else {
+        }
+        else {
             m_column_buffer.clear();
             m_column_buffer.resize(A_r().row_count());
             m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer);
@@ -714,10 +721,12 @@ void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned
                 m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
                 if (was_infeas != m_mpq_lar_core_solver.m_r_solver.inf_set_contains(j))
                     m_basic_columns_with_changed_cost.insert(j);
-        } else {
+            }
+            else {
                 m_mpq_lar_core_solver.m_r_solver.track_column_feasibility(j);
             }
-    } else {
+        }
+        else {
             numeric_pair<mpq> delta;
             if (m_mpq_lar_core_solver.m_r_solver.make_column_feasible(j, delta))
                 change_basic_columns_dependend_on_a_given_nb_column(j, delta);
@@ -803,7 +812,8 @@ void lar_solver::add_last_rows_to_lu(lp_primal_core_solver<K,L> & s) {
             if (f->m_refactor_counter + columns_to_replace.size() >= 200 || f->has_dense_submatrix()) {
                 delete f;
                 f = nullptr;
-        } else {
+            }
+            else {
                 f->add_last_rows_to_B(s.m_basis_heading, columns_to_replace);
             }
         }
@@ -1093,8 +1103,6 @@ void lar_solver::get_infeasibility_explanation(explanation& exp) const {
         lp_assert(explanation_is_correct(exp));
     }
 
-
-
     void lar_solver::get_infeasibility_explanation_for_inf_sign(
         explanation& exp,
         const vector<std::pair<mpq, unsigned>>& inf_row,
@@ -1115,41 +1123,47 @@ void lar_solver::get_infeasibility_explanation_for_inf_sign(
 
     // (x, y) != (x', y') => (x + delta*y) != (x' + delta*y')
     void lar_solver::get_model(std::unordered_map<var_index, mpq>& variable_values) const {
-    if (!(get_status() == lp_status::OPTIMAL || get_status() == lp_status::FEASIBLE) ||
-        !m_columns_with_changed_bounds.empty()) {
         variable_values.clear();
+        if (!init_model())
             return;
+
+        unsigned n = m_mpq_lar_core_solver.m_r_x.size();
+
+        for (unsigned j = 0; j < n; j++) 
+            variable_values[j] = get_value(column_index(j));
+
+        TRACE("lar_solver_model", tout << "delta = " << m_delta << "\nmodel:\n";
+               for (auto p : variable_values) tout << this->get_variable_name(p.first) << " = " << p.second << "\n";);
     }
 
+    bool lar_solver::init_model() const {
+        if (get_status() != lp_status::OPTIMAL && get_status() != lp_status::FEASIBLE)
+            return false;
+        if (!m_columns_with_changed_bounds.empty())
+            return false;
+
         lp_assert(m_mpq_lar_core_solver.m_r_solver.calc_current_x_is_feasible_include_non_basis());
-    variable_values.clear();
+        m_delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1));
-    mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1)); 
         unsigned j;
         unsigned n = m_mpq_lar_core_solver.m_r_x.size();
-    std::unordered_set<impq> set_of_different_pairs; 
-    std::unordered_set<mpq> set_of_different_singles;
         do {
-        set_of_different_pairs.clear();
+            m_set_of_different_pairs.clear();
-        set_of_different_singles.clear();
+            m_set_of_different_singles.clear();
             for (j = 0; j < n; j++) {
                 const numeric_pair<mpq>& rp = m_mpq_lar_core_solver.m_r_x[j];
-            set_of_different_pairs.insert(rp);
+                mpq x = rp.x + m_delta * rp.y;
-            mpq x = rp.x + delta * rp.y;
+                m_set_of_different_pairs.insert(rp);
-            set_of_different_singles.insert(x);
+                m_set_of_different_singles.insert(x);
-            if (set_of_different_pairs.size() != set_of_different_singles.size()) {
+                if (m_set_of_different_pairs.size() != m_set_of_different_singles.size()) {
-                delta /= mpq(2);
+                    m_delta /= mpq(2);
                     break;
                 }
-            variable_values[j] = x;
             }
         } while (j != n);
-    TRACE("lar_solver_model", tout << "delta = " << delta << "\nmodel:\n";
+        TRACE("lar_solver_model", tout << "delta = " << m_delta << "\nmodel:\n";);
-          for (auto p : variable_values ) {
+        return true;
-              tout << this->get_variable_name(p.first) << " = " << p.second << "\n";
-          });
     }
 
-
     void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq>& variable_values) const {
         mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1));
         for (unsigned i = 0; i < m_mpq_lar_core_solver.m_r_x.size(); i++) {
@@ -1158,6 +1172,30 @@ void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<var_in
         }
     }
 
+    mpq lar_solver::get_value(column_index const& j) const {
+        SASSERT(get_status() == lp_status::OPTIMAL || get_status() == lp_status::FEASIBLE);
+        numeric_pair<mpq> const& rp = get_column_value(j);
+        return rp.x + m_delta * rp.y;        
+    }
+
+    mpq lar_solver::get_value(tv const& t) const {
+        if (t.is_var())
+            get_value(t.column());
+        mpq r(0);
+        for (const auto& p : get_term(t)) 
+            r += p.coeff() * get_value(p.column());
+        return r;
+    }
+
+    impq lar_solver::get_ivalue(tv const& t) const {
+        if (t.is_var())
+            get_ivalue(t.column());
+        impq r;
+        for (const auto& p : get_term(t)) 
+            r += p.coeff() * get_ivalue(p.column());
+        return r;
+    }
+
     void lar_solver::get_rid_of_inf_eps() {
         bool y_is_zero = true;
         for (unsigned j = 0; j < number_of_vars(); j++) {
@@ -1218,9 +1256,11 @@ std::ostream& lar_solver::print_term(lar_term const& term, std::ostream & out) c
             mpq val = p.coeff();
             if (first) {
                 first = false;
-        } else if (is_pos(val)) {
+            }
+            else if (is_pos(val)) {
                 out << " + ";
-        } else {
+            }
+            else {
                 out << " - ";
                 val = -val;
             }
@@ -1258,7 +1298,8 @@ void lar_solver::fill_var_set_for_random_update(unsigned sz, var_index const * v
                 if (term_is_used_as_row(tv::unmask_term(var))) {
                     column_list.push_back(map_term_index_to_column_index(var));
                 }
-        } else {
+            }
+            else {
                 column_list.push_back(var);
             }
         }
@@ -1354,7 +1395,8 @@ void lar_solver::remove_last_column_from_basis_tableau(unsigned j) {
                 rslv.m_basis_heading[j_at_last_pos] = i;
             }
             rslv.m_basis.pop_back(); // remove j from the basis
-    } else {
+        }
+        else {
             int last_pos = static_cast<int>(rslv.m_nbasis.size()) - 1;
             lp_assert(last_pos >= 0);
             i = -1 - i;
@@ -1626,7 +1668,8 @@ bool lar_solver::term_coeffs_are_ok(const vector<std::pair<mpq, var_index>> & co
             if (!g_is_set) {
                 g_is_set = true;
                 g = p.first;
-        } else {
+            }
+            else {
                 g = gcd(g, p.first);
             }
         }
@@ -1770,7 +1813,8 @@ void lar_solver::register_in_fixed_var_table(unsigned j, unsigned & equal_to_j)
                 m_fixed_var_table_int.insert(key, j);
                 return;
             }
-    } else { // j is not integral column        
+        }
+        else { // j is not integral column        
             if (!m_fixed_var_table_real.find(key, k)) {
                 m_fixed_var_table_real.insert(key, j);
                 return;
@@ -1988,7 +2032,8 @@ void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstrain
         SASSERT(column_has_upper_bound(j));
         if (column_has_lower_bound(j)) {
             update_bound_with_ub_lb(j, kind, right_side, constraint_index);
-    } else {
+        }
+        else {
             update_bound_with_ub_no_lb(j, kind, right_side, constraint_index);
         }
     }
@@ -1997,7 +2042,8 @@ void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstr
         SASSERT(!column_has_upper_bound(j));
         if (column_has_lower_bound(j)) {
             update_bound_with_no_ub_lb(j, kind, right_side, constraint_index);
-    } else {
+        }
+        else {
             update_bound_with_no_ub_no_lb(j, kind, right_side, constraint_index);
         }
     }
@@ -2257,7 +2303,8 @@ void lar_solver::round_to_integer_solution() {
             if (del < -impq(mpq(1, 2))) {
                 del = impq(one_of_type<mpq>()) + del;
                 v = impq(ceil(v));
-        } else {
+            }
+            else {
                 v = flv;
             }
             m_incorrect_columns.insert(j);
@@ -2342,10 +2389,12 @@ void lar_solver::set_cut_strategy(unsigned cut_frequency) {
         if (cut_frequency < 4) {
             settings().m_int_gomory_cut_period = 2; // do it often
             settings().set_hnf_cut_period(4); // also create hnf cuts
-    } else if (cut_frequency == 4) { // enable all cuts and cube equally
+        }
+        else if (cut_frequency == 4) { // enable all cuts and cube equally
             settings().m_int_gomory_cut_period = 4;
             settings().set_hnf_cut_period(4);
-    } else {
+        }
+        else {
             // disable all heuristics except cube
             settings().m_int_gomory_cut_period = 10000000;
             settings().set_hnf_cut_period(100000000);
@@ -2360,7 +2409,8 @@ void lar_solver::register_normalized_term(const lar_term& t, lpvar j) {
         tout << ", normalized_t="; print_term_as_indices(normalized_t, tout) << "\n";);
         if (m_normalized_terms_to_columns.find(normalized_t) == m_normalized_terms_to_columns.end()) {
             m_normalized_terms_to_columns[normalized_t] = std::make_pair(a, j);
-    } else {
+        }
+        else {
             TRACE("lar_solver_terms", tout << "the term has been seen already\n";);
         }
     }

src/math/lp/lar_solver.h

@@ -137,10 +137,7 @@ class lar_solver : public column_namer {
 
     inline void clear_columns_with_changed_bounds() { m_columns_with_changed_bounds.clear(); }
     inline void increase_by_one_columns_with_changed_bounds() { m_columns_with_changed_bounds.increase_size_by_one(); }
-    inline void insert_to_columns_with_changed_bounds(unsigned j) {
+    inline void insert_to_columns_with_changed_bounds(unsigned j) { m_columns_with_changed_bounds.insert(j); }
-        m_columns_with_changed_bounds.insert(j);
-    }
-    
     
     void update_column_type_and_bound_check_on_equal(unsigned j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index, unsigned&);
     void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index);
@@ -289,6 +286,11 @@ class lar_solver : public column_namer {
     void collect_rounded_rows_to_fix();
     void register_normalized_term(const lar_term&, lpvar);
     void deregister_normalized_term(const lar_term&);
+
+    mutable std::unordered_set<impq> m_set_of_different_pairs;
+    mutable std::unordered_set<mpq>  m_set_of_different_singles;
+    mutable mpq m_delta;
+
 public:
     const map<mpq, unsigned, obj_hash<mpq>, default_eq<mpq>>& fixed_var_table_int() const {
         return m_fixed_var_table_int;
@@ -520,6 +522,10 @@ public:
     std::ostream& print_constraint_indices_only(const lar_base_constraint * c, std::ostream & out) const;
     std::ostream& print_implied_bound(const implied_bound& be, std::ostream & out) const;
     std::ostream& print_values(std::ostream& out) const;
+    bool init_model() const;
+    mpq get_value(column_index const& j) const;
+    mpq get_value(tv const& t) const;
+    impq get_ivalue(tv const& t) const;
     void get_model(std::unordered_map<var_index, mpq> & variable_values) const;
     void get_rid_of_inf_eps();
     void get_model_do_not_care_about_diff_vars(std::unordered_map<var_index, mpq> & variable_values) const;
@@ -615,6 +621,7 @@ public:
     inline const vector<unsigned> & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); }
     inline bool column_is_real(unsigned j) const { return !column_is_int(j); }	
     lp_status get_status() const;
+    bool has_changed_columns() const { return !m_columns_with_changed_bounds.empty();  }
     void set_status(lp_status s);
     lp_status solve();
     void fill_explanation_from_crossed_bounds_column(explanation & evidence) const;
@@ -631,8 +638,8 @@ public:
     inline const static_matrix<mpq, impq> & A_r() const { return m_mpq_lar_core_solver.m_r_A; }
     // columns
     bool column_is_int(column_index const& j) const { return column_is_int((unsigned)j); }
-    const impq& get_value(column_index const& j) const { return get_column_value(j); }
+    const impq& get_ivalue(column_index const& j) const { return get_column_value(j); }
-    const impq& get_column_value(unsigned j) const { return m_mpq_lar_core_solver.m_r_x[j]; }
+    const impq& get_column_value(column_index const& j) const { return m_mpq_lar_core_solver.m_r_x[j]; }
     inline
     var_index external_to_local(unsigned j) const {
         var_index local_j;

src/math/lp/lp_types.h

@@ -36,6 +36,20 @@ typedef unsigned lpvar;
 const lpvar null_lpvar = UINT_MAX;
 const constraint_index null_ci = UINT_MAX;
 
+class column_index {
+    unsigned m_index;
+    friend class lar_solver;
+    friend class lar_term;
+    friend nla::core;
+
+    operator unsigned() const { return m_index; }
+    
+public:
+    column_index(unsigned j): m_index(j) {}    
+    unsigned index() const { return m_index; }
+    bool is_null() const { return m_index == null_lpvar; }
+};
+
 
 // index that comes from term or variable.
 class tv {
@@ -49,6 +63,7 @@ public:
 
     // retrieve the identifier associated with tv
     unsigned id() const { return unmask_term(m_index); }
+    column_index column() const { SASSERT(is_var()); return column_index(id()); }
 
     // retrieve the raw index.
     unsigned index() const { return m_index; }
@@ -74,20 +89,6 @@ public:
 
 };
 
-class column_index {
-    unsigned m_index;
-    friend class lar_solver;
-    friend class lar_term;
-    friend nla::core;
-
-    operator unsigned() const { return m_index; }
-    
-public:
-    column_index(unsigned j): m_index(j) {}    
-    unsigned index() const { return m_index; }
-    bool is_null() const { return m_index == null_lpvar; }
-};
-
 }
 
 inline std::ostream& operator<<(std::ostream& out, lp::tv const& t) {

src/sat/smt/euf_solver.cpp

@@ -231,7 +231,6 @@ namespace euf {
         return true;
     }
 
-
     void solver::get_antecedents(literal l, constraint& j, literal_vector& r, bool probing) {
         expr* e = nullptr;
         euf::enode* n = nullptr;

src/smt/theory_lra.cpp

@@ -1474,6 +1474,7 @@ public:
             }
             else {
                 auto I = m_variable_values.find(t2.index());
+                std::cout << (I == E) << "\n";
                 if (I != E)
                     result += I->second * coeff;
             }
@@ -1619,7 +1620,7 @@ public:
         IF_VERBOSE(12, verbose_stream() << "final-check " << lp().get_status() << "\n");
         lbool is_sat = l_true;
         SASSERT(lp().ax_is_correct());
-        if (lp().get_status() != lp::lp_status::OPTIMAL) {
+        if (lp().get_status() != lp::lp_status::OPTIMAL || lp().has_changed_columns()) {
             is_sat = make_feasible();
         }
         final_check_status st = FC_DONE;