Remove non feasible costs

src/math/lp/lar_core_solver_def.h

@@ -40,7 +40,6 @@ void lar_core_solver::prefix_r() {
     if (m_r_solver.m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows) {
         m_r_solver.m_costs.resize(m_r_solver.m_n());
         m_r_solver.m_d.resize(m_r_solver.m_n());
-        m_r_solver.set_using_infeas_costs(true);
     }
 }
 

src/math/lp/lar_solver.cpp

@@ -191,9 +191,11 @@ namespace lp {
             stats().m_max_rows = A_r().row_count();
         if (strategy_is_undecided())
             decide_on_strategy_and_adjust_initial_state();
-
+        auto strategy_was =     settings().simplex_strategy();
+        settings().set_simplex_strategy(simplex_strategy_enum::tableau_rows);
         m_mpq_lar_core_solver.m_r_solver.m_look_for_feasible_solution_only = true;
         auto ret = solve();
+        settings().set_simplex_strategy(strategy_was);
         return ret;
     }
 
@@ -355,7 +357,6 @@ namespace lp {
         m_basic_columns_with_changed_cost.resize(m_mpq_lar_core_solver.m_r_x.size());
         move_non_basic_columns_to_bounds(false);
         auto& rslv = m_mpq_lar_core_solver.m_r_solver;
-        rslv.set_using_infeas_costs(false);
         lp_assert(costs_are_zeros_for_r_solver());
         lp_assert(reduced_costs_are_zeroes_for_r_solver());
         rslv.m_costs.resize(A_r().column_count(), zero_of_type<mpq>());
@@ -490,11 +491,11 @@ namespace lp {
     lp_status lar_solver::maximize_term(unsigned j_or_term,
         impq& term_max) {
         TRACE("lar_solver", print_values(tout););
+
         lar_term term = get_term_to_maximize(j_or_term);
         if (term.is_empty()) {
             return lp_status::UNBOUNDED;
         }
-
         impq prev_value;
         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()) {
@@ -710,14 +711,6 @@ namespace lp {
     void lar_solver::update_x_and_inf_costs_for_columns_with_changed_bounds_tableau() {
         for (auto j : m_columns_with_changed_bounds)
             update_x_and_inf_costs_for_column_with_changed_bounds(j);
-
-        if (tableau_with_costs()) {
-            if (m_mpq_lar_core_solver.m_r_solver.using_infeas_costs()) {
-                for (unsigned j : m_basic_columns_with_changed_cost)
-                    m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j);
-                lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
-            }
-        }
     }
 
 
@@ -1344,14 +1337,6 @@ namespace lp {
         for (unsigned j : became_feas)
             m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j);
 
-
-        if (use_tableau_costs()) {
-            for (unsigned j : became_feas)
-                m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j);
-            for (unsigned j : basic_columns_with_changed_cost)
-                m_mpq_lar_core_solver.m_r_solver.update_inf_cost_for_column_tableau(j);
-            lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
-        }
     }
 
     bool lar_solver::model_is_int_feasible() const {

src/math/lp/lp_core_solver_base.cpp

@@ -65,8 +65,6 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::pivot_column_tableau(un
 template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::transpose_rows_tableau(unsigned int, unsigned int);
 template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::inf_set_is_correct() const;
 template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::inf_set_is_correct() const;
-template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::infeasibility_costs_are_correct() const;
-template bool lp::lp_core_solver_base<lp::mpq, lp::mpq >::infeasibility_costs_are_correct() const;
 template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_from_basis(unsigned int);
 
 

src/math/lp/lp_core_solver_base.h

@@ -63,15 +63,12 @@ class lp_core_solver_base {
     bool current_x_is_infeasible() const { return m_inf_set.size() != 0; }
 private:
     u_set m_inf_set;
-    bool m_using_infeas_costs;
 public:
     const u_set& inf_set() const { return m_inf_set; }
     u_set& inf_set() { return m_inf_set; }
     void inf_set_increase_size_by_one() { m_inf_set.increase_size_by_one(); }
     bool inf_set_contains(unsigned j) const { return m_inf_set.contains(j); }
-    unsigned inf_set_size() const { return m_inf_set.size(); }
-    bool using_infeas_costs() const { return m_using_infeas_costs; }
-    void set_using_infeas_costs(bool val)  { m_using_infeas_costs = val; }
+    unsigned inf_set_size() const { return m_inf_set.size(); }    
     indexed_vector<T>     m_pivot_row; // this is the real pivot row of the simplex tableu
     static_matrix<T, X> & m_A; // the matrix A
     // vector<X> const &           m_b; // the right side
@@ -198,11 +195,7 @@ class lp_core_solver_base {
         if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows)
             return true;
         CASSERT("check_static_matrix", m_A.is_correct());
-        if (m_using_infeas_costs) {
-            if (infeasibility_costs_are_correct() == false) {
-                return false;
-            }
-        }
+
 
         unsigned n = m_A.column_count();
         for (unsigned j = 0; j < n; j++) {
@@ -236,8 +229,6 @@ class lp_core_solver_base {
         return below_bound(m_x[p], m_lower_bounds[p]);
     }
 
-    bool infeasibility_costs_are_correct() const;
-    bool infeasibility_cost_is_correct_for_column(unsigned j) const;
 
     bool x_above_lower_bound(unsigned p) const {
         return above_bound(m_x[p], m_lower_bounds[p]);

src/math/lp/lp_core_solver_base_def.h

@@ -43,7 +43,6 @@ lp_core_solver_base(static_matrix<T, X> & A,
     m_iters_with_no_cost_growing(0),
     m_status(lp_status::FEASIBLE),
     m_inf_set(A.column_count()),
-    m_using_infeas_costs(false),
     m_pivot_row(A.column_count()),
     m_A(A),
     m_basis(basis),
@@ -94,8 +93,6 @@ pivot_to_reduced_costs_tableau(unsigned i, unsigned j) {
 
 
 
-
-
 // template <typename T, typename X> void lp_core_solver_base<T, X>::
 // update_index_of_ed() {
 //     m_index_of_ed.clear();
@@ -434,57 +431,4 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::remove_from_ba
 }
 
 
-template <typename T, typename X> bool 
-lp_core_solver_base<T, X>::infeasibility_costs_are_correct() const {
-    if (! this->m_using_infeas_costs)
-        return true;
-    lp_assert(costs_on_nbasis_are_zeros());
-    for (unsigned j :this->m_basis) {
-        if (!infeasibility_cost_is_correct_for_column(j)) {
-            TRACE("lar_solver", tout << "incorrect cost for column " << j << std::endl;);
-            return false;
-        }
-        if (!is_zero(m_d[j])) {
-            TRACE("lar_solver", tout << "non zero inf cost for basis j = " << j << std::endl;);
-            return false;
-        }
-    }
-    return true;
-}
-
-template <typename T, typename X> bool
-lp_core_solver_base<T, X>::infeasibility_cost_is_correct_for_column(unsigned j)  const {
-    T r =  -one_of_type<T>();
-
-    switch (this->m_column_types[j]) {
-    case column_type::fixed:
-    case column_type::boxed:
-        if (this->x_above_upper_bound(j)) {
-            return (this->m_costs[j] == r);
-        }
-        if (this->x_below_low_bound(j)) {
-            return (this->m_costs[j] == -r);
-        }
-        return is_zero(this->m_costs[j]);
-
-    case column_type::lower_bound:
-        if (this->x_below_low_bound(j)) {
-            return this->m_costs[j] == -r;
-        }
-        return is_zero(this->m_costs[j]);
-
-    case column_type::upper_bound:
-        if (this->x_above_upper_bound(j)) {
-            return this->m_costs[j] == r;
-        }
-        return is_zero(this->m_costs[j]);
-    case column_type::free_column:
-        return is_zero(this->m_costs[j]);
-    default:
-        UNREACHABLE();
-        return true;
-    }
-}
-
-
 }

src/math/lp/lp_primal_core_solver.cpp

@@ -33,6 +33,5 @@ template unsigned lp_primal_core_solver<mpq, mpq>::solve();
 template unsigned lp_primal_core_solver<mpq, numeric_pair<mpq> >::solve();
 template bool lp::lp_primal_core_solver<lp::mpq, lp::mpq>::update_basis_and_x_tableau(int, int, lp::mpq const&);
 template bool lp::lp_primal_core_solver<lp::mpq, lp::numeric_pair<lp::mpq> >::update_basis_and_x_tableau(int, int, lp::numeric_pair<lp::mpq> const&);
-template void lp::lp_primal_core_solver<rational, lp::numeric_pair<rational> >::update_inf_cost_for_column_tableau(unsigned);
 
 }

src/math/lp/lp_primal_core_solver.h

@@ -334,28 +334,12 @@ class lp_primal_core_solver : public lp_core_solver_base<T, X> {
   void update_x_tableau_rows(unsigned entering, unsigned leaving,
                              const X &delta) {
     this->add_delta_to_x(entering, delta);
-    if (!this->using_infeas_costs()) {
-      for (const auto &c : this->m_A.m_columns[entering]) {
-        if (leaving != this->m_basis[c.var()]) {
-          this->add_delta_to_x_and_track_feasibility(
-              this->m_basis[c.var()], -delta * this->m_A.get_val(c));
-        }
-      }
-    } else { // using_infeas_costs() == true
-      lp_assert(this->column_is_feasible(entering));
-      lp_assert(this->m_costs[entering] == zero_of_type<T>());
-      // m_d[entering] can change because of the cost change for basic columns.
-      for (const auto &c : this->m_A.m_columns[entering]) {
-        unsigned j = this->m_basis[c.var()];
-        if (j != leaving)
-          this->add_delta_to_x(j, -delta * this->m_A.get_val(c));
-        update_inf_cost_for_column_tableau(j);
-        if (is_zero(this->m_costs[j]))
-          this->remove_column_from_inf_set(j);
-        else
-          this->insert_column_into_inf_set(j);
-      }
-    }
+    for (const auto &c : this->m_A.m_columns[entering]) {
+       if (leaving != this->m_basis[c.var()]) {
+         this->add_delta_to_x_and_track_feasibility(
+         this->m_basis[c.var()], -delta * this->m_A.get_val(c));
+       }
+    }    
   }
 
   void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) {
@@ -437,7 +421,6 @@ class lp_primal_core_solver : public lp_core_solver_base<T, X> {
   }
 
   void decide_on_status_when_cannot_find_entering() {
-    lp_assert(!need_to_switch_costs());
     this->set_status(this->current_x_is_feasible() ? lp_status::OPTIMAL
                                                    : lp_status::INFEASIBLE);
   }
@@ -628,11 +611,6 @@ class lp_primal_core_solver : public lp_core_solver_base<T, X> {
 
   bool column_is_benefitial_for_entering_basis(unsigned j) const;
   void init_infeasibility_costs();
-
-  void init_infeasibility_cost_for_column(unsigned j);
-  T get_infeasibility_cost_for_column(unsigned j) const;
-  void init_infeasibility_costs_for_changed_basis_only();
-
   void print_column(unsigned j, std::ostream &out);
 
   void print_bound_info_and_x(unsigned j, std::ostream &out);
@@ -652,8 +630,6 @@ class lp_primal_core_solver : public lp_core_solver_base<T, X> {
 
   void init_run_tableau();
   void update_x_tableau(unsigned entering, const X &delta);
-  void update_inf_cost_for_column_tableau(unsigned j);
-
   // the delta is between the old and the new cost (old - new)
   void update_reduced_cost_for_basic_column_cost_change(const T &delta,
                                                         unsigned j) {

src/math/lp/lp_primal_core_solver_def.h

@@ -258,122 +258,10 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::find_feas
     solve();
 }
 
-template <typename T, typename X>
-void lp_primal_core_solver<T, X>::init_infeasibility_costs_for_changed_basis_only() {
-    for (unsigned i :  this->m_ed.m_index)
-        init_infeasibility_cost_for_column(this->m_basis[i]);
-    this->set_using_infeas_costs(true);
-}
-
 
-template <typename T, typename X>
-void lp_primal_core_solver<T, X>::init_infeasibility_costs() {
-    lp_assert(this->m_x.size() >= this->m_n());
-    lp_assert(this->m_column_types.size() >= this->m_n());
-    for (unsigned j = this->m_n(); j--;)
-        init_infeasibility_cost_for_column(j);
-    this->set_using_infeas_costs(true);
-}
 
-template <typename T, typename X> T
-lp_primal_core_solver<T, X>::get_infeasibility_cost_for_column(unsigned j) const {
-    if (this->m_basis_heading[j] < 0) {
-        return zero_of_type<T>();
-    }
-    T ret;
-    // j is a basis column
-    switch (this->m_column_types[j]) {
-    case column_type::fixed:
-    case column_type::boxed:
-        if (this->x_above_upper_bound(j)) {
-            ret = 1;
-        } else if (this->x_below_low_bound(j)) {
-            ret = -1;
-        } else {
-            ret = numeric_traits<T>::zero();
-        }
-        break;
-    case column_type::lower_bound:
-        if (this->x_below_low_bound(j)) {
-            ret = -1;
-        } else {
-            ret = numeric_traits<T>::zero();
-        }
-        break;
-    case column_type::upper_bound:
-        if (this->x_above_upper_bound(j)) {
-            ret = 1;
-        } else {
-            ret = numeric_traits<T>::zero();
-        }
-        break;
-    case column_type::free_column:
-        ret = numeric_traits<T>::zero();
-        break;
-    default:
-        UNREACHABLE();
-        ret = numeric_traits<T>::zero(); // does not matter
-        break;
-    }
-
-    ret = - ret;
-
-    return ret;
-}
 
 
-// changed m_inf_set too!
-template <typename T, typename X> void
-lp_primal_core_solver<T, X>::init_infeasibility_cost_for_column(unsigned j) {
-
-    if (this->m_basis_heading[j] < 0) {
-        this->m_costs[j] = numeric_traits<T>::zero();
-        this->remove_column_from_inf_set(j);
-        return;
-    }
-    // j is a basis column
-    switch (this->m_column_types[j]) {
-    case column_type::fixed:
-    case column_type::boxed:
-        if (this->x_above_upper_bound(j)) {
-            this->m_costs[j] = 1;
-        } else if (this->x_below_low_bound(j)) {
-            this->m_costs[j] = -1;
-        } else {
-            this->m_costs[j] = numeric_traits<T>::zero();
-        }
-        break;
-    case column_type::lower_bound:
-        if (this->x_below_low_bound(j)) {
-            this->m_costs[j] = -1;
-        } else {
-            this->m_costs[j] = numeric_traits<T>::zero();
-        }
-        break;
-    case column_type::upper_bound:
-        if (this->x_above_upper_bound(j)) {
-            this->m_costs[j] = 1;
-        } else {
-            this->m_costs[j] = numeric_traits<T>::zero();
-        }
-        break;
-    case column_type::free_column:
-        this->m_costs[j] = numeric_traits<T>::zero();
-        break;
-    default:
-        UNREACHABLE();
-        break;
-    }
-
-    if (numeric_traits<T>::is_zero(this->m_costs[j])) {
-        this->remove_column_from_inf_set(j);
-    } else {
-        this->insert_column_into_inf_set(j);
-    }
-    this->m_costs[j] = - this->m_costs[j];
-
-}
-
 
 template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column(unsigned j, std::ostream & out) {
     out << this->column_name(j) << " " <<   j << " " << column_type_to_string(this->m_column_type[j]) << " x = " << this->m_x[j] << " " << "c = " << this->m_costs[j];;