Merge pull request #2180 from levnach/Prover

snap variables to bounds when maximizing terms

src/util/lp/int_solver.cpp

@@ -263,7 +263,7 @@ lia_move int_solver::find_cube() {
     if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
         TRACE("cube", tout << "cannot find a feasiblie solution";);
         _sp.pop();
-        move_non_basic_columns_to_bounds();
+        m_lar_solver->move_non_basic_columns_to_bounds();
         find_feasible_solution();
         // it can happen that we found an integer solution here
         return !m_lar_solver->r_basis_has_inf_int()? lia_move::sat: lia_move::undef;
@@ -296,7 +296,7 @@ lia_move int_solver::gomory_cut() {
     if ((m_number_of_calls) % settings().m_int_gomory_cut_period != 0)
         return lia_move::undef;
 
-    if (move_non_basic_columns_to_bounds()) {
+    if (m_lar_solver->move_non_basic_columns_to_bounds()) {
 #if Z3DEBUG 
         lp_status st =
 #endif
@@ -441,53 +441,6 @@ int int_solver::find_any_inf_int_column_basis_first() {
    return find_inf_int_nbasis_column();
 }
 
-bool int_solver::move_non_basic_column_to_bounds(unsigned j) {
-    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
-    auto & val = lcs.m_r_x[j];
-    switch (lcs.m_column_types()[j]) {
-    case column_type::boxed:
-        if (val != lcs.m_r_lower_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) {
-            if (random() % 2 == 0)
-                set_value_for_nbasic_column(j, lcs.m_r_lower_bounds()[j]);
-            else
-                set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
-            return true;
-        }
-        break;
-    case column_type::lower_bound:
-        if (val != lcs.m_r_lower_bounds()[j]) {
-            set_value_for_nbasic_column(j, lcs.m_r_lower_bounds()[j]);
-            return true;
-        }
-        break;
-    case column_type::upper_bound:
-        if (val != lcs.m_r_upper_bounds()[j]) {
-            set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
-            return true;
-        }
-        break;
-    default:
-        if (is_int(j) && !val.is_int()) {
-            set_value_for_nbasic_column(j, impq(floor(val)));
-            return true;
-        }
-        break;
-    }
-    return false;
-}
-
-bool int_solver::move_non_basic_columns_to_bounds() {
-    auto & lcs = m_lar_solver->m_mpq_lar_core_solver;
-    bool change = false;
-    for (unsigned j : lcs.m_r_nbasis) {
-        if (move_non_basic_column_to_bounds(j))
-            change = true;
-    }
-
-    if (settings().simplex_strategy() == simplex_strategy_enum::tableau_costs)
-        m_lar_solver->update_x_and_inf_costs_for_columns_with_changed_bounds_tableau();
-    return change;
-}
 
 void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val) {
     lp_assert(!is_base(j));
@@ -498,13 +451,6 @@ void int_solver::set_value_for_nbasic_column_ignore_old_values(unsigned j, const
 }
 
 
-void int_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
-    lp_assert(!is_base(j));
-    auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x[j];
-    auto delta = new_val - x;
-    x = new_val;
-    m_lar_solver->change_basic_columns_dependend_on_a_given_nb_column(j, delta);
-}
 
 void int_solver::patch_nbasic_column(unsigned j, bool patch_only_int_vals) {
     auto & lcs = m_lar_solver->m_mpq_lar_core_solver; 
@@ -537,7 +483,7 @@ void int_solver::patch_nbasic_column(unsigned j, bool patch_only_int_vals) {
         if (inf_u || l <= u) {
             TRACE("patch_int",
                   tout << "patching with l: " << l << '\n';);
-            set_value_for_nbasic_column(j, l);
+            m_lar_solver->set_value_for_nbasic_column(j, l);
         }
         else {
             TRACE("patch_int",
@@ -546,12 +492,12 @@ void int_solver::patch_nbasic_column(unsigned j, bool patch_only_int_vals) {
     }
     else if (!inf_u) {
         u = m_is_one ? floor(u) : m * floor(u / m);
-        set_value_for_nbasic_column(j, u);
+        m_lar_solver->set_value_for_nbasic_column(j, u);
         TRACE("patch_int",
               tout << "patching with u: " << u << '\n';);
     }
     else {
-        set_value_for_nbasic_column(j, impq(0));
+        m_lar_solver->set_value_for_nbasic_column(j, impq(0));
         TRACE("patch_int",
               tout << "patching with 0\n";);
     }

src/util/lp/int_solver.h

@@ -55,7 +55,6 @@ public:
     explanation const& get_explanation() const { return m_ex; }
     bool is_upper() const { return m_upper; }
 
-    bool move_non_basic_column_to_bounds(unsigned j);
     bool is_base(unsigned j) const;
     bool is_real(unsigned j) const;
     const impq & lower_bound(unsigned j) const;
@@ -95,7 +94,6 @@ private:
     bool is_fixed(unsigned j) const;
     bool is_free(unsigned j) const;
     bool value_is_int(unsigned j) const;
-    void set_value_for_nbasic_column(unsigned j, const impq & new_val);
     void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val);
     bool non_basic_columns_are_at_bounds() const;
     bool is_feasible() const;
@@ -108,7 +106,6 @@ private:
     int get_kth_inf_int(unsigned) const;
     lp_settings& settings();
     const lp_settings& settings() const;
-    bool move_non_basic_columns_to_bounds();
     void branch_infeasible_int_var(unsigned);
     lia_move mk_gomory_cut(unsigned inf_col, const row_strip<mpq>& row);
     lia_move proceed_with_gomory_cut(unsigned j);

src/util/lp/lar_solver.cpp

@@ -270,9 +270,9 @@ void lar_solver::propagate_bounds_for_touched_rows(bound_propagator & bp) {
     }
 }
 
-lp_status lar_solver::get_status() const { return m_status;}
+lp_status lar_solver::get_status() const { return m_status; }
 
-void lar_solver::set_status(lp_status s) {m_status = s;}
+void lar_solver::set_status(lp_status s) { m_status = s; }
 
 lp_status lar_solver::find_feasible_solution() {
     m_settings.st().m_make_feasible++;
@@ -453,9 +453,10 @@ void lar_solver::set_costs_to_zero(const lar_term& term) {
     lp_assert(costs_are_zeros_for_r_solver());
 }
 
-void lar_solver::prepare_costs_for_r_solver(const lar_term & term) {
+void lar_solver::prepare_costs_for_r_solver(const lar_term & term) {        
-        
     TRACE("lar_solver", print_term(term, tout << "prepare: ") << "\n";);
+    if (move_non_basic_columns_to_bounds())
+        find_feasible_solution();
     auto & rslv = m_mpq_lar_core_solver.m_r_solver;
     rslv.m_using_infeas_costs = false;
     lp_assert(costs_are_zeros_for_r_solver());
@@ -471,13 +472,71 @@ void lar_solver::prepare_costs_for_r_solver(const lar_term & term) {
     }
     lp_assert(rslv.reduced_costs_are_correct_tableau());
 }
+
+bool lar_solver::move_non_basic_columns_to_bounds() {
+    auto & lcs = m_mpq_lar_core_solver;
+    bool change = false;
+    for (unsigned j : lcs.m_r_nbasis) {
+        if (move_non_basic_column_to_bounds(j))
+            change = true;
+    }
+
+    if (settings().simplex_strategy() == simplex_strategy_enum::tableau_costs)
+        update_x_and_inf_costs_for_columns_with_changed_bounds_tableau();
+    return change;
+}
+
+bool lar_solver::move_non_basic_column_to_bounds(unsigned j) {
+    auto & lcs = m_mpq_lar_core_solver;
+    auto & val = lcs.m_r_x[j];
+    switch (lcs.m_column_types()[j]) {
+    case column_type::boxed:
+        if (val != lcs.m_r_lower_bounds()[j] && val != lcs.m_r_upper_bounds()[j]) {
+            if (random() % 2 == 0)
+                set_value_for_nbasic_column(j, lcs.m_r_lower_bounds()[j]);
+            else
+                set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
+            return true;
+        }
+        break;
+    case column_type::lower_bound:
+        if (val != lcs.m_r_lower_bounds()[j]) {
+            set_value_for_nbasic_column(j, lcs.m_r_lower_bounds()[j]);
+            return true;
+        }
+        break;
+    case column_type::upper_bound:
+        if (val != lcs.m_r_upper_bounds()[j]) {
+            set_value_for_nbasic_column(j, lcs.m_r_upper_bounds()[j]);
+            return true;
+        }
+        break;
+    default:
+        if (is_int(j) && !val.is_int()) {
+            set_value_for_nbasic_column(j, impq(floor(val)));
+            return true;
+        }
+        break;
+    }
+    return false;
+}
+
+void lar_solver::set_value_for_nbasic_column(unsigned j, const impq & new_val) {
+    lp_assert(!is_base(j));
+    auto & x = m_mpq_lar_core_solver.m_r_x[j];
+    auto delta = new_val - x;
+    x = new_val;
+    change_basic_columns_dependend_on_a_given_nb_column(j, delta);
+}
+
     
 bool lar_solver::maximize_term_on_corrected_r_solver(lar_term & term,
                                                      impq &term_max) {
     settings().backup_costs = false;
     bool ret = false;
-    TRACE("lar_solver", print_term(term, tout << "maximize: ") << "\n"; print_constraints(tout););
+    TRACE("lar_solver", print_term(term, tout << "maximize: ") << "\n"; print_constraints(tout); tout << ", strategy = " << (int)settings().simplex_strategy() << "\n";);
     switch (settings().simplex_strategy()) {
+        
     case simplex_strategy_enum::tableau_rows:
         prepare_costs_for_r_solver(term);
         settings().simplex_strategy() = simplex_strategy_enum::tableau_costs;

src/util/lp/lar_solver.h

@@ -359,7 +359,14 @@ public:
     void detect_rows_with_changed_bounds_for_column(unsigned j);
     
     void detect_rows_with_changed_bounds();
+    inline bool is_base(unsigned j) const {
+        return m_mpq_lar_core_solver.m_r_heading[j] >= 0;
+    }
 
+    bool move_non_basic_columns_to_bounds();
+
+    bool move_non_basic_column_to_bounds(unsigned j);
+    void set_value_for_nbasic_column(unsigned j, const impq & new_val);
     void update_x_and_inf_costs_for_columns_with_changed_bounds();
 
     void update_x_and_inf_costs_for_columns_with_changed_bounds_tableau();

src/util/lp/lp_primal_core_solver.h

@@ -869,7 +869,8 @@ public:
         for (unsigned j : this->m_basis)
             if(!basis_column_is_set_correctly(j))
                 return false;
-        return true;
+        
+        return this->m_basis_heading.size() == this->m_A.column_count() && this->m_basis.size() == this->m_A.row_count();
     }
 
     void init_run_tableau();

src/util/lp/lp_primal_core_solver_def.h

@@ -151,8 +151,9 @@ template <typename T, typename X>
 bool lp_primal_core_solver<T, X>::column_is_benefitial_for_entering_basis_precise(unsigned j) const {
     lp_assert (numeric_traits<T>::precise());
     if (this->m_using_infeas_costs && this->m_settings.use_breakpoints_in_feasibility_search)
-        return column_is_benefitial_for_entering_on_breakpoints(j);
+        return column_is_benefitial_for_entering_on_breakpoints(j);    
     const T& dj = this->m_d[j];
+    TRACE("lar_solver", tout << "dj=" << dj << "\n";); 
     switch (this->m_column_types[j]) {
     case column_type::fixed:  break;
     case column_type::free_column: