1 /*++ 2 Copyright (c) 2014 Microsoft Corporation 3 4 Module Name: 5 6 simplex.h 7 8 Abstract: 9 10 Multi-precision simplex tableau. 11 12 - It uses code from theory_arith where applicable. 13 14 - It is detached from the theory class and ASTs. 15 16 - It uses non-shared mpz/mpq's avoiding global locks and operations on rationals. 17 18 - It follows the same sparse tableau layout (no LU yet). 19 20 - It does not include features for non-linear arithmetic. 21 22 - Branch/bound/cuts is external. 23 24 Author: 25 26 Nikolaj Bjorner (nbjorner) 2014-01-15 27 28 Notes: 29 30 --*/ 31 32 #pragma once 33 34 #include "math/simplex/sparse_matrix.h" 35 #include "util/mpq_inf.h" 36 #include "util/heap.h" 37 #include "util/lbool.h" 38 #include "util/uint_set.h" 39 40 namespace simplex { 41 42 template<typename Ext> 43 class simplex { 44 45 typedef unsigned var_t; 46 typedef typename Ext::eps_numeral eps_numeral; 47 typedef typename Ext::numeral numeral; 48 typedef typename Ext::manager manager; 49 typedef typename Ext::eps_manager eps_manager; 50 typedef typename Ext::scoped_numeral scoped_numeral; 51 typedef _scoped_numeral<eps_manager> scoped_eps_numeral; 52 typedef _scoped_numeral_vector<eps_manager> scoped_eps_numeral_vector; 53 typedef sparse_matrix<Ext> matrix; 54 struct var_lt { operatorvar_lt55 bool operator()(var_t v1, var_t v2) const { return v1 < v2; } 56 }; 57 typedef heap<var_lt> var_heap; 58 59 struct stats { 60 unsigned m_num_pivots; 61 unsigned m_num_infeasible; 62 unsigned m_num_checks; statsstats63 stats() { reset(); } resetstats64 void reset() { 65 memset(this, 0, sizeof(*this)); 66 } 67 }; 68 69 enum pivot_strategy_t { 70 S_BLAND, 71 S_GREATEST_ERROR, 72 S_LEAST_ERROR, 73 S_DEFAULT 74 }; 75 76 struct var_info { 77 unsigned m_base2row:29; 78 unsigned m_is_base:1; 79 unsigned m_lower_valid:1; 80 unsigned m_upper_valid:1; 81 eps_numeral m_value; 82 eps_numeral m_lower; 83 eps_numeral m_upper; 84 numeral m_base_coeff; var_infovar_info85 var_info(): 86 m_base2row(0), 87 m_is_base(false), 88 m_lower_valid(false), 89 m_upper_valid(false) 90 {} 91 }; 92 93 static const var_t null_var; 94 reslimit& m_limit; 95 mutable manager m; 96 mutable eps_manager em; 97 mutable matrix M; 98 unsigned m_max_iterations; 99 var_heap m_to_patch; 100 vector<var_info> m_vars; 101 svector<var_t> m_row2base; 102 bool m_bland; 103 unsigned m_blands_rule_threshold; 104 random_gen m_random; 105 uint_set m_left_basis; 106 unsigned m_infeasible_var; 107 unsigned_vector m_base_vars; 108 stats m_stats; 109 110 public: simplex(reslimit & lim)111 simplex(reslimit& lim): 112 m_limit(lim), 113 M(m), 114 m_max_iterations(UINT_MAX), 115 m_to_patch(1024), 116 m_bland(false), 117 m_blands_rule_threshold(1000) {} 118 119 ~simplex(); 120 121 typedef typename matrix::row row; 122 typedef typename matrix::row_iterator row_iterator; 123 typedef typename matrix::col_iterator col_iterator; 124 125 void ensure_var(var_t v); 126 row add_row(var_t base, unsigned num_vars, var_t const* vars, numeral const* coeffs); 127 row get_infeasible_row(); get_base_var(row const & r)128 var_t get_base_var(row const& r) const { return m_row2base[r.id()]; } get_base_coeff(row const & r)129 numeral const& get_base_coeff(row const& r) const { return m_vars[m_row2base[r.id()]].m_base_coeff; } 130 void del_row(var_t base_var); 131 void set_lower(var_t var, eps_numeral const& b); 132 void set_upper(var_t var, eps_numeral const& b); get_lower(var_t var,scoped_eps_numeral & b)133 void get_lower(var_t var, scoped_eps_numeral& b) const { b = m_vars[var].m_lower; } get_upper(var_t var,scoped_eps_numeral & b)134 void get_upper(var_t var, scoped_eps_numeral& b) const { b = m_vars[var].m_upper; } get_lower(var_t var)135 eps_numeral const& get_lower(var_t var) const { return m_vars[var].m_lower; } get_upper(var_t var)136 eps_numeral const& get_upper(var_t var) const { return m_vars[var].m_upper; } 137 bool above_lower(var_t var, eps_numeral const& b) const; 138 bool below_upper(var_t var, eps_numeral const& b) const; 139 bool below_lower(var_t v) const; 140 bool above_upper(var_t v) const; lower_valid(var_t var)141 bool lower_valid(var_t var) const { return m_vars[var].m_lower_valid; } upper_valid(var_t var)142 bool upper_valid(var_t var) const { return m_vars[var].m_upper_valid; } 143 void unset_lower(var_t var); 144 void unset_upper(var_t var); 145 void set_value(var_t var, eps_numeral const& b); set_max_iterations(unsigned n)146 void set_max_iterations(unsigned n) { m_max_iterations = n; } 147 void reset(); 148 lbool make_feasible(); 149 lbool minimize(var_t var); 150 eps_numeral const& get_value(var_t v); 151 void display(std::ostream& out) const; 152 void display_row(std::ostream& out, row const& r, bool values = true); 153 get_num_vars()154 unsigned get_num_vars() const { return m_vars.size(); } 155 row_begin(row const & r)156 row_iterator row_begin(row const& r) { return M.row_begin(r); } row_end(row const & r)157 row_iterator row_end(row const& r) { return M.row_end(r); } 158 159 void collect_statistics(::statistics & st) const; 160 161 private: 162 163 void del_row(row const& r); 164 var_t select_var_to_fix(); 165 pivot_strategy_t pivot_strategy(); select_smallest_var()166 var_t select_smallest_var() { return m_to_patch.empty()?null_var:m_to_patch.erase_min(); } 167 var_t select_error_var(bool least); 168 void check_blands_rule(var_t v, unsigned& num_repeated); 169 bool make_var_feasible(var_t x_i); 170 void update_and_pivot(var_t x_i, var_t x_j, numeral const& a_ij, eps_numeral const& new_value); 171 void update_value(var_t v, eps_numeral const& delta); 172 void update_value_core(var_t v, eps_numeral const& delta); 173 void pivot(var_t x_i, var_t x_j, numeral const& a_ij); 174 void move_to_bound(var_t x, bool to_lower); 175 var_t select_pivot(var_t x_i, bool is_below, scoped_numeral& out_a_ij); 176 var_t select_pivot_blands(var_t x_i, bool is_below, scoped_numeral& out_a_ij); 177 var_t select_pivot_core(var_t x_i, bool is_below, scoped_numeral& out_a_ij); 178 int get_num_non_free_dep_vars(var_t x_j, int best_so_far); 179 180 var_t pick_var_to_leave(var_t x_j, bool is_pos, 181 scoped_eps_numeral& gain, scoped_numeral& new_a_ij, bool& inc); 182 183 184 void select_pivot_primal(var_t v, var_t& x_i, var_t& x_j, scoped_numeral& a_ij, bool& inc_x_i, bool& inc_x_j); 185 186 187 bool at_lower(var_t v) const; 188 bool at_upper(var_t v) const; 189 bool above_lower(var_t v) const; 190 bool below_upper(var_t v) const; outside_bounds(var_t v)191 bool outside_bounds(var_t v) const { return below_lower(v) || above_upper(v); } is_free(var_t v)192 bool is_free(var_t v) const { return !m_vars[v].m_lower_valid && !m_vars[v].m_upper_valid; } is_non_free(var_t v)193 bool is_non_free(var_t v) const { return !is_free(v); } is_base(var_t x)194 bool is_base(var_t x) const { return m_vars[x].m_is_base; } 195 void add_patch(var_t v); 196 197 bool well_formed() const; 198 bool well_formed_row(row const& r) const; 199 bool is_feasible() const; 200 }; 201 202 void ensure_rational_solution(simplex<mpq_ext>& s); 203 }; 204 205